diff options
author | schnetter <schnetter@578cdeb0-5ea1-4b81-8215-5a3b8777ee0b> | 2006-02-19 00:31:54 +0000 |
---|---|---|
committer | schnetter <schnetter@578cdeb0-5ea1-4b81-8215-5a3b8777ee0b> | 2006-02-19 00:31:54 +0000 |
commit | c4896403254a6de126501a42f9cb046538a44727 (patch) | |
tree | 5044935d014dab7bca78235a52a54f64dad45402 /src/RK45.c | |
parent | a7134900411d8adb1cb1f4d6962dfcf892983326 (diff) |
Implement RK45 Cash-Karp integrator in MoL.
This integrator is similar to the existing RK45 integrator. It also
supports adaptive time stepping, but uses slightly different
coefficients. The Numerical Recipes say that it has "slightly better
error properties".
git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/MoL/trunk@110 578cdeb0-5ea1-4b81-8215-5a3b8777ee0b
Diffstat (limited to 'src/RK45.c')
-rw-r--r-- | src/RK45.c | 61 |
1 files changed, 54 insertions, 7 deletions
@@ -67,11 +67,9 @@ void MoL_RK45Add(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, scratchstep, alphaindex, scratchindex; + CCTK_INT index, var, scratchstep; CCTK_INT totalsize; CCTK_REAL * restrict UpdateVar; @@ -84,7 +82,8 @@ void MoL_RK45Add(CCTK_ARGUMENTS) CCTK_REAL beta, gamma, gamma_error; - static const CCTK_REAL beta_array[5][5] = { + /* Fehlberg coefficients */ + static const CCTK_REAL beta_array_F[5][5] = { { 0.25, 0.0, 0.0, 0.0, 0.0 }, { 0.09375, 0.28125, 0.0, 0.0, 0.0 }, { 0.879380974055530268548020027310, -3.27719617660446062812926718252, 3.32089212562585343650432407829, 0.0, 0.0 }, @@ -92,7 +91,7 @@ void MoL_RK45Add(CCTK_ARGUMENTS) { -0.296296296296296296296296296296, 2.0, -1.38167641325536062378167641326, 0.452972709551656920077972709552, -0.275 } }; - static const CCTK_REAL gamma_array[6] = + static const CCTK_REAL gamma_array_F[6] = { 0.118518518518518518518518518519, 0.0, 0.518986354775828460038986354776, @@ -101,7 +100,7 @@ void MoL_RK45Add(CCTK_ARGUMENTS) 0.0363636363636363636363636363636 }; - static const CCTK_REAL gammastar_array[6] = + static const CCTK_REAL gammastar_array_F[6] = { 0.115740740740740740740740740741, 0.0, 0.548927875243664717348927875244, @@ -110,6 +109,54 @@ void MoL_RK45Add(CCTK_ARGUMENTS) 0.0 }; + /* Cash-Karp coefficients */ + static const CCTK_REAL beta_array_CK[5][5] = { + { 1.0/5.0, 0.0, 0.0, 0.0, 0.0, }, + { 3.0/40.0, 9.0/40.0, 0.0, 0.0, 0.0, }, + { 3.0/10.0, -9.0/10.0, 6.0/5.0, 0.0, 0.0, }, + { -11.0/54.0, 5.0/2.0, -70.0/27.0, 35.0/27.0, 0.0, }, + { 1631.0/55296.0, 175.0/512.0, 575.0/13824.0, 44275.0/110592.0, 253.0/4096.0, }, + }; + + static const CCTK_REAL gamma_array_CK[6] = { + 37.0/378.0, + 0.0, + 250.0/621.0, + 125.0/594.0, + 0.0, + 512.0/1771.0, + }; + + static const CCTK_REAL gammastar_array_CK[6] = { + 2825.0/27648.0, + 0.0, + 18575.0/48384.0, + 13525.0/55296.0, + 277.0/14336.0, + 1.0/4.0, + }; + + const CCTK_REAL (* restrict beta_array)[5]; + const CCTK_REAL * restrict gamma_array; + const CCTK_REAL * restrict gammastar_array; + + if (CCTK_EQUALS(ODE_Method, "RK45")) + { + beta_array = beta_array_F; + gamma_array = gamma_array_F; + gammastar_array = gammastar_array_F; + } + else if (CCTK_EQUALS(ODE_Method, "RK45CK")) + { + beta_array = beta_array_CK; + gamma_array = gamma_array_CK; + gammastar_array = gammastar_array_CK; + } + else + { + CCTK_WARN (0, "internal error"); + } + totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { |