* - 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:
*
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,
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);
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);
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;
};
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 &
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));
};
};
- time_stepper.differential_components = [&]() -> VectorType&
+ time_stepper.differential_components = [&]() -> IndexSet
{
- return diff;
+ return complete_index_set(2);
};
}
* A = [ 0 , -1; k^2, 0 ]
*
* y_0 = 0, k
- * y_0' = 0, 0
+ * y_0' = k, 0
*
* The exact solution is
*
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));
};
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;
};
{
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;
{
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);
-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
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