aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--interface.ccl4
-rw-r--r--param.ccl1
-rw-r--r--schedule.ccl9
-rw-r--r--src/IndexArrays.c4
-rw-r--r--src/ParamCheck.c8
-rw-r--r--src/RK45.c236
-rw-r--r--src/make.code.defn1
7 files changed, 263 insertions, 0 deletions
diff --git a/interface.ccl b/interface.ccl
index ac744ed..f9f45f5 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -220,3 +220,7 @@ CCTK_REAL ArraySandRScratchSpace TYPE = ARRAY DIM = 1 SIZE = MoL_Max_Evolved_Arr
#CCTK_COMPLEX ComplexArrayScratchSpace[MoL_Num_Scratch_Levels] TYPE = ARRAY DIM = 1 SIZE = MoL_Max_Evolved_ComplexArray_Size Timelevels = 1
#CCTK_COMPLEX ComplexArraySandRScratchSpace TYPE = ARRAY DIM = 1 SIZE = MoL_Max_Evolved_ComplexArray_Size Timelevels = 1
+
+#Error vector for RK45
+
+CCTK_REAL ErrorEstimate[MoL_Num_Evolved_Vars] TYPE=GF Timelevels = 1 tags='Prolongation="None"'
diff --git a/param.ccl b/param.ccl
index e982acb..a65216a 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"
+ "RK45" :: "RK45 with error estimation"
} "ICN"
KEYWORD Generic_Type "If using the generic method, which sort"
diff --git a/schedule.ccl b/schedule.ccl
index 8f2fcf6..43dc1f9 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -513,6 +513,15 @@ else if (CCTK_Equals(ODE_Method,"RK3"))
LANG: C
} "Updates calculated with the efficient Runge-Kutta 3 method"
}
+else if (CCTK_Equals(ODE_Method,"RK45"))
+{
+ STORAGE: ErrorEstimate
+
+ schedule MoL_RK45Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep
+ {
+ LANG: C
+ } "Updates calculated with the Runge-Kutta 45 method"
+}
else if (CCTK_Equals(ODE_Method,"ICN"))
{
schedule MoL_ICNAdd AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep
diff --git a/src/IndexArrays.c b/src/IndexArrays.c
index 096a745..60cd140 100644
--- a/src/IndexArrays.c
+++ b/src/IndexArrays.c
@@ -299,6 +299,10 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS)
{
sprintf(infoline, "Runge-Kutta 3");
}
+ else if (CCTK_EQUALS(ODE_Method,"RK45"))
+ {
+ sprintf(infoline, "Runge-Kutta 45");
+ }
else if (CCTK_EQUALS(ODE_Method,"ICN"))
{
sprintf(infoline, "Iterative Crank Nicholson with %i iterations",
diff --git a/src/ParamCheck.c b/src/ParamCheck.c
index 5316c9f..bd91e5e 100644
--- a/src/ParamCheck.c
+++ b/src/ParamCheck.c
@@ -129,6 +129,14 @@ int MoL_ParamCheck(CCTK_ARGUMENTS)
}
+ if ( (CCTK_Equals(ODE_Method, "RK45")) &&
+ ( !((MoL_Intermediate_Steps == 6)||(MoL_Num_Scratch_Levels > 5)) ) )
+ {
+ CCTK_PARAMWARN("When using the RK45 evolver the "
+ "number of intermediate steps must be 6 "
+ "and the number of scratch levels at least 6.");
+ }
+
return 0;
}
diff --git a/src/RK45.c b/src/RK45.c
new file mode 100644
index 0000000..41a0462
--- /dev/null
+++ b/src/RK45.c
@@ -0,0 +1,236 @@
+ /*@@
+ @file RK45.c
+ @date Sun May 26 03:47:15 2002
+ @author Ian Hawke
+ @desc
+ RK45 following Forsythe, Malcolm and Moler
+ (Computer Methods for Mathematical Computations).
+ @enddesc
+ @version $Header$
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "ExternalVariables.h"
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusBase_MoL_RK45_c);
+
+/********************************************************************
+ ********************* Local Data Types ***********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ ***************** Scheduled Routine Prototypes *********************
+ ********************************************************************/
+
+void MoL_RK45Add(CCTK_ARGUMENTS);
+
+/********************************************************************
+ ********************* Other Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Local Data *****************************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* External Routines **********************
+ ********************************************************************/
+
+ /*@@
+ @routine MoL_RK45Add
+ @date Sun May 26 03:50:44 2002
+ @author Ian Hawke
+ @desc
+ Performs a single step of a Runge-Kutta 45 type time
+ integration, storing the error estimate.
+ @enddesc
+ @calls
+ @calledby
+ @history
+
+ @endhistory
+
+@@*/
+
+void MoL_RK45Add(CCTK_ARGUMENTS)
+{
+
+ DECLARE_CCTK_ARGUMENTS;
+ DECLARE_CCTK_PARAMETERS;
+
+ cGroupDynamicData arraydata;
+ CCTK_INT groupindex, ierr;
+ CCTK_INT arraytotalsize, arraydim;
+
+ CCTK_INT index, var, scratchstep, alphaindex, scratchindex;
+ CCTK_INT totalsize;
+
+ CCTK_REAL *UpdateVar;
+ CCTK_REAL *RHSVar;
+ CCTK_REAL *ScratchVar;
+ CCTK_REAL *ErrorVar;
+ CCTK_REAL *OldVar;
+
+ CCTK_INT arrayscratchlocation;
+
+ CCTK_REAL beta, gamma, gamma_error;
+
+ static const CCTK_REAL beta_array[5][5] = {
+ { 0.25, 0.0, 0.0, 0.0, 0.0 },
+ { 0.09375, 0.28125, 0.0, 0.0, 0.0 },
+ { 0.879380974055530268548020027310, -3.27719617660446062812926718252, 3.32089212562585343650432407829, 0.0, 0.0 },
+ { 2.03240740740740740740740740741, -8.0, 7.17348927875243664717348927875, -0.205896686159844054580896686160, 0.0 },
+ { -0.296296296296296296296296296296, 2.0, -1.38167641325536062378167641326, 0.452972709551656920077972709552, -0.275 }
+ };
+
+ static const CCTK_REAL gamma_array[6] =
+ { 0.118518518518518518518518518519,
+ 0.0,
+ 0.518986354775828460038986354776,
+ 0.506131490342016657806131490342,
+ -0.18,
+ 0.0363636363636363636363636363636
+ };
+
+ static const CCTK_REAL gammastar_array[6] =
+ { 0.115740740740740740740740740741,
+ 0.0,
+ 0.548927875243664717348927875244,
+ 0.535331384015594541910331384016,
+ -0.2,
+ 0.0
+ };
+
+ totalsize = 1;
+ for (arraydim = 0; arraydim < cctk_dim; arraydim++)
+ {
+ totalsize *= cctk_lsh[arraydim];
+ }
+
+ /* Real GFs */
+
+ /* First store (dt times) the rhs in the scratch array. */
+
+ for (var = 0; var < MoLNumEvolvedVariables; var++)
+ {
+ UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0,
+ EvolvedVariableIndex[var]);
+ RHSVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0,
+ RHSVariableIndex[var]);
+ ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
+ CCTK_FirstVarIndex("MOL::SCRATCHSPACE")
+ + var
+ + MoL_Num_Evolved_Vars *
+ (MoL_Intermediate_Steps -
+ (*MoL_Intermediate_Step)));
+ for (index = 0; index < totalsize; index++)
+ {
+ ScratchVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac *
+ RHSVar[index];
+ }
+ }
+
+
+ for (var = 0; var < MoLNumEvolvedVariables; var++)
+ {
+ OldVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 1,
+ EvolvedVariableIndex[var]);
+ UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0,
+ EvolvedVariableIndex[var]);
+ RHSVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0,
+ CCTK_FirstVarIndex("MOL::SCRATCHSPACE")
+ + var
+ + MoL_Num_Evolved_Vars *
+ (MoL_Intermediate_Steps -
+ (*MoL_Intermediate_Step)));
+ ErrorVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
+ CCTK_FirstVarIndex("MOL::ERRORESTIMATE")
+ + var);
+
+ for (index = 0; index < totalsize; index++)
+ {
+ UpdateVar[index] = OldVar[index];
+ ErrorVar[index] = 0;
+ }
+
+ if (*MoL_Intermediate_Step - 1)
+ {
+
+ for (scratchstep = 0;
+ scratchstep < MoL_Intermediate_Steps - (*MoL_Intermediate_Step) + 1;
+ scratchstep++)
+ {
+
+ ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
+ CCTK_FirstVarIndex("MOL::SCRATCHSPACE")
+ + var
+ + MoL_Num_Evolved_Vars * scratchstep);
+
+ beta = beta_array[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)][scratchstep];
+
+ if ( (beta > MoL_Tiny)||(beta < -MoL_Tiny) )
+ {
+ for (index = 0; index < totalsize; index++)
+ {
+ UpdateVar[index] += beta * ScratchVar[index];
+ }
+ }
+
+ }
+
+ }
+ else
+ {
+ for (scratchstep = 0; scratchstep < 6; scratchstep++)
+ {
+
+ ScratchVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
+ CCTK_FirstVarIndex("MOL::SCRATCHSPACE")
+ + var
+ + MoL_Num_Evolved_Vars * scratchstep);
+
+ gamma = gamma_array[scratchstep];
+ gamma_error = gamma - gammastar_array[scratchstep];
+
+ if ( (gamma > MoL_Tiny)||(gamma < -MoL_Tiny) )
+ {
+ for (index = 0; index < totalsize; index++)
+ {
+ UpdateVar[index] += gamma * ScratchVar[index];
+ ErrorVar[index] += gamma_error * ScratchVar[index];
+ }
+ }
+ }
+
+ }
+
+ }
+
+ /* Real arrays */
+
+ arrayscratchlocation = 0;
+
+ for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
+ {
+
+ CCTK_WARN(0, "Ian has been too lazy to write the RK45 routine "
+ "out for array variables. Better send him an email...");
+
+ }
+
+ return;
+}
+
+/********************************************************************
+ ********************* Local Routines *************************
+ ********************************************************************/
+
diff --git a/src/make.code.defn b/src/make.code.defn
index adfbe27..544c7bd 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -11,6 +11,7 @@ SRCS = ChangeType.c \
ParamCheck.c \
RK2.c \
RK3.c \
+ RK45.c \
Registration.c \
RKCoefficients.c \
RHSNaNCheck.c \
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-08 19:12:37 +0000