SSATestCase.cc

Go to the documentation of this file.

// Copyright (C) 2009--2026 Ed Bueler, Constantine Khroulev, and David Maxwell
//
// 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 "pism/stressbalance/ssa/tests/SSATestCase.hh"
#include "pism/stressbalance/ssa/SSAFD_SNES.hh"
#include "pism/stressbalance/StressBalance.hh"
#include "pism/util/Context.hh"
#include "pism/util/Interpolation1D.hh"
#include "pism/util/io/File.hh"
#include "pism/util/pism_options.hh"
#include "pism/util/pism_utilities.hh"
 
#include "pism/stressbalance/ssa/SSAFD.hh"
#include "pism/stressbalance/ssa/SSAFEM.hh"
#include "pism/util/Logger.hh"
#include "pism/util/Time.hh"
#include "pism/util/io/SynchronousOutputWriter.hh"
#include "pism/util/io/IO_Flags.hh"
#include "pism/util/io/io_helpers.hh"
#include <vector>
 
namespace pism {
namespace stressbalance {
 
SSATestCase::SSATestCase(std::shared_ptr<SSA> ssa)
  : m_grid(ssa->grid()),
    m_ctx(m_grid->ctx()),
    m_config(m_ctx->config()),
    m_sys(m_ctx->unit_system()),
    m_tauc(m_grid, "tauc"),
    m_ice_enthalpy(m_grid, "enthalpy", array::WITH_GHOSTS, m_grid->z(), 1),
    m_bc_values(m_grid, "_bc"), // u_bc and v_bc
    m_bc_mask(m_grid, "bc_mask"),
    m_geometry(m_grid),
    m_ssa(ssa) {
 
  m_bc_mask.set_interpolation_type(NEAREST);
 
  // yield stress for basal till (plastic or pseudo-plastic model)
  m_tauc.metadata(0)
      .long_name("yield stress for basal till (plastic or pseudo-plastic model)")
      .units("Pa");
 
  // enthalpy
  m_ice_enthalpy.metadata(0)
      .long_name("ice enthalpy (includes sensible heat, latent heat, pressure)")
      .units("J kg^-1");
 
  // dirichlet boundary condition (FIXME: perhaps unused!)
  m_bc_values.metadata(0)
      .long_name("X-component of the SSA velocity boundary conditions")
      .units("m s^-1")
      .output_units("m year^-1");
  m_bc_values.metadata(1)
      .long_name("Y-component of the SSA velocity boundary conditions")
      .units("m s^-1")
      .output_units("m year^-1");
 
  units::System::Ptr sys  = m_grid->ctx()->unit_system();
  double fill_value =
      units::convert(sys, m_config->get_number("output.fill_value"), "m year^-1", "m second^-1");
 
  auto large_number = units::convert(m_sys, 1e6, "m year^-1", "m second^-1");
 
  m_bc_values.metadata(0)["valid_range"] = { -large_number, large_number };
  m_bc_values.metadata(0)["_FillValue"]  = { fill_value };
 
  m_bc_values.metadata(1)["valid_range"] = { -large_number, large_number };
  m_bc_values.metadata(1)["_FillValue"]  = { fill_value };
 
  m_bc_values.set(fill_value);
 
  // Dirichlet B.C. mask
  m_bc_mask.metadata().long_name("grounded_dragging_floating integer mask");
 
  m_bc_mask.metadata()["flag_values"]   = { 0.0, 1.0 };
  m_bc_mask.metadata()["flag_meanings"] = "no_data ssa.dirichlet_bc_location";
}
 
std::shared_ptr<SSA> SSATestCase::solver(std::shared_ptr<Grid> grid, const std::string &method) {
  if (method == "fem") {
    return std::make_shared<SSAFEM>(grid);
  }
  if (method == "fd") {
    return std::make_shared<SSAFD>(grid, false);
  }
  return  std::make_shared<SSAFD_SNES>(grid, false);
}
 
//! Initialize the test case at the start of a run
void SSATestCase::init() {
 
  m_ssa->init();
 
  // Allow the subclass to set the coefficients.
  initializeSSACoefficients();
}
 
//! Solve the SSA
void SSATestCase::run() {
  m_ctx->log()->message(2, "* Solving the SSA stress balance ...\n");
 
  m_geometry.ensure_consistency(m_config->get_number("stress_balance.ice_free_thickness_standard"));
 
  Inputs inputs;
  inputs.water_column_pressure = nullptr;
  inputs.geometry              = &m_geometry;
  inputs.enthalpy              = &m_ice_enthalpy;
  inputs.basal_yield_stress    = &m_tauc;
  inputs.bc_mask               = &m_bc_mask;
  inputs.bc_values             = &m_bc_values;
 
  bool full_update = true;
  m_ssa->update(inputs, full_update);
}
 
//! Report on the generated solution
void SSATestCase::report(const std::string &testname) {
 
  m_ctx->log()->message(3, m_ssa->stdout_report());
 
  double maxvecerr = 0.0, avvecerr = 0.0, avuerr = 0.0, avverr = 0.0, maxuerr = 0.0, maxverr = 0.0;
  double gmaxvecerr = 0.0, gavvecerr = 0.0, gavuerr = 0.0, gavverr = 0.0, gmaxuerr = 0.0,
         gmaxverr = 0.0;
 
  if (m_config->get_flag("basal_resistance.pseudo_plastic.enabled") &&
      m_config->get_number("basal_resistance.pseudo_plastic.q") != 1.0) {
    m_ctx->log()->message(1, "WARNING: numerical errors not valid for pseudo-plastic till\n");
  }
  m_ctx->log()->message(1, "NUMERICAL ERRORS in velocity relative to exact solution:\n");
 
  const array::Vector &vel_ssa = m_ssa->velocity();
 
  array::AccessScope list{ &vel_ssa };
 
  double exactvelmax = 0, gexactvelmax = 0;
  for (auto p : m_grid->points()) {
    const int i = p.i(), j = p.j();
 
    double uexact, vexact;
    double myx = m_grid->x(i), myy = m_grid->y(j);
 
    exactSolution(i, j, myx, myy, &uexact, &vexact);
 
    double exactnormsq = sqrt(uexact * uexact + vexact * vexact);
    exactvelmax        = std::max(exactnormsq, exactvelmax);
 
    // compute maximum errors
    const double uerr   = fabs(vel_ssa(i, j).u - uexact);
    const double verr   = fabs(vel_ssa(i, j).v - vexact);
    avuerr              = avuerr + uerr;
    avverr              = avverr + verr;
    maxuerr             = std::max(maxuerr, uerr);
    maxverr             = std::max(maxverr, verr);
    const double vecerr = sqrt(uerr * uerr + verr * verr);
    maxvecerr           = std::max(maxvecerr, vecerr);
    avvecerr            = avvecerr + vecerr;
  }
 
  unsigned int N = (m_grid->Mx() * m_grid->My());
 
  gexactvelmax = GlobalMax(m_grid->com, exactvelmax);
  gmaxuerr     = GlobalMax(m_grid->com, maxuerr);
  gmaxverr     = GlobalMax(m_grid->com, maxverr);
  gavuerr      = GlobalSum(m_grid->com, avuerr);
  gavuerr      = gavuerr / N;
  gavverr      = GlobalSum(m_grid->com, avverr);
  gavverr      = gavverr / N;
  gmaxvecerr   = GlobalMax(m_grid->com, maxvecerr);
  gavvecerr    = GlobalSum(m_grid->com, avvecerr);
  gavvecerr    = gavvecerr / N;
 
  using pism::units::convert;
 
  m_ctx->log()->message(
      1, "velocity  :  maxvector   prcntavvec      maxu      maxv       avu       avv\n");
  m_ctx->log()->message(1, "           %11.4f%13.5f%10.4f%10.4f%10.4f%10.4f\n",
                        convert(m_sys, gmaxvecerr, "m second-1", "m year-1"),
                        (gavvecerr / gexactvelmax) * 100.0,
                        convert(m_sys, gmaxuerr, "m second-1", "m year-1"),
                        convert(m_sys, gmaxverr, "m second-1", "m year-1"),
                        convert(m_sys, gavuerr, "m second-1", "m year-1"),
                        convert(m_sys, gavverr, "m second-1", "m year-1"));
 
  m_ctx->log()->message(1, "NUM ERRORS DONE\n");
 
  report_netcdf(testname, convert(m_sys, gmaxvecerr, "m second-1", "m year-1"),
                (gavvecerr / gexactvelmax) * 100.0,
                convert(m_sys, gmaxuerr, "m second-1", "m year-1"),
                convert(m_sys, gmaxverr, "m second-1", "m year-1"),
                convert(m_sys, gavuerr, "m second-1", "m year-1"),
                convert(m_sys, gavverr, "m second-1", "m year-1"));
}
 
namespace details {
static void write(const File &file, const char *name,
                  const char *units,
                  const char *long_name,
                  unsigned int start, double value,
                  io::Type type = io::PISM_DOUBLE) {
  file.define_dimension("N", io::PISM_UNLIMITED);
  file.define_variable(name, type, { "N" });
  file.write_attribute(name, "units", units);
  file.write_attribute(name, "long_name", long_name);
  file.write_variable(name, { start }, { 1 }, &value);
}
} // namespace details
 
void SSATestCase::report_netcdf(const std::string &testname, double max_vector, double rel_vector,
                                double max_u, double max_v, double avg_u, double avg_v) {
  auto sys = m_grid->ctx()->unit_system();
 
  VariableMetadata global_attributes("PISM_GLOBAL", sys);
  global_attributes["source"] = std::string("PISM ") + pism::revision;
 
  options::String filename("-report_file", "NetCDF error report file");
 
  if (not filename.is_set()) {
    return;
  }
 
  m_ctx->log()->message(2, "Also writing errors to '%s'...\n", filename->c_str());
 
  bool append = options::Bool("-append", "Append the NetCDF error report");
 
  io::Mode mode = io::PISM_READWRITE;
  if (not append) {
    mode = io::PISM_READWRITE_MOVE;
  }
 
  // Find the number of records in this file:
  File file(m_grid->com, filename, io::PISM_NETCDF3, mode); // OK to use NetCDF3.
  size_t start = static_cast<size_t>(file.dimension_length("N"));
 
  for (const auto &arg : global_attributes.all_strings()) {
    file.write_attribute("PISM_GLOBAL", arg.first, arg.second);
  }
 
  details::write(file, "N", "", "", start, (double)(start + 1));
 
  details::write(file, "dx", "meters", "", start, m_grid->dx());
 
  details::write(file, "dy", "meters", "", start, m_grid->dy());
 
  details::write(file, "test", "", "", start, (double)testname[0], io::PISM_INT);
 
  details::write(file, "max_velocity", "m year^-1", "maximum ice velocity magnitude error", start,
                 max_vector);
 
  details::write(file, "relative_velocity", "percent", "relative ice velocity magnitude error",
                 start, rel_vector);
 
  details::write(file, "maximum_u", "m year^-1",
                 "maximum error in the X-component of the ice velocity", start, max_u);
 
  details::write(file, "maximum_v", "m year^-1",
                 "maximum error in the Y-component of the ice velocity", start, max_v);
 
  details::write(file, "average_u", "m year^-1",
                 "average error in the X-component of the ice velocity", start, avg_u);
 
  details::write(file, "average_v", "m year^-1",
                 "average error in the Y-component of the ice velocity", start, avg_v);
 
  file.close();
}
 
void SSATestCase::exactSolution(int /*i*/, int /*j*/, double /*x*/, double /*y*/, double *u,
                                double *v) {
  *u = 0;
  *v = 0;
}
 
//! Save the computation and data to a file.
void SSATestCase::write(const std::string &filename) {
  auto writer = std::make_shared<SynchronousOutputWriter>(m_grid->com, *m_config);
  writer->initialize({}, true);
 
  // Write results to an output file:
  OutputFile file(writer, filename);
 
  // Write results to an output file:
  file.define_variable(m_ctx->time()->metadata());
  file.append_time(0.0);
 
  const array::Array *variables[] = {
    &m_geometry.ice_surface_elevation, &m_geometry.ice_thickness, &m_bc_mask,   &m_tauc,
    &m_geometry.bed_elevation,         &m_ice_enthalpy,           &m_bc_values, &m_ssa->velocity()
  };
 
  // define
  for (const auto *v : variables) {
    for (const auto &m : v->all_metadata()) {
      file.define_variable(m);
    }
  }
 
  // write
  for (const auto *v : variables) {
    v->write(file);
  }
 
  // write all easily available diagnostics:
  {
    std::vector<std::shared_ptr<array::Array>> vecs;
    for (auto &pair : m_ssa->spatial_diagnostics()) {
      try {
        vecs.push_back(pair.second->compute());
      } catch (RuntimeError &e) {
        // ignore errors: some diagnostics may need variables not present in this context
      }
    }
 
    // define
    for (auto &v : vecs) {
      for (auto &m : v->all_metadata()) {
        file.define_variable(m);
      }
    }
 
    // write
    for (auto &v : vecs) {
      v->write(file);
    }
  }
 
  array::Vector tmp(m_grid, "_exact");
  tmp.metadata(0)
      .long_name("X-component of the SSA exact solution")
      .units("m s^-1");
  tmp.metadata(1)
      .long_name("Y-component of the SSA exact solution")
      .units("m s^-1");
 
  array::AccessScope list(tmp);
  for (auto p : m_grid->points()) {
    const int i = p.i(), j = p.j();
 
    exactSolution(i, j, m_grid->x(i), m_grid->y(j),
                  &(tmp(i,j).u), &(tmp(i,j).v));
  }
 
  // define
  for (auto &m : tmp.all_metadata()) {
    file.define_variable(m);
  }
  // write
  tmp.write(file);
 
  tmp.metadata(0)
    .set_name("u_error")
    .long_name("X-component of the error (exact - computed)")
    .units("m s^-1");
  tmp.metadata(1)
    .set_name("v_error")
    .long_name("Y-component of the error (exact - computed)")
    .units("m s^-1");
 
  tmp.add(-1.0, m_ssa->velocity());
  // define
  for (auto &m : tmp.all_metadata()) {
    file.define_variable(m);
  }
  tmp.write(file);
 
  array::Scalar error_mag(m_grid, "error_mag");
  error_mag.metadata(0).long_name("magnitude of the error").units("m s^-1");
  array::compute_magnitude(tmp, error_mag);
 
  // define
  for (auto &m : error_mag.all_metadata()) {
    file.define_variable(m);
  }
  // write
  error_mag.write(file);
 
  file.close();
}
 
} // end of namespace stressbalance
} // end of namespace pism