aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorhawke <hawke@578cdeb0-5ea1-4b81-8215-5a3b8777ee0b>2003-07-22 07:07:14 +0000
committerhawke <hawke@578cdeb0-5ea1-4b81-8215-5a3b8777ee0b>2003-07-22 07:07:14 +0000
commite3096d84ffb4b909977066f158db78eb45f4a80b (patch)
tree9e41096faa3de1f2af27660866aa5e4eea4cffed
parent06f6fb8bd438528724715988336a8b8ef86192ce (diff)
Add two methods:
RK3. The optimized version of the TVD RK3 solver. Requires no scratch space so is about as efficient as ICN, but third order. Generic method from a parameter table. By specifying the number of intermediate steps and the alpha and beta arrays, create your own method at parameter time. Not well (or at all) documented because it doesn't seem to work correctly at the moment. Some tidying of extraneous code as well. git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/MoL/trunk@29 578cdeb0-5ea1-4b81-8215-5a3b8777ee0b
-rw-r--r--param.ccl11
-rw-r--r--schedule.ccl7
-rw-r--r--src/ChangeType.c223
-rw-r--r--src/Counter.c2
-rw-r--r--src/GenericRK.c8
-rw-r--r--src/IndexArrays.c9
-rw-r--r--src/InitialCopy.c5
-rw-r--r--src/ParamCheck.c40
-rw-r--r--src/RK3.c348
-rw-r--r--src/RKCoefficients.c58
-rw-r--r--src/Registration.c32
-rw-r--r--src/SetTime.c23
-rw-r--r--src/make.code.defn1
13 files changed, 565 insertions, 202 deletions
diff --git a/param.ccl b/param.ccl
index 0bf2bcc..0aeae83 100644
--- a/param.ccl
+++ b/param.ccl
@@ -86,12 +86,14 @@ KEYWORD ODE_Method "The ODE method use by MoL to do time integration"
"ICN" :: "Iterative Crank Nicholson"
"ICN-avg" :: "Iterative Crank Nicholson with averaging"
"RK2" :: "Efficient RK2"
+ "RK3" :: "Efficient RK3"
} "ICN"
KEYWORD Generic_Type "If using the generic method, which sort"
{
"RK" :: "One of the standard TVD Runga-Kutta methods"
"ICN" :: "Iterative Crank Nicholson as a generic method"
+ "Table" :: "Given from the generic method descriptor parameter"
} "RK"
CCTK_INT MoL_Intermediate_Steps "Number of intermediate steps taken by the ODE method"
@@ -111,3 +113,12 @@ CCTK_REAL MoL_Tiny "Effective local machine zero; required by generic solvers"
BOOLEAN initial_data_is_crap "If the initial data routine fails to set up the previous time levels, copy the current backwards"
{
} "no"
+
+# The default for this parameter corresponds to generic RK2
+STRING Generic_Method_Descriptor "A string used to create a table containing the description of the generic method"
+{
+ ".*" :: "Should contain the Alpha and Beta arrays, and the number of intermediate steps"
+} "GenericIntermediateSteps = 2 \
+ GenericAlphaCoeffs = { 1.0 0.0 0.5 0.5 } \
+ GenericBetaCoeffs = { 1.0 0.5 }"
+
diff --git a/schedule.ccl b/schedule.ccl
index 4321dcf..fe74af9 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -259,6 +259,13 @@ else if (CCTK_Equals(ODE_Method,"RK2"))
LANG: C
} "Updates calculated with the efficient Runge-Kutta 2 method"
}
+else if (CCTK_Equals(ODE_Method,"RK3"))
+{
+ schedule MoL_RK3Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep
+ {
+ LANG: C
+ } "Updates calculated with the efficient Runge-Kutta 3 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/ChangeType.c b/src/ChangeType.c
index 9361f3a..2f1d926 100644
--- a/src/ChangeType.c
+++ b/src/ChangeType.c
@@ -5,7 +5,7 @@
@desc
The external functions called (via function aliasing) by physics
thorns to tell MoL that they want these GFs to be treated as a
- different type to the original declaration
+ different type to the original declaration.
@enddesc
@version $Header$
@@*/
@@ -48,33 +48,6 @@ CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex);
CCTK_INT MoL_ChangeToNone(CCTK_INT RemoveIndex);
-CCTK_INT MoL_ChangeToEvolvedComplex(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex);
-
-CCTK_INT MoL_ChangeToConstrainedComplex(CCTK_INT ConstrainedIndex);
-
-CCTK_INT MoL_ChangeToSaveAndRestoreComplex(CCTK_INT SandRIndex);
-
-CCTK_INT MoL_ChangeToNoneComplex(CCTK_INT RemoveIndex);
-
-
-CCTK_INT MoL_ChangeToEvolvedArray(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex);
-
-CCTK_INT MoL_ChangeToConstrainedArray(CCTK_INT ConstrainedIndex);
-
-CCTK_INT MoL_ChangeToSaveAndRestoreArray(CCTK_INT SandRIndex);
-
-CCTK_INT MoL_ChangeToNoneArray(CCTK_INT RemoveIndex);
-
-
-CCTK_INT MoL_ChangeToEvolvedComplexArray(CCTK_INT EvolvedIndex,
- CCTK_INT RHSIndex);
-
-CCTK_INT MoL_ChangeToConstrainedComplexArray(CCTK_INT ConstrainedIndex);
-
-CCTK_INT MoL_ChangeToSaveAndRestoreComplexArray(CCTK_INT SandRIndex);
-
-CCTK_INT MoL_ChangeToNoneComplexArray(CCTK_INT RemoveIndex);
-
/********************************************************************
********************* Local Data *****************************
********************************************************************/
@@ -145,7 +118,7 @@ CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
case MOL_UNKNOWN_VARTYPE:
{
- timelevs = CCTK_NumTimeLevelsFromVarI(EvolvedIndex);
+ timelevs = CCTK_MaxTimeLevelsVI(EvolvedIndex);
if (timelevs < 1)
{
CCTK_VWarn(1,__LINE__,__FILE__,"MoL","Warning whilst trying to change variable index %i to evolved.", EvolvedIndex);
@@ -156,7 +129,7 @@ CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
CCTK_VWarn(1,__LINE__,__FILE__,"MoL","Warning whilst trying to change variable index %i to evolved.", EvolvedIndex);
CCTK_WARN(0, "The index passed only has a single timelevel.");
}
- timelevs = CCTK_NumTimeLevelsFromVarI(RHSIndex);
+ timelevs = CCTK_MaxTimeLevelsVI(RHSIndex);
if (timelevs < 1) {
CCTK_VWarn(1,__LINE__,__FILE__,"MoL","Warning whilst trying to change variable index %i to evolved (RHS GF).", RHSIndex);
CCTK_WARN(0, "The index passed does not correspond to a GF.");
@@ -182,7 +155,7 @@ CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
case MOL_CONSTRAINED_VARTYPE:
{
- timelevs = CCTK_NumTimeLevelsFromVarI(RHSIndex);
+ timelevs = CCTK_MaxTimeLevelsVI(RHSIndex);
if (timelevs < 1) {
CCTK_VWarn(1,__LINE__,__FILE__,"MoL","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.");
@@ -206,7 +179,7 @@ CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
case MOL_SANDR_VARTYPE:
{
- timelevs = CCTK_NumTimeLevelsFromVarI(RHSIndex);
+ timelevs = CCTK_MaxTimeLevelsVI(RHSIndex);
if (timelevs < 1) {
CCTK_VWarn(1,__LINE__,__FILE__,"MoL","Warning whilst trying to change variable index %i to evolved type from save and restore.", RHSIndex);
CCTK_WARN(0, "The RHS index passed does not correspond to a GF.");
@@ -263,6 +236,18 @@ CCTK_INT MoL_ChangeToConstrained(CCTK_INT ConstrainedIndex)
vartype = 0;
usedindex = -1;
+
+ timelevs = CCTK_MaxTimeLevelsVI(ConstrainedIndex);
+ if (timelevs < 1)
+ {
+ CCTK_VWarn(1,__LINE__,__FILE__,"MoL","Warning whilst trying to change variable index %i to evolved.", ConstrainedIndex);
+ CCTK_WARN(0, "The index passed does not correspond to a GF.");
+ }
+ else if (timelevs == 1)
+ {
+ CCTK_VWarn(1,__LINE__,__FILE__,"MoL","Warning whilst trying to change variable index %i to evolved.", ConstrainedIndex);
+ CCTK_WARN(0, "The index passed only has a single timelevel.");
+ }
for (index = 0; (index < MoLNumEvolvedVariables)&&(!vartype); index++)
{
@@ -300,17 +285,6 @@ CCTK_INT MoL_ChangeToConstrained(CCTK_INT ConstrainedIndex)
case MOL_UNKNOWN_VARTYPE:
{
- timelevs = CCTK_NumTimeLevelsFromVarI(ConstrainedIndex);
- if (timelevs < 1)
- {
- CCTK_VWarn(1,__LINE__,__FILE__,"MoL","Warning whilst trying to change variable index %i to evolved.", ConstrainedIndex);
- CCTK_WARN(0, "The index passed does not correspond to a GF.");
- }
- else if (timelevs == 1)
- {
- CCTK_VWarn(1,__LINE__,__FILE__,"MoL","Warning whilst trying to change variable index %i to evolved.", ConstrainedIndex);
- CCTK_WARN(0, "The index passed only has a single timelevel.");
- }
ConstrainedVariableIndex[MoLNumConstrainedVariables] = ConstrainedIndex;
MoLNumConstrainedVariables++;
break;
@@ -319,6 +293,14 @@ CCTK_INT MoL_ChangeToConstrained(CCTK_INT ConstrainedIndex)
case MOL_EVOLVED_VARTYPE:
{
+ for (index = usedindex; index < MoLNumEvolvedVariables - 1; index++)
+ {
+ EvolvedVariableIndex[index] = EvolvedVariableIndex[index+1];
+ RHSVariableIndex[index] = RHSVariableIndex[index+1];
+ }
+ MoLNumEvolvedVariables--;
+ ConstrainedVariableIndex[index] = ConstrainedIndex;
+ MoLNumConstrainedVariables++;
break;
}
@@ -378,6 +360,13 @@ CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex)
vartype = 0;
usedindex = -1;
+
+ timelevs = CCTK_MaxTimeLevelsVI(SandRIndex);
+ if (timelevs < 1)
+ {
+ CCTK_VWarn(1,__LINE__,__FILE__,"MoL","Warning whilst trying to change variable index %i to save and restore.", SandRIndex);
+ CCTK_WARN(0, "The index passed does not correspond to a GF.");
+ }
for (index = 0; (index < MoLNumEvolvedVariables)&&(!vartype); index++)
{
@@ -415,12 +404,6 @@ CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex)
case MOL_UNKNOWN_VARTYPE:
{
- timelevs = CCTK_NumTimeLevelsFromVarI(SandRIndex);
- if (timelevs < 1)
- {
- CCTK_VWarn(1,__LINE__,__FILE__,"MoL","Warning whilst trying to change variable index %i to save and restore.", SandRIndex);
- CCTK_WARN(0, "The index passed does not correspond to a GF.");
- }
SandRVariableIndex[MoLNumSandRVariables] = SandRIndex;
MoLNumSandRVariables++;
break;
@@ -429,7 +412,14 @@ CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex)
case MOL_EVOLVED_VARTYPE:
{
-
+ for (index = usedindex; index < MoLNumEvolvedVariables - 1; index++)
+ {
+ EvolvedVariableIndex[index] = EvolvedVariableIndex[index+1];
+ RHSVariableIndex[index] = RHSVariableIndex[index+1];
+ }
+ MoLNumEvolvedVariables--;
+ SandRVariableIndex[MoLNumSandRVariables] = SandRIndex;
+ MoLNumSandRVariables++;
break;
}
@@ -441,7 +431,7 @@ CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex)
ConstrainedVariableIndex[index] = ConstrainedVariableIndex[index+1];
}
MoLNumConstrainedVariables--;
- SandRVariableIndex[index] = SandRIndex;
+ SandRVariableIndex[MoLNumSandRVariables] = SandRIndex;
MoLNumSandRVariables++;
break;
}
@@ -574,137 +564,6 @@ CCTK_INT MoL_ChangeToNone(CCTK_INT RemoveIndex)
}
-/*
- Currently no-op test functions.
-*/
-
-
-CCTK_INT MoL_ChangeToEvolvedComplex(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
-{
- CCTK_INT dummy;
- dummy = EvolvedIndex;
- dummy = RHSIndex;
- return 0;
-}
-
-CCTK_INT MoL_ChangeToConstrainedComplex(CCTK_INT ConstrainedIndex)
-{
- CCTK_INT dummy;
- dummy = ConstrainedIndex;
- return 0;
-}
-
-CCTK_INT MoL_ChangeToSaveAndRestoreComplex(CCTK_INT SandRIndex)
-{
- CCTK_INT dummy;
- dummy = SandRIndex;
- return 0;
-}
-
-CCTK_INT MoL_ChangeToNoneComplex(CCTK_INT RemoveIndex)
-{
- CCTK_INT dummy;
- dummy = RemoveIndex;
- return 0;
-}
-
-
-CCTK_INT MoL_ChangeToEvolvedArray(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
-{
- CCTK_INT dummy;
- dummy = EvolvedIndex;
- dummy = RHSIndex;
- return 0;
-}
-
-CCTK_INT MoL_ChangeToConstrainedArray(CCTK_INT ConstrainedIndex)
-{
- CCTK_INT dummy;
- dummy = ConstrainedIndex;
- return 0;
-}
-
-CCTK_INT MoL_ChangeToSaveAndRestoreArray(CCTK_INT SandRIndex)
-{
- CCTK_INT dummy;
- dummy = SandRIndex;
- return 0;
-}
-
-CCTK_INT MoL_ChangeToNoneArray(CCTK_INT RemoveIndex)
-{
- CCTK_INT dummy;
- dummy = RemoveIndex;
- return 0;
-}
-
-
-CCTK_INT MoL_ChangeToEvolvedComplexArray(CCTK_INT EvolvedIndex,
- CCTK_INT RHSIndex)
-{
- CCTK_INT dummy;
- dummy = EvolvedIndex;
- dummy = RHSIndex;
- return 0;
-}
-
-CCTK_INT MoL_ChangeToConstrainedComplexArray(CCTK_INT ConstrainedIndex)
-{
- CCTK_INT dummy;
- dummy = ConstrainedIndex;
- return 0;
-}
-
-CCTK_INT MoL_ChangeToSaveAndRestoreComplexArray(CCTK_INT SandRIndex)
-{
- CCTK_INT dummy;
- dummy = SandRIndex;
- return 0;
-}
-
-CCTK_INT MoL_ChangeToNoneComplexArray(CCTK_INT RemoveIndex)
-{
- CCTK_INT dummy;
- dummy = RemoveIndex;
- return 0;
-}
-
-/*
- Fortran wrappers for the above functions.
- Should be replaced by using function aliasing eventually.
- */
-
-/*
-void CCTK_FCALL CCTK_FNAME(MoL_ChangeToEvolved)(int *ierr,
- CCTK_INT *EvolvedIndex,
- CCTK_INT *RHSIndex)
-{
- *ierr = MoL_ChangeToEvolved(*EvolvedIndex, *RHSIndex);
- return;
-}
-
-void CCTK_FCALL CCTK_FNAME(MoL_ChangeToConstrained)(int *ierr,
- CCTK_INT *ConstrainedIndex)
-{
- *ierr = MoL_ChangeToConstrained(*ConstrainedIndex);
- return;
-}
-
-void CCTK_FCALL CCTK_FNAME(MoL_ChangeToSaveAndRestore)(int *ierr,
- CCTK_INT *EvolvedIndex)
-{
- *ierr = MoL_ChangeToSaveAndRestore(*EvolvedIndex);
- return;
-}
-
-void CCTK_FCALL CCTK_FNAME(MoL_ChangeToNone)(int *ierr,
- CCTK_INT *EvolvedIndex)
-{
- *ierr = MoL_ChangeToNone(*EvolvedIndex);
- return;
-}
-*/
-
/********************************************************************
********************* Local Routines *************************
********************************************************************/
diff --git a/src/Counter.c b/src/Counter.c
index 2d7f5cd..d44278f 100644
--- a/src/Counter.c
+++ b/src/Counter.c
@@ -66,8 +66,6 @@ int MoL_SetCounter(CCTK_ARGUMENTS)
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
-
-/* CCTK_REAL *Var; */
*MoL_Intermediate_Step = MoL_Intermediate_Steps;
diff --git a/src/GenericRK.c b/src/GenericRK.c
index 60edf15..19ffd0b 100644
--- a/src/GenericRK.c
+++ b/src/GenericRK.c
@@ -150,6 +150,10 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS)
if (scratchstep)
{
/*
+ The following would work if all drivers considered
+ a vector group to have contiguous storage.
+ */
+ /*
ScratchVar = &ScratchSpace[(var * MoL_Num_Scratch_Levels +
scratchindex) * totalsize];
*/
@@ -191,6 +195,10 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS)
UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0,
EvolvedVariableIndex[var]);
/*
+ The following would work if all drivers considered
+ a vector group to have contiguous storage.
+ */
+ /*
ScratchVar = &ScratchSpace[(var * MoL_Num_Scratch_Levels +
MoL_Intermediate_Steps -
(*MoL_Intermediate_Step)) * totalsize];
diff --git a/src/IndexArrays.c b/src/IndexArrays.c
index daf0d68..8fbded3 100644
--- a/src/IndexArrays.c
+++ b/src/IndexArrays.c
@@ -232,6 +232,11 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS)
{
sprintf(infoline, "Generic Runge-Kutta %i",MoL_Intermediate_Steps);
}
+ else if (CCTK_EQUALS(Generic_Type,"Table"))
+ {
+ sprintf(infoline, "Generic method, options:\n %s\n",
+ Generic_Method_Descriptor);
+ }
else
{
CCTK_WARN(0, "Generic_Type not recognized!");
@@ -241,6 +246,10 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS)
{
sprintf(infoline, "Runge-Kutta 2");
}
+ else if (CCTK_EQUALS(ODE_Method,"RK3"))
+ {
+ sprintf(infoline, "Runge-Kutta 3");
+ }
else if (CCTK_EQUALS(ODE_Method,"ICN"))
{
sprintf(infoline, "Iterative Crank Nicholson with %i iterations",
diff --git a/src/InitialCopy.c b/src/InitialCopy.c
index 1980d5f..1092ddb 100644
--- a/src/InitialCopy.c
+++ b/src/InitialCopy.c
@@ -422,7 +422,10 @@ void MoL_FillAllLevels(CCTK_ARGUMENTS)
{
CurrentVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
EvolvedVariableIndex[var]);
- for (level = 1; level < CCTK_NumTimeLevelsFromVarI(EvolvedVariableIndex[var]); level++)
+ for (level = 1;
+ level < CCTK_ActiveTimeLevelsVI(cctkGH,
+ EvolvedVariableIndex[var]);
+ level++)
{
PreviousVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, level,
EvolvedVariableIndex[var]);
diff --git a/src/ParamCheck.c b/src/ParamCheck.c
index 2ade916..01e6cbd 100644
--- a/src/ParamCheck.c
+++ b/src/ParamCheck.c
@@ -12,6 +12,9 @@
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
static const char *rcsid = "$Header$";
CCTK_FILEVERSION(AlphaThorns_MoL_ParamCheck_c);
@@ -62,19 +65,54 @@ int MoL_ParamCheck(CCTK_ARGUMENTS)
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
+ CCTK_INT options_table, ierr, GenericIntermediateSteps;
+
if (CCTK_Equals(ODE_Method, "Generic"))
{
if (MoL_Num_Scratch_Levels < MoL_Intermediate_Steps - 1)
{
CCTK_PARAMWARN("When using a generic solver the number of scratch levels must be at least the number of intermediate steps - 1");
}
+ if (CCTK_Equals(Generic_Type, "Table"))
+ {
+ options_table =
+ Util_TableCreateFromString(Generic_Method_Descriptor);
+ if (options_table < 0)
+ {
+ CCTK_WARN(0, "Failed to create table from Generic_Method_Descriptor!");
+ }
+ ierr = Util_TableGetInt(options_table,
+ &GenericIntermediateSteps,
+ "GenericIntermediateSteps");
+ if (ierr < 1)
+ {
+ if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+ {
+ CCTK_WARN(0, "When using the generic table options you must set \"GenericIntermediateSteps\" in the options table");
+ }
+ else
+ {
+ CCTK_WARN(0, "Table error - check with Ian.");
+ }
+ }
+ if (MoL_Intermediate_Steps != GenericIntermediateSteps)
+ {
+ CCTK_PARAMWARN("The number of intermediate steps must equal the number specified in the table options!");
+ }
+ ierr = Util_TableDestroy(options_table);
+ }
}
if ( (CCTK_Equals(ODE_Method, "RK2"))&&(!(MoL_Intermediate_Steps == 2)) )
{
CCTK_PARAMWARN("When using the efficient RK2 evolver the number of intermediate steps must be 2");
}
-
+
+ if ( (CCTK_Equals(ODE_Method, "RK3"))&&(!(MoL_Intermediate_Steps == 3)) )
+ {
+ CCTK_PARAMWARN("When using the efficient RK3 evolver the number of intermediate steps must be 3");
+ }
+
return 0;
}
diff --git a/src/RK3.c b/src/RK3.c
new file mode 100644
index 0000000..b636698
--- /dev/null
+++ b/src/RK3.c
@@ -0,0 +1,348 @@
+ /*@@
+ @file RK3.c
+ @date Tue Jul 22 00:38:47 2003
+ @author Ian Hawke
+ @desc
+ A specialized third order Runge-Kutta time integrator. This is
+ the integrator that Shu refers to as the optimal TVD third
+ order method (see reference in documentation).
+ @enddesc
+ @version $Header$
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "ExternalVariables.h"
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(AlphaThorns_MoL_RK3_c);
+
+/********************************************************************
+ ********************* Local Data Types ***********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ ***************** Scheduled Routine Prototypes *********************
+ ********************************************************************/
+
+void MoL_RK3Add(CCTK_ARGUMENTS);
+
+/********************************************************************
+ ********************* Other Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Local Data *****************************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* External Routines **********************
+ ********************************************************************/
+
+ /*@@
+ @routine MoL_RK3Add
+ @date Tue Jul 22 00:39:55 2003
+ @author Ian Hawke
+ @desc
+ Performs third order Runge-Kutta time integration.
+ @enddesc
+ @calls
+ @calledby
+ @history
+
+ @endhistory
+
+ @@*/
+
+void MoL_RK3Add(CCTK_ARGUMENTS)
+{
+
+ DECLARE_CCTK_ARGUMENTS;
+ DECLARE_CCTK_PARAMETERS;
+
+ cGroupDynamicData arraydata;
+ CCTK_INT groupindex, ierr;
+ CCTK_INT arraytotalsize, arraydim;
+
+ CCTK_INT index, var;
+ CCTK_INT totalsize;
+ CCTK_REAL *OldVar;
+ CCTK_REAL *UpdateVar;
+ CCTK_REAL *RHSVar;
+
+ /* FIXME */
+
+#ifdef MOLDOESCOMPLEX
+
+ CCTK_COMPLEX *OldComplexVar;
+ CCTK_COMPLEX *UpdateComplexVar;
+ CCTK_COMPLEX *RHSComplexVar;
+ CCTK_COMPLEX Complex_Delta_Time = CCTK_Cmplx(CCTK_DELTA_TIME, 0);
+ CCTK_COMPLEX Complex_Half = CCTK_Cmplx(0.5, 0);
+ CCTK_COMPLEX Complex_Third = CCTK_CmplxDiv(CCTK_Cmplx(1.0, 0),
+ CCTK_Cmplx(3.0,0));
+ CCTK_COMPLEX Complex_TwoThird = CCTK_CmplxDiv(CCTK_Cmplx(2.0, 0),
+ CCTK_Cmplx(3.0,0));
+ CCTK_COMPLEX Complex_Quarter = CCTK_Cmplx(0.25, 0);
+ CCTK_COMPLEX Complex_ThreeQuarter = CCTK_Cmplx(0.75, 0);
+
+#endif
+
+#ifdef MOLDEBUG
+ printf("Inside RK3.\nStep %d.\nRefinement %d.\nTimestep %g.\nSpacestep %g.\nTime %g\n",
+ MoL_Intermediate_Steps - *MoL_Intermediate_Step + 1,
+ *cctk_levfac,
+ CCTK_DELTA_TIME,
+ CCTK_DELTA_SPACE(0),
+ cctk_time);
+#endif
+
+ totalsize = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];
+
+ switch (*MoL_Intermediate_Step)
+ {
+
+ case 3:
+ {
+ for (var = 0; var < MoLNumEvolvedVariables; var++)
+ {
+ UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
+ EvolvedVariableIndex[var]);
+ RHSVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
+ RHSVariableIndex[var]);
+
+ for (index = 0; index < totalsize; index++)
+ {
+ UpdateVar[index] += CCTK_DELTA_TIME * RHSVar[index];
+ }
+ }
+
+ for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
+ {
+ UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
+ EvolvedArrayVariableIndex[var]);
+ RHSVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
+ RHSArrayVariableIndex[var]);
+ groupindex = CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]);
+ ierr = CCTK_GroupDynamicData(cctkGH, groupindex,
+ &arraydata);
+ if (ierr)
+ {
+ CCTK_VWarn(0, __LINE__, __FILE__, "MoL",
+ "The driver does not return group information for group '%s'.",
+ CCTK_GroupName(groupindex));
+ }
+ arraytotalsize = 1;
+ for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
+ {
+ arraytotalsize *= arraydata.lsh[arraydim];
+ }
+
+ for (index = 0; index < arraytotalsize; index++)
+ {
+ UpdateVar[index] += CCTK_DELTA_TIME * RHSVar[index];
+ }
+ }
+
+ /* FIXME */
+
+#ifdef MOLDOESCOMPLEX
+
+ for (var = 0; var < MoLNumEvolvedComplexVariables; var++)
+ {
+ UpdateComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0,
+ EvolvedComplexVariableIndex[var]);
+ RHSComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0,
+ RHSComplexVariableIndex[var]);
+
+ for (index = 0; index < totalsize; index++)
+ {
+ UpdateComplexVar[index] = CCTK_CmplxAdd(UpdateComplexVar[index],
+ CCTK_CmplxMul(Complex_Delta_Time,
+ RHSComplexVar[index]));
+ }
+ }
+
+#endif
+
+ break;
+ }
+
+ case 2:
+ {
+ 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,
+ RHSVariableIndex[var]);
+
+ for (index = 0; index < totalsize; index++)
+ {
+ UpdateVar[index] = 0.25 * (3*OldVar[index] +
+ UpdateVar[index]) +
+ CCTK_DELTA_TIME * RHSVar[index];
+ }
+ }
+
+ for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
+ {
+ OldVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 1,
+ EvolvedArrayVariableIndex[var]);
+ UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
+ EvolvedArrayVariableIndex[var]);
+ RHSVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
+ RHSArrayVariableIndex[var]);
+ groupindex = CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]);
+ ierr = CCTK_GroupDynamicData(cctkGH, groupindex,
+ &arraydata);
+ if (ierr)
+ {
+ CCTK_VWarn(0, __LINE__, __FILE__, "MoL",
+ "The driver does not return group information for group '%s'.",
+ CCTK_GroupName(groupindex));
+ }
+ arraytotalsize = 1;
+ for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
+ {
+ arraytotalsize *= arraydata.lsh[arraydim];
+ }
+
+ for (index = 0; index < arraytotalsize; index++)
+ {
+ UpdateVar[index] = 0.25*(3*OldVar[index] + UpdateVar[index]) +
+ CCTK_DELTA_TIME * RHSVar[index];
+ }
+ }
+
+ /* FIXME */
+
+#ifdef MOLDOESCOMPLEX
+
+ for (var = 0; var < MoLNumEvolvedComplexVariables; var++)
+ {
+ OldComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 1,
+ EvolvedComplexVariableIndex[var]);
+ UpdateComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0,
+ EvolvedComplexVariableIndex[var]);
+ RHSComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0,
+ RHSComplexVariableIndex[var]);
+
+ for (index = 0; index < totalsize; index++)
+ {
+ UpdateComplexVar[index] =
+ CCTK_CmplxAdd(CCTK_CmplxAdd(CCTK_CmplxMul(Complex_ThreeQuarter,
+ OldComplexVar[index]),
+ CCTK_CmplxMul(ComplexQuarter,
+ UpdateComplexVar[index])),
+ CCTK_CmplxMul(CCTK_CmplxMul(Complex_Delta_Time,
+ RHSComplexVar[index])));
+ }
+ }
+
+#endif
+
+ break;
+ }
+ case 1:
+ {
+ 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,
+ RHSVariableIndex[var]);
+
+ for (index = 0; index < totalsize; index++)
+ {
+ UpdateVar[index] = (OldVar[index] + 2*UpdateVar[index])/3
+ + CCTK_DELTA_TIME * RHSVar[index];
+ }
+ }
+
+ for (var = 0; var < MoLNumEvolvedArrayVariables; var++)
+ {
+ OldVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 1,
+ EvolvedArrayVariableIndex[var]);
+ UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
+ EvolvedArrayVariableIndex[var]);
+ RHSVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0,
+ RHSArrayVariableIndex[var]);
+
+ groupindex = CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]);
+ ierr = CCTK_GroupDynamicData(cctkGH, groupindex,
+ &arraydata);
+ if (ierr)
+ {
+ CCTK_VWarn(0, __LINE__, __FILE__, "MoL",
+ "The driver does not return group information for group '%s'.",
+ CCTK_GroupName(groupindex));
+ }
+ arraytotalsize = 1;
+ for (arraydim = 0; arraydim < arraydata.dim; arraydim++)
+ {
+ arraytotalsize *= arraydata.lsh[arraydim];
+ }
+
+ for (index = 0; index < arraytotalsize; index++)
+ {
+ UpdateVar[index] = (OldVar[index] + 2*UpdateVar[index])/3
+ + CCTK_DELTA_TIME * RHSVar[index];
+ }
+ }
+
+ /* FIXME */
+
+#ifdef MOLDOESCOMPLEX
+
+ for (var = 0; var < MoLNumEvolvedComplexVariables; var++)
+ {
+ OldComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 1,
+ EvolvedComplexVariableIndex[var]);
+ UpdateComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0,
+ EvolvedComplexVariableIndex[var]);
+ RHSComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0,
+ RHSComplexVariableIndex[var]);
+
+ for (index = 0; index < totalsize; index++)
+ {
+ UpdateComplexVar[index] =
+ CCTK_CmplxAdd(CCTK_CmplxAdd(CCTK_CmplxMul(Complex_Third,
+ OldComplexVar[index]),
+ CCTK_CmplxMul(ComplexTwoThird,
+ UpdateComplexVar[index])),
+ CCTK_CmplxMul(Complex_Delta_Time,
+ RHSComplexVar[index]));
+ }
+ }
+
+#endif
+
+ break;
+ }
+ default:
+ {
+ CCTK_WARN(0, "RK3 expects MoL_Intermediate_Step to be in [1,3]. This should be caught at ParamCheck - bug Ian!");
+ break;
+ }
+
+ }
+
+ return;
+
+}
+
+/********************************************************************
+ ********************* Local Routines *************************
+ ********************************************************************/
diff --git a/src/RKCoefficients.c b/src/RKCoefficients.c
index 635df7a..cae00b2 100644
--- a/src/RKCoefficients.c
+++ b/src/RKCoefficients.c
@@ -14,6 +14,9 @@
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
static const char *rcsid = "$Header$";
CCTK_FILEVERSION(AlphaThorns_MoL_RKCoefficients_c);
@@ -68,6 +71,8 @@ int MoL_SetupRKCoefficients(CCTK_ARGUMENTS)
CCTK_INT i, j;
+ CCTK_INT ierr, options_table;
+
if (CCTK_Equals(Generic_Type,"ICN"))
{
for (i = 0; i < MoL_Intermediate_Steps; i++)
@@ -145,6 +150,59 @@ int MoL_SetupRKCoefficients(CCTK_ARGUMENTS)
CCTK_WARN(0, "RKCoefficients cannot do generic RK methods with MoL_Intermediate_Steps greater than 4");
}
}
+ else if (CCTK_Equals(Generic_Type,"Table"))
+ {
+ if (MoL_Num_Scratch_Levels < MoL_Intermediate_Steps - 1)
+ {
+ CCTK_WARN(0, "For generic methods, MoL_Num_Scratch_Levels should be at least MoL_Intermediate_Steps - 1");
+ }
+ options_table =
+ Util_TableCreateFromString(Generic_Method_Descriptor);
+ if (options_table < 0)
+ {
+ CCTK_WARN(0, "Failed to create table from Generic_Method_Descriptor!");
+ }
+ ierr = Util_TableGetRealArray(options_table,
+ (MoL_Num_Scratch_Levels + 1) *
+ MoL_Intermediate_Steps,
+ RKAlphaCoefficients,
+ "GenericAlphaCoeffs");
+ if (ierr < (MoL_Num_Scratch_Levels + 1) * MoL_Intermediate_Steps )
+ {
+ if (ierr >= 0)
+ {
+ CCTK_WARN(0, "Insufficient elements in the specified GenericAlphaCoeffs array");
+ }
+ else if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+ {
+ CCTK_WARN(0, "When using the generic table options you must set \"GenericAlphaCoeffs\" in the options table");
+ }
+ else
+ {
+ CCTK_WARN(0, "Table error - check with Ian.");
+ }
+ }
+ ierr = Util_TableGetRealArray(options_table,
+ MoL_Intermediate_Steps,
+ RKBetaCoefficients,
+ "GenericBetaCoeffs");
+ if (ierr < MoL_Intermediate_Steps)
+ {
+ if (ierr >= 0)
+ {
+ CCTK_WARN(0, "Insufficient elements in the specified GenericBetaCoeffs array");
+ }
+ else if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+ {
+ CCTK_WARN(0, "When using the generic table options you must set \"GenericBetaCoeffs\" in the options table");
+ }
+ else
+ {
+ CCTK_WARN(0, "Table error - check with Ian.");
+ }
+ }
+ ierr = Util_TableDestroy(options_table);
+ }
else
{
CCTK_WARN(0, "RKCoefficients does not recognize the value of Generic_Type");
diff --git a/src/Registration.c b/src/Registration.c
index 2dd4767..2d851c3 100644
--- a/src/Registration.c
+++ b/src/Registration.c
@@ -735,8 +735,8 @@ CCTK_INT MoL_RegisterEvolvedReal(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
CCTK_VarName(RHSIndex));
}
- numtimelevs1 = CCTK_NumTimeLevelsFromVarI(EvolvedIndex);
- numtimelevs2 = CCTK_NumTimeLevelsFromVarI(RHSIndex);
+ numtimelevs1 = CCTK_MaxTimeLevelsVI(EvolvedIndex);
+ numtimelevs2 = CCTK_MaxTimeLevelsVI(RHSIndex);
if ( (numtimelevs1 < 0) || (numtimelevs2 < 0) )
{
@@ -864,7 +864,7 @@ CCTK_INT MoL_RegisterConstrainedReal(CCTK_INT ConstrainedIndex)
CCTK_VarName(ConstrainedIndex));
}
- numtimelevs = CCTK_NumTimeLevelsFromVarI(ConstrainedIndex);
+ numtimelevs = CCTK_MaxTimeLevelsVI(ConstrainedIndex);
if (numtimelevs < 1) {
@@ -962,7 +962,7 @@ CCTK_INT MoL_RegisterSaveAndRestoreReal(CCTK_INT SandRIndex)
CCTK_VarName(SandRIndex));
}
- numtimelevs = CCTK_NumTimeLevelsFromVarI(SandRIndex);
+ numtimelevs = CCTK_MaxTimeLevelsVI(SandRIndex);
if (numtimelevs < 1) {
@@ -1164,8 +1164,8 @@ CCTK_INT MoL_RegisterEvolvedComplex(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
CCTK_VarName(RHSIndex));
}
- numtimelevs1 = CCTK_NumTimeLevelsFromVarI(EvolvedIndex);
- numtimelevs2 = CCTK_NumTimeLevelsFromVarI(RHSIndex);
+ numtimelevs1 = CCTK_MaxTimeLevelsVI(EvolvedIndex);
+ numtimelevs2 = CCTK_MaxTimeLevelsVI(RHSIndex);
if ( (numtimelevs1 < 0) || (numtimelevs2 < 0) )
{
@@ -1279,7 +1279,7 @@ CCTK_INT MoL_RegisterConstrainedComplex(CCTK_INT ConstrainedIndex)
CCTK_VarName(ConstrainedIndex));
}
- numtimelevs = CCTK_NumTimeLevelsFromVarI(ConstrainedIndex);
+ numtimelevs = CCTK_MaxTimeLevelsVI(ConstrainedIndex);
if (numtimelevs < 1) {
@@ -1373,7 +1373,7 @@ CCTK_INT MoL_RegisterSaveAndRestoreComplex(CCTK_INT SandRIndex)
CCTK_VarName(SandRIndex));
}
- numtimelevs = CCTK_NumTimeLevelsFromVarI(SandRIndex);
+ numtimelevs = CCTK_MaxTimeLevelsVI(SandRIndex);
if (numtimelevs < 1) {
@@ -1574,8 +1574,8 @@ CCTK_INT MoL_RegisterEvolvedArray(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
CCTK_VarName(RHSIndex));
}
- numtimelevs1 = CCTK_NumTimeLevelsFromVarI(EvolvedIndex);
- numtimelevs2 = CCTK_NumTimeLevelsFromVarI(RHSIndex);
+ numtimelevs1 = CCTK_MaxTimeLevelsVI(EvolvedIndex);
+ numtimelevs2 = CCTK_MaxTimeLevelsVI(RHSIndex);
if ( (numtimelevs1 < 0) || (numtimelevs2 < 0) )
{
@@ -1683,7 +1683,7 @@ CCTK_INT MoL_RegisterConstrainedArray(CCTK_INT ConstrainedIndex)
CCTK_VarName(ConstrainedIndex));
}
- numtimelevs = CCTK_NumTimeLevelsFromVarI(ConstrainedIndex);
+ numtimelevs = CCTK_MaxTimeLevelsVI(ConstrainedIndex);
if (numtimelevs < 1) {
@@ -1771,7 +1771,7 @@ CCTK_INT MoL_RegisterSaveAndRestoreArray(CCTK_INT SandRIndex)
CCTK_VarName(SandRIndex));
}
- numtimelevs = CCTK_NumTimeLevelsFromVarI(SandRIndex);
+ numtimelevs = CCTK_MaxTimeLevelsVI(SandRIndex);
if (numtimelevs < 1) {
@@ -1975,8 +1975,8 @@ CCTK_INT MoL_RegisterEvolvedComplexArray(CCTK_INT EvolvedIndex,
CCTK_VarName(RHSIndex));
}
- numtimelevs1 = CCTK_NumTimeLevelsFromVarI(EvolvedIndex);
- numtimelevs2 = CCTK_NumTimeLevelsFromVarI(RHSIndex);
+ numtimelevs1 = CCTK_MaxTimeLevelsVI(EvolvedIndex);
+ numtimelevs2 = CCTK_MaxTimeLevelsVI(RHSIndex);
if ( (numtimelevs1 < 0) || (numtimelevs2 < 0) )
{
@@ -2090,7 +2090,7 @@ CCTK_INT MoL_RegisterConstrainedComplexArray(CCTK_INT ConstrainedIndex)
CCTK_VarName(ConstrainedIndex));
}
- numtimelevs = CCTK_NumTimeLevelsFromVarI(ConstrainedIndex);
+ numtimelevs = CCTK_MaxTimeLevelsVI(ConstrainedIndex);
if (numtimelevs < 1) {
@@ -2184,7 +2184,7 @@ CCTK_INT MoL_RegisterSaveAndRestoreComplexArray(CCTK_INT SandRIndex)
CCTK_VarName(SandRIndex));
}
- numtimelevs = CCTK_NumTimeLevelsFromVarI(SandRIndex);
+ numtimelevs = CCTK_MaxTimeLevelsVI(SandRIndex);
if (numtimelevs < 1) {
diff --git a/src/SetTime.c b/src/SetTime.c
index ba76ffe..5c2a213 100644
--- a/src/SetTime.c
+++ b/src/SetTime.c
@@ -184,6 +184,18 @@ int MoL_ResetTime(CCTK_ARGUMENTS)
cctkGH->cctk_time = (*Original_Time);
}
}
+ else if (CCTK_EQUALS(ODE_Method,"RK3"))
+ {
+ if (*MoL_Intermediate_Step == 2)
+ {
+ cctkGH->cctk_time = (*Original_Time);
+ }
+ else if (*MoL_Intermediate_Step == 1)
+ {
+ cctkGH->cctk_time = (*Original_Time) +
+ 0.5*(*Original_Delta_Time)/cctkGH->cctk_timefac;
+ }
+ }
#ifdef MOLDEBUG
printf("MoL has once more reset t (%d): %f.\n", *MoL_Intermediate_Step, cctkGH->cctk_time);
fflush(stdout);
@@ -250,6 +262,17 @@ int MoL_ResetDeltaTime(CCTK_ARGUMENTS)
cctkGH->cctk_delta_time = 0.5*(*Original_Delta_Time);
}
}
+ else if (CCTK_EQUALS(ODE_Method,"RK3"))
+ {
+ if (*MoL_Intermediate_Step == 2)
+ {
+ cctkGH->cctk_delta_time = 0.25*(*Original_Delta_Time);
+ }
+ else if (*MoL_Intermediate_Step == 1)
+ {
+ cctkGH->cctk_delta_time = 2.0/3.0*(*Original_Delta_Time);
+ }
+ }
#ifdef MOLDEBUG
printf("MoL has once more reset dt (%d): %f.\n",
*MoL_Intermediate_Step,
diff --git a/src/make.code.defn b/src/make.code.defn
index 1afaae5..48a4349 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -10,6 +10,7 @@ SRCS = ChangeType.c \
InitialCopy.c \
ParamCheck.c \
RK2.c \
+ RK3.c \
Registration.c \
RKCoefficients.c \
SandR.c \