CarpetMG: Snapshot of current state

darcs-hash:20041208222153-891bb-d9ac7d30cfb9ed9eda22377bcb72317b349cbc73.gz

6 files changed, 172 insertions, 243 deletions

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<CCTK_REAL> 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<CCTK_INT> 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<options.factors.size(); ++n) {
         options.factors.at(n) = 1.0;
       }
-      icnt = (Util_TableGetRealArray
-              (options_table, options.factors.size(),
-               &options.factors.front(), "factors"));
+      int const 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");
+      options.mgiters = 10;
+      int const 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");
+      int const 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");
+      int const icnt
+        = Util_TableGetInt (options_table, &options.poststeps, "poststeps");
       assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY);
     }
     
@@ -471,34 +487,11 @@ namespace CarpetMG {
              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,
@@ -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<cGH *>(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<cGH *>(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<CCTK_INT> 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<dim; ++d) {
       // TODO: correct this for solving on grid arrays instead of grid
@@ -845,10 +827,11 @@ namespace CarpetMG {
                ? (static_cast<CCTK_REAL const *>
                   (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);