Minor cleanup inspired by clang-tidy

src/energy/enthSystem.cc

@@ -1,4 +1,4 @@
-// Copyright (C) 2009-2018, 2020 Andreas Aschwanden and Ed Bueler and Constantine Khroulev
+// Copyright (C) 2009-2018, 2020, 2021 Andreas Aschwanden and Ed Bueler and Constantine Khroulev
 //
 // This file is part of PISM.
 //

src/stressbalance/blatter/Blatter.cc

@@ -89,7 +89,7 @@ void Blatter::compute_node_type(double min_thickness) {
       for (int k = 0; k < fem::q1::n_chi; ++k) {
         int ii, jj;
         E.local_to_global(k, ii, jj);
-        node_type(ii, jj) += interior;
+        node_type(ii, jj) += static_cast<double>(interior);
       }
     }
   }
@@ -189,7 +189,7 @@ bool Blatter::grounding_line(const double *F) {
  * over this face when computing lateral boundary conditions.
  */
 bool Blatter::partially_submerged_face(int face, const double *z, const double *sea_level) {
-  auto nodes = fem::q13d::incident_nodes[face];
+  const auto *nodes = fem::q13d::incident_nodes[face];
 
   // number of nodes per face
   int N = 4;
@@ -225,7 +225,7 @@ bool Blatter::marine_boundary(int face,
                               const int *node_type,
                               const double *ice_bottom,
                               const double *sea_level) {
-  auto nodes = fem::q13d::incident_nodes[face];
+  const auto *nodes = fem::q13d::incident_nodes[face];
 
   // number of nodes per face
   int N = 4;
@@ -449,8 +449,8 @@ PetscErrorCode Blatter::setup(DM pism_da, Periodicity periodicity, int Mz,
       stencil_width = 1;
 
     DMBoundaryType
-      bx = periodicity & X_PERIODIC ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE,
-      by = periodicity & Y_PERIODIC ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE,
+      bx = (periodicity & X_PERIODIC) != 0 ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE,
+      by = (periodicity & Y_PERIODIC) != 0 ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE,
       bz = DM_BOUNDARY_NONE;
 
     ierr = DMDACreate3d(comm,
@@ -512,7 +512,7 @@ PetscErrorCode Blatter::setup(DM pism_da, Periodicity periodicity, int Mz,
 
     PetscBool ksp_use_ew = PETSC_FALSE;
     ierr = SNESKSPGetUseEW(m_snes, &ksp_use_ew); CHKERRQ(ierr);
-    m_ksp_use_ew = ksp_use_ew;
+    m_ksp_use_ew = (ksp_use_ew != 0U);
   }
 
   return 0;
@@ -570,7 +570,7 @@ void Blatter::init_2d_parameters(const Inputs &inputs) {
  */
 void Blatter::init_ice_hardness(const Inputs &inputs, const petsc::DM &da) {
 
-  auto enthalpy = inputs.enthalpy;
+  const auto *enthalpy = inputs.enthalpy;
   // PISM's vertical grid:
   const auto &zlevels = enthalpy->levels();
   auto Mz = zlevels.size();
@@ -639,15 +639,15 @@ void Blatter::nodal_parameter_values(const fem::Q1Element3 &element,
 
     auto p = P[I.j][I.i];
 
-    node_type[n]        = p.node_type;
+    node_type[n]        = static_cast<int>(p.node_type);
     bottom_elevation[n] = p.bed;
     ice_thickness[n]    = p.thickness;
 
-    if (surface_elevation) {
+    if (surface_elevation != nullptr) {
       surface_elevation[n] = p.bed + p.thickness;
     }
 
-    if (sea_level) {
+    if (sea_level != nullptr) {
       sea_level[n] = p.sea_level;
     }
   }
@@ -714,7 +714,7 @@ void Blatter::report_mesh_info() {
 
       int node_type[4];
       for (int k = 0; k < 4; ++k) {
-        node_type[k] = P[k].node_type;
+        node_type[k] = static_cast<int>(P[k].node_type);
       }
 
       if (exterior_element(node_type)) {
@@ -1103,8 +1103,8 @@ void Blatter::copy_solution() {
   for (Points p(*m_grid); p; p.next()) {
     const int i = p.i(), j = p.j();
 
-    auto u = m_u_sigma->get_column(i, j);
-    auto v = m_v_sigma->get_column(i, j);
+    auto *u = m_u_sigma->get_column(i, j);
+    auto *v = m_v_sigma->get_column(i, j);
 
     for (int k = 0; k < Mz; ++k) {
       u[k] = x[j][i][k].u;      // STORAGE_ORDER
@@ -1143,8 +1143,8 @@ void Blatter::set_initial_guess(const IceModelVec3 &u_sigma,
   for (Points p(*m_grid); p; p.next()) {
     const int i = p.i(), j = p.j();
 
-    auto u = u_sigma.get_column(i, j);
-    auto v = v_sigma.get_column(i, j);
+    const auto *u = u_sigma.get_column(i, j);
+    const auto *v = v_sigma.get_column(i, j);
 
     for (int k = 0; k < Mz; ++k) {
       x[j][i][k].u = u[k];      // STORAGE_ORDER

src/stressbalance/blatter/Blatter.hh

@@ -118,11 +118,11 @@ protected:
 
   void write_model_state_impl(const File &output) const;
 
-  bool exterior_element(const int *node_type);
+  static bool exterior_element(const int *node_type);
 
-  bool grounding_line(const double *F);
+  static bool grounding_line(const double *F);
 
-  bool partially_submerged_face(int face, const double *z, const double *sea_level);
+  static bool partially_submerged_face(int face, const double *z, const double *sea_level);
 
   void compute_node_type(double min_thickness);
 
@@ -160,7 +160,7 @@ protected:
                               const Vector2 *u_nodal,
                               double K[2 * fem::q13d::n_chi][2 * fem::q13d::n_chi]);
 
-  void compute_residual(DMDALocalInfo *info, const Vector2 ***xg, Vector2 ***yg);
+  void compute_residual(DMDALocalInfo *info, const Vector2 ***X, Vector2 ***R);
 
   void residual_dirichlet(const DMDALocalInfo &info,
                           Parameters **P,
@@ -198,11 +198,11 @@ protected:
   static PetscErrorCode jacobian_callback(DMDALocalInfo *info,
                                           const Vector2 ***x,
                                           Mat A, Mat J,
-                                          Blatter *data);
+                                          Blatter *solver);
 
   static PetscErrorCode function_callback(DMDALocalInfo *info,
                                           const Vector2 ***x, Vector2 ***f,
-                                          Blatter *data);
+                                          Blatter *solver);
 
   virtual void init_2d_parameters(const Inputs &inputs);
 

src/stressbalance/blatter/BlatterMod.cc

@@ -1,4 +1,4 @@
-/* Copyright (C) 2020 PISM Authors
+/* Copyright (C) 2020, 2021 PISM Authors
  *
  * This file is part of PISM.
  *
@@ -85,8 +85,8 @@ void BlatterMod::transfer(const IceModelVec2S &ice_thickness) {
   for (Points p(*m_grid); p; p.next()) {
     const int i = p.i(), j = p.j();
 
-    auto u = m_u.get_column(i, j);
-    auto v = m_v.get_column(i, j);
+    auto *u = m_u.get_column(i, j);
+    auto *v = m_v.get_column(i, j);
 
     double H = ice_thickness(i, j);
 
@@ -113,7 +113,7 @@ void BlatterMod::transfer(const IceModelVec2S &ice_thickness) {
 void BlatterMod::compute_max_diffusivity(const IceModelVec2V &velocity,
                                          const IceModelVec2S &ice_thickness,
                                          const IceModelVec2S &surface) {
-  double eps = 1e-3;
+  const double eps = 1e-3;
   double
     dx = m_grid->dx(),
     dy = m_grid->dy();

src/stressbalance/blatter/jacobian.cc

@@ -33,8 +33,8 @@ namespace stressbalance {
  * Computes the Jacobian contribution of the "main" part of the Blatter system.
  */
 void Blatter::jacobian_f(const fem::Q1Element3 &element,
-                         const Vector2 *velocity,
-                         const double *hardness,
+                         const Vector2 *u_nodal,
+                         const double *B_nodal,
                          double K[16][16]) {
   int Nk = fem::q13d::n_chi;
 
@@ -46,8 +46,8 @@ void Blatter::jacobian_f(const fem::Q1Element3 &element,
 
   double *B = m_work[0];
 
-  element.evaluate(velocity, u, u_x, u_y, u_z);
-  element.evaluate(hardness, B);
+  element.evaluate(u_nodal, u, u_x, u_y, u_z);
+  element.evaluate(B_nodal, B);
 
   // loop over all quadrature points
   for (int q = 0; q < element.n_pts(); ++q) {
@@ -122,7 +122,7 @@ void Blatter::jacobian_f(const fem::Q1Element3 &element,
 void Blatter::jacobian_basal(const fem::Q1Element3Face &face,
                              const double *tauc_nodal,
                              const double *f_nodal,
-                             const Vector2 *velocity,
+                             const Vector2 *u_nodal,
                              double K[16][16]) {
   int Nk = fem::q13d::n_chi;
 
@@ -132,7 +132,7 @@ void Blatter::jacobian_basal(const fem::Q1Element3Face &face,
     *tauc       = m_work[0],
     *floatation = m_work[1];
 
-  face.evaluate(velocity, u);
+  face.evaluate(u_nodal, u);
   face.evaluate(tauc_nodal, tauc);
   face.evaluate(f_nodal, floatation);
 
@@ -253,7 +253,7 @@ void Blatter::compute_jacobian(DMDALocalInfo *petsc_info,
   DataAccess<double***> hardness(info.da, 3, GHOSTED);
 
   IceModelVec::AccessList list(m_parameters);
-  auto P = m_parameters.array();
+  auto *P = m_parameters.array();
 
   // loop over all the elements that have at least one owned node
   for (int j = info.gys; j < info.gys + info.gym - 1; j++) {
@@ -320,8 +320,7 @@ void Blatter::compute_jacobian(DMDALocalInfo *petsc_info,
 
           fem::Q1Element3Face *face = grounding_line(floatation) ? &m_face100 : &m_face4;
 
-          // face 4 is the bottom face in fem::q13d::incident_nodes
-          face->reset(4, z);
+          face->reset(fem::q13d::FACE_BOTTOM, z);
 
           jacobian_basal(*face, basal_yield_stress, floatation, velocity, K);
         }

src/stressbalance/blatter/residual.cc

@@ -348,7 +348,7 @@ void Blatter::compute_residual(DMDALocalInfo *petsc_info,
   DataAccess<double***> ice_hardness(info.da, 3, GHOSTED);
 
   IceModelVec::AccessList list(m_parameters);
-  auto P = m_parameters.array();
+  auto *P = m_parameters.array();
 
   // Compute the residual at Dirichlet nodes and set it to zero elsewhere.
   residual_dirichlet(info, P, X, R);
@@ -423,8 +423,7 @@ void Blatter::compute_residual(DMDALocalInfo *petsc_info,
 
           // use an N*N-point equally-spaced quadrature at grounding lines
           fem::Q1Element3Face *face = grounding_line(floatation) ? &m_face100 : &m_face4;
-          // face 4 is the bottom face in fem::q13d::incident_nodes
-          face->reset(4, z);
+          face->reset(fem::q13d::FACE_BOTTOM, z);
 
           residual_basal(element, *face, basal_yield_stress, floatation, velocity, R_nodal);
         }
@@ -444,8 +443,7 @@ void Blatter::compute_residual(DMDALocalInfo *petsc_info,
 
         // top boundary (verification tests only)
         if (k == info.mz - 2) {
-          // face 5 is the top face in fem::q13d::incident_nodes
-          m_face4.reset(5, z);
+          m_face4.reset(fem::q13d::FACE_TOP, z);
 
           residual_surface(element, m_face4, R_nodal);
         }

src/stressbalance/sia/SIAFD.cc

@@ -41,7 +41,7 @@ namespace pism {
 namespace stressbalance {
 
 SIAFD::SIAFD(IceGrid::ConstPtr g)
-  : SSB_Modifier(g),
+  : SSB_Modifier(std::move(g)),
     m_stencil_width(m_config->get_number("grid.max_stencil_width")),
     m_work_2d_0(m_grid, "work_vector_2d_0", WITH_GHOSTS, m_stencil_width),
     m_work_2d_1(m_grid, "work_vector_2d_1", WITH_GHOSTS, m_stencil_width),
@@ -766,7 +766,9 @@ void SIAFD::compute_diffusivity(bool full_update,
                                   "too rough.\n"
                                   "Increase stress_balance.sia.max_diffusivity to suppress this message.", m_D_max);
 
-  } else if (high_diffusivity_counter > 0) {
+  }
+
+  if (high_diffusivity_counter > 0) {
     // This can happen only if stress_balance.sia.limit_diffusivity is true and this
     // limiting mechanism was active (high_diffusivity_counter is incremented only if
     // limit_diffusivity is true).
@@ -955,16 +957,16 @@ void SIAFD::compute_3d_horizontal_velocity(const Geometry &geometry,
 }
 
 //! Determine if `accumulation_time` corresponds to an interglacial period.
-bool SIAFD::interglacial(double accumulation_time) {
+bool SIAFD::interglacial(double accumulation_time) const {
   if (accumulation_time < m_eemian_start) {
     return false;
-  } else if (accumulation_time < m_eemian_end) {
-    return true;
-  } else if (accumulation_time < m_holocene_start) {
-    return false;
-  } else {
+  }
+
+  if (accumulation_time < m_eemian_end) {
     return true;
   }
+
+  return (accumulation_time >= m_holocene_start);
 }
 
 const IceModelVec2Stag& SIAFD::surface_gradient_x() const {

src/stressbalance/sia/SIAFD.hh

@@ -101,7 +101,7 @@ protected:
 
   virtual void compute_I(const Geometry &geometry);
 
-  bool interglacial(double accumulation_time);
+  bool interglacial(double accumulation_time) const;
 
   const unsigned int m_stencil_width;
 

src/util/ColumnInterpolation.cc

@@ -1,4 +1,4 @@
-/* Copyright (C) 2014, 2015 PISM Authors
+/* Copyright (C) 2014, 2015, 2021 PISM Authors
  *
  * This file is part of PISM.
  *
@@ -184,7 +184,9 @@ static std::vector<unsigned int> init_interpolation_indexes(const std::vector<do
     if (z_output[k] <= z_input.front()) {
       result[k] = 0;
       continue;
-    } else if (z_output[k] >= z_input.back()) {
+    }
+
+    if (z_output[k] >= z_input.back()) {
       result[k] = z_input.size() - 1;
       continue;
     }
@@ -216,11 +218,8 @@ void ColumnInterpolation::init_interpolation() {
     dz_max = std::max(dz, dz_max);
   }
 
-  if (fabs(dz_max - dz_min) <= 1.0e-8) {
-    m_use_linear_interpolation = true;
-  } else {
-    m_use_linear_interpolation = false;
-  }
+  const double eps = 1.0e-8;
+  m_use_linear_interpolation = (fabs(dz_max - dz_min) <= eps);
 
   // initialize quadratic interpolation constants
   if (not m_use_linear_interpolation) {