Add the Blatter stress balance to the stress balance factory

Now we can choose it in pismr and such by using "-stress_balance blatter".
pism · Oct 14, 2020 · 7b691a4 · 7b691a4
1 parent fd9a5c2
commit 7b691a4
Show file tree

Hide file tree

Showing 7 changed files with 41 additions and 18 deletions.
diff --git a/src/pism_config.cdl b/src/pism_config.cdl
@@ -2079,7 +2079,7 @@ netcdf pism_config {
     pism_config:stress_balance.ice_free_thickness_standard_units = "meters";
 
     pism_config:stress_balance.model = "sia";
-    pism_config:stress_balance.model_choices = "none,prescribed_sliding,weertman_sliding,sia,ssa,prescribed_sliding+sia,weertman_sliding+sia,ssa+sia";
+    pism_config:stress_balance.model_choices = "none,prescribed_sliding,weertman_sliding,sia,ssa,prescribed_sliding+sia,weertman_sliding+sia,ssa+sia,blatter";
     pism_config:stress_balance.model_doc = "Stress balance model";
     pism_config:stress_balance.model_option = "stress_balance";
     pism_config:stress_balance.model_type = "keyword";
@@ -2280,6 +2280,11 @@ netcdf pism_config {
     pism_config:stress_balance.blatter.Mz = 5;
     pism_config:stress_balance.blatter.Mz_doc = "; Number of vertical grid levels in the ice (in the Blatter solver).";
 
+    pism_config:stress_balance.blatter.n_levels_units = "count";
+    pism_config:stress_balance.blatter.n_levels_type = "integer";
+    pism_config:stress_balance.blatter.n_levels = 1;
+    pism_config:stress_balance.blatter.n_levels_doc = "Number of multigrid levels to pad the grid for";
+
     pism_config:stress_balance.ssa.strength_extension.constant_nu = 9.48680701906572e+14;
     pism_config:stress_balance.ssa.strength_extension.constant_nu_doc = "The SSA is made elliptic by use of a constant value for the product of viscosity (nu) and thickness (H).  This value for nu comes from hardness (bar B)=1.9e8 `Pa s^{1/3}` :cite:`MacAyealetal` and a typical strain rate of 0.001 year-1:  `\\nu = (\\bar B) / (2 \\cdot 0.001^{2/3})`.  Compare the value of 9.45e14 Pa s = 30 MPa year in :cite:`Ritzetal2001`.";
     pism_config:stress_balance.ssa.strength_extension.constant_nu_type = "number";

diff --git a/src/stressbalance/blatter/Blatter.cc b/src/stressbalance/blatter/Blatter.cc
@@ -725,22 +725,21 @@ PetscErrorCode Blatter::jacobian_callback(DMDALocalInfo *info,
  * @param[in] grid PISM's grid.
  * @param[in] Mz number of vertical levels
  * @param[in] n_levels maximum number of grid levels to use
+ * @param[in] coarsening_factor grid coarsening factor
  */
-Blatter::Blatter(IceGrid::ConstPtr grid, int Mz, int n_levels)
+Blatter::Blatter(IceGrid::ConstPtr grid, int Mz, int n_levels, int coarsening_factor)
   : ShallowStressBalance(grid) {
 
   auto pism_da = grid->get_dm(1, 0);
 
-  int coarsening_factor = m_config->get_number("stress_balance.blatter.coarsening_factor");
-
   int ierr = setup(*pism_da, Mz, n_levels, coarsening_factor);
   if (ierr != 0) {
     throw RuntimeError(PISM_ERROR_LOCATION,
                        "Failed to allocate a Blatter solver instance");
   }
 
   {
-    int mz = Mz + grid_padding(Mz, n_levels);
+    int mz = Mz + grid_padding(Mz, coarsening_factor, n_levels);
     std::vector<double> sigma(mz);
     double dz = 1.0 / (mz - 1);
     for (int i = 0; i < mz; ++i) {
@@ -867,7 +866,7 @@ PetscErrorCode Blatter::setup(DM pism_da, int Mz, int n_levels, int coarsening_f
     {
       // x direction
       {
-        int pad_x = grid_padding(Mx, n_levels);
+        int pad_x = grid_padding(Mx, coarsening_factor, n_levels);
 
         new_lx[Nx - 1] += pad_x;
         Mx             += pad_x;
@@ -876,15 +875,15 @@ PetscErrorCode Blatter::setup(DM pism_da, int Mz, int n_levels, int coarsening_f
 
       // y direction
       {
-        int pad_y = grid_padding(My, n_levels);
+        int pad_y = grid_padding(My, coarsening_factor, n_levels);
 
         new_ly[Ny - 1] += pad_y;
         My             += pad_y;
         y_max          += pad_y * m_grid->dy();
       }
 
       // z direction
-      Mz += grid_padding(Mz, n_levels);
+      Mz += grid_padding(Mz, coarsening_factor, n_levels);
     }
 
     ierr = DMDACreate3d(PETSC_COMM_WORLD,

diff --git a/src/stressbalance/blatter/Blatter.hh b/src/stressbalance/blatter/Blatter.hh
@@ -32,7 +32,7 @@ namespace stressbalance {
 
 class Blatter : public ShallowStressBalance {
 public:
-  Blatter(IceGrid::ConstPtr grid, int Mz, int n_levels);
+  Blatter(IceGrid::ConstPtr grid, int Mz, int n_levels, int coarsening_factor);
   virtual ~Blatter();
 
   void update(const Inputs &inputs, bool);

diff --git a/src/stressbalance/blatter/Poisson3.cc b/src/stressbalance/blatter/Poisson3.cc
@@ -479,8 +479,10 @@ Poisson3::Poisson3(IceGrid::ConstPtr grid, int Mz, int n_levels)
                        "Failed to allocate a Poisson3 instance");
   }
 
+  int coarsening_factor = 2;
+
   {
-    int mz = Mz + grid_padding(Mz, n_levels);
+    int mz = Mz + grid_padding(Mz, coarsening_factor, n_levels);
     std::vector<double> sigma(mz);
     double dz = 1.0 / (mz - 1);
     for (int i = 0; i < mz; ++i) {
@@ -589,6 +591,8 @@ PetscErrorCode Poisson3::setup(DM pism_da, int Mz, int n_levels) {
       }
     }
 
+    int coarsening_factor = 2;
+
     double
       x_max = m_grid->Lx(),
       x_min = -x_max,
@@ -599,7 +603,7 @@ PetscErrorCode Poisson3::setup(DM pism_da, int Mz, int n_levels) {
     {
       // x direction
       {
-        int pad_x = grid_padding(Mx, n_levels);
+        int pad_x = grid_padding(Mx, coarsening_factor, n_levels);
 
         new_lx[Nx - 1] += pad_x;
         Mx             += pad_x;
@@ -608,15 +612,15 @@ PetscErrorCode Poisson3::setup(DM pism_da, int Mz, int n_levels) {
 
       // y direction
       {
-        int pad_y = grid_padding(My, n_levels);
+        int pad_y = grid_padding(My, coarsening_factor, n_levels);
 
         new_ly[Ny - 1] += pad_y;
         My             += pad_y;
         y_max          += pad_y * m_grid->dy();
       }
 
       // z direction
-      Mz += grid_padding(Mz, n_levels);
+      Mz += grid_padding(Mz, coarsening_factor, n_levels);
     }
 
     ierr = DMDACreate3d(PETSC_COMM_WORLD,

diff --git a/src/stressbalance/blatter/grid_hierarchy.cc b/src/stressbalance/blatter/grid_hierarchy.cc
@@ -51,17 +51,18 @@ double grid_z(double b, double H, int Mz, int k) {
  * Compute the padding needed to allow for `n_levels` of coarsening.
  *
  * @param[in] N number of grid points (nodes)
+ * @param[in] factor coarsening factor
  * @param[in] n_levels number of coarsening levels
  *
  * @return padding amount
  */
-int grid_padding(int N, int n_levels) {
+int grid_padding(int N, int factor, int n_levels) {
   // number of spaces
   int k = N - 1;
   int C = 1;
   for (int n = 0; n < n_levels; ++n) {
-    C *= 2;
-    k = (k % 2 ? k + 1: k) / 2;
+    C *= factor;
+    k = (k % factor ? k + (factor - (k % factor)): k) / factor;
   }
   return (C * k + 1) - N;
 }

diff --git a/src/stressbalance/blatter/grid_hierarchy.hh b/src/stressbalance/blatter/grid_hierarchy.hh
@@ -56,7 +56,7 @@ double grid_xy(double min, double delta, int k);
 
 double grid_z(double b, double H, int Mz, int k);
 
-int grid_padding(int N, int n_levels);
+int grid_padding(int N, int coarsening_factor, int n_levels);
 
 DMDALocalInfo grid_transpose(const DMDALocalInfo &input);
 

diff --git a/src/stressbalance/factory.cc b/src/stressbalance/factory.cc
@@ -25,6 +25,9 @@
 #include "SSB_Modifier.hh"
 #include "pism/regional/SSAFD_Regional.hh"
 #include "pism/regional/SIAFD_Regional.hh"
+#include "blatter/Blatter.hh"
+#include "blatter/BlatterMod.hh"
+
 #include "pism/util/pism_utilities.hh"
 #include "pism/util/Context.hh"
 #include "pism/stressbalance/ssa/SSAFEM.hh"
@@ -35,17 +38,28 @@ namespace stressbalance {
 std::shared_ptr<StressBalance> create(const std::string &model,
                                       IceGrid::ConstPtr grid,
                                       bool regional) {
-  std::shared_ptr<ShallowStressBalance> sliding;
 
   auto config = grid->ctx()->config();
 
+  if (model == "blatter") {
+    int Mz = config->get_number("stress_balance.blatter.Mz");
+    int n_levels = config->get_number("stress_balance.blatter.n_levels");
+    int coarsening_factor = config->get_number("stress_balance.blatter.coarsening_factor");
+
+    std::shared_ptr<Blatter> blatter(new Blatter(grid, Mz, n_levels, coarsening_factor));
+    std::shared_ptr<BlatterMod> mod(new BlatterMod(blatter));
+
+    return std::shared_ptr<StressBalance>(new StressBalance(grid, blatter, mod));
+  }
+
   SSAFactory SSA;
   if (config->get_string("stress_balance.ssa.method") == "fd") {
     SSA = regional ? SSAFD_RegionalFactory : SSAFDFactory;
   } else {
     SSA = SSAFEMFactory;
   }
 
+  std::shared_ptr<ShallowStressBalance> sliding;
   if (member(model, {"none", "sia"})) {
     sliding.reset(new ZeroSliding(grid));
   } else if (member(model, {"prescribed_sliding", "prescribed_sliding+sia"})) {