More work on CFBC verification

src/stressbalance/blatter/residual.cc

@@ -194,7 +194,7 @@ void Blatter::residual_lateral(const fem::Q1Element3 &element,
 
     double
       ice_depth   = s[q] - z[q],
-      p_ice       = m_rho_ice_g *   ice_depth,
+      p_ice       = m_rho_ice_g * ice_depth,
       water_depth = std::max(sea_level[q] - z[q], 0.0),
       p_ocean     = m_rho_ocean_g * water_depth;
 

src/stressbalance/blatter/verification/BlatterTestCFBC.cc

@@ -23,6 +23,7 @@
 #include "pism/stressbalance/StressBalance.hh"
 #include "pism/geometry/Geometry.hh"
 #include "pism/stressbalance/blatter/util/DataAccess.hh"
+#include "manufactured_solutions.hh"
 
 namespace pism {
 namespace stressbalance {
@@ -40,6 +41,8 @@ BlatterTestCFBC::BlatterTestCFBC(IceGrid::ConstPtr grid, int Mz, int n_levels,
   m_rho_i = m_config->get_number("constants.ice.density");
 
   m_rho_w = m_config->get_number("constants.sea_water.density");
+
+  m_L = 2.0 * m_grid->Lx();
 }
 
 bool BlatterTestCFBC::dirichlet_node(const DMDALocalInfo &info,
@@ -49,13 +52,9 @@ bool BlatterTestCFBC::dirichlet_node(const DMDALocalInfo &info,
 }
 
 Vector2 BlatterTestCFBC::u_bc(double x, double y, double z) {
-  (void) x;
   (void) y;
-  (void) z;
-
-  double u = 1.0 / (2.0 * m_B) * x * z * m_g * (m_rho_w - m_rho_i);
 
-  return {u, 0.0};
+  return blatter_xz_cfbc_exact(x, z, m_B, m_L, m_rho_i, m_rho_w, m_g);
 }
 
 void BlatterTestCFBC::residual_source_term(const fem::Q1Element3 &element,
@@ -85,13 +84,13 @@ void BlatterTestCFBC::residual_source_term(const fem::Q1Element3 &element,
   for (int q = 0; q < element.n_pts(); ++q) {
     auto W = element.weight(q);
 
-    double F = m_g * z[q] * (m_rho_w - m_rho_i);
+    auto F = blatter_xz_cfbc_source(x[q], z[q], m_L, m_rho_i, m_rho_w, m_g);
 
     // loop over all test functions
     for (int t = 0; t < element.n_chi(); ++t) {
       const auto &psi = element.chi(q, t);
 
-      residual[t].u += W * psi.val * F;
+      residual[t] += W * psi.val * F;
     }
   }
 }
@@ -114,13 +113,13 @@ void BlatterTestCFBC::residual_surface(const fem::Q1Element3 &element,
   for (int q = 0; q < face.n_pts(); ++q) {
     auto W = face.weight(q);
 
-    double F = (1.0 / 8.0) * m_g * pow(x[q], 2) * (m_rho_w - m_rho_i);
+    auto F = -1.0 * blatter_xz_cfbc_surface(x[q], m_L, m_rho_i, m_rho_w, m_g);
 
     // loop over all test functions
     for (int t = 0; t < element.n_chi(); ++t) {
       auto psi = face.chi(q, t);
 
-      residual[t].u += W * psi * F;
+      residual[t] += W * psi * F;
     }
   }
 }
@@ -151,16 +150,15 @@ void BlatterTestCFBC::residual_basal(const fem::Q1Element3 &element,
   for (int q = 0; q < face.n_pts(); ++q) {
     auto W = face.weight(q);
 
-    double F = - (1.0 / 8.0) * m_g * pow(x[q], 2) * (m_rho_w - m_rho_i);
+    auto F = -1.0 * blatter_xz_cfbc_base(x[q], m_L, m_rho_i, m_rho_w, m_g);
 
     // loop over all test functions
     for (int t = 0; t < element.n_chi(); ++t) {
       auto psi = face.chi(q, t);
 
-      residual[t].u += W * psi * F;
+      residual[t] += W * psi * F;
     }
   }
-
 }
 
 void BlatterTestCFBC::jacobian_basal(const fem::Q1Element3Face &face,

src/stressbalance/blatter/verification/BlatterTestCFBC.hh

@@ -26,7 +26,7 @@ namespace pism {
 namespace stressbalance {
 
 /*!
- * Implements Dirichlet BC for an X-Z flow line setup testing the implementation of lateral (ice margin) boundary conditions.
+ * Implements an X-Z flow line setup testing the implementation of lateral (ice margin) boundary conditions.
  */
 class BlatterTestCFBC : public Blatter {
 public:
@@ -63,6 +63,7 @@ private:
   double m_g;
   double m_rho_i;
   double m_rho_w;
+  double m_L;
 };
 
 } // end of namespace stressbalance

src/stressbalance/blatter/verification/manufactured_solutions.cc

@@ -84,15 +84,15 @@ Vector2 blatter_xz_source_surface(double x, double z, double A, double rho, doub
   };
 }
 
-Vector2 blatter_xz_cfbc_exact(double x, double z, double B, double L, double rho_i, double rho_w)
+Vector2 blatter_xz_cfbc_exact(double x, double z, double B, double L, double rho_i, double rho_w, double g)
 {
   return {
     (1.0/2.0)*L*g*z*(rho_i - rho_w)*sin(M_PI*x/L)/(M_PI*B),
     0
   };
 }
 
-Vector2 blatter_xz_cfbc_source(double x, double z, double B, double L, double rho_i, double rho_w)
+Vector2 blatter_xz_cfbc_source(double x, double z, double L, double rho_i, double rho_w, double g)
 {
   double x0 = M_PI/L;
   return {
@@ -101,15 +101,15 @@ Vector2 blatter_xz_cfbc_source(double x, double z, double B, double L, double rh
   };
 }
 
-Vector2 blatter_xz_cfbc_surface(double x, double z, double B, double L, double rho_i, double rho_w)
+Vector2 blatter_xz_cfbc_surface(double x, double L, double rho_i, double rho_w, double g)
 {
   return {
     (1.0/4.0)*L*g*(rho_i - rho_w)*sin(M_PI*x/L)/M_PI,
     0
   };
 }
 
-Vector2 blatter_xz_cfbc_base(double x, double z, double B, double L, double rho_i, double rho_w)
+Vector2 blatter_xz_cfbc_base(double x, double L, double rho_i, double rho_w, double g)
 {
   return {
     -1.0/4.0*L*g*(rho_i - rho_w)*sin(M_PI*x/L)/M_PI,

src/stressbalance/blatter/verification/manufactured_solutions.hh

@@ -15,12 +15,12 @@ Vector2 blatter_xz_source_bed(double x, double z, double A, double rho, double g
 
 Vector2 blatter_xz_source_surface(double x, double z, double A, double rho, double g, double s_0, double alpha, double H, double beta);
 
-Vector2 blatter_xz_cfbc_exact(double x, double z, double B, double L, double rho_i, double rho_w);
+Vector2 blatter_xz_cfbc_exact(double x, double z, double B, double L, double rho_i, double rho_w, double g);
 
-Vector2 blatter_xz_cfbc_source(double x, double z, double B, double L, double rho_i, double rho_w);
+Vector2 blatter_xz_cfbc_source(double x, double z, double L, double rho_i, double rho_w, double g);
 
-Vector2 blatter_xz_cfbc_surface(double x, double z, double B, double L, double rho_i, double rho_w);
+Vector2 blatter_xz_cfbc_surface(double x, double L, double rho_i, double rho_w, double g);
 
-Vector2 blatter_xz_cfbc_base(double x, double z, double B, double L, double rho_i, double rho_w);
+Vector2 blatter_xz_cfbc_base(double x, double L, double rho_i, double rho_w, double g);
 
 } // end of namespace pism

src/stressbalance/blatter/verification/test_xz_cfbc.py

@@ -37,20 +37,22 @@ def lateral_bc():
     return (eta(u0, v0, n) * M(u0, v0).row(0) * N)[0].subs({nx: 1, ny: 0, nz : 0, x : L})
 
 def print_code(header=False):
-    args = ["x", "z", "B", "L", "rho_i", "rho_w"]
+    args = ["x", "z", "B", "L", "rho_i", "rho_w", "g"]
+    source_args = ["x", "z", "L", "rho_i", "rho_w", "g"]
+    surface_args = ["x", "L", "rho_i", "rho_w", "g"]
     if header:
         declare(name="blatter_xz_cfbc_exact", args=args)
-        declare(name="blatter_xz_cfbc_source", args=args)
-        declare(name="blatter_xz_cfbc_surface", args=args)
-        declare(name="blatter_xz_cfbc_base", args=args)
+        declare(name="blatter_xz_cfbc_source", args=source_args)
+        declare(name="blatter_xz_cfbc_surface", args=surface_args)
+        declare(name="blatter_xz_cfbc_base", args=surface_args)
         return
 
     u0, v0 = exact()
     define(u0, v0, name="blatter_xz_cfbc_exact", args=args)
 
     f_u, f_v = source_term(eta(u0, v0, n), u0, v0)
-    define(f_u, f_v, name="blatter_xz_cfbc_source", args=args)
+    define(f_u, f_v, name="blatter_xz_cfbc_source", args=source_args)
 
     f_s = surface_bc()
-    define(f_s, S(0), name="blatter_xz_cfbc_surface", args=args)
-    define(-f_s, S(0), name="blatter_xz_cfbc_base", args=args)
+    define(f_s, S(0), name="blatter_xz_cfbc_surface", args=surface_args)
+    define(-f_s, S(0), name="blatter_xz_cfbc_base", args=surface_args)

test/blatter_verification.py

@@ -15,6 +15,9 @@
 
 import numpy as np
 
+# Keep petsc4py from suppressing error messages
+PISM.PETSc.Sys.popErrorHandler()
+
 ctx = PISM.Context()
 config = ctx.config
 
@@ -30,9 +33,9 @@
 config.set_number("constants.ice.density", 910.0)
 config.set_number("constants.standard_gravity", 9.81)
 
-def expt(xs, ys):
+def expt(ns, errors):
     "Compute the convergence rate using a polynomial fit."
-    return -np.polyfit(np.log(xs), np.log(ys), 1)[0]
+    return -np.polyfit(np.log(ns), np.log(errors), 1)[0]
 
 class TestXY(TestCase):
     """2D (x-y) verification test using a manufactured solution.
@@ -249,7 +252,8 @@ def setUp(self):
         self.opt = PISM.PETSc.Options()
 
         self.opts = {"-snes_monitor": "",
-                     "-pc_type": "mg"
+                     "-pc_type": "mg",
+                     "-pc_mg_levels": "1", # added here to make tearDown clean it up
                      }
 
         for k, v in self.opts.items():
@@ -453,9 +457,6 @@ def plot(self):
 class TestCFBC(TestCase):
     """Constant viscosity 2D (x-z) verification test checking the implementation of CFBC.
 
-    u = 1/4 * x * z * g * (rho_w - rho_i)
-    v = 0
-
     on a square 0 <= x <= 1, -1 <= z <= 0 with sea level at z = 0 (fully submerged).
 
     Dirichet BC at x = 0, periodic in the Y direction, lateral BC at x = 1. No basal drag.
@@ -465,21 +466,15 @@ def setUp(self):
         "Set PETSc options"
 
         self.H = 1e3
+        self.L = 1.0
 
         self.opt = PISM.PETSc.Options()
 
-        self.opts = {"-snes_monitor": "",
-                     "-ksp_monitor": "",
-                     "-pc_type": "mg"
-                     }
+        self.opts = {"-snes_monitor": ""}
 
         for k, v in self.opts.items():
             self.opt.setValue(k, v)
 
-        # the magnitude of the regularization constant affects the accuracy near the
-        # "dome"
-        config.set_number("flow_law.Schoof_regularizing_velocity", 1e-5)
-
         # constant viscocity
         n = 1.0
         config.set_number("stress_balance.blatter.Glen_exponent", n)
@@ -502,7 +497,7 @@ def inputs(self, N):
         P = PISM.GridParameters(config)
 
         # Domain: [0, 1] * [-dx, dx] * [-1, 0]
-        Lx = 0.5
+        Lx = 0.5 * self.L
         Mx = N
         dx = (2 * Lx) / (Mx - 1)
 
@@ -530,7 +525,8 @@ def inputs(self, N):
         enthalpy.set(1e5)
 
         yield_stress = PISM.IceModelVec2S(grid, "tauc", PISM.WITHOUT_GHOSTS)
-
+        # this value is not important: we use a compensatory term at the base instead of
+        # the sliding law
         yield_stress.set(0.0)
 
         geometry.bed_elevation.set(-self.H)
@@ -551,32 +547,29 @@ def exact_solution(self, grid, bed, Z):
 
         rho_i = config.get_number("constants.ice.density")
         rho_w = config.get_number("constants.sea_water.density")
-        g = config.get_number("constants.standard_gravity")
+        g     = config.get_number("constants.standard_gravity")
 
         u = np.zeros_like(Z)
         with PISM.vec.Access([bed, exact]):
             for (i, j) in grid.points():
                 x = grid.x(i)
                 for k, z_sigma in enumerate(Z):
                     z = bed[i, j] + self.H * z_sigma
-                    u[k] = U = 1 / (4 * self.B) * x**2 * z * g * (rho_w - rho_i)
+                    u[k] = PISM.blatter_xz_cfbc_exact(x, z, self.B, self.L, rho_i, rho_w, g).u
                 exact.set_column(i, j, u)
 
         return exact
 
-    def error_norm(self, N, n_mg):
+    def error_norm(self, N):
         "Return the infinity norm of errors for the u component."
         geometry, enthalpy, yield_stress = self.inputs(N)
 
-        # set the number of multigrid levels
-        self.opt.setValue("-pc_mg_levels", n_mg)
-
         grid = enthalpy.grid()
 
         blatter_Mz = N
-        # do not pad the vertical grid
+        # no geometric multigrid here
         n_levels = 0
-        coarsening_factor = 4
+        coarsening_factor = 1
 
         model = PISM.BlatterTestCFBC(grid, blatter_Mz, n_levels, coarsening_factor)
 
@@ -605,6 +598,7 @@ def error_norm(self, N, n_mg):
 
         # compute the error
         u_error = PISM.IceModelVec3(grid, "error", PISM.WITHOUT_GHOSTS, u_model_z)
+        u_error.set_attrs("exact", "error", "m / s", "m / year", "", 0)
         u_error.copy_from(u_exact)
         u_error.add(-1.0, model.velocity_u_sigma())
 
@@ -617,11 +611,8 @@ def test(self):
 
         # refinement path
         Ns = [11, 21]
-        # number of MG levels to use for a particular grid size, assuming the coarsening
-        # factor of 4
-        mg_levels = [1, 2]
 
-        norms = [self.error_norm(N, n_mg) for (N, n_mg) in zip(Ns, mg_levels)]
+        norms = [self.error_norm(N) for N in Ns]
 
         # Compute the exponent for the convergence rate
         expt_u = expt(Ns, norms)
@@ -633,25 +624,24 @@ def test(self):
 
     def plot(self):
         Ns = [11, 21, 41, 81]
-        mg_levels = [1, 2, 3, 4]
 
         try:
             self.setUp()
-            norms = [self.error_norm(N, n_mg) for (N, n_mg) in zip(Ns, mg_levels)]
+            norms = [self.error_norm(N) for N in Ns]
         finally:
             self.tearDown()
 
-        # the domain is 100km long and 1km high
-        dxs = 100e3 / (np.array(Ns) - 1)
+        dxs = self.L / (np.array(Ns) - 1)
 
         p = np.polyfit(np.log(dxs), np.log(norms), 1)
         fit = np.exp(np.polyval(p, np.log(dxs)))
 
         from bokeh.plotting import figure, show, output_file
 
-        output_file("test-xz.html")
-        f = figure(title = "Blatter-Pattyn stress balance: verification test XZ",
-                   x_axis_type="log", y_axis_type="log", x_axis_label="dx", y_axis_label="max error")
+        output_file("test-xz-cfbc.html")
+        f = figure(title = "Blatter-Pattyn stress balance: verification test XZ-CFBC",
+                   x_axis_type="log", y_axis_type="log", x_axis_label="dx",
+                   y_axis_label="max error")
         f.scatter(dxs, norms)
         f.line(dxs, norms, legend_label="max error", line_width=2)
         f.line(dxs, fit, line_dash="dashed", line_color="red", line_width=2,
@@ -662,10 +652,7 @@ def plot(self):
 
 if __name__ == "__main__":
 
-    # TestXY().plot()
-    # TestXZ().plot()
-
-    t = TestCFBC()
-    t.setUp()
-    print(t.error_norm(41, 2))
-    t.tearDown()
+    for test in [TestXY(),
+                 TestXZ(),
+                 TestCFBC()]:
+        test.plot()