diff options
author | Erik Schnetter <schnetter@cct.lsu.edu> | 2005-12-22 19:55:00 +0000 |
---|---|---|
committer | Erik Schnetter <schnetter@cct.lsu.edu> | 2005-12-22 19:55:00 +0000 |
commit | 2453e3b19d421340c616bb706852f422f4e5da95 (patch) | |
tree | 9a017bddca5b26189a03551fc3372b64ed5caf57 /CarpetDev/CarpetJacobi | |
parent | 92fe2ab05a26d77228a8b76ca1a446cb0fd044c8 (diff) |
CarpetJacobi: Correct time level cycling
When the solution has multiple time levels, the solution needs to be
copyied from the past to the current time level after time level
cycling.
Add some debugging code, disabled by default.
Correct handling of restrict qualifier.
darcs-hash:20051222195540-dae7b-5c276af06aed498dbbea3f23f589641377c2bfa7.gz
Diffstat (limited to 'CarpetDev/CarpetJacobi')
-rw-r--r-- | CarpetDev/CarpetJacobi/src/Jacobi.cc | 60 |
1 files changed, 44 insertions, 16 deletions
diff --git a/CarpetDev/CarpetJacobi/src/Jacobi.cc b/CarpetDev/CarpetJacobi/src/Jacobi.cc index 4e326d74c..7a01b0f49 100644 --- a/CarpetDev/CarpetJacobi/src/Jacobi.cc +++ b/CarpetDev/CarpetJacobi/src/Jacobi.cc @@ -1,6 +1,7 @@ #include <cassert> #include <cmath> #include <cstdio> +#include <cstring> #include <ctime> #include <iostream> #include <vector> @@ -28,12 +29,17 @@ namespace Carpet { +#undef restrict #ifdef CCTK_CXX_RESTRICT # define restrict CCTK_CXX_RESTRICT #endif +#define DEBUG 0 + + + namespace CarpetJacobi { using namespace std; @@ -233,14 +239,14 @@ namespace CarpetJacobi { global_time = 1.0 * iter / ipow (maxval (maxspacereflevelfact), reflevelpower); * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time; -// cout << "CJ global iter " << iter << " time " << global_time << flush << endl; + if (DEBUG) cout << "CJ global iter " << iter << " time " << global_time << flush << endl; // Calculate residual, smooth, and apply boundary conditions ierr = 0; BEGIN_REFLEVEL_LOOP(cctkGH) { const int do_every - = ipow (mglevelfact * maxval (maxspacereflevelfact - / spacereflevelfact), + = ipow (mglevelfact + * maxval (maxspacereflevelfact / spacereflevelfact), reflevelpower); if (reflevel <= solve_level && (iter-istep) % do_every == 0) { @@ -252,7 +258,24 @@ namespace CarpetJacobi { = (1.0 * (iter - istep + do_every) / ipow (maxval (maxspacereflevelfact), reflevelpower)); CycleTimeLevels (cctkGH); -// cout << "CJ residual iter " << iter << " reflevel " << reflevel << " time " << cctkGH->cctk_time << flush << endl; + if (DEBUG) cout << "CJ residual iter " << iter << " reflevel " << reflevel << " time " << cctkGH->cctk_time << flush << endl; + + // Advance time levels + BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { + BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { + for (int n=0; n<nvars; ++n) { + if (CCTK_ActiveTimeLevelsVI(cctkGH, var[n]) > 1) { + size_t npoints = 1; + for (int d=0; d<cctkGH->cctk_dim; ++d) { + npoints *= cctkGH->cctk_lsh[d]; + } + memcpy (CCTK_VarDataPtrI (cctkGH, 0, var[n]), + CCTK_VarDataPtrI (cctkGH, 1, var[n]), + npoints * sizeof(CCTK_REAL)); + } + } // for n + } END_LOCAL_COMPONENT_LOOP; + } END_MAP_LOOP; // Calculate residual BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { @@ -264,14 +287,16 @@ namespace CarpetJacobi { } END_MAP_LOOP; // Smooth and apply boundary conditions - CCTK_REAL levfac = 0; - for (int d=0; d<dim; ++d) { - levfac += 1 / - ipow (cctkGH->cctk_delta_space[d] / cctkGH->cctk_levfac[d], 2); - } - levfac = 1 / levfac; - BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { + + CCTK_REAL levfac = 0; + for (int d=0; d<dim; ++d) { + CCTK_REAL const delta + = cctkGH->cctk_delta_space[d] / cctkGH->cctk_levfac[d]; + levfac += 1 / ipow (delta, 2); + } + levfac = 1 / levfac; + BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { if (! ierr) { @@ -284,6 +309,7 @@ namespace CarpetJacobi { assert (resptr); assert (resptr != varptr); const CCTK_REAL fac = factor * factors[n] * levfac; + if (DEBUG) cout << "CJ smoothing fac=" << fac << endl; assert (cctkGH->cctk_dim == 3); for (int k=cctkGH->cctk_nghostzones[2]; k<cctkGH->cctk_lsh[2]-cctkGH->cctk_nghostzones[2]; @@ -306,7 +332,7 @@ namespace CarpetJacobi { ierr = applybnds (cctkGH, options_table, userdata); #endif - } + } // if ! ierr } END_LOCAL_COMPONENT_LOOP; } END_MAP_LOOP; @@ -317,7 +343,7 @@ namespace CarpetJacobi { #endif - } + } // if do_every } END_REFLEVEL_LOOP; if (ierr) { @@ -334,12 +360,12 @@ namespace CarpetJacobi { ierr = 0; BEGIN_REVERSE_REFLEVEL_LOOP(cctkGH) { const int do_every - = ipow (mglevelfact * maxval (maxspacereflevelfact - / spacereflevelfact), + = ipow (mglevelfact + * maxval (maxspacereflevelfact / spacereflevelfact), reflevelpower); if (reflevel <= solve_level && iter % do_every == 0) { -// cout << "CJ restrict iter " << iter << " reflevel " << reflevel << " time " << cctkGH->cctk_time << flush << endl; + if (DEBUG) cout << "CJ restrict iter " << iter << " reflevel " << reflevel << " time " << cctkGH->cctk_time << flush << endl; if (! ierr) { if (reflevel < solve_level) { @@ -369,6 +395,7 @@ namespace CarpetJacobi { = (const CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, res[n]); assert (resptr); const CCTK_REAL fac = factors[n]; + if (DEBUG) cout << "CJ norm fac=" << fac << endl; assert (cctkGH->cctk_dim == 3); for (int k=cctkGH->cctk_nghostzones[2]; k<cctkGH->cctk_lsh[2]-cctkGH->cctk_nghostzones[2]; @@ -386,6 +413,7 @@ namespace CarpetJacobi { } } } + if (DEBUG) cout << "CJ norm=" << norm_l2 << endl; } // for n } |