From a4ed7cfd73db6ef8418ffb47f4237cc529f8c710 Mon Sep 17 00:00:00 2001 From: eschnett Date: Tue, 22 Jan 2013 21:00:28 +0000 Subject: MoL Update New integrator Euler. This is an explicit, first-order method, mostly useful for debugging. New integrators AB (Adams-Bashforth) with various orders. These are explicit integrators using several past timelevels to provide higher-order integration with a single RHS evaluation each. Introduce LinearCombination, a generic routine to calculate linear combinations of grid functions. This simplifies existing code, and can be overloaded if MoL should run on a device (e.g. with OpenCL). Replace cctk_lsh with cctk_ash where necessary. git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/MoL/trunk@190 578cdeb0-5ea1-4b81-8215-5a3b8777ee0b --- interface.ccl | 26 ++++- param.ccl | 28 ++++- schedule.ccl | 16 ++- src/AB.c | 234 ++++++++++++++++++++++++++++++++++++++++ src/Euler.c | 116 ++++++++++++++++++++ src/GenericRK.c | 2 +- src/ICN.c | 8 +- src/IndexArrays.c | 37 +++++++ src/InitialCopy.c | 50 +++------ src/Operators.c | 312 +++++++++++++++++++++++++++++++++++++++++++++++++++++ src/Operators.h | 15 +++ src/ParamCheck.c | 12 +++ src/RK2-MR-2_1.c | 2 +- src/RK2.c | 163 +++++++++------------------- src/RK3.c | 250 +++++++++++++----------------------------- src/RK4-MR-2_1.c | 2 +- src/RK4-RK2.c | 2 +- src/RK4.c | 111 ++++++++----------- src/RK45.c | 2 +- src/RK65.c | 2 +- src/RK87.c | 172 +++++++++++++---------------- src/SandR.c | 2 +- src/SetTime.c | 14 +++ src/make.code.defn | 15 +-- 24 files changed, 1073 insertions(+), 520 deletions(-) create mode 100644 src/AB.c create mode 100644 src/Euler.c create mode 100644 src/Operators.c create mode 100644 src/Operators.h diff --git a/interface.ccl b/interface.ccl index d2d14f9..94a4bb7 100644 --- a/interface.ccl +++ b/interface.ccl @@ -206,10 +206,28 @@ PROVIDES FUNCTION MoLNumIntegratorSubsteps WITH MoL_NumIntegratorSubsteps \ # PROVIDES FUNCTION MoLChangeToNoneComplexArray WITH \ # MoL_ChangeToNoneComplexArray LANGUAGE C -################################################## -### Functions provided by devices ### -################################################## - +############################################################## +### Aliased functions that can be provided by other thorns ### +### to override low-level grid variable operations. ### +############################################################## + +CCTK_INT FUNCTION LinearCombination \ + (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN var, \ + CCTK_REAL IN scale, \ + CCTK_INT ARRAY IN srcs, \ + CCTK_INT ARRAY IN tls, \ + CCTK_REAL ARRAY IN facts, \ + CCTK_INT IN nsrcs) +USES FUNCTION LinearCombination + +CCTK_INT FUNCTION Accelerator_NotifyDataModified \ + (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT ARRAY IN variables, \ + CCTK_INT ARRAY IN tls, \ + CCTK_INT IN nvariables, \ + CCTK_INT IN on_device) +USES FUNCTION Accelerator_NotifyDataModified private: diff --git a/param.ccl b/param.ccl index 43a1fd4..cce525f 100644 --- a/param.ccl +++ b/param.ccl @@ -91,6 +91,7 @@ KEYWORD ODE_Method "The ODE method use by MoL to do time integration" "Generic" :: "Generic Shu-Osher Runga-Kutta type" "ICN" :: "Iterative Crank Nicholson" "ICN-avg" :: "Iterative Crank Nicholson with averaging" + "Euler" :: "Euler" "RK2" :: "Efficient RK2" "RK3" :: "Efficient RK3" "RK4" :: "Efficient RK4" @@ -98,6 +99,7 @@ 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" + "AB" :: "Adams-Bashforth" "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" @@ -120,6 +122,19 @@ BOOLEAN ICN_avg_swapped "Use swapped averages in ICN method?" { } "no" +KEYWORD AB_Type "If using the the AB method, which sort" +{ + "1" :: "same as forward Euler" + "2" :: "second order" + "3" :: "third order" + "4" :: "fourth order" + "5" :: "fifth order" +} "1" + +BOOLEAN AB_initially_reduce_order "Reduce order of accuracy initially so that no past timelevels of initial data are required" +{ +} "yes" + CCTK_INT MoL_Intermediate_Steps "Number of intermediate steps taken by the ODE method" { 1:* :: "Anything greater than 1" @@ -152,6 +167,11 @@ BOOLEAN copy_ID_after_MoL_PostStep \ { } "no" +BOOLEAN set_ID_boundaries "Should boundaries be overwritten (via synchronization, prolongation, boundary conditions) by MoL?" +{ +} "yes" + + # 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" @@ -221,6 +241,8 @@ KEYWORD verbose "How verbose should MoL be?" "extreme" :: "Everything you never wanted to know" } "normal" -BOOLEAN set_ID_boundaries "Should boundaries be overwritten (via synchronization, prolongation, boundary conditions) by MoL?" -{ -} "yes" + + +shares: Cactus + +USES CCTK_REAL cctk_initial_time diff --git a/schedule.ccl b/schedule.ccl index 48283fd..21896e3 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -551,6 +551,13 @@ if (CCTK_Equals(ODE_Method,"Generic")) LANG: C } "Updates calculated with a generic method" } +else if (CCTK_Equals(ODE_Method,"Euler")) +{ + schedule MoL_EulerAdd AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep + { + LANG: C + } "Updates calculated with the Euler method" +} else if (CCTK_Equals(ODE_Method,"RK2")) { schedule MoL_RK2Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep @@ -613,6 +620,13 @@ else if (CCTK_Equals(ODE_Method,"ICN-avg")) LANG: C } "Updates calculated with the averaging ICN method" } +else if (CCTK_Equals(ODE_Method,"AB")) +{ + schedule MoL_ABAdd AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep + { + LANG: C + } "Updates calculated with the Adams-Bashforth" +} 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 @@ -765,7 +779,7 @@ schedule MoL_ResetDeltaTime IN MoL_Step AFTER (MoL_PostStep MoL_PostStepModify) ### variables to their original state. ### ################################################## -schedule MoL_RestoreSandR IN MoL_Evolution AFTER MoL_FinishLoop +schedule MoL_RestoreSandR IN MoL_Evolution AFTER (MoL_ReduceAdaptiveError MoL_FinishLoop) { LANG: C } "Restoring the Save and Restore variables to the original state" diff --git a/src/AB.c b/src/AB.c new file mode 100644 index 0000000..ee7fdad --- /dev/null +++ b/src/AB.c @@ -0,0 +1,234 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include +#include + +#include "ExternalVariables.h" + + + +/* Coefficients taken from + * http://en.wikipedia.org/wiki/Linear_multistep_method, which cites + * (Hairer, Nørsett & Wanner 1993, §III.1; Butcher 2003, p. 103). + * + * Hairer, Ernst; Nørsett, Syvert Paul; Wanner, Gerhard (1993), + * Solving ordinary differential equations I: Nonstiff problems (2nd + * ed.), Berlin: Springer Verlag, ISBN 978-3-540-56670-0. + * + * Butcher, John C. (2003), Numerical Methods for Ordinary + * Differential Equations, John Wiley, ISBN 978-0-471-96758-3. + */ + +/* The following Mathematic expression (see + * http://en.wikipedia.org/wiki/Linear_multistep_method) calculates + * Adams-Bashforth coefficients for arbitrary orders: + * + * b = Table[ + * Table[(-1)^j/(j! (s - j - 1)!) Integrate[ + * Product[If[j == i, 1, u + i], {i, 0, s - 1}], {u, 0, 1}], {j, 0, + * s - 1}], {s, 1, 8}] + * + * where s is the desired order. The coefficients up to order 8 are: + * + * {1} + * + * {3/2, -(1/2)} + * + * {23/12, -(4/3), 5/12} + * + * {55/24, -(59/24), 37/24, -(3/8)} + * + * {1901/720, -(1387/360), 109/30, -(637/360), 251/720} + * + * {4277/1440, -(2641/480), 4991/720, -(3649/720), 959/480, -(95/288)} + * + * {198721/60480, -(18637/2520), 235183/20160, -(10754/945), + * 135713/20160, -(5603/2520), 19087/60480} + * + * {16083/4480, -(1152169/120960), 242653/13440, -(296053/13440), + * 2102243/120960, -(115747/13440), 32863/13440, -(5257/17280)} + */ + + + +static void order1 (CCTK_REAL* restrict const UpdateVar, + CCTK_REAL const* restrict *restrict const RHSVars, + CCTK_REAL const dt, + int const totalsize) +{ + CCTK_REAL const* restrict const RHSVar0 = RHSVars[0]; +#pragma omp parallel for + for (int index = 0; index < totalsize; index++) { + UpdateVar[index] += dt * RHSVar0[index]; + } +} + +static void order2 (CCTK_REAL* restrict const UpdateVar, + CCTK_REAL const* restrict *restrict const RHSVars, + CCTK_REAL const dt, + int const totalsize) +{ + CCTK_REAL const* restrict const RHSVar0 = RHSVars[0]; + CCTK_REAL const* restrict const RHSVar1 = RHSVars[1]; + CCTK_REAL const f0 = + (3.0/2.0) * dt; + CCTK_REAL const f1 = - (1.0/2.0) * dt; +#pragma omp parallel for + for (int index = 0; index < totalsize; index++) { + UpdateVar[index] += f0 * RHSVar0[index] + f1 * RHSVar1[index]; + } +} + +static void order3 (CCTK_REAL* restrict const UpdateVar, + CCTK_REAL const* restrict *restrict const RHSVars, + CCTK_REAL const dt, + int const totalsize) +{ + CCTK_REAL const* restrict const RHSVar0 = RHSVars[0]; + CCTK_REAL const* restrict const RHSVar1 = RHSVars[1]; + CCTK_REAL const* restrict const RHSVar2 = RHSVars[2]; + CCTK_REAL const f0 = + (23.0/12.0) * dt; + CCTK_REAL const f1 = - ( 4.0/ 3.0) * dt; + CCTK_REAL const f2 = + ( 5.0/12.0) * dt; +#pragma omp parallel for + for (int index = 0; index < totalsize; index++) { + UpdateVar[index] += + f0 * RHSVar0[index] + f1 * RHSVar1[index] + f2 * RHSVar2[index]; + } +} + +static void order4 (CCTK_REAL* restrict const UpdateVar, + CCTK_REAL const* restrict *restrict const RHSVars, + CCTK_REAL const dt, + int const totalsize) +{ + CCTK_REAL const* restrict const RHSVar0 = RHSVars[0]; + CCTK_REAL const* restrict const RHSVar1 = RHSVars[1]; + CCTK_REAL const* restrict const RHSVar2 = RHSVars[2]; + CCTK_REAL const* restrict const RHSVar3 = RHSVars[3]; + CCTK_REAL const f0 = + (55.0/24.0) * dt; + CCTK_REAL const f1 = - (59.0/24.0) * dt; + CCTK_REAL const f2 = + (37.0/24.0) * dt; + CCTK_REAL const f3 = - ( 3.0/ 8.0) * dt; +#pragma omp parallel for + for (int index = 0; index < totalsize; index++) { + UpdateVar[index] += + f0 * RHSVar0[index] + f1 * RHSVar1[index] + f2 * RHSVar2[index] + + f3 * RHSVar3[index]; + } +} + +static void order5 (CCTK_REAL* restrict const UpdateVar, + CCTK_REAL const* restrict *restrict const RHSVars, + CCTK_REAL const dt, + int const totalsize) +{ + CCTK_REAL const* restrict const RHSVar0 = RHSVars[0]; + CCTK_REAL const* restrict const RHSVar1 = RHSVars[1]; + CCTK_REAL const* restrict const RHSVar2 = RHSVars[2]; + CCTK_REAL const* restrict const RHSVar3 = RHSVars[3]; + CCTK_REAL const* restrict const RHSVar4 = RHSVars[4]; + CCTK_REAL const f0 = + (1901.0/720.0) * dt; + CCTK_REAL const f1 = - (1387.0/360.0) * dt; + CCTK_REAL const f2 = + ( 109.0/ 30.0) * dt; + CCTK_REAL const f3 = - ( 637.0/360.0) * dt; + CCTK_REAL const f4 = + ( 251.0/720.0) * dt; +#pragma omp parallel for + for (int index = 0; index < totalsize; index++) { + UpdateVar[index] += + f0 * RHSVar0[index] + f1 * RHSVar1[index] + f2 * RHSVar2[index] + + f3 * RHSVar3[index] + f4 * RHSVar4[index]; + } +} + +/* Array of function pointers */ +static +void (* const orders[]) (CCTK_REAL* restrict const UpdateVar, + CCTK_REAL const* restrict *restrict const RHSVars, + CCTK_REAL const dt, + int const totalsize) += { order1, order2, order3, order4, order5 }; +static int const max_order = sizeof orders / sizeof *orders; + + + +void MoL_ABAdd(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + CCTK_REAL const dt = CCTK_DELTA_TIME; + + /* Determine the order of accuracy */ + int order; + if (CCTK_EQUALS(AB_Type,"1")) { + order = 1 ; + } else if (CCTK_EQUALS(AB_Type,"2")) { + order = 2; + } else if (CCTK_EQUALS(AB_Type,"3")) { + order = 3; + } else if (CCTK_EQUALS(AB_Type,"4")) { + order = 4; + } else if (CCTK_EQUALS(AB_Type,"5")) { + order = 5; + } else { + abort(); + } + if (AB_initially_reduce_order) { + /* Reduce the order for the first time steps */ + /* int const iteration = cctk_iteration; */ + int const iteration = 1 + floor((cctk_time - cctk_initial_time) / dt + 0.5); + if (order > iteration) { + order = iteration; + CCTK_VInfo (CCTK_THORNSTRING, + "Reducing Adams-Bashforth order to %d", order); + } + } +/* printf ("MoL AB: iter=%d, order=%d\n", cctk_iteration, order); */ + assert (order >= 1 && order <= max_order); + + int totalsize = 1; + for (int arraydim = 0; arraydim < cctk_dim; arraydim++) { + totalsize *= cctk_ash[arraydim]; + } + + for (int var = 0; var < MoLNumEvolvedVariables; var++) { + CCTK_REAL* restrict const UpdateVar = + CCTK_VarDataPtrI(cctkGH, 0, EvolvedVariableIndex[var]); + CCTK_REAL const* restrict RHSVars[order]; + for (int tl = 0; tl < order; tl++) { + RHSVars[tl] = CCTK_VarDataPtrI(cctkGH, tl, RHSVariableIndex[var]); + } + + /* Add RHS */ + (orders[order-1]) (UpdateVar, RHSVars, dt, totalsize); + } /* var */ + + for (int var = 0; var < MoLNumEvolvedArrayVariables; var++) { + CCTK_REAL* restrict const UpdateVar = + CCTK_VarDataPtrI(cctkGH, 0, EvolvedArrayVariableIndex[var]); + CCTK_REAL const* restrict RHSVars[order]; + for (int tl = 0; tl < order; tl++) { + RHSVars[tl] = CCTK_VarDataPtrI(cctkGH, tl, RHSArrayVariableIndex[var]); + } + + int const groupindex = + CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]); + cGroupDynamicData arraydata; + int const ierr = CCTK_GroupDynamicData(cctkGH, groupindex, &arraydata); + if (ierr) { + CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "The driver does not return group information for group '%s'.", + CCTK_GroupName(groupindex)); + } + int arraytotalsize = 1; + for (int arraydim = 0; arraydim < arraydata.dim; arraydim++) { + arraytotalsize *= arraydata.ash[arraydim]; + } + + /* Add RHS */ + (orders[order-1]) (UpdateVar, RHSVars, dt, arraytotalsize); + } /* var */ + +} diff --git a/src/Euler.c b/src/Euler.c new file mode 100644 index 0000000..0033cd1 --- /dev/null +++ b/src/Euler.c @@ -0,0 +1,116 @@ + /*@@ + @file Euler.c + @date 2012-12-05 + @author Erik Schnetter + @desc + An explicit Euler time integrator. This method is unstable in most + cases, but is sometimes useful for debugging. + @enddesc + @version $Header$ + @@*/ + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include "ExternalVariables.h" +#include "Operators.h" + +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusBase_MoL_Euler_c); + +/******************************************************************** + ********************* Local Data Types *********************** + ********************************************************************/ + +/******************************************************************** + ********************* Local Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ***************** Scheduled Routine Prototypes ********************* + ********************************************************************/ + +void MoL_EulerAdd(CCTK_ARGUMENTS); + +/******************************************************************** + ********************* Other Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ********************* Local Data ***************************** + ********************************************************************/ + +/******************************************************************** + ********************* External Routines ********************** + ********************************************************************/ + + /*@@ + @routine MoL_EulerAdd + @date 2012-12-05 + @author Erik Schnetter + @desc + Performs first order Euler time integration. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +void MoL_EulerAdd(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + int totalsize = 1; + for (int arraydim = 0; arraydim < cctk_dim; arraydim++) + { + totalsize *= cctk_ash[arraydim]; + } + + switch (*MoL_Intermediate_Step) + { + case 1: + { + for (int var = 0; var < MoLNumEvolvedVariables; var++) + { + int const nsrcs = 2; + CCTK_INT const srcs[] = + {EvolvedVariableIndex[var], RHSVariableIndex[var]}; + CCTK_INT const tls[] = {1, 0}; + CCTK_REAL const facts[] = {1.0, CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedVariableIndex[var], 0.0, + srcs, tls, facts, nsrcs); + } + + for (int var = 0; var < MoLNumEvolvedArrayVariables; var++) + { + int const nsrcs = 2; + CCTK_INT const srcs[] = + {EvolvedArrayVariableIndex[var], RHSArrayVariableIndex[var]}; + CCTK_INT const tls[] = {1, 0}; + CCTK_REAL const facts[] = {1.0, CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedArrayVariableIndex[var], 0.0, + srcs, tls, facts, nsrcs); + } + + break; + } + default: + { + CCTK_WARN(CCTK_WARN_ABORT, + "Euler expects MoL_Intermediate_Step to be in [1,1]. " + "This should be caught at ParamCheck - bug Ian!"); + } + } +} + +/******************************************************************** + ********************* Local Routines ************************* + ********************************************************************/ diff --git a/src/GenericRK.c b/src/GenericRK.c index 4a7b854..6c88ca6 100644 --- a/src/GenericRK.c +++ b/src/GenericRK.c @@ -99,7 +99,7 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS) totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { - totalsize *= cctk_lsh[arraydim]; + totalsize *= cctk_ash[arraydim]; } if (scratchspace_firstindex == -99) diff --git a/src/ICN.c b/src/ICN.c index fb71783..af074e3 100644 --- a/src/ICN.c +++ b/src/ICN.c @@ -105,7 +105,7 @@ void MoL_ICNAdd(CCTK_ARGUMENTS) totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { - totalsize *= cctk_lsh[arraydim]; + totalsize *= cctk_ash[arraydim]; } #ifdef MOLDEBUG @@ -149,7 +149,7 @@ void MoL_ICNAdd(CCTK_ARGUMENTS) arraytotalsize = 1; for (arraydim = 0; arraydim < arraydata.dim; arraydim++) { - arraytotalsize *= arraydata.lsh[arraydim]; + arraytotalsize *= arraydata.ash[arraydim]; } /* CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, */ @@ -256,7 +256,7 @@ void MoL_ICNAverage(CCTK_ARGUMENTS) totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { - totalsize *= cctk_lsh[arraydim]; + totalsize *= cctk_ash[arraydim]; } #ifdef MOLDEBUG @@ -298,7 +298,7 @@ void MoL_ICNAverage(CCTK_ARGUMENTS) arraytotalsize = 1; for (arraydim = 0; arraydim < arraydata.dim; arraydim++) { - arraytotalsize *= arraydata.lsh[arraydim]; + arraytotalsize *= arraydata.ash[arraydim]; } #pragma omp parallel for diff --git a/src/IndexArrays.c b/src/IndexArrays.c index 485a4bd..65a4dac 100644 --- a/src/IndexArrays.c +++ b/src/IndexArrays.c @@ -12,6 +12,7 @@ #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" +#include "util_String.h" #include #include @@ -313,6 +314,10 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS) CCTK_WARN(0, "Generic_Type not recognized!"); } } + else if (CCTK_EQUALS(ODE_Method,"Euler")) + { + sprintf(infoline, "Euler"); + } else if (CCTK_EQUALS(ODE_Method,"RK2")) { sprintf(infoline, "Runge-Kutta 2"); @@ -352,6 +357,38 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS) "Averaging iterative Crank Nicholson with %i iterations", MoL_Intermediate_Steps); } + else if (CCTK_EQUALS(ODE_Method,"AB")) + { + if (CCTK_EQUALS(AB_Type,"1")) + { + Util_snprintf(infoline, sizeof infoline, + "Adams-Bashforth of order 1"); + } + else if (CCTK_EQUALS(AB_Type,"2")) + { + Util_snprintf(infoline, sizeof infoline, + "Adams-Bashforth of order 2"); + } + else if (CCTK_EQUALS(AB_Type,"3")) + { + Util_snprintf(infoline, sizeof infoline, + "Adams-Bashforth of order 3"); + } + else if (CCTK_EQUALS(AB_Type,"4")) + { + Util_snprintf(infoline, sizeof infoline, + "Adams-Bashforth of order 4"); + } + else if (CCTK_EQUALS(AB_Type,"5")) + { + Util_snprintf(infoline, sizeof infoline, + "Adams-Bashforth of order 5"); + } + else + { + CCTK_WARN(0, "AB_Type not recognized!"); + } + } else if (CCTK_EQUALS(ODE_Method,"RK2-MR-2:1")) { sprintf(infoline, "Multi-rate 2:1 Runge-Kutta 2"); diff --git a/src/InitialCopy.c b/src/InitialCopy.c index b421908..5cef13c 100644 --- a/src/InitialCopy.c +++ b/src/InitialCopy.c @@ -20,6 +20,7 @@ #include "cctk_Parameters.h" #include "ExternalVariables.h" +#include "Operators.h" static const char *rcsid = "$Header$"; @@ -104,11 +105,15 @@ void MoL_InitialCopy(CCTK_ARGUMENTS) totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { - totalsize *= cctk_lsh[arraydim]; + totalsize *= cctk_ash[arraydim]; } for (var = 0; var < MoLNumEvolvedVariables; var++) { + const int nsrc = 1; + const int srcs[1] = {EvolvedVariableIndex[var]}; + const int tls[1] = {1}; + const CCTK_REAL facts[1] = {1.0}; StorageOn = CCTK_QueryGroupStorageI(cctkGH, CCTK_GroupIndexFromVarI(EvolvedVariableIndex[var])); @@ -130,20 +135,8 @@ void MoL_InitialCopy(CCTK_ARGUMENTS) CCTK_WARN(0, "The grid function does not have storage assigned."); } - PreviousVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 1, - EvolvedVariableIndex[var]); - CurrentVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedVariableIndex[var]); - if (PreviousVar && CurrentVar) - { - memcpy(CurrentVar, PreviousVar, totalsize * sizeof(CCTK_REAL)); - } - else - { - CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,"Null pointer for variable %s", - CCTK_VarName(EvolvedVariableIndex[var])); - } - + MoL_LinearCombination(cctkGH, EvolvedVariableIndex[var], 0.0, + srcs, tls, facts, nsrc); } /* Set up the array sizes */ @@ -186,7 +179,7 @@ void MoL_InitialCopy(CCTK_ARGUMENTS) arraytotalsize = 1; for (arraydim = 0; arraydim < arraydata.dim; arraydim++) { - arraytotalsize *= arraydata.lsh[arraydim]; + arraytotalsize *= arraydata.ash[arraydim]; } ArrayScratchSizes[var] = arraytotalsize; @@ -526,12 +519,11 @@ void MoL_InitRHS(CCTK_ARGUMENTS) totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { - totalsize *= cctk_lsh[arraydim]; + totalsize *= cctk_ash[arraydim]; } for (var = 0; var < MoLNumEvolvedVariables; var++) { - StorageOn = CCTK_QueryGroupStorageI(cctkGH, CCTK_GroupIndexFromVarI(RHSVariableIndex[var])); @@ -552,22 +544,8 @@ void MoL_InitRHS(CCTK_ARGUMENTS) CCTK_WARN(0, "The grid function does not have storage assigned."); } - RHSVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, - RHSVariableIndex[var]); - if (RHSVar) - { -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - RHSVar[index] = 0; - } - } - else - { - CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,"Null pointer for variable %s", - CCTK_VarName(RHSVariableIndex[var])); - } - + MoL_LinearCombination(cctkGH, RHSVariableIndex[var], 0.0, + NULL, NULL, NULL, 0); } for (var = 0; var < MoLNumEvolvedArrayVariables; var++) @@ -588,7 +566,7 @@ void MoL_InitRHS(CCTK_ARGUMENTS) arraytotalsize = 1; for (arraydim = 0; arraydim < arraydata.dim; arraydim++) { - arraytotalsize *= arraydata.lsh[arraydim]; + arraytotalsize *= arraydata.ash[arraydim]; } if (arraytotalsize) @@ -690,7 +668,7 @@ void MoL_FillAllLevels(CCTK_ARGUMENTS) totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { - totalsize *= cctk_lsh[arraydim]; + totalsize *= cctk_ash[arraydim]; } for (var = 0; var < MoLNumEvolvedVariables; var++) diff --git a/src/Operators.c b/src/Operators.c new file mode 100644 index 0000000..e66312e --- /dev/null +++ b/src/Operators.c @@ -0,0 +1,312 @@ +#include "Operators.h" +#include +#include +#include + +/* These are MoL's low-level operators. If they are overloaded as + aliased functions, these aliased functions are called; otherwise, a + default implementation is used. */ + +/* The aliased functions should never be called from MoL's time + integrators, because they may not exist. Instead, the operators + defined here (with the MoL_ prefix) should be used. */ + +static +void +error_no_storage(int const var, int const tl) +{ + char *const fullname = CCTK_FullName(var); + CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Variable %s, timelevel %d has no storage", fullname, tl); + free(fullname); + return; +} + + + +/* Some common special cases to improve performance */ +static +void +op_real_set_0(CCTK_REAL *restrict const varptr, + ptrdiff_t const npoints) +{ +#pragma omp parallel for + for (ptrdiff_t i=0; i + +CCTK_INT +MoL_LinearCombination(cGH const *const cctkGH, + CCTK_INT const var, + CCTK_REAL const scale, + CCTK_INT const srcs[], + CCTK_INT const tls[], + CCTK_REAL const facts[], + CCTK_INT const nsrcs); + +#endif // #ifndef OPERATORS_H diff --git a/src/ParamCheck.c b/src/ParamCheck.c index b80aa1b..96494bf 100644 --- a/src/ParamCheck.c +++ b/src/ParamCheck.c @@ -116,6 +116,12 @@ void MoL_ParamCheck(CCTK_ARGUMENTS) } } + if ( (CCTK_Equals(ODE_Method, "Euler"))&&(!(MoL_Intermediate_Steps == 1)) ) + { + CCTK_PARAMWARN("When using the Euler evolver the " + "number of intermediate steps must be 1"); + } + if ( (CCTK_Equals(ODE_Method, "RK2"))&&(!(MoL_Intermediate_Steps == 2)) ) { CCTK_PARAMWARN("When using the efficient RK2 evolver the " @@ -160,6 +166,12 @@ void MoL_ParamCheck(CCTK_ARGUMENTS) "and the number of scratch levels at least 13."); } + if ( (CCTK_Equals(ODE_Method, "AB"))&&(!(MoL_Intermediate_Steps == 1)) ) + { + CCTK_PARAMWARN("When using the Adams-Bashforth evolver the " + "number of intermediate steps must be 1"); + } + 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 " diff --git a/src/RK2-MR-2_1.c b/src/RK2-MR-2_1.c index 8eaabd2..9179cdd 100644 --- a/src/RK2-MR-2_1.c +++ b/src/RK2-MR-2_1.c @@ -99,7 +99,7 @@ CCTK_WARN(0, "not implemented"); totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { - totalsize *= cctk_lsh[arraydim]; + totalsize *= cctk_ash[arraydim]; } if (scratchspace_firstindex == -99) diff --git a/src/RK2.c b/src/RK2.c index 58dd03f..ea4f12c 100644 --- a/src/RK2.c +++ b/src/RK2.c @@ -16,6 +16,7 @@ #include "cctk_Parameters.h" #include "ExternalVariables.h" +#include "Operators.h" static const char *rcsid = "$Header$"; @@ -68,15 +69,10 @@ void MoL_RK2Add(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; + CCTK_INT var; CCTK_INT totalsize; - CCTK_REAL const * restrict OldVar; - CCTK_REAL * restrict UpdateVar; - CCTK_REAL const * restrict RHSVar; /* FIXME */ @@ -103,7 +99,7 @@ void MoL_RK2Add(CCTK_ARGUMENTS) totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { - totalsize *= cctk_lsh[arraydim]; + totalsize *= cctk_ash[arraydim]; } switch (*MoL_Intermediate_Step) @@ -113,45 +109,24 @@ void MoL_RK2Add(CCTK_ARGUMENTS) { for (var = 0; var < MoLNumEvolvedVariables; var++) { - UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedVariableIndex[var]); - RHSVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 0, - RHSVariableIndex[var]); - -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - UpdateVar[index] += CCTK_DELTA_TIME * RHSVar[index]; - } + int const nsrcs = 1; + CCTK_INT const srcs[] = {RHSVariableIndex[var]}; + CCTK_INT const tls[] = {0}; + CCTK_REAL const facts[] = {CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedVariableIndex[var], 1.0, + srcs, tls, facts, nsrcs); } for (var = 0; var < MoLNumEvolvedArrayVariables; var++) { - UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedArrayVariableIndex[var]); - RHSVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 0, - RHSArrayVariableIndex[var]); - groupindex = CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]); - ierr = CCTK_GroupDynamicData(cctkGH, groupindex, - &arraydata); - if (ierr) - { - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "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]; - } - -#pragma omp parallel for - for (index = 0; index < arraytotalsize; index++) - { - UpdateVar[index] += CCTK_DELTA_TIME * RHSVar[index]; - } + int const nsrcs = 1; + CCTK_INT const srcs[] = {RHSArrayVariableIndex[var]}; + CCTK_INT const tls[] = {0}; + CCTK_REAL const facts[] = {CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedArrayVariableIndex[var], 1.0, + srcs, tls, facts, nsrcs); } /* FIXME */ @@ -160,18 +135,13 @@ void MoL_RK2Add(CCTK_ARGUMENTS) for (var = 0; var < MoLNumEvolvedComplexVariables; var++) { - UpdateComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedComplexVariableIndex[var]); - RHSComplexVar = (CCTK_COMPLEX const*)CCTK_VarDataPtrI(cctkGH, 0, - RHSComplexVariableIndex[var]); - -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - UpdateComplexVar[index] = CCTK_CmplxAdd(UpdateComplexVar[index], - CCTK_CmplxMul(Complex_Delta_Time, - RHSComplexVar[index])); - } + int const nsrcs = 1; + CCTK_INT const srcs[] = {RHSComplexVariableIndex[var]}; + CCTK_INT const tls[] = {0}; + CCTK_REAL const facts[] = {CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedComplexVariableIndex[var], 1.0, + srcs, tls, facts, nsrcs); } #endif @@ -182,52 +152,26 @@ void MoL_RK2Add(CCTK_ARGUMENTS) { for (var = 0; var < MoLNumEvolvedVariables; var++) { - OldVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 1, - EvolvedVariableIndex[var]); - UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedVariableIndex[var]); - RHSVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 0, - RHSVariableIndex[var]); - -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - UpdateVar[index] = 0.5 * (OldVar[index] + UpdateVar[index]) + - CCTK_DELTA_TIME * RHSVar[index]; - } + int const nsrcs = 2; + CCTK_INT const srcs[] = + {EvolvedVariableIndex[var], RHSVariableIndex[var]}; + CCTK_INT const tls[] = {1, 0}; + CCTK_REAL const facts[] = {0.5, CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedVariableIndex[var], 0.5, + srcs, tls, facts, nsrcs); } for (var = 0; var < MoLNumEvolvedArrayVariables; var++) { - OldVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 1, - EvolvedArrayVariableIndex[var]); - UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedArrayVariableIndex[var]); - RHSVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 0, - RHSArrayVariableIndex[var]); - - groupindex = CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]); - ierr = CCTK_GroupDynamicData(cctkGH, groupindex, - &arraydata); - if (ierr) - { - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "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]; - } - -#pragma omp parallel for - for (index = 0; index < arraytotalsize; index++) - { - UpdateVar[index] = 0.5 * (OldVar[index] + UpdateVar[index]) + - CCTK_DELTA_TIME * RHSVar[index]; - } + int const nsrcs = 2; + CCTK_INT const srcs[] = + {EvolvedArrayVariableIndex[var], RHSArrayVariableIndex[var]}; + CCTK_INT const tls[] = {1, 0}; + CCTK_REAL const facts[] = {0.5, CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedArrayVariableIndex[var], 0.5, + srcs, tls, facts, nsrcs); } /* FIXME */ @@ -236,23 +180,14 @@ void MoL_RK2Add(CCTK_ARGUMENTS) for (var = 0; var < MoLNumEvolvedComplexVariables; var++) { - OldComplexVar = (CCTK_COMPLEX const*)CCTK_VarDataPtrI(cctkGH, 1, - EvolvedComplexVariableIndex[var]); - UpdateComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedComplexVariableIndex[var]); - RHSComplexVar = (CCTK_COMPLEX const*)CCTK_VarDataPtrI(cctkGH, 0, - RHSComplexVariableIndex[var]); - -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - UpdateComplexVar[index] = - CCTK_CmplxAdd(CCTK_CmplxMul(Complex_Half, - (CCTK_CmplxAdd(OldComplexVar[index], - UpdateComplexVar[index]))), - CCTK_CmplxMul(Complex_Delta_Time, - RHSComplexVar[index])); - } + int const nsrcs = 2; + CCTK_INT const srcs[] = + {EvolvedComplexVariableIndex[var], RHSComplexVariableIndex[var]}; + CCTK_INT const tls[] = {1, 0}; + CCTK_REAL const facts[] = {0.5, CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedComplexVariableIndex[var], 0.5, + srcs, tls, facts, nsrcs); } #endif diff --git a/src/RK3.c b/src/RK3.c index f3b5e9c..38a55ab 100644 --- a/src/RK3.c +++ b/src/RK3.c @@ -15,6 +15,7 @@ #include "cctk_Parameters.h" #include "ExternalVariables.h" +#include "Operators.h" static const char *rcsid = "$Header$"; @@ -67,15 +68,10 @@ void MoL_RK3Add(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; + CCTK_INT var; CCTK_INT totalsize; - CCTK_REAL const * restrict OldVar; - CCTK_REAL * restrict UpdateVar; - CCTK_REAL const * restrict RHSVar; /* FIXME */ @@ -106,7 +102,7 @@ void MoL_RK3Add(CCTK_ARGUMENTS) totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { - totalsize *= cctk_lsh[arraydim]; + totalsize *= cctk_ash[arraydim]; } switch (*MoL_Intermediate_Step) @@ -116,45 +112,24 @@ void MoL_RK3Add(CCTK_ARGUMENTS) { for (var = 0; var < MoLNumEvolvedVariables; var++) { - UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedVariableIndex[var]); - RHSVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, - RHSVariableIndex[var]); - -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - UpdateVar[index] += CCTK_DELTA_TIME * RHSVar[index]; - } + int const nsrcs = 1; + CCTK_INT const srcs[] = {RHSVariableIndex[var]}; + CCTK_INT const tls[] = {0}; + CCTK_REAL const facts[] = {CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedVariableIndex[var], 1.0, + srcs, tls, facts, nsrcs); } 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__, CCTK_THORNSTRING, - "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]; - } - -#pragma omp parallel for - for (index = 0; index < arraytotalsize; index++) - { - UpdateVar[index] += CCTK_DELTA_TIME * RHSVar[index]; - } + int const nsrcs = 1; + CCTK_INT const srcs[] = {RHSArrayVariableIndex[var]}; + CCTK_INT const tls[] = {0}; + CCTK_REAL const facts[] = {CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedArrayVariableIndex[var], 1.0, + srcs, tls, facts, nsrcs); } /* FIXME */ @@ -163,18 +138,13 @@ void MoL_RK3Add(CCTK_ARGUMENTS) for (var = 0; var < MoLNumEvolvedComplexVariables; var++) { - UpdateComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedComplexVariableIndex[var]); - RHSComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0, - RHSComplexVariableIndex[var]); - -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - UpdateComplexVar[index] = CCTK_CmplxAdd(UpdateComplexVar[index], - CCTK_CmplxMul(Complex_Delta_Time, - RHSComplexVar[index])); - } + int const nsrcs = 1; + CCTK_INT const srcs[] = {RHSComplexVariableIndex[var]}; + CCTK_INT const tls[] = {0}; + CCTK_REAL const facts[] = {CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedComplexVariableIndex[var], 1.0, + srcs, tls, facts, nsrcs); } #endif @@ -186,52 +156,26 @@ void MoL_RK3Add(CCTK_ARGUMENTS) { for (var = 0; var < MoLNumEvolvedVariables; var++) { - OldVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 1, - EvolvedVariableIndex[var]); - UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedVariableIndex[var]); - RHSVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 0, - RHSVariableIndex[var]); - -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - UpdateVar[index] = 0.25 * (3*OldVar[index] + - UpdateVar[index]) + - CCTK_DELTA_TIME * RHSVar[index]; - } + int const nsrcs = 2; + CCTK_INT const srcs[] = + {EvolvedVariableIndex[var], RHSVariableIndex[var]}; + CCTK_INT const tls[] = {1, 0}; + CCTK_REAL const facts[] = {0.75, CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedVariableIndex[var], 0.25, + srcs, tls, facts, nsrcs); } for (var = 0; var < MoLNumEvolvedArrayVariables; var++) { - OldVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 1, - EvolvedArrayVariableIndex[var]); - UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedArrayVariableIndex[var]); - RHSVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 0, - RHSArrayVariableIndex[var]); - groupindex = CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]); - ierr = CCTK_GroupDynamicData(cctkGH, groupindex, - &arraydata); - if (ierr) - { - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "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]; - } - -#pragma omp parallel for - for (index = 0; index < arraytotalsize; index++) - { - UpdateVar[index] = 0.25*(3*OldVar[index] + UpdateVar[index]) + - CCTK_DELTA_TIME * RHSVar[index]; - } + int const nsrcs = 2; + CCTK_INT const srcs[] = + {EvolvedArrayVariableIndex[var], RHSArrayVariableIndex[var]}; + CCTK_INT const tls[] = {1, 0}; + CCTK_REAL const facts[] = {0.75, CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedArrayVariableIndex[var], 0.25, + srcs, tls, facts, nsrcs); } /* FIXME */ @@ -240,24 +184,14 @@ void MoL_RK3Add(CCTK_ARGUMENTS) for (var = 0; var < MoLNumEvolvedComplexVariables; var++) { - OldComplexVar = (CCTK_COMPLEX const*)CCTK_VarDataPtrI(cctkGH, 1, - EvolvedComplexVariableIndex[var]); - UpdateComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedComplexVariableIndex[var]); - RHSComplexVar = (CCTK_COMPLEX const*)CCTK_VarDataPtrI(cctkGH, 0, - RHSComplexVariableIndex[var]); - -#pragma omp parallel for - 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]))); - } + int const nsrcs = 2; + CCTK_INT const srcs[] = + {EvolvedComplexVariableIndex[var], RHSComplexVariableIndex[var]}; + CCTK_INT const tls[] = {1, 0}; + CCTK_REAL const facts[] = {0.75, CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedComplexVariableIndex[var], 0.25, + srcs, tls, facts, nsrcs); } #endif @@ -266,56 +200,28 @@ void MoL_RK3Add(CCTK_ARGUMENTS) } case 1: { - const CCTK_REAL one_third = 1.0 / 3.0; - for (var = 0; var < MoLNumEvolvedVariables; var++) { - OldVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 1, - EvolvedVariableIndex[var]); - UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedVariableIndex[var]); - RHSVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 0, - RHSVariableIndex[var]); - -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - UpdateVar[index] = (OldVar[index] + 2*UpdateVar[index]) * one_third - + CCTK_DELTA_TIME * RHSVar[index]; - } + int const nsrcs = 2; + CCTK_INT const srcs[] = + {EvolvedVariableIndex[var], RHSVariableIndex[var]}; + CCTK_INT const tls[] = {1, 0}; + CCTK_REAL const facts[] = {1.0/3.0, CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedVariableIndex[var], 2.0/3.0, + srcs, tls, facts, nsrcs); } for (var = 0; var < MoLNumEvolvedArrayVariables; var++) { - OldVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 1, - EvolvedArrayVariableIndex[var]); - UpdateVar = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedArrayVariableIndex[var]); - RHSVar = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, 0, - RHSArrayVariableIndex[var]); - - groupindex = CCTK_GroupIndexFromVarI(EvolvedArrayVariableIndex[var]); - ierr = CCTK_GroupDynamicData(cctkGH, groupindex, - &arraydata); - if (ierr) - { - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "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]; - } - -#pragma omp parallel for - for (index = 0; index < arraytotalsize; index++) - { - UpdateVar[index] = (OldVar[index] + 2*UpdateVar[index]) * one_third - + CCTK_DELTA_TIME * RHSVar[index]; - } + int const nsrcs = 2; + CCTK_INT const srcs[] = + {EvolvedArrayVariableIndex[var], RHSArrayVariableIndex[var]}; + CCTK_INT const tls[] = {1, 0}; + CCTK_REAL const facts[] = {1.0/3.0, CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedArrayVariableIndex[var], 2.0/3.0, + srcs, tls, facts, nsrcs); } /* FIXME */ @@ -324,24 +230,14 @@ void MoL_RK3Add(CCTK_ARGUMENTS) for (var = 0; var < MoLNumEvolvedComplexVariables; var++) { - OldComplexVar = (CCTK_COMPLEX const*)CCTK_VarDataPtrI(cctkGH, 1, - EvolvedComplexVariableIndex[var]); - UpdateComplexVar = (CCTK_COMPLEX*)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedComplexVariableIndex[var]); - RHSComplexVar = (CCTK_COMPLEX const*)CCTK_VarDataPtrI(cctkGH, 0, - RHSComplexVariableIndex[var]); - -#pragma omp parallel for - 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])); - } + int const nsrcs = 2; + CCTK_INT const srcs[] = + {EvolvedComplexVariableIndex[var], RHSComplexVariableIndex[var]}; + CCTK_INT const tls[] = {1, 0}; + CCTK_REAL const facts[] = {1.0/3.0, CCTK_DELTA_TIME}; + MoL_LinearCombination(cctkGH, + EvolvedComplexVariableIndex[var], 2.0/3.0, + srcs, tls, facts, nsrcs); } #endif diff --git a/src/RK4-MR-2_1.c b/src/RK4-MR-2_1.c index 815a41e..0a48cb0 100644 --- a/src/RK4-MR-2_1.c +++ b/src/RK4-MR-2_1.c @@ -99,7 +99,7 @@ CCTK_WARN(0, "not implemented"); totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { - totalsize *= cctk_lsh[arraydim]; + totalsize *= cctk_ash[arraydim]; } if (scratchspace_firstindex == -99) diff --git a/src/RK4-RK2.c b/src/RK4-RK2.c index 6e5b85e..472418e 100644 --- a/src/RK4-RK2.c +++ b/src/RK4-RK2.c @@ -54,7 +54,7 @@ void MoL_RK4_RK2_Add(CCTK_ARGUMENTS) int const step = MoL_Intermediate_Steps - *MoL_Intermediate_Step; int totalsize = 1; - for (int d=0; dcctk_timefac; diff --git a/src/RK4.c b/src/RK4.c index 1d790ab..4e72389 100644 --- a/src/RK4.c +++ b/src/RK4.c @@ -13,8 +13,8 @@ #include "cctk_Arguments.h" #include "cctk_Parameters.h" -#include #include "ExternalVariables.h" +#include "Operators.h" static const char *rcsid = "$Header$"; @@ -99,14 +99,14 @@ CCTK_WARN(0, "not implemented"); /* Keep a running total of alpha as we perform the substeps, so that we know the "real" alpha (including round-off errors) when we calculate the final result. */ - static CCTK_REAL sum_alpha; - /* the last MoL intermediate step that was summed */ - static CCTK_INT last_sum_step; - + CCTK_REAL const time_rhs = 1.0; + CCTK_REAL const old_time = 0.0; + static CCTK_REAL time, scratch_time; + totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { - totalsize *= cctk_lsh[arraydim]; + totalsize *= cctk_ash[arraydim]; } if (scratchspace_firstindex == -99) @@ -131,21 +131,14 @@ CCTK_WARN(0, "not implemented"); case 3: alpha = 1.0; beta = 1.0 / 6.0; + break; } - /* Initialise alpha before the first intermediate step */ if (MoL_Intermediate_Steps == (*MoL_Intermediate_Step)) { - sum_alpha = 0.0; - last_sum_step = -1; - } - /* Add alpha once per intermediate step */ - if (last_sum_step != (*MoL_Intermediate_Step)) - { - sum_alpha += alpha; - last_sum_step = (*MoL_Intermediate_Step); + time = 0.0; } - + /* FIXME */ #ifdef MOLDOESCOMPLEX @@ -159,71 +152,53 @@ CCTK_WARN(0, "not implemented"); 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 generic RK. Variable %d (%s). RHS %d (%s). beta %g.\n", - EvolvedVariableIndex[var], - CCTK_VarName(EvolvedVariableIndex[var]), - RHSVariableIndex[var], - CCTK_VarName(RHSVariableIndex[var]), - beta); -#endif - -#pragma omp parallel for - for (index = 0; index < totalsize; index++) { - UpdateVar[index] = OldVar[index] + - (*Original_Delta_Time) / cctkGH->cctk_timefac * beta * 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 + int const nsrcs = 2; + CCTK_INT const srcs[] = + {EvolvedVariableIndex[var], RHSVariableIndex[var]}; + CCTK_INT const tls[] = {1, 0}; + CCTK_REAL const facts[] = + {1.0, (*Original_Delta_Time) / cctkGH->cctk_timefac * beta}; + MoL_LinearCombination(cctkGH, + EvolvedVariableIndex[var], 0.0, + srcs, tls, facts, nsrcs); + time = facts[0] * old_time + facts[1] * time_rhs; } - ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, - scratchspace_firstindex - + var ); - + /* scratch storage */ if ((*MoL_Intermediate_Step) == MoL_Intermediate_Steps) { -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - ScratchVar[index] = 0; - } + MoL_LinearCombination(cctkGH, + scratchspace_firstindex + var, 0.0, + NULL, NULL, NULL, 0); + scratch_time = 0.0; } if ((*MoL_Intermediate_Step)>1) { -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - ScratchVar[index] += alpha * UpdateVar[index]; - } + int const nsrcs = 1; + CCTK_INT const srcs[] = {EvolvedVariableIndex[var]}; + CCTK_INT const tls[] = {0}; + CCTK_REAL const facts[] = {alpha}; + MoL_LinearCombination(cctkGH, + scratchspace_firstindex + var, 1.0, + srcs, tls, facts, nsrcs); + scratch_time += facts[0] * time; } else { -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - UpdateVar[index] += ScratchVar[index] - (sum_alpha - 1.0) * OldVar[index]; - } + int const nsrcs = 2; + CCTK_INT const srcs[] = + {scratchspace_firstindex + var, EvolvedVariableIndex[var]}; + CCTK_INT const tls[] = {0, 1}; + CCTK_REAL const facts[] = {1.0, -4.0/3.0}; + MoL_LinearCombination(cctkGH, + EvolvedVariableIndex[var], 1.0, + srcs, tls, facts, nsrcs); + time += facts[0] * scratch_time + facts[1] * old_time; } - } - /* Complex GFs */ /* FIXME */ @@ -259,7 +234,7 @@ CCTK_WARN(0, "not done"); arraytotalsize = 1; for (arraydim = 0; arraydim < arraydata.dim; arraydim++) { - arraytotalsize *= arraydata.lsh[arraydim]; + arraytotalsize *= arraydata.ash[arraydim]; } ScratchVar = &ArrayScratchSpace[arrayscratchlocation]; @@ -293,7 +268,7 @@ CCTK_WARN(0, "not done"); #pragma omp parallel for for (index = 0; index < arraytotalsize; index++) { - UpdateVar[index] += ScratchVar[index] - (sum_alpha - 1.0) * OldVar[index]; + UpdateVar[index] += ScratchVar[index] - 4.0 / 3.0 * OldVar[index]; } } arrayscratchlocation += arraytotalsize; diff --git a/src/RK45.c b/src/RK45.c index ea547fd..2e520da 100644 --- a/src/RK45.c +++ b/src/RK45.c @@ -159,7 +159,7 @@ void MoL_RK45Add(CCTK_ARGUMENTS) int totalsize = 1; for (int arraydim = 0; arraydim < cctk_dim; arraydim++) { - totalsize *= cctk_lsh[arraydim]; + totalsize *= cctk_ash[arraydim]; } /* Real GFs */ diff --git a/src/RK65.c b/src/RK65.c index f333741..e5dc629 100644 --- a/src/RK65.c +++ b/src/RK65.c @@ -115,7 +115,7 @@ void MoL_RK65Add(CCTK_ARGUMENTS) totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { - totalsize *= cctk_lsh[arraydim]; + totalsize *= cctk_ash[arraydim]; } /* Real GFs */ diff --git a/src/RK87.c b/src/RK87.c index dee29da..68c7ff2 100644 --- a/src/RK87.c +++ b/src/RK87.c @@ -14,6 +14,7 @@ #include "cctk_Parameters.h" #include "ExternalVariables.h" +#include "Operators.h" static const char *rcsid = "$Header$"; @@ -70,17 +71,9 @@ void MoL_RK87Add(CCTK_ARGUMENTS) CCTK_INT arraydim; static CCTK_INT scratchspace_firstindex = -99; - CCTK_INT index, var, scratchstep; + CCTK_INT var, scratchstep; CCTK_INT totalsize; - CCTK_REAL * restrict UpdateVar; - CCTK_REAL const * restrict RHSVar; - CCTK_REAL * restrict ScratchVar; - CCTK_REAL * restrict ErrorVar; - CCTK_REAL const * restrict OldVar; - - CCTK_REAL beta, gamma, gamma_error; - static const CCTK_REAL beta_array[12][12] = { { 1.0/18.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, { 1.0/48.0, 1.0/16.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, @@ -131,7 +124,7 @@ void MoL_RK87Add(CCTK_ARGUMENTS) totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { - totalsize *= cctk_lsh[arraydim]; + totalsize *= cctk_ash[arraydim]; } if (scratchspace_firstindex == -99) @@ -145,122 +138,101 @@ void MoL_RK87Add(CCTK_ARGUMENTS) for (var = 0; var < MoLNumEvolvedVariables; var++) { - const CCTK_REAL tmp = (*Original_Delta_Time) / cctkGH->cctk_timefac; - - UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedVariableIndex[var]); - RHSVar = (CCTK_REAL const *)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))); -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - ScratchVar[index] = tmp * RHSVar[index]; - } - + int const step = MoL_Intermediate_Steps - (*MoL_Intermediate_Step); + int const scratchvarindex = + scratchspace_firstindex + MoL_Num_Evolved_Vars * step + var; + + int const nsrcs = 1; + CCTK_INT const srcs[] = {RHSVariableIndex[var]}; + CCTK_INT const tls[] = {0}; + CCTK_REAL const facts[] = {(*Original_Delta_Time) / cctkGH->cctk_timefac}; + MoL_LinearCombination(cctkGH, + scratchvarindex, 1.0, + srcs, tls, facts, nsrcs); } - for (var = 0; var < MoLNumEvolvedVariables; var++) { - OldVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 1, - EvolvedVariableIndex[var]); - UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, - EvolvedVariableIndex[var]); - RHSVar = (CCTK_REAL const *)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); + int const error_firstindex = CCTK_FirstVarIndex("MOL::ERRORESTIMATE"); + int const step = MoL_Intermediate_Steps - (*MoL_Intermediate_Step); if (*MoL_Intermediate_Step - 1) - { - -#pragma omp parallel for - for (index = 0; index < totalsize; index++) + { + int const nsrcs = step + 2; + CCTK_INT srcs[nsrcs]; + CCTK_INT tls[nsrcs]; + CCTK_REAL facts[nsrcs]; + srcs[0] = EvolvedVariableIndex[var]; + tls[0] = 1; + facts[0] = 1.0; + + for (scratchstep = 0; scratchstep < step + 1; scratchstep++) { - UpdateVar[index] = OldVar[index]; + int const scratchvarindex = + scratchspace_firstindex + MoL_Num_Evolved_Vars * scratchstep + var; + srcs [scratchstep+1] = scratchvarindex; + tls [scratchstep+1] = 0; + facts[scratchstep+1] = beta_array[step][scratchstep]; } - for (scratchstep = 0; - scratchstep < MoL_Intermediate_Steps - (*MoL_Intermediate_Step) + 1; - scratchstep++) + MoL_LinearCombination(cctkGH, + EvolvedVariableIndex[var], 0.0, + srcs, tls, facts, nsrcs); + } + else + { { + int const nsrcs = 14; + CCTK_INT srcs[nsrcs]; + CCTK_INT tls[nsrcs]; + CCTK_REAL facts[nsrcs]; + srcs[0] = EvolvedVariableIndex[var]; + tls[0] = 1; + facts[0] = 1.0; - 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 (scratchstep = 0; scratchstep < 13; scratchstep++) { -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - UpdateVar[index] += beta * ScratchVar[index]; - } + int const scratchvarindex = + scratchspace_firstindex + MoL_Num_Evolved_Vars * scratchstep + var; + srcs [scratchstep+1] = scratchvarindex; + tls [scratchstep+1] = 0; + facts[scratchstep+1] = gamma_array[scratchstep]; } - - } - - } - else - { -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - UpdateVar[index] = OldVar[index]; - ErrorVar[index] = 0; + MoL_LinearCombination(cctkGH, + EvolvedVariableIndex[var], 0.0, + srcs, tls, facts, nsrcs); } - - for (scratchstep = 0; scratchstep < 13; scratchstep++) + { + int const nsrcs = 13; + CCTK_INT srcs[nsrcs]; + CCTK_INT tls[nsrcs]; + CCTK_REAL facts[nsrcs]; - 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 (scratchstep = 0; scratchstep < 13; scratchstep++) { -#pragma omp parallel for - for (index = 0; index < totalsize; index++) - { - UpdateVar[index] += gamma * ScratchVar[index]; - ErrorVar[index] += gamma_error * ScratchVar[index]; - } + int const scratchvarindex = + scratchspace_firstindex + MoL_Num_Evolved_Vars * scratchstep + var; + srcs [scratchstep] = scratchvarindex; + tls [scratchstep] = 0; + facts[scratchstep] = + gamma_array[scratchstep] - gammastar_array[scratchstep]; } + MoL_LinearCombination(cctkGH, + error_firstindex + var, 0.0, + srcs, tls, facts, nsrcs); } - + } - + } /* Real arrays */ - for (var = 0; var < MoLNumEvolvedArrayVariables; var++) - { - - CCTK_WARN(0, "Peter has been too lazy to write the RK87 routine " - "out for array variables. Better send him an email..."); - - } + CCTK_WARN(0, "Peter has been too lazy to write the RK87 routine " + "out for array variables. Better send him an email..."); return; } diff --git a/src/SandR.c b/src/SandR.c index e644c53..577a8fe 100644 --- a/src/SandR.c +++ b/src/SandR.c @@ -91,7 +91,7 @@ void MoL_RestoreSandR(CCTK_ARGUMENTS) totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { - totalsize *= cctk_lsh[arraydim]; + totalsize *= cctk_ash[arraydim]; } for (var = 0; var < MoLNumSandRVariables; var++) diff --git a/src/SetTime.c b/src/SetTime.c index edef89f..b60f3bb 100644 --- a/src/SetTime.c +++ b/src/SetTime.c @@ -234,6 +234,13 @@ void MoL_ResetTime(CCTK_ARGUMENTS) cctkGH->cctk_time = previous_times[MoL_Intermediate_Steps - *MoL_Intermediate_Step]; } + else if (CCTK_EQUALS(ODE_Method,"Euler")) + { + if (*MoL_Intermediate_Step == 1) + { + cctkGH->cctk_time = (*Original_Time); + } + } else if (CCTK_EQUALS(ODE_Method,"RK2")) { if (*MoL_Intermediate_Step == 1) @@ -433,6 +440,13 @@ void MoL_ResetDeltaTime(CCTK_ARGUMENTS) (*MoL_Intermediate_Step)] * (*Original_Delta_Time); } + else if (CCTK_EQUALS(ODE_Method,"Euler")) + { + if (*MoL_Intermediate_Step == 1) + { + cctkGH->cctk_delta_time = (*Original_Delta_Time); + } + } else if (CCTK_EQUALS(ODE_Method,"RK2")) { if (*MoL_Intermediate_Step == 1) diff --git a/src/make.code.defn b/src/make.code.defn index 80e1b56..f11c36b 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -2,14 +2,20 @@ # $Header$ # Source files in this directory -SRCS = ChangeType.c \ +SRCS = AB.c \ + ChangeType.c \ Counter.c \ + Euler.c \ GenericRK.c \ ICN.c \ IndexArrays.c \ InitialCopy.c \ + Operators.c \ ParamCheck.c \ RK2.c \ + RK2-MR-2_1.c \ + RK4-MR-2_1.c \ + RK4-RK2.c \ RK3.c \ RK4.c \ RK45.c \ @@ -21,11 +27,8 @@ SRCS = ChangeType.c \ SandR.c \ SetTime.c \ Startup.c \ - StepSize.c \ - RK2-MR-2_1.c \ - RK4-MR-2_1.c \ - RK4-RK2.c - + StepSize.c + # Subdirectories containing source files SUBDIRS = -- cgit v1.2.3