9#include "cantera/numerics/sundials_headers.h"
17#if CT_SUNDIALS_VERSION >= 60
18 return N_VNew_Serial(
static_cast<sd_size_t
>(N), context.get());
20 return N_VNew_Serial(
static_cast<sd_size_t
>(N));
47static int ida_rhs(realtype t, N_Vector y, N_Vector ydot, N_Vector r,
void* f_data)
50 return f->
evalDaeNoThrow(t, NV_DATA_S(y), NV_DATA_S(ydot), NV_DATA_S(r));
56static void ida_err(
int error_code,
const char* module,
57 const char* function,
char* msg,
void* eh_data)
71IdasIntegrator::~IdasIntegrator()
77 N_VDestroy_Serial(
m_y);
86 N_VDestroy_Serial(m_constraints);
92 return NV_Ith_S(
m_y, k);
97 return NV_DATA_S(
m_y);
110 for (
size_t i=0; i<n; i++) {
139 checkError(flag,
"setMaxOrder",
"IDASetMaxOrd");
148 int flag = IDASetMaxStep(
m_ida_mem, hmax);
149 checkError(flag,
"setMaxStepSize",
"IDASetMaxStep");
180 int flag = IDAGetNumSteps(
m_ida_mem, &val);
181 checkError(flag,
"solverStats",
"IDAGetNumSteps");
182 stats[
"steps"] = val;
184 stats[
"res_evals"] = val;
186 stats[
"lin_solve_setups"] = val;
188 stats[
"err_tests_fails"] = val;
191 IDAGetNumNonlinSolvIters(
m_ida_mem, &val);
192 stats[
"nonlinear_iters"] = val;
193 IDAGetNumNonlinSolvConvFails(
m_ida_mem, &val);
194 stats[
"nonlinear_conv_fails"] = val;
198void IdasIntegrator::setMaxNonlinIterations(
int n)
203 checkError(flag,
"setMaxNonlinIterations",
"IDASetMaxNonlinIters");
207void IdasIntegrator::setMaxNonlinConvFailures(
int n)
212 checkError(flag,
"setMaxNonlinConvFailures",
"IDASetMaxConvFails");
216void IdasIntegrator::includeAlgebraicInErrorTest(
bool yesno)
221 checkError(flag,
"inclAlgebraicInErrorTest",
"IDASetSuppressAlg");
235 N_VDestroy_Serial(
m_y);
241 N_VDestroy_Serial(
m_ydot);
249 "not enough absolute tolerance values specified.");
253 N_VDestroy_Serial(m_constraints);
267 #if CT_SUNDIALS_VERSION >= 60
273 throw CanteraError(
"IdasIntegrator::initialize",
"IDACreate failed.");
277 if (flag != IDA_SUCCESS) {
278 if (flag == IDA_MEM_FAIL) {
280 "Memory allocation failed.");
281 }
else if (flag == IDA_ILL_INPUT) {
283 "Illegal value for IDAInit input argument.");
285 throw CanteraError(
"IdasIntegrator::initialize",
"IDAInit failed.");
292 flag = IDASetId(
m_ida_mem, m_constraints);
297 checkError(flag,
"initialize",
"IDASVtolerances");
300 checkError(flag,
"initialize",
"IDASStolerances");
304 checkError(flag,
"initialize",
"IDASetUserData");
307 throw CanteraError(
"IdasIntegrator::initialize",
"Sensitivity analysis "
308 "for DAE systems is not fully implemented");
312 checkError(flag,
"initialize",
"IDASetSensParams");
317void IdasIntegrator::reinitialize(
double t0,
FuncEval& func)
327 checkError(result,
"reinitialize",
"IDAReInit");
334 sd_size_t N =
static_cast<sd_size_t
>(
m_neq);
335 SUNLinSolFree((SUNLinearSolver)
m_linsol);
337 #if CT_SUNDIALS_VERSION >= 60
342 #if CT_SUNDIALS_VERSION >= 60
347 #if CT_SUNDIALS_USE_LAPACK
348 #if CT_SUNDIALS_VERSION >= 60
355 #if CT_SUNDIALS_VERSION >= 60
362 #if CT_SUNDIALS_VERSION >= 60
369 }
else if (
m_type ==
"GMRES") {
370 #if CT_SUNDIALS_VERSION >= 60
373 #elif CT_SUNDIALS_VERSION >= 40
382 "unsupported linear solver flag '{}'",
m_type);
391 checkError(flag,
"applyOptions",
"IDASetMaxOrd");
401 checkError(flag,
"applyOptions",
"IDASetMaxNonlinIters");
405 checkError(flag,
"applyOptions",
"IDASetMaxConvFails");
409 checkError(flag,
"applyOptions",
"IDASetSuppressAlg");
413void IdasIntegrator::sensInit(
double t0,
FuncEval& func)
419 #if CT_SUNDIALS_VERSION >= 60
420 m_yS = N_VCloneVectorArray(
static_cast<int>(
m_np), y);
422 m_yS = N_VCloneVectorArray_Serial(
static_cast<int>(
m_np), y);
424 for (
size_t n = 0; n <
m_np; n++) {
425 N_VConst(0.0,
m_yS[n]);
427 N_VDestroy_Serial(y);
429 #if CT_SUNDIALS_VERSION >= 60
430 m_ySdot = N_VCloneVectorArray(
static_cast<int>(
m_np), ydot);
432 m_ySdot = N_VCloneVectorArray_Serial(
static_cast<int>(
m_np), ydot);
434 for (
size_t n = 0; n <
m_np; n++) {
438 int flag = IDASensInit(
m_ida_mem,
static_cast<sd_size_t
>(
m_np),
442 vector<double> atol(
m_np);
443 for (
size_t n = 0; n <
m_np; n++) {
449 checkError(flag,
"sensInit",
"IDASensSStolerances");
456 }
else if (tout <
m_time) {
458 "Cannot integrate backwards in time.\n"
459 "Requested time {} < current time {}",
466 "Maximum number of timesteps ({}) taken without reaching output "
467 "time ({}).\nCurrent integrator time: {}",
471 if (flag != IDA_SUCCESS) {
473 if (!f_errs.empty()) {
474 f_errs =
"Exceptions caught during RHS evaluation:\n" + f_errs;
477 "IDA error encountered. Error code: {}\n{}\n"
479 "Components with largest weighted error estimates:\n{}",
492 if (flag != IDA_SUCCESS) {
494 if (!f_errs.empty()) {
495 f_errs =
"Exceptions caught during RHS evaluation:\n" + f_errs;
498 "IDA error encountered. Error code: {}\n{}\n"
500 "Components with largest weighted error estimates:\n{}",
508double IdasIntegrator::sensitivity(
size_t k,
size_t p)
516 checkError(flag,
"sensitivity",
"IDAGetSens");
521 throw CanteraError(
"IdasIntegrator::sensitivity",
522 "sensitivity: k out of range ({})", k);
525 throw CanteraError(
"IdasIntegrator::sensitivity",
526 "sensitivity: p out of range ({})", p);
528 return NV_Ith_S(
m_yS[p],k);
538 vector<tuple<double, double, size_t>> weightedErrors;
539 for (
size_t i=0; i<
m_neq; i++) {
540 double err = NV_Ith_S(errs, i) * NV_Ith_S(errw, i);
541 weightedErrors.emplace_back(-abs(err), err, i);
546 N = std::min(N,
static_cast<int>(
m_neq));
547 sort(weightedErrors.begin(), weightedErrors.end());
548 fmt::memory_buffer s;
549 for (
int i=0; i<N; i++) {
550 fmt_append(s,
"{}: {}\n", get<2>(weightedErrors[i]), get<1>(weightedErrors[i]));
556 const string& idaMethod)
const
558 if (flag == IDA_SUCCESS) {
560 }
else if (flag == IDA_MEM_NULL) {
562 "IDAS integrator is not initialized");
564 const char* flagname = IDAGetReturnFlagName(flag);
566 "{} returned error code {} ({}):\n{}",
575 throw CanteraError(
"IdasIntegrator::setMethod",
"unknown method");
Header file for class IdasIntegrator.
A map of string keys to values whose type can vary at runtime.
Base class for exceptions thrown by Cantera classes.
Virtual base class for ODE/DAE right-hand-side function evaluators.
virtual void getStateDae(double *y, double *ydot)
Fill in the vectors y and ydot with the current state of the system.
vector< double > m_paramScales
Scaling factors for each sensitivity parameter.
int evalDaeNoThrow(double t, double *y, double *ydot, double *residual)
Evaluate the right-hand side using return code to indicate status.
void clearErrors()
Clear any previously-stored suppressed errors.
virtual size_t neq() const =0
Number of equations.
virtual void getConstraints(double *constraints)
Given a vector of length neq(), mark which variables should be considered algebraic constraints.
virtual size_t nparams() const
Number of sensitivity parameters.
string getErrors() const
Return a string containing the text of any suppressed errors.
vector< double > m_sens_params
Values for the problem parameters for which sensitivities are computed This is the array which is per...
Wrapper for Sundials IDAS solver.
double step(double tout) override
Integrate the system of equations.
void setMaxStepSize(double hmax) override
Set the maximum step size.
double * solution() override
The current value of the solution of the system of equations.
double m_t0
The start time for the integrator.
N_Vector m_y
The current system state.
void * m_linsol
Sundials linear solver object.
int m_maxNonlinConvFails
Maximum number of nonlinear convergence failures.
double m_init_step
Initial IDA step size.
void checkError(long flag, const string &ctMethod, const string &idaMethod) const
Check whether an IDAS method indicated an error.
double m_abstolsens
Scalar absolute tolerance for sensitivities.
size_t m_nabs
! Number of variables for which absolute tolerances were provided
void setMaxOrder(int n) override
Set the maximum integration order that will be used.
double m_time
The current integrator time, corresponding to m_y.
int m_maxNonlinIters
Maximum number of nonlinear solver iterations at one solution.
bool m_sens_ok
Indicates whether the sensitivities stored in m_yS and m_ySdot have been updated for the current inte...
int maxSteps() override
Returns the maximum number of time-steps the integrator can take before reaching the next output time...
SundialsContext m_sundials_ctx
SUNContext object for Sundials>=6.0.
FuncEval * m_func
Object implementing the DAE residual function .
double m_abstols
Scalar absolute tolerance.
double m_hmax
Maximum time step size. Zero means infinity.
int m_itol
Flag indicating whether scalar (IDA_SS) or vector (IDA_SV) absolute tolerances are being used.
int m_maxErrTestFails
Maximum number of error test failures in attempting one step.
int m_maxord
Maximum order allowed for the BDF method.
void setSensitivityTolerances(double reltol, double abstol) override
Set the sensitivity error tolerances.
void applyOptions()
Applies user-specified options to the underlying IDAS solver.
void integrate(double tout) override
Integrate the system of equations.
void setTolerances(double reltol, size_t n, double *abstol) override
Set error tolerances.
void setLinearSolverType(const string &linearSolverType) override
Set the linear solver type.
double m_reltol
Relative tolerance for all state variables.
string m_error_message
Error message information provide by IDAS.
double m_reltolsens
Scalar relative tolerance for sensitivities.
N_Vector * m_yS
Sensitivities of y, size m_np by m_neq.
void setMaxErrTestFails(int n) override
Set the maximum permissible number of error test failures.
void setMaxSteps(int nmax) override
Set the maximum number of time-steps the integrator can take before reaching the next output time.
size_t m_neq
Number of equations/variables in the system.
AnyMap solverStats() const override
Get solver stats from integrator.
string m_type
The linear solver type.
size_t m_np
Number of sensitivity parameters.
void * m_linsol_matrix
matrix used by Sundials
N_Vector m_ydot
The time derivatives of the system state.
N_Vector * m_ySdot
Sensitivities of ydot, size m_np by m_neq.
double m_tInteg
The latest time reached by the integrator. May be greater than m_time.
void initialize(double t0, FuncEval &func) override
Initialize the integrator for a new problem.
int m_maxsteps
Maximum number of internal steps taken in a call to integrate()
string getErrorInfo(int N)
Returns a string listing the weighted error estimates associated with each solution component.
IdasIntegrator()
Constructor.
void setMethod(MethodType t) override
Set the solution method.
void * m_ida_mem
Pointer to the IDA memory for the problem.
bool m_setSuppressAlg
If true, the algebraic variables don't contribute to error tolerances.
N_Vector m_abstol
Absolute tolerances for each state variable.
virtual string linearSolverType() const
Return the integrator problem type.
virtual int lastOrder() const
Order used during the last solution step.
A wrapper for managing a SUNContext object, need for Sundials >= 6.0.
void fmt_append(fmt::memory_buffer &b, Args... args)
Versions 6.2.0 and 6.2.1 of fmtlib do not include this define before they include windows....
Namespace for the Cantera kernel.
static void ida_err(int error_code, const char *module, const char *function, char *msg, void *eh_data)
Function called by IDA when an error is encountered instead of writing to stdout.
MethodType
Specifies the method used to integrate the system of equations.
@ BDF_Method
Backward Differentiation.
static int ida_rhs(realtype t, N_Vector y, N_Vector ydot, N_Vector r, void *f_data)
Function called by IDA to evaluate the residual, given y and ydot.
Contains declarations for string manipulation functions within Cantera.