Implement ISMIP-HOM experiments

src/pythonbindings/pism_blatter.i

@@ -5,6 +5,9 @@ pism_class(pism::stressbalance::Blatter,
 pism_class(pism::stressbalance::BlatterTest1,
            "pism/stressbalance/blatter/verification/BlatterTest1.hh")
 
+pism_class(pism::stressbalance::BlatterISMIPHOM,
+           "pism/stressbalance/blatter/ismip-hom/BlatterISMIPHOM.hh")
+
 /* BlatterMod has to be wrapped after Blatter*/
 pism_class(pism::stressbalance::BlatterMod,
            "pism/stressbalance/blatter/BlatterMod.hh")

src/stressbalance/blatter/CMakeLists.txt

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

src/stressbalance/blatter/ismip-hom/BlatterISMIPHOM.cc

@@ -0,0 +1,140 @@
+/* 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 "BlatterISMIPHOM.hh"
+
+#include "pism/util/node_types.hh"
+
+namespace pism {
+namespace stressbalance {
+
+static double A_surface(double x, double y, double L) {
+  (void) y;
+  (void) L;
+  double alpha = 0.5 * (M_PI / 180.0); // 0.5 degrees
+  return -x * tan(alpha);
+}
+
+static double A_bed(double x, double y, double L) {
+  double omega = 2 * M_PI / L;
+  return A_surface(x, y, L) - 1000.0 + 500.0 * sin(omega * x) * sin(omega * y);
+}
+
+static double B_bed(double x, double y, double L) {
+  (void) y;
+  double omega = 2 * M_PI / L;
+  return A_surface(x, y, L) - 1000.0 + 500.0 * sin(omega * x);
+}
+
+static double C_surface(double x, double y, double L) {
+  (void) y;
+  (void) L;
+  double alpha = 0.1 * (M_PI / 180.0); // 0.1 degrees
+  return -x * tan(alpha);
+}
+
+static double C_bed(double x, double y, double L) {
+  return C_surface(x, y, L) - 1000.0;
+}
+
+BlatterISMIPHOM::BlatterISMIPHOM(IceGrid::ConstPtr grid,
+                                 int Mz,
+                                 int n_levels,
+                                 int coarsening_factor,
+                                 ISMIPHOMTest test)
+  : Blatter(grid, Mz, n_levels, coarsening_factor),
+    m_test(test),
+    m_L(2.0 * grid->Lx()) {
+
+  char testname[] = "ABCD";
+  m_log->message(2, "Running ISMIP-HOM Experiment %c (L = %d km)...\n",
+                 testname[m_test],
+                 (int)(m_L * 1e-3));
+
+  switch (m_test) {
+  case HOM_A:
+    m_b = A_bed;
+    m_s = A_surface;
+    break;
+  case HOM_B:
+    m_b = B_bed;
+    // test B surface is the same as for test A
+    m_s = A_surface;
+    break;
+  case HOM_C:
+    m_b = C_bed;
+    m_s = C_surface;
+    break;
+  case HOM_D:
+    // test D geometry is the same as for test C
+    m_b = C_bed;
+    m_s = C_surface;
+    break;
+  default:
+    m_b = nullptr;
+    m_s = nullptr;
+    break;
+  }
+}
+
+void BlatterISMIPHOM::nodal_parameter_values(const fem::Q1Element3 &element,
+                                             Parameters **P,
+                                             int i,
+                                             int j,
+                                             int *node_type,
+                                             double *bottom_elevation,
+                                             double *ice_thickness,
+                                             double *surface_elevation,
+                                             double *sea_level) const {
+  (void) P;
+
+  // This method is called before we get a chance to "reset" the current element, so the
+  // element argument does not yet know about its physical coordinates. This is why we
+  // have to compute x and y "by hand".
+  double
+    x_min = m_grid->x0() - m_grid->Lx(),
+    y_min = m_grid->y0() - m_grid->Ly(),
+    dx    = m_grid->dx(),
+    dy    = m_grid->dy();
+
+  for (int n = 0; n < fem::q13d::n_chi; ++n) {
+    auto I = element.local_to_global(i, j, 0, n);
+
+    double
+      x = x_min + I.i * dx,
+      y = y_min + I.j * dy;
+
+    node_type[n]        = NODE_INTERIOR;
+    bottom_elevation[n] = m_b(x, y, m_L);
+    ice_thickness[n]    = m_s(x, y, m_L) - m_b(x, y, m_L);
+
+    if (surface_elevation) {
+      surface_elevation[n] = m_s(x, y, m_L);
+    }
+
+    if (sea_level) {
+      // set sea level low enough to ensure that all ice is grounded
+      sea_level[n] = bottom_elevation[n] - 1.0;
+    }
+  }
+}
+
+
+} // end of namespace stressbalance
+} // end of namespace pism

src/stressbalance/blatter/ismip-hom/BlatterISMIPHOM.hh

@@ -0,0 +1,63 @@
+/* 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 BLATTERISMIPHOM_H
+#define BLATTERISMIPHOM_H
+
+
+#include "pism/stressbalance/blatter/Blatter.hh"
+
+namespace pism {
+namespace stressbalance {
+
+enum ISMIPHOMTest {HOM_A, HOM_B, HOM_C, HOM_D};
+
+/*!
+ * This class implements periodic geometry experiments from the ISMIP-HOM
+ * inter-comparison.
+ */
+class BlatterISMIPHOM : public Blatter {
+public:
+  BlatterISMIPHOM(IceGrid::ConstPtr grid, int Mz, int n_levels, int coarsening_factor,
+                  ISMIPHOMTest test);
+
+protected:
+  void nodal_parameter_values(const fem::Q1Element3 &element,
+                              Parameters **P,
+                              int i,
+                              int j,
+                              int *node_type,
+                              double *bottom,
+                              double *thickness,
+                              double *surface,
+                              double *sea_level) const;
+  ISMIPHOMTest m_test;
+
+  typedef double (*geometry)(double x, double y, double L);
+
+  geometry m_b, m_s;
+
+  // domain length scale
+  double m_L;
+};
+
+} // end of namespace stressbalance
+} // end of namespace pism
+
+#endif /* BLATTERISMIPHOM_H */

test/ismip-hom/Makefile

@@ -0,0 +1,10 @@
+all: pism ismiphom plot
+
+pism:
+	mpiexec -n 2 python3 run-ismiphom.py -Mx 101 -My 101 -snes_monitor -ksp_monitor -pc_type mg -pc_mg_levels 2 -stress_balance.blatter.coarsening_factor 4 -mg_coarse_pc_type gamg
+
+ismiphom:
+	python3 convert-ismiphom.py
+
+plot:
+	python3 plot-ismiphom.py

test/ismip-hom/README.rst

@@ -0,0 +1,26 @@
+.. default-role:: literal
+
+This directory contains scripts that run ISMIP-HOM_ experiments A-D and plot results.
+
+In short, download and unpack the ISMIP-HOM_ supplement, then run:
+
+.. code::
+
+   python3 run-ismiphom.py
+
+   ln -s path/to/ismip_all .
+
+   python3 convert-ismiphom.py
+   python3 plot-ismiphom.py
+
+The script `run-ismiphom.py` uses PISM's Python bindings to run the Blatter solver in
+PISM. See the top of this script for details.
+
+The script `convert-ismiphom.py` reads submitted model results (see the supplement to
+ISMIP-HOM_), samples them along the line (x, 0.25), and saves to files. To run this
+script, place the `ismip_all` directory (or a symlink to it) in this directory.
+
+The script `plot-ismiphom.py` uses Bokeh_ and data processed by `convert-ismiphom.py`.
+
+.. _Bokeh: https://bokeh.org/_
+.. _ISMIP-HOM: https://tc.copernicus.org/articles/2/95/2008/

test/ismip-hom/convert-ismiphom.py

@@ -0,0 +1,89 @@
+#!/usr/bin/env python3
+
+"""Reads ISMIP-HOM data, interpolates to the same grid, and saves to files using NumPy's
+savez_compressed().
+
+"""
+
+import glob
+import re
+import os
+import numpy as np
+from scipy.interpolate import interp1d, interp2d, griddata
+
+def vx_bd(filename, xs):
+    "Load x-velocity for flowline experiments (B, D) and interpolate onto a given grid."
+    data = np.loadtxt(filename)
+    x = data[:, 0]
+    v = data[:, 1]
+    # Fix some data using the fact that velocities are periodic in x:
+    if np.isnan(v[0]):
+        v[0] = v[-1]
+    return interp1d(x, v, fill_value="extrapolate")(xs)
+
+def vx_ac(filename, xs):
+    "Load x-velocity for 3D experiments (A, C) and interpolate onto a given grid."
+    data = np.loadtxt(filename)
+    x = data[:, 0]
+    y = data[:, 1]
+    v = data[:, 2]
+
+    # Remove points with missing data
+    xx = []
+    yy = []
+    vv = []
+    for X, Y, Z in zip(x, y, v):
+        if not np.isnan(Z):
+            xx.append(X)
+            yy.append(Y)
+            vv.append(Z)
+
+    ys = np.zeros_like(xs) + 0.25
+
+    # method="linear" and "cubic" cause segfaults
+    result = griddata((xx, yy), vv, (xs, ys), method="nearest")
+
+    return result
+
+def sample(files, exp, length_scale, xs, func):
+    "Sample x-velocity along the flow for plotting."
+    result = {}
+
+    for filename in files:
+        pattern = ".*/([a-z0-9]+){exp}{length}(_surf)?\\.txt".format(exp=exp, length=length_scale)
+        m = re.match(pattern, filename)
+
+        if m is not None:
+            model = m.group(1)
+            result[model] = func(filename, xs)
+        else:
+            pass
+
+    return result
+
+if __name__ == "__main__":
+    ismip_prefix = "./ismip_all/"
+
+    files = glob.glob(ismip_prefix + "**/*.txt", recursive=True)
+
+    # 401 is the highest grid resolution in ISMIP-HOM data
+    N_samples=401
+    xs = np.linspace(0, 1, N_samples)
+
+    for ex in "abcd":
+        print("Experiment: ", ex)
+
+        sampler = vx_bd if ex in "bd" else vx_ac
+
+        for length_scale in ["005", "010", "020", "040", "080", "160"]:
+            print("  Length scale: ", length_scale)
+
+            data = sample(files, ex, length_scale, xs, sampler)
+
+            print("Models: ", sorted(list(data.keys())))
+
+            data["x"] = xs
+
+            filename = "ismip-hom-{}-{}".format(ex, length_scale)
+
+            np.savez_compressed(filename, **data)