diff options
Diffstat (limited to 'CarpetDev/CarpetJacobi/src')
-rw-r--r-- | CarpetDev/CarpetJacobi/src/Jacobi.cc | 694 |
1 files changed, 338 insertions, 356 deletions
diff --git a/CarpetDev/CarpetJacobi/src/Jacobi.cc b/CarpetDev/CarpetJacobi/src/Jacobi.cc index c326b66cd..a4686af0a 100644 --- a/CarpetDev/CarpetJacobi/src/Jacobi.cc +++ b/CarpetDev/CarpetJacobi/src/Jacobi.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/src/Jacobi.cc,v 1.3 2004/01/25 14:57:29 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/src/Jacobi.cc,v 1.1 2003/09/02 14:35:58 schnetter Exp $ #include <cassert> #include <cmath> @@ -14,18 +14,14 @@ #include "util_ErrorCodes.h" #include "util_Table.h" -#include "TATelliptic.h" - #include "carpet.hh" - -#include "th.hh" +#include "TATelliptic.h" namespace Carpet { // TODO: fix this - void CycleTimeLevels (const cGH* cctkGH); - void Restrict (const cGH* cctkGH); + void Restrict (const cGH* cgh); }; @@ -52,6 +48,134 @@ namespace CarpetJacobi { + namespace common { + + const int * var; + const int * res; + int nvars; + int options_table; + calcfunc calcres; + calcfunc applybnds; + void * userdata; + + CCTK_REAL factor; + CCTK_REAL * factors; + + CCTK_REAL norm_count; + CCTK_REAL norm_l2; + + int ierr; + + void call_calcres (cGH * const cctkGH) + { + if (ierr) return; + ierr = calcres (cctkGH, options_table, userdata); + } + + void norm (cGH * const cctkGH) + { + int ierr; + + cGroup groupdata; + ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[0]), &groupdata); + assert (!ierr); + cGroupDynamicData groupdyndata; + ierr = CCTK_GroupDynamicData (cctkGH, CCTK_GroupIndexFromVarI(var[0]), + &groupdyndata); + assert (!ierr); + const int thedim = groupdata.dim; + assert (thedim>=0 && thedim<=dim); + assert (thedim == groupdyndata.dim); + int lsh[dim], nghostzones[dim]; + for (int d=0; d<thedim; ++d) { + lsh[d] = groupdyndata.lsh[d]; + nghostzones[d] = groupdyndata.nghostzones[d]; + } + for (int d=thedim; d<dim; ++d) { + lsh[d] = 1; + nghostzones[d] = 0; + } + for (int d=0; d<dim; ++d) { + assert (lsh[d]>=0); + assert (nghostzones[d] >= 0 && 2*nghostzones[d] <= lsh[d]); + } + + for (int n=0; n<nvars; ++n) { + const CCTK_REAL * restrict resptr + = (const CCTK_REAL *)CCTK_VarDataPtrI (cctkGH, 0, res[n]); + assert (resptr); + const CCTK_REAL fac = factors[n]; + for (int k=nghostzones[2]; k<lsh[2]-nghostzones[2]; ++k) { + for (int j=nghostzones[1]; j<lsh[1]-nghostzones[1]; ++j) { + for (int i=nghostzones[0]; i<lsh[0]-nghostzones[0]; ++i) { + const int ind = i + lsh[0] * (j + lsh[1] * k); + ++norm_count; + // TODO: scale the norm by the resolution? + norm_l2 += pow(fac * resptr[ind], 2); + } + } + } + } + } + + void smooth (cGH * const cctkGH) + { + int ierr; + + cGroup groupdata; + ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[0]), &groupdata); + assert (!ierr); + cGroupDynamicData groupdyndata; + ierr = CCTK_GroupDynamicData (cctkGH, CCTK_GroupIndexFromVarI(var[0]), + &groupdyndata); + assert (!ierr); + const int thedim = groupdata.dim; + assert (thedim>=0 && thedim<=dim); + assert (thedim == groupdyndata.dim); + int lsh[dim], nghostzones[dim]; + for (int d=0; d<thedim; ++d) { + lsh[d] = groupdyndata.lsh[d]; + nghostzones[d] = groupdyndata.nghostzones[d]; + } + for (int d=thedim; d<dim; ++d) { + lsh[d] = 1; + nghostzones[d] = 0; + } + for (int d=0; d<dim; ++d) { + assert (lsh[d]>=0); + assert (nghostzones[d] >= 0 && 2*nghostzones[d] <= lsh[d]); + } + + for (int n=0; n<nvars; ++n) { + CCTK_REAL * restrict varptr + = (CCTK_REAL *)CCTK_VarDataPtrI (cctkGH, 0, var[n]); + assert (varptr); + const CCTK_REAL * restrict resptr + = (const CCTK_REAL *)CCTK_VarDataPtrI (cctkGH, 0, res[n]); + assert (resptr); + assert (resptr != varptr); + const CCTK_REAL fac = factor * factors[n]; + for (int k=nghostzones[2]; k<lsh[2]-nghostzones[2]; ++k) { + for (int j=nghostzones[1]; j<lsh[1]-nghostzones[1]; ++j) { + for (int i=nghostzones[0]; i<lsh[0]-nghostzones[0]; ++i) { + const int ind = i + lsh[0] * (j + lsh[1] * k); + varptr[ind] += fac * resptr[ind]; + } + } + } + } + } + + void call_applybnds (cGH * const cctkGH) + { + if (ierr) return; + ierr = applybnds (cctkGH, options_table, userdata); + } + + } // namespace common + + + void CarpetJacobi_register () { int const ierr = TATelliptic_RegisterSolver @@ -70,6 +194,7 @@ namespace CarpetJacobi { const calcfunc applybnds, void * const userdata) { + DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; // Error control @@ -87,19 +212,19 @@ namespace CarpetJacobi { assert (res); for (int n=0; n<nvars; ++n) { assert (var[n] >= 0 && var[n] < CCTK_NumVars()); - assert (CCTK_GroupTypeFromVarI(var[n]) == CCTK_GF); - assert (CCTK_GroupDimFromVarI(var[n]) == dim); + assert (CCTK_GroupTypeFromVarI(var[n]) == CCTK_GF + || CCTK_GroupTypeFromVarI(var[n]) == CCTK_ARRAY); + assert (CCTK_GroupDimFromVarI(var[n]) <= dim); assert (CCTK_VarTypeI(var[n]) == CCTK_VARIABLE_REAL); - assert (CCTK_QueryGroupStorageI(cctkGH, - CCTK_GroupIndexFromVarI(var[n]))); + assert (CCTK_QueryGroupStorageI(cctkGH, CCTK_GroupIndexFromVarI(var[n]))); } for (int n=0; n<nvars; ++n) { assert (res[n] >= 0 && res[n] < CCTK_NumVars()); - assert (CCTK_GroupTypeFromVarI(res[n]) == CCTK_GF); - assert (CCTK_GroupDimFromVarI(res[n]) == dim); + assert (CCTK_GroupTypeFromVarI(res[n]) == CCTK_GF + || CCTK_GroupTypeFromVarI(res[n]) == CCTK_ARRAY); + assert (CCTK_GroupDimFromVarI(res[n]) <= dim); assert (CCTK_VarTypeI(res[n]) == CCTK_VARIABLE_REAL); - assert (CCTK_QueryGroupStorageI(cctkGH, - CCTK_GroupIndexFromVarI(res[n]))); + assert (CCTK_QueryGroupStorageI(cctkGH, CCTK_GroupIndexFromVarI(res[n]))); } for (int n=0; n<nvars; ++n) { assert (var[n] != res[n]); @@ -111,7 +236,47 @@ namespace CarpetJacobi { } } +#if 0 + // Get domain description + cGroup groupdata; + ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[0]), &groupdata); + assert (!ierr); + cGroupDynamicData groupdyndata; + ierr = CCTK_GroupDynamicData (cctkGH, CCTK_GroupIndexFromVarI(var[0]), + &groupdyndata); + assert (!ierr); + const int thedim = groupdata.dim; + assert (thedim>=0 && thedim<=dim); + assert (thedim == groupdyndata.dim); + int lsh[dim], nghostzones[dim]; + for (int d=0; d<thedim; ++d) { + lsh[d] = groupdyndata.lsh[d]; + nghostzones[d] = groupdyndata.nghostzones[d]; + } + for (int d=thedim; d<dim; ++d) { + lsh[d] = 1; + nghostzones[d] = 0; + } + for (int d=0; d<dim; ++d) { + assert (lsh[d]>=0); + assert (nghostzones[d] >= 0 && 2*nghostzones[d] <= lsh[d]); + } + // Check all variables + for (int n=0; n<nvars; ++n) { + ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[n]), &groupdata); + assert (!ierr); + ierr = CCTK_GroupDynamicData + (cctkGH, CCTK_GroupIndexFromVarI(var[n]), &groupdyndata); + assert (!ierr); + assert (groupdata.dim == thedim); + assert (groupdyndata.dim == thedim); + for (int d=0; d<thedim; ++d) { + assert (groupdyndata.lsh[d] == lsh[d]); + assert (groupdyndata.nghostzones[d] == nghostzones[d]); + } + } +#endif assert (options_table >= 0); @@ -145,8 +310,7 @@ namespace CarpetJacobi { assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY); ierr = Util_TableGetReal (options_table, &deceleration, "deceleration"); assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY); - ierr = Util_TableGetRealArray (options_table, - nvars, &factors[0], "factors"); + ierr = Util_TableGetRealArray (options_table, nvars, &factors[0], "factors"); assert (ierr==nvars || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY); assert (maxiters >= 0); @@ -160,375 +324,193 @@ namespace CarpetJacobi { assert (factors[n] != 0); } + // Check state + // (want global mode, or at least level mode) + assert (reflevel==-1 || component==-1); + const int saved_reflevel = reflevel; + const int saved_mglevel = mglevel; + set_mglevel ((cGH *)cctkGH, -1); + set_reflevel ((cGH *)cctkGH, -1); + // Fill common block + common::var = var; + common::res = res; + common::nvars = nvars; + common::options_table = options_table; + common::calcres = calcres; + common::applybnds = applybnds; + common::userdata = userdata; + common::factor = factor; + common::factors = &factors[0]; - assert (is_global_mode() || is_level_mode()); - - // The level to solve on - const int solve_level = is_level_mode() ? reflevel : reflevels-1; + // Init statistics + CCTK_REAL norm2 = HUGE_VAL; + CCTK_REAL speed = 0.0; + CCTK_REAL throttle = 1.0; + bool did_hit_min = false; + bool did_hit_max = false; + const time_t starttime = time(0); + time_t nexttime = starttime; - if (solve_level < solve_minlevel) { - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, - "Did not solve because the level is too coarse."); - } - Util_TableSetInt (options_table, 0, "iters"); - Util_TableSetReal (options_table, -1.0, "error"); - - // Return as if the solving had not converged - return -1; + if (verbose || veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, "Solving through adaptive Jacobi iterations"); } -#warning "TODO: assert that all levels are at the same time" - - // Switch to global mode - BEGIN_GLOBAL_MODE(cctkGH) { + int iter; + for (iter=0; iter<=maxiters; ++iter) { - // Reset the initial data time - global_time = 0; - delta_time = 1; - * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time; - * const_cast<CCTK_REAL *> (& cctkGH->cctk_delta_time) = delta_time; - BEGIN_REFLEVEL_LOOP(cctkGH) { - * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time; - for (int m=0; m<maps; ++m) { - vtt[m]->set_time (reflevel, mglevel, 0); - vtt[m]->set_delta - (reflevel, mglevel, 1.0 / ipow(reflevelfact, reflevelpower)); + // Calculate residual + common::ierr = 0; + CallLocalFunction ((cGH *)cctkGH, common::call_calcres); + ierr = common::ierr; + if (ierr != 0) { + if (verbose || veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, "Residual calculation reported error; aborting solver."); } - } END_REFLEVEL_LOOP; + Util_TableSetInt (options_table, iter, "iters"); + Util_TableDeleteKey (options_table, "error"); + goto done; + } + // Save old norm + const CCTK_REAL old_norm2 = norm2; + // Calculate norm + common::norm_count = 0; + common::norm_l2 = 0; + CallLocalFunction ((cGH *)cctkGH, common::norm); + const int sum_handle = CCTK_ReductionArrayHandle ("sum"); + assert (sum_handle >= 0); + CCTK_REAL reduce_in[2], reduce_out[2]; + reduce_in[0] = common::norm_count; + reduce_in[1] = common::norm_l2; + ierr = CCTK_ReduceLocArrayToArray1D (cctkGH, -1, sum_handle, + reduce_in, reduce_out, 2, + CCTK_VARIABLE_REAL); + norm2 = sqrt(reduce_out[1] / reduce_out[0]); - // Init statistics - vector<CCTK_REAL> norm_counts (solve_level+1); - vector<CCTK_REAL> norm_l2s (solve_level+1); - CCTK_REAL norm2 = HUGE_VAL; - CCTK_REAL speed = 0.0; - CCTK_REAL throttle = 1.0; - bool did_hit_min = false; - bool did_hit_max = false; - const time_t starttime = time(0); - time_t nexttime = starttime; + // Calculate convergence speed + speed = old_norm2 / norm2; + // Log output if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, - "Solving through adaptive Jacobi iterations"); - } - - const int icoarsestep - = ipow (mglevelfact * maxreflevelfact, reflevelpower); - const int istep - = ipow (mglevelfact * maxreflevelfact / ipow(reffact, solve_level), - reflevelpower); - int iter = istep; - for (;; iter+=istep) { - - global_time = 1.0 * iter / ipow(maxreflevelfact, reflevelpower); - * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time; -// 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 * (maxreflevelfact/reflevelfact), - reflevelpower); - if (reflevel <= solve_level && (iter-istep) % do_every == 0) { - - // Advance time - for (int m=0; m<maps; ++m) { - vtt[m]->advance_time (reflevel, mglevel); - } - * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) - = (1.0 * (iter - istep + do_every) - / ipow(maxreflevelfact, reflevelpower)); - CycleTimeLevels (cctkGH); -// cout << "CJ residual iter " << iter << " reflevel " << reflevel << " time " << cctkGH->cctk_time << flush << endl; - - // Calculate residual - BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { - BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { - if (! ierr) { - ierr = calcres (cctkGH, options_table, userdata); - } - } END_LOCAL_COMPONENT_LOOP; - } END_MAP_LOOP; - - // Smooth and apply boundary conditions - CCTK_REAL levfac = 0; - for (int d=0; d<dim; ++d) { - levfac += 1 / - pow (cctkGH->cctk_delta_space[d] / cctkGH->cctk_levfac[d], 2); - } - levfac = 1 / levfac; - - BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { - BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { - if (! ierr) { - - for (int n=0; n<nvars; ++n) { - CCTK_REAL * restrict varptr - = (CCTK_REAL *)CCTK_VarDataPtrI (cctkGH, 0, var[n]); - assert (varptr); - const CCTK_REAL * restrict resptr - = (const CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, res[n]); - assert (resptr); - assert (resptr != varptr); - const CCTK_REAL fac = factor * factors[n] * levfac; - assert (cctkGH->cctk_dim == 3); - for (int k=cctkGH->cctk_nghostzones[2]; - k<cctkGH->cctk_lsh[2]-cctkGH->cctk_nghostzones[2]; - ++k) { - for (int j=cctkGH->cctk_nghostzones[1]; - j<cctkGH->cctk_lsh[1]-cctkGH->cctk_nghostzones[1]; - ++j) { - for (int i=cctkGH->cctk_nghostzones[0]; - i<cctkGH->cctk_lsh[0]-cctkGH->cctk_nghostzones[0]; - ++i) { - const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k); - varptr[ind] += fac * resptr[ind]; - } - } - } - } // for n - - // Apply boundary conditions - ierr = applybnds (cctkGH, options_table, userdata); - - } - } END_LOCAL_COMPONENT_LOOP; - } END_MAP_LOOP; - - } - } END_REFLEVEL_LOOP; - - if (ierr) { - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, - "Residual calculation or boundary conditions reported error; aborting solver."); - } - Util_TableSetInt (options_table, iter, "iters"); - Util_TableDeleteKey (options_table, "error"); - goto done; - } - - // Restrict and calculate norm - ierr = 0; - BEGIN_REVERSE_REFLEVEL_LOOP(cctkGH) { - const int do_every - = ipow (mglevelfact * (maxreflevelfact/reflevelfact), - reflevelpower); - if (reflevel <= solve_level && iter % do_every == 0) { - -// cout << "CJ restrict iter " << iter << " reflevel " << reflevel << " time " << cctkGH->cctk_time << flush << endl; - - if (! ierr) { - if (reflevel < solve_level) { - Restrict (cctkGH); - } - } - - BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { - BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { - if (! ierr) { - ierr = applybnds (cctkGH, options_table, userdata); - } - } END_LOCAL_COMPONENT_LOOP; - } END_MAP_LOOP; - - // Initialise norm - CCTK_REAL norm_count = 0; - CCTK_REAL norm_l2 = 0; - - BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { - BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { - if (! ierr) { - - // Calculate norm - for (int n=0; n<nvars; ++n) { - const CCTK_REAL * restrict resptr - = (const CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, res[n]); - assert (resptr); - const CCTK_REAL fac = factors[n]; - assert (cctkGH->cctk_dim == 3); - for (int k=cctkGH->cctk_nghostzones[2]; - k<cctkGH->cctk_lsh[2]-cctkGH->cctk_nghostzones[2]; - ++k) { - for (int j=cctkGH->cctk_nghostzones[1]; - j<cctkGH->cctk_lsh[1]-cctkGH->cctk_nghostzones[1]; - ++j) { - for (int i=cctkGH->cctk_nghostzones[0]; - i<cctkGH->cctk_lsh[0]-cctkGH->cctk_nghostzones[0]; - ++i) { - const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k); - ++norm_count; - // TODO: scale the norm by the resolution? - norm_l2 += pow(fac * resptr[ind], 2); - } - } - } - } // for n - - } - } END_LOCAL_COMPONENT_LOOP; - } END_MAP_LOOP; - - const int sum_handle = CCTK_ReductionArrayHandle ("sum"); - assert (sum_handle >= 0); - CCTK_REAL reduce_in[2], reduce_out[2]; - reduce_in[0] = norm_count; - reduce_in[1] = norm_l2; - ierr = CCTK_ReduceLocArrayToArray1D - (cctkGH, -1, sum_handle, reduce_in, reduce_out, 2, - CCTK_VARIABLE_REAL); - norm_counts.at(reflevel) = reduce_out[0]; - norm_l2s.at(reflevel) = reduce_out[1]; - -#warning "TODO" -#if 0 - CCTK_OutputVarAs (cctkGH, "WaveToyFO::phi", "phi-ell"); - CCTK_OutputVarAs (cctkGH, "IDSWTEsimpleFO::residual", "residual-ell"); -#endif - - } - } END_REVERSE_REFLEVEL_LOOP; - - if (ierr) { - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, - "Residual calculation or boundary conditions reported error; aborting solver."); - } - Util_TableSetInt (options_table, iter, "iters"); - Util_TableDeleteKey (options_table, "error"); - goto done; - } - - // Save old norm - const CCTK_REAL old_norm2 = norm2; - - // Calculate new norm - CCTK_REAL norm_count = 0; - CCTK_REAL norm_l2 = 0; - for (int rl=0; rl<=solve_level; ++rl) { - norm_count += norm_counts[rl]; - norm_l2 += norm_l2s[rl]; - } - norm2 = sqrt(norm_l2 / norm_count); - - // Calculate convergence speed - speed = old_norm2 / norm2; - - // Log output - if (verbose || veryverbose) { - const time_t currenttime = time(0); - if ((iter % icoarsestep == 0 && iter >= maxiters) + const time_t currenttime = time(0); + if (iter == maxiters #if HAVE_ISNAN - || isnan(norm2) + || isnan(norm2) #endif - || norm2 <= minerror - || (iter % outevery == 0 && currenttime >= nexttime)) { - if (verbose || veryverbose) { - if (did_hit_min) { - CCTK_VInfo (CCTK_THORNSTRING, - "Convergence factor became too small: artificially kept at minfactor (may lead to instability)"); - } - if (did_hit_max) { - CCTK_VInfo (CCTK_THORNSTRING, - "Convergence factor became too large: artificially kept at maxfactor (may lead to slower convergence)"); - } - did_hit_min = false; - did_hit_max = false; - } - if (veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, - "Iteration %d, time %g, factor %g, throttle %g, residual %g, speed %g", - iter, double(currenttime-starttime), - double(factor), double(throttle), double(norm2), - double(speed)); - } else { + || norm2 <= minerror + || (iter % outevery == 0 && currenttime >= nexttime)) { + if (verbose || veryverbose) { + if (did_hit_min) { CCTK_VInfo (CCTK_THORNSTRING, - "Iteration %d, time %g, residual %g", - iter, double(currenttime-starttime), double(norm2)); + "Convergence factor became too small: artificially kept at minfactor (may lead to instability)"); } - if (outeveryseconds > 0) { - while (nexttime <= currenttime) nexttime += outeveryseconds; + if (did_hit_max) { + CCTK_VInfo (CCTK_THORNSTRING, + "Convergence factor became too large: artificially kept at maxfactor (may lead to slower convergence)"); } + did_hit_min = false; + did_hit_max = false; } - } - -#if HAVE_ISNAN - // Exit if things go bad - if (isnan(norm2)) { - if (verbose || veryverbose) { + if (veryverbose) { CCTK_VInfo (CCTK_THORNSTRING, - "Encountered NaN. Aborting solver."); + "Iteration %d, time %g, factor %g, throttle %g, residual %g, speed %g", + iter, (double)(currenttime-starttime), (double)factor, (double)throttle, (double)norm2, (double)speed); + } else { + CCTK_VInfo (CCTK_THORNSTRING, + "Iteration %d, time %g, residual %g", + iter, (double)(currenttime-starttime), (double)norm2); } - break; - } -#endif - - // Return when desired accuracy is reached - if (norm2 <= minerror) { - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, "Done solving."); + if (outeveryseconds > 0) { + while (nexttime <= currenttime) nexttime += outeveryseconds; } - Util_TableSetInt (options_table, iter, "iters"); - Util_TableSetReal (options_table, norm2, "error"); - ierr = 0; - goto done; - } - - // Exit if the maximum number of iterations has been reached - if (iter % icoarsestep == 0 && iter >= maxiters) { - break; - } - - // Adjust speed - if (speed > 1.0) { - throttle = acceleration; - } else { - throttle = deceleration; - } - factor *= throttle; - if (factor < minfactor) { - did_hit_min = true; - factor = minfactor; } - if (factor > maxfactor) { - did_hit_max = true; - factor = maxfactor; + } + + if (iter == maxiters) break; + +#if HAVE_ISNAN + // Exit if things go bad + if (isnan(norm2)) { + if (verbose || veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, "Encountered NaN. Aborting solver."); } - - } // for iter + break; + } +#endif - // Did not solve - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, "Did not converge."); + // Return when desired accuracy is reached + if (norm2 <= minerror) { + if (verbose || veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, "Done solving."); + } + Util_TableSetInt (options_table, iter, "iters"); + Util_TableSetReal (options_table, norm2, "error"); + ierr = 0; + goto done; } - Util_TableSetInt (options_table, iter, "iters"); - Util_TableSetReal (options_table, norm2, "error"); - ierr = -1; - - done: ; + // Adjust speed + if (speed > 1.0) { + throttle = acceleration; + } else { + throttle = deceleration; + } + factor *= throttle; + if (factor < minfactor) { + did_hit_min = true; + factor = minfactor; + } + if (factor > maxfactor) { + did_hit_max = true; + factor = maxfactor; + } + // Smooth + CallLocalFunction ((cGH *)cctkGH, common::smooth); - // Reset the initial time -#warning "TODO: reset the initial time a bit more intelligently" - global_time = 0; - delta_time = 1; - * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time; - * const_cast<CCTK_REAL *> (& cctkGH->cctk_delta_time) = delta_time; - BEGIN_REFLEVEL_LOOP(cctkGH) { - * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time; - for (int m=0; m<maps; ++m) { - vtt[m]->set_time (reflevel, mglevel, 0); - vtt[m]->set_delta (reflevel, mglevel, 1.0 / reflevelfact); + // Apply boundary conditions + common::ierr = 0; + CallLocalFunction ((cGH *)cctkGH, common::call_applybnds); + ierr = common::ierr; + if (ierr != 0) { + if (verbose || veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, "Boundary conditions reported error; aborting solver."); } - } END_REFLEVEL_LOOP; + Util_TableSetInt (options_table, iter, "iters"); + Util_TableDeleteKey (options_table, "error"); + goto done; + } - } END_GLOBAL_MODE; +#if 0 + // Restrict + BEGIN_REVERSE_REFLEVEL_LOOP(cctkGH) { + BEGIN_MGLEVEL_LOOP(cctkGH) { + Restrict (cctkGH); + // TODO: do something here +// CCTK_ScheduleTraverse ("CCTK_POSTRESTRICT", (cGH *)cctkGH, CallFunction); + CallLocalFunction ((cGH *)cctkGH, common::call_applybnds); + assert (! common::ierr); + } END_MGLEVEL_LOOP; + } END_REVERSE_REFLEVEL_LOOP; +#endif + + } // for iter + + // Did not solve + if (verbose || veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, "Did not converge."); + } + Util_TableSetInt (options_table, iter, "iters"); + Util_TableSetReal (options_table, norm2, "error"); + ierr = -1; + + done: + + // Restore state + set_reflevel ((cGH *)cctkGH, saved_reflevel); + set_mglevel ((cGH *)cctkGH, saved_mglevel); return ierr; } |