diff options
author | eschnett <eschnett@578cdeb0-5ea1-4b81-8215-5a3b8777ee0b> | 2011-03-25 21:44:37 +0000 |
---|---|---|
committer | eschnett <eschnett@578cdeb0-5ea1-4b81-8215-5a3b8777ee0b> | 2011-03-25 21:44:37 +0000 |
commit | 42b20a1be5e2d6baf633cd68adff4991e869c154 (patch) | |
tree | 521fabb118f7f65a6deabb290f04ccd5e4284d7e | |
parent | 61a33fbee9c1b55997fe7f2e0b8c9c25f4bca9ba (diff) |
Improve roundoff accuracy of efficient RK4
I just encountered a weird case where evolving Minkowski with the
efficient RK4 implementation was not static. It turns out that the
reason was floating-point round-off error in RK4. During each of the
substeps, RK4 adds terms with factors of 1/3 and 2/3, and for the
final step, adds a term with a factor of 4/3. Due to round-off, these
terms do not cancel exactly, so that the time-evolved lapse is
slightly different from 1 (although all the lapse RHS terms are 0).
The same goes for the diagonal terms of the metric.
The attached patch keeps track of the round-off error when adding
these terms, leading to a more accurate integration result.
The problem described above is not always visible, but also depends on
optimisation settings.
git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/MoL/trunk@145 578cdeb0-5ea1-4b81-8215-5a3b8777ee0b
-rw-r--r-- | src/RK4.c | 19 |
1 files changed, 15 insertions, 4 deletions
@@ -86,8 +86,8 @@ CCTK_WARN(0, "not implemented"); #endif static CCTK_INT scratchspace_firstindex = -99; - CCTK_INT index, var, scratchstep, alphaindex, scratchindex; - CCTK_INT totalsize, singlearraysize; + CCTK_INT index, var; + CCTK_INT totalsize; CCTK_REAL alpha, beta; CCTK_REAL * restrict UpdateVar; CCTK_REAL * restrict OldVar; @@ -96,6 +96,11 @@ CCTK_WARN(0, "not implemented"); CCTK_INT arrayscratchlocation; + /* 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; + totalsize = 1; for (arraydim = 0; arraydim < cctk_dim; arraydim++) { @@ -126,6 +131,12 @@ CCTK_WARN(0, "not implemented"); beta = 1.0 / 6.0; } + if (MoL_Intermediate_Steps == (*MoL_Intermediate_Step)) + { + sum_alpha = 0.0; + } + sum_alpha += alpha; + /* FIXME */ #ifdef MOLDOESCOMPLEX @@ -197,7 +208,7 @@ CCTK_WARN(0, "not implemented"); #pragma omp parallel for for (index = 0; index < totalsize; index++) { - UpdateVar[index] += ScratchVar[index] - 4.0 / 3.0 * OldVar[index]; + UpdateVar[index] += ScratchVar[index] - (sum_alpha - 1.0) * OldVar[index]; } } @@ -276,7 +287,7 @@ CCTK_WARN(0, "not done"); UpdateVar[index] += ScratchVar[index] - 4.0 / 3.0 * OldVar[index]; } } - arrayscratchlocation += arraytotalsize; + arrayscratchlocation += arraytotalsize; } return; |