]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Brusselator test with jacobian.
authorLuca Heltai <luca.heltai@sissa.it>
Wed, 27 Sep 2017 12:33:48 +0000 (14:33 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Wed, 27 Sep 2017 12:36:18 +0000 (14:36 +0200)
tests/sundials/harmonic_oscillator_05.cc [new file with mode: 0644]
tests/sundials/harmonic_oscillator_05.output [new file with mode: 0644]
tests/sundials/harmonic_oscillator_05.prm [new file with mode: 0644]
tests/sundials/harmonic_oscillator_06.cc [new file with mode: 0644]
tests/sundials/harmonic_oscillator_06.output [new file with mode: 0644]
tests/sundials/harmonic_oscillator_06.prm [new file with mode: 0644]

diff --git a/tests/sundials/harmonic_oscillator_05.cc b/tests/sundials/harmonic_oscillator_05.cc
new file mode 100644 (file)
index 0000000..27d2f81
--- /dev/null
@@ -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 <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;
+}
diff --git a/tests/sundials/harmonic_oscillator_05.output b/tests/sundials/harmonic_oscillator_05.output
new file mode 100644 (file)
index 0000000..507044c
--- /dev/null
@@ -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 (file)
index 0000000..55a91a0
--- /dev/null
@@ -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 (file)
index 0000000..c58e56c
--- /dev/null
@@ -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 <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;
+}
diff --git a/tests/sundials/harmonic_oscillator_06.output b/tests/sundials/harmonic_oscillator_06.output
new file mode 100644 (file)
index 0000000..507044c
--- /dev/null
@@ -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 (file)
index 0000000..55a91a0
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.