diff options
| -rw-r--r-- | param.ccl | 2 | ||||
| -rw-r--r-- | schedule.ccl | 18 | ||||
| -rw-r--r-- | src/IndexArrays.c | 8 | ||||
| -rw-r--r-- | src/ParamCheck.c | 20 | ||||
| -rw-r--r-- | src/RK65.c | 249 | ||||
| -rw-r--r-- | src/RK87.c | 272 | ||||
| -rw-r--r-- | src/SetTime.c | 59 | ||||
| -rw-r--r-- | src/StepSize.c | 24 | ||||
| -rw-r--r-- | src/make.code.defn | 2 |
9 files changed, 649 insertions, 5 deletions
@@ -88,6 +88,8 @@ KEYWORD ODE_Method "The ODE method use by MoL to do time integration" "RK2" :: "Efficient RK2" "RK3" :: "Efficient RK3" "RK45" :: "RK45 with error estimation" + "RK65" :: "RK65 with error estimation" + "RK87" :: "RK87 with error estimation" } "ICN" KEYWORD Generic_Type "If using the generic method, which sort" diff --git a/schedule.ccl b/schedule.ccl index add994d..3c5b82a 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -567,6 +567,24 @@ else if (CCTK_Equals(ODE_Method,"RK45")) LANG: C } "Updates calculated with the Runge-Kutta 45 method" } +else if (CCTK_Equals(ODE_Method,"RK65")) +{ + STORAGE: ErrorEstimate ErrorScalars + + schedule MoL_RK65Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep + { + LANG: C + } "Updates calculated with the Runge-Kutta 65 method" +} +else if (CCTK_Equals(ODE_Method,"RK87")) +{ + STORAGE: ErrorEstimate ErrorScalars + + schedule MoL_RK87Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep + { + LANG: C + } "Updates calculated with the Runge-Kutta 87 method" +} else if (CCTK_Equals(ODE_Method,"ICN")) { schedule MoL_ICNAdd AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep diff --git a/src/IndexArrays.c b/src/IndexArrays.c index 59f3950..505b1ab 100644 --- a/src/IndexArrays.c +++ b/src/IndexArrays.c @@ -308,6 +308,14 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS) { sprintf(infoline, "Runge-Kutta 45"); } + else if (CCTK_EQUALS(ODE_Method,"RK65")) + { + sprintf(infoline, "Runge-Kutta 65"); + } + else if (CCTK_EQUALS(ODE_Method,"RK87")) + { + sprintf(infoline, "Runge-Kutta 87"); + } else if (CCTK_EQUALS(ODE_Method,"ICN")) { sprintf(infoline, "Iterative Crank Nicholson with %i iterations", diff --git a/src/ParamCheck.c b/src/ParamCheck.c index 3bd1c21..272117b 100644 --- a/src/ParamCheck.c +++ b/src/ParamCheck.c @@ -136,15 +136,31 @@ int MoL_ParamCheck(CCTK_ARGUMENTS) "and the number of scratch levels at least 6."); } + if ( (CCTK_Equals(ODE_Method, "RK65")) && + ( !((MoL_Intermediate_Steps == 8)||(MoL_Num_Scratch_Levels > 7)) ) ) + { + CCTK_PARAMWARN("When using the RK65 evolver the " + "number of intermediate steps must be 8 " + "and the number of scratch levels at least 8."); + } + + if ( (CCTK_Equals(ODE_Method, "RK87")) && + ( !((MoL_Intermediate_Steps == 13)||(MoL_Num_Scratch_Levels > 12)) ) ) + { + CCTK_PARAMWARN("When using the RK87 evolver the " + "number of intermediate steps must be 13 " + "and the number of scratch levels at least 13."); + } + if (adaptive_stepsize) { - if (CCTK_Equals(ODE_Method, "RK45")) + if (CCTK_Equals(ODE_Method, "RK45")||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 solver"); + CCTK_PARAMWARN("Adaptive time step sizes are only possible with the RK45, RK65 and RK87 solvers"); } } diff --git a/src/RK65.c b/src/RK65.c new file mode 100644 index 0000000..22bd26b --- /dev/null +++ b/src/RK65.c @@ -0,0 +1,249 @@ + /*@@ + @file RK65.c + @date Sun May 26 03:47:15 2002 + @author Peter Diener (based on RK45.c by Ian Hawke) + @desc + RK65 following P. J. Prince and J. R. Dormand + Journal of Computational and Applied Mathematics, volume 7, no 1, 1981 + @enddesc + @version $Header$ + @@*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include "ExternalVariables.h" + +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusBase_MoL_RK65_c); + +/******************************************************************** + ********************* Local Data Types *********************** + ********************************************************************/ + +/******************************************************************** + ********************* Local Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ***************** Scheduled Routine Prototypes ********************* + ********************************************************************/ + +void MoL_RK65Add(CCTK_ARGUMENTS); + +/******************************************************************** + ********************* Other Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ********************* Local Data ***************************** + ********************************************************************/ + +/******************************************************************** + ********************* External Routines ********************** + ********************************************************************/ + + /*@@ + @routine MoL_RK65Add + @date Sun May 26 03:50:44 2002 + @author Peter Diener (based on MoL_RK45Add by Ian Hawke) + @desc + Performs a single step of a Runge-Kutta 65 type time + integration, storing the error estimate. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +void MoL_RK65Add(CCTK_ARGUMENTS) +{ + + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + cGroupDynamicData arraydata; + CCTK_INT groupindex, ierr; + CCTK_INT arraytotalsize, arraydim; + + CCTK_INT index, var, scratchstep, alphaindex, scratchindex; + CCTK_INT totalsize; + + CCTK_REAL * restrict UpdateVar; + CCTK_REAL const * restrict RHSVar; + CCTK_REAL * restrict ScratchVar; + CCTK_REAL * restrict ErrorVar; + CCTK_REAL const * restrict OldVar; + + CCTK_INT arrayscratchlocation; + + CCTK_REAL beta, gamma, gamma_error; + + static const CCTK_REAL beta_array[7][7] = { + { 1.0/10.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, + { -2.0/81.0, 20.0/81.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, + { 615.0/1372.0, -270.0/343.0, 1053.0/1372.0, 0.0, 0.0, 0.0, 0.0 }, + { 3243.0/5500.0, -54.0/55.0, 50949.0/71500.0, 4998.0/17875.0, 0.0, 0.0, 0.0 }, + { -26492.0/37125.0, 72.0/55.0, 2808.0/23375.0, -24206.0/37125.0, 338.0/459.0, 0.0, 0.0 }, + { 5561.0/2376.0, -35.0/11.0, -24117.0/31603.0, 899983.0/200772.0, -5225.0/1836.0, 3925.0/4056.0, 0.0 }, + { 465467.0/266112.0, -2945.0/1232.0, -5610201.0/14158144.0, 10513573.0/3212352.0, -424325.0/205632.0, 376225.0/454272.0, 0.0 } + }; + + static const CCTK_REAL gamma_array[8] = + { 61.0/864.0, + 0.0, + 98415.0/321776.0, + 16807.0/146016.0, + 1375.0/7344.0, + 1375.0/5408.0, + -37.0/1120.0, + 1.0/10.0 + }; + + static const CCTK_REAL gammastar_array[8] = + { 821.0/10800.0, + 0.0, + 19683.0/71825.0, + 175273.0/912600.0, + 395.0/3672.0, + 785.0/2704.0, + 3.0/50.0, + 0.0 + }; + + totalsize = 1; + for (arraydim = 0; arraydim < cctk_dim; arraydim++) + { + totalsize *= cctk_lsh[arraydim]; + } + + /* Real GFs */ + + /* First store (dt times) the rhs in the scratch array. */ + + for (var = 0; var < MoLNumEvolvedVariables; var++) + { + const CCTK_REAL tmp = (*Original_Delta_Time) / cctkGH->cctk_timefac; + + UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, + EvolvedVariableIndex[var]); + RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, + RHSVariableIndex[var]); + ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, + CCTK_FirstVarIndex("MOL::SCRATCHSPACE") + + var + + MoL_Num_Evolved_Vars * + (MoL_Intermediate_Steps - + (*MoL_Intermediate_Step))); + for (index = 0; index < totalsize; index++) + { + ScratchVar[index] = tmp * RHSVar[index]; + } + } + + + for (var = 0; var < MoLNumEvolvedVariables; var++) + { + OldVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 1, + EvolvedVariableIndex[var]); + UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, + EvolvedVariableIndex[var]); + RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, + CCTK_FirstVarIndex("MOL::SCRATCHSPACE") + + var + + MoL_Num_Evolved_Vars * + (MoL_Intermediate_Steps - + (*MoL_Intermediate_Step))); + ErrorVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, + CCTK_FirstVarIndex("MOL::ERRORESTIMATE") + + var); + + if (*MoL_Intermediate_Step - 1) + { + + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] = OldVar[index]; + } + + for (scratchstep = 0; + scratchstep < MoL_Intermediate_Steps - (*MoL_Intermediate_Step) + 1; + scratchstep++) + { + + ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, + CCTK_FirstVarIndex("MOL::SCRATCHSPACE") + + var + + MoL_Num_Evolved_Vars * scratchstep); + + beta = beta_array[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)][scratchstep]; + + if ( (beta > MoL_Tiny)||(beta < -MoL_Tiny) ) + { + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += beta * ScratchVar[index]; + } + } + + } + + } + else + { + + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] = OldVar[index]; + ErrorVar[index] = 0; + } + + for (scratchstep = 0; scratchstep < 8; scratchstep++) + { + + ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, + CCTK_FirstVarIndex("MOL::SCRATCHSPACE") + + var + + MoL_Num_Evolved_Vars * scratchstep); + + gamma = gamma_array[scratchstep]; + gamma_error = gamma - gammastar_array[scratchstep]; + + if ( (gamma > MoL_Tiny)||(gamma < -MoL_Tiny) ) + { + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += gamma * ScratchVar[index]; + ErrorVar[index] += gamma_error * ScratchVar[index]; + } + } + } + + } + + } + + /* Real arrays */ + + arrayscratchlocation = 0; + + for (var = 0; var < MoLNumEvolvedArrayVariables; var++) + { + + CCTK_WARN(0, "Peter has been too lazy to write the RK65 routine " + "out for array variables. Better send him an email..."); + + } + + return; +} + +/******************************************************************** + ********************* Local Routines ************************* + ********************************************************************/ + diff --git a/src/RK87.c b/src/RK87.c new file mode 100644 index 0000000..3ce19b6 --- /dev/null +++ b/src/RK87.c @@ -0,0 +1,272 @@ + /*@@ + @file RK87.c + @date Sun May 26 03:47:15 2002 + @author Peter Diener (based on RK45.c by Ian Hawke) + @desc + RK87 following P. J. Prince and J. R. Dormand + Journal of Computational and Applied Mathematics, volume 7, no 1, 1981 + @enddesc + @version $Header$ + @@*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include "ExternalVariables.h" + +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusBase_MoL_RK87_c); + +/******************************************************************** + ********************* Local Data Types *********************** + ********************************************************************/ + +/******************************************************************** + ********************* Local Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ***************** Scheduled Routine Prototypes ********************* + ********************************************************************/ + +void MoL_RK87Add(CCTK_ARGUMENTS); + +/******************************************************************** + ********************* Other Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ********************* Local Data ***************************** + ********************************************************************/ + +/******************************************************************** + ********************* External Routines ********************** + ********************************************************************/ + + /*@@ + @routine MoL_RK87Add + @date Sun May 26 03:50:44 2002 + @author Peter Diener (based on MoL_RK45Add by Ian Hawke) + @desc + Performs a single step of a Runge-Kutta 87 type time + integration, storing the error estimate. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +void MoL_RK87Add(CCTK_ARGUMENTS) +{ + + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + cGroupDynamicData arraydata; + CCTK_INT groupindex, ierr, i, j; + CCTK_INT arraytotalsize, arraydim; + + static CCTK_INT scratchspace_firstindex = -99; + CCTK_INT index, var, scratchstep, alphaindex, scratchindex; + CCTK_INT totalsize; + + CCTK_REAL * restrict UpdateVar; + CCTK_REAL const * restrict RHSVar; + CCTK_REAL * restrict ScratchVar; + CCTK_REAL * restrict ErrorVar; + CCTK_REAL const * restrict OldVar; + + CCTK_INT arrayscratchlocation; + + CCTK_REAL beta, gamma, gamma_error; + + static const CCTK_REAL beta_array[12][12] = { + { 1.0/18.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, + { 1.0/48.0, 1.0/16.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, + { 1.0/32.0, 0.0, 3.0/32.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, + { 5.0/16.0, 0.0, -75.0/64.0, 75.0/64.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, + { 3.0/80.0, 0.0, 0.0, 3.0/16.0, 3.0/20.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, + { 29443841.0/614563906.0, 0.0, 0.0, 77736538.0/692538347.0, -28693883.0/1125000000.0, 23124283.0/1800000000.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, + { 16016141.0/946692911.0, 0.0, 0.0, 61564180.0/158732637.0, 22789713.0/633445777.0, 545815736.0/2771057229.0, -180193667.0/1043307555.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, + { 39632708.0/573591083.0, 0.0, 0.0, -433636366.0/683701615.0, -421739975.0/2616292301.0, 100302831.0/723423059.0, 790204164.0/839813087.0, 800635310.0/3783071287.0, 0.0, 0.0, 0.0, 0.0 }, + { 246121993.0/1340847787.0, 0.0, 0.0, -37695042795.0/15268766246.0, -309121744.0/1061227803.0, -12992083.0/490766935.0, 6005943493.0/2108947869.0, 393006217.0/1396673457.0, 123872331.0/1001029789.0, 0.0, 0.0, 0.0 }, + { -1028468189.0/846180014.0, 0.0, 0.0, 8478235783.0/508512852.0, 1311729495.0/1432422823.0, -10304129995.0/1701304382.0, -48777925059.0/3047939560.0, 15336726248.0/1032824649.0, -45442868181.0/3398467696.0, 3065993473.0/597172653.0, 0.0, 0.0 }, + { 185892177.0/718116043.0, 0.0, 0.0, -3185094517.0/667107341.0, -477755414.0/1098053517.0, -703635378.0/230739211.0, 5731566787.0/1027545527.0, 5232866602.0/850066563.0, -4093664535.0/808688257.0, 3962137247.0/1805957418.0, 65686358.0/487910083.0, 0.0 }, + { 403863854.0/491063109.0, 0.0, 0.0, -5068492393.0/434740067.0, -411421997.0/543043805.0, 652783627.0/914296604.0, 11173962825.0/925320556.0, -13158990841.0/6184727034.0, 3936647629.0/1978049680.0, -160528059.0/685178525.0, 248638103.0/1413531060.0, 0.0 } + }; + + static const CCTK_REAL gamma_array[13] = + { 14005451.0/335480064.0, + 0.0, + 0.0, + 0.0, + 0.0, + -59238493.0/1068277825.0, + 181606767.0/758867731.0, + 561292985.0/797845732.0, + -1041891430.0/1371343529.0, + 760417239.0/1151165299.0, + 118820643.0/751138087.0, + -528747749.0/2220607170.0, + 1.0/4.0 + }; + + static const CCTK_REAL gammastar_array[13] = + { 13451932.0/455176623.0, + 0.0, + 0.0, + 0.0, + 0.0, + -808719846.0/976000145.0, + 1757004468.0/5645159321.0, + 656045339.0/265891186.0, + -3867574721.0/1518517206.0, + 465885868.0/322736535.0, + 53011238.0/667516719.0, + 2.0/45.0, + 0.0 + }; + + totalsize = 1; + for (arraydim = 0; arraydim < cctk_dim; arraydim++) + { + totalsize *= cctk_lsh[arraydim]; + } + + if (scratchspace_firstindex == -99) + { + scratchspace_firstindex = CCTK_FirstVarIndex("MOL::SCRATCHSPACE"); + } + + /* Real GFs */ + + /* First store (dt times) the rhs in the scratch array. */ + + for (var = 0; var < MoLNumEvolvedVariables; var++) + { + const CCTK_REAL tmp = (*Original_Delta_Time) / cctkGH->cctk_timefac; + + UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, + EvolvedVariableIndex[var]); + RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, + RHSVariableIndex[var]); + ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, + CCTK_FirstVarIndex("MOL::SCRATCHSPACE") + + var + + MoL_Num_Evolved_Vars * + (MoL_Intermediate_Steps - + (*MoL_Intermediate_Step))); + for (index = 0; index < totalsize; index++) + { + ScratchVar[index] = tmp * RHSVar[index]; + } + + } + + + for (var = 0; var < MoLNumEvolvedVariables; var++) + { + OldVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 1, + EvolvedVariableIndex[var]); + UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, + EvolvedVariableIndex[var]); + RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, + CCTK_FirstVarIndex("MOL::SCRATCHSPACE") + + var + + MoL_Num_Evolved_Vars * + (MoL_Intermediate_Steps - + (*MoL_Intermediate_Step))); + ErrorVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, + CCTK_FirstVarIndex("MOL::ERRORESTIMATE") + + var); + + if (*MoL_Intermediate_Step - 1) + { + + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] = OldVar[index]; + } + + for (scratchstep = 0; + scratchstep < MoL_Intermediate_Steps - (*MoL_Intermediate_Step) + 1; + scratchstep++) + { + + ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, + CCTK_FirstVarIndex("MOL::SCRATCHSPACE") + + var + + MoL_Num_Evolved_Vars * scratchstep); + + beta = beta_array[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)][scratchstep]; + + if ( (beta > MoL_Tiny)||(beta < -MoL_Tiny) ) + { + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += beta * ScratchVar[index]; + } + } + + } + + } + else + { + + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] = OldVar[index]; + ErrorVar[index] = 0; + } + + for (scratchstep = 0; scratchstep < 13; scratchstep++) + { + + ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, + CCTK_FirstVarIndex("MOL::SCRATCHSPACE") + + var + + MoL_Num_Evolved_Vars * scratchstep); + + gamma = gamma_array[scratchstep]; + gamma_error = gamma - gammastar_array[scratchstep]; + + if ( (gamma > MoL_Tiny)||(gamma < -MoL_Tiny) ) + { + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += gamma * ScratchVar[index]; + ErrorVar[index] += gamma_error * ScratchVar[index]; + } + } + + } + + } + + } + + /* Real arrays */ + + arrayscratchlocation = 0; + + for (var = 0; var < MoLNumEvolvedArrayVariables; var++) + { + + CCTK_WARN(0, "Peter has been too lazy to write the RK87 routine " + "out for array variables. Better send him an email..."); + + } + + return; +} + +/******************************************************************** + ********************* Local Routines ************************* + ********************************************************************/ + diff --git a/src/SetTime.c b/src/SetTime.c index af044f4..1d1596d 100644 --- a/src/SetTime.c +++ b/src/SetTime.c @@ -57,6 +57,35 @@ static const CCTK_REAL alpha_array[6] = { 1.0/2.0, }; +/* RK65 coefficients */ +static const CCTK_REAL alpha_array65[8] = { + 0.0, + 1.0/10.0, + 2.0/9.0, + 3.0/7.0, + 3.0/5.0, + 4.0/5.0, + 1.0, + 1.0 +}; + +/* RK87 coefficients */ + static const CCTK_REAL alpha_array87[13] = { + 0.0, + 1.0/18.0, + 1.0/12.0, + 1.0/8.0, + 5.0/16.0, + 3.0/8.0, + 59.0/400.0, + 93.0/200.0, + 5490023248.0/9719169821.0, + 13.0/20.0, + 1201146811.0/1299019798.0, + 1.0, + 1.0 +}; + /******************************************************************** ********************* External Routines ********************** ********************************************************************/ @@ -216,6 +245,22 @@ int MoL_ResetTime(CCTK_ARGUMENTS) + ((alpha_array[substep] - 1) * (* Original_Delta_Time) / cctkGH->cctk_timefac)); } + else if (CCTK_EQUALS(ODE_Method,"RK65")) + { + const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step); + cctkGH->cctk_time + = ((* Original_Time) + + ((alpha_array65[substep] - 1) + * (* Original_Delta_Time) / cctkGH->cctk_timefac)); + } + else if (CCTK_EQUALS(ODE_Method,"RK87")) + { + const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step); + cctkGH->cctk_time + = ((* Original_Time) + + ((alpha_array87[substep] - 1) + * (* Original_Delta_Time) / cctkGH->cctk_timefac)); + } #ifdef MOLDEBUG printf("MoL has once more reset t (%d): %f.\n", *MoL_Intermediate_Step, cctkGH->cctk_time); @@ -301,6 +346,20 @@ int MoL_ResetDeltaTime(CCTK_ARGUMENTS) = ((alpha_array[substep + 1] - alpha_array[substep]) * (* Original_Delta_Time)); } + else if (CCTK_EQUALS(ODE_Method,"RK65")) + { + const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step); + cctkGH->cctk_delta_time + = ((alpha_array65[substep + 1] - alpha_array65[substep]) + * (* Original_Delta_Time)); + } + else if (CCTK_EQUALS(ODE_Method,"RK87")) + { + const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step); + cctkGH->cctk_delta_time + = ((alpha_array87[substep + 1] - alpha_array87[substep]) + * (* Original_Delta_Time)); + } #ifdef MOLDEBUG printf("MoL has once more reset dt (%d): %f.\n", *MoL_Intermediate_Step, diff --git a/src/StepSize.c b/src/StepSize.c index 03d6fb9..cb5499e 100644 --- a/src/StepSize.c +++ b/src/StepSize.c @@ -213,7 +213,7 @@ void MoL_ReduceAdaptiveError(CCTK_ARGUMENTS) int redop; CCTK_REAL red_local[2], red_global[2]; - + CCTK_REAL p1, p2; int ierr; /* Get global result over all processors */ @@ -235,6 +235,24 @@ void MoL_ReduceAdaptiveError(CCTK_ARGUMENTS) CCTK_VInfo (CCTK_THORNSTRING, "Integration accuracy quotient is %g", (double)*Error); } + if ( CCTK_EQUALS(ODE_Method,"RK45") ) + { + p1 = 0.2; + p2 = 0.25; + } + + if ( CCTK_EQUALS(ODE_Method,"RK65") ) + { + p1 = 1.0/6.0; + p2 = 0.2; + } + + if ( CCTK_EQUALS(ODE_Method,"RK87") ) + { + p1 = 0.125; + p2 = 1.0/7.0; + } + /* Decide whether to accept this step */ *MoL_Stepsize_Bad = *Error > 1; @@ -249,7 +267,7 @@ void MoL_ReduceAdaptiveError(CCTK_ARGUMENTS) } cctkGH->cctk_delta_time - = safety_factor * (*Original_Delta_Time) / pow(*Error, 0.2); + = safety_factor * (*Original_Delta_Time) / pow(*Error, p1); if (! CCTK_EQUALS(verbose, "none")) { CCTK_VInfo (CCTK_THORNSTRING, "Setting time step to %g", (double)cctkGH->cctk_delta_time); @@ -271,7 +289,7 @@ void MoL_ReduceAdaptiveError(CCTK_ARGUMENTS) { /* The error is acceptable; estimate the next step size */ *EstimatedDt = - safety_factor * (*Original_Delta_Time) / pow(*Error, 0.25); + safety_factor * (*Original_Delta_Time) / pow(*Error, p2); if (*EstimatedDt > (*Original_Delta_Time) * maximum_increase) { diff --git a/src/make.code.defn b/src/make.code.defn index 7683c09..f62902e 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -12,6 +12,8 @@ SRCS = ChangeType.c \ RK2.c \ RK3.c \ RK45.c \ + RK65.c \ + RK87.c \ Registration.c \ RKCoefficients.c \ RHSNaNCheck.c \ |