# include <deal.II/lac/petsc_vector.h>
# endif
-# include <deal.II/sundials/copy.h>
+# include <deal.II/sundials/n_vector.h>
# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
# include <arkode/arkode_impl.h>
Assert(user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(user_data);
- GrowingVectorMemory<VectorType> mem;
-
- typename VectorMemory<VectorType>::Pointer src_yy(mem);
- solver.reinit_vector(*src_yy);
-
- typename VectorMemory<VectorType>::Pointer dst_yp(mem);
- solver.reinit_vector(*dst_yp);
-
- copy(*src_yy, yy);
- int err = solver.explicit_function(tt, *src_yy, *dst_yp);
+ auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
+ auto *dst_yp = internal::unwrap_nvector<VectorType>(yp);
- copy(yp, *dst_yp);
-
- return err;
+ return solver.explicit_function(tt, *src_yy, *dst_yp);
}
Assert(user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(user_data);
- GrowingVectorMemory<VectorType> mem;
- typename VectorMemory<VectorType>::Pointer src_yy(mem);
- solver.reinit_vector(*src_yy);
+ auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
+ auto *dst_yp = internal::unwrap_nvector<VectorType>(yp);
- typename VectorMemory<VectorType>::Pointer dst_yp(mem);
- solver.reinit_vector(*dst_yp);
-
- copy(*src_yy, yy);
-
- int err = solver.implicit_function(tt, *src_yy, *dst_yp);
-
- copy(yp, *dst_yp);
-
- return err;
+ return solver.implicit_function(tt, *src_yy, *dst_yp);
}
Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
- GrowingVectorMemory<VectorType> mem;
-
- typename VectorMemory<VectorType>::Pointer src_ypred(mem);
- solver.reinit_vector(*src_ypred);
-
- typename VectorMemory<VectorType>::Pointer src_fpred(mem);
- solver.reinit_vector(*src_fpred);
-
- copy(*src_ypred, ypred);
- copy(*src_fpred, fpred);
+ auto *src_ypred = internal::unwrap_nvector_const<VectorType>(ypred);
+ auto *src_fpred = internal::unwrap_nvector_const<VectorType>(fpred);
// avoid reinterpret_cast
bool jcurPtr_tmp = false;
int err = solver.setup_jacobian(convfail,
Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
- GrowingVectorMemory<VectorType> mem;
-
- typename VectorMemory<VectorType>::Pointer src(mem);
- solver.reinit_vector(*src);
- typename VectorMemory<VectorType>::Pointer src_ycur(mem);
- solver.reinit_vector(*src_ycur);
+ auto *dst = internal::unwrap_nvector<VectorType>(b);
+ auto *src_ycur = internal::unwrap_nvector_const<VectorType>(ycur);
+ auto *src_fcur = internal::unwrap_nvector_const<VectorType>(fcur);
- typename VectorMemory<VectorType>::Pointer src_fcur(mem);
- solver.reinit_vector(*src_fcur);
-
- typename VectorMemory<VectorType>::Pointer dst(mem);
- solver.reinit_vector(*dst);
-
- copy(*src, b);
- copy(*src_ycur, ycur);
- copy(*src_fcur, fcur);
+ // make a temporary copy to work on in the user call
+ VectorType src = *dst;
int err = solver.solve_jacobian_system(arkode_mem->ark_tn,
arkode_mem->ark_gamma,
*src_ycur,
*src_fcur,
- *src,
+ src,
*dst);
- copy(b, *dst);
+
return err;
}
Assert(arkode_mem->ark_user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
- GrowingVectorMemory<VectorType> mem;
-
- typename VectorMemory<VectorType>::Pointer src(mem);
- solver.reinit_vector(*src);
- typename VectorMemory<VectorType>::Pointer dst(mem);
- solver.reinit_vector(*dst);
+ auto *dst = internal::unwrap_nvector<VectorType>(b);
- copy(*src, b);
+ // make a temporary copy to work on in the user call
+ VectorType src = *dst;
- int err = solver.solve_mass_system(*src, *dst);
- copy(b, *dst);
+ int err = solver.solve_mass_system(src, *dst);
return err;
}
*static_cast<ARKode<VectorType> *>(user_data);
GrowingVectorMemory<VectorType> mem;
- typename VectorMemory<VectorType>::Pointer src_v(mem);
- solver.reinit_vector(*src_v);
- typename VectorMemory<VectorType>::Pointer dst_Jv(mem);
- solver.reinit_vector(*dst_Jv);
- typename VectorMemory<VectorType>::Pointer src_y(mem);
- solver.reinit_vector(*src_y);
- typename VectorMemory<VectorType>::Pointer src_fy(mem);
- solver.reinit_vector(*src_fy);
- copy(*src_v, v);
- copy(*src_y, y);
- copy(*src_fy, fy);
-
- int err =
- solver.jacobian_times_vector(*src_v, *dst_Jv, t, *src_y, *src_fy);
- copy(Jv, *dst_Jv);
- return err;
+ auto *src_v = internal::unwrap_nvector_const<VectorType>(v);
+ auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
+ auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
+
+ auto *dst_Jv = internal::unwrap_nvector<VectorType>(Jv);
+
+ return solver.jacobian_times_vector(*src_v, *dst_Jv, t, *src_y, *src_fy);
}
Assert(user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(user_data);
- GrowingVectorMemory<VectorType> mem;
- typename VectorMemory<VectorType>::Pointer src_y(mem);
- solver.reinit_vector(*src_y);
- typename VectorMemory<VectorType>::Pointer src_fy(mem);
- solver.reinit_vector(*src_fy);
- copy(*src_y, y);
- copy(*src_fy, fy);
- int err = solver.jacobian_times_setup(t, *src_y, *src_fy);
- return err;
+
+ auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
+ auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
+
+ return solver.jacobian_times_setup(t, *src_y, *src_fy);
}
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(user_data);
- GrowingVectorMemory<VectorType> mem;
- typename VectorMemory<VectorType>::Pointer src_y(mem);
- solver.reinit_vector(*src_y);
- typename VectorMemory<VectorType>::Pointer src_fy(mem);
- solver.reinit_vector(*src_fy);
- typename VectorMemory<VectorType>::Pointer src_r(mem);
- solver.reinit_vector(*src_r);
- typename VectorMemory<VectorType>::Pointer dst_z(mem);
- solver.reinit_vector(*dst_z);
- copy(*src_y, y);
- copy(*src_fy, fy);
- copy(*src_r, r);
-
- int status = solver.jacobian_preconditioner_solve(
+ auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
+ auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
+ auto *src_r = internal::unwrap_nvector_const<VectorType>(r);
+
+ auto *dst_z = internal::unwrap_nvector<VectorType>(z);
+
+ return solver.jacobian_preconditioner_solve(
t, *src_y, *src_fy, *src_r, *dst_z, gamma, delta, lr);
- copy(z, *dst_z);
- return status;
}
Assert(user_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(user_data);
- GrowingVectorMemory<VectorType> mem;
- typename VectorMemory<VectorType>::Pointer src_y(mem);
- solver.reinit_vector(*src_y);
- typename VectorMemory<VectorType>::Pointer src_fy(mem);
- solver.reinit_vector(*src_fy);
- copy(*src_y, y);
- copy(*src_fy, fy);
+
+ auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
+ auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
+
return solver.jacobian_preconditioner_setup(
t, *src_y, *src_fy, jok, *jcurPtr, gamma);
}
Assert(mtimes_data != nullptr, ExcInternalError());
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(mtimes_data);
- GrowingVectorMemory<VectorType> mem;
- typename VectorMemory<VectorType>::Pointer src_v(mem);
- solver.reinit_vector(*src_v);
- typename VectorMemory<VectorType>::Pointer dst_Mv(mem);
- solver.reinit_vector(*dst_Mv);
- copy(*src_v, v);
- int err = solver.mass_times_vector(t, *src_v, *dst_Mv);
- copy(Mv, *dst_Mv);
+ auto *src_v = internal::unwrap_nvector_const<VectorType>(v);
+ auto *dst_Mv = internal::unwrap_nvector<VectorType>(Mv);
- return err;
+ return solver.mass_times_vector(t, *src_v, *dst_Mv);
}
ARKode<VectorType> &solver =
*static_cast<ARKode<VectorType> *>(user_data);
- GrowingVectorMemory<VectorType> mem;
- typename VectorMemory<VectorType>::Pointer src_r(mem);
- solver.reinit_vector(*src_r);
- typename VectorMemory<VectorType>::Pointer dst_z(mem);
- solver.reinit_vector(*dst_z);
- copy(*src_r, r);
-
- int status =
- solver.mass_preconditioner_solve(t, *src_r, *dst_z, delta, lr);
- copy(z, *dst_z);
- return status;
+ auto *src_r = internal::unwrap_nvector_const<VectorType>(r);
+ auto *dst_z = internal::unwrap_nvector<VectorType>(z);
+
+ return solver.mass_preconditioner_solve(t, *src_r, *dst_z, delta, lr);
}
auto content = access_content<VectorType>(LS);
ARKode<VectorType> &solver = *(content->solver);
- GrowingVectorMemory<VectorType> mem;
- typename VectorMemory<VectorType>::Pointer rhs(mem);
- solver.reinit_vector(*rhs);
- typename VectorMemory<VectorType>::Pointer dst(mem);
- solver.reinit_vector(*dst);
+ auto *src_b = internal::unwrap_nvector_const<VectorType>(b);
+ auto *dst_x = internal::unwrap_nvector<VectorType>(x);
- copy(*rhs, b);
-
- SundialsOperator<VectorType> op(solver,
+ SundialsOperator<VectorType> op(solver,
content->A_data,
content->a_times_fn);
+
SundialsPreconditioner<VectorType> preconditioner(
solver, content->P_data, content->preconditioner_solve, tol);
- int err = content->lsolve(op, preconditioner, *dst, *rhs, tol);
- copy(x, *dst);
- return err;
+ return content->lsolve(op, preconditioner, *dst_x, *src_b, tol);
}
const MPI_Comm & mpi_comm)
: data(data)
, arkode_mem(nullptr)
- , yy(nullptr)
- , abs_tolls(nullptr)
, communicator(is_serial_vector<VectorType>::value ?
MPI_COMM_SELF :
Utilities::MPI::duplicate_communicator(mpi_comm))
unsigned int
ARKode<VectorType>::solve_ode(VectorType &solution)
{
- const unsigned int system_size = solution.size();
-
- // The solution is stored in solution. Here we take only a view of it.
-# ifdef DEAL_II_WITH_MPI
- if (is_serial_vector<VectorType>::value == false)
- {
- const IndexSet is = solution.locally_owned_elements();
- const std::size_t local_system_size = is.n_elements();
-
- yy = N_VNew_Parallel(communicator, local_system_size, system_size);
-
- abs_tolls =
- N_VNew_Parallel(communicator, local_system_size, system_size);
- }
- else
-# endif
- {
- Assert(is_serial_vector<VectorType>::value,
- ExcInternalError(
- "Trying to use a serial code with a parallel vector."));
- yy = N_VNew_Serial(system_size);
- abs_tolls = N_VNew_Serial(system_size);
- }
-
DiscreteTime time(data.initial_time,
data.final_time,
data.initial_step_size);
(void)status;
AssertARKode(status);
- copy(solution, yy);
-
time.set_next_step_size(actual_next_time - time.get_current_time());
time.advance_time();
time.get_step_number());
}
- // Free the vectors which are no longer used.
-# ifdef DEAL_II_WITH_MPI
- if (is_serial_vector<VectorType>::value == false)
- {
- N_VDestroy_Parallel(yy);
- N_VDestroy_Parallel(abs_tolls);
- }
- else
-# endif
- {
- N_VDestroy_Serial(yy);
- N_VDestroy_Serial(abs_tolls);
- }
-
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,
- const VectorType &solution)
+ ARKode<VectorType>::reset(const double current_time,
+ const double current_time_step,
+ VectorType & solution)
{
- unsigned int system_size;
-
if (arkode_mem)
ARKodeFree(&arkode_mem);
arkode_mem = ARKodeCreate();
- // Free the vectors which are no longer used.
- if (yy)
- {
- free_vector(yy);
- free_vector(abs_tolls);
- }
-
int status;
(void)status;
- system_size = solution.size();
-# ifdef DEAL_II_WITH_MPI
- if (is_serial_vector<VectorType>::value == false)
- {
- const IndexSet is = solution.locally_owned_elements();
- const std::size_t local_system_size = is.n_elements();
-
- yy = N_VNew_Parallel(communicator, local_system_size, system_size);
-
- abs_tolls =
- N_VNew_Parallel(communicator, local_system_size, system_size);
- }
- else
-# endif
- {
- yy = N_VNew_Serial(system_size);
- abs_tolls = N_VNew_Serial(system_size);
- }
-
- copy(yy, solution);
Assert(explicit_function || implicit_function,
ExcFunctionNotProvided("explicit_function || implicit_function"));
+ // 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);
+
status = ARKodeInit(
arkode_mem,
explicit_function ? &t_arkode_explicit_function<VectorType> : nullptr,
if (get_local_tolerances)
{
- copy(abs_tolls, get_local_tolerances());
+ abs_tolls = make_nvector_view(get_local_tolerances());
status =
ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tolls);
AssertARKode(status);
void
ARKode<VectorType>::reset(const double current_time,
const double current_time_step,
- const VectorType &solution)
+ VectorType &solution)
{
if (arkode_mem)
ARKStepFree(&arkode_mem);
- // Free the vectors which are no longer used.
- if (yy)
- {
- free_vector(yy);
- free_vector(abs_tolls);
- }
-
int status;
(void)status;
- yy = create_vector(solution);
- abs_tolls = create_vector(solution);
- copy(yy, solution);
+
+ // 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);
Assert(explicit_function || implicit_function,
ExcFunctionNotProvided("explicit_function || implicit_function"));
current_time,
yy);
+ Assert(arkode_mem != nullptr, ExcInternalError());
+
if (get_local_tolerances)
{
- copy(abs_tolls, get_local_tolerances());
+ abs_tolls = internal::make_nvector_view(get_local_tolerances());
status =
ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tolls);
AssertARKode(status);
}
else
{
- N_Vector y_template = create_vector(solution);
+ auto y_template = internal::make_nvector_view(solution);
SUNNonlinearSolver fixed_point_solver =
SUNNonlinSol_FixedPoint(y_template,
data.anderson_acceleration_subspace);
- free_vector(y_template);
status = ARKStepSetNonlinearSolver(arkode_mem, fixed_point_solver);
AssertARKode(status);
}
- template <typename VectorType>
- N_Vector
- ARKode<VectorType>::create_vector(const VectorType &template_vector) const
- {
- N_Vector new_vector;
- auto system_size = template_vector.size();
-
-# ifdef DEAL_II_WITH_MPI
- if (is_serial_vector<VectorType>::value == false)
- {
- new_vector =
- N_VNew_Parallel(communicator,
- template_vector.locally_owned_elements().n_elements(),
- system_size);
- }
- else
-# endif
- {
- new_vector = N_VNew_Serial(system_size);
- }
- return new_vector;
- }
-
- template <typename VectorType>
- void
- ARKode<VectorType>::free_vector(N_Vector vector) const
- {
-# ifdef DEAL_II_WITH_MPI
- if (is_serial_vector<VectorType>::value == false)
- N_VDestroy_Parallel(vector);
- else
-# endif
- N_VDestroy_Serial(vector);
- }
-
-
-
template <typename VectorType>
void *
ARKode<VectorType>::get_arkode_memory() const
SundialsOperator<VectorType>::vmult(VectorType & dst,
const VectorType &src) const
{
- N_Vector sun_dst = solver.create_vector(dst);
- N_Vector sun_src = solver.create_vector(src);
- copy(sun_src, src);
- int status = a_times_fn(A_data, sun_src, sun_dst);
+ auto sun_dst = internal::make_nvector_view(dst);
+ auto sun_src = internal::make_nvector_view(src);
+ int status = a_times_fn(A_data, sun_src, sun_dst);
(void)status;
AssertARKode(status);
-
- copy(dst, sun_dst);
-
- solver.free_vector(sun_dst);
- solver.free_vector(sun_src);
}
return;
}
- N_Vector sun_dst = solver.create_vector(dst);
- N_Vector sun_src = solver.create_vector(src);
- copy(sun_src, src);
+ auto sun_dst = internal::make_nvector_view(dst);
+ auto sun_src = internal::make_nvector_view(src);
// for custom preconditioners no distinction between left and right
// preconditioning is made
int status =
p_solve_fn(P_data, sun_src, sun_dst, tol, 0 /*precondition_type*/);
(void)status;
AssertARKode(status);
-
- copy(dst, sun_dst);
-
- solver.free_vector(sun_dst);
- solver.free_vector(sun_src);
}
# endif