diff options
| -rw-r--r-- | interface.ccl | 10 | ||||
| -rw-r--r-- | param.ccl | 38 | ||||
| -rw-r--r-- | schedule.ccl | 89 | ||||
| -rw-r--r-- | src/ParamCheck.c | 13 | ||||
| -rw-r--r-- | src/StepSize.c | 347 | ||||
| -rw-r--r-- | src/make.code.defn | 3 |
6 files changed, 480 insertions, 20 deletions
diff --git a/interface.ccl b/interface.ccl index f9f45f5..40e1c79 100644 --- a/interface.ccl +++ b/interface.ccl @@ -194,6 +194,7 @@ CCTK_REAL RKBetaCoefficients TYPE=ARRAY DIM=1 SIZE=MoL_Intermediate_Steps DISTRI CCTK_INT MoL_Counters TYPE = SCALAR { MoL_Intermediate_Step + MoL_Stepsize_Bad } "The counter for the time integration method" CCTK_REAL MoL_Original_Time TYPE = SCALAR @@ -221,6 +222,13 @@ CCTK_REAL ArraySandRScratchSpace TYPE = ARRAY DIM = 1 SIZE = MoL_Max_Evolved_Arr #CCTK_COMPLEX ComplexArraySandRScratchSpace TYPE = ARRAY DIM = 1 SIZE = MoL_Max_Evolved_ComplexArray_Size Timelevels = 1 -#Error vector for RK45 +#Error vector and scalars for RK45 CCTK_REAL ErrorEstimate[MoL_Num_Evolved_Vars] TYPE=GF Timelevels = 1 tags='Prolongation="None"' + +CCTK_REAL ErrorScalars TYPE=SCALAR +{ + Error + Count + EstimatedDt +} "Global error estimate" @@ -132,6 +132,44 @@ BOOLEAN disable_prolongation "If Mesh refinement is enabled should we use buffer { } "yes" + + +BOOLEAN adaptive_stepsize "Choose the time step size adaptively" +{ +} "no" + +REAL maximum_absolute_error "Maximum allowed absolute error for adaptive stepsize control" +{ + (0.0:*) :: "" +} 1.0e-6 + +REAL maximum_relative_error "Maximum allowed relative error for adaptive stepsize control" +{ + (0.0:*) :: "" +} 1.0e-6 + +REAL RHS_error_weight "Weight of the RHS in the relative error calculation" +{ + 0.0:* :: "should be between zero and one" +} 1.0 + +REAL safety_factor "Safety factor for stepsize control" +{ + (0.0:*) :: "should be less than one" +} 0.9 + +REAL maximum_decrease "Maximum stepsize decrease factor" +{ + (0.0:*) :: "should be larger than one" +} 10.0 + +REAL maximum_increase "Maximum stepsize increase factor" +{ + (0.0:*) :: "should be larger than one" +} 5.0 + + + KEYWORD verbose "How verbose should MoL be?" { "none" :: "No output at all (not implemented)" diff --git a/schedule.ccl b/schedule.ccl index 43dc1f9..4e0ae19 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -130,6 +130,22 @@ if (initial_data_is_crap) } "A bad routine. Fills all previous timelevels with data copied from the current." } +########################################## +### Initialise the step size control ### +########################################## + +schedule MoL_StartLoop AT Evol BEFORE MoL_Evolution +{ + LANG: C + OPTIONS: LEVEL +} "Initialise the step size control" + +schedule MoL_StartLoop AT Initial +{ + LANG: C + OPTIONS: LEVEL +} "Initialise the step size control" + ###################################################### ### The evolution step. This is almost a self ### ### contained EVOL step with PRE and POST steps ### @@ -149,7 +165,7 @@ if (MoL_Num_Scratch_Levels) { if (MoL_Max_Evolved_Array_Size) { - schedule GROUP MoL_Evolution AT Evol + schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad { STORAGE: ScratchSpace STORAGE: SandRScratchSpace @@ -163,7 +179,7 @@ if (MoL_Num_Scratch_Levels) } else { - schedule GROUP MoL_Evolution AT Evol + schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad { STORAGE: ScratchSpace STORAGE: SandRScratchSpace @@ -178,7 +194,7 @@ if (MoL_Num_Scratch_Levels) { if (MoL_Max_Evolved_Array_Size) { - schedule GROUP MoL_Evolution AT Evol + schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad { STORAGE: ScratchSpace STORAGE: ArrayScratchSpace @@ -191,7 +207,7 @@ if (MoL_Num_Scratch_Levels) } else { - schedule GROUP MoL_Evolution AT Evol + schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad { STORAGE: ScratchSpace # STORAGE: ComplexScratchSpace @@ -208,7 +224,7 @@ if (MoL_Num_Scratch_Levels) { if (MoL_Max_Evolved_Array_Size) { - schedule GROUP MoL_Evolution AT Evol + schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad { STORAGE: SandRScratchSpace STORAGE: ArrayScratchSpace @@ -221,7 +237,7 @@ if (MoL_Num_Scratch_Levels) } else { - schedule GROUP MoL_Evolution AT Evol + schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad { STORAGE: SandRScratchSpace # STORAGE: ComplexScratchSpace @@ -235,7 +251,7 @@ if (MoL_Num_Scratch_Levels) { if (MoL_Max_Evolved_Array_Size) { - schedule GROUP MoL_Evolution AT Evol + schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad { STORAGE: ArrayScratchSpace STORAGE: ArraySandRScratchSpace @@ -247,7 +263,7 @@ if (MoL_Num_Scratch_Levels) } else { - schedule GROUP MoL_Evolution AT Evol + schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad { # STORAGE: ComplexScratchSpace # STORAGE: ComplexSandRScratchSpace @@ -266,7 +282,7 @@ else { if (MoL_Max_Evolved_Array_Size) { - schedule GROUP MoL_Evolution AT Evol + schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad { STORAGE: SandRScratchSpace STORAGE: ArraySandRScratchSpace @@ -278,7 +294,7 @@ else } else { - schedule GROUP MoL_Evolution AT Evol + schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad { STORAGE: SandRScratchSpace # STORAGE: ComplexScratchSpace @@ -292,7 +308,7 @@ else { if (MoL_Max_Evolved_Array_Size) { - schedule GROUP MoL_Evolution AT Evol + schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad { STORAGE: ArraySandRScratchSpace # STORAGE: ComplexScratchSpace @@ -303,7 +319,7 @@ else } else { - schedule GROUP MoL_Evolution AT Evol + schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad { # STORAGE: ComplexScratchSpace # STORAGE: ComplexSandRScratchSpace @@ -319,7 +335,7 @@ else { if (MoL_Max_Evolved_Array_Size) { - schedule GROUP MoL_Evolution AT Evol + schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad { STORAGE: SandRScratchSpace STORAGE: ArraySandRScratchSpace @@ -331,7 +347,7 @@ else } else { - schedule GROUP MoL_Evolution AT Evol + schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad { STORAGE: SandRScratchSpace # STORAGE: ComplexScratchSpace @@ -345,7 +361,7 @@ else { if (MoL_Max_Evolved_Array_Size) { - schedule GROUP MoL_Evolution AT Evol + schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad { STORAGE: ArraySandRScratchSpace # STORAGE: ComplexScratchSpace @@ -356,7 +372,7 @@ else } else { - schedule GROUP MoL_Evolution AT Evol + schedule GROUP MoL_Evolution AT Evol WHILE MoL::MoL_Stepsize_Bad { # STORAGE: ComplexScratchSpace # STORAGE: ComplexSandRScratchSpace @@ -515,7 +531,7 @@ else if (CCTK_Equals(ODE_Method,"RK3")) } else if (CCTK_Equals(ODE_Method,"RK45")) { - STORAGE: ErrorEstimate + STORAGE: ErrorEstimate ErrorScalars schedule MoL_RK45Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep { @@ -595,6 +611,45 @@ schedule MoL_RestoreSandR IN MoL_Evolution AFTER MoL_PostStep LANG: C } "Restoring the Save and Restore variables to the original state" +################################################### +### Loop until the step size was small enough ### +################################################### + +if (adaptive_stepsize) +{ + # Adaptive step size control + schedule MoL_InitAdaptiveError IN MoL_Evolution AFTER MoL_PostStep + { + LANG: C + OPTIONS: LEVEL + } "Control the step size: initialize error check variables" + + schedule MoL_FindAdaptiveError IN MoL_Evolution AFTER MoL_InitAdaptiveError + { + LANG: C + } "Control the step size: compute error check variables" + + schedule MoL_ReduceAdaptiveError IN MoL_Evolution AFTER MoL_FindAdaptiveError + { + LANG: C + OPTIONS: LEVEL + } "Control the step size: reduce error check variables" + + schedule MoL_SetEstimatedDt AT POSTSTEP + { + LANG: C + OPTIONS: LEVEL + } "Control the step size: set the new timestep" +} +else +{ + schedule MoL_FinishLoop IN MoL_Evolution AFTER MoL_PostStep + { + LANG: C + OPTIONS: LEVEL + } "Control the step size" +} + ################################################################ ### At the end (but before driver terminate to avoid those ### ### irritating segfaults) free the index arrays. ### diff --git a/src/ParamCheck.c b/src/ParamCheck.c index bd91e5e..3bd1c21 100644 --- a/src/ParamCheck.c +++ b/src/ParamCheck.c @@ -128,7 +128,6 @@ int MoL_ParamCheck(CCTK_ARGUMENTS) "number of intermediate steps must be 3"); } - if ( (CCTK_Equals(ODE_Method, "RK45")) && ( !((MoL_Intermediate_Steps == 6)||(MoL_Num_Scratch_Levels > 5)) ) ) { @@ -137,6 +136,18 @@ int MoL_ParamCheck(CCTK_ARGUMENTS) "and the number of scratch levels at least 6."); } + if (adaptive_stepsize) + { + if (CCTK_Equals(ODE_Method, "RK45")) + { + /* everything is fine, do nothing */ + } + else + { + CCTK_PARAMWARN("Adaptive time step sizes are only possible with the RK45 solver"); + } + } + return 0; } diff --git a/src/StepSize.c b/src/StepSize.c new file mode 100644 index 0000000..61266fc --- /dev/null +++ b/src/StepSize.c @@ -0,0 +1,347 @@ + /*@@ + @file StepSize.c + @date Tue Sep 07 2004 + @author Erik Schnetter + @desc + Control the time step size. + @enddesc + @version $Header$ + @@*/ + +#include <assert.h> +#include <math.h> +#include <stdlib.h> + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include "ExternalVariables.h" + +static const char *rcsid = "$Header$"; + +CCTK_FILEVERSION(CactusBase_MoL_StepSize_c); + +/******************************************************************** + ********************* Local Data Types *********************** + ********************************************************************/ + +/******************************************************************** + ********************* Local Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ***************** Scheduled Routine Prototypes ********************* + ********************************************************************/ + +void MoL_StartLoop(CCTK_ARGUMENTS); + +void MoL_InitAdaptiveError(CCTK_ARGUMENTS); +void MoL_FindAdaptiveError(CCTK_ARGUMENTS); +void MoL_ReduceAdaptiveError(CCTK_ARGUMENTS); + +void MoL_FinishLoop(CCTK_ARGUMENTS); + +/******************************************************************** + ********************* Other Routine Prototypes ********************* + ********************************************************************/ + +/******************************************************************** + ********************* Local Data ***************************** + ********************************************************************/ + +/******************************************************************** + ********************* External Routines ********************** + ********************************************************************/ + + /*@@ + @routine MoL_StartLoop + @date Tue Sep 07 2004 + @author Erik Schnetter + @desc + Start the step size control loop, so that at least one iteration is done. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +void +MoL_StartLoop(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + *MoL_Stepsize_Bad = 1; + + if (adaptive_stepsize) + { + *EstimatedDt = cctkGH->cctk_delta_time; + } + +} + + /*@@ + @routine MoL_InitAdaptiveError + @date Tue Sep 07 2004 + @author Erik Schnetter + @desc + Initialize error counters for adaptive stepsize control + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +static inline CCTK_REAL +square (CCTK_REAL const x) +{ + return x * x; +} + +void MoL_InitAdaptiveError(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + /* Initialise global error */ + *Error = 0; + *Count = 0; +} + + /*@@ + @routine MoL_FindAdaptiveError + @date Thu Jan 27 10:22:26 2005 + @author Erik Schnetter + @desc + Compute the error local to this component/patch/... + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +void MoL_FindAdaptiveError(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + CCTK_REAL const * restrict UpdateVar; + CCTK_REAL const * restrict RHSVar; + CCTK_REAL const * restrict ErrorVar; + + int imin[3], imax[3]; + + int var; + int i, j, k; + int d; + + int ierr; + + assert (cctk_dim <= 3); + for (d = 0; d < cctk_dim; d++) + { + imin[d] = cctk_bbox[2*d] ? 0 : cctk_nghostzones[d]; + imax[d] = cctk_lsh[d] - (cctk_bbox[2*d+1] ? 0 : cctk_nghostzones[d]); + } + + /* Calculate absolute error */ + for (var = 0; var < MoLNumEvolvedVariables; var++) + { + + UpdateVar = CCTK_VarDataPtrI(cctkGH, 0, EvolvedVariableIndex[var]); + RHSVar = CCTK_VarDataPtrI(cctkGH, 0, RHSVariableIndex[var]); + ErrorVar + = CCTK_VarDataPtrI(cctkGH, 0, + CCTK_FirstVarIndex("MOL::ERRORESTIMATE") + var); + + assert (cctk_dim == 3); + for (k = imin[2]; k < imax[2]; k++) + { + for (j = imin[1]; j < imax[1]; j++) + { + for (i = imin[0]; i < imax[0]; i++) + { + int const index = CCTK_GFINDEX3D(cctkGH, i, j, k); + CCTK_REAL const scale + = (square(maximum_absolute_error) + + square(maximum_relative_error * UpdateVar[index]) + + square(maximum_relative_error * RHS_error_weight + * (*Original_Delta_Time) * RHSVar[index])); + *Error += square(ErrorVar[index]) / scale; + } + } + } + + } /* for var */ + + *Count + += (MoLNumEvolvedVariables + * (imax[0] - imin[0]) * (imax[1] - imin[1]) * (imax[0] - imin[0])); +} + + /*@@ + @routine MoL_ReduceAdaptiveError + @date Thu Jan 27 10:23:14 2005 + @author Erik Schnetter + @desc + Find the global error estimate. + Change the timestep based on the error estimate. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +void MoL_ReduceAdaptiveError(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + int redop; + + CCTK_REAL red_local[2], red_global[2]; + + int ierr; + + /* Get global result over all processors */ + redop = CCTK_ReductionHandle ("sum"); + assert (redop >= 0); + + red_local[0] = *Error; + red_local[1] = *Count; + ierr = CCTK_ReduceLocArrayToArray1D + (cctkGH, -1, redop, red_local, red_global, 2, CCTK_VARIABLE_REAL); + assert (ierr == 0); + *Error = red_global[0]; + *Count = red_global[1]; + + /* Calculate L2-norm */ + *Error = sqrt(*Error / *Count); + if (! CCTK_EQUALS(verbose, "none")) + { + CCTK_VInfo (CCTK_THORNSTRING, "Integration accuracy quotient is %g", (double)*Error); + } + + /* Decide whether to accept this step */ + *MoL_Stepsize_Bad = *Error > 1; + + if (*MoL_Stepsize_Bad) + { + /* The error is too large; reject the time step and reduce the + step size */ + cctkGH->cctk_time -= cctkGH->cctk_delta_time; + if (! CCTK_EQUALS(verbose, "none")) + { + CCTK_VInfo (CCTK_THORNSTRING, "*** REJECTING TIME STEP ***"); + } + + cctkGH->cctk_delta_time + = safety_factor * (*Original_Delta_Time) / pow(*Error, 0.2); + if (! CCTK_EQUALS(verbose, "none")) + { + CCTK_VInfo (CCTK_THORNSTRING, "Setting time step to %g", (double)cctkGH->cctk_delta_time); + } + + if (cctkGH->cctk_delta_time < (*Original_Delta_Time) / maximum_decrease) + { + /* No more than a factor of 10 decrease */ + cctkGH->cctk_delta_time = (*Original_Delta_Time) / maximum_decrease; + if (! CCTK_EQUALS(verbose, "none")) + { + CCTK_VInfo (CCTK_THORNSTRING, " Time step reduction too large; clamping time step to %g", (double)cctkGH->cctk_delta_time); + } + } + + cctkGH->cctk_time += cctkGH->cctk_delta_time; + } + else + { + /* The error is acceptable; estimate the next step size */ + *EstimatedDt = + safety_factor * (*Original_Delta_Time) / pow(*Error, 0.25); + + if (*EstimatedDt > (*Original_Delta_Time) * maximum_increase) + { + /* No more than a factor of 5 increase */ + *EstimatedDt = (*Original_Delta_Time) * maximum_increase; + if (! CCTK_EQUALS(verbose, "none")) + { + CCTK_VInfo (CCTK_THORNSTRING, " Time step increase too large; clamping time step to %g", (double)(*EstimatedDt)); + } + } + } +} + + /*@@ + @routine MoL_SetEstimatedDt + @date Thu Jan 27 14:05:08 2005 + @author Ian Hawke + @desc + Actually set the timestep in PostStep. + This avoids problems when the timestep is changed in the middle of the + evolution loop. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +void MoL_SetEstimatedDt(CCTK_ARGUMENTS) +{ + + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + cctkGH->cctk_delta_time = *EstimatedDt; + + if (! CCTK_EQUALS(verbose, "none")) + { + CCTK_VInfo (CCTK_THORNSTRING, "Setting time step to %g", + (double)cctkGH->cctk_delta_time); + } + +} + + /*@@ + @routine MoL_FinishLoop + @date Thu Jan 27 10:24:19 2005 + @author Erik Schnetter + @desc + Loop control if adaptive timestepping is not used. + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +void MoL_FinishLoop(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + /* Keep time step size unchanged */ + *MoL_Stepsize_Bad = 0; +} + +/******************************************************************** + ********************* Local Routines ************************* + ********************************************************************/ diff --git a/src/make.code.defn b/src/make.code.defn index 544c7bd..7683c09 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -17,7 +17,8 @@ SRCS = ChangeType.c \ RHSNaNCheck.c \ SandR.c \ SetTime.c \ - Startup.c + Startup.c \ + StepSize.c # Subdirectories containing source files SUBDIRS = |