diff options
Diffstat (limited to 'CarpetDev/CarpetJacobi/src/Jacobi.cc')
| -rw-r--r-- | CarpetDev/CarpetJacobi/src/Jacobi.cc | 698 |
1 files changed, 356 insertions, 342 deletions
diff --git a/CarpetDev/CarpetJacobi/src/Jacobi.cc b/CarpetDev/CarpetJacobi/src/Jacobi.cc index ad673b342..c326b66cd 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.2 2003/11/19 14:05:14 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/src/Jacobi.cc,v 1.3 2004/01/25 14:57:29 schnetter Exp $ #include <cassert> #include <cmath> @@ -14,14 +14,18 @@ #include "util_ErrorCodes.h" #include "util_Table.h" -#include "carpet.hh" #include "TATelliptic.h" +#include "carpet.hh" + +#include "th.hh" + namespace Carpet { // TODO: fix this - void Restrict (const cGH* cgh); + void CycleTimeLevels (const cGH* cctkGH); + void Restrict (const cGH* cctkGH); }; @@ -48,134 +52,6 @@ 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 @@ -194,7 +70,6 @@ namespace CarpetJacobi { const calcfunc applybnds, void * const userdata) { - DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; // Error control @@ -212,19 +87,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 - || CCTK_GroupTypeFromVarI(var[n]) == CCTK_ARRAY); - assert (CCTK_GroupDimFromVarI(var[n]) <= dim); + assert (CCTK_GroupTypeFromVarI(var[n]) == CCTK_GF); + 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 - || CCTK_GroupTypeFromVarI(res[n]) == CCTK_ARRAY); - assert (CCTK_GroupDimFromVarI(res[n]) <= dim); + assert (CCTK_GroupTypeFromVarI(res[n]) == CCTK_GF); + 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]); @@ -236,47 +111,7 @@ 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); @@ -310,7 +145,8 @@ 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); @@ -324,197 +160,375 @@ 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; - if (reflevel!=-1) { - 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]; - // 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; + assert (is_global_mode() || is_level_mode()); + + // The level to solve on + const int solve_level = is_level_mode() ? reflevel : reflevels-1; - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, "Solving through adaptive Jacobi iterations"); + 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; } - int iter; - for (iter=0; iter<=maxiters; ++iter) { +#warning "TODO: assert that all levels are at the same time" + + // Switch to global mode + BEGIN_GLOBAL_MODE(cctkGH) { - // 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."); + // 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)); } - Util_TableSetInt (options_table, iter, "iters"); - Util_TableDeleteKey (options_table, "error"); - goto done; - } + } END_REFLEVEL_LOOP; - // 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]); - // Calculate convergence speed - speed = old_norm2 / norm2; + // 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; - // Log output if (verbose || veryverbose) { - const time_t currenttime = time(0); - if (iter == maxiters -#if HAVE_ISNAN - || 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)"); + 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); } - if (did_hit_max) { - CCTK_VInfo (CCTK_THORNSTRING, - "Convergence factor became too large: artificially kept at maxfactor (may lead to slower convergence)"); + * 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); } - did_hit_min = false; - did_hit_max = false; + 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; + } - if (veryverbose) { + } END_REFLEVEL_LOOP; + + if (ierr) { + if (verbose || 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 { + "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, - "Iteration %d, time %g, residual %g", - iter, (double)(currenttime-starttime), (double)norm2); + "Residual calculation or boundary conditions reported error; aborting solver."); } - if (outeveryseconds > 0) { - while (nexttime <= currenttime) nexttime += outeveryseconds; + 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) +#if HAVE_ISNAN + || 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 { + CCTK_VInfo (CCTK_THORNSTRING, + "Iteration %d, time %g, residual %g", + iter, double(currenttime-starttime), double(norm2)); + } + if (outeveryseconds > 0) { + while (nexttime <= currenttime) nexttime += outeveryseconds; + } } } - } - - 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."); + // Exit if things go bad + if (isnan(norm2)) { + if (verbose || veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "Encountered NaN. Aborting solver."); + } + break; } - break; - } #endif - - // Return when desired accuracy is reached - if (norm2 <= minerror) { - if (verbose || veryverbose) { - CCTK_VInfo (CCTK_THORNSTRING, "Done solving."); + + // 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 = 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; + } + + } // for iter - // 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; + // 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: ; - // Smooth - CallLocalFunction ((cGH *)cctkGH, common::smooth); - // 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."); - } - Util_TableSetInt (options_table, iter, "iters"); - Util_TableDeleteKey (options_table, "error"); - goto done; - } -#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 + // 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); + } + } END_REFLEVEL_LOOP; - } // 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 - if (reflevel!=saved_reflevel) { - set_reflevel ((cGH *)cctkGH, saved_reflevel); - set_mglevel ((cGH *)cctkGH, saved_mglevel); - } + } END_GLOBAL_MODE; return ierr; } |