Commit work on the second verification test

src/stressbalance/blatter/CMakeLists.txt

@@ -5,6 +5,7 @@ add_library (blatter OBJECT
   BlatterMod.cc
   util/grid_hierarchy.cc
   verification/BlatterTest1.cc
+  verification/BlatterTest2.cc
   verification/manufactured_solutions.cc
   ismip-hom/BlatterISMIPHOM.cc
   )

src/stressbalance/blatter/verification/BlatterTest2.cc

@@ -0,0 +1,111 @@
+/* Copyright (C) 2020 PISM Authors
+ *
+ * This file is part of PISM.
+ *
+ * PISM is free software; you can redistribute it and/or modify it under the
+ * terms of the GNU General Public License as published by the Free Software
+ * Foundation; either version 3 of the License, or (at your option) any later
+ * version.
+ *
+ * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PISM; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+#include "BlatterTest2.hh"
+
+#include "pism/rheology/IsothermalGlen.hh"
+
+#include "manufactured_solutions.hh"
+
+namespace pism {
+namespace stressbalance {
+
+BlatterTest2::BlatterTest2(IceGrid::ConstPtr grid, int Mz, int n_levels, int coarsening_factor)
+  : Blatter(grid, Mz, n_levels, coarsening_factor) {
+
+  // use the isothermal Glen flow law
+  m_flow_law.reset(new rheology::IsothermalGlen("stress_balance.blatter.", *m_config, m_EC));
+
+  // make sure we use the same Glen exponent
+  assert(m_flow_law->exponent() == 3.0);
+
+  // make sure the grid is periodic in the Y direction
+  assert(m_grid->periodicity() == Y_PERIODIC);
+
+  // store constant hardness (enthalpy and pressure values are irrelevant)
+  m_B = m_flow_law->hardness(1e5, 0);
+
+  // surface elevation parameter (1/m)
+  m_alpha = 4e-8;
+
+  // surface elevation parameter (m)
+  m_s0 = 2000.0;
+
+  // ice density times g
+  m_rhog = (m_config->get_number("constants.ice.density") *
+            m_config->get_number("constants.standard_gravity"));
+
+  // ice thickness (m)
+  m_H0 = 1000.0;
+}
+
+bool BlatterTest2::neumann_bc_face(int face, const int *node_type) {
+  (void) face;
+  (void) node_type;
+
+  return false;
+}
+
+bool BlatterTest2::dirichlet_node(const DMDALocalInfo &info, const fem::Element3::GlobalIndex& I) {
+  // use Dirichlet BC on the whole map-plane boundary
+  return (I.i == 0 or I.i == info.mx - 1);
+}
+
+Vector2 BlatterTest2::u_bc(double x, double y, double z) {
+  (void) y;
+
+  return exact_xz(x, z, m_B, m_rhog, m_s0, m_alpha, m_H0);
+}
+
+void BlatterTest2::residual_source_term(const fem::Q1Element3 &element,
+                                        const double *surface,
+                                        Vector2 *residual) {
+  (void) surface;
+
+  // compute x and y coordinates of quadrature points
+  const int Nq = element.n_pts();
+  std::vector<double> x(Nq), z(Nq);
+  {
+    const int Nk = fem::q13d::n_chi;
+    double x_nodal[Nk], z_nodal[Nk];
+
+    for (int n = 0; n < Nk; ++n) {
+      x_nodal[n] = element.x(n);
+      z_nodal[n] = element.z(n);
+    }
+
+    element.evaluate(x_nodal, x.data());
+    element.evaluate(z_nodal, z.data());
+  }
+
+  // loop over all quadrature points
+  for (int q = 0; q < Nq; ++q) {
+    auto W = element.weight(q);
+
+    // loop over all test functions
+    for (int t = 0; t < element.n_chi(); ++t) {
+      const auto &psi = element.chi(q, t);
+
+      residual[t] += W * psi.val * source_xz(x[q], z[q], m_B, m_rhog, m_s0, m_alpha, m_H0);
+    }
+  }
+}
+
+} // end of namespace stressbalance
+} // end of namespace pism

src/stressbalance/blatter/verification/BlatterTest2.hh

@@ -0,0 +1,75 @@
+/* Copyright (C) 2020 PISM Authors
+ *
+ * This file is part of PISM.
+ *
+ * PISM is free software; you can redistribute it and/or modify it under the
+ * terms of the GNU General Public License as published by the Free Software
+ * Foundation; either version 3 of the License, or (at your option) any later
+ * version.
+ *
+ * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PISM; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+#ifndef BLATTERTEST2_H
+#define BLATTERTEST2_H
+
+#include "pism/stressbalance/blatter/Blatter.hh"
+
+namespace pism {
+namespace stressbalance {
+
+/*!
+ * Implements Dirichlet BC and the source term for a verification test from Tezaur et al,
+ * 2015 (), section 4.2.
+ *
+ * Domain: [-50km, 50km] * [-1, 1] * [b, s].
+ *
+ * Here b is the bed elevation, s is the surface elevation.
+ *
+ * s = b + H0, H0 = 1000m.
+ *
+ * Dirichlet BC are imposed at all the nodes along the lateral boundary.
+ *
+ * Natural boundary conditions are used on the top boundary.
+ *
+ * The basal boundary condition includes the sliding condition and a correction for the
+ * chosen manufactured solution.
+ *
+ */
+class BlatterTest2 : public Blatter {
+public:
+  BlatterTest2(IceGrid::ConstPtr grid, int Mz, int n_levels, int coarsening_factor);
+
+protected:
+  bool neumann_bc_face(int face, const int *node_type);
+
+  bool dirichlet_node(const DMDALocalInfo &info, const fem::Element3::GlobalIndex& I);
+
+  Vector2 u_bc(double x, double y, double z);
+
+  void residual_source_term(const fem::Q1Element3 &element,
+                            const double *surface,
+                            Vector2 *residual);
+  //! constant ice hardness
+  double m_B;
+
+  double m_alpha;
+
+  double m_s0;
+
+  double m_H0;
+
+  double m_rhog;
+};
+
+} // end of namespace stressbalance
+} // end of namespace pism
+
+#endif /* BLATTERTEST2_H */