From 08dcb39eb8e3d0bc92e222159d153c2552dceae0 Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Wed, 8 Dec 2004 22:21:00 +0000 Subject: CarpetMG: Snapshot of current state darcs-hash:20041208222153-891bb-d9ac7d30cfb9ed9eda22377bcb72317b349cbc73.gz --- CarpetDev/CarpetMG/interface.ccl | 9 ++ CarpetDev/CarpetMG/par/charge-ml2.par | 15 ++- CarpetDev/CarpetMG/par/charge-ml4.par | 13 ++- CarpetDev/CarpetMG/par/charge.par | 121 -------------------- CarpetDev/CarpetMG/src/README | 52 +++++++++ CarpetDev/CarpetMG/src/mg.cc | 205 ++++++++++++++++------------------ 6 files changed, 172 insertions(+), 243 deletions(-) delete mode 100644 CarpetDev/CarpetMG/par/charge.par (limited to 'CarpetDev/CarpetMG') diff --git a/CarpetDev/CarpetMG/interface.ccl b/CarpetDev/CarpetMG/interface.ccl index 1023297ee..8d757bc39 100644 --- a/CarpetDev/CarpetMG/interface.ccl +++ b/CarpetDev/CarpetMG/interface.ccl @@ -12,3 +12,12 @@ USES INCLUDE HEADER: carpet.hh CCTK_INT FUNCTION \ SymmetryTableHandleForGrid (CCTK_POINTER_TO_CONST IN cctkGH) REQUIRES FUNCTION SymmetryTableHandleForGrid + +CCTK_INT FUNCTION \ + Boundary_SelectVarForBCI (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN faces, \ + CCTK_INT IN boundary_width, \ + CCTK_INT IN table_handle, \ + CCTK_INT IN var_index, \ + CCTK_STRING IN bc_name) +REQUIRES FUNCTION Boundary_SelectVarForBCI diff --git a/CarpetDev/CarpetMG/par/charge-ml2.par b/CarpetDev/CarpetMG/par/charge-ml2.par index 5ea07bc51..325e1f7df 100644 --- a/CarpetDev/CarpetMG/par/charge-ml2.par +++ b/CarpetDev/CarpetMG/par/charge-ml2.par @@ -37,7 +37,7 @@ Carpet::init_3_timelevels = yes ActiveThorns = "TATelliptic CarpetMG" CarpetMG::verbose = yes -#CarpetMG::veryverbose = yes +CarpetMG::veryverbose = yes @@ -93,7 +93,7 @@ CarpetRegrid::coordinates = " ActiveThorns = "WaveToyC" -WaveToyC::bound = radiation +WaveToyC::bound = radiation @@ -101,10 +101,13 @@ ActiveThorns = "IDScalarWave IDSWTEcarpet" IDScalarWave::initial_data = charge-TATelliptic-carpet -IDSWTEcarpet::solver = CarpetMG -IDSWTEcarpet::options = "" -IDSWTEcarpet::radius = 5.5 -IDSWTEcarpet::charge = 1.0 +IDSWTEcarpet::solver = CarpetMG +#IDSWTEcarpet::options = "presteps=4 poststeps=4" +IDSWTEcarpet::options = "presteps=4 poststeps=4 mgiters=10 maxiters=1000" +IDSWTEcarpet::shape = Gauss +IDSWTEcarpet::radius = 5.5 +IDSWTEcarpet::charge = 1.0 +IDSWTEcarpet::boundary = scalar diff --git a/CarpetDev/CarpetMG/par/charge-ml4.par b/CarpetDev/CarpetMG/par/charge-ml4.par index 78767862b..04804f0a6 100644 --- a/CarpetDev/CarpetMG/par/charge-ml4.par +++ b/CarpetDev/CarpetMG/par/charge-ml4.par @@ -37,7 +37,7 @@ Carpet::init_3_timelevels = yes ActiveThorns = "TATelliptic CarpetMG" CarpetMG::verbose = yes -#CarpetMG::veryverbose = yes +CarpetMG::veryverbose = yes @@ -103,10 +103,13 @@ ActiveThorns = "IDScalarWave IDSWTEcarpet" IDScalarWave::initial_data = charge-TATelliptic-carpet -IDSWTEcarpet::solver = CarpetMG -IDSWTEcarpet::options = "mgiters=10" -IDSWTEcarpet::radius = 5.5 -IDSWTEcarpet::charge = 1.0 +IDSWTEcarpet::solver = CarpetMG +#IDSWTEcarpet::options = "presteps=4 poststeps=4" +IDSWTEcarpet::options = "presteps=4 poststeps=4 mgiters=10 maxiters=1000" +IDSWTEcarpet::shape = Gauss +IDSWTEcarpet::radius = 5.5 +IDSWTEcarpet::charge = 1.0 +IDSWTEcarpet::boundary = scalar diff --git a/CarpetDev/CarpetMG/par/charge.par b/CarpetDev/CarpetMG/par/charge.par deleted file mode 100644 index 675b2ecfb..000000000 --- a/CarpetDev/CarpetMG/par/charge.par +++ /dev/null @@ -1,121 +0,0 @@ -!DESC "Charged sphere initial data, solved with CarpetMG" - -Cactus::cctk_full_warnings = yes -Cactus::cctk_timer_output = full - -Cactus::cctk_itlast = 5120 - - - -ActiveThorns = "NaNCatcher" - - - -ActiveThorns = "IOUtil" - -IO::out_dir = $parfile - - - -ActiveThorns = "Carpet CarpetLib CarpetInterp CarpetReduce CarpetSlab" - -Carpet::domain_from_coordbase = yes -Carpet::max_refinement_levels = 10 - -Carpet::ghost_size = 2 -Carpet::buffer_width = 4 - -Carpet::prolongation_order_space = 3 -Carpet::prolongation_order_time = 2 - -Carpet::convergence_level = 0 - -Carpet::init_3_timelevels = yes - - - -ActiveThorns = "TATelliptic CarpetMG" - -CarpetMG::verbose = yes -#CarpetMG::veryverbose = yes - - - -ActiveThorns = "Boundary CartGrid3D CoordBase ReflectionSymmetry SymBase" - -CoordBase::domainsize = minmax -CoordBase::spacing = numcells - -CoordBase::xmin = 0.0 -CoordBase::ymin = 0.0 -CoordBase::zmin = 0.0 -CoordBase::xmax = 20.0 -CoordBase::ymax = 20.0 -CoordBase::zmax = 20.0 - -CoordBase::ncells_x = 40 -CoordBase::ncells_y = 40 -CoordBase::ncells_z = 40 - -CoordBase::boundary_size_x_lower = 2 -CoordBase::boundary_size_y_lower = 2 -CoordBase::boundary_size_z_lower = 2 - -CoordBase::boundary_shiftout_x_lower = 1 -CoordBase::boundary_shiftout_y_lower = 1 -CoordBase::boundary_shiftout_z_lower = 1 - -CartGrid3D::type = coordbase -CartGrid3D::avoid_origin = no - -ReflectionSymmetry::reflection_x = yes -ReflectionSymmetry::reflection_y = yes -ReflectionSymmetry::reflection_z = yes - -ReflectionSymmetry::avoid_origin_x = no -ReflectionSymmetry::avoid_origin_y = no -ReflectionSymmetry::avoid_origin_z = no - - - -ActiveThorns = "WaveToyC" - -WaveToyC::bound = radiation - - - -ActiveThorns = "IDScalarWave IDSWTEcarpet" - -IDScalarWave::initial_data = charge-TATelliptic-carpet - -IDSWTEcarpet::solver = CarpetMG -IDSWTEcarpet::options = "" -IDSWTEcarpet::radius = 5.5 -IDSWTEcarpet::charge = 1.0 - - - -ActiveThorns = "Time" - -Time::dtfac = 0.5 - - - -ActiveThorns = "IOBasic" - -IOBasic::outInfo_every = 1 -IOBasic::outInfo_vars = "wavetoy::phi" - - - -ActiveThorns = "CarpetIOScalar" - -IOScalar::outScalar_every = 1 -IOScalar::outScalar_vars = "wavetoy::phi" - - - -ActiveThorns = "CarpetIOASCII" - -IOASCII::out1D_every = 1 -IOASCII::out1D_vars = "wavetoy::phi" diff --git a/CarpetDev/CarpetMG/src/README b/CarpetDev/CarpetMG/src/README index dd65e6254..2da9fe325 100644 --- a/CarpetDev/CarpetMG/src/README +++ b/CarpetDev/CarpetMG/src/README @@ -1,3 +1,55 @@ +Old documentation: + +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 + + + +******************************************************************************* + + + +New documentation: + equation: L(varc) = rhsc diff --git a/CarpetDev/CarpetMG/src/mg.cc b/CarpetDev/CarpetMG/src/mg.cc index ab1132758..aba03edd2 100644 --- a/CarpetDev/CarpetMG/src/mg.cc +++ b/CarpetDev/CarpetMG/src/mg.cc @@ -28,6 +28,7 @@ namespace CarpetMG { struct options_t { + // Width of outer boundary CCTK_INT bndwidth; @@ -38,17 +39,19 @@ namespace CarpetMG { // Coarsest level on which the full multigrid algorithm should start CCTK_INT firstlevel; - // Maximum number of direct solver iterations - CCTK_INT maxiters; + // Direct solver + CCTK_INT maxiters; // Maximum number of iterations + CCTK_REAL sor_factor; // SOR factor // Weight factors for the equations CCTK_REAL factor; vector factors; - // Maximum number of multigrid iterations (per level) and smoothing steps - CCTK_INT mgiters; - CCTK_INT presteps; - CCTK_INT poststeps; + // Multigrid iterations + CCTK_INT mgiters; // maximum number of iterations (per level) + CCTK_INT presteps; // number of presmoothing steps + CCTK_INT poststeps; // number of postsmoothing steps + }; @@ -111,7 +114,8 @@ namespace CarpetMG { vector const & wgt, options_t const & options, CCTK_REAL const old_error, - CCTK_REAL & error); + CCTK_REAL & error, + bool const use_sor = false); void norm (cGH const * restrict const cctkGH, @@ -338,19 +342,22 @@ namespace CarpetMG { // 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); + { + int const icnt + = (Util_TableGetIntArray + (symtable, 6, + options.symmetry_handle, "symmetry_handle")); + assert (icnt == 6); + } + { + int const icnt + = (Util_TableGetIntArray + (symtable, 6, + options.symmetry_zone_width, "symmetry_zone_width")); + assert (icnt == 6); + } } // Get full multigrid information @@ -368,40 +375,49 @@ namespace CarpetMG { = Util_TableGetInt (options_table, &options.maxiters, "maxiters"); assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY); } + { + options.sor_factor = 1.0; + int const icnt + = Util_TableGetReal (options_table, &options.sor_factor, "sor_factor"); + 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"); + int const 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= 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, @@ -515,8 +508,8 @@ namespace CarpetMG { if (verbose) { CCTK_VInfo (CCTK_THORNSTRING, - "%*s[%d] beginning multigrid solve", - indwidth(), "", reflevel); + "%*s[%d] beginning multigrid solve, desired residual %g", + indwidth(), "", reflevel, double(minerror)); } // Loop until converged @@ -564,70 +557,40 @@ namespace CarpetMG { // 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); + assert (aux.empty()); int const coarse_reflevel = reflevel - 1; BEGIN_GLOBAL_MODE (cctkGH) { enter_level_mode (const_cast(cctkGH), coarse_reflevel); try { + // Restrict variable restrict_var (cctkGH, var, options); boundary (cctkGH, options_table, applybnds, userdata); copy (cctkGH, sav, var, options); + + // Restrict residual 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); + CCTK_REAL coarse_error; + norm (cctkGH, res, rhs, options, coarse_error); + if (veryverbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "%*s[%d] iteration %d, initial coarse residual %g", + indwidth(), "", reflevel, + iter, double(coarse_error)); + } + // Recurse - CCTK_REAL const coarse_minerror = error / ipow(reflevelfact, dim); + CCTK_REAL const coarse_minerror = error / ipow(reffact, dim); int const ierr = multigrid (cctkGH, var, res, rhs, sav, wgt, aux, @@ -649,11 +612,29 @@ namespace CarpetMG { leave_level_mode (const_cast(cctkGH)); } END_GLOBAL_MODE; +#warning "TODO" + // save old solution + copy (cctkGH, sav, var, options); + zero (cctkGH, res, options); prolongate_var (cctkGH, res, options); add (cctkGH, var, res, options); boundary (cctkGH, options_table, applybnds, userdata); + CCTK_REAL const old_error = error; + + try { + residual (cctkGH, options_table, calcres, userdata); + } catch (char const *) { + assert (0); + } + norm (cctkGH, res, rhs, options, error); + if (error > old_error) { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Residual increased during recursion at level %d from %g to %g", + reflevel, double(old_error), double(error)); + } + if (veryverbose) { CCTK_VInfo (CCTK_THORNSTRING, "%*s[%d] iteration %d, after recursion, residual %g", @@ -733,8 +714,8 @@ namespace CarpetMG { if (verbose) { CCTK_VInfo (CCTK_THORNSTRING, - "%*s[%d] beginning direct solve", - indwidth(), "", reflevel); + "%*s[%d] beginning direct solve, desired residual %g", + indwidth(), "", reflevel, double(minerror)); } assert (is_level_mode()); @@ -759,7 +740,7 @@ namespace CarpetMG { residual (cctkGH, options_table, calcres, userdata); CCTK_REAL const old_error = error; - smooth (cctkGH, var, res, rhs, wgt, options, old_error, error); + smooth (cctkGH, var, res, rhs, wgt, options, old_error, error, true); boundary (cctkGH, options_table, applybnds, userdata); if (veryverbose) { @@ -799,7 +780,8 @@ namespace CarpetMG { vector const & wgt, options_t const & options, CCTK_REAL const old_error, - CCTK_REAL & error) + CCTK_REAL & error, + bool const use_sor) { DECLARE_CCTK_ARGUMENTS; @@ -810,7 +792,7 @@ namespace CarpetMG { // Calculate grid spacings - CCTK_REAL dxinv2 = 0; + CCTK_REAL dxinv2 = 0.0; CCTK_REAL dx[dim]; for (int d=0; d (CCTK_VarDataPtrI (cctkGH, 0, wgt.at(n)))) : NULL); - assert (wgt.size() == 0 or wgtptr); + assert (wgt.empty() or wgtptr); CCTK_REAL const fac - = options.factor * (wgtptr ? 1.0 : options.factors.at(n)); + = ((use_sor ? options.sor_factor : options.factor) + * (wgtptr ? 1.0 : options.factors.at(n) * mdxinv2inv)); int imin[3], imax[3]; interior_shape (cctkGH, options, imin, imax); @@ -859,9 +842,9 @@ namespace CarpetMG { 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; + CCTK_REAL const w = wgtptr ? fac / wgtptr[ind] : fac; - varptr[ind] -= fac * w * diff; + varptr[ind] -= w * diff; ++ count; error2 += pow(diff, 2); -- cgit v1.2.3