From 30fd78b3839c87d443f0b8a997922cb3729d143b Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Wed, 27 Sep 2017 14:33:48 +0200 Subject: [PATCH] Brusselator test with jacobian. --- tests/sundials/harmonic_oscillator_05.cc | 133 +++++++++++++++++ tests/sundials/harmonic_oscillator_05.output | 102 +++++++++++++ tests/sundials/harmonic_oscillator_05.prm | 15 ++ tests/sundials/harmonic_oscillator_06.cc | 145 +++++++++++++++++++ tests/sundials/harmonic_oscillator_06.output | 102 +++++++++++++ tests/sundials/harmonic_oscillator_06.prm | 15 ++ 6 files changed, 512 insertions(+) create mode 100644 tests/sundials/harmonic_oscillator_05.cc create mode 100644 tests/sundials/harmonic_oscillator_05.output create mode 100644 tests/sundials/harmonic_oscillator_05.prm create mode 100644 tests/sundials/harmonic_oscillator_06.cc create mode 100644 tests/sundials/harmonic_oscillator_06.output create mode 100644 tests/sundials/harmonic_oscillator_06.prm diff --git a/tests/sundials/harmonic_oscillator_05.cc b/tests/sundials/harmonic_oscillator_05.cc new file mode 100644 index 0000000000..27d2f81b78 --- /dev/null +++ b/tests/sundials/harmonic_oscillator_05.cc @@ -0,0 +1,133 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2015 by the deal2lkit authors +// +// This file is part of the deal2lkit library. +// +// The deal2lkit library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal2lkit distribution. +// +//----------------------------------------------------------- + +#include "../tests.h" +#include +#include +#include +#include + + +// Test implicit-explicit time stepper. Only solve_jacobian_system. +// Brusselator benchmark + +/** + * This test problem is called "brusselator", and is a typical benchmark for + * ODE solvers. This problem has 3 dependent variables u, v and w, that depend + * on the independent variable t via the IVP system + * + * du/dt = a − (w + 1)u + v u^2 + * dv/dt = w u − v u^2 + * dw/dt = (b − w)/eps -w u + * + * We integrate over the interval 0 ≤ t ≤ 10, with the initial conditions + * + * u(0) = 3.9, v(0) = 1.1, w(0) = 2.8, + * + * and parameters + * + * a = 1.2, b = 2.5, and eps = 10−5 + * + * The implicit part only contains the stiff part of the problem (the part with + * eps in right hand side of the third equation). + */ +int main (int argc, char **argv) +{ + std::ofstream out("output"); + + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int); + + typedef Vector VectorType; + + ParameterHandler prm; + SUNDIALS::ARKode::AdditionalData data; + data.add_parameters(prm); + + if (false) + { + std::ofstream ofile(SOURCE_DIR "/harmonic_oscillator_05.prm"); + prm.print_parameters(ofile, ParameterHandler::ShortText); + ofile.close(); + } + + std::ifstream ifile(SOURCE_DIR "/harmonic_oscillator_05.prm"); + prm.parse_input(ifile); + + SUNDIALS::ARKode ode(MPI_COMM_WORLD,data); + + ode.reinit_vector = [&] (VectorType&v) + { + // Three independent variables + v.reinit(3); + }; + + // Parameters + double u0 = 3.9, v0 = 1.1, w0 = 2.8, a = 1.2, b= 2.5, eps = 1e-5; + // Explicit jacobian. + FullMatrix J(3,3); + + ode.implicit_function = [&] (double, + const VectorType &y, + VectorType &ydot) -> int + { + ydot[0] = 0; + ydot[1] = 0; + ydot[2] = (b-y[2])/eps; + return 0; + }; + + + ode.explicit_function = [&] (double, + const VectorType &y, + VectorType &ydot) -> int + { + ydot[0] = a-(y[2]+1)*y[0]+y[1]*y[0]*y[0]; + ydot[1] = y[2]*y[0]-y[1]*y[0]*y[0]; + ydot[2] = -y[2]*y[0]; + return 0; + }; + + ode.solve_jacobian_system = [&] (const double t, + const double gamma, + const VectorType &, + const VectorType &, + const VectorType &src, + VectorType &dst) -> int + { + J = 0; + J(0,0) = 1; + J(1,1) = 1; + J(2,2) = 1 + gamma/eps; + J.gauss_jordan(); + J.vmult(dst, src); + return 0; + }; + + ode.output_step = [&](const double t, + const VectorType &sol, + const unsigned int step_number) -> int + { + out << t << " " + << sol[0] << " " << sol[1] << " " << sol[2] << std::endl; + return 0; + }; + + Vector y(3); + y[0] = u0; + y[1] = v0; + y[2] = w0; + ode.solve_ode(y); + return 0; +} diff --git a/tests/sundials/harmonic_oscillator_05.output b/tests/sundials/harmonic_oscillator_05.output new file mode 100644 index 0000000000..507044c3b2 --- /dev/null +++ b/tests/sundials/harmonic_oscillator_05.output @@ -0,0 +1,102 @@ +0 3.9 1.1 2.8 +0.1 3.99995 0.720431 2.49989 +0.2 3.78894 0.660967 2.49996 +0.3 3.52275 0.681501 2.49984 +0.4 3.26217 0.72296 2.49995 +0.5 3.01935 0.771864 2.49995 +0.6 2.79594 0.82467 2.49988 +0.7 2.59116 0.880241 2.49998 +0.8 2.4037 0.938097 2.49993 +0.9 2.23217 0.997959 2.49992 +1 2.07527 1.0596 2.49999 +1.1 1.93182 1.12281 2.49992 +1.2 1.80073 1.18737 2.49995 +1.3 1.68106 1.25304 2.49999 +1.4 1.57194 1.31959 2.49993 +1.5 1.47261 1.38677 2.49997 +1.6 1.38237 1.45434 2.49999 +1.7 1.3006 1.52203 2.49994 +1.8 1.2267 1.58963 2.49998 +1.9 1.16015 1.65689 2.49998 +2 1.10046 1.72361 2.49996 +2.1 1.04715 1.78959 2.49998 +2.2 0.999783 1.85466 2.49997 +2.3 0.957942 1.91866 2.49997 +2.4 0.921233 1.98145 2.49999 +2.5 0.889283 2.04291 2.49997 +2.6 0.861743 2.10294 2.49998 +2.7 0.838291 2.16142 2.49999 +2.8 0.818628 2.21827 2.49997 +2.9 0.802485 2.27338 2.49998 +3 0.78962 2.32667 2.49998 +3.1 0.779824 2.37802 2.49998 +3.2 0.772919 2.42731 2.49998 +3.3 0.768758 2.47441 2.49999 +3.4 0.767229 2.51916 2.49998 +3.5 0.768257 2.56138 2.5 +3.6 0.7718 2.60085 2.49995 +3.7 0.777859 2.63733 2.50005 +3.8 0.786474 2.67052 2.49999 +3.9 0.797733 2.70008 2.49988 +4 0.811773 2.72558 2.5001 +4.1 0.828789 2.74657 2.5001 +4.2 0.84904 2.76245 2.50009 +4.3 0.872862 2.77257 2.50012 +4.4 0.900673 2.77611 2.49986 +4.5 0.932987 2.77216 2.50009 +4.6 0.970423 2.75959 2.49999 +4.7 1.01371 2.73716 2.50002 +4.8 1.06366 2.70339 2.50006 +4.9 1.12115 2.65673 2.49999 +5 1.187 2.59554 2.49993 +5.1 1.26179 2.51839 2.50005 +5.2 1.34554 2.42435 2.49986 +5.3 1.43719 2.31362 2.50002 +5.4 1.53404 2.18823 2.50001 +5.5 1.6313 2.05269 2.4999 +5.6 1.72208 1.91414 2.49987 +5.7 1.79854 1.7815 2.49991 +5.8 1.85371 1.66352 2.49995 +5.9 1.88345 1.56671 2.49996 +6 1.88731 1.4941 2.49995 +6.1 1.86803 1.44543 2.49992 +6.2 1.83015 1.41827 2.49989 +6.3 1.77867 1.40922 2.50001 +6.4 1.71817 1.41482 2.50005 +6.5 1.65242 1.43201 2.50005 +6.6 1.58434 1.45824 2.50005 +6.7 1.51607 1.49149 2.49991 +6.8 1.44916 1.53017 2.4999 +6.9 1.38465 1.57301 2.50006 +7 1.32327 1.61902 2.49995 +7.1 1.26548 1.6674 2.49989 +7.2 1.21156 1.7175 2.50005 +7.3 1.16164 1.76879 2.49998 +7.4 1.11578 1.82081 2.49991 +7.5 1.07396 1.87319 2.50002 +7.6 1.03609 1.92558 2.5 +7.7 1.00208 1.97772 2.49993 +7.8 0.971784 2.02935 2.49999 +7.9 0.945072 2.08025 2.5 +8 0.921787 2.13022 2.49995 +8.1 0.901781 2.17907 2.49997 +8.2 0.884906 2.22664 2.5 +8.3 0.871025 2.27275 2.49997 +8.4 0.860013 2.31723 2.49997 +8.5 0.851761 2.35992 2.49999 +8.6 0.84618 2.40062 2.49998 +8.7 0.843199 2.43916 2.49998 +8.8 0.842774 2.47531 2.49999 +8.9 0.844885 2.50883 2.49997 +9 0.849543 2.53948 2.50004 +9.1 0.856788 2.56694 2.49994 +9.2 0.866696 2.59088 2.49999 +9.3 0.879382 2.61091 2.50009 +9.4 0.895002 2.6266 2.49992 +9.5 0.913755 2.63743 2.49985 +9.6 0.935893 2.64284 2.49987 +9.7 0.961713 2.64217 2.50002 +9.8 0.991567 2.63469 2.49996 +9.9 1.02585 2.61958 2.50009 +10 1.06497 2.59596 2.49997 +10 1.06497 2.59596 2.49997 diff --git a/tests/sundials/harmonic_oscillator_05.prm b/tests/sundials/harmonic_oscillator_05.prm new file mode 100644 index 0000000000..55a91a0b7c --- /dev/null +++ b/tests/sundials/harmonic_oscillator_05.prm @@ -0,0 +1,15 @@ +set Final time = 10. +set Initial time = 0. +set Time interval between each output = 0.1 +subsection Error control + set Absolute error tolerance = 0.000001 + set Relative error tolerance = 0.000010 +end +subsection Running parameters + set Implicit function is linear = true + set Implicit function is time independent = true + set Initial step size = 0.010000 + set Maximum number of nonlinear iterations = 10 + set Maximum order of ARK = 5 + set Minimum step size = 0.000001 +end diff --git a/tests/sundials/harmonic_oscillator_06.cc b/tests/sundials/harmonic_oscillator_06.cc new file mode 100644 index 0000000000..c58e56c0c1 --- /dev/null +++ b/tests/sundials/harmonic_oscillator_06.cc @@ -0,0 +1,145 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2015 by the deal2lkit authors +// +// This file is part of the deal2lkit library. +// +// The deal2lkit library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal2lkit distribution. +// +//----------------------------------------------------------- + +#include "../tests.h" +#include +#include +#include +#include + + +// Test implicit-explicit time stepper. Both setup and solve_jacobian_system. +// Brusselator benchmark + +/** + * This test problem is called "brusselator", and is a typical benchmark for + * ODE solvers. This problem has 3 dependent variables u, v and w, that depend + * on the independent variable t via the IVP system + * + * du/dt = a − (w + 1)u + v u^2 + * dv/dt = w u − v u^2 + * dw/dt = (b − w)/eps -w u + * + * We integrate over the interval 0 ≤ t ≤ 10, with the initial conditions + * + * u(0) = 3.9, v(0) = 1.1, w(0) = 2.8, + * + * and parameters + * + * a = 1.2, b = 2.5, and eps = 10−5 + * + * The implicit part only contains the stiff part of the problem (the part with + * eps in right hand side of the third equation). + */ +int main (int argc, char **argv) +{ + std::ofstream out("output"); + + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int); + + typedef Vector VectorType; + + ParameterHandler prm; + SUNDIALS::ARKode::AdditionalData data; + data.add_parameters(prm); + + if (false) + { + std::ofstream ofile(SOURCE_DIR "/harmonic_oscillator_06.prm"); + prm.print_parameters(ofile, ParameterHandler::ShortText); + ofile.close(); + } + + std::ifstream ifile(SOURCE_DIR "/harmonic_oscillator_06.prm"); + prm.parse_input(ifile); + + SUNDIALS::ARKode ode(MPI_COMM_WORLD,data); + + ode.reinit_vector = [&] (VectorType&v) + { + // Three independent variables + v.reinit(3); + }; + + // Parameters + double u0 = 3.9, v0 = 1.1, w0 = 2.8, a = 1.2, b= 2.5, eps = 1e-5; + // Explicit jacobian. + FullMatrix J(3,3); + + ode.implicit_function = [&] (double, + const VectorType &y, + VectorType &ydot) -> int + { + ydot[0] = 0; + ydot[1] = 0; + ydot[2] = (b-y[2])/eps; + return 0; + }; + + + ode.explicit_function = [&] (double, + const VectorType &y, + VectorType &ydot) -> int + { + ydot[0] = a-(y[2]+1)*y[0]+y[1]*y[0]*y[0]; + ydot[1] = y[2]*y[0]-y[1]*y[0]*y[0]; + ydot[2] = -y[2]*y[0]; + return 0; + }; + + + ode.setup_jacobian = [&] (const int, + const double, + const double gamma, + const VectorType &, + const VectorType &, + bool &j_is_current) -> int + { + J = 0; + J(0,0) = 1; + J(1,1) = 1; + J(2,2) = 1 + gamma/eps; + J.gauss_jordan(); + j_is_current = true; + return 0; + }; + + ode.solve_jacobian_system = [&] (const double t, + const double gamma, + const VectorType &, + const VectorType &, + const VectorType &src, + VectorType &dst) -> int + { + J.vmult(dst, src); + return 0; + }; + + ode.output_step = [&](const double t, + const VectorType &sol, + const unsigned int step_number) -> int + { + out << t << " " + << sol[0] << " " << sol[1] << " " << sol[2] << std::endl; + return 0; + }; + + Vector y(3); + y[0] = u0; + y[1] = v0; + y[2] = w0; + ode.solve_ode(y); + return 0; +} diff --git a/tests/sundials/harmonic_oscillator_06.output b/tests/sundials/harmonic_oscillator_06.output new file mode 100644 index 0000000000..507044c3b2 --- /dev/null +++ b/tests/sundials/harmonic_oscillator_06.output @@ -0,0 +1,102 @@ +0 3.9 1.1 2.8 +0.1 3.99995 0.720431 2.49989 +0.2 3.78894 0.660967 2.49996 +0.3 3.52275 0.681501 2.49984 +0.4 3.26217 0.72296 2.49995 +0.5 3.01935 0.771864 2.49995 +0.6 2.79594 0.82467 2.49988 +0.7 2.59116 0.880241 2.49998 +0.8 2.4037 0.938097 2.49993 +0.9 2.23217 0.997959 2.49992 +1 2.07527 1.0596 2.49999 +1.1 1.93182 1.12281 2.49992 +1.2 1.80073 1.18737 2.49995 +1.3 1.68106 1.25304 2.49999 +1.4 1.57194 1.31959 2.49993 +1.5 1.47261 1.38677 2.49997 +1.6 1.38237 1.45434 2.49999 +1.7 1.3006 1.52203 2.49994 +1.8 1.2267 1.58963 2.49998 +1.9 1.16015 1.65689 2.49998 +2 1.10046 1.72361 2.49996 +2.1 1.04715 1.78959 2.49998 +2.2 0.999783 1.85466 2.49997 +2.3 0.957942 1.91866 2.49997 +2.4 0.921233 1.98145 2.49999 +2.5 0.889283 2.04291 2.49997 +2.6 0.861743 2.10294 2.49998 +2.7 0.838291 2.16142 2.49999 +2.8 0.818628 2.21827 2.49997 +2.9 0.802485 2.27338 2.49998 +3 0.78962 2.32667 2.49998 +3.1 0.779824 2.37802 2.49998 +3.2 0.772919 2.42731 2.49998 +3.3 0.768758 2.47441 2.49999 +3.4 0.767229 2.51916 2.49998 +3.5 0.768257 2.56138 2.5 +3.6 0.7718 2.60085 2.49995 +3.7 0.777859 2.63733 2.50005 +3.8 0.786474 2.67052 2.49999 +3.9 0.797733 2.70008 2.49988 +4 0.811773 2.72558 2.5001 +4.1 0.828789 2.74657 2.5001 +4.2 0.84904 2.76245 2.50009 +4.3 0.872862 2.77257 2.50012 +4.4 0.900673 2.77611 2.49986 +4.5 0.932987 2.77216 2.50009 +4.6 0.970423 2.75959 2.49999 +4.7 1.01371 2.73716 2.50002 +4.8 1.06366 2.70339 2.50006 +4.9 1.12115 2.65673 2.49999 +5 1.187 2.59554 2.49993 +5.1 1.26179 2.51839 2.50005 +5.2 1.34554 2.42435 2.49986 +5.3 1.43719 2.31362 2.50002 +5.4 1.53404 2.18823 2.50001 +5.5 1.6313 2.05269 2.4999 +5.6 1.72208 1.91414 2.49987 +5.7 1.79854 1.7815 2.49991 +5.8 1.85371 1.66352 2.49995 +5.9 1.88345 1.56671 2.49996 +6 1.88731 1.4941 2.49995 +6.1 1.86803 1.44543 2.49992 +6.2 1.83015 1.41827 2.49989 +6.3 1.77867 1.40922 2.50001 +6.4 1.71817 1.41482 2.50005 +6.5 1.65242 1.43201 2.50005 +6.6 1.58434 1.45824 2.50005 +6.7 1.51607 1.49149 2.49991 +6.8 1.44916 1.53017 2.4999 +6.9 1.38465 1.57301 2.50006 +7 1.32327 1.61902 2.49995 +7.1 1.26548 1.6674 2.49989 +7.2 1.21156 1.7175 2.50005 +7.3 1.16164 1.76879 2.49998 +7.4 1.11578 1.82081 2.49991 +7.5 1.07396 1.87319 2.50002 +7.6 1.03609 1.92558 2.5 +7.7 1.00208 1.97772 2.49993 +7.8 0.971784 2.02935 2.49999 +7.9 0.945072 2.08025 2.5 +8 0.921787 2.13022 2.49995 +8.1 0.901781 2.17907 2.49997 +8.2 0.884906 2.22664 2.5 +8.3 0.871025 2.27275 2.49997 +8.4 0.860013 2.31723 2.49997 +8.5 0.851761 2.35992 2.49999 +8.6 0.84618 2.40062 2.49998 +8.7 0.843199 2.43916 2.49998 +8.8 0.842774 2.47531 2.49999 +8.9 0.844885 2.50883 2.49997 +9 0.849543 2.53948 2.50004 +9.1 0.856788 2.56694 2.49994 +9.2 0.866696 2.59088 2.49999 +9.3 0.879382 2.61091 2.50009 +9.4 0.895002 2.6266 2.49992 +9.5 0.913755 2.63743 2.49985 +9.6 0.935893 2.64284 2.49987 +9.7 0.961713 2.64217 2.50002 +9.8 0.991567 2.63469 2.49996 +9.9 1.02585 2.61958 2.50009 +10 1.06497 2.59596 2.49997 +10 1.06497 2.59596 2.49997 diff --git a/tests/sundials/harmonic_oscillator_06.prm b/tests/sundials/harmonic_oscillator_06.prm new file mode 100644 index 0000000000..55a91a0b7c --- /dev/null +++ b/tests/sundials/harmonic_oscillator_06.prm @@ -0,0 +1,15 @@ +set Final time = 10. +set Initial time = 0. +set Time interval between each output = 0.1 +subsection Error control + set Absolute error tolerance = 0.000001 + set Relative error tolerance = 0.000010 +end +subsection Running parameters + set Implicit function is linear = true + set Implicit function is time independent = true + set Initial step size = 0.010000 + set Maximum number of nonlinear iterations = 10 + set Maximum order of ARK = 5 + set Minimum step size = 0.000001 +end -- 2.39.5