diff options
author | Erik Schnetter <schnetter@aei.mpg.de> | 2004-11-24 06:12:00 +0000 |
---|---|---|
committer | Erik Schnetter <schnetter@aei.mpg.de> | 2004-11-24 06:12:00 +0000 |
commit | 7ec2e7dae5b17f1ec9989e8f1b2d024945a4fdc5 (patch) | |
tree | 927882f8b14a03cd05785bb1807d8bab1f665c78 /CarpetDev/CarpetMG/src/mg.cc | |
parent | 66a6c39b2267cdc4c17031fed70b8ba46202f119 (diff) |
CarpetMG: A multigrid solver for Carpet
darcs-hash:20041124061259-891bb-566d93667d53c7746838a83c4171f6874dd496f2.gz
Diffstat (limited to 'CarpetDev/CarpetMG/src/mg.cc')
-rw-r--r-- | CarpetDev/CarpetMG/src/mg.cc | 1302 |
1 files changed, 1302 insertions, 0 deletions
diff --git a/CarpetDev/CarpetMG/src/mg.cc b/CarpetDev/CarpetMG/src/mg.cc new file mode 100644 index 000000000..ab1132758 --- /dev/null +++ b/CarpetDev/CarpetMG/src/mg.cc @@ -0,0 +1,1302 @@ +#include <algorithm> +#include <cassert> +#include <climits> +#include <cmath> +#include <vector> + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include "util_ErrorCodes.h" +#include "util_Table.h" + +#include "carpet.hh" +#include "TATelliptic.h" + + + +namespace CarpetMG { + + using namespace std; + using namespace Carpet; + + + + char const * const solver_name = "CarpetMG"; + + + + struct options_t { + // Width of outer boundary + CCTK_INT bndwidth; + + // Symmetry information + CCTK_INT symmetry_handle[2*dim]; + CCTK_INT symmetry_zone_width[2*dim]; + + // Coarsest level on which the full multigrid algorithm should start + CCTK_INT firstlevel; + + // Maximum number of direct solver iterations + CCTK_INT maxiters; + + // Weight factors for the equations + CCTK_REAL factor; + vector<CCTK_REAL> factors; + + // Maximum number of multigrid iterations (per level) and smoothing steps + CCTK_INT mgiters; + CCTK_INT presteps; + CCTK_INT poststeps; + }; + + + + extern "C" { + + int + CarpetMG_register (); + + void + CarpetMG_paramcheck (CCTK_ARGUMENTS); + + } + + int + solve (cGH const * restrict const cctkGH, + int const * restrict const var, + int const * restrict const res, + int const nvar, + int const options_table, + calcfunc const calcres, + calcfunc const applybnds, + void * const userdata); + int + multigrid (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & var, + vector<CCTK_INT> const & res, + vector<CCTK_INT> const & rhs, + vector<CCTK_INT> const & sav, + vector<CCTK_INT> const & wgt, + vector<CCTK_INT> const & aux, + int const options_table, + calcfunc const calcres, + calcfunc const applybnds, + void * const userdata, + options_t const & options, + CCTK_REAL const minerror); + + int + direct_solve (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & var, + vector<CCTK_INT> const & res, + vector<CCTK_INT> const & rhs, + vector<CCTK_INT> const & wgt, + vector<CCTK_INT> const & aux, + int const options_table, + calcfunc const calcres, + calcfunc const applybnds, + void * const userdata, + options_t const & options, + CCTK_REAL const minerror); + + + + void + smooth (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & var, + vector<CCTK_INT> const & res, + vector<CCTK_INT> const & rhs, + vector<CCTK_INT> const & wgt, + options_t const & options, + CCTK_REAL const old_error, + CCTK_REAL & error); + + void + norm (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & res, + vector<CCTK_INT> const & rhs, + options_t const & options, + CCTK_REAL & error); + + + + void + residual (cGH const * restrict const cctkGH, + int const options_table, + calcfunc const calcres, + void * const userdata) + throw (char const *); + + void + boundary (cGH const * restrict const cctkGH, + int const options_table, + calcfunc const applybnds, + void * const userdata) + throw (char const *); + + + + void + restrict_var (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & var, + options_t const & options); + + void + prolongate_var (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & var, + options_t const & options); + + + + void + zero (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & dst, + options_t const & options); + + void + copy (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & dst, + vector<CCTK_INT> const & src, + options_t const & options); + + void + add (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & dst, + vector<CCTK_INT> const & src, + options_t const & options); + + void + subtract (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & dst, + vector<CCTK_INT> const & src, + options_t const & options); + + + + void + getvars (int const options_table, + char const * restrict const name, + vector<CCTK_INT> & vars); + + void + checkvar (int const vi, + vector<int> & usedvars); + + void + interior_shape (cGH const * restrict const cctkGH, + options_t const & options, + int * restrict const imin, + int * restrict const imax); + + + + int + indwidth (); + + + + // Register this solver + int + CarpetMG_register () + { + int const ierr = TATelliptic_RegisterSolver (solve, solver_name); + assert (! ierr); + + return 0; + } + + + + // Check parameters + void + CarpetMG_paramcheck (CCTK_ARGUMENTS) + { + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + if (CCTK_EQUALS (direct_solver, solver_name)) { + CCTK_VParamWarn (CCTK_THORNSTRING, + "It is not possible to use \"%s\" as direct solver for \"%s\" -- this would lead to an infinite recursion", + solver_name, solver_name); + } + } + + + + // Solve the system + int + solve (cGH const * restrict const cctkGH, + int const * restrict const var_, + int const * restrict const res_, + int const nvar, + int const options_table, + calcfunc const calcres, + calcfunc const applybnds, + void * const userdata) + { + DECLARE_CCTK_PARAMETERS; + + // Check arguments + assert (cctkGH); + assert (var_); + assert (res_); + assert (nvar >= 0); + assert (calcres); + assert (applybnds); + + + + if (verbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "%*s[%d] beginning to solve", + indwidth(), "", reflevel); + } + + + + // Solve the equation on this and some coarser levels, leaving all + // finer levels untouched. + assert (is_level_mode()); + + + + // Copy solution variables from arguments + vector<CCTK_INT> var (nvar); + for (int n=0; n<nvar; ++n) { + var.at(n) = var_[n]; + } + + // Copy residual variables from arguments + vector<CCTK_INT> res (nvar); + for (int n=0; n<nvar; ++n) { + res.at(n) = res_[n]; + } + + // Get RHS variables from options + vector<CCTK_INT> rhs (nvar); + { + int const icnt + = Util_TableGetIntArray (options_table, nvar, &rhs.front(), "rhs"); + assert (icnt == nvar); + } + + // Get save variables from options + vector<CCTK_INT> sav (nvar); + { + int const icnt + = Util_TableGetIntArray (options_table, nvar, &sav.front(), "sav"); + assert (icnt == nvar); + } + + // Get weight variables from options + vector<CCTK_INT> wgt; + getvars (options_table, "wgt", wgt); + int const nwgt = wgt.size(); + assert (nwgt == 0 or nwgt == nvar); + + // Get auxiliary variables from options + vector<CCTK_INT> aux; + getvars (options_table, "aux", aux); + int const naux = aux.size(); + assert (naux >= 0); + + + + // Check variables + { + vector<int> usedvars; + + for (int n=0; n<nvar; ++n) { + checkvar (var.at(n), usedvars); + checkvar (res.at(n), usedvars); + checkvar (rhs.at(n), usedvars); + checkvar (sav.at(n), usedvars); + } + + for (int n=0; n<nwgt; ++n) { + checkvar (wgt.at(n), usedvars); + } + + for (int n=0; n<naux; ++n) { + checkvar (aux.at(n), usedvars); + } + } + + + + options_t options; + + // Get outer boundary information + { + options.bndwidth = 1; + int const icnt + = Util_TableGetInt (options_table, &options.bndwidth, "bndwidth"); + assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY); + } + + // Get symmetry information + { + int icnt; + int const symtable = SymmetryTableHandleForGrid (cctkGH); + assert (symtable >= 0); + + icnt = (Util_TableGetIntArray + (symtable, 6, + options.symmetry_handle, "symmetry_handle")); + assert (icnt == 6); + + icnt = (Util_TableGetIntArray + (symtable, 6, + options.symmetry_zone_width, "symmetry_zone_width")); + assert (icnt == 6); + } + + // Get full multigrid information + { + options.firstlevel = 0; + int const icnt + = Util_TableGetInt (options_table, &options.firstlevel, "firstlevel"); + assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY); + } + + // Get direct solver information + { + options.maxiters = 100; + int const icnt + = Util_TableGetInt (options_table, &options.maxiters, "maxiters"); + assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY); + } + + // Get equation information + { + int icnt; + + options.factor = 0.5; + icnt = Util_TableGetReal (options_table, &options.factor, "factor"); + assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY); + + options.factors.resize (var.size()); + for (int n=0; n<options.factors.size(); ++n) { + options.factors.at(n) = 1.0; + } + icnt = (Util_TableGetRealArray + (options_table, options.factors.size(), + &options.factors.front(), "factors")); + assert (icnt == options.factors.size() + or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY); + } + + // Get multigrid information + { + int icnt; + + options.mgiters = 100; + icnt = Util_TableGetInt (options_table, &options.mgiters, "mgiters"); + assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY); + + options.presteps = 2; + icnt = Util_TableGetInt (options_table, &options.presteps, "presteps"); + assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY); + + options.poststeps = 2; + icnt = Util_TableGetInt (options_table, &options.poststeps, "poststeps"); + assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY); + } + + + + CCTK_REAL minerror = 1.0e-8; + { + int const icnt + = Util_TableGetReal (options_table, &minerror, "minerror"); + assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY); + } + + + + if (reflevel < options.firstlevel) return +1; + + + + // Initialise RHS to zero + BEGIN_GLOBAL_MODE (cctkGH) { + BEGIN_REFLEVEL_LOOP (cctkGH) { + zero (cctkGH, rhs, options); + } END_REFLEVEL_LOOP; + } END_GLOBAL_MODE; + + + + CCTK_REAL error; + int const ierr + = multigrid (cctkGH, + var, res, rhs, sav, wgt, aux, + options_table, + calcres, applybnds, userdata, + options, + minerror); + + + + if (verbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "%*s[%d] finished solving", + indwidth(), "", reflevel); + } + + return ierr; + } + + + + // Perform recursive multigrid cycles until converged + // Return values: + // positive: did not converge, continue + // zero: did converge, exit + // negative: error, abort + int + multigrid (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & var, + vector<CCTK_INT> const & res, + vector<CCTK_INT> const & rhs, + vector<CCTK_INT> const & sav, + vector<CCTK_INT> const & wgt, + vector<CCTK_INT> const & aux, + int const options_table, + calcfunc const calcres, + calcfunc const applybnds, + void * const userdata, + options_t const & options, + CCTK_REAL const minerror) + { +#if 0 + DECLARE_CCTK_ARGUMENTS; +#endif + DECLARE_CCTK_PARAMETERS; + + // Solve on this and some coarser levels + assert (is_level_mode()); + +#if 0 + // Calculate the domain size + int minsize = INT_MAX; + int size = 1; + int maxsize = INT_MAX; + for (d=0; d<dim; ++d) { + int const gsize = cctk_gsh[d]; + assert (gsize >= 0); + minsize = min (minsize, gsize); + + int const lsize = (cctk_lsh[d] + - (cctk_bbox[2*d ] ? 0 : cctk_nghostzones[d]) + - (cctk_bbox[2*d+1] ? 0 : cctk_nghostzones[d])); + assert (lsize >= 0); + assert (lsize <= maxsize); + size *= lsize; + maxsize /= lsize; + } +#endif + + // Solve directly when on the coarsest level + if (reflevel == 0) { + return direct_solve (cctkGH, + var, res, rhs, wgt, aux, + options_table, + calcres, applybnds, userdata, + options, + minerror); + } + + + + try { + + if (verbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "%*s[%d] beginning multigrid solve", + indwidth(), "", reflevel); + } + + // Loop until converged + CCTK_REAL error = HUGE_VAL; + int iter = 0; + while (error > minerror) { + + if (iter >= options.mgiters) { + // Did not converge + if (verbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "%*s[%d] finished multigrid solve after %d iterations, residual %g", + indwidth(), "", reflevel, iter, double(error)); + } + return 1; + } + ++iter; + + + + // Presmooth + { + int step = 0; + while (error > minerror) { + + ++step; + if (step > options.presteps) break; + + residual (cctkGH, options_table, calcres, userdata); + CCTK_REAL const old_error = error; + smooth (cctkGH, var, res, rhs, wgt, options, old_error, error); + boundary (cctkGH, options_table, applybnds, userdata); + + if (veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "%*s[%d] iteration %d, presmoothing step %d, residual %g", + indwidth(), "", reflevel, + iter, step, double(error)); + } + + } // while error > minerror + } + + + + // Restrict + + // solve L(u) = r + // set f := L(u) - r + // + // given: f = L(u ) - r + // wanted: 0 = L(u') - r + // known: R(f) - f is small + // assume: R(u'-u) is also small + // define: fc = R(f) + // uc = R(u) + // fc = L(uc) - rc (defines rc) + // solve: L(uc') = rc + // use: R(u'-u) = uc'-uc (defines u') + // yields: f' = L(u') - r + // where f' is smaller than f + // + // uc <- R(u) + // fc <- R(f) + // rc <- fc + // fc <- L(uc) - rc + // rc <- fc + // then solve L(uc) = rc + // + // the above is identical to: + // uc <- R(u) + // rc <- R(f) + // rc <- L(uc) - rc + // then solve L(uc) = rc + // + // now change the notation: + // varc <- R(var) + // resc-rhsc <- R(res-rhs) + // rhsc <- resc-rhsc + // resc-rhsc <- L(varc) - rhsc + // rhsc <- resc-rhsc + // then solve L(varc) = rhsc + // + // the above is identical to: + // varc <- R(var) + // rhsc <- R(res-rhs) + // rhsc <- L(varc) - rhsc + // then solve L(varc) = rhsc + + residual (cctkGH, options_table, calcres, userdata); + norm (cctkGH, res, rhs, options, error); + subtract (cctkGH, res, rhs, options); + // TODO: restrict and fixup aux + assert (aux.size() == 0); + + int const coarse_reflevel = reflevel - 1; + BEGIN_GLOBAL_MODE (cctkGH) { + enter_level_mode (const_cast<cGH *>(cctkGH), coarse_reflevel); + try { + + restrict_var (cctkGH, var, options); + boundary (cctkGH, options_table, applybnds, userdata); + copy (cctkGH, sav, var, options); + zero (cctkGH, res, options); + restrict_var (cctkGH, res, options); + copy (cctkGH, rhs, res, options); + residual (cctkGH, options_table, calcres, userdata); + subtract (cctkGH, rhs, res, options); + + // Recurse + CCTK_REAL const coarse_minerror = error / ipow(reflevelfact, dim); + int const ierr + = multigrid (cctkGH, + var, res, rhs, sav, wgt, aux, + options_table, + calcres, applybnds, userdata, + options, + coarse_minerror); + if (ierr < 0) throw "multigrid"; + + // Prolongate + copy (cctkGH, res, var, options); + subtract (cctkGH, res, sav, options); + + } catch (char const *) { + // TODO + assert (0); + } + + leave_level_mode (const_cast<cGH *>(cctkGH)); + } END_GLOBAL_MODE; + + zero (cctkGH, res, options); + prolongate_var (cctkGH, res, options); + add (cctkGH, var, res, options); + boundary (cctkGH, options_table, applybnds, userdata); + + if (veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "%*s[%d] iteration %d, after recursion, residual %g", + indwidth(), "", reflevel, + iter, double(error)); + } + + + + // Postsmooth + { + int step = 0; + while (error > minerror) { + + ++step; + if (step > options.poststeps) break; + + residual (cctkGH, options_table, calcres, userdata); + CCTK_REAL const old_error = error; + smooth (cctkGH, var, res, rhs, wgt, options, old_error, error); + boundary (cctkGH, options_table, applybnds, userdata); + + if (veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "%*s[%d] iteration %d, postsmoothing step %d, residual %g", + indwidth(), "", reflevel, + iter, step, double(error)); + } + + } // while error > minerror + } + + + + } // while error > minerror + + if (verbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "%*s[%d] finished multigrid solve after %d iterations, residual %g", + indwidth(), "", reflevel, iter, double(error)); + } + + // Everything went fine + return 0; + + } catch (char const *) { + + // There was an error + return -1; + + } + } + + + + // Solve directly + // This assumes that the grid covers the whole domain + // Return values: + // positive: did not converge, continue + // zero: did converge, exit + // negative: error, abort + int + direct_solve (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & var, + vector<CCTK_INT> const & res, + vector<CCTK_INT> const & rhs, + vector<CCTK_INT> const & wgt, + vector<CCTK_INT> const & aux, + int const options_table, + calcfunc const calcres, + calcfunc const applybnds, + void * const userdata, + options_t const & options, + CCTK_REAL const minerror) + { + DECLARE_CCTK_PARAMETERS; + + if (verbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "%*s[%d] beginning direct solve", + indwidth(), "", reflevel); + } + + assert (is_level_mode()); + + try { + + // Loop until converged + CCTK_REAL error = HUGE_VAL; + int iter = 0; + while (error > minerror) { + + if (iter >= options.maxiters) { + // Did not converge + if (verbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "%*s[%d] finished direct solve after %d iterations, residual is %g", + indwidth(), "", reflevel, iter, double(error)); + } + return 1; + } + ++iter; + + residual (cctkGH, options_table, calcres, userdata); + CCTK_REAL const old_error = error; + smooth (cctkGH, var, res, rhs, wgt, options, old_error, error); + boundary (cctkGH, options_table, applybnds, userdata); + + if (veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "%*s[%d] iteration %d, residual %g", + indwidth(), "", reflevel, + iter, double(error)); + } + + } // while error > min_error + + if (verbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "%*s[%d] finished direct solve after %d iterations, residual is %g", + indwidth(), "", reflevel, iter, double(error)); + } + + // Everything went fine + return 0; + + } catch (void *) { + + // There was an error + return -1; + + } + } + + + + // Smooth + void + smooth (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & var, + vector<CCTK_INT> const & res, + vector<CCTK_INT> const & rhs, + vector<CCTK_INT> const & wgt, + options_t const & options, + CCTK_REAL const old_error, + CCTK_REAL & error) + { + DECLARE_CCTK_ARGUMENTS; + + // Initialise errors + CCTK_REAL count = 0.0; + CCTK_REAL error2 = 0.0; + + + + // Calculate grid spacings + CCTK_REAL dxinv2 = 0; + CCTK_REAL dx[dim]; + for (int d=0; d<dim; ++d) { + // TODO: correct this for solving on grid arrays instead of grid + // functions + dx[d] = CCTK_DELTA_SPACE(d); + dxinv2 += 1.0 / pow(dx[d], 2); + } + CCTK_REAL const mdxinv2inv = 1.0 / (-2.0 * dxinv2); + + + + // Smooth and calculate errors + BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { + BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { + DECLARE_CCTK_ARGUMENTS; + for (int n=0; n<var.size(); ++n) { + + CCTK_REAL * restrict const varptr + = (static_cast<CCTK_REAL *> + (CCTK_VarDataPtrI (cctkGH, 0, var.at(n)))); + assert (varptr); + CCTK_REAL const * restrict const resptr + = (static_cast<CCTK_REAL const *> + (CCTK_VarDataPtrI (cctkGH, 0, res.at(n)))); + assert (resptr); + CCTK_REAL const * restrict const rhsptr + = (static_cast<CCTK_REAL const *> + (CCTK_VarDataPtrI (cctkGH, 0, rhs.at(n)))); + assert (rhsptr); + CCTK_REAL const * restrict const wgtptr + = (wgt.size() > 0 + ? (static_cast<CCTK_REAL const *> + (CCTK_VarDataPtrI (cctkGH, 0, wgt.at(n)))) + : NULL); + assert (wgt.size() == 0 or wgtptr); + + CCTK_REAL const fac + = options.factor * (wgtptr ? 1.0 : options.factors.at(n)); + + int imin[3], imax[3]; + interior_shape (cctkGH, options, imin, imax); + + for (int k=imin[2]; k<imax[2]; ++k) { + for (int j=imin[1]; j<imax[1]; ++j) { + for (int i=imin[0]; i<imax[0]; ++i) { + int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); + + CCTK_REAL const diff = resptr[ind] - rhsptr[ind]; + CCTK_REAL const w = wgtptr ? 1.0 / wgtptr[ind] : mdxinv2inv; + + varptr[ind] -= fac * w * diff; + + ++ count; + error2 += pow(diff, 2); + + } + } + } + } + } END_LOCAL_COMPONENT_LOOP; + } END_MAP_LOOP; + + + + // Reduce errors + int const sum_handle = CCTK_ReductionArrayHandle ("sum"); + assert (sum_handle >= 0); + CCTK_REAL reduce_in[2], reduce_out[2]; + reduce_in[0] = count; + reduce_in[1] = error2; + int const ierr = CCTK_ReduceLocArrayToArray1D + (cctkGH, -1, sum_handle, reduce_in, reduce_out, 2, CCTK_VARIABLE_REAL); + count = reduce_out[0]; + error2 = reduce_out[1]; + if (count > 0) { + error = sqrt(error2 / count); + } else { + error = 0.0; + } + + + + // Sanity check + if (error > old_error) { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Residual increased during smoothing at level %d from %g to %g", + reflevel, double(old_error), double(error)); + } + } + + + + // Calculate the norm of the residual without smoothing + void + norm (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & res, + vector<CCTK_INT> const & rhs, + options_t const & options, + CCTK_REAL & error) + { + DECLARE_CCTK_ARGUMENTS; + + // Initialise errors + CCTK_REAL count = 0.0; + CCTK_REAL error2 = 0.0; + + + + // Calculate errors + BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { + BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { + DECLARE_CCTK_ARGUMENTS; + for (int n=0; n<res.size(); ++n) { + + CCTK_REAL const * restrict const resptr + = (static_cast<CCTK_REAL const *> + (CCTK_VarDataPtrI (cctkGH, 0, res.at(n)))); + assert (resptr); + CCTK_REAL const * restrict const rhsptr + = (static_cast<CCTK_REAL const *> + (CCTK_VarDataPtrI (cctkGH, 0, rhs.at(n)))); + assert (rhsptr); + + int imin[3], imax[3]; + interior_shape (cctkGH, options, imin, imax); + + for (int k=imin[2]; k<imax[2]; ++k) { + for (int j=imin[1]; j<imax[1]; ++j) { + for (int i=imin[0]; i<imax[0]; ++i) { + int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); + + CCTK_REAL const diff = resptr[ind] - rhsptr[ind]; + + ++ count; + error2 += pow(diff, 2); + + } + } + } + } + } END_LOCAL_COMPONENT_LOOP; + } END_MAP_LOOP; + + + + // Reduce errors + int const sum_handle = CCTK_ReductionArrayHandle ("sum"); + assert (sum_handle >= 0); + CCTK_REAL reduce_in[2], reduce_out[2]; + reduce_in[0] = count; + reduce_in[1] = error2; + int const ierr = CCTK_ReduceLocArrayToArray1D + (cctkGH, -1, sum_handle, reduce_in, reduce_out, 2, CCTK_VARIABLE_REAL); + count = reduce_out[0]; + error2 = reduce_out[1]; + if (count > 0) { + error = sqrt(error2 / count); + } else { + error = 0.0; + } + } + + + + // Calculate the residual + void + residual (cGH const * restrict const cctkGH, + int const options_table, + calcfunc const calcres, + void * const userdata) + throw (char const *) + { + int ierr = 0; + 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; + if (ierr) throw "calcres"; + } + + + + // Apply the boundary condition + void + boundary (cGH const * restrict const cctkGH, + int const options_table, + calcfunc const applybnds, + void * const userdata) + throw (char const *) + { + int ierr = 0; + ierr = applybnds (cctkGH, options_table, userdata); + if (ierr) throw "applybnds"; + ierr = CallScheduleGroup (const_cast<cGH *>(cctkGH), "ApplyBCs"); + assert (! ierr); + } + + + + // Restrict the solution variables from the next finer level + void + restrict_var (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & var, + options_t const & options) + { + assert (reflevel < reflevels - 1); + int const tl = 0; + for (comm_state<dim> state; !state.done(); state.step()) { + for (int m=0; m<maps; ++m) { + const CCTK_REAL time = vtt.at(m)->time (tl, reflevel, mglevel); + for (int c=0; c<vhh.at(m)->components(reflevel); ++c) { + for (int n=0; n<var.size(); ++n) { + int const vi = var.at(n); + assert (vi >= 0); + int const gi = CCTK_GroupIndexFromVarI (vi); + assert (gi >= 0); + int const v0 = CCTK_FirstVarIndexI (gi); + assert (v0 >= 0); + arrdata.at(gi).at(m).data.at(vi-v0)->ref_restrict + (state, tl, reflevel, c, mglevel, time); + } // for n + } // for c + } // for m + } // for state + } + + + + // Prolongate the solution variables from the next coarser level + void + prolongate_var (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & var, + options_t const & options) + { + assert (reflevel > 0); + int const tl = 0; + for (comm_state<dim> state; !state.done(); state.step()) { + for (int m=0; m<maps; ++m) { + const CCTK_REAL time = vtt.at(m)->time (tl, reflevel, mglevel); + for (int c=0; c<vhh.at(m)->components(reflevel); ++c) { + for (int n=0; n<var.size(); ++n) { + int const vi = var.at(n); + assert (vi >= 0); + int const gi = CCTK_GroupIndexFromVarI (vi); + assert (gi >= 0); + int const v0 = CCTK_FirstVarIndexI (gi); + assert (v0 >= 0); + arrdata.at(gi).at(m).data.at(vi-v0)->ref_prolongate + (state, tl, reflevel, c, mglevel, time); + } // for n + } // for c + } // for m + } // for state + } + + + + void + zero (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & dst, + options_t const & options) + { + BEGIN_MAP_LOOP (cctkGH, CCTK_GF) { + BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) { + DECLARE_CCTK_ARGUMENTS; + for (int n=0; n<dst.size(); ++n) { + CCTK_REAL * restrict const dstptr + = (static_cast<CCTK_REAL *> + (CCTK_VarDataPtrI (cctkGH, 0, dst.at(n)))); + assert (dstptr); + for (int k=0; k<cctk_lsh[2]; ++k) { + for (int j=0; j<cctk_lsh[1]; ++j) { + for (int i=0; i<cctk_lsh[0]; ++i) { + int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); + dstptr[ind] = 0.0; + } + } + } + } + } END_LOCAL_COMPONENT_LOOP; + } END_MAP_LOOP; + } + + + + void + copy (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & dst, + vector<CCTK_INT> const & src, + options_t const & options) + { + assert (dst.size() == src.size()); + BEGIN_MAP_LOOP (cctkGH, CCTK_GF) { + BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) { + DECLARE_CCTK_ARGUMENTS; + for (int n=0; n<dst.size(); ++n) { + CCTK_REAL * restrict const dstptr + = (static_cast<CCTK_REAL *> + (CCTK_VarDataPtrI (cctkGH, 0, dst.at(n)))); + assert (dstptr); + CCTK_REAL const * restrict const srcptr + = (static_cast<CCTK_REAL *> + (CCTK_VarDataPtrI (cctkGH, 0, src.at(n)))); + assert (srcptr); + for (int k=0; k<cctk_lsh[2]; ++k) { + for (int j=0; j<cctk_lsh[1]; ++j) { + for (int i=0; i<cctk_lsh[0]; ++i) { + int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); + dstptr[ind] = srcptr[ind]; + } + } + } + } + } END_LOCAL_COMPONENT_LOOP; + } END_MAP_LOOP; + } + + + + void + add (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & dst, + vector<CCTK_INT> const & src, + options_t const & options) + { + assert (dst.size() == src.size()); + BEGIN_MAP_LOOP (cctkGH, CCTK_GF) { + BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) { + DECLARE_CCTK_ARGUMENTS; + for (int n=0; n<dst.size(); ++n) { + CCTK_REAL * restrict const dstptr + = (static_cast<CCTK_REAL *> + (CCTK_VarDataPtrI (cctkGH, 0, dst.at(n)))); + assert (dstptr); + CCTK_REAL const * restrict const srcptr + = (static_cast<CCTK_REAL *> + (CCTK_VarDataPtrI (cctkGH, 0, src.at(n)))); + assert (srcptr); + for (int k=0; k<cctk_lsh[2]; ++k) { + for (int j=0; j<cctk_lsh[1]; ++j) { + for (int i=0; i<cctk_lsh[0]; ++i) { + int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); + dstptr[ind] += srcptr[ind]; + } + } + } + } + } END_LOCAL_COMPONENT_LOOP; + } END_MAP_LOOP; + } + + + + void + subtract (cGH const * restrict const cctkGH, + vector<CCTK_INT> const & dst, + vector<CCTK_INT> const & src, + options_t const & options) + { + assert (dst.size() == src.size()); + BEGIN_MAP_LOOP (cctkGH, CCTK_GF) { + BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) { + DECLARE_CCTK_ARGUMENTS; + for (int n=0; n<dst.size(); ++n) { + CCTK_REAL * restrict const dstptr + = (static_cast<CCTK_REAL *> + (CCTK_VarDataPtrI (cctkGH, 0, dst.at(n)))); + assert (dstptr); + CCTK_REAL const * restrict const srcptr + = (static_cast<CCTK_REAL *> + (CCTK_VarDataPtrI (cctkGH, 0, src.at(n)))); + assert (srcptr); + for (int k=0; k<cctk_lsh[2]; ++k) { + for (int j=0; j<cctk_lsh[1]; ++j) { + for (int i=0; i<cctk_lsh[0]; ++i) { + int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); + dstptr[ind] -= srcptr[ind]; + } + } + } + } + } END_LOCAL_COMPONENT_LOOP; + } END_MAP_LOOP; + } + + + + // Get a variable list from an options table + void + getvars (int const options_table, + char const * restrict const name, + vector<CCTK_INT> & vars) + { + int const nvars = Util_TableGetIntArray (options_table, 0, NULL, name); + if (nvars == UTIL_ERROR_TABLE_NO_SUCH_KEY) { + vars.resize (0); + } else if (nvars >= 0) { + assert (nvars >= 0); + vars.resize (nvars); + int const icnt + = Util_TableGetIntArray (options_table, nvars, &vars.front(), name); + assert (icnt == nvars); + } else { + assert (0); + } + } + + + + // Check a grid variable index + void + checkvar (int const vi, + vector<int> & usedvars) + { + int ierr; + + assert (vi>=0 and vi<CCTK_NumVars()); + + int const gi = CCTK_GroupIndexFromVarI (vi); + assert (gi>=0); + + cGroup group; + ierr = CCTK_GroupData (gi, &group); + assert (! ierr); + + assert (group.grouptype == CCTK_GF); + assert (group.vartype == CCTK_VARIABLE_REAL); + assert (group.dim == 3); + + for (int n=0; n<usedvars.size(); ++n) { + assert (vi != usedvars.at(n)); + } + usedvars.push_back (vi); + } + + + + // Calculate the shape of the interior + void + interior_shape (cGH const * restrict const cctkGH, + options_t const & options, + int * restrict const imin, + int * restrict const imax) + { + DECLARE_CCTK_ARGUMENTS; + + for (int d=0; d<dim; ++d) { + + imin[d] = 0; + if (cctk_bbox[2*d]) { + if (options.symmetry_handle[2*d] >= 0) { + imin[d] += options.symmetry_zone_width[2*d]; + } else { + imin[d] += options.bndwidth; + } + } else { + imin[d] += cctk_nghostzones[d]; + } + + imax[d] = cctk_lsh[d]; + if (cctk_bbox[2*d+1]) { + if (options.symmetry_handle[2*d+1] >= 0) { + imax[d] -= options.symmetry_zone_width[2*d+1]; + } else { + imax[d] -= options.bndwidth; + } + } else { + imax[d] -= cctk_nghostzones[d]; + } + + assert (imin[d] <= imax[d]); + + } // for d + } + + + + // Determine indentation + int + indwidth () + { + if (reflevel == -1) return 0; + return 3 * (reflevels - reflevel - 1); + } + +} // namespace CarpetMG |