# include <sundials/sundials_linearsolver.h>
+# include <exception>
# include <functional>
# include <memory>
class LinearSolverWrapper
{
public:
- explicit LinearSolverWrapper(LinearSolveFunction<VectorType> lsolve
+ explicit LinearSolverWrapper(
+ const LinearSolveFunction<VectorType> &lsolve,
+ std::exception_ptr & pending_exception
# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
- ,
- SUNContext linsol_ctx
+ ,
+ SUNContext linsol_ctx
# endif
);
# endif
, mpi_communicator(mpi_comm)
, last_end_time(data.initial_time)
+ , pending_exception(nullptr)
{
set_functions_to_trigger_an_assert();
(void)status;
AssertARKode(status);
# endif
+
+ Assert(pending_exception == nullptr, ExcInternalError());
}
,
arkode_ctx
# endif
- );
+ ,
+ pending_exception);
sun_linear_solver = *linear_solver;
}
else
,
arkode_ctx
# endif
- );
+ ,
+ pending_exception);
sun_mass_linear_solver = *mass_solver;
}
else
# define AssertSundialsSolver(code) \
Assert(code >= 0, ExcSundialsSolverError(code))
+ namespace
+ {
+ /**
+ * A function that calls the function object given by its first argument
+ * with the set of arguments following at the end. If the call returns
+ * regularly, the current function returns zero to indicate success. If the
+ * call fails with an exception of type RecoverableUserCallbackError, then
+ * the current function returns 1 to indicate that the called function
+ * object thought the error it encountered is recoverable. If the call fails
+ * with any other exception, then the current function returns with an error
+ * code of -1. In each of the last two cases, the exception thrown by `f`
+ * is captured and `eptr` is set to the exception. In case of success,
+ * `eptr` is set to `nullptr`.
+ */
+ template <typename F, typename... Args>
+ int
+ call_and_possibly_capture_exception(const F & f,
+ std::exception_ptr &eptr,
+ Args &&...args)
+ {
+ // See whether there is already something in the exception pointer
+ // variable. This can only happen if we had previously had
+ // a recoverable exception, and the underlying library actually
+ // did recover successfully. In that case, we can abandon the
+ // exception previously thrown. If eptr contains anything other,
+ // then we really don't know how that could have happened, and
+ // should probably bail out:
+ if (eptr)
+ {
+ try
+ {
+ std::rethrow_exception(eptr);
+ }
+ catch (const RecoverableUserCallbackError &)
+ {
+ // ok, ignore, but reset the pointer
+ eptr = nullptr;
+ }
+ catch (...)
+ {
+ // uh oh:
+ AssertThrow(false, ExcInternalError());
+ }
+ }
+
+ // Call the function and if that succeeds, return zero:
+ try
+ {
+ f(std::forward<Args>(args)...);
+ eptr = nullptr;
+ return 0;
+ }
+ // If the call failed with a recoverable error, then
+ // ignore the exception for now (but store a pointer to it)
+ // and return a positive return value (+1). If the underlying
+ // implementation manages to recover
+ catch (const RecoverableUserCallbackError &)
+ {
+ eptr = std::current_exception();
+ return 1;
+ }
+ // For any other exception, capture the exception and
+ // return -1:
+ catch (const std::exception &)
+ {
+ eptr = std::current_exception();
+ return -1;
+ }
+ }
+ } // namespace
+
+
namespace internal
{
/**
template <typename VectorType>
struct LinearSolverContent
{
- LinearSolverContent()
+ LinearSolverContent(std::exception_ptr &pending_exception)
: a_times_fn(nullptr)
, preconditioner_setup(nullptr)
, preconditioner_solve(nullptr)
# endif
, P_data(nullptr)
, A_data(nullptr)
+ , pending_exception(pending_exception)
{}
ATimesFn a_times_fn;
void *P_data;
void *A_data;
+
+ /**
+ * A reference to a location where we can store exceptions, should they
+ * be thrown by a linear solver object.
+ */
+ std::exception_ptr &pending_exception;
};
} // namespace internal
N_Vector b,
realtype tol)
{
- auto content = access_content<VectorType>(LS);
+ LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
auto *src_b = unwrap_nvector_const<VectorType>(b);
auto *dst_x = unwrap_nvector<VectorType>(x);
# endif
tol);
- return content->lsolve(op, preconditioner, *dst_x, *src_b, tol);
+ return call_and_possibly_capture_exception(content->lsolve,
+ content->pending_exception,
+ op,
+ preconditioner,
+ *dst_x,
+ *src_b,
+ tol);
}
int
arkode_linsol_setup(SUNLinearSolver LS, SUNMatrix /*ignored*/)
{
- auto content = access_content<VectorType>(LS);
+ LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
+
if (content->preconditioner_setup)
return content->preconditioner_setup(content->P_data);
return 0;
int
arkode_linsol_set_a_times(SUNLinearSolver LS, void *A_data, ATimesFn ATimes)
{
- auto content = access_content<VectorType>(LS);
+ LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
+
content->A_data = A_data;
content->a_times_fn = ATimes;
return 0;
PSetupFn p_setup,
PSolveFn p_solve)
{
- auto content = access_content<VectorType>(LS);
+ LinearSolverContent<VectorType> *content = access_content<VectorType>(LS);
+
content->P_data = P_data;
content->preconditioner_setup = p_setup;
content->preconditioner_solve = p_solve;
template <typename VectorType>
internal::LinearSolverWrapper<VectorType>::LinearSolverWrapper(
- LinearSolveFunction<VectorType> lsolve
+ const LinearSolveFunction<VectorType> &lsolve,
+ std::exception_ptr & pending_exception
# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
,
SUNContext linsol_ctx
# endif
)
- : content(std::make_unique<LinearSolverContent<VectorType>>())
+ : content(
+ std::make_unique<LinearSolverContent<VectorType>>(pending_exception))
{
# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
sun_linear_solver = SUNLinSolNewEmpty(linsol_ctx);
double kappa = 1.0;
- ode.explicit_function =
- [&](double, const VectorType &y, VectorType &ydot) -> int {
+ ode.explicit_function = [&](double, const VectorType &y, VectorType &ydot) {
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;
- };
+ ode.output_step =
+ [&](const double t, const VectorType &sol, const unsigned int step_number) {
+ deallog << t << ' ' << sol[0] << ' ' << sol[1] << std::endl;
+ };
Vector<double> y(2);
y[0] = 0;
y[1] = kappa;
ode.solve_ode(y);
- return 0;
}
double kappa = 1.0;
- ode.implicit_function =
- [&](double, const VectorType &y, VectorType &ydot) -> int {
+ ode.implicit_function = [&](double, const VectorType &y, VectorType &ydot) {
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 {
- // limit the output to every 10th step and increase the precision to make
- // the test more robust
- if (step_number % 10 == 0)
- deallog << t << ' ' << std::setprecision(7) << sol[0] << ' ' << sol[1]
- << std::endl;
- return 0;
- };
+ ode.output_step =
+ [&](const double t, const VectorType &sol, const unsigned int step_number) {
+ // limit the output to every 10th step and increase the precision to make
+ // the test more robust
+ if (step_number % 10 == 0)
+ deallog << t << ' ' << std::setprecision(7) << sol[0] << ' ' << sol[1]
+ << std::endl;
+ };
Vector<double> y(2);
y[0] = 0;
y[1] = kappa;
ode.solve_ode(y);
- return 0;
}
// Parameters
double u0 = 3.9, v0 = 1.1, w0 = 2.8, a = 1.2, b = 2.5, eps = 1e-5;
- ode.implicit_function =
- [&](double, const VectorType &y, VectorType &ydot) -> int {
+ ode.implicit_function = [&](double, const VectorType &y, VectorType &ydot) {
ydot[0] = 0;
ydot[1] = 0;
ydot[2] = -y[2] / eps;
- return 0;
};
- ode.explicit_function =
- [&](double, const VectorType &y, VectorType &ydot) -> int {
+ ode.explicit_function = [&](double, const VectorType &y, VectorType &ydot) {
ydot[0] = a - (y[2] + 1) * y[0] + y[1] * y[0] * y[0];
ydot[1] = y[2] * y[0] - y[1] * y[0] * y[0];
ydot[2] = b / eps - y[2] * y[0];
- return 0;
};
- ode.output_step = [&](const double t,
- const VectorType & sol,
- const unsigned int step_number) -> int {
- // limit the output to every 10th step and increase the precision to make
- // the test more robust
- if (step_number % 10 == 0)
- deallog << t << ' ' << std::setprecision(10) << sol[0] << ' ' << sol[1]
- << ' ' << sol[2] << std::endl;
- return 0;
- };
+ ode.output_step =
+ [&](const double t, const VectorType &sol, const unsigned int step_number) {
+ // limit the output to every 10th step and increase the precision to make
+ // the test more robust
+ if (step_number % 10 == 0)
+ deallog << t << ' ' << std::setprecision(10) << sol[0] << ' ' << sol[1]
+ << ' ' << sol[2] << std::endl;
+ };
// This test, for reasons I don't fully understand, generates some output
// which varies between environments much more than the other ARKODE
// tests. Work around it by setting a fairly stringent maximum time step.
- ode.custom_setup = [&](void *arkode_mem) -> int {
+ ode.custom_setup = [&](void *arkode_mem) {
int ierr = ARKStepSetMinStep(arkode_mem, 1e-8);
AssertThrow(ierr == 0, ExcInternalError());
ierr = ARKStepSetMaxStep(arkode_mem, 1e-4);
AssertThrow(ierr == 0, ExcInternalError());
ierr = ARKStepSetMaxNumSteps(arkode_mem, 5000);
AssertThrow(ierr == 0, ExcInternalError());
- return 0;
};
VectorType y(3);
y[1] = v0;
y[2] = w0;
ode.solve_ode(y);
- return 0;
}
FullMatrix<double> J(3, 3);
J(2, 2) = -1.0 / eps;
- ode.implicit_function =
- [&](double, const VectorType &y, VectorType &ydot) -> int {
+ ode.implicit_function = [&](double, const VectorType &y, VectorType &ydot) {
ydot[0] = 0;
ydot[1] = 0;
ydot[2] = -y[2] / eps;
- return 0;
};
- ode.explicit_function =
- [&](double, const VectorType &y, VectorType &ydot) -> int {
+ ode.explicit_function = [&](double, const VectorType &y, VectorType &ydot) {
ydot[0] = a - (y[2] + 1) * y[0] + y[1] * y[0] * y[0];
ydot[1] = y[2] * y[0] - y[1] * y[0] * y[0];
ydot[2] = b / eps - y[2] * y[0];
- return 0;
};
ode.jacobian_times_vector = [&](const VectorType &src,
VectorType & dst,
double t,
const VectorType & /*y*/,
- const VectorType & /*fy*/) -> int {
+ const VectorType & /*fy*/) {
J.vmult(dst, src);
-
- return 0;
};
- ode.output_step = [&](const double t,
- const VectorType & sol,
- const unsigned int step_number) -> int {
- deallog << std::setprecision(16) << t << ' ' << sol[0] << ' ' << sol[1]
- << ' ' << sol[2] << std::endl;
- return 0;
- };
+ ode.output_step =
+ [&](const double t, const VectorType &sol, const unsigned int step_number) {
+ deallog << std::setprecision(16) << t << ' ' << sol[0] << ' ' << sol[1]
+ << ' ' << sol[2] << std::endl;
+ };
Vector<double> y(3);
y[0] = u0;
y[1] = v0;
y[2] = w0;
ode.solve_ode(y);
- return 0;
}
// Explicit jacobian.
FullMatrix<double> J(3, 3);
- ode.implicit_function =
- [&](double, const VectorType &y, VectorType &ydot) -> int {
+ ode.implicit_function = [&](double, const VectorType &y, VectorType &ydot) {
ydot[0] = 0;
ydot[1] = 0;
ydot[2] = -y[2] / eps;
- return 0;
};
- ode.explicit_function =
- [&](double, const VectorType &y, VectorType &ydot) -> int {
+ ode.explicit_function = [&](double, const VectorType &y, VectorType &ydot) {
ydot[0] = a - (y[2] + 1) * y[0] + y[1] * y[0] * y[0];
ydot[1] = y[2] * y[0] - y[1] * y[0] * y[0];
ydot[2] = b / eps - y[2] * y[0];
- return 0;
};
const double gamma,
const VectorType &,
const VectorType &,
- bool &j_is_current) -> int {
+ bool &j_is_current) {
J = 0;
J(0, 0) = 1;
J(1, 1) = 1;
J(2, 2) = 1 + gamma / eps;
J.gauss_jordan();
j_is_current = true;
- return 0;
};
ode.solve_jacobian_system = [&](const double t,
const VectorType &,
const VectorType &,
const VectorType &src,
- VectorType & dst) -> int {
- J.vmult(dst, src);
- return 0;
- };
+ VectorType & dst) { J.vmult(dst, src); };
- ode.output_step = [&](const double t,
- const VectorType & sol,
- const unsigned int step_number) -> int {
- deallog << t << ' ' << sol[0] << ' ' << sol[1] << ' ' << sol[2]
- << std::endl;
- return 0;
- };
+ ode.output_step =
+ [&](const double t, const VectorType &sol, const unsigned int step_number) {
+ deallog << t << ' ' << sol[0] << ' ' << sol[1] << ' ' << sol[2]
+ << std::endl;
+ };
Vector<double> y(3);
y[0] = u0;
y[1] = v0;
y[2] = w0;
ode.solve_ode(y);
- return 0;
}
// Explicit jacobian.
FullMatrix<double> J(3, 3);
- ode.implicit_function =
- [&](double, const VectorType &y, VectorType &ydot) -> int {
+ ode.implicit_function = [&](double, const VectorType &y, VectorType &ydot) {
ydot[0] = 0;
ydot[1] = 0;
ydot[2] = -y[2] / eps;
- return 0;
};
- ode.explicit_function =
- [&](double, const VectorType &y, VectorType &ydot) -> int {
+ ode.explicit_function = [&](double, const VectorType &y, VectorType &ydot) {
ydot[0] = a - (y[2] + 1) * y[0] + y[1] * y[0] * y[0];
ydot[1] = y[2] * y[0] - y[1] * y[0] * y[0];
ydot[2] = b / eps - y[2] * y[0];
- return 0;
};
ode.jacobian_times_setup =
- [&](realtype t, const VectorType &y, const VectorType &fy) -> int {
- J = 0;
- J(2, 2) = -1.0 / eps;
- return 0;
- };
+ [&](realtype t, const VectorType &y, const VectorType &fy) {
+ J = 0;
+ J(2, 2) = -1.0 / eps;
+ };
ode.jacobian_times_vector = [&](const VectorType &v,
VectorType & Jv,
double t,
const VectorType &y,
- const VectorType &fy) -> int {
- J.vmult(Jv, v);
- return 0;
- };
+ const VectorType &fy) { J.vmult(Jv, v); };
ode.solve_linearized_system =
[&](SUNDIALS::SundialsOperator<VectorType> &op,
SUNDIALS::SundialsPreconditioner<VectorType> &,
VectorType & x,
const VectorType &b,
- double tol) -> int {
- ReductionControl control;
- SolverCG<VectorType> solver_cg(control);
- solver_cg.solve(op, x, b, PreconditionIdentity());
- return 0;
- };
-
- ode.output_step = [&](const double t,
- const VectorType & sol,
- const unsigned int step_number) -> int {
- deallog << t << ' ' << sol[0] << ' ' << sol[1] << ' ' << sol[2]
- << std::endl;
- return 0;
- };
+ double tol) {
+ ReductionControl control;
+ SolverCG<VectorType> solver_cg(control);
+ solver_cg.solve(op, x, b, PreconditionIdentity());
+ };
+
+ ode.output_step =
+ [&](const double t, const VectorType &sol, const unsigned int step_number) {
+ deallog << t << ' ' << sol[0] << ' ' << sol[1] << ' ' << sol[2]
+ << std::endl;
+ };
// after 5.2.0 a special interpolation mode should be used for stiff problems
#if DEAL_II_SUNDIALS_VERSION_GTE(5, 2, 0)
y[1] = v0;
y[2] = w0;
ode.solve_ode(y);
- return 0;
}
// Explicit jacobian.
FullMatrix<double> J(3, 3);
- ode.implicit_function =
- [&](double, const VectorType &y, VectorType &ydot) -> int {
+ ode.implicit_function = [&](double, const VectorType &y, VectorType &ydot) {
ydot[0] = 0;
ydot[1] = 0;
ydot[2] = -y[2] / eps;
- return 0;
};
- ode.explicit_function =
- [&](double, const VectorType &y, VectorType &ydot) -> int {
+ ode.explicit_function = [&](double, const VectorType &y, VectorType &ydot) {
ydot[0] = a - (y[2] + 1) * y[0] + y[1] * y[0] * y[0];
ydot[1] = y[2] * y[0] - y[1] * y[0] * y[0];
ydot[2] = b / eps - y[2] * y[0];
- return 0;
};
ode.jacobian_times_setup =
- [&](realtype t, const VectorType &y, const VectorType &fy) -> int {
- J = 0;
- J(2, 2) = -1.0 / eps;
- return 0;
- };
+ [&](realtype t, const VectorType &y, const VectorType &fy) {
+ J = 0;
+ J(2, 2) = -1.0 / eps;
+ };
ode.jacobian_times_vector = [&](const VectorType &v,
VectorType & Jv,
double t,
const VectorType &y,
- const VectorType &fy) -> int {
- J.vmult(Jv, v);
- return 0;
- };
+ const VectorType &fy) { J.vmult(Jv, v); };
ode.solve_linearized_system =
[&](SUNDIALS::SundialsOperator<VectorType> & op,
SUNDIALS::SundialsPreconditioner<VectorType> &prec,
VectorType & x,
const VectorType & b,
- double tol) -> int {
- ReductionControl control;
- SolverCG<VectorType> solver_cg(control);
- solver_cg.solve(op, x, b, prec);
- return 0;
- };
+ double tol) {
+ ReductionControl control;
+ SolverCG<VectorType> solver_cg(control);
+ solver_cg.solve(op, x, b, prec);
+ };
ode.jacobian_preconditioner_setup = [&](double t,
const VectorType &y,
const VectorType &fy,
int jok,
int & jcur,
- double gamma) -> int {
+ double gamma) {
deallog << "jacobian_preconditioner_setup called\n";
- return 0;
};
ode.jacobian_preconditioner_solve = [&](double t,
VectorType & z,
double gamma,
double delta,
- int lr) -> int {
+ int lr) {
deallog << "jacobian_preconditioner_solve called\n";
z = r;
- return 0;
};
- ode.output_step = [&](const double t,
- const VectorType & sol,
- const unsigned int step_number) -> int {
- deallog << t << ' ' << sol[0] << ' ' << sol[1] << ' ' << sol[2]
- << std::endl;
- return 0;
- };
+ ode.output_step =
+ [&](const double t, const VectorType &sol, const unsigned int step_number) {
+ deallog << t << ' ' << sol[0] << ' ' << sol[1] << ' ' << sol[2]
+ << std::endl;
+ };
// after 5.2.0 a special interpolation mode should be used for stiff problems
#if DEAL_II_SUNDIALS_VERSION_GTE(5, 2, 0)
y[1] = v0;
y[2] = w0;
ode.solve_ode(y);
- return 0;
}
M(0, 0) = M(1, 1) = 2.0 / 3;
M(1, 0) = M(0, 1) = 1.0 / 3;
- ode.implicit_function =
- [&](double, const VectorType &y, VectorType &ydot) -> int {
+ ode.implicit_function = [&](double, const VectorType &y, VectorType &ydot) {
K.vmult(ydot, y);
- return 0;
};
- ode.explicit_function =
- [&](double, const VectorType &y, VectorType &ydot) -> int {
+ ode.explicit_function = [&](double, const VectorType &y, VectorType &ydot) {
ydot[0] = 1;
ydot[1] = 2;
- return 0;
};
ode.jacobian_times_vector = [&](const VectorType &v,
VectorType & Jv,
double t,
const VectorType &y,
- const VectorType &fy) -> int {
- K.vmult(Jv, v);
- return 0;
- };
+ const VectorType &fy) { K.vmult(Jv, v); };
const auto solve_function =
[&](SUNDIALS::SundialsOperator<VectorType> & op,
SUNDIALS::SundialsPreconditioner<VectorType> &prec,
VectorType & x,
const VectorType & b,
- double tol) -> int {
- ReductionControl control;
- SolverCG<VectorType> solver_cg(control);
- solver_cg.solve(op, x, b, prec);
- return 0;
- };
+ double tol) {
+ ReductionControl control;
+ SolverCG<VectorType> solver_cg(control);
+ solver_cg.solve(op, x, b, prec);
+ };
ode.solve_linearized_system = solve_function;
FullMatrix<double> M_inv(2, 2);
- ode.mass_preconditioner_solve = [&](double t,
- const VectorType &r,
- VectorType & z,
- double gamma,
- int lr) -> int {
- LogStream::Prefix prefix("mass_preconditioner_solve");
- deallog << "applied" << std::endl;
- M_inv.vmult(z, r);
- return 0;
- };
+ ode.mass_preconditioner_solve =
+ [&](double t, const VectorType &r, VectorType &z, double gamma, int lr) {
+ LogStream::Prefix prefix("mass_preconditioner_solve");
+ deallog << "applied" << std::endl;
+ M_inv.vmult(z, r);
+ };
- ode.mass_preconditioner_setup = [&](double t) -> int {
+ ode.mass_preconditioner_setup = [&](double t) {
LogStream::Prefix prefix("mass_preconditioner_setup");
deallog << "applied" << std::endl;
M_inv.invert(M);
- return 0;
};
- ode.mass_times_vector =
- [&](const double t, const VectorType &v, VectorType &Mv) -> int {
- M.vmult(Mv, v);
- return 0;
- };
+ ode.mass_times_vector = [&](const double t,
+ const VectorType &v,
+ VectorType & Mv) { M.vmult(Mv, v); };
- 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;
- };
+ ode.output_step =
+ [&](const double t, const VectorType &sol, const unsigned int step_number) {
+ deallog << t << ' ' << sol[0] << ' ' << sol[1] << std::endl;
+ };
Vector<double> y(2);
y[0] = 1;
y[1] = 0;
ode.solve_ode(y);
- return 0;
}
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 {
+ ode->explicit_function = [&](double, const VectorType &y, VectorType &ydot) {
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 << std::setprecision(16) << t << ' ' << sol[0] << ' ' << sol[1]
- << std::endl;
- return 0;
- };
+ ode->output_step =
+ [&](const double t, const VectorType &sol, const unsigned int step_number) {
+ deallog << std::setprecision(16) << t << ' ' << sol[0] << ' ' << sol[1]
+ << std::endl;
+ };
return ode;
}
ode->reset(0.0, 0.01, y0);
ode->solve_ode_incrementally(y, 2.0);
}
-
- return 0;
}
projects / dealii.git / commitdiff