]> https://gitweb.dealii.org/ - dealii.git/commitdiff
More sensible defaults, and document optional functions.
authorLuca Heltai <luca.heltai@sissa.it>
Wed, 6 Sep 2017 12:22:03 +0000 (14:22 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Fri, 15 Sep 2017 12:39:45 +0000 (14:39 +0200)
include/deal.II/sundials/ida_interface.h
source/sundials/ida_interface.cc
tests/quick_tests/sundials-ida.cc
tests/sundials/harmonic_oscillator_01.cc
tests/sundials/harmonic_oscillator_01.output
tests/sundials/harmonic_oscillator_01.prm

index f2264676565503fc81a663fe2242189c72983554..da37f7f587f79781ac16830271f7802e44bb6b47 100644 (file)
@@ -60,9 +60,15 @@ namespace SUNDIALS
    *  - residual;
    *  - setup_jacobian;
    *  - solve_jacobian_system;
-   *  - output_step;
+   *
+   * Optionally, also the following functions could be implemented. By default
+   * they do nothing, or are not required. If you call the constructor in a way
+   * that requires any of the not-implemented functions, an Assertion will be
+   * thrown.
    *  - solver_should_restart;
-   *  - differential_components.
+   *  - differential_components;
+   *  - output_step;
+   *  - get_local_tolerances;
    *
    * Citing from the SUNDIALS documentation:
    *
@@ -180,8 +186,8 @@ namespace SUNDIALS
                  const unsigned int &max_order = 5,
                  const double &output_period = .1,
                  const bool &ignore_algebraic_terms_for_errors = true,
-                 const std::string &ic_type = "use_y_diff",
-                 const std::string &reset_type = "use_y_dot",
+                 const std::string &ic_type = "none",
+                 const std::string &reset_type = "none",
                  const double &ic_alpha = .33,
                  const unsigned int &ic_max_iter = 5,
                  const unsigned int &max_non_linear_iterations = 10,
index 785b3d56b5df2166663d5f8dd926d1aca6f844db..cca559a2f2ccba28b0fc487e50dafa00344012d9 100644 (file)
@@ -408,7 +408,6 @@ namespace SUNDIALS
 
     copy(yy, solution);
     copy(yp, solution_dot);
-    copy(diff_id, differential_components());
 
     status = IDAInit(ida_mem, t_dae_residual<VectorType>, current_time, yy, yp);
     AssertIDA(status);
@@ -429,8 +428,18 @@ namespace SUNDIALS
     status = IDASetUserData(ida_mem, (void *) this);
     AssertIDA(status);
 
-    status = IDASetId(ida_mem, diff_id);
-    AssertIDA(status);
+    if (ic_type == "use_y_diff" || reset_type == "use_y_diff" || ignore_algebraic_terms_for_errors)
+      {
+        VectorType diff_comp_vector(solution);
+        diff_comp_vector = 0.0;
+        auto dc = differential_components();
+        for (auto i = dc.begin(); i != dc.end(); ++i)
+          diff_comp_vector[*i] = 1.0;
+
+        copy(diff_id, diff_comp_vector);
+        status = IDASetId(ida_mem, diff_id);
+        AssertIDA(status);
+      }
 
     status = IDASetSuppressAlg(ida_mem, ignore_algebraic_terms_for_errors);
     AssertIDA(status);
@@ -508,9 +517,9 @@ namespace SUNDIALS
   void IDAInterface<VectorType>::set_functions_to_trigger_an_assert()
   {
 
-    create_new_vector = []() ->std::shared_ptr<VectorType>
+    create_new_vector = []() ->std::unique_ptr<VectorType>
     {
-      std::shared_ptr<VectorType> p;
+      std::unique_ptr<VectorType> p;
       AssertThrow(false, ExcFunctionNotProvided("create_new_vector"));
       return p;
     };
@@ -548,23 +557,21 @@ namespace SUNDIALS
                      const VectorType &,
                      const unsigned int)
     {
-      AssertThrow(false, ExcFunctionNotProvided("output_step"));
+      return;
     };
 
     solver_should_restart = [](const double,
                                VectorType &,
                                VectorType &) ->bool
     {
-      bool ret=false;
-      AssertThrow(false, ExcFunctionNotProvided("solver_should_restart"));
-      return ret;
+      return false;
     };
 
-    differential_components = []() ->VectorType &
+    differential_components = []() ->IndexSet
     {
-      std::shared_ptr<VectorType> y;
+      IndexSet i;
       AssertThrow(false, ExcFunctionNotProvided("differential_components"));
-      return *y;
+      return i;
     };
 
     get_local_tolerances = []() ->VectorType &
index 5f83751c2fb05465ac18d0390db67f85f465a307..2f7c20c58a1339043ad9ed3eba5852cf657ab067 100644 (file)
@@ -68,9 +68,9 @@ public:
     diff[0] = 1.0;
     diff[1] = 1.0;
 
-    time_stepper.create_new_vector = [&] () -> std::shared_ptr<Vector<double> >
+    time_stepper.create_new_vector = [&] () -> std::unique_ptr<Vector<double> >
     {
-      return std::shared_ptr<Vector<double>>(new Vector<double>(2));
+      return std::unique_ptr<Vector<double>>(new Vector<double>(2));
     };
 
 
@@ -127,9 +127,9 @@ public:
     };
 
 
-    time_stepper.differential_components = [&]() -> VectorType&
+    time_stepper.differential_components = [&]() -> IndexSet
     {
-      return diff;
+      return complete_index_set(2);
     };
   }
 
index 4c9f1bec27fc23f56b3644d4e65de8ee8293e635..8bded1027b44ef6d5ce8d5bd7c495075b915e98f 100644 (file)
@@ -39,7 +39,7 @@
  * A = [ 0 , -1; k^2, 0 ]
  *
  * y_0  = 0, k
- * y_0' = 0, 0
+ * y_0' = k, 0
  *
  * The exact solution is
  *
@@ -58,19 +58,15 @@ public:
   HarmonicOscillator(double _kappa=1.0) :
     y(2),
     y_dot(2),
-    diff(2),
     J(2,2),
     A(2,2),
     Jinv(2,2),
     kappa(_kappa),
     out("output")
   {
-    diff[0] = 1.0;
-    diff[1] = 1.0;
-
-    time_stepper.create_new_vector = [&] () -> std::shared_ptr<Vector<double> >
+    time_stepper.create_new_vector = [&] () -> std::unique_ptr<Vector<double> >
     {
-      return std::shared_ptr<Vector<double>>(new Vector<double>(2));
+      return std::unique_ptr<Vector<double>>(new Vector<double>(2));
     };
 
 
@@ -116,7 +112,7 @@ public:
                                    const unsigned int step_number) -> int
     {
       out << t << " "
-      << sol << " " << sol_dot;
+      << sol[0] << " " << sol[1] << " " << sol_dot[0] << " " << sol_dot[1] << std::endl;
       return 0;
     };
 
@@ -126,24 +122,18 @@ public:
     {
       return false;
     };
-
-
-    time_stepper.differential_components = [&]() -> VectorType&
-    {
-      return diff;
-    };
   }
 
   void run()
   {
     y[1] = kappa;
+    y_dot[0] = kappa;
     time_stepper.solve_dae(y,y_dot);
   }
   SUNDIALS::IDAInterface<Vector<double> >  time_stepper;
 private:
   Vector<double> y;
   Vector<double> y_dot;
-  Vector<double> diff;
   FullMatrix<double> J;
   FullMatrix<double> A;
   FullMatrix<double> Jinv;
@@ -157,7 +147,7 @@ int main (int argc, char **argv)
 {
   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int);
 
-  HarmonicOscillator ode(2*numbers::PI);
+  HarmonicOscillator ode(1.0);
   ParameterHandler prm;
   ode.time_stepper.add_parameters(prm);
 
index 0f0f792c21e5bbe0a63799e1b977de7a8becccc0..7241ab086fc104c0f04af39b13176de1d09636b3 100644 (file)
@@ -1,42 +1,33 @@
-0 0.000e+00 6.283e+00 
- 6.283e+00 -2.481e-05 
-0.05 3.090e-01 5.976e+00 
- 5.976e+00 -1.220e+01 
-0.1 5.878e-01 5.083e+00 
- 5.083e+00 -2.321e+01 
-0.15 8.090e-01 3.693e+00 
- 3.693e+00 -3.194e+01 
-0.2 9.511e-01 1.942e+00 
- 1.941e+00 -3.755e+01 
-0.25 1.000e+00 -5.322e-05 
- -1.559e-04 -3.948e+01 
-0.3 9.510e-01 -1.942e+00 
- -1.942e+00 -3.755e+01 
-0.35 8.090e-01 -3.693e+00 
- -3.693e+00 -3.194e+01 
-0.4 5.878e-01 -5.083e+00 
- -5.083e+00 -2.320e+01 
-0.45 3.090e-01 -5.975e+00 
- -5.976e+00 -1.220e+01 
-0.5 -2.180e-05 -6.283e+00 
- -6.283e+00 1.846e-03 
-0.55 -3.090e-01 -5.975e+00 
- -5.975e+00 1.220e+01 
-0.6 -5.878e-01 -5.083e+00 
- -5.083e+00 2.321e+01 
-0.65 -8.090e-01 -3.693e+00 
- -3.693e+00 3.194e+01 
-0.7 -9.510e-01 -1.941e+00 
- -1.941e+00 3.754e+01 
-0.75 -9.999e-01 2.282e-04 
- 3.915e-04 3.948e+01 
-0.8 -9.510e-01 1.942e+00 
- 1.942e+00 3.754e+01 
-0.85 -8.089e-01 3.693e+00 
- 3.693e+00 3.193e+01 
-0.9 -5.877e-01 5.083e+00 
- 5.083e+00 2.320e+01 
-0.95 -3.090e-01 5.975e+00 
- 5.975e+00 1.220e+01 
-1 3.380e-05 6.283e+00 
- 6.283e+00 -1.802e-03 
+0 0 1 1 0
+0.2 0.198666 0.98005 0.980056 -0.198643
+0.4 0.389414 0.921051 0.921063 -0.3894
+0.6 0.56464 0.825325 0.825346 -0.564654
+0.8 0.717354 0.696694 0.696713 -0.717368
+1 0.841469 0.540288 0.540295 -0.841478
+1.2 0.932034 0.362343 0.362338 -0.932036
+1.4 0.985441 0.169953 0.16995 -0.985443
+1.6 0.999562 -0.0292115 -0.0292147 -0.999563
+1.8 0.973833 -0.227212 -0.227216 -0.973833
+2 0.90928 -0.416153 -0.416157 -0.90928
+2.2 0.808478 -0.588504 -0.588507 -0.808476
+2.4 0.675444 -0.737392 -0.737395 -0.675442
+2.6 0.515482 -0.856883 -0.856886 -0.515479
+2.8 0.33497 -0.942212 -0.942214 -0.334967
+3 0.141104 -0.989978 -0.98998 -0.141101
+3.2 -0.058387 -0.998277 -0.998278 0.0583908
+3.4 -0.25555 -0.966777 -0.966777 0.255554
+3.6 -0.442525 -0.896735 -0.896735 0.442529
+3.8 -0.611857 -0.790944 -0.790942 0.611861
+4 -0.756797 -0.65362 -0.653617 0.7568
+4.2 -0.871565 -0.490238 -0.490235 0.871568
+4.4 -0.951586 -0.307312 -0.307309 0.951588
+4.6 -0.993671 -0.112135 -0.112131 0.993672
+4.8 -0.996141 0.0875123 0.0875162 0.996141
+5 -0.958897 0.28367 0.283674 0.958897
+5.2 -0.883426 0.468519 0.468523 0.883425
+5.4 -0.772735 0.634689 0.634693 0.772733
+5.6 -0.631238 0.775556 0.775559 0.631236
+5.8 -0.464576 0.885504 0.885506 0.464573
+6 -0.279392 0.960149 0.960151 0.279389
+6.2 -0.0830705 0.996516 0.996517 0.0830668
+6.3 0.0168302 0.99983 0.99983 -0.0168291
index f91eb0de80196aa0522794fd4058c433a8917903..81776b7c1935eb976549f661a2758d4856d473df 100644 (file)
@@ -1,16 +1,16 @@
 set Absolute error tolerance                      = 0.000001
-set Final time                                    = 1.000000
-set Ignore algebraic terms for error computations = true
+set Final time                                    = 6.3
+set Ignore algebraic terms for error computations = false
 set Initial condition Newton max iterations       = 5
 set Initial condition Newton parameter            = 0.330000
-set Initial condition type                        = use_y_diff
-set Initial condition type after restart          = use_y_dot
-set Initial step size                             = 0.000100
-set Initial time                                  = 0.000000
+set Initial condition type                        = none
+set Initial condition type after restart          = none
+set Initial step size                             = 0.1
+set Initial time                                  = 0
 set Maximum number of nonlinear iterations        = 10
 set Maximum order of BDF                          = 5
 set Minimum step size                             = 0.000001
-set Relative error tolerance                      = 0.000010
+set Relative error tolerance                      = 0.00001
 set Show output of time steps                     = true
-set Time interval between each output             = 0.05
+set Time interval between each output             = 0.2
 set Use local tolerances                          = false

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.