diff options
author | rhaas <rhaas@578cdeb0-5ea1-4b81-8215-5a3b8777ee0b> | 2012-08-02 16:34:52 +0000 |
---|---|---|
committer | rhaas <rhaas@578cdeb0-5ea1-4b81-8215-5a3b8777ee0b> | 2012-08-02 16:34:52 +0000 |
commit | 26067015fd74e5079d09a589fc3bab9e7fd13f22 (patch) | |
tree | a3939c89d8d4643ebf529722c1edfc8b8cbbe2a6 /src/RK4-RK2.c | |
parent | 7923fe9d90ae7b78450708403bd7d451e0227d15 (diff) |
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
Diffstat (limited to 'src/RK4-RK2.c')
-rw-r--r-- | src/RK4-RK2.c | 342 |
1 files changed, 342 insertions, 0 deletions
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 <stdio.h> +#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; +} |