diff options
author | eschnett <> | 2001-03-01 11:40:00 +0000 |
---|---|---|
committer | eschnett <> | 2001-03-01 11:40:00 +0000 |
commit | 310f0ea48d18866b773136aed11200b6eda6378b (patch) | |
tree | 445d3e34ce8b89812994b6614f7bc9f4acbc7fe2 /CarpetDev/CarpetJacobi/src |
Initial revision
darcs-hash:20010301114010-f6438-12fb8a9ffcc80e86c0a97e37b5b0dae0dbc59b79.gz
Diffstat (limited to 'CarpetDev/CarpetJacobi/src')
-rw-r--r-- | CarpetDev/CarpetJacobi/src/Jacobi.cc | 536 | ||||
-rw-r--r-- | CarpetDev/CarpetJacobi/src/make.code.defn | 9 |
2 files changed, 545 insertions, 0 deletions
diff --git a/CarpetDev/CarpetJacobi/src/Jacobi.cc b/CarpetDev/CarpetJacobi/src/Jacobi.cc new file mode 100644 index 000000000..c326b66cd --- /dev/null +++ b/CarpetDev/CarpetJacobi/src/Jacobi.cc @@ -0,0 +1,536 @@ +// $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> +#include <cstdio> +#include <ctime> +#include <iostream> +#include <vector> + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include "util_ErrorCodes.h" +#include "util_Table.h" + +#include "TATelliptic.h" + +#include "carpet.hh" + +#include "th.hh" + + + +namespace Carpet { + // TODO: fix this + void CycleTimeLevels (const cGH* cctkGH); + void Restrict (const cGH* cctkGH); +}; + + + +namespace CarpetJacobi { + + using namespace std; + using namespace Carpet; + + + + extern "C" { + void CarpetJacobi_register (); + } + + int CarpetJacobi_solve (const cGH * const cctkGH, + const int * const var, + const int * const res, + const int nvars, + const int options_table, + const calcfunc calcres, + const calcfunc applybnds, + void * const userdata); + + + + void CarpetJacobi_register () + { + int const ierr = TATelliptic_RegisterSolver + (CarpetJacobi_solve, "CarpetJacobi"); + assert (!ierr); + } + + + + int CarpetJacobi_solve (const cGH * const cctkGH, + const int * const var, + const int * const res, + const int nvars, + const int options_table, + const calcfunc calcres, + const calcfunc applybnds, + void * const userdata) + { + DECLARE_CCTK_PARAMETERS; + + // Error control + int ierr; + + if (! CCTK_IsThornActive(CCTK_THORNSTRING)) { + CCTK_WARN (0, "Thorn " CCTK_THORNSTRING " has not been activated. It is therefore not possible to call CarpetJacobi_solve."); + } + + // Check arguments + assert (cctkGH); + + assert (nvars > 0); + assert (var); + 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_VarTypeI(var[n]) == CCTK_VARIABLE_REAL); + 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_VarTypeI(res[n]) == CCTK_VARIABLE_REAL); + assert (CCTK_QueryGroupStorageI(cctkGH, + CCTK_GroupIndexFromVarI(res[n]))); + } + for (int n=0; n<nvars; ++n) { + assert (var[n] != res[n]); + for (int nn=0; nn<n; ++nn) { + assert (var[nn] != var[n]); + assert (var[nn] != res[n]); + assert (res[nn] != var[n]); + assert (res[nn] != res[n]); + } + } + + + + assert (options_table >= 0); + + assert (calcres); + assert (applybnds); + + // Get options from options table + CCTK_INT maxiters = 10000; + CCTK_REAL minerror = 1e-8; + CCTK_REAL factor = 1.0e-2; // slow start + CCTK_REAL minfactor = 1.0e-8; + CCTK_REAL maxfactor = 1.0; + CCTK_REAL acceleration = 1.2; // slow acceleration + CCTK_REAL deceleration = 0.1; // emergency brake + vector<CCTK_REAL> factors (nvars); + for (int n=0; n<nvars; ++n) { + factors[n] = 1.0; + } + + ierr = Util_TableGetInt (options_table, &maxiters, "maxiters"); + assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY); + ierr = Util_TableGetReal (options_table, &minerror, "minerror"); + assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY); + ierr = Util_TableGetReal (options_table, &factor, "factor"); + assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY); + ierr = Util_TableGetReal (options_table, &minfactor, "minfactor"); + assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY); + ierr = Util_TableGetReal (options_table, &maxfactor, "maxfactor"); + assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY); + ierr = Util_TableGetReal (options_table, &acceleration, "acceleration"); + 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"); + assert (ierr==nvars || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY); + + assert (maxiters >= 0); + assert (minerror >= 0); + assert (factor > 0); + assert (minfactor > 0 && minfactor <= factor); + assert (maxfactor >= factor); + assert (acceleration >= 1); + assert (deceleration > 0 && deceleration <= 1); + for (int n=0; n<nvars; ++n) { + assert (factors[n] != 0); + } + + + + assert (is_global_mode() || is_level_mode()); + + // The level to solve on + const int solve_level = is_level_mode() ? reflevel : reflevels-1; + + 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; + } + +#warning "TODO: assert that all levels are at the same time" + + // Switch to global mode + BEGIN_GLOBAL_MODE(cctkGH) { + + // 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)); + } + } END_REFLEVEL_LOOP; + + + + // 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; + + 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) +#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 HAVE_ISNAN + // Exit if things go bad + if (isnan(norm2)) { + if (verbose || veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "Encountered NaN. Aborting solver."); + } + break; + } +#endif + + // 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; + } + + // 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 + + // 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: ; + + + + // 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; + + } END_GLOBAL_MODE; + + return ierr; + } + +} // namespace CarpetJacobi diff --git a/CarpetDev/CarpetJacobi/src/make.code.defn b/CarpetDev/CarpetJacobi/src/make.code.defn new file mode 100644 index 000000000..efd1ad4c4 --- /dev/null +++ b/CarpetDev/CarpetJacobi/src/make.code.defn @@ -0,0 +1,9 @@ +# Main make.code.defn file for thorn CarpetJacobi +# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/src/make.code.defn,v 1.1 2003/09/02 14:35:58 schnetter Exp $ + +# Source files in this directory +SRCS = Jacobi.cc + +# Subdirectories containing source files +SUBDIRS = + |