Skip to content

Commit

Permalink
add plotting files
Browse files Browse the repository at this point in the history
  • Loading branch information
@balos1
balos1 committed Apr 18, 2024
1 parent 6f91a8b commit 9204c82
Show file tree
Hide file tree
Showing 10 changed files with 403 additions and 125 deletions.
4 changes: 2 additions & 2 deletions examples/arkode/C_serial/ark_brusselator.c
View file
Expand Up @@ -112,7 +112,7 @@ int main(void)
if (check_flag(&flag, "SUNContext_Create", 1)) { return 1; }

/* set up the test problem according to the desired test */
if (test == 1)
if (test == 0)
{
u0 = SUN_RCONST(3.9);
v0 = SUN_RCONST(1.1);
Expand All @@ -121,7 +121,7 @@ int main(void)
b = SUN_RCONST(2.5);
ep = SUN_RCONST(1.0e-5);
}
else if (test == 3)
else if (test == 1)
{
u0 = SUN_RCONST(3.0);
v0 = SUN_RCONST(3.0);
Expand Down
179 changes: 109 additions & 70 deletions examples/cvode/serial/cvBrusselator.c
View file
Expand Up @@ -50,26 +50,41 @@
* (dense direct, GMRES with the Jacobian computed by difference quotients,
* or GMRES with analytical Jacobian), and the reactor type.
*
* ./cv_brusselator [solver_type] [reactor_type]
* ./cv_brusselator [reactor_type] [nonlinear_solver_type] [linear_solver_type]
*
* Options:
* num_batches <int>
* solver_type:
* 0 - SUNDIALS GMRES with difference quotients Jacobian
* 1 - SUNDIALS GMRES with analytical Jacobian
* 2 - SUNDIALS dense direct linear solver with analytical Jacobian
* reactor_type:
* 0 - Reactor 0
* 1 - Reactor 1
* 2 - Reactor 2
* nonlinear_linear_solver_type:
* 0 - Newton
* 1 - Fixedpoint
* 2 - Automatic switching (via GS algorthim) between Newton and Fixedpoint, start with Newton
* 3 - Automatic switching (via GS algorthim) between Newton and Fixedpoint, start with Fixedpoint
* 12 - Automatic switching (via LSODKR algorthim) between Newton and Fixedpoint, start with Newton
* 13 - Automatic switching (via LSODKR algorthim) between Newton and Fixedpoint, start with Fixedpoint
* linear_linear_solver_type:
* 0 - SUNDIALS GMRES with difference quotients Jtv
* 1 - SUNDIALS dense linear solver with analytic Jacobian
* --------------------------------------------------------------------------*/

#include <stdio.h>
#include <cvode/cvode.h>
#include <nvector/nvector_serial.h>
#include <sunlinsol/sunlinsol_spgmr.h>
#include <sunlinsol/sunlinsol_dense.h>
#include <sundials/sundials_core.h>

#if defined(SUNDIALS_EXTENDED_PRECISION)
#define GSYM "Lg"
#define ESYM "Le"
#define FSYM "Lf"
#else
#define GSYM "g"
#define ESYM "e"
#define FSYM "f"
#endif

/* Functions Called by the Solver */
static int f(sunrealtype t, N_Vector y, N_Vector ydot, void *user_data);
Expand All @@ -78,10 +93,7 @@ static int Jac(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J,
void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);

/* Private function to output results */
static void PrintOutput(sunrealtype t, sunrealtype y1, sunrealtype y2, sunrealtype y3);

/* Private function to print final statistics */
static void PrintFinalStats(void *cvode_mem, SUNLinearSolver LS);
static void PrintOutput(FILE* fp, sunrealtype t, sunrealtype y1, sunrealtype y2, sunrealtype y3);

/* Private function to check function return values */
static int check_retval(void *returnvalue, const char *funcname, int opt);
Expand All @@ -100,38 +112,47 @@ typedef struct {

int main(int argc, char *argv[])
{
const sunrealtype T0 = SUN_RCONST(0.0); /* initial time */
const sunrealtype Tf = SUN_RCONST(10.0); /* final time */
const sunrealtype dTout = SUN_RCONST(1.0); /* time between outputs */
const int Nt = (int) ceil(Tf/dTout); /* number of output times */
const sunrealtype reltol = 1.0e-4; /* relative integrator tolerance */
const sunrealtype T0 = SUN_RCONST(0.0); /* initial time */
const sunrealtype Tf = SUN_RCONST(10.0); /* final time */
// const sunrealtype Tf = SUN_RCONST(0.1); /* final time */
// const sunrealtype dTout = SUN_RCONST(1.0); /* time between outputs */
const sunrealtype dTout = SUN_RCONST(0.01); /* time between outputs */
const int Nt = (int) ceil(Tf/dTout); /* number of output times */
const sunrealtype reltol = 1.0e-6; /* relative integrator tolerance */
const sunrealtype abstol = 1.0e-10; /* absolute integrator tolerance */

int retval;
N_Vector y, abstol;
SUNContext sunctx;
N_Vector y;
SUNMatrix A;
SUNLinearSolver LS;
void *cvode_mem;
FILE* outputfp;

y = abstol = NULL;
y = NULL;
A = NULL;
LS = NULL;
cvode_mem = NULL;
outputfp = fopen("cvBrusselator_solution.txt", "w+");

/* Create the SUNDIALS context */
SUNContext sunctx;
SUNContext_Create(SUN_COMM_NULL, &sunctx);

/* Set defaults */
int solver_type = 0;
int nonlinear_solver_type = 0;
int linear_solver_type = 0;
int reactor_type = 2;

/* Parse command line arguments and setup UserData */
int argi = 0;
// if (argc > 1) {
// solver_type = atoi(argv[++argi]);
// }
if (argc > 2) {
if (argc > 1) {
reactor_type = atoi(argv[++argi]);
}
if (argc > 2) {
nonlinear_solver_type = atoi(argv[++argi]);
}
if (argc > 3) {
linear_solver_type = atoi(argv[++argi]);
}

/* Set the Reaction parameters according to reactor_type */
UserData udata;
Expand All @@ -149,7 +170,7 @@ int main(int argc, char *argv[])
udata.a = SUN_RCONST(0.5);
udata.b = SUN_RCONST(3.0);
udata.ep = SUN_RCONST(5.0e-4);
} else if (reactor_type == 2) {
} else {
udata.u0 = SUN_RCONST(1.2);
udata.v0 = SUN_RCONST(3.1);
udata.w0 = SUN_RCONST(3.0);
Expand All @@ -160,25 +181,22 @@ int main(int argc, char *argv[])

/* Create for initial conditions the and absolute tolerance vector */
y = N_VNew_Serial(3, sunctx);
abstol = N_VClone(y);

/* Initialize y */
sunrealtype* ydata = N_VGetArrayPointer(y);
ydata[0] = udata.u0;
ydata[1] = udata.v0;
ydata[2] = udata.w0;

/* Set the vector absolute tolerance */
sunrealtype* abstol_data = N_VGetArrayPointer(abstol);
abstol_data[0] = 1.0e-10;
abstol_data[1] = 1.0e-10;
abstol_data[2] = 1.0e-10;

/* Call CVodeCreate to create the solver memory and specify the
* Backward Differentiation Formula */
cvode_mem = CVodeCreate(CV_BDF, sunctx);
if (check_retval((void *)cvode_mem, "CVodeCreate", 0)) { return(1); }

/* Must be called before CVodeInit */
retval = CVodeSetNonlinearSolverAlgorithm(cvode_mem, nonlinear_solver_type);
if (check_retval(&retval, "CVodeSetNonlinearSolverAlgorithm", 1)) { return(1); }

/* 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. */
Expand All @@ -191,28 +209,51 @@ int main(int argc, char *argv[])

/* Call CVodeSVtolerances to specify the scalar relative tolerance
* and vector absolute tolerances */
retval = CVodeSVtolerances(cvode_mem, reltol, abstol);
retval = CVodeSStolerances(cvode_mem, reltol, abstol);
if (check_retval(&retval, "CVodeSVtolerances", 1)) { return(1); }

/* Create the SUNLinearSolver object for use by CVode */
if (solver_type == 0) {
if (linear_solver_type == 0) {
LS = SUNLinSol_SPGMR(y, SUN_PREC_NONE, 0, sunctx);
if(check_retval((void *)LS, "SUNLinSol_SPGMR", 0)) return(1);

/* Call CVodeSetLinearSolver to attach the matrix and linear solver to CVode */
retval = CVodeSetLinearSolver(cvode_mem, LS, NULL);
if(check_retval(&retval, "CVodeSetLinearSolver", 1)) return(1);
} else if (linear_solver_type == 1) {
A = SUNDenseMatrix(3, 3, sunctx);
if (check_retval((void*)A, "SUNDenseMatrix", 0)) { return 1; }

LS = SUNLinSol_Dense(y, A, sunctx);
if (check_retval((void*)LS, "SUNLinSol_Dense", 0)) { return 1; }

/* Call CVodeSetLinearSolver to attach the matrix and linear solver to CVode */
retval = CVodeSetLinearSolver(cvode_mem, LS, A);
if(check_retval(&retval, "CVodeSetLinearSolver", 1)) return(1);

retval = CVodeSetJacFn(cvode_mem, Jac);
if (check_retval(&retval, "CVodeSetJacFn", 1)) { return 1; }
}

/* In loop, call CVode, print results, and reactor_type for error.
Break out of loop when preset output times have been reached. */
printf(" \nBrusselator problem\n\n");
/* Initial problem output */
printf("\nBrusselator ODE test problem:\n");
printf(" initial conditions: u0 = %" GSYM ", v0 = %" GSYM
", w0 = %" GSYM "\n",
udata.u0, udata.v0, udata.w0);
printf(" problem parameters: a = %" GSYM ", b = %" GSYM ", ep = %" GSYM
"\n",
udata.a, udata.b, udata.ep);
printf(" reltol = %.1" ESYM ", abstol = %.1" ESYM "\n", reltol, abstol);
printf(" reactor_type = %d, nonlinear_solver_type = %d, linear_solver_type = %d\n\n", reactor_type, nonlinear_solver_type, linear_solver_type);

sunrealtype t = T0;
sunrealtype tout = T0+dTout;
PrintOutput(outputfp, t, ydata[0], ydata[1], ydata[2]);
for (int iout=0; iout<Nt; iout++) {
retval = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
PrintOutput(t, ydata[0], ydata[1], ydata[2]);
PrintOutput(outputfp, t, ydata[0], ydata[1], ydata[2]);
if (check_retval(&retval, "CVode", 1)) break;
if (retval == CV_SUCCESS) {
tout += dTout;
Expand All @@ -224,9 +265,8 @@ int main(int argc, char *argv[])
printf("\nFinal Statistics:\n");
CVodePrintAllStats(cvode_mem, stdout, SUN_OUTPUTFORMAT_TABLE);

/* Free y and abstol vectors */
/* Free state vector */
N_VDestroy(y);
N_VDestroy(abstol);

/* Free integrator memory */
CVodeFree(&cvode_mem);
Expand All @@ -237,6 +277,8 @@ int main(int argc, char *argv[])
/* Free the matrix memory */
if (A) { SUNMatDestroy(A); }

fclose(outputfp);

SUNContext_Free(&sunctx);

return(0);
Expand All @@ -263,7 +305,7 @@ int f(sunrealtype t, N_Vector y, N_Vector ydot, void *user_data)
ydata = N_VGetArrayPointer(y);
ydotdata = N_VGetArrayPointer(ydot);

a = udata->a; b = udata->b, ep = udata->ep;
a = udata->a; b = udata->b; ep = udata->ep;

u = ydata[0];
v = ydata[1];
Expand All @@ -278,7 +320,6 @@ int f(sunrealtype t, N_Vector y, N_Vector ydot, void *user_data)

/*
* Jacobian routine. Compute J(t,y) = df/dy.
* This is done on the GPU.
*/

int Jac(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J,
Expand All @@ -289,7 +330,7 @@ int Jac(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J,
sunrealtype u, v, w, a, b, ep;

ydata = N_VGetArrayPointer(y);
// Jdata = ?
Jdata = SUNDenseMatrix_Data(J);

a = udata->a; b = udata->b, ep = udata->ep;

Expand All @@ -299,19 +340,19 @@ int Jac(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J,
w = ydata[2];

/* first col of block */
Jdata[0] = -(w+1.0) + 2.0*u*v;
Jdata[1] = u*u;
Jdata[2] = -u;
Jdata[0] = -(w+1.0) + 2.0*u*v;
Jdata[1] = u*u;
Jdata[2] = -u;

/* second col of block */
Jdata[3] = u*u;
Jdata[4] = -u*u;
Jdata[5] = u;
Jdata[3] = w - 2.0*u*v;
Jdata[4] = -u*u;
Jdata[5] = u;

/* third col of block */
Jdata[6] = -w;
Jdata[7] = 0.0;
Jdata[8] = -1.0/ep - u;
Jdata[6] = -w;
Jdata[7] = 0.0;
Jdata[8] = -1.0/ep - u;

return(0);
}
Expand All @@ -322,14 +363,15 @@ int Jac(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix J,
*-------------------------------
*/

void PrintOutput(sunrealtype t, sunrealtype y1, sunrealtype y2, sunrealtype y3)
void PrintOutput(FILE* fp, sunrealtype t, sunrealtype y1, sunrealtype y2, sunrealtype y3)
{
#if defined(SUNDIALS_DOUBLE_PRECISION)
printf("At t = %0.4e y =%14.6e %14.6e %14.6e\n", t, y1, y2, y3);
fprintf(stdout, "At t = %0.4e y =%14.6e %14.6e %14.6e\n", t, y1, y2, y3);
fprintf(fp, "%0.4e %14.6e %14.6e %14.6e\n", t, y1, y2, y3);
#else
printf("At t = %0.4e y =%14.6e %14.6e %14.6e\n", t, y1, y2, y3);
fprintf(fp, "%0.4e %14.6e %14.6e %14.6e\n", t, y1, y2, y3);
fprintf(stdout, "At t = %0.4e y =%14.6e %14.6e %14.6e\n", t, y1, y2, y3);
#endif

return;
}

Expand All @@ -339,33 +381,30 @@ void PrintOutput(sunrealtype t, sunrealtype y1, sunrealtype y2, sunrealtype y3)
* returned NULL pointer
* opt == 1 means SUNDIALS function returns an integer value so check if
* retval < 0
* opt == 2 means function allocates memory so check if returned
* NULL pointer
*/

int check_retval(void *returnvalue, const char *funcname, int opt)
{
int *retval;

/* Check if SUNDIALS function returned NULL pointer - no memory allocated */
if (opt == 0 && returnvalue == NULL) {
fprintf(stderr, "\nSUNDIALS_ERROR: %s(failed - returned NULL pointer\n\n",
if (opt == 0 && returnvalue == NULL)
{
fprintf(stderr, "\nSUNDIALS ERROR: %s(failed - returned NULL pointer\n\n",
funcname);
return(1); }

return(1);
}
/* Check if retval < 0 */
else if (opt == 1) {
else if (opt == 1)
{
retval = (int *) returnvalue;
if (*retval < 0) {
fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with retval = %d\n\n",
if (*retval < 0)
{
fprintf(stderr, "\nSUNDIALS ERROR: %s() failed with retval = %d\n\n",
funcname, *retval);
return(1); }}

/* Check if function returned NULL pointer - no memory allocated */
else if (opt == 2 && returnvalue == NULL) {
fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
funcname);
return(1); }
return(1);
}
}

return(0);
}

0 comments on commit 9204c82

Please sign in to comment.