From 90c037e4501ad733aa6ad0baeebe80062ba79a37 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Wed, 6 Sep 2017 14:22:03 +0200 Subject: [PATCH] More sensible defaults, and document optional functions. --- include/deal.II/sundials/ida_interface.h | 14 ++-- source/sundials/ida_interface.cc | 31 ++++---- tests/quick_tests/sundials-ida.cc | 8 +-- tests/sundials/harmonic_oscillator_01.cc | 22 ++---- tests/sundials/harmonic_oscillator_01.output | 75 +++++++++----------- tests/sundials/harmonic_oscillator_01.prm | 16 ++--- 6 files changed, 80 insertions(+), 86 deletions(-) diff --git a/include/deal.II/sundials/ida_interface.h b/include/deal.II/sundials/ida_interface.h index f226467656..da37f7f587 100644 --- a/include/deal.II/sundials/ida_interface.h +++ b/include/deal.II/sundials/ida_interface.h @@ -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, diff --git a/source/sundials/ida_interface.cc b/source/sundials/ida_interface.cc index 785b3d56b5..cca559a2f2 100644 --- a/source/sundials/ida_interface.cc +++ b/source/sundials/ida_interface.cc @@ -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, 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::set_functions_to_trigger_an_assert() { - create_new_vector = []() ->std::shared_ptr + create_new_vector = []() ->std::unique_ptr { - std::shared_ptr p; + std::unique_ptr 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 y; + IndexSet i; AssertThrow(false, ExcFunctionNotProvided("differential_components")); - return *y; + return i; }; get_local_tolerances = []() ->VectorType & diff --git a/tests/quick_tests/sundials-ida.cc b/tests/quick_tests/sundials-ida.cc index 5f83751c2f..2f7c20c58a 100644 --- a/tests/quick_tests/sundials-ida.cc +++ b/tests/quick_tests/sundials-ida.cc @@ -68,9 +68,9 @@ public: diff[0] = 1.0; diff[1] = 1.0; - time_stepper.create_new_vector = [&] () -> std::shared_ptr > + time_stepper.create_new_vector = [&] () -> std::unique_ptr > { - return std::shared_ptr>(new Vector(2)); + return std::unique_ptr>(new Vector(2)); }; @@ -127,9 +127,9 @@ public: }; - time_stepper.differential_components = [&]() -> VectorType& + time_stepper.differential_components = [&]() -> IndexSet { - return diff; + return complete_index_set(2); }; } diff --git a/tests/sundials/harmonic_oscillator_01.cc b/tests/sundials/harmonic_oscillator_01.cc index 4c9f1bec27..8bded1027b 100644 --- a/tests/sundials/harmonic_oscillator_01.cc +++ b/tests/sundials/harmonic_oscillator_01.cc @@ -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 > + time_stepper.create_new_vector = [&] () -> std::unique_ptr > { - return std::shared_ptr>(new Vector(2)); + return std::unique_ptr>(new Vector(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 > time_stepper; private: Vector y; Vector y_dot; - Vector diff; FullMatrix J; FullMatrix A; FullMatrix 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); diff --git a/tests/sundials/harmonic_oscillator_01.output b/tests/sundials/harmonic_oscillator_01.output index 0f0f792c21..7241ab086f 100644 --- a/tests/sundials/harmonic_oscillator_01.output +++ b/tests/sundials/harmonic_oscillator_01.output @@ -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 diff --git a/tests/sundials/harmonic_oscillator_01.prm b/tests/sundials/harmonic_oscillator_01.prm index f91eb0de80..81776b7c19 100644 --- a/tests/sundials/harmonic_oscillator_01.prm +++ b/tests/sundials/harmonic_oscillator_01.prm @@ -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 -- 2.39.5