--- /dev/null
+//-----------------------------------------------------------
+//
+// 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 <deal.II/sundials/arkode.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+
+// 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<double> VectorType;
+
+ ParameterHandler prm;
+ SUNDIALS::ARKode<VectorType>::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<VectorType> 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<double> 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<double> y(3);
+ y[0] = u0;
+ y[1] = v0;
+ y[2] = w0;
+ ode.solve_ode(y);
+ return 0;
+}
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+//-----------------------------------------------------------
+//
+// 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 <deal.II/sundials/arkode.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+
+// 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<double> VectorType;
+
+ ParameterHandler prm;
+ SUNDIALS::ARKode<VectorType>::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<VectorType> 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<double> 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<double> y(3);
+ y[0] = u0;
+ y[1] = v0;
+ y[2] = w0;
+ ode.solve_ode(y);
+ return 0;
+}
--- /dev/null
+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
--- /dev/null
+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