aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/IndexArrays.c6
-rw-r--r--src/ParamCheck.c8
-rw-r--r--src/RK45.c61
-rw-r--r--src/SetTime.c44
4 files changed, 103 insertions, 16 deletions
diff --git a/src/IndexArrays.c b/src/IndexArrays.c
index 505b1ab..5e59734 100644
--- a/src/IndexArrays.c
+++ b/src/IndexArrays.c
@@ -306,7 +306,11 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS)
}
else if (CCTK_EQUALS(ODE_Method,"RK45"))
{
- sprintf(infoline, "Runge-Kutta 45");
+ sprintf(infoline, "Runge-Kutta 45 (Fehlberg)");
+ }
+ else if (CCTK_EQUALS(ODE_Method,"RK45CK"))
+ {
+ sprintf(infoline, "Runge-Kutta 45 (Cash-Karp)");
}
else if (CCTK_EQUALS(ODE_Method,"RK65"))
{
diff --git a/src/ParamCheck.c b/src/ParamCheck.c
index 272117b..c2346b0 100644
--- a/src/ParamCheck.c
+++ b/src/ParamCheck.c
@@ -128,10 +128,10 @@ int MoL_ParamCheck(CCTK_ARGUMENTS)
"number of intermediate steps must be 3");
}
- if ( (CCTK_Equals(ODE_Method, "RK45")) &&
+ if ( (CCTK_Equals(ODE_Method, "RK45") || CCTK_Equals(ODE_Method, "RK45CK")) &&
( !((MoL_Intermediate_Steps == 6)||(MoL_Num_Scratch_Levels > 5)) ) )
{
- CCTK_PARAMWARN("When using the RK45 evolver the "
+ CCTK_PARAMWARN("When using the RK45 or RK45CK evolver, the "
"number of intermediate steps must be 6 "
"and the number of scratch levels at least 6.");
}
@@ -154,13 +154,13 @@ int MoL_ParamCheck(CCTK_ARGUMENTS)
if (adaptive_stepsize)
{
- if (CCTK_Equals(ODE_Method, "RK45")||CCTK_Equals(ODE_Method, "RK65")||CCTK_Equals(ODE_Method, "RK87"))
+ if (CCTK_Equals(ODE_Method, "RK45")||CCTK_Equals(ODE_Method, "RK45CK")||CCTK_Equals(ODE_Method, "RK65")||CCTK_Equals(ODE_Method, "RK87"))
{
/* everything is fine, do nothing */
}
else
{
- CCTK_PARAMWARN("Adaptive time step sizes are only possible with the RK45, RK65 and RK87 solvers");
+ CCTK_PARAMWARN("Adaptive time step sizes are only possible with the RK45, RK45CK, RK65, and RK87 solvers");
}
}
diff --git a/src/RK45.c b/src/RK45.c
index f92159a..3ca2456 100644
--- a/src/RK45.c
+++ b/src/RK45.c
@@ -67,11 +67,9 @@ void MoL_RK45Add(CCTK_ARGUMENTS)
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
- cGroupDynamicData arraydata;
- CCTK_INT groupindex, ierr;
- CCTK_INT arraytotalsize, arraydim;
+ CCTK_INT arraydim;
- CCTK_INT index, var, scratchstep, alphaindex, scratchindex;
+ CCTK_INT index, var, scratchstep;
CCTK_INT totalsize;
CCTK_REAL * restrict UpdateVar;
@@ -84,7 +82,8 @@ void MoL_RK45Add(CCTK_ARGUMENTS)
CCTK_REAL beta, gamma, gamma_error;
- static const CCTK_REAL beta_array[5][5] = {
+ /* Fehlberg coefficients */
+ static const CCTK_REAL beta_array_F[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 },
@@ -92,7 +91,7 @@ void MoL_RK45Add(CCTK_ARGUMENTS)
{ -0.296296296296296296296296296296, 2.0, -1.38167641325536062378167641326, 0.452972709551656920077972709552, -0.275 }
};
- static const CCTK_REAL gamma_array[6] =
+ static const CCTK_REAL gamma_array_F[6] =
{ 0.118518518518518518518518518519,
0.0,
0.518986354775828460038986354776,
@@ -101,7 +100,7 @@ void MoL_RK45Add(CCTK_ARGUMENTS)
0.0363636363636363636363636363636
};
- static const CCTK_REAL gammastar_array[6] =
+ static const CCTK_REAL gammastar_array_F[6] =
{ 0.115740740740740740740740740741,
0.0,
0.548927875243664717348927875244,
@@ -110,6 +109,54 @@ void MoL_RK45Add(CCTK_ARGUMENTS)
0.0
};
+ /* Cash-Karp coefficients */
+ static const CCTK_REAL beta_array_CK[5][5] = {
+ { 1.0/5.0, 0.0, 0.0, 0.0, 0.0, },
+ { 3.0/40.0, 9.0/40.0, 0.0, 0.0, 0.0, },
+ { 3.0/10.0, -9.0/10.0, 6.0/5.0, 0.0, 0.0, },
+ { -11.0/54.0, 5.0/2.0, -70.0/27.0, 35.0/27.0, 0.0, },
+ { 1631.0/55296.0, 175.0/512.0, 575.0/13824.0, 44275.0/110592.0, 253.0/4096.0, },
+ };
+
+ static const CCTK_REAL gamma_array_CK[6] = {
+ 37.0/378.0,
+ 0.0,
+ 250.0/621.0,
+ 125.0/594.0,
+ 0.0,
+ 512.0/1771.0,
+ };
+
+ static const CCTK_REAL gammastar_array_CK[6] = {
+ 2825.0/27648.0,
+ 0.0,
+ 18575.0/48384.0,
+ 13525.0/55296.0,
+ 277.0/14336.0,
+ 1.0/4.0,
+ };
+
+ const CCTK_REAL (* restrict beta_array)[5];
+ const CCTK_REAL * restrict gamma_array;
+ const CCTK_REAL * restrict gammastar_array;
+
+ if (CCTK_EQUALS(ODE_Method, "RK45"))
+ {
+ beta_array = beta_array_F;
+ gamma_array = gamma_array_F;
+ gammastar_array = gammastar_array_F;
+ }
+ else if (CCTK_EQUALS(ODE_Method, "RK45CK"))
+ {
+ beta_array = beta_array_CK;
+ gamma_array = gamma_array_CK;
+ gammastar_array = gammastar_array_CK;
+ }
+ else
+ {
+ CCTK_WARN (0, "internal error");
+ }
+
totalsize = 1;
for (arraydim = 0; arraydim < cctk_dim; arraydim++)
{
diff --git a/src/SetTime.c b/src/SetTime.c
index 1d1596d..3e9ae3f 100644
--- a/src/SetTime.c
+++ b/src/SetTime.c
@@ -47,8 +47,8 @@ int MoL_ResetDeltaTime(CCTK_ARGUMENTS);
********************* Local Data *****************************
********************************************************************/
-/* RK45 coefficients */
-static const CCTK_REAL alpha_array[6] = {
+/* RK45 Fehlberg coefficients */
+static const CCTK_REAL alpha_array_F[6] = {
0.0,
1.0/4.0,
3.0/8.0,
@@ -57,6 +57,16 @@ static const CCTK_REAL alpha_array[6] = {
1.0/2.0,
};
+/* RK45 Cash-Karp coefficients */
+static const CCTK_REAL alpha_array_CK[6] = {
+ 0.0,
+ 1.0/5.0,
+ 3.0/10.0,
+ 3.0/5.0,
+ 1.0,
+ 7.0/8.0,
+};
+
/* RK65 coefficients */
static const CCTK_REAL alpha_array65[8] = {
0.0,
@@ -237,9 +247,22 @@ int MoL_ResetTime(CCTK_ARGUMENTS)
0.5*(*Original_Delta_Time)/cctkGH->cctk_timefac;
}
}
- else if (CCTK_EQUALS(ODE_Method,"RK45"))
+ else if (CCTK_EQUALS(ODE_Method,"RK45") || CCTK_EQUALS(ODE_Method,"RK45CK"))
{
const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step);
+ const CCTK_REAL * alpha_array;
+ if (CCTK_EQUALS(ODE_Method, "RK45"))
+ {
+ alpha_array = alpha_array_F;
+ }
+ else if (CCTK_EQUALS(ODE_Method, "RK45CK"))
+ {
+ alpha_array = alpha_array_CK;
+ }
+ else
+ {
+ CCTK_WARN (0, "internal error");
+ }
cctkGH->cctk_time
= ((* Original_Time)
+ ((alpha_array[substep] - 1)
@@ -339,9 +362,22 @@ int MoL_ResetDeltaTime(CCTK_ARGUMENTS)
cctkGH->cctk_delta_time = 2.0/3.0*(*Original_Delta_Time);
}
}
- else if (CCTK_EQUALS(ODE_Method,"RK45"))
+ else if (CCTK_EQUALS(ODE_Method,"RK45") || CCTK_EQUALS(ODE_Method,"RK45CK"))
{
const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step);
+ const CCTK_REAL * alpha_array;
+ if (CCTK_EQUALS(ODE_Method, "RK45"))
+ {
+ alpha_array = alpha_array_F;
+ }
+ else if (CCTK_EQUALS(ODE_Method, "RK45CK"))
+ {
+ alpha_array = alpha_array_CK;
+ }
+ else
+ {
+ CCTK_WARN (0, "internal error");
+ }
cctkGH->cctk_delta_time
= ((alpha_array[substep + 1] - alpha_array[substep])
* (* Original_Delta_Time));
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-20 20:24:51 +0000