From 26067015fd74e5079d09a589fc3bab9e7fd13f22 Mon Sep 17 00:00:00 2001 From: rhaas Date: Thu, 2 Aug 2012 16:34:52 +0000 Subject: MoL: add Multirate capabilities. This add three new multirate RK schemes to MoL. Flags indicate whether it is time to execute slow RHS computation. For instance, in the RK4-RK2 scheme, there are 4 substeps in total, but the RK2 RHS are only evaluated in the very first and in the very last step of the four substeps. From: Christian Reisswig, minor changes by Roland Haas git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/MoL/trunk@175 578cdeb0-5ea1-4b81-8215-5a3b8777ee0b --- interface.ccl | 27 +++ param.ccl | 9 + schedule.ccl | 30 ++- src/ChangeType.c | 409 +++++++++++++++++++++++++++++++++++++++ src/Counter.c | 73 ++++++- src/ExternalVariables.h | 3 + src/IndexArrays.c | 47 +++++ src/MoL.h | 3 + src/MoLFunctions.h | 3 + src/ParamCheck.c | 18 ++ src/RK2-MR-2_1.c | 358 ++++++++++++++++++++++++++++++++++ src/RK4-MR-2_1.c | 502 ++++++++++++++++++++++++++++++++++++++++++++++++ src/RK4-RK2.c | 342 +++++++++++++++++++++++++++++++++ src/Registration.c | 446 +++++++++++++++++++++++++++++++++++++++++- src/SetTime.c | 57 ++++++ src/Startup.c | 3 + src/StepSize.c | 4 + src/make.code.defn | 7 +- 18 files changed, 2336 insertions(+), 5 deletions(-) create mode 100644 src/RK2-MR-2_1.c create mode 100644 src/RK4-MR-2_1.c create mode 100644 src/RK4-RK2.c diff --git a/interface.ccl b/interface.ccl index fa72e0c..d2d14f9 100644 --- a/interface.ccl +++ b/interface.ccl @@ -34,14 +34,21 @@ USES FUNCTION EnableProlongating CCTK_INT FUNCTION MoLRegisterEvolved(CCTK_INT IN EvolvedIndex, \ CCTK_INT IN RHSIndex) +CCTK_INT FUNCTION MoLRegisterEvolvedSlow(CCTK_INT IN EvolvedIndex, \ + CCTK_INT IN RHSIndexSlow) + CCTK_INT FUNCTION MoLRegisterConstrained(CCTK_INT IN ConstrainedIndex) CCTK_INT FUNCTION MoLRegisterSaveAndRestore(CCTK_INT IN SandRIndex) CCTK_INT FUNCTION MoLRegisterEvolvedGroup(CCTK_INT IN EvolvedIndex, \ CCTK_INT IN RHSIndex) +CCTK_INT FUNCTION MoLRegisterEvolvedGroupSlow(CCTK_INT IN EvolvedIndex, \ + CCTK_INT IN RHSIndexSlow) CCTK_INT FUNCTION MoLRegisterConstrainedGroup(CCTK_INT IN ConstrainedIndex) CCTK_INT FUNCTION MoLRegisterSaveAndRestoreGroup(CCTK_INT IN SandRIndex) CCTK_INT FUNCTION MoLChangeToEvolved(CCTK_INT IN EvolvedIndex, \ CCTK_INT IN RHSIndex) +CCTK_INT FUNCTION MoLChangeToEvolvedSlow(CCTK_INT IN EvolvedIndex, \ + CCTK_INT IN RHSIndexSlow) CCTK_INT FUNCTION MoLChangeToConstrained(CCTK_INT IN ConstrainedIndex) CCTK_INT FUNCTION MoLChangeToSaveAndRestore(CCTK_INT IN SandRIndex) CCTK_INT FUNCTION MoLChangeToNone(CCTK_INT IN RemoveIndex) @@ -49,17 +56,21 @@ CCTK_INT FUNCTION MoLQueryEvolvedRHS(CCTK_INT IN EvolvedIndex) CCTK_INT FUNCTION MoLNumIntegratorSubsteps() PROVIDES FUNCTION MoLRegisterEvolved WITH MoL_RegisterEvolved LANGUAGE C +PROVIDES FUNCTION MoLRegisterEvolvedSlow WITH MoL_RegisterEvolvedSlow LANGUAGE C PROVIDES FUNCTION MoLRegisterConstrained WITH MoL_RegisterConstrained \ LANGUAGE C PROVIDES FUNCTION MoLRegisterSaveAndRestore WITH MoL_RegisterSaveAndRestore \ LANGUAGE C PROVIDES FUNCTION MoLRegisterEvolvedGroup WITH MoL_RegisterEvolvedGroup \ LANGUAGE C +PROVIDES FUNCTION MoLRegisterEvolvedGroupSlow WITH MoL_RegisterEvolvedGroupSlow \ + LANGUAGE C PROVIDES FUNCTION MoLRegisterConstrainedGroup WITH \ MoL_RegisterConstrainedGroup LANGUAGE C PROVIDES FUNCTION MoLRegisterSaveAndRestoreGroup WITH \ MoL_RegisterSaveAndRestoreGroup LANGUAGE C PROVIDES FUNCTION MoLChangeToEvolved WITH MoL_ChangeToEvolved LANGUAGE C +PROVIDES FUNCTION MoLChangeToEvolvedSlow WITH MoL_ChangeToEvolvedSlow LANGUAGE C PROVIDES FUNCTION MoLChangeToConstrained WITH MoL_ChangeToConstrained \ LANGUAGE C PROVIDES FUNCTION MoLChangeToSaveAndRestore WITH MoL_ChangeToSaveAndRestore \ @@ -226,6 +237,17 @@ CCTK_INT MoL_Counters \ { MoL_Intermediate_Step MoL_Stepsize_Bad + + # A flag indicating whether it is time for slow RHS evaluation. + # Oustide the MoL loop, it is guaranteed to be 1. + # It is only zero for certain MoL substeps when multirate methods are used. + MoL_SlowStep + + # A flag indicating whether it is time for slow post step computations (e.g. applying BCs) + # Oustide the MoL loop, it is guaranteed to be 1. + # It is only zero for certain MoL substeps when multirate methods are used. + MoL_SlowPostStep + } "The counter for the time integration method" CCTK_REAL MoL_Original_Time \ @@ -241,6 +263,11 @@ CCTK_REAL ScratchSpace[MoL_Num_Evolved_Vars*MoL_Num_Scratch_Levels] \ Timelevels = 1 \ TAGS = 'Prolongation="None" Checkpoint="no"' +CCTK_REAL ScratchSpaceSlow[MoL_Num_Evolved_Vars_Slow*MoL_Num_Scratch_Levels] \ + TYPE = GF \ + Timelevels = 1 \ + TAGS = 'Prolongation="None" Checkpoint="no"' + CCTK_REAL SandRScratchSpace[MoL_Num_SaveAndRestore_Vars] \ TYPE = GF \ Timelevels = 1 \ diff --git a/param.ccl b/param.ccl index ba43b51..e82dd07 100644 --- a/param.ccl +++ b/param.ccl @@ -8,6 +8,12 @@ CCTK_INT MoL_Num_Evolved_Vars "The maximum number of variables to be evolved by 0:* :: "Anything non negative. Added to by other thorns." } 0 +CCTK_INT MoL_Num_Evolved_Vars_Slow "The maximum number of 'slow' variables to be evolved by MoL" ACCUMULATOR = (x+y) +{ + 0:* :: "Anything non negative. Added to by other thorns." +} 0 + + CCTK_INT MoL_Num_Constrained_Vars "The maximum number of constrained variables with timelevels that MoL needs to know about" ACCUMULATOR = (x+y) { 0:* :: "Anything non negative. Added to by other thorns." @@ -92,6 +98,9 @@ KEYWORD ODE_Method "The ODE method use by MoL to do time integration" "RK45CK" :: "RK45CK (Cash-Karp) with error estimation" "RK65" :: "RK65 with error estimation" "RK87" :: "RK87 with error estimation" + "RK2-MR-2:1" :: "2nd order 2:1 multirate RK scheme based on RK2 due to Schlegel et al 2009" + "RK4-MR-2:1" :: "3rd order 2:1 multirate RK scheme based on RK43 due to Schlegel et al 2009" + "RK4-RK2" :: "RK4 as fast method and RK2 as slow method" } "ICN" KEYWORD Generic_Type "If using the generic method, which sort" diff --git a/schedule.ccl b/schedule.ccl index f0aed52..8dc65f9 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -399,7 +399,13 @@ else } } } - + + +if (MoL_Num_Evolved_Vars_Slow) +{ + STORAGE: ScratchSpaceSlow +} + ###################################################### ### StartStep contains the routines that just do ### ### internal MoL stuff; setting the counter, and ### @@ -607,6 +613,28 @@ else if (CCTK_Equals(ODE_Method,"ICN-avg")) LANG: C } "Updates calculated with the averaging ICN method" } +else if (CCTK_Equals(ODE_Method,"RK2-MR-2:1")) +{ + schedule MoL_RK2_MR_2_1_Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep + { + LANG: C + } "Updates calculated with the multirate Runge-Kutta 2 method" +} +else if (CCTK_Equals(ODE_Method,"RK4-MR-2:1")) +{ + schedule MoL_RK4_MR_2_1_Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep + { + LANG: C + } "Updates calculated with the multirate Runge-Kutta 4 method" +} +else if (CCTK_Equals(ODE_Method,"RK4-RK2")) +{ + schedule MoL_RK4_RK2_Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep + { + LANG: C + } "Updates calculated with the multirate RK4/RK2 method" +} + ################################################## ### Physics thorns can apply boundaries and ### diff --git a/src/ChangeType.c b/src/ChangeType.c index e0bc333..9c1db4a 100644 --- a/src/ChangeType.c +++ b/src/ChangeType.c @@ -27,6 +27,7 @@ CCTK_FILEVERSION(CactusBase_MoL_ChangeType_c); #define MOL_EVOLVED_VARTYPE 1 #define MOL_CONSTRAINED_VARTYPE 2 #define MOL_SANDR_VARTYPE 3 +#define MOL_EVOLVEDSLOW_VARTYPE 4 /******************************************************************** ********************* Local Routine Prototypes ********************* @@ -42,6 +43,8 @@ CCTK_FILEVERSION(CactusBase_MoL_ChangeType_c); CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex); +CCTK_INT MoL_ChangeToEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndexSlow); + CCTK_INT MoL_ChangeToConstrained(CCTK_INT ConstrainedIndex); CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex); @@ -95,6 +98,16 @@ CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex) } } + for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!vartype); index++) + { + vartype = MOL_EVOLVEDSLOW_VARTYPE * + (EvolvedVariableIndexSlow[index] == EvolvedIndex); + if (vartype) + { + usedindex = index; + } + } + for (index = 0; (index < MoLNumConstrainedVariables)&&(!vartype); index++) { vartype = MOL_CONSTRAINED_VARTYPE * @@ -174,6 +187,46 @@ CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex) break; } + case MOL_EVOLVEDSLOW_VARTYPE: + + { + timelevs = CCTK_MaxTimeLevelsVI(RHSIndex); + if (timelevs < 1) { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "Warning whilst trying to change variable " + "index %i to evolved type from constrained.", RHSIndex); + CCTK_WARN(0, "The RHS index passed does not correspond to a GF."); + } + + if (!(MoLNumEvolvedVariables < MoL_Num_Evolved_Vars)) + { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "Warning whilst trying to change variable " + "index %i (%s) to evolved.", + EvolvedIndex, CCTK_VarName(EvolvedIndex)); + CCTK_WARN(0, "When changing type there are more variables " + "than the accumulator parameter " + "MoL_Num_Evolved_Vars allows. Check that " + "you are accumulating onto this parameter correctly"); + } + + for (index = usedindex; index < MoLNumEvolvedVariablesSlow - 1; + index++) + { + EvolvedVariableIndexSlow[index] = EvolvedVariableIndexSlow[index+1]; + } + MoLNumEvolvedVariablesSlow--; + EvolvedVariableIndex[MoLNumEvolvedVariables] = EvolvedIndex; + RHSVariableIndex[MoLNumEvolvedVariables] = RHSIndex; + MoLNumEvolvedVariables++; +#ifdef MOLDEBUG + printf("Changing (constrained): vars %d var %d (%s).\n", + MoLNumEvolvedVariables, EvolvedIndex, + CCTK_VarName(EvolvedIndex)); +#endif + break; + } + case MOL_CONSTRAINED_VARTYPE: { @@ -265,6 +318,266 @@ CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex) return 0; +} + + /*@@ + @routine MoL_ChangeToEvolvedSlow + @date Thu May 30 16:45:30 2002 + @author Ian Hawke, Roland Haas + @desc + Changes a variable to evolved slow type. Checks to see which type it was + before and does the bookkeeping on the index arrays. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +CCTK_INT MoL_ChangeToEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex) +{ + + DECLARE_CCTK_PARAMETERS; + + CCTK_INT index, usedindex; + CCTK_INT vartype; /* See the defines at the top of file */ + CCTK_INT timelevs; + + vartype = 0; + usedindex = -1; + + for (index = 0; (index < MoLNumEvolvedVariables)&&(!vartype); index++) + { + vartype = MOL_EVOLVED_VARTYPE * + (EvolvedVariableIndex[index] == EvolvedIndex); + if (vartype) + { + usedindex = index; + } + } + + for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!vartype); index++) + { + vartype = MOL_EVOLVEDSLOW_VARTYPE * + (EvolvedVariableIndexSlow[index] == EvolvedIndex); + if (vartype) + { + usedindex = index; + } + } + + for (index = 0; (index < MoLNumConstrainedVariables)&&(!vartype); index++) + { + vartype = MOL_CONSTRAINED_VARTYPE * + (ConstrainedVariableIndex[index] == EvolvedIndex); + if (vartype) + { + usedindex = index; + } + } + + for (index = 0; (index < MoLNumSandRVariables)&&(!vartype); index++) + { + vartype = MOL_SANDR_VARTYPE * + (SandRVariableIndex[index] == EvolvedIndex); + if (vartype) + { + usedindex = index; + } + } + + switch (vartype) + { + + case MOL_UNKNOWN_VARTYPE: + + { + timelevs = CCTK_MaxTimeLevelsVI(EvolvedIndex); + if (timelevs < 1) + { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "Warning whilst trying to change variable " + "index %i to slow evolved.", EvolvedIndex); + CCTK_WARN(0, "The index passed does not correspond to a GF."); + } + else if (timelevs == 1) + { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "Warning whilst trying to change variable " + "index %i to slow evolved.", EvolvedIndex); + CCTK_WARN(0, "The index passed only has a single timelevel."); + } + timelevs = CCTK_MaxTimeLevelsVI(RHSIndex); + if (timelevs < 1) { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "Warning whilst trying to change variable " + "index %i to slow evolved (RHS GF).", RHSIndex); + CCTK_WARN(0, "The index passed does not correspond to a GF."); + } + + if (!(MoLNumEvolvedVariablesSlow < MoL_Num_Evolved_Vars_Slow)) + { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "Warning whilst trying to change variable " + "index %i (%s) to slow evolved.", + EvolvedIndex, CCTK_VarName(EvolvedIndex)); + CCTK_WARN(0, "When changing type there are more variables " + "than the accumulator parameter " + "MoL_Num_Evolved_Vars_Slow allows. Check that " + "you are accumulating onto this parameter correctly"); + } + + EvolvedVariableIndexSlow[MoLNumEvolvedVariablesSlow] = EvolvedIndex; + RHSVariableIndexSlow[MoLNumEvolvedVariablesSlow] = RHSIndex; + MoLNumEvolvedVariablesSlow++; +#ifdef MOLDEBUG + printf("Changing (unknown): vars %d var %d (%s).\n", + MoLNumEvolvedVariablesSlow, EvolvedIndex, + CCTK_VarName(EvolvedIndex)); +#endif + break; + } + + case MOL_EVOLVED_VARTYPE: + + { + timelevs = CCTK_MaxTimeLevelsVI(RHSIndex); + if (timelevs < 1) { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "Warning whilst trying to change variable " + "index %i to slow evolved type from constrained.", RHSIndex); + CCTK_WARN(0, "The RHS index passed does not correspond to a GF."); + } + + if (!(MoLNumEvolvedVariablesSlow < MoL_Num_Evolved_Vars_Slow)) + { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "Warning whilst trying to change variable " + "index %i (%s) to slow evolved.", + EvolvedIndex, CCTK_VarName(EvolvedIndex)); + CCTK_WARN(0, "When changing type there are more variables " + "than the accumulator parameter " + "MoL_Num_Evolved_Vars_Slow allows. Check that " + "you are accumulating onto this parameter correctly"); + } + + for (index = usedindex; index < MoLNumEvolvedVariables - 1; + index++) + { + EvolvedVariableIndex[index] = EvolvedVariableIndex[index+1]; + } + MoLNumEvolvedVariables--; + EvolvedVariableIndexSlow[MoLNumEvolvedVariablesSlow] = EvolvedIndex; + RHSVariableIndexSlow[MoLNumEvolvedVariablesSlow] = RHSIndex; + MoLNumEvolvedVariablesSlow++; +#ifdef MOLDEBUG + printf("Changing (constrained): vars %d var %d (%s).\n", + MoLNumEvolvedVariablesSlow, EvolvedIndex, + CCTK_VarName(EvolvedIndex)); +#endif + break; + } + + case MOL_EVOLVEDSLOW_VARTYPE: + + { + RHSVariableIndex[usedindex] = RHSIndex; + break; + } + + case MOL_CONSTRAINED_VARTYPE: + + { + timelevs = CCTK_MaxTimeLevelsVI(RHSIndex); + if (timelevs < 1) { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "Warning whilst trying to change variable " + "index %i to slow evolved type from constrained.", RHSIndex); + CCTK_WARN(0, "The RHS index passed does not correspond to a GF."); + } + + if (!(MoLNumEvolvedVariablesSlow < MoL_Num_Evolved_Vars_Slow)) + { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "Warning whilst trying to change variable " + "index %i (%s) to slow evolved.", + EvolvedIndex, CCTK_VarName(EvolvedIndex)); + CCTK_WARN(0, "When changing type there are more variables " + "than the accumulator parameter " + "MoL_Num_Evolved_Vars_Slow allows. Check that " + "you are accumulating onto this parameter correctly"); + } + + for (index = usedindex; index < MoLNumConstrainedVariables - 1; + index++) + { + ConstrainedVariableIndex[index] = ConstrainedVariableIndex[index+1]; + } + MoLNumConstrainedVariables--; + EvolvedVariableIndexSlow[MoLNumEvolvedVariablesSlow] = EvolvedIndex; + RHSVariableIndexSlow[MoLNumEvolvedVariablesSlow] = RHSIndex; + MoLNumEvolvedVariablesSlow++; +#ifdef MOLDEBUG + printf("Changing (constrained): vars %d var %d (%s).\n", + MoLNumEvolvedVariables, EvolvedIndex, + CCTK_VarName(EvolvedIndex)); +#endif + break; + } + + case MOL_SANDR_VARTYPE: + + { + timelevs = CCTK_MaxTimeLevelsVI(RHSIndex); + if (timelevs < 1) { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "Warning whilst trying to change variable " + "index %i to slow evolved type from save and " + "restore.", RHSIndex); + CCTK_WARN(0, "The RHS index passed does not correspond to a GF."); + } + + if (!(MoLNumEvolvedVariablesSlow < MoL_Num_Evolved_Vars_Slow)) + { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "Warning whilst trying to change variable " + "index %i (%s) to slow evolved.", + EvolvedIndex, CCTK_VarName(EvolvedIndex)); + CCTK_WARN(0, "When changing type there are more variables " + "than the accumulator parameter " + "MoL_Num_Evolved_Vars_Slow allows. Check that " + "you are accumulating onto this parameter correctly"); + } + + for (index = usedindex; index < MoLNumSandRVariables - 1; index++) + { + SandRVariableIndex[index] = SandRVariableIndex[index+1]; + } + MoLNumSandRVariables--; + EvolvedVariableIndexSlow[MoLNumEvolvedVariablesSlow] = EvolvedIndex; + RHSVariableIndexSlow[MoLNumEvolvedVariablesSlow] = RHSIndex; + MoLNumEvolvedVariablesSlow++; +#ifdef MOLDEBUG + printf("Changing (SandR): vars %d var %d (%s).\n", + MoLNumEvolvedVariablesSlow, EvolvedIndex, + CCTK_VarName(EvolvedIndex)); +#endif + break; + } + + default: + + { + CCTK_WARN(0, "Something is seriously wrong in ChangeType.c! " + "Case out of range in switch statement."); + } + + } + + return 0; + } /*@@ @@ -321,6 +634,16 @@ CCTK_INT MoL_ChangeToConstrained(CCTK_INT ConstrainedIndex) } } + for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!vartype); index++) + { + vartype = MOL_EVOLVEDSLOW_VARTYPE * + (EvolvedVariableIndexSlow[index] == ConstrainedIndex); + if (vartype) + { + usedindex = index; + } + } + for (index = 0; (index < MoLNumConstrainedVariables)&&(!vartype); index++) { vartype = MOL_CONSTRAINED_VARTYPE * @@ -397,6 +720,34 @@ CCTK_INT MoL_ChangeToConstrained(CCTK_INT ConstrainedIndex) break; } + case MOL_EVOLVEDSLOW_VARTYPE: + + { + for (index = usedindex; index < MoLNumEvolvedVariablesSlow - 1; index++) + { + EvolvedVariableIndex[index] = EvolvedVariableIndex[index+1]; + RHSVariableIndex[index] = RHSVariableIndex[index+1]; + } + MoLNumEvolvedVariablesSlow--; + if (timelevs > 1) + { + if ( !(MoLNumConstrainedVariables < MoL_Num_Constrained_Vars)) + { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "Warning whilst trying to change variable " + "index %i (%s) to constrained.", + ConstrainedIndex, CCTK_VarName(ConstrainedIndex)); + CCTK_WARN(0, "When changing type there are more variables " + "than the accumulator parameter " + "MoL_Num_Constrained_Vars allows. Check that " + "you are accumulating onto this parameter correctly"); + } + ConstrainedVariableIndex[MoLNumConstrainedVariables] = ConstrainedIndex; + MoLNumConstrainedVariables++; + } + break; + } + case MOL_CONSTRAINED_VARTYPE: { @@ -489,6 +840,16 @@ CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex) } } + for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!vartype); index++) + { + vartype = MOL_EVOLVEDSLOW_VARTYPE * + (EvolvedVariableIndexSlow[index] == SandRIndex); + if (vartype) + { + usedindex = index; + } + } + for (index = 0; (index < MoLNumConstrainedVariables)&&(!vartype); index++) { vartype = MOL_CONSTRAINED_VARTYPE * @@ -558,6 +919,32 @@ CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex) break; } + case MOL_EVOLVEDSLOW_VARTYPE: + + { + + if (!(MoLNumSandRVariables < MoL_Num_SaveAndRestore_Vars)) + { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "Warning whilst trying to change variable " + "index %i (%s) to save and restore.", + SandRIndex, CCTK_VarName(SandRIndex)); + CCTK_WARN(0, "When changing type there are more variables " + "than the accumulator parameter " + "MoL_Num_SaveAndRestore_Vars allows. Check that " + "you are accumulating onto this parameter correctly"); + } + for (index = usedindex; index < MoLNumEvolvedVariablesSlow - 1; index++) + { + EvolvedVariableIndex[index] = EvolvedVariableIndex[index+1]; + RHSVariableIndex[index] = RHSVariableIndex[index+1]; + } + MoLNumEvolvedVariablesSlow--; + SandRVariableIndex[MoLNumSandRVariables] = SandRIndex; + MoLNumSandRVariables++; + break; + } + case MOL_CONSTRAINED_VARTYPE: { @@ -639,6 +1026,16 @@ CCTK_INT MoL_ChangeToNone(CCTK_INT RemoveIndex) } } + for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!vartype); index++) + { + vartype = MOL_EVOLVEDSLOW_VARTYPE * + (EvolvedVariableIndexSlow[index] == RemoveIndex); + if (vartype) + { + usedindex = index; + } + } + for (index = 0; (index < MoLNumConstrainedVariables)&&(!vartype); index++) { vartype = MOL_CONSTRAINED_VARTYPE * @@ -680,6 +1077,18 @@ CCTK_INT MoL_ChangeToNone(CCTK_INT RemoveIndex) break; } + case MOL_EVOLVEDSLOW_VARTYPE: + + { + for (index = usedindex; index < MoLNumEvolvedVariablesSlow - 1; index++) + { + EvolvedVariableIndex[index] = EvolvedVariableIndex[index+1]; + RHSVariableIndex[index] = RHSVariableIndex[index+1]; + } + MoLNumEvolvedVariablesSlow--; + break; + } + case MOL_CONSTRAINED_VARTYPE: { diff --git a/src/Counter.c b/src/Counter.c index 44515de..0093f1f 100644 --- a/src/Counter.c +++ b/src/Counter.c @@ -69,7 +69,13 @@ void MoL_SetCounter(CCTK_ARGUMENTS) *MoL_Intermediate_Step = MoL_Intermediate_Steps; - + // Always indicate execution of slow steps. + // If multirate methods are used, they are set to zero + // for certain substeps. + // In any case, outside the MoL loop, these scalars are guaranteed to be 1 + *MoL_SlowStep = 1; + *MoL_SlowPostStep = 1; + /* #ifdef HAVE_CARPET */ if ((*MoL_Intermediate_Step)) { @@ -114,6 +120,71 @@ void MoL_DecrementCounter(CCTK_ARGUMENTS) (*MoL_Intermediate_Step) --; + + if (CCTK_EQUALS(ODE_Method, "RK2-MR-2:1")) + { + switch (MoL_Intermediate_Steps - *MoL_Intermediate_Step) + { + case 2: + case 4: + if (*MoL_SlowStep) + *MoL_SlowPostStep = 1; + else + *MoL_SlowPostStep = 0; + // don't do slow step + *MoL_SlowStep = 0; + break; + default: + if (!(*MoL_SlowStep)) + *MoL_SlowPostStep = 0; + else + *MoL_SlowPostStep = 1; + // do a slow step! + *MoL_SlowStep = 1; + } + } + + if (CCTK_EQUALS(ODE_Method, "RK4-MR-2:1")) + { + switch (MoL_Intermediate_Steps - *MoL_Intermediate_Step) + { + case 2: + case 4: + case 7: + case 9: + // don't do slow step + *MoL_SlowStep = 0; + break; + default: + // do a slow step! + *MoL_SlowStep = 1; + } + } + + if (CCTK_EQUALS(ODE_Method, "RK4-RK2")) + { + switch (MoL_Intermediate_Steps - *MoL_Intermediate_Step) + { + case 0: + // do a slow step! + *MoL_SlowStep = 1; + // don't sync! + *MoL_SlowPostStep = 0; + break; + case 1: + case 2: + // don't do slow step + *MoL_SlowStep = 0; + *MoL_SlowPostStep = 0; + break; + case 3: + // do a slow step! + *MoL_SlowStep = 1; + // sync! + *MoL_SlowPostStep = 1; + } + } + /* #ifdef HAVE_CARPET */ if (! (*MoL_Intermediate_Step)) { diff --git a/src/ExternalVariables.h b/src/ExternalVariables.h index 54a5129..5f8c1e6 100644 --- a/src/ExternalVariables.h +++ b/src/ExternalVariables.h @@ -16,12 +16,15 @@ extern CCTK_INT *EvolvedVariableIndex; +extern CCTK_INT *EvolvedVariableIndexSlow; extern CCTK_INT *RHSVariableIndex; +extern CCTK_INT *RHSVariableIndexSlow; extern CCTK_INT *ConstrainedVariableIndex; extern CCTK_INT *SandRVariableIndex; extern CCTK_INT MoLNumEvolvedVariables; +extern CCTK_INT MoLNumEvolvedVariablesSlow; extern CCTK_INT MoLNumConstrainedVariables; extern CCTK_INT MoLNumSandRVariables; diff --git a/src/IndexArrays.c b/src/IndexArrays.c index d3ceccb..485a4bd 100644 --- a/src/IndexArrays.c +++ b/src/IndexArrays.c @@ -110,6 +110,23 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS) } } + if (MoL_Num_Evolved_Vars_Slow) + { + EvolvedVariableIndexSlow = (CCTK_INT *)malloc(MoL_Num_Evolved_Vars_Slow * + sizeof(CCTK_INT)); + if (!EvolvedVariableIndexSlow) + { + CCTK_WARN(0,"Failed to allocate the slow evolved variable index array"); + } + + RHSVariableIndexSlow = (CCTK_INT *)malloc(MoL_Num_Evolved_Vars_Slow * + sizeof(CCTK_INT)); + if (!RHSVariableIndexSlow) + { + CCTK_WARN(0,"Failed to allocate the slow RHS variable index array"); + } + } + if (MoL_Num_Constrained_Vars) { ConstrainedVariableIndex = (CCTK_INT *)malloc(MoL_Num_Constrained_Vars * @@ -335,6 +352,18 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS) "Averaging iterative Crank Nicholson with %i iterations", MoL_Intermediate_Steps); } + else if (CCTK_EQUALS(ODE_Method,"RK2-MR-2:1")) + { + sprintf(infoline, "Multi-rate 2:1 Runge-Kutta 2"); + } + else if (CCTK_EQUALS(ODE_Method,"RK4-MR-2:1")) + { + sprintf(infoline, "Multi-rate 2:1 Runge-Kutta 4"); + } + else if (CCTK_EQUALS(ODE_Method,"RK4-RK2")) + { + sprintf(infoline, "Multi-rate 2:1 Runge-Kutta 4 and Runge-Kutta 2"); + } else { CCTK_WARN(0, "ODE_Method not recognized!"); @@ -345,6 +374,12 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS) free(infoline); infoline = NULL; + // These scalars must be 1 oustide of the MoL loop. + // They will only be zero for certain substeps when multirate methods are used. + // Otherwise, they are guaranteed to always be ONE. + *MoL_SlowPostStep = 1; + *MoL_SlowStep = 1; + return; } @@ -381,6 +416,18 @@ void MoL_FreeIndexArrays(CCTK_ARGUMENTS) RHSVariableIndex = NULL; } + if (EvolvedVariableIndexSlow) + { + free(EvolvedVariableIndexSlow); + EvolvedVariableIndexSlow = NULL; + } + + if (RHSVariableIndexSlow) + { + free(RHSVariableIndexSlow); + RHSVariableIndexSlow = NULL; + } + if (ConstrainedVariableIndex) { free(ConstrainedVariableIndex); diff --git a/src/MoL.h b/src/MoL.h index 2135870..707c7e2 100644 --- a/src/MoL.h +++ b/src/MoL.h @@ -20,6 +20,7 @@ extern "C" { #ifdef CCODE CCTK_INT MoL_RegisterEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex); +CCTK_INT MoL_RegisterEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndexSlow); CCTK_INT MoL_RegisterConstrained(CCTK_INT ConstrainedIndex); CCTK_INT MoL_RegisterSaveAndRestore(CCTK_INT SandRIndex); CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex); @@ -28,6 +29,8 @@ CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex); CCTK_INT MoL_ChangeToNone(CCTK_INT RemoveIndex); CCTK_INT MoL_RegisterEvolvedGroup(CCTK_INT EvolvedGroupIndex, CCTK_INT RHSGroupIndex); +CCTK_INT MoL_RegisterEvolvedGroupSlow(CCTK_INT EvolvedGroupIndexSlow, + CCTK_INT RHSGroupIndexSlow); CCTK_INT MoL_RegisterConstrainedGroup(CCTK_INT ConstrainedGroupIndex); CCTK_INT MoL_RegisterSaveAndRestoreGroup(CCTK_INT SandRGroupIndex); diff --git a/src/MoLFunctions.h b/src/MoLFunctions.h index 9f99e73..6d92bf0 100644 --- a/src/MoLFunctions.h +++ b/src/MoLFunctions.h @@ -13,6 +13,7 @@ #define MOL_FUNCTIONS_H CCTK_INT MoL_RegisterEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex); +CCTK_INT MoL_RegisterEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndexSlow); CCTK_INT MoL_RegisterConstrained(CCTK_INT ConstrainedIndex); CCTK_INT MoL_RegisterSaveAndRestore(CCTK_INT SandRIndex); CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex); @@ -21,6 +22,8 @@ CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex); CCTK_INT MoL_ChangeToNone(CCTK_INT RemoveIndex); CCTK_INT MoL_RegisterEvolvedGroup(CCTK_INT EvolvedGroupIndex, CCTK_INT RHSGroupIndex); +CCTK_INT MoL_RegisterEvolvedGroupSlow(CCTK_INT EvolvedGroupIndexSlow, + CCTK_INT RHSGroupIndexSlow); CCTK_INT MoL_RegisterConstrainedGroup(CCTK_INT ConstrainedGroupIndex); CCTK_INT MoL_RegisterSaveAndRestoreGroup(CCTK_INT SandRGroupIndex); diff --git a/src/ParamCheck.c b/src/ParamCheck.c index 0257eb2..5ae51d6 100644 --- a/src/ParamCheck.c +++ b/src/ParamCheck.c @@ -159,6 +159,24 @@ void MoL_ParamCheck(CCTK_ARGUMENTS) "number of intermediate steps must be 13 " "and the number of scratch levels at least 13."); } + + if ( (CCTK_Equals(ODE_Method, "RK2-MR-2:1"))&&( !((MoL_Intermediate_Steps == 5) || (MoL_Num_Scratch_Levels > 4))) ) + { + CCTK_PARAMWARN("When using the multirate 2-1 RK2 evolver the " + "number of intermediate steps must be 5 and the number of scratch levels at least 5"); + } + + if ( (CCTK_Equals(ODE_Method, "RK4-MR-2:1"))&&( !((MoL_Intermediate_Steps == 10) || (MoL_Num_Scratch_Levels > 9))) ) + { + CCTK_PARAMWARN("When using the multirate 2-1 RK4 evolver the " + "number of intermediate steps must be 10 and the number of scratch levels at least 10"); + } + + if ( (CCTK_Equals(ODE_Method, "RK4-RK2"))&&( !((MoL_Intermediate_Steps == 4) || (MoL_Num_Scratch_Levels > 3))) ) + { + CCTK_PARAMWARN("When using the multirate RK4-RK2 evolver the " + "number of intermediate steps must be 4 and the number of scratch levels at least 3"); + } if (adaptive_stepsize) { diff --git a/src/RK2-MR-2_1.c b/src/RK2-MR-2_1.c new file mode 100644 index 0000000..8eaabd2 --- /dev/null +++ b/src/RK2-MR-2_1.c @@ -0,0 +1,358 @@ + /*@@ + @file RK2-MR-2_1.c + @date 2012-03-25 + @author Christian Reisswig + @desc + A routine to perform RK2 2:1 MR evolution. Mostly copied from + genericRK.c + @enddesc + @version $Header$ + @@*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include +#include "ExternalVariables.h" + +//#define MOLDEBUG + +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusBase_MoL_RK2_MR_2_1_c); + +/******************************************************************** + ********************* Local Data Types *********************** + ********************************************************************/ + +/******************************************************************** + ********************* Local Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ***************** Scheduled Routine Prototypes ********************* + ********************************************************************/ + +void MoL_RK2_MR_2_1_Add(CCTK_ARGUMENTS); + +/******************************************************************** + ********************* Other Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ********************* Local Data ***************************** + ********************************************************************/ + +/******************************************************************** + ********************* External Routines ********************** + ********************************************************************/ + + /*@@ + @routine MoL_RK2_MR_2_1_Add + @date + @author + @desc + Performs a single step of a RK2_MR_2_1 type time + integration. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +void MoL_RK2_MR_2_1_Add(CCTK_ARGUMENTS) +{ + + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + CCTK_INT 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; + static CCTK_INT scratchspace_firstindex_slow = -99; + CCTK_INT index, var, scratchstep; + CCTK_INT totalsize; + CCTK_REAL alpha[4], beta[5]; + CCTK_REAL alpha_slow[4], beta_slow[5]; + CCTK_REAL * restrict UpdateVar; + CCTK_REAL * restrict OldVar; + CCTK_REAL const * restrict RHSVar; + CCTK_REAL * restrict ScratchVar; + + totalsize = 1; + for (arraydim = 0; arraydim < cctk_dim; arraydim++) + { + totalsize *= cctk_lsh[arraydim]; + } + + if (scratchspace_firstindex == -99) + { + scratchspace_firstindex = CCTK_FirstVarIndex("MOL::SCRATCHSPACE"); + } + + if (scratchspace_firstindex_slow == -99) + { + scratchspace_firstindex_slow = CCTK_FirstVarIndex("MOL::SCRATCHSPACESLOW"); + } + + switch (MoL_Intermediate_Steps - (*MoL_Intermediate_Step)) + { + case 0: + alpha[0] = 1.0/2.0; + alpha[1] = 0; + alpha[2] = 0; + alpha[3] = 0; + alpha_slow[0] = 1.0/2.0; + alpha_slow[1] = 0; + alpha_slow[2] = 0; + alpha_slow[3] = 0; + break; + case 1: + alpha[0] = 1.0/4.0; + alpha[1] = 1.0/4.0; + alpha[2] = 0; + alpha[3] = 0; + alpha_slow[0] = 1.0/2.0; + alpha_slow[1] = 0; + alpha_slow[2] = 0; + alpha_slow[3] = 0; + break; + case 2: + alpha[0] = 1.0/4.0; + alpha[1] = 1.0/4.0; + alpha[2] = 1.0/2.0; + alpha[3] = 0.0; + alpha_slow[0] = 1.0; + alpha_slow[1] = 0; + alpha_slow[2] = 0; + alpha_slow[3] = 0; + break; + case 3: + alpha[0] = 1.0/4.0; + alpha[1] = 1.0/4.0; + alpha[2] = 1.0/4.0; + alpha[3] = 1.0/4.0; + alpha_slow[0] = 1.0; + alpha_slow[1] = 0; + alpha_slow[2] = 0; + alpha_slow[3] = 0; + } + + beta[0] = 1.0/4.0; + beta[1] = 1.0/4.0; + beta[2] = 1.0/4.0; + beta[3] = 1.0/4.0; + beta[4] = 0.0; + + beta_slow[0] = 1.0/2.0; + beta_slow[1] = 0.0; + beta_slow[2] = 0.0; + beta_slow[3] = 0.0; + beta_slow[4] = 1.0/2.0; + + + /* FIXME */ + + + /* Real GFs, the "fast" part */ + + 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 multirate RK. Variable %d (%s). RHS %d (%s). beta %g.\n", + EvolvedVariableIndex[var], + CCTK_VarName(EvolvedVariableIndex[var]), + RHSVariableIndex[var], + CCTK_VarName(RHSVariableIndex[var]), + beta[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)]); +#endif + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex + + var + + MoLNumEvolvedVariables * (MoL_Intermediate_Steps - (*MoL_Intermediate_Step))); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + ScratchVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac * 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 + } + + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] = OldVar[index]; + } + + if ((*MoL_Intermediate_Step)>1) + { + //printf("Step %d \n", MoL_Intermediate_Steps - (*MoL_Intermediate_Step)); + for (scratchstep = 0; scratchstep <= MoL_Intermediate_Steps - (*MoL_Intermediate_Step); scratchstep++) + { + + //printf("Scratch Step %d, alpha %g \n", scratchstep, alpha[scratchstep]); + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex + + var + + MoLNumEvolvedVariables * scratchstep); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += alpha[scratchstep] * ScratchVar[index]; + } + } + } + else + { + //printf("Final Step!\n"); + + for (scratchstep = 0; scratchstep < MoL_Intermediate_Steps; scratchstep++) + { + + //printf("Scratch Step %d, beta %g \n", scratchstep, beta[scratchstep]); + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex + + var + + MoLNumEvolvedVariables * scratchstep); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += beta[scratchstep] * ScratchVar[index]; + } + } + } + + } + + + for (var = 0; var < MoLNumEvolvedVariablesSlow; var++) + { + + UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, + EvolvedVariableIndexSlow[var]); + OldVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 1, + EvolvedVariableIndexSlow[var]); + RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, + RHSVariableIndexSlow[var]); +/* #define MOLDEBUG 1 */ +#ifdef MOLDEBUG + printf("In multirate RK. SLOW Variable %d (%s). RHS %d (%s). beta %g.\n", + EvolvedVariableIndexSlow[var], + CCTK_VarName(EvolvedVariableIndexSlow[var]), + RHSVariableIndexSlow[var], + CCTK_VarName(RHSVariableIndexSlow[var]), + beta[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)]); +#endif + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex_slow + + var + + MoLNumEvolvedVariablesSlow * (MoL_Intermediate_Steps - (*MoL_Intermediate_Step))); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + ScratchVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac * RHSVar[index]; + +#ifdef MOLDEBUG + if (CCTK_EQUALS(verbose,"extreme")) + { + printf("SLOW 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 + } + + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] = OldVar[index]; + } + + if ((*MoL_Intermediate_Step)>1) + { + //printf("Step %d \n", MoL_Intermediate_Steps - (*MoL_Intermediate_Step)); + for (scratchstep = 0; scratchstep <= /*MoL_Intermediate_Steps - (*MoL_Intermediate_Step)*/ 0; scratchstep++) + { + + //printf("Scratch Step %d, alpha %g \n", scratchstep, alpha_slow[scratchstep]); + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex_slow + + var + + MoLNumEvolvedVariablesSlow * scratchstep); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += alpha_slow[scratchstep] * ScratchVar[index]; + } + } + } + else + { + //printf("Final Step!\n"); + + for (scratchstep = 0; scratchstep < MoL_Intermediate_Steps; scratchstep+=4) + { + + //printf("Scratch Step %d, beta %g \n", scratchstep, beta_slow[scratchstep]); + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex_slow + + var + + MoLNumEvolvedVariablesSlow * scratchstep); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += beta_slow[scratchstep] * ScratchVar[index]; + } + } + } + + } + + return; +} diff --git a/src/RK4-MR-2_1.c b/src/RK4-MR-2_1.c new file mode 100644 index 0000000..815a41e --- /dev/null +++ b/src/RK4-MR-2_1.c @@ -0,0 +1,502 @@ + /*@@ + @file RK4-MR-2_1.c + @date 2012-03-25 + @author Christian Reisswig + @desc + A routine to perform 3rd order 2:1 multirate RK evolution. Mostly copied + from genericRK.c + @enddesc + @version $Header$ + @@*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include +#include "ExternalVariables.h" + +//#define MOLDEBUG + +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusBase_MoL_RK4_MR_2_1_c); + +/******************************************************************** + ********************* Local Data Types *********************** + ********************************************************************/ + +/******************************************************************** + ********************* Local Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ***************** Scheduled Routine Prototypes ********************* + ********************************************************************/ + +void MoL_RK4_MR_2_1_Add(CCTK_ARGUMENTS); + +/******************************************************************** + ********************* Other Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ********************* Local Data ***************************** + ********************************************************************/ + +/******************************************************************** + ********************* External Routines ********************** + ********************************************************************/ + + /*@@ + @routine MoL_RK4_MR_2_1_Add + @date + @author + @desc + Performs a single step of a RK43 2:1 type time + integration. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +void MoL_RK4_MR_2_1_Add(CCTK_ARGUMENTS) +{ + + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + CCTK_INT 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; + static CCTK_INT scratchspace_firstindex_slow = -99; + CCTK_INT index, var, scratchstep; + CCTK_INT totalsize; + CCTK_REAL alpha[9], beta[10]; + CCTK_REAL alpha_slow[9], beta_slow[10]; + CCTK_REAL * restrict UpdateVar; + CCTK_REAL * restrict OldVar; + CCTK_REAL const * restrict RHSVar; + CCTK_REAL * restrict ScratchVar; + + totalsize = 1; + for (arraydim = 0; arraydim < cctk_dim; arraydim++) + { + totalsize *= cctk_lsh[arraydim]; + } + + if (scratchspace_firstindex == -99) + { + scratchspace_firstindex = CCTK_FirstVarIndex("MOL::SCRATCHSPACE"); + } + + if (scratchspace_firstindex_slow == -99) + { + scratchspace_firstindex_slow = CCTK_FirstVarIndex("MOL::SCRATCHSPACESLOW"); + } + + switch (MoL_Intermediate_Steps - (*MoL_Intermediate_Step)) + { + case 0: + alpha[0] = 1.0/4.0; + alpha[1] = 0; + alpha[2] = 0; + alpha[3] = 0; + alpha[4] = 0; + alpha[5] = 0; + alpha[6] = 0; + alpha[7] = 0; + alpha[8] = 0; + alpha_slow[0] = 1.0/4.0; + alpha_slow[1] = 0; + alpha_slow[2] = 0; + alpha_slow[3] = 0; + alpha_slow[4] = 0.0; + alpha_slow[5] = 0; + alpha_slow[6] = 0; + alpha_slow[7] = 0; + alpha_slow[8] = 0; + break; + case 1: + alpha[0] = -1.0/12.0; + alpha[1] = 1.0/3.0; + alpha[2] = 0; + alpha[3] = 0; + alpha[4] = 0.0; + alpha[5] = 0.0; + alpha[6] = 0; + alpha[7] = 0; + alpha[8] = 0; + alpha_slow[0] = 1.0/4.0; + alpha_slow[1] = 0; + alpha_slow[2] = 0; + alpha_slow[3] = 0; + alpha_slow[4] = 0.0; + alpha_slow[5] = 0; + alpha_slow[6] = 0; + alpha_slow[7] = 0; + alpha_slow[8] = 0; + break; + case 2: + alpha[0] = 1.0/6.0; + alpha[1] = -1.0/6.0; + alpha[2] = 1.0/2.0; + alpha[3] = 0.0; + alpha[4] = 0.0; + alpha[5] = 0.0; + alpha[6] = 0.0; + alpha[7] = 0.0; + alpha[8] = 0.0; + alpha_slow[0] = 1.0/2.0; + alpha_slow[1] = 0; + alpha_slow[2] = 0; + alpha_slow[3] = 0; + alpha_slow[4] = 0.0; + alpha_slow[5] = 0; + alpha_slow[6] = 0; + alpha_slow[7] = 0; + alpha_slow[8] = 0; + break; + case 3: + alpha[0] = 1.0/12.0; + alpha[1] = 1.0/6.0; + alpha[2] = 1.0/6.0; + alpha[3] = 1.0/12.0; + alpha[4] = 0.0; + alpha[5] = 0.0; + alpha[6] = 0.0; + alpha[7] = 0.0; + alpha[8] = 0.0; + alpha_slow[0] = 1.0/2.0; + alpha_slow[1] = 0; + alpha_slow[2] = 0; + alpha_slow[3] = 0; + alpha_slow[4] = 0.0; + alpha_slow[5] = 0; + alpha_slow[6] = 0; + alpha_slow[7] = 0; + alpha_slow[8] = 0; + case 4: + alpha[0] = 1.0/12.0; + alpha[1] = 1.0/6.0; + alpha[2] = 1.0/6.0; + alpha[3] = 1.0/6.0; + alpha[4] = 1.0/12.0; + alpha[5] = 0.0; + alpha[6] = 0.0; + alpha[7] = 0.0; + alpha[8] = 0.0; + alpha_slow[0] = -1.0/6.0; + alpha_slow[1] = 0; + alpha_slow[2] = 0; + alpha_slow[3] = 0; + alpha_slow[4] = 2.0/3.0; + alpha_slow[5] = 0; + alpha_slow[6] = 0; + alpha_slow[7] = 0; + alpha_slow[7] = 0; + case 5: + alpha[0] = 1.0/12.0; + alpha[1] = 1.0/6.0; + alpha[2] = 1.0/6.0; + alpha[3] = 1.0/12.0; + alpha[4] = 0.0; + alpha[5] = 1.0/4.0; + alpha[6] = 0.0; + alpha[7] = 0.0; + alpha[8] = 0.0; + alpha_slow[0] = 1.0/12.0; + alpha_slow[1] = 0; + alpha_slow[2] = 0; + alpha_slow[3] = 0; + alpha_slow[4] = 1.0/6.0; + alpha_slow[5] = 1.0/2.0; + alpha_slow[6] = 0; + alpha_slow[7] = 0; + alpha_slow[7] = 0; + case 6: + alpha[0] = 1.0/12.0; + alpha[1] = 1.0/6.0; + alpha[2] = 1.0/6.0; + alpha[3] = 1.0/12.0; + alpha[4] = 0.0; + alpha[5] = -1.0/12.0; + alpha[6] = 1.0/3.0; + alpha[7] = 0.0; + alpha[8] = 0.0; + alpha_slow[0] = 1.0/12.0; + alpha_slow[1] = 0; + alpha_slow[2] = 0; + alpha_slow[3] = 0; + alpha_slow[4] = 1.0/6.0; + alpha_slow[5] = 1.0/2.0; + alpha_slow[6] = 0; + alpha_slow[7] = 0; + alpha_slow[7] = 0; + case 7: + alpha[0] = 1.0/12.0; + alpha[1] = 1.0/6.0; + alpha[2] = 1.0/6.0; + alpha[3] = 1.0/12.0; + alpha[4] = 0.0; + alpha[5] = 1.0/6.0; + alpha[6] = -1.0/6.0; + alpha[7] = 1.0/2.0; + alpha[8] = 0.0; + alpha_slow[0] = 1.0/3.0; + alpha_slow[1] = 0; + alpha_slow[2] = 0; + alpha_slow[3] = 0; + alpha_slow[4] = -1.0/3.0; + alpha_slow[5] = 1.0; + alpha_slow[6] = 0; + alpha_slow[7] = 0; + alpha_slow[7] = 0; + case 8: + alpha[0] = 1.0/12.0; + alpha[1] = 1.0/6.0; + alpha[2] = 1.0/6.0; + alpha[3] = 1.0/12.0; + alpha[4] = 0.0; + alpha[5] = 1.0/12.0; + alpha[6] = 1.0/6.0; + alpha[7] = 1.0/6.0; + alpha[8] = 1.0/12.0; + alpha_slow[0] = 1.0/3.0; + alpha_slow[1] = 0; + alpha_slow[2] = 0; + alpha_slow[3] = 0; + alpha_slow[4] = -1.0/3.0; + alpha_slow[5] = 1.0; + alpha_slow[6] = 0; + alpha_slow[7] = 0; + alpha_slow[7] = 0; + } + + beta[0] = 1.0/12.0; + beta[1] = 1.0/6.0; + beta[2] = 1.0/6.0; + beta[3] = 1.0/12.0; + beta[4] = 0.0; + beta[5] = 1.0/12.0; + beta[6] = 1.0/6.0; + beta[7] = 1.0/6.0; + beta[8] = 1.0/12.0; + beta[9] = 0.0; + + beta_slow[0] = 1.0/6.0; + beta_slow[1] = 0.0; + beta_slow[2] = 0.0; + beta_slow[3] = 0.0; + beta_slow[4] = 1.0/3.0; + beta_slow[5] = 1.0/3.0; + beta_slow[6] = 0.0; + beta_slow[7] = 0.0; + beta_slow[8] = 0.0; + beta_slow[9] = 1.0/6.0; + + /* FIXME */ + + + /* Real GFs, the "fast" part */ + + 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 multirate RK. Variable %d (%s). RHS %d (%s). beta %g.\n", + EvolvedVariableIndex[var], + CCTK_VarName(EvolvedVariableIndex[var]), + RHSVariableIndex[var], + CCTK_VarName(RHSVariableIndex[var]), + beta[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)]); +#endif + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex + + var + + MoLNumEvolvedVariables * (MoL_Intermediate_Steps - (*MoL_Intermediate_Step))); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + ScratchVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac * 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 + } + + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] = OldVar[index]; + } + + if ((*MoL_Intermediate_Step)>1) + { + //printf("Step %d \n", MoL_Intermediate_Steps - (*MoL_Intermediate_Step)); + for (scratchstep = 0; scratchstep <= MoL_Intermediate_Steps - (*MoL_Intermediate_Step); scratchstep++) + { + + //printf("Scratch Step %d, alpha %g \n", scratchstep, alpha[scratchstep]); + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex + + var + + MoLNumEvolvedVariables * scratchstep); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += alpha[scratchstep] * ScratchVar[index]; + } + } + } + else + { + //printf("Final Step!\n"); + + for (scratchstep = 0; scratchstep < MoL_Intermediate_Steps; scratchstep++) + { + + //printf("Scratch Step %d, beta %g \n", scratchstep, beta[scratchstep]); + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex + + var + + MoLNumEvolvedVariables * scratchstep); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += beta[scratchstep] * ScratchVar[index]; + } + } + } + + } + + + for (var = 0; var < MoLNumEvolvedVariablesSlow; var++) + { + + UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, + EvolvedVariableIndexSlow[var]); + OldVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 1, + EvolvedVariableIndexSlow[var]); + RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, + RHSVariableIndexSlow[var]); +/* #define MOLDEBUG 1 */ +#ifdef MOLDEBUG + printf("In multirate RK. SLOW Variable %d (%s). RHS %d (%s). beta %g.\n", + EvolvedVariableIndexSlow[var], + CCTK_VarName(EvolvedVariableIndexSlow[var]), + RHSVariableIndexSlow[var], + CCTK_VarName(RHSVariableIndexSlow[var]), + beta[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)]); +#endif + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex_slow + + var + + MoLNumEvolvedVariablesSlow * (MoL_Intermediate_Steps - (*MoL_Intermediate_Step))); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + ScratchVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac * RHSVar[index]; + +#ifdef MOLDEBUG + if (CCTK_EQUALS(verbose,"extreme")) + { + printf("SLOW 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 + } + + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] = OldVar[index]; + } + + if ((*MoL_Intermediate_Step)>1) + { + //printf("Step %d \n", MoL_Intermediate_Steps - (*MoL_Intermediate_Step)); + for (scratchstep = 0; scratchstep <= MoL_Intermediate_Steps - (*MoL_Intermediate_Step); scratchstep++) + { + + //printf("Scratch Step %d, alpha %g \n", scratchstep, alpha[scratchstep]); + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex_slow + + var + + MoLNumEvolvedVariablesSlow * scratchstep); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += alpha_slow[scratchstep] * ScratchVar[index]; + } + } + } + else + { + //printf("Final Step!\n"); + + for (scratchstep = 0; scratchstep < MoL_Intermediate_Steps; scratchstep++) + { + + //printf("Scratch Step %d, beta %g \n", scratchstep, beta[scratchstep]); + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex_slow + + var + + MoLNumEvolvedVariablesSlow * scratchstep); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += beta_slow[scratchstep] * ScratchVar[index]; + } + } + } + + } + + return; +} diff --git a/src/RK4-RK2.c b/src/RK4-RK2.c new file mode 100644 index 0000000..40e84fb --- /dev/null +++ b/src/RK4-RK2.c @@ -0,0 +1,342 @@ + /*@@ + @file RK4-RK2.c + @date 2012-03-25 + @author Christian Reisswig + @desc + A routine to perform homegrown RK4RK2 evolution. Mostly copied from + genericRK.c + @enddesc + @version $Header$ + @@*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include +#include "ExternalVariables.h" + +//#define MOLDEBUG + +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusBase_MoL_RK4_RK2_c); + +/******************************************************************** + ********************* Local Data Types *********************** + ********************************************************************/ + +/******************************************************************** + ********************* Local Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ***************** Scheduled Routine Prototypes ********************* + ********************************************************************/ + +void MoL_RK4_RK2_Add(CCTK_ARGUMENTS); + +/******************************************************************** + ********************* Other Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ********************* Local Data ***************************** + ********************************************************************/ + +/******************************************************************** + ********************* External Routines ********************** + ********************************************************************/ + + /*@@ + @routine MoL_RK4_RK2_Add + @date + @author + @desc + Performs a single step of a RK4_RK2 type time + integration. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +void MoL_RK4_RK2_Add(CCTK_ARGUMENTS) +{ + + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + CCTK_INT 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; + static CCTK_INT scratchspace_firstindex_slow = -99; + CCTK_INT index, var, scratchstep; + CCTK_INT totalsize; + CCTK_REAL alpha[3], beta[4]; + CCTK_REAL alpha_slow[3], beta_slow[4]; + CCTK_REAL * restrict UpdateVar; + CCTK_REAL * restrict OldVar; + CCTK_REAL const * restrict RHSVar; + CCTK_REAL * restrict ScratchVar; + + totalsize = 1; + for (arraydim = 0; arraydim < cctk_dim; arraydim++) + { + totalsize *= cctk_lsh[arraydim]; + } + + if (scratchspace_firstindex == -99) + { + scratchspace_firstindex = CCTK_FirstVarIndex("MOL::SCRATCHSPACE"); + } + + if (scratchspace_firstindex_slow == -99) + { + scratchspace_firstindex_slow = CCTK_FirstVarIndex("MOL::SCRATCHSPACESLOW"); + } + + switch (MoL_Intermediate_Steps - (*MoL_Intermediate_Step)) + { + case 0: + alpha[0] = 1.0/2.0; + alpha[1] = 0; + alpha[2] = 0; + alpha_slow[0] = 0; + alpha_slow[1] = 0; + alpha_slow[2] = 0; + break; + case 1: + alpha[0] = 0.0; + alpha[1] = 1.0/2.0; + alpha[2] = 0; + alpha_slow[0] = 0; + alpha_slow[1] = 0; + alpha_slow[2] = 0; + break; + case 2: + alpha[0] = 0.0; + alpha[1] = 0.0; + alpha[2] = 1.0; + alpha_slow[0] = 1.0; + alpha_slow[1] = 0; + alpha_slow[2] = 0; + break; + } + + beta[0] = 1.0/3.0; + beta[1] = 1.0/6.0; + beta[2] = 1.0/6.0; + beta[3] = 1.0/3.0; + + beta_slow[0] = 1.0/2.0; + beta_slow[1] = 0.0; + beta_slow[2] = 0.0; + beta_slow[3] = 1.0/2.0; + + + /* FIXME */ + + + /* Real GFs, the "fast" part */ + + 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 multirate RK. Variable %d (%s). RHS %d (%s). beta %g.\n", + EvolvedVariableIndex[var], + CCTK_VarName(EvolvedVariableIndex[var]), + RHSVariableIndex[var], + CCTK_VarName(RHSVariableIndex[var]), + beta[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)]); +#endif + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex + + var + + MoLNumEvolvedVariables * (MoL_Intermediate_Steps - (*MoL_Intermediate_Step))); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + ScratchVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac * 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 + } + + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] = OldVar[index]; + } + + + if ((*MoL_Intermediate_Step)>1) + { + //printf("Step %d \n", MoL_Intermediate_Steps - (*MoL_Intermediate_Step)); + for (scratchstep = 0; scratchstep <= MoL_Intermediate_Steps - (*MoL_Intermediate_Step); scratchstep++) + { + + //printf("Scratch Step %d, alpha %g \n", scratchstep, alpha[scratchstep]); + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex + + var + + MoLNumEvolvedVariables * scratchstep); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += alpha[scratchstep] * ScratchVar[index]; + } + } + } + else + { + //printf("Final Step!\n"); + + for (scratchstep = 0; scratchstep < MoL_Intermediate_Steps; scratchstep++) + { + + //printf("Scratch Step %d, beta %g \n", scratchstep, beta[scratchstep]); + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex + + var + + MoLNumEvolvedVariables * scratchstep); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += beta[scratchstep] * ScratchVar[index]; + } + } + } + + } + + + for (var = 0; var < MoLNumEvolvedVariablesSlow; var++) + { + + UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, + EvolvedVariableIndexSlow[var]); + OldVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 1, + EvolvedVariableIndexSlow[var]); + RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, + RHSVariableIndexSlow[var]); +/* #define MOLDEBUG 1 */ +#ifdef MOLDEBUG + printf("In multirate RK. SLOW Variable %d (%s). RHS %d (%s). beta %g.\n", + EvolvedVariableIndexSlow[var], + CCTK_VarName(EvolvedVariableIndexSlow[var]), + RHSVariableIndexSlow[var], + CCTK_VarName(RHSVariableIndexSlow[var]), + beta[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)]); +#endif + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex_slow + + var + + MoLNumEvolvedVariablesSlow * (MoL_Intermediate_Steps - (*MoL_Intermediate_Step))); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + ScratchVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac * RHSVar[index]; + +#ifdef MOLDEBUG + if (CCTK_EQUALS(verbose,"extreme")) + { + printf("SLOW 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 + } + + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] = OldVar[index]; + } + + if ((*MoL_Intermediate_Step)>1) + { + //printf("Step %d \n", MoL_Intermediate_Steps - (*MoL_Intermediate_Step)); + for (scratchstep = 0; scratchstep <= /*MoL_Intermediate_Steps - (*MoL_Intermediate_Step)*/ 0; scratchstep++) + { + + //printf("Scratch Step %d, alpha %g \n", scratchstep, alpha_slow[scratchstep]); + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex_slow + + var + + MoLNumEvolvedVariablesSlow * scratchstep); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += alpha_slow[scratchstep] * ScratchVar[index]; + } + } + } + else + { + //printf("Final Step!\n"); + + for (scratchstep = 0; scratchstep < MoL_Intermediate_Steps; scratchstep++) + { + + //printf("Scratch Step %d, beta %g \n", scratchstep, beta_slow[scratchstep]); + + ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, + scratchspace_firstindex_slow + + var + + MoLNumEvolvedVariablesSlow * scratchstep); + +#pragma omp parallel for + for (index = 0; index < totalsize; index++) + { + UpdateVar[index] += beta_slow[scratchstep] * ScratchVar[index]; + } + } + } + + } + + return; +} diff --git a/src/Registration.c b/src/Registration.c index 06cf953..5f58d36 100644 --- a/src/Registration.c +++ b/src/Registration.c @@ -40,6 +40,8 @@ void MoL_SetScheduleStatus(CCTK_ARGUMENTS); CCTK_INT MoL_RegisterEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex); +CCTK_INT MoL_RegisterEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndexSlow); + CCTK_INT MoL_RegisterConstrained(CCTK_INT ConstrainedIndex); CCTK_INT MoL_RegisterSaveAndRestore(CCTK_INT SandRIndex); @@ -52,9 +54,10 @@ CCTK_INT MoL_RegisterConstrainedGroup(CCTK_INT ConstrainedGroupIndex); CCTK_INT MoL_RegisterSaveAndRestoreGroup(CCTK_INT SandRGroupIndex); - CCTK_INT MoL_RegisterEvolvedReal(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex); +CCTK_INT MoL_RegisterEvolvedRealSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndexSlow); + CCTK_INT MoL_RegisterConstrainedReal(CCTK_INT ConstrainedIndex); CCTK_INT MoL_RegisterSaveAndRestoreReal(CCTK_INT SandRIndex); @@ -62,6 +65,10 @@ CCTK_INT MoL_RegisterSaveAndRestoreReal(CCTK_INT SandRIndex); CCTK_INT MoL_RegisterEvolvedRealGroup(CCTK_INT EvolvedGroupIndex, CCTK_INT RHSGroupIndex); +CCTK_INT MoL_RegisterEvolvedRealGroupSlow(CCTK_INT EvolvedGroupIndex, + CCTK_INT RHSGroupIndexSlow); + + CCTK_INT MoL_RegisterConstrainedRealGroup(CCTK_INT ConstrainedGroupIndex); CCTK_INT MoL_RegisterSaveAndRestoreRealGroup(CCTK_INT SandRGroupIndex); @@ -269,6 +276,109 @@ CCTK_INT MoL_RegisterEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex) } + + /*@@ + @routine MoL_RegisterEvolvedSlow + @date Thu May 30 11:36:59 2002 + @author Ian Hawke + @desc + Given the index of the GF to be evolved and the RHS GF, it stores + the indexes for later use together with various error checking. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +CCTK_INT MoL_RegisterEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndexSlow) +{ + + CCTK_INT retval, GroupType; + + retval = 0; + + if (!ScheduleStatus) + { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "MoL registration routine called too early!\n" + "Trying to register variable '%s'," + "Please ensure that all calls to MoL registration routines " + "occur within the \"MoL_Register\" timebin.", + CCTK_VarName(EvolvedIndex)); + retval++; + } + + GroupType = CCTK_GroupTypeFromVarI(EvolvedIndex); + if (GroupType < 0) + { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Evolved index %i is not a real variable index.", + EvolvedIndex); + retval++; + } + + if (!retval) + { + switch (CCTK_VarTypeI(EvolvedIndex)) + { + case CCTK_VARIABLE_REAL: + { + switch (GroupType) + { + case CCTK_GF: + { + retval += + MoL_RegisterEvolvedRealSlow(EvolvedIndex, + RHSIndexSlow); + break; + } + case CCTK_ARRAY: + case CCTK_SCALAR: + { + CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING, + "The variable '%s' is not supported for multirate methods", + CCTK_VarName(EvolvedIndex)); + retval++; + break; + } + default: + { + CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING, + "The variable '%s' is not supported for multirate methods", + CCTK_VarName(EvolvedIndex)); + retval++; + break; + } + } + break; + } + case CCTK_VARIABLE_COMPLEX: + { + CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING, + "The COMPLEX variable '%s' is not supported for multirate methods", + CCTK_VarName(EvolvedIndex)); + break; + } + default: + { + CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING, + "The variable '%s' is neither REAL nor COMPLEX.", + CCTK_VarName(EvolvedIndex)); + retval++; + break; + } + } + } + + return retval; + +} + + + /*@@ @routine MoL_RegisterConstrained @date Thu May 30 12:35:58 2002 @@ -606,6 +716,91 @@ CCTK_INT MoL_RegisterEvolvedGroup(CCTK_INT EvolvedGroupIndex, } +CCTK_INT MoL_RegisterEvolvedGroupSlow(CCTK_INT EvolvedGroupIndex, + CCTK_INT RHSGroupIndex2) +{ + + CCTK_INT retval, GroupFirstVar; + + retval = 0; + + if (!ScheduleStatus) + { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "MoL registration routine called too early!\n" + "Trying to register group '%s'," + "Please ensure that all calls to MoL registration routines " + "occur within the \"MoL_Register\" timebin.", + CCTK_GroupName(EvolvedGroupIndex)); + retval++; + } + + GroupFirstVar = CCTK_FirstVarIndexI(EvolvedGroupIndex); + if (GroupFirstVar < 0) + { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Evolved group index %i is not a real group index.", + EvolvedGroupIndex); + retval++; + } + + switch (CCTK_VarTypeI(CCTK_FirstVarIndexI(EvolvedGroupIndex))) + { + case CCTK_VARIABLE_REAL: + { + switch (CCTK_GroupTypeI(EvolvedGroupIndex)) + { + case CCTK_GF: + { + retval += + MoL_RegisterEvolvedRealGroupSlow(EvolvedGroupIndex, + RHSGroupIndex2); + break; + } + case CCTK_ARRAY: + case CCTK_SCALAR: + { + CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING, + "The group '%s' is not supported by multirate RK", + CCTK_GroupName(EvolvedGroupIndex)); + retval++; + break; + } + default: + { + CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING, + "The group '%s' is not supported by multirate RK", + CCTK_GroupName(EvolvedGroupIndex)); + retval++; + break; + } + } + break; + } + case CCTK_VARIABLE_COMPLEX: + { + CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING, + "The COMPLEX group '%s' is not supported by multirate RK", + CCTK_GroupName(EvolvedGroupIndex)); + retval++; + break; + } + default: + { + CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING, + "The group '%s' is neither REAL nor COMPLEX.", + CCTK_GroupName(EvolvedGroupIndex)); + retval++; + break; + } + } + + return retval; + +} + + + CCTK_INT MoL_RegisterConstrainedGroup(CCTK_INT ConstrainedGroupIndex) { @@ -992,6 +1187,198 @@ CCTK_INT MoL_RegisterEvolvedReal(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex) } + + /*@@ + @routine MoL_RegisterEvolvedSlow + @date Thu May 30 11:36:59 2002 + @author Ian Hawke + @desc + Given the index of the GF to be evolved and the RHS GF, it stores + the indexes for later use together with various error checking. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +CCTK_INT MoL_RegisterEvolvedRealSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndexSlow) +{ + + DECLARE_CCTK_PARAMETERS; + + CCTK_INT /* ierr, */ index, varused, numtimelevs1, numtimelevs2; + +#ifdef MOLDEBUG + printf("Arrived in MoLRegisterEvolvedSlow \n"); + printf("The indexes are %d and %d.\n",EvolvedIndex, RHSIndexSlow); + printf("These correspond to variables %s and %s.\n", + CCTK_VarName(EvolvedIndex),CCTK_VarName(RHSIndexSlow)); + printf("The pointer to EvolvedVariableIndex: %p\n", + EvolvedVariableIndex); +#endif + + if (!(CCTK_GroupTypeFromVarI(EvolvedIndex)==CCTK_GF)) + { + CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING, + "The variable %s is not a GF and so should " + "not be registered with MoLRegisterEvolved.", + CCTK_VarName(EvolvedIndex)); + } + + if (!(CCTK_GroupTypeFromVarI(RHSIndexSlow)==CCTK_GF)) + { + CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING, + "The rhs variable %s is not a GF and so should " + "not be registered with MoLRegisterEvolved.", + CCTK_VarName(RHSIndexSlow)); + } + + if (!(CCTK_VarTypeI(EvolvedIndex)==CCTK_VARIABLE_REAL)) + { + CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING, + "The variable %s is not of type CCTK_REAL and so " + "should not be registered with MoLRegisterEvolved.", + CCTK_VarName(EvolvedIndex)); + } + + if (!(CCTK_VarTypeI(RHSIndexSlow)==CCTK_VARIABLE_REAL)) + { + CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING, + "The rhs variable %s is not of type CCTK_REAL and " + "so should not be registered with MoLRegisterEvolved.", + CCTK_VarName(RHSIndexSlow)); + } + + numtimelevs1 = CCTK_MaxTimeLevelsVI(EvolvedIndex); + numtimelevs2 = CCTK_MaxTimeLevelsVI(RHSIndexSlow); + + if ( (numtimelevs1 < 0) || (numtimelevs2 < 0) ) + { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "Warning for variable index %i", EvolvedIndex); + CCTK_WARN(0, "The index passed does not correspond to a GF."); + } + + if (numtimelevs1 < 2) + { + CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING, + "Warning for variable index %i name %s", + EvolvedIndex, CCTK_VarName(EvolvedIndex)); + CCTK_WARN(0, "The GF passed only has one timelevel. " + "It must have at least two."); + } + + varused = 0; + + for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!varused); index++) + { + varused = (EvolvedIndex == EvolvedVariableIndexSlow[index]); +#ifdef MOLDEBUG + printf("Registering %d. Checking index %d which is %d\n",EvolvedIndex, + index,EvolvedVariableIndexSlow[index]); +#endif + } + + if (varused) + { + + CCTK_VWarn(2,__LINE__,__FILE__,CCTK_THORNSTRING, + "The GF %s has already been registered " + "as an evolved variable with RHS GF %s. " + "The attempt to register with RHS GF %s will be ignored", + CCTK_VarName(EvolvedIndex), + CCTK_VarName(RHSVariableIndexSlow[index-1]), + CCTK_VarName(RHSIndexSlow)); + + } + else + { + + if (MoLNumEvolvedVariablesSlow+1 > MoL_Num_Evolved_Vars_Slow) + { + CCTK_WARN(0,"You have tried to register more evolved " + "variables than the accumulator parameter " + "MoL_Num_Evolved_Variables allows. Check that " + "you are accumulating onto this parameter correctly"); + } + + EvolvedVariableIndexSlow[MoLNumEvolvedVariablesSlow] = EvolvedIndex; + RHSVariableIndexSlow[MoLNumEvolvedVariablesSlow] = RHSIndexSlow; + + MoLNumEvolvedVariablesSlow++; + +#ifdef MOLDEBUG + printf("The max number is now %d. Just added %d (%s).\n", + MoLNumEvolvedVariablesSlow, EvolvedIndex, + CCTK_VarName(EvolvedIndex)); +#endif + + } + + varused = -1; + + for (index = 0; (index < MoLNumConstrainedVariables)&&(!(varused+1)); + index++) + { + if (EvolvedIndex == ConstrainedVariableIndex[index]) + { + varused = index; + } + + } + + if ((varused+1)) + { + for (index = varused; index < MoLNumConstrainedVariables-1; index++) + { + ConstrainedVariableIndex[index] = ConstrainedVariableIndex[index+1]; + } + MoLNumConstrainedVariables--; + } + + varused = -1; + + for (index = 0; (index < MoLNumSandRVariables)&&(!(varused+1)); index++) + { + if (EvolvedIndex == SandRVariableIndex[index]) + { + varused = index; + } + +#ifdef MOLDEBUG + printf("Checking SandR var %d. Index %d (evolvedindex %d).\n", + index, SandRVariableIndex[index], EvolvedIndex); +#endif + + } + + if ((varused+1)) + { + for (index = varused; index < MoLNumSandRVariables-1; index++) + { + SandRVariableIndex[index] = SandRVariableIndex[index+1]; + +#ifdef MOLDEBUG + printf("The registered evolved variable was SandR." + " Now index %d is %d (%s).\n", + index, SandRVariableIndex[index], + CCTK_VarName(SandRVariableIndex[index])); +#endif + + } + MoLNumSandRVariables--; + } + + return 0; + +} + + + + /*@@ @routine MoL_RegisterConstrained @date Thu May 30 12:35:58 2002 @@ -1256,6 +1643,51 @@ CCTK_INT MoL_RegisterEvolvedRealGroup(CCTK_INT EvolvedGroupIndex, } +CCTK_INT MoL_RegisterEvolvedRealGroupSlow(CCTK_INT EvolvedGroupIndex, + CCTK_INT RHSGroupIndexSlow) +{ + + CCTK_INT EvolvedGroupFirstVar, RHSGroupFirstVar, GroupNumVars, retval; + CCTK_INT EvolvedVar, RHSVar; + + EvolvedGroupFirstVar = CCTK_FirstVarIndexI(EvolvedGroupIndex); + if (EvolvedGroupFirstVar < 0) + { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Evolved group index %i is not a real group index.", + EvolvedGroupIndex); + } + + RHSGroupFirstVar = CCTK_FirstVarIndexI(RHSGroupIndexSlow); + if (RHSGroupFirstVar < 0) + { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "RHS group index %d is not a real group index.", + RHSGroupIndexSlow); + } + + GroupNumVars = CCTK_NumVarsInGroupI(EvolvedGroupIndex); + if (CCTK_NumVarsInGroupI(RHSGroupIndexSlow) != GroupNumVars) + { + CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, + "There are a different number of variables in " + "evolved group %d and RHS group %d.", + EvolvedGroupIndex, RHSGroupIndexSlow); + } + + retval = 0; + + for (EvolvedVar = EvolvedGroupFirstVar, RHSVar = RHSGroupFirstVar; + EvolvedVar < EvolvedGroupFirstVar + GroupNumVars; + EvolvedVar++, RHSVar++) + { + retval += MoL_RegisterEvolvedRealSlow(EvolvedVar, RHSVar); + } + + return retval; +} + + CCTK_INT MoL_RegisterConstrainedRealGroup(CCTK_INT ConstrainedGroupIndex) { @@ -2741,6 +3173,7 @@ CCTK_INT MoL_RegisterSaveAndRestoreComplexArrayGroup(CCTK_INT SandRGroupIndex) @desc Given the index of a registered evolved GF, returns the corresponding RHS GF. + Checks both regular and slow evolved variables. @enddesc @calls @calledby @@ -2769,6 +3202,17 @@ CCTK_INT MoL_QueryEvolvedRHS(CCTK_INT EvolvedIndex) } } + // If we haven't found any registered variable that matches, + // we also try to check whether this has been registered with the slow multirate sector! + + for (index = 0; index < MoLNumEvolvedVariablesSlow; index++) + { + if (EvolvedVariableIndexSlow[index] == EvolvedIndex) + { + return RHSVariableIndexSlow[index]; + } + } + return -1; } diff --git a/src/SetTime.c b/src/SetTime.c index c04cc0b..32e893c 100644 --- a/src/SetTime.c +++ b/src/SetTime.c @@ -298,6 +298,63 @@ void MoL_ResetTime(CCTK_ARGUMENTS) + ((alpha_array87[substep] - 1) * (* Original_Delta_Time) / cctkGH->cctk_timefac)); } + else if (CCTK_EQUALS(ODE_Method,"RK2-MR-2:1")) + { + 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,"RK4-MR-2:1")) + { + 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 *= 3.0/4.0; + break; + case 3: + case 4: + case 5: + dt *= 1.0/2.0; + break; + case 6: + case 7: + dt *= 1.0/4.0; + break; + default: + dt = 0; + } + cctkGH->cctk_time = (*Original_Time) - dt; + } + else if (CCTK_EQUALS(ODE_Method,"RK4-RK2")) + { + 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; + } #ifdef MOLDEBUG printf("MoL has once more reset t (%d): %f.\n", *MoL_Intermediate_Step, cctkGH->cctk_time); diff --git a/src/Startup.c b/src/Startup.c index 06ae44d..e5a2d5e 100644 --- a/src/Startup.c +++ b/src/Startup.c @@ -26,11 +26,14 @@ CCTK_FILEVERSION(CactusBase_MoL_Startup_c); ********************************************************************/ CCTK_INT *EvolvedVariableIndex = NULL; +CCTK_INT *EvolvedVariableIndexSlow = NULL; CCTK_INT *RHSVariableIndex = NULL; +CCTK_INT *RHSVariableIndexSlow = NULL; CCTK_INT *ConstrainedVariableIndex = NULL; CCTK_INT *SandRVariableIndex = NULL; CCTK_INT MoLNumEvolvedVariables = 0; +CCTK_INT MoLNumEvolvedVariablesSlow = 0; CCTK_INT MoLNumConstrainedVariables = 0; CCTK_INT MoLNumSandRVariables = 0; diff --git a/src/StepSize.c b/src/StepSize.c index cb5499e..b20c752 100644 --- a/src/StepSize.c +++ b/src/StepSize.c @@ -358,6 +358,10 @@ void MoL_FinishLoop(CCTK_ARGUMENTS) /* Keep time step size unchanged */ *MoL_Stepsize_Bad = 0; + + // Set these flags to ONE outside of MoL-loop! + *MoL_SlowPostStep = 1; + *MoL_SlowStep = 1; } /******************************************************************** diff --git a/src/make.code.defn b/src/make.code.defn index 926dbbb..80e1b56 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -21,8 +21,11 @@ SRCS = ChangeType.c \ SandR.c \ SetTime.c \ Startup.c \ - StepSize.c - + StepSize.c \ + RK2-MR-2_1.c \ + RK4-MR-2_1.c \ + RK4-RK2.c + # Subdirectories containing source files SUBDIRS = -- cgit v1.2.3