Merge branch 'ckhrulev/vec-cleanup' into dev

pism · Oct 6, 2021 · a109d85 · a109d85
2 parents c6d15fc + 1fcf547
commit a109d85
Show file tree

Hide file tree

Showing 87 changed files with 1,199 additions and 1,700 deletions.
diff --git a/examples/inverse/pismi.py b/examples/inverse/pismi.py
@@ -568,8 +568,8 @@ def run():
     r_mag.metadata().set_number("_FillValue", convert(-0.01, 'm/year', 'm/s'))
     r_mag.metadata().set_number("valid_min", 0.0)
 
-    r_mag.set_to_magnitude(residual)
-    r_mag.mask_by(vecs.land_ice_thickness)
+    PISM.compute_magnitude(residual, r_mag)
+    PISM.apply_mask(vecs.land_ice_thickness, 0.0, r_mag)
 
     vecs.add(residual, writing=True)
     vecs.add(r_mag, writing=True)

diff --git a/site-packages/PISM/invert/ssa.py b/site-packages/PISM/invert/ssa.py
@@ -495,13 +495,13 @@ def printTikhonovProgress(invssasolver, it, data):
 
     logMessage("design objective %.8g; weighted %.8g\n" % (designVal, designVal * dWeight))
     if 'grad_JTikhonov' in data:
-        logMessage("gradient: design %.8g state %.8g sum %.8g\n" % (data.grad_JDesign.norm(norm_type) * dWeight,
-                                                                    data.grad_JState.norm(norm_type) * sWeight,
-                                                                    data.grad_JTikhonov.norm(norm_type)))
+        logMessage("gradient: design %.8g state %.8g sum %.8g\n" % (data.grad_JDesign.norm(norm_type)[0] * dWeight,
+                                                                    data.grad_JState.norm(norm_type)[0] * sWeight,
+                                                                    data.grad_JTikhonov.norm(norm_type)[0]))
     else:
-        logMessage("gradient: design %.8g state %.8g; constraints: %.8g\n" % (data.grad_JDesign.norm(norm_type) * dWeight,
-                                                                              data.grad_JState.norm(norm_type) * sWeight,
-                                                                              data.constraints.norm(norm_type)))
+        logMessage("gradient: design %.8g state %.8g; constraints: %.8g\n" % (data.grad_JDesign.norm(norm_type)[0] * dWeight,
+                                                                              data.grad_JState.norm(norm_type)[0] * sWeight,
+                                                                              data.constraints.norm(norm_type)[0]))
     logMessage("tikhonov functional: %.8g\n" % (stateVal * sWeight + designVal * dWeight))
 
 

diff --git a/site-packages/PISM/sia.py b/site-packages/PISM/sia.py
@@ -1,4 +1,4 @@
-# Copyright (C) 2011, 2012, 2014, 2015, 2016, 2017, 2018, 2020 David Maxwell and Constantine Khroulev
+# Copyright (C) 2011, 2012, 2014, 2015, 2016, 2017, 2018, 2020, 2021 David Maxwell and Constantine Khroulev
 #
 # This file is part of PISM.
 #
@@ -62,12 +62,15 @@ def computeSIASurfaceVelocities(modeldata, siasolver=PISM.SIAFD):
     v = sia.velocity_v()
 
     vel_sia = PISM.model.create2dVelocityVec(grid, name="_sia", stencil_width=1)
-    tmp = PISM.IceModelVec2S(grid, 'tmp', PISM.WITHOUT_GHOSTS)
+    u_sia = PISM.IceModelVec2S(grid, 'u_sia', PISM.WITHOUT_GHOSTS)
+    v_sia = PISM.IceModelVec2S(grid, 'v_sia', PISM.WITHOUT_GHOSTS)
 
-    u.extract_surface(md.vecs.thk, tmp)
-    vel_sia.set_component(0, tmp)
+    PISM.extract_surface(u, md.vecs.thk, u_sia)
+    PISM.extract_surface(v, md.vecs.thk, v_sia)
 
-    v.extract_surface(md.vecs.thk, tmp)
-    vel_sia.set_component(1, tmp)
+    with PISM.vec.Access([vel_sia, u_sia, v_sia]):
+        for (i, j) in grid.points():
+            vel_sia[i, j].u = u_sia[i, j]
+            vel_sia[i, j].v = v_sia[i, j]
 
     return vel_sia
diff --git a/site-packages/PISM/ssa.py b/site-packages/PISM/ssa.py
@@ -141,8 +141,9 @@ def write(self, filename, append=False):
         vecs.add(vel_ssa)
 
         velbar_mag = model.createCBarVec(self.grid)
-        velbar_mag.set_to_magnitude(vel_ssa)
-        velbar_mag.mask_by(vecs.thk, util.convert(-0.01, "m/year", "m/second"))
+        PISM.compute_magnitude(vel_ssa, velbar_mag)
+        PISM.apply_mask(vecs.thk, util.convert(-0.01, "m/year", "m/second"),
+                     velbar_mag)
         vecs.add(velbar_mag)
 
         taud = PISM.SSA_taud(self.ssa).compute()

diff --git a/src/age/AgeModel.cc b/src/age/AgeModel.cc
@@ -166,7 +166,7 @@ void AgeModel::update(double t, double dt, const AgeModelInputs &inputs) {
   }
   loop.check();
 
-  m_work.update_ghosts(m_ice_age);
+  m_ice_age.copy_from(m_work);
 }
 
 const IceModelVec3 & AgeModel::age() const {

diff --git a/src/coupler/frontalmelt/DischargeGiven.cc b/src/coupler/frontalmelt/DischargeGiven.cc
@@ -154,7 +154,7 @@ void DischargeGiven::update_impl(const FrontalMeltInputs &inputs, double t, doub
     if (apply(cell_type, i, j) and cell_type.ice_free(i, j)) {
 
       auto R = m_frontal_melt_rate->star(i, j);
-      auto M = cell_type.int_star(i, j);
+      auto M = cell_type.star(i, j);
 
       int N = 0;
       double R_sum = 0.0;

diff --git a/src/coupler/frontalmelt/DischargeRouting.cc b/src/coupler/frontalmelt/DischargeRouting.cc
@@ -141,7 +141,7 @@ void DischargeRouting::update_impl(const FrontalMeltInputs &inputs, double t, do
     if (apply(cell_type, i, j) and cell_type.ice_free(i, j)) {
 
       auto R = m_frontal_melt_rate->star(i, j);
-      auto M = cell_type.int_star(i, j);
+      auto M = cell_type.star(i, j);
 
       int N = 0;
       double R_sum = 0.0;

diff --git a/src/coupler/frontalmelt/FrontalMelt.cc b/src/coupler/frontalmelt/FrontalMelt.cc
@@ -89,7 +89,7 @@ void FrontalMelt::compute_retreat_rate(const Geometry &geometry,
 
         auto H = ice_thickness.star(i, j);
         auto h = surface_elevation.star(i, j);
-        auto M = cell_type.int_star(i, j);
+        auto M = cell_type.star(i, j);
 
         double H_threshold = part_grid_threshold_thickness(M, H, h, bed);
 

diff --git a/src/coupler/ocean/Pico.cc b/src/coupler/ocean/Pico.cc
@@ -145,9 +145,11 @@ void Pico::init_impl(const Geometry &geometry) {
   m_geometry.init(geometry.cell_type);
 
   // FIXME: m_n_basins is a misnomer
-  m_n_basins = static_cast<int>(m_geometry.basin_mask().max()) + 1;
+  m_n_basins = static_cast<int>(max(m_geometry.basin_mask())) + 1;
 
-  m_log->message(4, "PICO basin min=%f, max=%f\n", m_geometry.basin_mask().min(), m_geometry.basin_mask().max());
+  m_log->message(4, "PICO basin min=%f, max=%f\n",
+                 min(m_geometry.basin_mask()),
+                 max(m_geometry.basin_mask()));
 
   PicoPhysics physics(*m_config);
 
@@ -210,7 +212,7 @@ static void extend_basal_melt_rates(const IceModelVec2CellType &cell_type,
 
     const int i = p.i(), j = p.j();
 
-    auto M = cell_type.int_box(i, j);
+    auto M = cell_type.box(i, j);
 
     bool potential_partially_filled_cell =
       ((M.ij == MASK_GROUNDED or M.ij == MASK_ICE_FREE_OCEAN) and
@@ -275,7 +277,7 @@ void Pico::update_impl(const Geometry &geometry, double t, double dt) {
   m_geometry.update(bed_elevation, cell_type);
 
   // FIXME: m_n_shelves is not really the number of shelves.
-  m_n_shelves = static_cast<int>(m_geometry.ice_shelf_mask().max()) + 1;
+  m_n_shelves = static_cast<int>(max(m_geometry.ice_shelf_mask())) + 1;
 
   // Physical part of PICO
   {
@@ -483,7 +485,7 @@ void Pico::set_ocean_input_fields(const PicoPhysics &physics,
 
       // find all basins b, in which the ice shelf s has a calving front with potential ocean water intrusion
       if (mask.as_int(i, j) == MASK_FLOATING) {
-        auto M = mask.int_star(i, j);
+        auto M = mask.star(i, j);
         if (M.n == MASK_ICE_FREE_OCEAN or
             M.e == MASK_ICE_FREE_OCEAN or
             M.s == MASK_ICE_FREE_OCEAN or

diff --git a/src/coupler/ocean/PicoGeometry.cc b/src/coupler/ocean/PicoGeometry.cc
@@ -81,7 +81,7 @@ void PicoGeometry::init(const IceModelVec2CellType &cell_type) {
 
   m_basin_mask.regrid(opt.filename, CRITICAL);
 
-  m_n_basins = static_cast<int>(m_basin_mask.max()) + 1;
+  m_n_basins = static_cast<int>(max(m_basin_mask)) + 1;
 
   m_n_basin_neighbors.resize(2*m_n_basins);
   get_basin_neighbors(cell_type, m_basin_mask, m_n_basin_neighbors);
@@ -121,7 +121,7 @@ void PicoGeometry::update(const IceModelVec2S &bed_elevation,
   // computing ice_shelf_mask and box_mask could be done at the same time
   {
     compute_ice_shelf_mask(m_ice_rises, m_lake_mask, m_ice_shelves);
-    auto n_shelves = static_cast<int>(m_ice_shelves.max()) + 1;
+    auto n_shelves = static_cast<int>(max(m_ice_shelves)) + 1;
 
     std::vector<int> cfs_in_basins_per_shelf(n_shelves*m_n_basins, 0);
     std::vector<int> most_shelf_cells_in_basin(n_shelves, 0);
@@ -524,8 +524,8 @@ void PicoGeometry::get_basin_neighbors(const IceModelVec2CellType &cell_type,
     int b = basin_mask.as_int(i, j);
     if (result[2 * b] == 0 or result[2 * b + 1] == 0) {
 
-      auto B = basin_mask.int_star(i, j);
-      auto M = cell_type.int_star(i, j);
+      auto B = basin_mask.star(i, j);
+      auto M = cell_type.star(i, j);
       int bn = 0; // neighbor basin id
 
       if (M.ij == MASK_ICE_FREE_OCEAN and
@@ -593,7 +593,7 @@ void PicoGeometry::identify_calving_front_connection(const IceModelVec2CellType
       n_shelf_cells_per_basin[sb]++;
 
       if (cell_type.as_int(i, j) == MASK_FLOATING) {
-        auto M = cell_type.int_star(i, j);
+        auto M = cell_type.star(i, j);
         if (M.n == MASK_ICE_FREE_OCEAN or
             M.e == MASK_ICE_FREE_OCEAN or
             M.s == MASK_ICE_FREE_OCEAN or
@@ -774,7 +774,7 @@ void PicoGeometry::compute_distances_cf(const IceModelVec2Int &ocean_mask,
         // if this is an ice shelf cell (or an ice rise) or a hole in an ice shelf
 
         // label the shelf cells adjacent to the ice front with 1
-        auto M = ocean_mask.int_star(i, j);
+        auto M = ocean_mask.star(i, j);
 
         if (M.n == 2 or M.e == 2 or M.s == 2 or M.w == 2) {
           // i.e. there is a neighboring open ocean cell
@@ -830,7 +830,7 @@ void eikonal_equation(IceModelVec2Int &mask) {
 
       if (mask.as_int(i, j) == 0) {
 
-        auto R = mask.int_star(i, j);
+        auto R = mask.star(i, j);
 
         if (R.ij == 0 and
             (R.n == current_label or R.s == current_label or

diff --git a/src/energy/BTU_Full.cc b/src/energy/BTU_Full.cc
@@ -1,4 +1,4 @@
-/* Copyright (C) 2016, 2017, 2018, 2019, 2020 PISM Authors
+/* Copyright (C) 2016, 2017, 2018, 2019, 2020, 2021 PISM Authors
  *
  * This file is part of PISM.
  *
@@ -57,20 +57,24 @@ BTU_Full::BTU_Full(IceGrid::ConstPtr g, const BTUGrid &grid)
     m_Mbz = grid.Mbz;
     m_Lbz = grid.Lbz;
 
-    std::map<std::string, std::string> attrs;
-    attrs["units"]     = "m";
-    attrs["long_name"] = "Z-coordinate in bedrock";
-    attrs["axis"]      = "Z";
-    attrs["positive"]  = "up";
-
     std::vector<double> z(m_Mbz);
     double dz = m_Lbz / (m_Mbz - 1);
     for (unsigned int k = 0; k < m_Mbz; ++k) {
       z[k] = -m_Lbz + k * dz;
     }
     z.back() = 0.0;
 
-    m_temp.reset(new IceModelVec3(m_grid, "litho_temp", "zb", z, attrs));
+    m_temp.reset(new IceModelVec3(m_grid, "litho_temp", WITHOUT_GHOSTS, z));
+    m_temp->metadata(0).z().set_name("zb");
+
+    std::map<std::string, std::string> z_attrs =
+      {{"units", "m"},
+       {"long_name", "Z-coordinate in bedrock"},
+       {"axis", "Z"},
+       {"positive", "up"}};
+    for (const auto &z_attr : z_attrs) {
+      m_temp->metadata(0).z().set_string(z_attr.first, z_attr.second);
+    }
 
     m_temp->set_attrs("model_state",
                       "lithosphere (bedrock) temperature, in BTU_Full",

diff --git a/src/energy/EnergyModel.cc b/src/energy/EnergyModel.cc
@@ -279,7 +279,7 @@ void EnergyModel::update(double t, double dt, const Inputs &inputs) {
     // this call should fill m_work with new values of enthalpy
     this->update_impl(t, dt, inputs);
 
-    m_work.update_ghosts(m_ice_enthalpy);
+    m_ice_enthalpy.copy_from(m_work);
   }
   profiling.end("ice_energy");
 

diff --git a/src/energy/TemperatureModel.cc b/src/energy/TemperatureModel.cc
@@ -353,7 +353,7 @@ void TemperatureModel::update_impl(double t, double dt, const Inputs &inputs) {
   }
 
   // copy to m_ice_temperature, updating ghosts
-  m_work.update_ghosts(m_ice_temperature);
+  m_ice_temperature.copy_from(m_work);
 
   // Set ice enthalpy in place. EnergyModel::update will scatter ghosts
   compute_enthalpy_cold(m_work, ice_thickness, m_work);

diff --git a/src/energy/btutest.cc b/src/energy/btutest.cc
@@ -202,8 +202,8 @@ int main(int argc, char *argv[]) {
 
     // compute numerical error
     heat_flux_at_ice_base.shift(-FF);
-    double max_error = heat_flux_at_ice_base.norm(NORM_INFINITY);
-    double avg_error = heat_flux_at_ice_base.norm(NORM_1);
+    double max_error = heat_flux_at_ice_base.norm(NORM_INFINITY)[0];
+    double avg_error = heat_flux_at_ice_base.norm(NORM_1)[0];
     heat_flux_at_ice_base.shift(+FF); // shift it back for writing
     avg_error /= (grid->Mx() * grid->My());
     log->message(2,

diff --git a/src/fracturedensity/FractureDensity.cc b/src/fracturedensity/FractureDensity.cc
@@ -37,11 +37,10 @@ FractureDensity::FractureDensity(IceGrid::ConstPtr grid,
     m_age_new(grid, "new_fracture_age", WITHOUT_GHOSTS),
     m_toughness(grid, "fracture_toughness", WITHOUT_GHOSTS),
     m_strain_rates(grid, "strain_rates", WITHOUT_GHOSTS,
-                   2,           // stencil width
-                   2),          // dof
-    m_deviatoric_stresses(grid, "sigma", WITHOUT_GHOSTS,
-                          0, // stencil width
-                          3), // dof
+                   2,           // dof
+                   2),          // stencil width
+    m_deviatoric_stresses(grid, "sigma",
+                          WITHOUT_GHOSTS, 3),
     m_velocity(grid, "ghosted_velocity", WITH_GHOSTS, 1),
     m_flow_law(flow_law) {
 
@@ -474,8 +473,8 @@ void FractureDensity::update(double dt,
     }
   }
 
-  A_new.update_ghosts(A);
-  D_new.update_ghosts(D);
+  A.copy_from(A_new);
+  D.copy_from(D_new);
 }
 
 DiagnosticList FractureDensity::diagnostics_impl() const {

diff --git a/src/fracturedensity/FractureDensity.hh b/src/fracturedensity/FractureDensity.hh
@@ -1,4 +1,4 @@
-/* Copyright (C) 2019 PISM Authors
+/* Copyright (C) 2019, 2021 PISM Authors
  *
  * This file is part of PISM.
  *
@@ -26,7 +26,6 @@
 namespace pism {
 
 class IceModelVec2S;
-class IceModelVec2;
 class Geometry;
 
 class FractureDensity : public Component {
@@ -69,10 +68,10 @@ private:
   IceModelVec2S m_toughness;
 
   //! major and minor principal components of horizontal strain-rate tensor (temporary storage)
-  IceModelVec2 m_strain_rates;
+  IceModelVec3 m_strain_rates;
 
   //! components of horizontal stress tensor along axes and shear stress (temporary storage)
-  IceModelVec2 m_deviatoric_stresses;
+  IceModelVec3 m_deviatoric_stresses;
 
   //! Ghosted copy of the ice velocity
   IceModelVec2V m_velocity;

diff --git a/src/frontretreat/FrontRetreat.cc b/src/frontretreat/FrontRetreat.cc
@@ -192,7 +192,7 @@ void FrontRetreat::update_geometry(double dt,
         Href_old = Href(i, j);
 
       // Compute the number of floating neighbors and the neighbor-averaged ice thickness:
-      double H_threshold = part_grid_threshold_thickness(m_cell_type.int_star(i, j),
+      double H_threshold = part_grid_threshold_thickness(m_cell_type.star(i, j),
                                                          ice_thickness.star(i, j),
                                                          surface_elevation.star(i, j),
                                                          bed(i, j));
@@ -215,8 +215,8 @@ void FrontRetreat::update_geometry(double dt,
         int N = 0;
         {
           auto
-            M  = m_cell_type.int_star(i, j),
-            BC = bc_mask.int_star(i, j);
+            M  = m_cell_type.star(i, j),
+            BC = bc_mask.star(i, j);
 
           auto
             b  = bed.star(i, j),

diff --git a/src/frontretreat/calving/EigenCalving.cc b/src/frontretreat/calving/EigenCalving.cc
@@ -1,4 +1,4 @@
-/* Copyright (C) 2013, 2014, 2015, 2016, 2017, 2018, 2019 PISM Authors
+/* Copyright (C) 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2021 PISM Authors
  *
  * This file is part of PISM.
  *

diff --git a/src/frontretreat/calving/HayhurstCalving.cc b/src/frontretreat/calving/HayhurstCalving.cc
@@ -1,4 +1,4 @@
-/* Copyright (C) 2016, 2017, 2018, 2019 PISM Authors
+/* Copyright (C) 2016, 2017, 2018, 2019, 2021 PISM Authors
  *
  * This file is part of PISM.
  *
@@ -137,7 +137,7 @@ void HayhurstCalving::update(const IceModelVec2CellType &cell_type,
     if (cell_type.ice_free(i, j) and cell_type.next_to_ice(i, j) ) {
 
       auto R = m_calving_rate.star(i, j);
-      auto M = cell_type.int_star(i, j);
+      auto M = cell_type.star(i, j);
 
       int N = 0;
       double R_sum = 0.0;

diff --git a/src/frontretreat/calving/StressCalving.cc b/src/frontretreat/calving/StressCalving.cc
@@ -1,4 +1,4 @@
-/* Copyright (C) 2016, 2017, 2018, 2019, 2020 PISM Authors
+/* Copyright (C) 2016, 2017, 2018, 2019, 2020, 2021 PISM Authors
  *
  * This file is part of PISM.
  *
@@ -26,9 +26,8 @@ StressCalving::StressCalving(IceGrid::ConstPtr grid,
                              unsigned int stencil_width)
   : Component(grid),
     m_stencil_width(stencil_width),
-    m_strain_rates(m_grid, "strain_rates", WITH_GHOSTS,
-                   m_stencil_width,
-                   2 /* 2 components */),
+    m_strain_rates(m_grid, "strain_rates", WITH_GHOSTS, {"eigen1", "eigen2"},
+                   m_stencil_width),
     m_calving_rate(m_grid, "calving_rate", WITHOUT_GHOSTS),
     m_cell_type(m_grid, "cell_type", WITH_GHOSTS)
 {

diff --git a/src/frontretreat/calving/StressCalving.hh b/src/frontretreat/calving/StressCalving.hh
@@ -1,4 +1,4 @@
-/* Copyright (C) 2016, 2018, 2019 PISM Authors
+/* Copyright (C) 2016, 2018, 2019, 2021 PISM Authors
  *
  * This file is part of PISM.
  *
@@ -38,7 +38,7 @@ public:
 protected:
   const int m_stencil_width;
 
-  IceModelVec2 m_strain_rates;
+  IceModelVec3 m_strain_rates;
 
   IceModelVec2S m_calving_rate;