From 2c8448f84da0b9f4ae13cf1a02c71e0c626e4782 Mon Sep 17 00:00:00 2001 From: hawke Date: Thu, 1 Jul 2004 11:02:39 +0000 Subject: Add the Classic RK3 method (as a generic method, so use Generic_Type). Agrees with other RK3's to floating point round off (except at boundaries) for linear case. Uses more storage and is slower than standard RK3 so I don't recommend it. This showed up (so I fixed) a bug with the generic methods when used with Carpet. git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/MoL/trunk@72 578cdeb0-5ea1-4b81-8215-5a3b8777ee0b --- param.ccl | 1 + src/GenericRK.c | 87 +++++++++++++++++++++++++++++++++++----------------- src/IndexArrays.c | 4 +++ src/ParamCheck.c | 8 +++++ src/RKCoefficients.c | 35 +++++++++++++++++++-- src/SetTime.c | 6 ++-- 6 files changed, 108 insertions(+), 33 deletions(-) diff --git a/param.ccl b/param.ccl index 74752ad..e982acb 100644 --- a/param.ccl +++ b/param.ccl @@ -94,6 +94,7 @@ 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" + "Classic RK3" :: "Efficient RK3 - classical version" } "RK" CCTK_INT MoL_Intermediate_Steps "Number of intermediate steps taken by the ODE method" diff --git a/src/GenericRK.c b/src/GenericRK.c index 6c052af..942866b 100644 --- a/src/GenericRK.c +++ b/src/GenericRK.c @@ -108,7 +108,7 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS) #ifdef MOLDOESCOMPLEX - Complex_Delta_Time = CCTK_Cmplx((*Original_Delta_Time), 0); + Complex_Delta_Time = CCTK_Cmplx((*Original_Delta_Time) / cctkGH->cctk_timefac, 0); Complex_beta = CCTK_Cmplx(beta, 0); #endif @@ -122,22 +122,26 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS) EvolvedVariableIndex[var]); RHSVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, RHSVariableIndex[var]); -/* #define MOLDEBUG */ +/* #define MOLDEBUG 1 */ #ifdef MOLDEBUG - printf("In generic RK. Variable %d (%s). RHS %d (%s).\n", + 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])); + CCTK_VarName(RHSVariableIndex[var]), + beta); #endif for (index = 0; index < totalsize; index++) { - UpdateVar[index] = (*Original_Delta_Time) * beta * RHSVar[index]; + UpdateVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac * beta * RHSVar[index]; #ifdef MOLDEBUG - printf("Variable: %d. Index: %d. dt: %f. beta %f. RHS: %f. q: %f.\n", - var, index, (*Original_Delta_Time), beta, RHSVar[index], - UpdateVar[index]); + 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 } @@ -150,6 +154,15 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS) scratchindex = scratchstep - 1; alpha = RKAlphaCoefficients[alphaindex]; +#ifdef MOLDEBUG + printf("In generic RK. Variable %d (%s). RHS %d (%s). step %d. alpha %g.\n", + EvolvedVariableIndex[var], + CCTK_VarName(EvolvedVariableIndex[var]), + RHSVariableIndex[var], + CCTK_VarName(RHSVariableIndex[var]), + scratchstep, + alpha); +#endif if (scratchstep) { @@ -166,9 +179,12 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS) + var + MoL_Num_Evolved_Vars * scratchindex); #ifdef MOLDEBUG - printf("Reading from scratch space, initial address %ld index %d\n", - ScratchVar, (var * MoL_Num_Scratch_Levels + - scratchindex) * totalsize); + if (CCTK_EQUALS(verbose,"extreme")) + { + printf("Reading from scratch space, initial address %ld index %d\n", + ScratchVar, (var * MoL_Num_Scratch_Levels + + scratchindex) * totalsize); + } #endif } else @@ -183,10 +199,13 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS) { UpdateVar[index] += alpha * ScratchVar[index]; #ifdef MOLDEBUG - printf("Variable: %d. Index: %d. step: %d. " - "alpha: %f. Scratch: %f. q: %f.\n", - var, index, (*MoL_Intermediate_Step), alpha, - ScratchVar[index], UpdateVar[index]); + if (CCTK_EQUALS(verbose,"extreme")) + { + printf("Variable: %d. Index: %d. step: %d. " + "alpha: %f. Scratch: %f. q: %f.\n", + var, index, (*MoL_Intermediate_Step), alpha, + ScratchVar[index], UpdateVar[index]); + } #endif } } @@ -225,9 +244,12 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS) { ScratchVar[index] = UpdateVar[index]; #ifdef MOLDEBUG - printf("Variable: %d. Index: %d. step: %d. Scratch: %f.\n", - var, index, (*MoL_Intermediate_Step), ScratchVar[index]); - fflush(stdout); + if (CCTK_EQUALS(verbose,"extreme")) + { + printf("Variable: %d. Index: %d. step: %d. Scratch: %f.\n", + var, index, (*MoL_Intermediate_Step), ScratchVar[index]); + fflush(stdout); + } #endif } } @@ -344,7 +366,7 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS) for (index = 0; index < arraytotalsize; index++) { - UpdateVar[index] = (*Original_Delta_Time) * beta * RHSVar[index]; + UpdateVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac * beta * RHSVar[index]; } for (scratchstep = 0; @@ -363,9 +385,12 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS) (MoL_Max_Evolved_Array_Size+1) + arrayscratchlocation]; #ifdef MOLDEBUG - printf("Reading from scratch space, initial address %ld index %d\n", - ScratchVar, scratchindex*(MoL_Max_Evolved_Array_Size+1) + - arrayscratchlocation); + if (CCTK_EQUALS(verbose,"extreme")) + { + printf("Reading from scratch space, initial address %ld index %d\n", + ScratchVar, scratchindex*(MoL_Max_Evolved_Array_Size+1) + + arrayscratchlocation); + } #endif } else @@ -380,10 +405,13 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS) { UpdateVar[index] += alpha * ScratchVar[index]; #ifdef MOLDEBUG - printf("Variable: %d. Index: %d. step: %d. " - "alpha: %f. Scratch: %f. q: %f.\n", - var, index, (*MoL_Intermediate_Step), - alpha, ScratchVar[index], UpdateVar[index]); + if (CCTK_EQUALS(verbose,"extreme")) + { + printf("Variable: %d. Index: %d. step: %d. " + "alpha: %f. Scratch: %f. q: %f.\n", + var, index, (*MoL_Intermediate_Step), + alpha, ScratchVar[index], UpdateVar[index]); + } #endif } } @@ -433,8 +461,11 @@ void MoL_GenericRKAdd(CCTK_ARGUMENTS) { ScratchVar[index] = UpdateVar[index]; #ifdef MOLDEBUG - printf("Variable: %d. Index: %d. step: %d. Scratch: %f.\n", - var, index, (*MoL_Intermediate_Step), ScratchVar[index]); + if (CCTK_EQUALS(verbose,"extreme")) + { + printf("Variable: %d. Index: %d. step: %d. Scratch: %f.\n", + var, index, (*MoL_Intermediate_Step), ScratchVar[index]); + } #endif } arrayscratchlocation += arraytotalsize; diff --git a/src/IndexArrays.c b/src/IndexArrays.c index e0df297..096a745 100644 --- a/src/IndexArrays.c +++ b/src/IndexArrays.c @@ -282,6 +282,10 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS) sprintf(infoline, "Generic method, options:\n %s\n", Generic_Method_Descriptor); } + else if (CCTK_EQUALS(Generic_Type,"Classic RK3")) + { + sprintf(infoline, "Classic Runge-Kutta 3"); + } else { CCTK_WARN(0, "Generic_Type not recognized!"); diff --git a/src/ParamCheck.c b/src/ParamCheck.c index 47c573a..5316c9f 100644 --- a/src/ParamCheck.c +++ b/src/ParamCheck.c @@ -75,6 +75,13 @@ int MoL_ParamCheck(CCTK_ARGUMENTS) "of scratch levels must be at least the " "number of intermediate steps - 1"); } + if ( (CCTK_Equals(Generic_Type, "Classic RK3"))&& + ((!(MoL_Intermediate_Steps == 3))||(!(MoL_Num_Scratch_Levels > 1))) ) + { + CCTK_PARAMWARN("When using the classic RK3 evolver the " + "number of intermediate steps must be 3 " + "and the number of scratch levels at least 2"); + } if (CCTK_Equals(Generic_Type, "Table")) { options_table = @@ -121,6 +128,7 @@ int MoL_ParamCheck(CCTK_ARGUMENTS) "number of intermediate steps must be 3"); } + return 0; } diff --git a/src/RKCoefficients.c b/src/RKCoefficients.c index 905fb60..57ac38a 100644 --- a/src/RKCoefficients.c +++ b/src/RKCoefficients.c @@ -73,7 +73,36 @@ int MoL_SetupRKCoefficients(CCTK_ARGUMENTS) CCTK_INT ierr, options_table; - if (CCTK_Equals(Generic_Type,"ICN")) + if (CCTK_EQUALS(Generic_Type,"Classic RK3")) + { + if (MoL_Num_Scratch_Levels != 2) + { + CCTK_WARN(0, "For Classic RK3, MoL_Num_Scratch_Levels " + "should be at least 2"); + } + if (MoL_Intermediate_Steps != 3) + { + CCTK_WARN(0, "For Classic RK3, MoL_Intermediate_Steps " + "should be at least 3"); + } + for (i = 0; i < MoL_Intermediate_Steps; i++) + { + for (j = 0; j < MoL_Num_Scratch_Levels + 1; j++) + { + RKAlphaCoefficients[i * MoL_Intermediate_Steps + j] = 0.0; + } + RKBetaCoefficients[i] = 0.0; + } + RKAlphaCoefficients[0] = 1.0; + RKAlphaCoefficients[3] = 1.0; + RKAlphaCoefficients[6] = 1.0 / 9.0; + RKAlphaCoefficients[7] = 4.0 / 9.0; + RKAlphaCoefficients[8] = 4.0 / 9.0; + RKBetaCoefficients[0] = 0.5; + RKBetaCoefficients[1] = 0.75; + RKBetaCoefficients[2] = 4.0 / 9.0; + } + else if (CCTK_EQUALS(Generic_Type,"ICN")) { for (i = 0; i < MoL_Intermediate_Steps; i++) { @@ -92,7 +121,7 @@ int MoL_SetupRKCoefficients(CCTK_ARGUMENTS) } } } - else if (CCTK_Equals(Generic_Type,"RK")) + else if (CCTK_EQUALS(Generic_Type,"RK")) { if (MoL_Num_Scratch_Levels < MoL_Intermediate_Steps - 1) { @@ -152,7 +181,7 @@ int MoL_SetupRKCoefficients(CCTK_ARGUMENTS) "with MoL_Intermediate_Steps greater than 4"); } } - else if (CCTK_Equals(Generic_Type,"Table")) + else if (CCTK_EQUALS(Generic_Type,"Table")) { if (MoL_Num_Scratch_Levels < MoL_Intermediate_Steps - 1) { diff --git a/src/SetTime.c b/src/SetTime.c index 76a06a4..94380b1 100644 --- a/src/SetTime.c +++ b/src/SetTime.c @@ -19,6 +19,8 @@ static const char *rcsid = "$Header$"; CCTK_FILEVERSION(CactusBase_MoL_SetTime_c); +/* #define MOLDEBUG 1 */ + /******************************************************************** ********************* Local Data Types *********************** ********************************************************************/ @@ -192,7 +194,7 @@ int MoL_ResetTime(CCTK_ARGUMENTS) } else if (*MoL_Intermediate_Step == 1) { - cctkGH->cctk_time = (*Original_Time) + + cctkGH->cctk_time = (*Original_Time) - 0.5*(*Original_Delta_Time)/cctkGH->cctk_timefac; } } @@ -254,7 +256,7 @@ int MoL_ResetDeltaTime(CCTK_ARGUMENTS) { cctkGH->cctk_delta_time = RKBetaCoefficients[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)] * - (*Original_Delta_Time)/cctkGH->cctk_timefac; + (*Original_Delta_Time); } else if (CCTK_EQUALS(ODE_Method,"RK2")) { -- cgit v1.2.3