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 --- src/AB.c | 234 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 234 insertions(+) create mode 100644 src/AB.c (limited to 'src/AB.c') 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 */ + +} -- cgit v1.2.3