diff --git a/examples/arkode/C_serial/ark_damped_harmonic_symplectic.c b/examples/arkode/C_serial/ark_damped_harmonic_symplectic.c index 940e5be81d..7acbf5543b 100644 --- a/examples/arkode/C_serial/ark_damped_harmonic_symplectic.c +++ b/examples/arkode/C_serial/ark_damped_harmonic_symplectic.c @@ -74,28 +74,34 @@ static int check_retval(void* returnvalue, const char* funcname, int opt); int main(int argc, char* argv[]) { ProgramArgs args; - SUNContext sunctx = NULL; - N_Vector y = NULL; - sunrealtype* ydata = NULL; - sunrealtype tout = NAN; - sunrealtype tret = NAN; - void* arkode_mem = NULL; - int iout = 0; - int retval = 0; + SUNContext sunctx = NULL; + N_Vector y = NULL; + sunrealtype* ydata = NULL; + sunrealtype tout = NAN; + sunrealtype tret = NAN; + void* arkode_mem = NULL; + int iout = 0; + int retval = 0; + int order = 0; + int use_compsums = 0; + int num_output_times = 0; + sunrealtype Tf = SUN_RCONST(0.0); + sunrealtype dt = SUN_RCONST(0.0); + sunrealtype dTout = SUN_RCONST(0.0); + const sunrealtype T0 = SUN_RCONST(0.0); /* Parse the command line arguments */ if (ParseArgs(argc, argv, &args)) { return 1; }; /* Default integrator options */ - int order = args.order; - int use_compsums = args.use_compsums; - const int num_output_times = args.num_output_times; + order = args.order; + use_compsums = args.use_compsums; + num_output_times = args.num_output_times; /* Default problem parameters */ - const sunrealtype T0 = SUN_RCONST(0.0); - sunrealtype Tf = args.Tf; - sunrealtype dt = args.dt; - const sunrealtype dTout = (Tf - T0) / ((sunrealtype)num_output_times); + Tf = args.Tf; + dt = args.dt; + dTout = (Tf - T0) / ((sunrealtype)num_output_times); /* Create the SUNDIALS context object for this simulation */ retval = SUNContext_Create(NULL, &sunctx); @@ -206,6 +212,8 @@ int pdot(sunrealtype t, N_Vector yvec, N_Vector ydotvec, void* user_data) int ParseArgs(int argc, char* argv[], ProgramArgs* args) { + int argi = 0; + args->order = 4; args->num_output_times = 8; args->use_compsums = 0; @@ -213,7 +221,7 @@ int ParseArgs(int argc, char* argv[], ProgramArgs* args) args->Tf = SUN_RCONST(10.0) * PI; args->dt = SUN_RCONST(1e-3); - for (int argi = 1; argi < argc; argi++) + for (argi = 1; argi < argc; argi++) { if (!strcmp(argv[argi], "--order")) { diff --git a/examples/arkode/C_serial/ark_harmonic_symplectic.c b/examples/arkode/C_serial/ark_harmonic_symplectic.c index b9d5966db8..400fb73527 100644 --- a/examples/arkode/C_serial/ark_harmonic_symplectic.c +++ b/examples/arkode/C_serial/ark_harmonic_symplectic.c @@ -82,33 +82,39 @@ int main(int argc, char* argv[]) { ProgramArgs args; UserData udata; - SUNContext sunctx = NULL; - N_Vector y = NULL; - N_Vector solution = NULL; - sunrealtype* ydata = NULL; - sunrealtype tout = NAN; - sunrealtype tret = NAN; - sunrealtype err = NAN; - void* arkode_mem = NULL; - int iout = 0; - int retval = 0; + SUNContext sunctx = NULL; + N_Vector y = NULL; + N_Vector solution = NULL; + sunrealtype* ydata = NULL; + sunrealtype tout = NAN; + sunrealtype tret = NAN; + sunrealtype err = NAN; + void* arkode_mem = NULL; + int iout = 0; + int retval = 0; + int order = 0; + int use_compsums = 0; + int num_output_times = 0; + sunrealtype Tf = SUN_RCONST(0.0); + sunrealtype dt = SUN_RCONST(0.0); + sunrealtype dTout = SUN_RCONST(0.0); + const sunrealtype T0 = SUN_RCONST(0.0); + const sunrealtype A = SUN_RCONST(10.0); + const sunrealtype phi = SUN_RCONST(0.0); + const sunrealtype omega = SUN_RCONST(1.0); /* Parse the command line arguments */ if (ParseArgs(argc, argv, &args)) { return 1; }; - /* Default integrator options */ - int order = args.order; - int use_compsums = args.use_compsums; - const int num_output_times = args.num_output_times; + /* Default integrator options and problem parameters */ + order = args.order; + use_compsums = args.use_compsums; + num_output_times = args.num_output_times; + Tf = args.Tf; + dt = args.dt; + dTout = (Tf - T0) / ((sunrealtype)num_output_times); /* Default problem parameters */ - const sunrealtype T0 = SUN_RCONST(0.0); - sunrealtype Tf = args.Tf; - sunrealtype dt = args.dt; - const sunrealtype A = SUN_RCONST(10.0); - const sunrealtype phi = SUN_RCONST(0.0); - const sunrealtype omega = SUN_RCONST(1.0); - const sunrealtype dTout = (Tf - T0) / ((sunrealtype)num_output_times); /* Create the SUNDIALS context object for this simulation */ retval = SUNContext_Create(NULL, &sunctx); @@ -249,6 +255,8 @@ int vdot(sunrealtype t, N_Vector yvec, N_Vector ydotvec, void* user_data) int ParseArgs(int argc, char* argv[], ProgramArgs* args) { + int argi = 0; + args->order = 4; args->num_output_times = 8; args->use_compsums = 0; @@ -256,7 +264,7 @@ int ParseArgs(int argc, char* argv[], ProgramArgs* args) args->dt = SUN_RCONST(1e-3); args->Tf = SUN_RCONST(2.0) * PI; - for (int argi = 1; argi < argc; argi++) + for (argi = 1; argi < argc; argi++) { if (!strcmp(argv[argi], "--order")) { diff --git a/examples/arkode/C_serial/ark_kepler.c b/examples/arkode/C_serial/ark_kepler.c index 3a1c09614a..ae4346d455 100644 --- a/examples/arkode/C_serial/ark_kepler.c +++ b/examples/arkode/C_serial/ark_kepler.c @@ -153,6 +153,7 @@ int main(int argc, char* argv[]) } else { + int i = 0; /* Compute the order of accuracy of the method by testing it with different step sizes. */ sunrealtype acc_orders[NUM_DT]; @@ -167,14 +168,7 @@ int main(int argc, char* argv[]) sunrealtype ord_max_acc = 0, ord_max_conv = 0, ord_avg = 0, ord_est = 0; sunrealtype refine = SUN_RCONST(.5); sunrealtype dt = (expected_order >= 3) ? SUN_RCONST(1e-1) : SUN_RCONST(1e-3); - sunrealtype dts[NUM_DT] = {dt, - dt * refine, - dt * refine * refine, - dt * pow(refine, 3), - dt * pow(refine, 4), - dt * pow(refine, 5), - dt * pow(refine, 6), - dt * pow(refine, 7)}; + sunrealtype dts[NUM_DT]; /* Create a reference solution using 8th order ERK with a small time step */ const int old_step_mode = args.step_mode; @@ -197,8 +191,10 @@ int main(int argc, char* argv[]) args.stepper = old_stepper; args.method_name = old_method_name; + for (i = 0; i < NUM_DT; i++) { dts[i] = dt * pow(refine, i); } + /* Compute the error with various step sizes */ - for (int i = 0; i < NUM_DT; i++) + for (i = 0; i < NUM_DT; i++) { /* Set the dt to use for this solve */ args.dt = dts[i]; @@ -636,9 +632,10 @@ int ComputeConvergence(int num_dt, sunrealtype* orders, sunrealtype* ord_max, sunrealtype* ord_est) { /* Compute/print overall estimated convergence rate */ + int i = 0; sunrealtype det = 0; *ord_avg = 0, *ord_max = 0, *ord_est = 0; - for (int i = 1; i < num_dt; i++) + for (i = 1; i < num_dt; i++) { *ord_avg += orders[i - 1]; *ord_max = SUNMAX(*ord_max, orders[i - 1]); @@ -651,6 +648,8 @@ int ComputeConvergence(int num_dt, sunrealtype* orders, int ParseArgs(int argc, char* argv[], ProgramArgs* args) { + int argi = 0; + args->step_mode = 0; args->stepper = 0; args->method_name = NULL; @@ -662,7 +661,7 @@ int ParseArgs(int argc, char* argv[], ProgramArgs* args) args->check_order = 0; args->num_output_times = 50; - for (int argi = 1; argi < argc; argi++) + for (argi = 1; argi < argc; argi++) { if (!strcmp(argv[argi], "--step-mode")) { diff --git a/src/arkode/arkode.c b/src/arkode/arkode.c index e0964ba0f0..fb5142b78c 100644 --- a/src/arkode/arkode.c +++ b/src/arkode/arkode.c @@ -27,6 +27,7 @@ #include "arkode_impl.h" #include "arkode_interp_impl.h" +#include "sundials/sundials_config.h" #include #include @@ -2333,7 +2334,7 @@ int arkYddNorm(ARKodeMem ark_mem, realtype hg, realtype *yddnrm) return(ARK_SUCCESS); } -inline static +SUNDIALS_STATIC_INLINE void compensatedSum(sunrealtype base, sunrealtype inc, sunrealtype *sum, sunrealtype *error) { sunrealtype err = *error; diff --git a/src/arkode/arkode_sprk.c b/src/arkode/arkode_sprk.c index e161699f2d..f4e54c2531 100644 --- a/src/arkode/arkode_sprk.c +++ b/src/arkode/arkode_sprk.c @@ -138,17 +138,19 @@ ARKodeSPRKStorage ARKodeSymplecticMcLachlan2() ARKodeSPRKStorage ARKodeSymplecticMcLachlan3() { + sunrealtype w = 0.0; + sunrealtype y = 0.0; + sunrealtype z = 0.0; ARKodeSPRKStorage sprk_storage = ARKodeSPRKStorage_Alloc(3); sprk_storage->q = 3; sprk_storage->stages = 3; - sunrealtype z = - -SUNRpowerR((SUN_RCONST(2.0) / SUN_RCONST(27.0)) - - SUN_RCONST(1.0) / (SUN_RCONST(9.0) * SUNRsqrt(3.0)), - SUN_RCONST(1.0) / SUN_RCONST(3.0)); - sunrealtype w = -SUN_RCONST(2.0) / SUN_RCONST(3.0) + - SUN_RCONST(1.0) / (SUN_RCONST(9.0) * z) + z; - sunrealtype y = (SUN_RCONST(1.0) + w * w) / SUN_RCONST(4.0); + z = -SUNRpowerR((SUN_RCONST(2.0) / SUN_RCONST(27.0)) - + SUN_RCONST(1.0) / (SUN_RCONST(9.0) * SUNRsqrt(3.0)), + SUN_RCONST(1.0) / SUN_RCONST(3.0)); + w = -SUN_RCONST(2.0) / SUN_RCONST(3.0) + + SUN_RCONST(1.0) / (SUN_RCONST(9.0) * z) + z; + y = (SUN_RCONST(1.0) + w * w) / SUN_RCONST(4.0); sprk_storage->a[0] = SUNRsqrt(SUN_RCONST(1.0) / (SUN_RCONST(9.0) * y) - w / SUN_RCONST(2.0) + SUNRsqrt(y)) - SUN_RCONST(1.0) / (SUN_RCONST(3.0) * SUNRsqrt(y)); @@ -211,16 +213,16 @@ ARKodeSPRKStorage ARKodeSymplecticYoshida6() ARKodeSPRKStorage sprk_storage = ARKodeSPRKStorage_Alloc(8); sprk_storage->q = 6; sprk_storage->stages = 8; - sprk_storage->a[0] = SUN_RCONST(0.7845136104775572638194976338663498757768); - sprk_storage->a[1] = SUN_RCONST(0.2355732133593581336847931829785346016865); - sprk_storage->a[2] = -SUN_RCONST(1.177679984178871006946415680964315734639); - sprk_storage->a[3] = SUN_RCONST(1.315186320683911218884249728238862514352); - sprk_storage->a[4] = sprk_storage->a[2]; - sprk_storage->a[5] = sprk_storage->a[1]; - sprk_storage->a[6] = sprk_storage->a[0]; - sprk_storage->a[7] = SUN_RCONST(0.0); - sprk_storage->ahat[0] = sprk_storage->a[0] / SUN_RCONST(2.0); - sprk_storage->ahat[1] = (sprk_storage->a[0] + sprk_storage->a[1]) / + sprk_storage->a[0] = SUN_RCONST(0.7845136104775572638194976338663498757768); + sprk_storage->a[1] = SUN_RCONST(0.2355732133593581336847931829785346016865); + sprk_storage->a[2] = -SUN_RCONST(1.177679984178871006946415680964315734639); + sprk_storage->a[3] = SUN_RCONST(1.315186320683911218884249728238862514352); + sprk_storage->a[4] = sprk_storage->a[2]; + sprk_storage->a[5] = sprk_storage->a[1]; + sprk_storage->a[6] = sprk_storage->a[0]; + sprk_storage->a[7] = SUN_RCONST(0.0); + sprk_storage->ahat[0] = sprk_storage->a[0] / SUN_RCONST(2.0); + sprk_storage->ahat[1] = (sprk_storage->a[0] + sprk_storage->a[1]) / SUN_RCONST(2.0); sprk_storage->ahat[2] = (sprk_storage->a[1] + sprk_storage->a[2]) / SUN_RCONST(2.0); @@ -238,7 +240,7 @@ ARKodeSPRKStorage ARKodeSymplecticYoshida6() (Original) Suzuki, M., & Umeno, K. (1993). Higher-order decomposition theory of exponential operators and its applications to QMC and nonlinear dynamics. - Computer simulation studies in condensed-matter physics VI, 74-86. + Computer simulation studies in condensed-matter physics VI, 74-86. https://doi.org/10.1007/978-3-642-78448-4_7 McLachlan, R.I.: On the Numerical Integration of Ordinary Differential @@ -252,24 +254,24 @@ ARKodeSPRKStorage ARKodeSymplecticSuzukiUmeno816() ARKodeSPRKStorage sprk_storage = ARKodeSPRKStorage_Alloc(16); sprk_storage->q = 8; sprk_storage->stages = 16; - sprk_storage->a[0] = SUN_RCONST(0.7416703643506129534482278017838063156035); - sprk_storage->a[1] = -SUN_RCONST(0.4091008258000315939973000958935634173099); - sprk_storage->a[2] = SUN_RCONST(0.1907547102962383799538762564503716627355); - sprk_storage->a[3] = -SUN_RCONST(0.5738624711160822666563877266355357421595); - sprk_storage->a[4] = SUN_RCONST(0.2990641813036559238444635406886029882258); - sprk_storage->a[5] = SUN_RCONST(0.3346249182452981837849579798821822886337); - sprk_storage->a[6] = SUN_RCONST(0.3152930923967665966320566638110024309941); - sprk_storage->a[7] = -SUN_RCONST(0.7968879393529163540197888401737330534463); - sprk_storage->a[8] = sprk_storage->a[6]; - sprk_storage->a[9] = sprk_storage->a[5]; - sprk_storage->a[10] = sprk_storage->a[4]; - sprk_storage->a[11] = sprk_storage->a[3]; - sprk_storage->a[12] = sprk_storage->a[2]; - sprk_storage->a[13] = sprk_storage->a[1]; - sprk_storage->a[14] = sprk_storage->a[0]; - sprk_storage->a[15] = SUN_RCONST(0.0); - sprk_storage->ahat[0] = sprk_storage->a[0] / SUN_RCONST(2.0); - sprk_storage->ahat[1] = (sprk_storage->a[0] + sprk_storage->a[1]) / + sprk_storage->a[0] = SUN_RCONST(0.7416703643506129534482278017838063156035); + sprk_storage->a[1] = -SUN_RCONST(0.4091008258000315939973000958935634173099); + sprk_storage->a[2] = SUN_RCONST(0.1907547102962383799538762564503716627355); + sprk_storage->a[3] = -SUN_RCONST(0.5738624711160822666563877266355357421595); + sprk_storage->a[4] = SUN_RCONST(0.2990641813036559238444635406886029882258); + sprk_storage->a[5] = SUN_RCONST(0.3346249182452981837849579798821822886337); + sprk_storage->a[6] = SUN_RCONST(0.3152930923967665966320566638110024309941); + sprk_storage->a[7] = -SUN_RCONST(0.7968879393529163540197888401737330534463); + sprk_storage->a[8] = sprk_storage->a[6]; + sprk_storage->a[9] = sprk_storage->a[5]; + sprk_storage->a[10] = sprk_storage->a[4]; + sprk_storage->a[11] = sprk_storage->a[3]; + sprk_storage->a[12] = sprk_storage->a[2]; + sprk_storage->a[13] = sprk_storage->a[1]; + sprk_storage->a[14] = sprk_storage->a[0]; + sprk_storage->a[15] = SUN_RCONST(0.0); + sprk_storage->ahat[0] = sprk_storage->a[0] / SUN_RCONST(2.0); + sprk_storage->ahat[1] = (sprk_storage->a[0] + sprk_storage->a[1]) / SUN_RCONST(2.0); sprk_storage->ahat[2] = (sprk_storage->a[1] + sprk_storage->a[2]) / SUN_RCONST(2.0); @@ -428,8 +430,7 @@ ARKodeSPRKStorage ARKodeSPRKStorage_Load(ARKODE_SPRKMethodID id) case ARKODE_SPRK_MCLACHLAN_2_2: return ARKodeSymplecticMcLachlan2(); case ARKODE_SPRK_MCLACHLAN_3_3: return ARKodeSymplecticMcLachlan3(); case ARKODE_SPRK_MCLACHLAN_4_4: return ARKodeSymplecticMcLachlan4(); - case ARKODE_SPRK_CANDY_ROZMUS_4_4: - return ARKodeSymplecticCandyRozmus4(); + case ARKODE_SPRK_CANDY_ROZMUS_4_4: return ARKodeSymplecticCandyRozmus4(); case ARKODE_SPRK_MCLACHLAN_5_6: return ARKodeSymplecticMcLachlan5(); case ARKODE_SPRK_YOSHIDA_6_8: return ARKodeSymplecticYoshida6(); case ARKODE_SPRK_SUZUKI_UMENO_8_16: return ARKodeSymplecticSuzukiUmeno816();