diff options
-rw-r--r-- | src/SetTime.c | 25 |
1 files changed, 25 insertions, 0 deletions
diff --git a/src/SetTime.c b/src/SetTime.c index 94380b1..af044f4 100644 --- a/src/SetTime.c +++ b/src/SetTime.c @@ -47,6 +47,16 @@ int MoL_ResetDeltaTime(CCTK_ARGUMENTS); ********************* Local Data ***************************** ********************************************************************/ +/* RK45 coefficients */ +static const CCTK_REAL alpha_array[6] = { + 0.0, + 1.0/4.0, + 3.0/8.0, + 12.0/13.0, + 1.0, + 1.0/2.0, +}; + /******************************************************************** ********************* External Routines ********************** ********************************************************************/ @@ -198,6 +208,14 @@ int MoL_ResetTime(CCTK_ARGUMENTS) 0.5*(*Original_Delta_Time)/cctkGH->cctk_timefac; } } + else if (CCTK_EQUALS(ODE_Method,"RK45")) + { + const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step); + cctkGH->cctk_time + = ((* Original_Time) + + ((alpha_array[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); @@ -276,6 +294,13 @@ int MoL_ResetDeltaTime(CCTK_ARGUMENTS) cctkGH->cctk_delta_time = 2.0/3.0*(*Original_Delta_Time); } } + else if (CCTK_EQUALS(ODE_Method,"RK45")) + { + const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step); + cctkGH->cctk_delta_time + = ((alpha_array[substep + 1] - alpha_array[substep]) + * (* Original_Delta_Time)); + } #ifdef MOLDEBUG printf("MoL has once more reset dt (%d): %f.\n", *MoL_Intermediate_Step, |