diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/IndexArrays.c | 6 | ||||
-rw-r--r-- | src/ParamCheck.c | 8 | ||||
-rw-r--r-- | src/RK45.c | 61 | ||||
-rw-r--r-- | src/SetTime.c | 44 |
4 files changed, 103 insertions, 16 deletions
diff --git a/src/IndexArrays.c b/src/IndexArrays.c index 505b1ab..5e59734 100644 --- a/src/IndexArrays.c +++ b/src/IndexArrays.c @@ -306,7 +306,11 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS) } else if (CCTK_EQUALS(ODE_Method,"RK45")) { - sprintf(infoline, "Runge-Kutta 45"); + sprintf(infoline, "Runge-Kutta 45 (Fehlberg)"); + } + else if (CCTK_EQUALS(ODE_Method,"RK45CK")) + { + sprintf(infoline, "Runge-Kutta 45 (Cash-Karp)"); } else if (CCTK_EQUALS(ODE_Method,"RK65")) { diff --git a/src/ParamCheck.c b/src/ParamCheck.c index 272117b..c2346b0 100644 --- a/src/ParamCheck.c +++ b/src/ParamCheck.c @@ -128,10 +128,10 @@ int MoL_ParamCheck(CCTK_ARGUMENTS) "number of intermediate steps must be 3"); } - if ( (CCTK_Equals(ODE_Method, "RK45")) && + if ( (CCTK_Equals(ODE_Method, "RK45") || CCTK_Equals(ODE_Method, "RK45CK")) && ( !((MoL_Intermediate_Steps == 6)||(MoL_Num_Scratch_Levels > 5)) ) ) { - CCTK_PARAMWARN("When using the RK45 evolver the " + CCTK_PARAMWARN("When using the RK45 or RK45CK evolver, the " "number of intermediate steps must be 6 " "and the number of scratch levels at least 6."); } @@ -154,13 +154,13 @@ int MoL_ParamCheck(CCTK_ARGUMENTS) if (adaptive_stepsize) { - if (CCTK_Equals(ODE_Method, "RK45")||CCTK_Equals(ODE_Method, "RK65")||CCTK_Equals(ODE_Method, "RK87")) + if (CCTK_Equals(ODE_Method, "RK45")||CCTK_Equals(ODE_Method, "RK45CK")||CCTK_Equals(ODE_Method, "RK65")||CCTK_Equals(ODE_Method, "RK87")) { /* everything is fine, do nothing */ } else { - CCTK_PARAMWARN("Adaptive time step sizes are only possible with the RK45, RK65 and RK87 solvers"); + CCTK_PARAMWARN("Adaptive time step sizes are only possible with the RK45, RK45CK, RK65, and RK87 solvers"); } } @@ -67,11 +67,9 @@ void MoL_RK45Add(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; - cGroupDynamicData arraydata; - CCTK_INT groupindex, ierr; - CCTK_INT arraytotalsize, arraydim; + CCTK_INT arraydim; - CCTK_INT index, var, scratchstep, alphaindex, scratchindex; + CCTK_INT index, var, scratchstep; CCTK_INT totalsize; CCTK_REAL * restrict UpdateVar; @@ -84,7 +82,8 @@ void MoL_RK45Add(CCTK_ARGUMENTS) CCTK_REAL beta, gamma, gamma_error; - static const CCTK_REAL beta_array[5][5] = { + /* Fehlberg coefficients */ + static const CCTK_REAL beta_array_F[5][5] = { { 0.25, 0.0, 0.0, 0.0, 0.0 }, { 0.09375, 0.28125, 0.0, 0.0, 0.0 }, { 0.879380974055530268548020027310, -3.27719617660446062812926718252, 3.32089212562585343650432407829, 0.0, 0.0 }, @@ -92,7 +91,7 @@ void MoL_RK45Add(CCTK_ARGUMENTS) { -0.296296296296296296296296296296, 2.0, -1.38167641325536062378167641326, 0.452972709551656920077972709552, -0.275 } }; - static const CCTK_REAL gamma_array[6] = + static const CCTK_REAL gamma_array_F[6] = { 0.118518518518518518518518518519, 0.0, 0.518986354775828460038986354776, @@ -101,7 +100,7 @@ void MoL_RK45Add(CCTK_ARGUMENTS) 0.0363636363636363636363636363636 }; - static const CCTK_REAL gammastar_array[6] = + static const CCTK_REAL gammastar_array_F[6] = { 0.115740740740740740740740740741, 0.0, 0.548927875243664717348927875244, @@ -110,6 +109,54 @@ void MoL_RK45Add(CCTK_ARGUMENTS) 0.0 }; + /* Cash-Karp coefficients */ + static const CCTK_REAL beta_array_CK[5][5] = { + { 1.0/5.0, 0.0, 0.0, 0.0, 0.0, }, + { 3.0/40.0, 9.0/40.0, 0.0, 0.0, 0.0, }, + { 3.0/10.0, -9.0/10.0, 6.0/5.0, 0.0, 0.0, }, + { -11.0/54.0, 5.0/2.0, -70.0/27.0, 35.0/27.0, 0.0, }, + { 1631.0/55296.0, 175.0/512.0, 575.0/13824.0, 44275.0/110592.0, 253.0/4096.0, }, + }; + + static const CCTK_REAL gamma_array_CK[6] = { + 37.0/378.0, + 0.0, + 250.0/621.0, + 125.0/594.0, + 0.0, + 512.0/1771.0, + }; + + static const CCTK_REAL gammastar_array_CK[6] = { + 2825.0/27648.0, + 0.0, + 18575.0/48384.0, + 13525.0/55296.0, + 277.0/14336.0, + 1.0/4.0, + }; + + const CCTK_REAL (* restrict beta_array)[5]; + const CCTK_REAL * restrict gamma_array; + const CCTK_REAL * restrict gammastar_array; + + if (CCTK_EQUALS(ODE_Method, "RK45")) + { + beta_array = beta_array_F; + gamma_array = gamma_array_F; + gammastar_array = gammastar_array_F; + } + else if (CCTK_EQUALS(ODE_Method, "RK45CK")) + { + beta_array = beta_array_CK; + gamma_array = gamma_array_CK; + gammastar_array = gammastar_array_CK; + } + else + { + CCTK_WARN (0, "internal error"); + } + totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { diff --git a/src/SetTime.c b/src/SetTime.c index 1d1596d..3e9ae3f 100644 --- a/src/SetTime.c +++ b/src/SetTime.c @@ -47,8 +47,8 @@ int MoL_ResetDeltaTime(CCTK_ARGUMENTS); ********************* Local Data ***************************** ********************************************************************/ -/* RK45 coefficients */ -static const CCTK_REAL alpha_array[6] = { +/* RK45 Fehlberg coefficients */ +static const CCTK_REAL alpha_array_F[6] = { 0.0, 1.0/4.0, 3.0/8.0, @@ -57,6 +57,16 @@ static const CCTK_REAL alpha_array[6] = { 1.0/2.0, }; +/* RK45 Cash-Karp coefficients */ +static const CCTK_REAL alpha_array_CK[6] = { + 0.0, + 1.0/5.0, + 3.0/10.0, + 3.0/5.0, + 1.0, + 7.0/8.0, +}; + /* RK65 coefficients */ static const CCTK_REAL alpha_array65[8] = { 0.0, @@ -237,9 +247,22 @@ int MoL_ResetTime(CCTK_ARGUMENTS) 0.5*(*Original_Delta_Time)/cctkGH->cctk_timefac; } } - else if (CCTK_EQUALS(ODE_Method,"RK45")) + else if (CCTK_EQUALS(ODE_Method,"RK45") || CCTK_EQUALS(ODE_Method,"RK45CK")) { const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step); + const CCTK_REAL * alpha_array; + if (CCTK_EQUALS(ODE_Method, "RK45")) + { + alpha_array = alpha_array_F; + } + else if (CCTK_EQUALS(ODE_Method, "RK45CK")) + { + alpha_array = alpha_array_CK; + } + else + { + CCTK_WARN (0, "internal error"); + } cctkGH->cctk_time = ((* Original_Time) + ((alpha_array[substep] - 1) @@ -339,9 +362,22 @@ int MoL_ResetDeltaTime(CCTK_ARGUMENTS) cctkGH->cctk_delta_time = 2.0/3.0*(*Original_Delta_Time); } } - else if (CCTK_EQUALS(ODE_Method,"RK45")) + else if (CCTK_EQUALS(ODE_Method,"RK45") || CCTK_EQUALS(ODE_Method,"RK45CK")) { const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step); + const CCTK_REAL * alpha_array; + if (CCTK_EQUALS(ODE_Method, "RK45")) + { + alpha_array = alpha_array_F; + } + else if (CCTK_EQUALS(ODE_Method, "RK45CK")) + { + alpha_array = alpha_array_CK; + } + else + { + CCTK_WARN (0, "internal error"); + } cctkGH->cctk_delta_time = ((alpha_array[substep + 1] - alpha_array[substep]) * (* Original_Delta_Time)); |