Skip to content

Commit

Permalink
update cvode version to 6.2.0
Browse files Browse the repository at this point in the history
  • Loading branch information
@munan
munan committed May 6, 2024
1 parent 6c52580 commit e8fa1da
Show file tree
Hide file tree
Showing 10 changed files with 793 additions and 787 deletions.
94 changes: 55 additions & 39 deletions cvodeDense.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,31 +1,39 @@
#include "cvodeDense.h"

CvodeDense::CvodeDense(Ode &ode,
CvodeDense::CvodeDense(Ode &ode,
const double reltol, const double *abstol,
const bool userJac)
:ode_(ode),
kDimen(ode.Dimen()),
reltol_(reltol),
abstol_(NULL),
cvode_mem_(NULL),
sunctx_(NULL),
dense_matrix_(NULL),
dense_ls_(NULL),
t_solve_(0.),
t_solve_eq_(0.),
nst_last_(0),
nst_eq_(0) {
int flag;
// Create the SUNDIALS context
flag = SUNContext_Create(NULL, &sunctx_);
CheckFlag(&flag, "SUNContext_Create", 1);
//create vector y
ode_.y_ = N_VNew_Serial(kDimen, sunctx_);
CheckFlag(static_cast<void *>(ode_.y_), "N_VNew_Serial", 0);

/* -----------Initialize absolute value vector----------*/
abstol_ = N_VNew_Serial(kDimen);
CheckFlag((void *)abstol_, "N_VNew_Serial", 0);
abstol_ = N_VNew_Serial(kDimen, sunctx_);
CheckFlag(static_cast<void *>(abstol_), "N_VNew_Serial", 0);
for (int i=0; i<kDimen; i++) {
NV_Ith_S(abstol_, i) = abstol[i];
}

/* ----------------Intialize solver---------------------*/
/* Call CVodeCreate to create the solver memory and specify the
/* Call CVodeCreate to create the solver memory and specify the
* Backward Differentiation Formula and the use of a Newton iteration */
cvode_mem_ = CVodeCreate(CV_BDF);
cvode_mem_ = CVodeCreate(CV_BDF, sunctx_);
CheckFlag((void *)cvode_mem_, "CVodeCreate", 0);

/* Set the user data pointer to ode*/
Expand All @@ -35,7 +43,7 @@ CvodeDense::CvodeDense(Ode &ode,
/* Call CVodeInit to initialize the integrator memory and specify the
* user's right hand side function in y'=f(t,y), the inital time T0, and
* the initial dependent variable vector y. */
flag = CVodeInit(cvode_mem_, ode_.wrapRHS,
flag = CVodeInit(cvode_mem_, ode_.wrapRHS,
ode_.t_, ode_.y_);
CheckFlag(&flag, "CVodeInit", 1);

Expand All @@ -49,21 +57,21 @@ CvodeDense::CvodeDense(Ode &ode,
//CheckFlag(&flag, "CVDense", 1);

/* Create dense SUNMatrix for use in linear solves */
dense_matrix_ = SUNDenseMatrix(kDimen, kDimen);
dense_matrix_ = SUNDenseMatrix(kDimen, kDimen, sunctx_);
CheckFlag((void *)dense_matrix_, "SUNDenseMatrix", 0);

/* Create dense SUNLinearSolver object for use by CVode */
dense_ls_ = SUNDenseLinearSolver(ode_.y_, dense_matrix_);
dense_ls_ = SUNLinSol_Dense(ode_.y_, dense_matrix_, sunctx_);
CheckFlag((void *)dense_ls_, "SUNDenseLinearSolver", 0);

/* Call CVDlsSetLinearSolver to attach the matrix and linear solver to CVode */
flag = CVDlsSetLinearSolver(cvode_mem_, dense_ls_, dense_matrix_);
CheckFlag(&flag, "CVDlsSetLinearSolver", 1);
/* Call CVodeSetLinearSolver to attach the matrix and linear solver to CVode */
flag = CVodeSetLinearSolver(cvode_mem_, dense_ls_, dense_matrix_);
CheckFlag(&flag, "CVodeSetLinearSolver", 1);

/* Set the Jacobian routine to Jac (user-supplied) */
if (userJac) {
flag = CVDlsSetJacFn(cvode_mem_, ode_.wrapJac);
CheckFlag(&flag, "CVDlsSetJacFn", 1);
flag = CVodeSetJacFn(cvode_mem_, ode_.wrapJac);
CheckFlag(&flag, "CVodeSetJacFn", 1);
}
}

Expand Down Expand Up @@ -97,10 +105,18 @@ void CvodeDense::SetInitStep(const double hinit) {


CvodeDense::~CvodeDense() {
/* destroy vector y*/
N_VDestroy(ode_.y_);
/*Free abstol_ vector*/
N_VDestroy_Serial(abstol_);
N_VDestroy(abstol_);
/* Free integrator memory */
CVodeFree(&cvode_mem_);
// Free linear solver memory
SUNLinSolFree(dense_ls_);
// Free matrix memory
SUNMatDestroy(dense_matrix_);
// Free SUNDIALS context
SUNContext_Free(&sunctx_);
}

void CvodeDense::Solve(const double tfinal) {
Expand All @@ -111,16 +127,16 @@ void CvodeDense::Solve(const double tfinal) {
CV_NORMAL);
CheckFlag(&flag, "CVode", 3);
t = clock() - t;
t_solve_ = ((float)t)/CLOCKS_PER_SEC; /*record the run time*/
t_solve_ = ((float)t)/CLOCKS_PER_SEC; /*record the run time*/
nst_last_ = GetNstep() - nst;
/*floor small negative abundances to avoid numerical problems*/
/*TODO: adjust species abundances to conserve the elements*/
for (int i=0; i<kDimen; i++) {
if ( NV_Ith_S(ode_.y_, i) < 0. ) {
printf("Warning: correcting negative aboundances of species ");
//printf("%d from %0.4e to abstol_[i] = %0.4e.\n", i,
//printf("%d from %0.4e to abstol_[i] = %0.4e.\n", i,
// NV_Ith_S(ode_.y_, i), NV_Ith_S(abstol_, i));
printf("%d from %0.4e to zero.\n", i,
printf("%d from %0.4e to zero.\n", i,
NV_Ith_S(ode_.y_, i));
//NV_Ith_S(ode_.y_, i) = NV_Ith_S(abstol_, i);
NV_Ith_S(ode_.y_, i) = 0.;
Expand All @@ -139,7 +155,7 @@ void CvodeDense::SolveEq(const double tolfac, const double tmax,
}

/* Guess the time toward equalibrium.
* To make a guess, we compute dx/dt at the initial conditions and at a point
* To make a guess, we compute dx/dt at the initial conditions and at a point
* slightly perturbed from them, and then use the most restrictive timestep we find.
* We also require the minmum time of 10yr=3.16e8s and maximum time of
* tmax/10 */
Expand All @@ -151,13 +167,13 @@ void CvodeDense::SolveEq(const double tolfac, const double tmax,
double dt1 = tevol_max;
double dt2 = tevol_max;
double tevol = 0.;
N_Vector ydot = N_VNew_Serial(kDimen);
N_Vector ypert = N_VNew_Serial(kDimen);
N_Vector ydot = N_VNew_Serial(kDimen, sunctx_);
N_Vector ypert = N_VNew_Serial(kDimen, sunctx_);
/* get the minimum component of dt1, restrict to be smaller than tevol_max*/
ode_.RHS(ode_.t_, ode_.y_, ydot);
for (int i=0; i<kDimen; i++) {
dt1 = std::min( dt1, fabs(
( NV_Ith_S(ode_.y_, i) + NV_Ith_S(abstol_, i) )
dt1 = std::min( dt1, fabs(
( NV_Ith_S(ode_.y_, i) + NV_Ith_S(abstol_, i) )
/ ( NV_Ith_S(ydot, i) + small__ ) )
);
}
Expand All @@ -169,8 +185,8 @@ void CvodeDense::SolveEq(const double tolfac, const double tmax,
* than tevol_max*/
ode_.RHS(ode_.t_+dt1, ypert, ydot);
for (int i=0; i<kDimen; i++) {
dt2 = std::min( dt2, fabs(
( NV_Ith_S(ypert, i) + NV_Ith_S(abstol_, i) )
dt2 = std::min( dt2, fabs(
( NV_Ith_S(ypert, i) + NV_Ith_S(abstol_, i) )
/ ( NV_Ith_S(ydot, i) + small__ ) )
);
}
Expand All @@ -183,11 +199,11 @@ void CvodeDense::SolveEq(const double tolfac, const double tmax,
}

/*destroy created vectors*/
N_VDestroy_Serial(ydot);
N_VDestroy_Serial(ypert);
N_VDestroy(ydot);
N_VDestroy(ypert);

/*find equalibrium.*/
N_Vector yprev = N_VNew_Serial(kDimen);
N_Vector yprev = N_VNew_Serial(kDimen, sunctx_);
double err[kDimen];
double maxerr;
clock_t t = clock();
Expand All @@ -196,7 +212,7 @@ void CvodeDense::SolveEq(const double tolfac, const double tmax,
/*reinitialize ode*/
flag = CVodeReInit(cvode_mem_, ode_.t_, ode_.y_);
CheckFlag(&flag, "CVodeReInit", 1);

/*copy y to yprev*/
for (int i=0; i<kDimen; i++) {
NV_Ith_S(yprev, i) = NV_Ith_S(ode_.y_, i);
Expand All @@ -205,7 +221,7 @@ void CvodeDense::SolveEq(const double tolfac, const double tmax,
Solve(ode_.t_ + tevol);
/*Calculate the relative differences of y from last timestep*/
for (int i=0; i<kDimen; i++) {
err[i] = fabs( NV_Ith_S(ode_.y_, i) / ( NV_Ith_S(yprev, i) + small__ )
err[i] = fabs( NV_Ith_S(ode_.y_, i) / ( NV_Ith_S(yprev, i) + small__ )
- 1. );
/*not consider aboundances below 10 times absolute tolerance*/
if ( NV_Ith_S(ode_.y_, i) < 10*NV_Ith_S(abstol_, i) ) {
Expand All @@ -219,27 +235,27 @@ void CvodeDense::SolveEq(const double tolfac, const double tmax,
ode_.t_ - tevol, ode_.t_, maxerr);
}
if ( maxerr < reltol_*tolfac ) {
N_VDestroy_Serial(yprev);
N_VDestroy(yprev);
t = clock() - t;
t_solve_eq_ = ((float)t)/CLOCKS_PER_SEC; /*record the run time*/
t_solve_eq_ = ((float)t)/CLOCKS_PER_SEC; /*record the run time*/
nst_eq_ = GetNstep() - nst;
return;
}
if (tevol > tevol_max) {
N_VDestroy_Serial(yprev);
N_VDestroy(yprev);
t = clock() - t;
t_solve_eq_ = ((float)t)/CLOCKS_PER_SEC; /*record the run time*/
t_solve_eq_ = ((float)t)/CLOCKS_PER_SEC; /*record the run time*/
nst_eq_ = GetNstep() - nst;
printf(
"Warning: not reaching equalibrium in tevol_max, residual=%0.4e.\n",
"Warning: not reaching equalibrium in tevol_max, residual=%0.4e.\n",
maxerr);
return;
}
tevol *= 2.;
}

}

void CvodeDense::PrintStats() const {
long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf, nge;
int flag;
Expand All @@ -257,10 +273,10 @@ void CvodeDense::PrintStats() const {
flag = CVodeGetNumNonlinSolvConvFails(cvode_mem_, &ncfn);
CheckFlag(&flag, "CVodeGetNumNonlinSolvConvFails", 1);

flag = CVDlsGetNumJacEvals(cvode_mem_, &nje);
CheckFlag(&flag, "CVDlsGetNumJacEvals", 1);
flag = CVDlsGetNumRhsEvals(cvode_mem_, &nfeLS);
CheckFlag(&flag, "CVDlsGetNumRhsEvals", 1);
flag = CVodeGetNumJacEvals(cvode_mem_, &nje);
CheckFlag(&flag, "CVodeGetNumJacEvals", 1);
flag = CVodeGetNumLinRhsEvals(cvode_mem_, &nfeLS);
CheckFlag(&flag, "CVodeGetNumLinRhsEvals", 1);

flag = CVodeGetNumGEvals(cvode_mem_, &nge);
CheckFlag(&flag, "CVodeGetNumGEvals", 1);
Expand Down
11 changes: 6 additions & 5 deletions cvodeDense.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,11 +15,11 @@
#include <time.h> /*clock_t, clock, CLOCKS_PER_SEC*/
#include <algorithm> /*std::max, std::max_element*/
#include "ode.h"
#include "sundial.h" /*Ith, IJth macro and CheckFlag function*/
#include "sundial.h" /*CheckFlag function*/

class CvodeDense {
public:
CvodeDense(Ode &ode,
CvodeDense(Ode &ode,
const double reltol, const double *abstol,
const bool userJac=false);
~CvodeDense();
Expand All @@ -39,7 +39,7 @@ class CvodeDense {
* tolfac: conservative factor of tolerance, default 10. The real torlerance
* is set as the relative torlerance times the tolfac.
* tmax: maximum time for reaching equalibrium, default 3.16e13 s or * * 1Myr.*/
void SolveEq(const double tolfac=10, const double tmax=3.16e13,
void SolveEq(const double tolfac=10, const double tmax=3.16e13,
const bool verbose=false, const double tmin=3.16e10);
void PrintStats() const;
/*write t_solve_ and t_solve_eq to file*/
Expand All @@ -62,14 +62,15 @@ class CvodeDense {
const int kDimen;
realtype reltol_;
N_Vector abstol_;
SUNContext sunctx_;
SUNMatrix dense_matrix_;
SUNLinearSolver dense_ls_;
void *cvode_mem_;
/*record the run time and number of steps for the last time calling CVode
/*record the run time and number of steps for the last time calling CVode
* solver in Solve function.*/
double t_solve_;
long int nst_last_;
/*recored the run time and number of steps for the last time calling SolveEq
/*recored the run time and number of steps for the last time calling SolveEq
* to evolve to equalibrium.*/
double t_solve_eq_;
long int nst_eq_;
Expand Down
5 changes: 0 additions & 5 deletions gow17.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -291,9 +291,6 @@ gow17::gow17()
fCplusCR_(1.),
fHeplusgr_(1.)
{
y_ = N_VNew_Serial(kDimen); /* Create serial vector y*/
CheckFlag((void *)y_, "N_VNew_Serial", 0);

/*Default values of physical parameters*/
Zg_ = 1.;
Zd_ = 1.;
Expand All @@ -315,8 +312,6 @@ gow17::gow17()
}

gow17::~gow17() {
/* destroy serial vector y*/
N_VDestroy_Serial(y_);
}

/*-------------------------RHS and Jac-----------------------------*/
Expand Down
6 changes: 3 additions & 3 deletions ode.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
#include <nvector/nvector_serial.h> /* N_Vector type*/
#include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix */
#include <stdio.h>
#include "sundial.h" /*Ith IJth macro and CheckFlag function*/
#include "sundial.h" /*CheckFlag function*/
class Ode {
friend class CvodeDense;
public:
Expand All @@ -25,7 +25,7 @@ class Ode {
*Note that user_data is always set to pointer to Ode class.*/
virtual void WriteStat (FILE *pf) const;
static int wrapJac(const realtype t,
const N_Vector y, const N_Vector fy,
const N_Vector y, const N_Vector fy,
SUNMatrix J, void *user_data,
N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
static int wrapRHS(const realtype t, const N_Vector y,
Expand All @@ -41,7 +41,7 @@ class Ode {
N_Vector ydot) = 0;
/*Jacobian*/
virtual int Jac(const realtype t,
const N_Vector y, const N_Vector fy,
const N_Vector y, const N_Vector fy,
SUNMatrix J, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) = 0;
protected:
/*status of ODE. Note the the ODE solver will change this. */
Expand Down
2 changes: 1 addition & 1 deletion out_example_simple/prof000000.dat
View file

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion out_example_simple/prof000001.dat
View file

Large diffs are not rendered by default.

Loading

0 comments on commit e8fa1da

Please sign in to comment.