Close
#include <stdio.h>
#include <cvodes/cvodes.h>
#include <cvodes/cvodes_dense.h>
#include <nvector/nvector_serial.h>
#define Ith(v,i) NV_Ith_S(v,i-1)
#define NEQ 3
#define Y1 RCONST(1.0)
#define Y2 RCONST(0.0)
#define Y3 RCONST(0.0)
#define RTOL RCONST(1.0e-4)
#define ATOL1 RCONST(7.0e-7)
#define ATOL2 RCONST(1.0e-12)
#define ATOL3 RCONST(1.0e-5)
#define T0 RCONST(0.0)
#define T1 RCONST(0.4)
#define TMULT RCONST(10.0)
#define NOUT 12
static int f(realtype t, N_Vector y, N_Vector ydot, void *f_data);
static void PrintOutput(realtype t, realtype y1, realtype y2, realtype y3);
static void PrintFinalStats(void *cvode_mem);
int main()
{
realtype reltol, t, tout;
N_Vector y, abstol;
void *cvode_mem;
int iout;
booleantype check_negative;
y = N_VNew_Serial(NEQ);
abstol = N_VNew_Serial(NEQ);
Ith(y,1) = Y1;
Ith(y,2) = Y2;
Ith(y,3) = Y3;
reltol = RTOL;
Ith(abstol,1) = ATOL1;
Ith(abstol,2) = ATOL2;
Ith(abstol,3) = ATOL3;
cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
CVodeMalloc(cvode_mem, f, T0, y, CV_SV, reltol, abstol);
CVodeSetFdata(cvode_mem, &check_negative);
CVDense(cvode_mem, NEQ);
printf("Ignore negative solution components\n\n");
check_negative = FALSE;
iout = 0; tout = T1;
while(1) {
CVode(cvode_mem, tout, y, &t, CV_NORMAL);
PrintOutput(t, Ith(y,1), Ith(y,2), Ith(y,3));
iout++;
tout *= TMULT;
if (iout == NOUT) break;
}
PrintFinalStats(cvode_mem);
printf("Intercept negative solution components\n\n");
check_negative = TRUE;
Ith(y,1) = Y1;
Ith(y,2) = Y2;
Ith(y,3) = Y3;
CVodeReInit(cvode_mem, f, T0, y, CV_SV, reltol, abstol);
iout = 0; tout = T1;
while(1) {
CVode(cvode_mem, tout, y, &t, CV_NORMAL);
PrintOutput(t, Ith(y,1), Ith(y,2), Ith(y,3));
iout++;
tout *= TMULT;
if (iout == NOUT) break;
}
PrintFinalStats(cvode_mem);
N_VDestroy_Serial(y);
CVodeFree(&cvode_mem);
return(0);
}
static int f(realtype t, N_Vector y, N_Vector ydot, void *f_data)
{
realtype y1, y2, y3, yd1, yd3;
booleantype *check_negative;
check_negative = (booleantype *)f_data;
y1 = Ith(y,1); y2 = Ith(y,2); y3 = Ith(y,3);
if ( *check_negative && (y1<0 || y2<0 || y3<0) )
return(1);
Ith(ydot,1) = yd1 = RCONST(-0.04)*y1 + RCONST(1.0e4)*y2*y3;
Ith(ydot,3) = yd3 = RCONST(3.0e7)*y2*y2;
Ith(ydot,2) = -yd1 - yd3;
return(0);
}
static void PrintOutput(realtype t, realtype y1, realtype y2, realtype y3)
{
printf("At t = %0.4le y =%14.6le %14.6le %14.6le\n", t, y1, y2, y3);
return;
}
static void PrintFinalStats(void *cvode_mem)
{
long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf;
CVodeGetNumSteps(cvode_mem, &nst);
CVodeGetNumRhsEvals(cvode_mem, &nfe);
CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
CVodeGetNumErrTestFails(cvode_mem, &netf);
CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
CVDenseGetNumJacEvals(cvode_mem, &nje);
CVDenseGetNumRhsEvals(cvode_mem, &nfeLS);
printf("\nFinal Statistics:\n");
printf("nst = %-6ld nfe = %-6ld nsetups = %-6ld nfeLS = %-6ld nje = %ld\n",
nst, nfe, nsetups, nfeLS, nje);
printf("nni = %-6ld ncfn = %-6ld netf = %-6ld\n \n",
nni, ncfn, netf);
}
Close