# ifdef DEAL_II_WITH_MPI
# include <nvector/nvector_parallel.h>
# endif
+# include <deal.II/base/discrete_time.h>
+
# include <deal.II/sundials/n_vector.h>
# include <deal.II/sundials/sunlinsol_wrapper.h>
/**
* Integrate the initial value problem. This function returns the final
* number of computed steps.
+ *
+ * @param solution On input, this vector contains the initial condition. On
+ * output, it contains the solution at the final time.
*/
unsigned int
solve_ode(VectorType &solution);
+ /**
+ * Integrate the initial value problem. Compared to the function above, this
+ * function allows to specify an @p intermediate_time for the next solution.
+ * Repeated calls of this function must use monotonously increasing values
+ * for @p intermediate_time. The last solution state is saved internally
+ * along with the @p intermediate_time and will be reused as initial
+ * condition for the next call.
+ *
+ * Users may find this function useful when integrating ARKode into an outer
+ * time loop of their own, especially when output_step() is too restrictive.
+ *
+ * @note @p intermediate_time may be larger than AdditionalData::final_time,
+ * which is ignored by this function.
+ *
+ * @param solution The final solution. If the solver restarts, either
+ * because it is the first ever solve or the flag @p reset_solver is
+ * set, the vector is also used as initial condition.
+ * @param intermediate_time The time for the incremental solution step. Must
+ * be greater than the last time that was used in a previous call to this
+ * function.
+ * @param reset_solver Optional flag to recreate all internal objects which
+ * may be desirable for spatial adaptivity methods. If set to `true`,
+ * reset() is called before solving the ODE, which sets @p solution as
+ * initial condition. This will *not* reset the stored time from previous
+ * calls to this function.
+ */
+ unsigned int
+ solve_ode_incrementally(VectorType & solution,
+ const double intermediate_time,
+ const bool reset_solver = false);
+
/**
* Clear internal memory and start with clean objects. This function is
* called when the simulation starts and when the user returns true to a
* a different function to solver_should_restart() that performs all mesh
* changes, transfers the solution to the new mesh, and returns true.
*
- * @param[in] t The new starting time
- * @param[in] h The new starting time step
- * @param[in,out] y The new initial solution
+ * @param t The new starting time
+ * @param h The new starting time step
+ * @param y The new initial solution
*/
void
- reset(const double t, const double h, VectorType &y);
+ reset(const double t, const double h, const VectorType &y);
/**
* Provides user access to the internally used ARKODE memory.
<< "Please provide an implementation for the function \""
<< arg1 << "\"");
+ /**
+ * Internal routine to call ARKode repeatedly.
+ */
+ int
+ do_evolve_time(VectorType & solution,
+ dealii::DiscreteTime &time,
+ const bool do_reset);
+
# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
/**
/**
* Set up the solver and preconditioner for a non-identity mass matrix in
* the ARKODE memory object based on the user-specified functions.
+ * @param solution The solution vector which is used as a template to create
+ * new vectors.
*/
void
- setup_mass_solver();
+ setup_mass_solver(const VectorType &solution);
# endif
*/
void *arkode_mem;
- /**
- * ARKode solution vector.
- */
- internal::NVectorView<VectorType> yy;
-
- /**
- * ARKode absolute tolerances vector.
- */
- internal::NVectorView<VectorType> abs_tolls;
-
/**
* MPI communicator. SUNDIALS solver runs happily in
* parallel. Note that if the library is compiled without MPI
MPI_Comm communicator;
/**
- * Memory pool of vectors.
+ * The final time in the last call to solve_ode().
*/
- GrowingVectorMemory<VectorType> mem;
+ double last_end_time;
# if DEAL_II_SUNDIALS_VERSION_GTE(4, 0, 0)
std::unique_ptr<internal::LinearSolverWrapper<VectorType>> linear_solver;
Assert(user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(user_data);
- GrowingVectorMemory<VectorType> mem;
auto *src_v = internal::unwrap_nvector_const<VectorType>(v);
auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
, communicator(is_serial_vector<VectorType>::value ?
MPI_COMM_SELF :
Utilities::MPI::duplicate_communicator(mpi_comm))
+ , last_end_time(data.initial_time)
{
set_functions_to_trigger_an_assert();
}
DiscreteTime time(data.initial_time,
data.final_time,
data.initial_step_size);
- reset(time.get_current_time(), time.get_next_step_size(), solution);
- if (output_step)
- output_step(time.get_current_time(), solution, time.get_step_number());
+ return do_evolve_time(solution, time, /* force_solver_restart = */ true);
+ }
+
- while (time.is_at_end() == false)
+
+ template <typename VectorType>
+ unsigned int
+ ARKode<VectorType>::solve_ode_incrementally(VectorType & solution,
+ const double intermediate_time,
+ const bool reset_solver)
+ {
+ AssertThrow(
+ intermediate_time > last_end_time,
+ dealii::ExcMessage(
+ "The requested intermediate time is smaller than the last requested "
+ "intermediate time."));
+
+ const bool do_reset = reset_solver || arkode_mem == nullptr;
+ DiscreteTime time(last_end_time, intermediate_time, data.initial_step_size);
+ return do_evolve_time(solution, time, do_reset);
+ }
+
+
+
+ template <typename VectorType>
+ int
+ ARKode<VectorType>::do_evolve_time(VectorType & solution,
+ DiscreteTime &time,
+ const bool do_reset)
+ {
+ if (do_reset)
+ {
+ reset(time.get_current_time(), time.get_next_step_size(), solution);
+ if (output_step)
+ output_step(time.get_current_time(),
+ solution,
+ time.get_step_number());
+ }
+
+ auto solution_nvector = internal::make_nvector_view(solution);
+
+ while (!time.is_at_end())
{
time.set_desired_next_step_size(data.output_period);
double actual_next_time;
# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
- const auto status = SundialsARKode(
- arkode_mem, time.get_next_time(), yy, &actual_next_time, ARK_NORMAL);
+ const auto status = SundialsARKode(arkode_mem,
+ time.get_next_time(),
+ solution_nvector,
+ &actual_next_time,
+ ARK_NORMAL);
# else
- const auto status = ARKStepEvolve(
- arkode_mem, time.get_next_time(), yy, &actual_next_time, ARK_NORMAL);
+ const auto status = ARKStepEvolve(arkode_mem,
+ time.get_next_time(),
+ solution_nvector,
+ &actual_next_time,
+ ARK_NORMAL);
# endif
(void)status;
AssertARKode(status);
solution,
time.get_step_number());
}
-
+ last_end_time = time.get_current_time();
return time.get_step_number();
}
+
+
# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
template <typename VectorType>
void
- ARKode<VectorType>::reset(const double current_time,
- const double current_time_step,
- VectorType & solution)
+ ARKode<VectorType>::reset(const double current_time,
+ const double current_time_step,
+ const VectorType &solution)
{
+ last_end_time = current_time;
if (arkode_mem)
ARKodeFree(&arkode_mem);
// just a view on the memory in solution, all write operations on yy by
// ARKODE will automatically be mirrored to solution
- yy = internal::make_nvector_view(solution);
+ auto initial_condition_nvector = internal::make_nvector_view(solution);
status = ARKodeInit(
arkode_mem,
explicit_function ? &t_arkode_explicit_function<VectorType> : nullptr,
implicit_function ? &t_arkode_implicit_function<VectorType> : nullptr,
current_time,
- yy);
+ initial_condition_nvector);
AssertARKode(status);
if (get_local_tolerances)
{
- abs_tolls = make_nvector_view(get_local_tolerances());
+ const auto abs_tols = make_nvector_view(get_local_tolerances());
status =
- ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tolls);
+ ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tols);
AssertARKode(status);
}
else
void
ARKode<VectorType>::reset(const double current_time,
const double current_time_step,
- VectorType &solution)
+ const VectorType &solution)
{
+ last_end_time = current_time;
if (arkode_mem)
ARKStepFree(&arkode_mem);
// just a view on the memory in solution, all write operations on yy by
// ARKODE will automatically be mirrored to solution
- yy = internal::make_nvector_view(solution);
+ auto initial_condition_nvector = internal::make_nvector_view(solution);
Assert(explicit_function || implicit_function,
ExcFunctionNotProvided("explicit_function || implicit_function"));
explicit_function ? &t_arkode_explicit_function<VectorType> : nullptr,
implicit_function ? &t_arkode_implicit_function<VectorType> : nullptr,
current_time,
- yy);
+ initial_condition_nvector);
Assert(arkode_mem != nullptr, ExcInternalError());
if (get_local_tolerances)
{
- abs_tolls = internal::make_nvector_view(get_local_tolerances());
+ const auto abs_tols =
+ internal::make_nvector_view(get_local_tolerances());
status =
- ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tolls);
+ ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tols);
AssertARKode(status);
}
else
setup_system_solver(solution);
- setup_mass_solver();
+ setup_mass_solver(solution);
status = ARKStepSetInitStep(arkode_mem, current_time_step);
AssertARKode(status);
{
// use default solver from SUNDIALS
// TODO give user options
+ auto y_template = internal::make_nvector_view(solution);
sun_linear_solver =
- SUNLinSol_SPGMR(yy,
+ SUNLinSol_SPGMR(y_template,
PREC_NONE,
0 /*krylov subvectors, 0 uses default*/);
}
template <typename VectorType>
void
- ARKode<VectorType>::setup_mass_solver()
+ ARKode<VectorType>::setup_mass_solver(const VectorType &solution)
{
int status;
(void)status;
}
else
{
+ auto y_template = internal::make_nvector_view(solution);
sun_mass_linear_solver =
- SUNLinSol_SPGMR(yy,
+ SUNLinSol_SPGMR(y_template,
PREC_NONE,
0 /*krylov subvectors, 0 uses default*/);
}
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II 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.md at
+// the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+
+#include <deal.II/base/parameter_handler.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/sundials/arkode.h>
+
+#include "../tests.h"
+
+// Test repeated calls to solve_ode with harmonic oscillator problem
+
+using VectorType = Vector<double>;
+
+double kappa = 1.0;
+
+std::unique_ptr<SUNDIALS::ARKode<VectorType>>
+create_solver()
+{
+ ParameterHandler prm;
+ SUNDIALS::ARKode<VectorType>::AdditionalData data;
+ data.add_parameters(prm);
+
+ std::ifstream ifile(SOURCE_DIR "/arkode_repeated_solve.prm");
+ prm.parse_input(ifile);
+
+ auto ode = std::make_unique<SUNDIALS::ARKode<VectorType>>(data);
+
+ // will yield analytic solution y[0] = sin(kappa*t); y[1] = kappa*cos(kappa*t)
+ ode->explicit_function =
+ [&](double, const VectorType &y, VectorType &ydot) -> int {
+ ydot[0] = y[1];
+ ydot[1] = -kappa * kappa * y[0];
+ return 0;
+ };
+
+ ode->output_step = [&](const double t,
+ const VectorType & sol,
+ const unsigned int step_number) -> int {
+ deallog << t << " " << sol[0] << " " << sol[1] << std::endl;
+ return 0;
+ };
+ return ode;
+}
+
+int
+main(int argc, char **argv)
+{
+ initlog();
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(
+ argc, argv, numbers::invalid_unsigned_int);
+
+ auto ode = create_solver();
+ {
+ Vector<double> y(2);
+ y[0] = 0;
+ y[1] = kappa;
+
+ // solve the ode in one go
+ ode->solve_ode(y);
+
+ deallog << "First solve done" << std::endl;
+
+ // solve again -> will call reset automatically
+ y[0] = 0;
+ y[1] = kappa;
+ ode->solve_ode(y);
+ deallog << "Second solve done" << std::endl;
+
+ // Solve for a time greater than the originally requested final time
+ ode->solve_ode_incrementally(y, 3.0);
+ deallog << "Solve further than final time done" << std::endl;
+ }
+
+ // solve in two steps with new object
+ ode = create_solver();
+
+ {
+ Vector<double> y(2);
+ y[0] = 0;
+ y[1] = kappa;
+
+ ode->solve_ode_incrementally(y, 1.0);
+ deallog << "Split solve - part 1 done" << std::endl;
+ }
+
+ {
+ // N.B. initial conditions are not necessary as the internal state is reused
+ Vector<double> y(2);
+
+ ode->solve_ode_incrementally(y, 2.0);
+ deallog << "Split solve - part 2 done" << std::endl;
+ }
+
+ {
+ Vector<double> y(2);
+ y[0] = 0;
+ y[1] = kappa;
+ // solving with single-argument call reruns from start
+ ode->solve_ode(y);
+ deallog << "Rerun with single argument call done" << std::endl;
+ }
+
+ {
+ Vector<double> y0(2), y(2);
+ y0[0] = 0;
+ y0[1] = kappa;
+ // manually reset solver
+ ode->reset(0.0, 0.01, y0);
+ ode->solve_ode_incrementally(y, 2.0);
+ }
+
+ return 0;
+}
--- /dev/null
+set Final time = 2.0
+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 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
+
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.100000 0.0998334 0.995004
+DEAL::0.200000 0.198669 0.980067
+DEAL::0.300000 0.295521 0.955336
+DEAL::0.400000 0.389420 0.921060
+DEAL::0.500000 0.479426 0.877583
+DEAL::0.600000 0.564643 0.825335
+DEAL::0.700000 0.644219 0.764841
+DEAL::0.800000 0.717356 0.696707
+DEAL::0.900000 0.783330 0.621606
+DEAL::1.00000 0.841478 0.540292
+DEAL::1.10000 0.891214 0.453587
+DEAL::1.20000 0.932042 0.362356
+DEAL::1.30000 0.963560 0.267498
+DEAL::1.40000 0.985452 0.169958
+DEAL::1.50000 0.997498 0.0707265
+DEAL::1.60000 0.999576 -0.0292043
+DEAL::1.70000 0.991667 -0.128845
+DEAL::1.80000 0.973849 -0.227205
+DEAL::1.90000 0.946302 -0.323292
+DEAL::2.00000 0.909300 -0.416148
+DEAL::First solve done
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.100000 0.0998334 0.995004
+DEAL::0.200000 0.198669 0.980067
+DEAL::0.300000 0.295521 0.955336
+DEAL::0.400000 0.389420 0.921060
+DEAL::0.500000 0.479426 0.877583
+DEAL::0.600000 0.564643 0.825335
+DEAL::0.700000 0.644219 0.764841
+DEAL::0.800000 0.717356 0.696707
+DEAL::0.900000 0.783330 0.621606
+DEAL::1.00000 0.841478 0.540292
+DEAL::1.10000 0.891214 0.453587
+DEAL::1.20000 0.932042 0.362356
+DEAL::1.30000 0.963560 0.267498
+DEAL::1.40000 0.985452 0.169958
+DEAL::1.50000 0.997498 0.0707265
+DEAL::1.60000 0.999576 -0.0292043
+DEAL::1.70000 0.991667 -0.128845
+DEAL::1.80000 0.973849 -0.227205
+DEAL::1.90000 0.946302 -0.323292
+DEAL::2.00000 0.909300 -0.416148
+DEAL::Second solve done
+DEAL::2.10000 0.863210 -0.504849
+DEAL::2.20000 0.808497 -0.588504
+DEAL::2.30000 0.745707 -0.666278
+DEAL::2.40000 0.675463 -0.737397
+DEAL::2.50000 0.598472 -0.801147
+DEAL::2.60000 0.515502 -0.856891
+DEAL::2.70000 0.427379 -0.904075
+DEAL::2.80000 0.334987 -0.942226
+DEAL::2.90000 0.239250 -0.970961
+DEAL::3.00000 0.141118 -0.989995
+DEAL::Solve further than final time done
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.100000 0.0998334 0.995004
+DEAL::0.200000 0.198669 0.980067
+DEAL::0.300000 0.295521 0.955336
+DEAL::0.400000 0.389420 0.921060
+DEAL::0.500000 0.479426 0.877583
+DEAL::0.600000 0.564643 0.825335
+DEAL::0.700000 0.644219 0.764841
+DEAL::0.800000 0.717356 0.696707
+DEAL::0.900000 0.783330 0.621606
+DEAL::1.00000 0.841478 0.540292
+DEAL::Split solve - part 1 done
+DEAL::1.10000 0.891214 0.453587
+DEAL::1.20000 0.932042 0.362356
+DEAL::1.30000 0.963560 0.267498
+DEAL::1.40000 0.985452 0.169958
+DEAL::1.50000 0.997498 0.0707265
+DEAL::1.60000 0.999576 -0.0292043
+DEAL::1.70000 0.991667 -0.128845
+DEAL::1.80000 0.973849 -0.227205
+DEAL::1.90000 0.946302 -0.323292
+DEAL::2.00000 0.909300 -0.416148
+DEAL::Split solve - part 2 done
+DEAL::0.00000 0.00000 1.00000
+DEAL::0.100000 0.0998334 0.995004
+DEAL::0.200000 0.198669 0.980067
+DEAL::0.300000 0.295521 0.955336
+DEAL::0.400000 0.389420 0.921060
+DEAL::0.500000 0.479426 0.877583
+DEAL::0.600000 0.564643 0.825335
+DEAL::0.700000 0.644219 0.764841
+DEAL::0.800000 0.717356 0.696707
+DEAL::0.900000 0.783330 0.621606
+DEAL::1.00000 0.841478 0.540292
+DEAL::1.10000 0.891214 0.453587
+DEAL::1.20000 0.932042 0.362356
+DEAL::1.30000 0.963560 0.267498
+DEAL::1.40000 0.985452 0.169958
+DEAL::1.50000 0.997498 0.0707265
+DEAL::1.60000 0.999576 -0.0292043
+DEAL::1.70000 0.991667 -0.128845
+DEAL::1.80000 0.973849 -0.227205
+DEAL::1.90000 0.946302 -0.323292
+DEAL::2.00000 0.909300 -0.416148
+DEAL::Rerun with single argument call done
+DEAL::0.100000 0.0998334 0.995004
+DEAL::0.200000 0.198669 0.980067
+DEAL::0.300000 0.295521 0.955336
+DEAL::0.400000 0.389420 0.921060
+DEAL::0.500000 0.479426 0.877583
+DEAL::0.600000 0.564643 0.825335
+DEAL::0.700000 0.644219 0.764841
+DEAL::0.800000 0.717356 0.696707
+DEAL::0.900000 0.783330 0.621606
+DEAL::1.00000 0.841478 0.540292
+DEAL::1.10000 0.891214 0.453587
+DEAL::1.20000 0.932042 0.362356
+DEAL::1.30000 0.963560 0.267498
+DEAL::1.40000 0.985452 0.169958
+DEAL::1.50000 0.997498 0.0707265
+DEAL::1.60000 0.999576 -0.0292043
+DEAL::1.70000 0.991667 -0.128845
+DEAL::1.80000 0.973849 -0.227205
+DEAL::1.90000 0.946302 -0.323292
+DEAL::2.00000 0.909300 -0.416148
projects / dealii.git / commitdiff