From eeb8bdf88f9b3de750e95bdb98ae78a8f9b4ad49 Mon Sep 17 00:00:00 2001 From: hawke Date: Fri, 28 Jul 2006 09:43:58 +0000 Subject: Efficient RK4, as provided by Yosef Zlochower. git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/MoL/trunk@113 578cdeb0-5ea1-4b81-8215-5a3b8777ee0b --- param.ccl | 1 + schedule.ccl | 30 ++++++ src/IndexArrays.c | 4 + src/ParamCheck.c | 8 ++ src/RK4.c | 296 +++++++++++++++++++++++++++++++++++++++++++++++++++++ src/SetTime.c | 16 +++ src/make.code.defn | 1 + 7 files changed, 356 insertions(+) create mode 100644 src/RK4.c diff --git a/param.ccl b/param.ccl index ea83504..bd8e86f 100644 --- a/param.ccl +++ b/param.ccl @@ -87,6 +87,7 @@ KEYWORD ODE_Method "The ODE method use by MoL to do time integration" "ICN-avg" :: "Iterative Crank Nicholson with averaging" "RK2" :: "Efficient RK2" "RK3" :: "Efficient RK3" + "RK4" :: "Efficient RK4" "RK45" :: "RK45 (Fehlberg) with error estimation" "RK45CK" :: "RK45CK (Cash-Karp) with error estimation" "RK65" :: "RK65 with error estimation" diff --git a/schedule.ccl b/schedule.ccl index aff477d..b54905c 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -558,6 +558,13 @@ else if (CCTK_Equals(ODE_Method,"RK3")) LANG: C } "Updates calculated with the efficient Runge-Kutta 3 method" } +else if (CCTK_Equals(ODE_Method,"RK4")) +{ + schedule MoL_RK4Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep + { + LANG: C + } "Updates calculated with the efficient Runge-Kutta 4 method" +} else if (CCTK_Equals(ODE_Method,"RK45") || CCTK_Equals(ODE_Method,"RK45CK")) { STORAGE: ErrorEstimate ErrorScalars @@ -625,6 +632,29 @@ schedule GROUP MoL_PostStep AT PostRestrict { } "Ensure that everything is correct after restriction" +############################################################ +### Additional boundary condition bins as requested by ### +### Yosef Zlochower. ### +############################################################ + +# schedule GROUP MoL_OldBdry_Wrap IN MoL_PostStep +# { +# } "Wrapper group, do not schedule directly into this group" +# +# schedule GROUP MoL_OldStyleBoundaries in MoL_OldBdry_Wrap +# { +# } "Place old style boundary routines here" +# +# schedule MoL_OldBdry_SetDt IN MoL_OldBdry_Wrap BEFORE MoL_OldStyleBoundaries +# { +# LANGUAGE: C +# } "Store and change dt" +# +# schedule MoL_OldBdry_ResetDt IN MoL_OldBdry_Wrap AFTER MoL_OldStyleBoundaries +# { +# LANGUAGE: C +# } "Reset dt" + ################################################# ### Final internal MoL stuff; decrement the ### ### counter, change time and timestep ### diff --git a/src/IndexArrays.c b/src/IndexArrays.c index 5e59734..e30d169 100644 --- a/src/IndexArrays.c +++ b/src/IndexArrays.c @@ -304,6 +304,10 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS) { sprintf(infoline, "Runge-Kutta 3"); } + else if (CCTK_EQUALS(ODE_Method,"RK4")) + { + sprintf(infoline, "Runge-Kutta 4"); + } else if (CCTK_EQUALS(ODE_Method,"RK45")) { sprintf(infoline, "Runge-Kutta 45 (Fehlberg)"); diff --git a/src/ParamCheck.c b/src/ParamCheck.c index c2346b0..112c530 100644 --- a/src/ParamCheck.c +++ b/src/ParamCheck.c @@ -128,6 +128,14 @@ int MoL_ParamCheck(CCTK_ARGUMENTS) "number of intermediate steps must be 3"); } + if ( (CCTK_Equals(ODE_Method, "RK4"))&&(!(MoL_Intermediate_Steps == 4)) + && (!(MoL_Num_Scratch_Levels > 1))) + { + CCTK_PARAMWARN("When using the efficient RK4 evolver the " + "number of intermediate steps must be 4, and" + " the number of scratch levels at least 1"); + } + if ( (CCTK_Equals(ODE_Method, "RK45") || CCTK_Equals(ODE_Method, "RK45CK")) && ( !((MoL_Intermediate_Steps == 6)||(MoL_Num_Scratch_Levels > 5)) ) ) { diff --git a/src/RK4.c b/src/RK4.c new file mode 100644 index 0000000..940402b --- /dev/null +++ b/src/RK4.c @@ -0,0 +1,296 @@ + /*@@ + @file RK4.c + @date Fri July 14, 2006 + @author Yosef Zlochower + @desc + A routine to perform RK4 evolution. Mostly copied from + genericRK.c + @enddesc + @version $Header$ + @@*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include +#include "ExternalVariables.h" + +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusBase_MoL_RK4_c); + +/******************************************************************** + ********************* Local Data Types *********************** + ********************************************************************/ + +/******************************************************************** + ********************* Local Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ***************** Scheduled Routine Prototypes ********************* + ********************************************************************/ + +void MoL_RK4Add(CCTK_ARGUMENTS); + +/******************************************************************** + ********************* Other Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ********************* Local Data ***************************** + ********************************************************************/ + +/******************************************************************** + ********************* External Routines ********************** + ********************************************************************/ + + /*@@ + @routine MoL_RK4Add + @date + @author + @desc + Performs a single step of a RK4 type time + integration. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +void MoL_RK4Add(CCTK_ARGUMENTS) +{ + + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + cGroupDynamicData arraydata; + CCTK_INT groupindex, ierr; + CCTK_INT arraytotalsize, arraydim; + + /* FIXME */ + +#ifdef MOLDOESCOMPLEX + +CCTK_WARN(0, "not implemented"); + CCTK_COMPLEX Complex_alpha, Complex_beta, Complex_Delta_Time; + CCTK_COMPLEX * restrict OldComplexVar; + CCTK_COMPLEX * restrict UpdateComplexVar; + CCTK_COMPLEX const * restrict RHSComplexVar; + CCTK_COMPLEX * restrict ScratchComplexVar; + +#endif + + static CCTK_INT scratchspace_firstindex = -99; + CCTK_INT index, var, scratchstep, alphaindex, scratchindex; + CCTK_INT totalsize, singlearraysize; + CCTK_REAL alpha, beta; + CCTK_REAL * restrict UpdateVar; + CCTK_REAL * restrict OldVar; + CCTK_REAL const * restrict RHSVar; + CCTK_REAL * restrict ScratchVar; + + CCTK_INT arrayscratchlocation; + + totalsize = 1; + for (arraydim = 0; arraydim < cctk_dim; arraydim++) + { + totalsize *= cctk_lsh[arraydim]; + } + + if (scratchspace_firstindex == -99) + { + scratchspace_firstindex = CCTK_FirstVarIndex("MOL::SCRATCHSPACE"); + } + + switch (MoL_Intermediate_Steps - (*MoL_Intermediate_Step)) + { + case 0: + alpha = 1.0 / 3.0; + beta = 0.5; + break; + case 1: + alpha = 2.0 / 3.0; + beta = 0.5; + break; + case 2: + alpha = 1.0 / 3.0; + beta = 1.0; + break; + case 3: + alpha = 1.0; + beta = 1.0 / 6.0; + } + + /* FIXME */ + +#ifdef MOLDOESCOMPLEX + + Complex_Delta_Time = CCTK_Cmplx((*Original_Delta_Time) / cctkGH->cctk_timefac, 0); + Complex_beta = CCTK_Cmplx(beta, 0); + +#endif + + /* Real GFs */ + + for (var = 0; var < MoLNumEvolvedVariables; var++) + { + + UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, + EvolvedVariableIndex[var]); + OldVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 1, + EvolvedVariableIndex[var]); + RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, + RHSVariableIndex[var]); +/* #define MOLDEBUG 1 */ +#ifdef MOLDEBUG + printf("In generic RK. Variable %d (%s). RHS %d (%s). beta %g.\n", + EvolvedVariableIndex[var], + CCTK_VarName(EvolvedVariableIndex[var]), + RHSVariableIndex[var], + CCTK_VarName(RHSVariableIndex[var]), + beta); +#endif + + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] = OldVar[index] + + (*Original_Delta_Time) / cctkGH->cctk_timefac * beta * RHSVar[index]; +#ifdef MOLDEBUG + if (CCTK_EQUALS(verbose,"extreme")) + { + printf("Variable: %d. Index: %d. dt: %f. beta %f. RHS: %f. q: %f.\n", + var, index, (*Original_Delta_Time) / cctkGH->cctk_timefac, beta, RHSVar[index], + UpdateVar[index]); + } +#endif + } + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex + + var ); + + /* scratch storage */ + if ((*MoL_Intermediate_Step) == MoL_Intermediate_Steps) + { + for (index = 0; index < totalsize; index++) + { + ScratchVar[index] = 0; + } + } + + if ((*MoL_Intermediate_Step)>1) + { + for (index = 0; index < totalsize; index++) + { + ScratchVar[index] += alpha * UpdateVar[index]; + } + } + else + { + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += ScratchVar[index] - 4.0 / 3.0 * OldVar[index]; + } + } + + } + + + /* Complex GFs */ + + /* FIXME */ + +#ifdef MOLDOESCOMPLEX +CCTK_WARN(0, "not done"); +#endif + + /* Real arrays */ + + arrayscratchlocation = 0; + + arraytotalsize = MoL_Max_Evolved_Array_Size; + if (MoL_Num_Scratch_Levels) + { + groupindex = CCTK_GroupIndex("MOL::ARRAYSCRATCHSPACE"); + ierr = CCTK_GroupDynamicData(cctkGH, groupindex, + &arraydata); + if (ierr) + { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "The driver does not return group information " + "for group '%s' (%d).", + CCTK_GroupName(groupindex), ierr); + } + arraytotalsize = 1; + for (arraydim = 0; arraydim < arraydata.dim; arraydim++) + { + arraytotalsize *= arraydata.lsh[arraydim]; + } + } + singlearraysize = arraytotalsize / MoLNumEvolvedArrayVariables; + + for (var = 0; var < MoLNumEvolvedArrayVariables; var++) + { + + UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, + EvolvedArrayVariableIndex[var]); + OldVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 1, + EvolvedArrayVariableIndex[var]); + RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, + RHSArrayVariableIndex[var]); + + groupindex = CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]); + ierr = CCTK_GroupDynamicData(cctkGH, groupindex, + &arraydata); + if (ierr) + { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "The driver does not return group information " + "for group '%s'.", + CCTK_GroupName(groupindex)); + } + arraytotalsize = 1; + for (arraydim = 0; arraydim < arraydata.dim; arraydim++) + { + arraytotalsize *= arraydata.lsh[arraydim]; + } + + ScratchVar = &ArrayScratchSpace[arrayscratchlocation]; + + for (index = 0; index < arraytotalsize; index++) + { + UpdateVar[index] = OldVar[index] + + (*Original_Delta_Time) / cctkGH->cctk_timefac * beta * RHSVar[index]; + } + + if ((*MoL_Intermediate_Step) == MoL_Intermediate_Steps) + { + for (index = 0; index < arraytotalsize; index++) + { + ScratchVar[index] = 0; + } + } + + if ((*MoL_Intermediate_Step)>1) + { + for (index = 0; index < arraytotalsize; index++) + { + ScratchVar[index] += alpha * UpdateVar[index]; + } + } + else + { + for (index = 0; index < arraytotalsize; index++) + { + UpdateVar[index] += ScratchVar[index] - 4.0 / 3.0 * OldVar[index]; + } + } + arrayscratchlocation += arraytotalsize; + } + + return; +} diff --git a/src/SetTime.c b/src/SetTime.c index 3e9ae3f..5149eb2 100644 --- a/src/SetTime.c +++ b/src/SetTime.c @@ -247,6 +247,22 @@ int MoL_ResetTime(CCTK_ARGUMENTS) 0.5*(*Original_Delta_Time)/cctkGH->cctk_timefac; } } + else if (CCTK_EQUALS(ODE_Method,"RK4")) + { + const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step); + + CCTK_REAL dt = (*Original_Delta_Time)/cctkGH->cctk_timefac; + switch (substep) + { + case 1: + case 2: + dt *= 0.5; + break; + default: + dt = 0; + } + cctkGH->cctk_time = (*Original_Time) - dt; + } else if (CCTK_EQUALS(ODE_Method,"RK45") || CCTK_EQUALS(ODE_Method,"RK45CK")) { const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step); diff --git a/src/make.code.defn b/src/make.code.defn index f62902e..926dbbb 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -11,6 +11,7 @@ SRCS = ChangeType.c \ ParamCheck.c \ RK2.c \ RK3.c \ + RK4.c \ RK45.c \ RK65.c \ RK87.c \ -- cgit v1.2.3