src/base/stressbalance/SSAFD.cc

Go to the documentation of this file.

// Copyright (C) 2004--2011 Constantine Khroulev, Ed Bueler and Jed Brown
//
// 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 2 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 "SSAFD.hh"
#include "Mask.hh"

SSA *SSAFDFactory(IceGrid &g, IceBasalResistancePlasticLaw &b,
                IceFlowLaw &i, EnthalpyConverter &ec,
                const NCConfigVariable &c)
{
  return new SSAFD(g,b,i,ec,c);
}



PetscErrorCode SSAFD::allocate_fd() {
  PetscErrorCode ierr;

  // note SSADA and SSAX are allocated in SSA::allocate()
  ierr = VecDuplicate(SSAX, &SSARHS); CHKERRQ(ierr);

  ierr = DAGetMatrix(SSADA, MATAIJ, &SSAStiffnessMatrix); CHKERRQ(ierr);

  ierr = KSPCreate(grid.com, &SSAKSP); CHKERRQ(ierr);
  // the default PC type somehow is ILU, which now fails (?) while block jacobi
  //   seems to work; runtime options can override (see test J in vfnow.py)
  PC pc;
  ierr = KSPGetPC(SSAKSP,&pc); CHKERRQ(ierr);
  ierr = PCSetType(pc,PCBJACOBI); CHKERRQ(ierr);
  ierr = KSPSetFromOptions(SSAKSP); CHKERRQ(ierr);

  const PetscScalar power = 1.0 / ice.exponent();
  char unitstr[TEMPORARY_STRING_LENGTH];
  snprintf(unitstr, sizeof(unitstr), "Pa s%f", power);
  ierr = hardness.create(grid, "hardness", false); CHKERRQ(ierr);
  ierr = hardness.set_attrs("diagnostic",
                            "vertically-averaged ice hardness",
                            unitstr, ""); CHKERRQ(ierr);

  ierr = nuH.create(grid, "nuH", true); CHKERRQ(ierr);
  ierr = nuH.set_attrs("internal",
                       "ice thickness times effective viscosity",
                       "Pa s m", ""); CHKERRQ(ierr);

  ierr = nuH_old.create(grid, "nuH_old", true); CHKERRQ(ierr);
  ierr = nuH_old.set_attrs("internal",
                           "ice thickness times effective viscosity (before an update)",
                           "Pa s m", ""); CHKERRQ(ierr);

  scaling = 1.0e9;  // comparable to typical beta for an ice stream;

  // The nuH viewer:
  view_nuh = false;
  nuh_viewer_size = 300;
  nuh_viewer = PETSC_NULL;
  return 0;
}

PetscErrorCode SSAFD::deallocate_fd() {
  PetscErrorCode ierr;

  if (SSAKSP != PETSC_NULL) {
    ierr = KSPDestroy(SSAKSP); CHKERRQ(ierr);
  }

  if (SSAStiffnessMatrix != PETSC_NULL) {
    ierr = MatDestroy(SSAStiffnessMatrix); CHKERRQ(ierr);
  }

  if (SSARHS != PETSC_NULL) {
    ierr = VecDestroy(SSARHS); CHKERRQ(ierr);
  }

  return 0;
}

PetscErrorCode SSAFD::init(PISMVars &vars) {
  PetscErrorCode ierr;
  ierr = SSA::init(vars); CHKERRQ(ierr);
  ierr = verbPrintf(2,grid.com,
                    "  [using the KSP-based finite difference implementation]\n"); CHKERRQ(ierr);

  ierr = PetscOptionsBegin(grid.com, "", "SSAFD options", ""); CHKERRQ(ierr);
  {
    bool flag;
    ierr = PISMOptionsInt("-ssa_nuh_viewer_size", "nuH viewer size",
                          nuh_viewer_size, flag); CHKERRQ(ierr);
    ierr = PISMOptionsIsSet("-ssa_view_nuh", "Enable the SSAFD nuH runtime viewer",
                            view_nuh); CHKERRQ(ierr);
  }
  ierr = PetscOptionsEnd(); CHKERRQ(ierr);

  if (config.get_flag("calving_front_stress_boundary_condition")) {
    ierr = verbPrintf(2,grid.com,
      "  using PISM-PIK calving-front stress boundary condition ...\n"); CHKERRQ(ierr);
  }

  return 0;
}


PetscErrorCode SSAFD::assemble_rhs(Vec rhs) {
  PetscErrorCode ierr;
  PISMVector2 **rhs_uv;
  const double dx = grid.dx, dy = grid.dy;

  Mask M;

  const double standard_gravity = config.get("standard_gravity"),
    ocean_rho = config.get("sea_water_density");
  const bool use_cfbc = config.get_flag("calving_front_stress_boundary_condition");

  ierr = VecSet(rhs, 0.0); CHKERRQ(ierr);

  // get driving stress components
  ierr = compute_driving_stress(taud); CHKERRQ(ierr);

  ierr = taud.begin_access(); CHKERRQ(ierr);
  ierr = DAVecGetArray(SSADA, rhs, &rhs_uv); CHKERRQ(ierr);

  if (vel_bc && bc_locations) {
    ierr = vel_bc->begin_access(); CHKERRQ(ierr);
    ierr = bc_locations->begin_access(); CHKERRQ(ierr);
  }

  if (use_cfbc) {
    ierr = thickness->begin_access(); CHKERRQ(ierr);
    ierr = bed->begin_access(); CHKERRQ(ierr);
    ierr = mask->begin_access(); CHKERRQ(ierr);
  }

  for (PetscInt i = grid.xs; i < grid.xs + grid.xm; ++i) {
    for (PetscInt j = grid.ys; j < grid.ys + grid.ym; ++j) {

      if (vel_bc && (bc_locations->as_int(i, j) == 1)) {
        rhs_uv[i][j].u = scaling * (*vel_bc)(i, j).u;
        rhs_uv[i][j].v = scaling * (*vel_bc)(i, j).v;
        continue;
      }


      if (use_cfbc) {
        PetscScalar H_ij = (*thickness)(i,j);
        PetscInt M_ij = mask->as_int(i,j),
          M_e = mask->as_int(i + 1,j),
          M_w = mask->as_int(i - 1,j),
          M_n = mask->as_int(i,j + 1),
          M_s = mask->as_int(i,j - 1);

        // Note: this sets velocities at both ice-free ocean and ice-free
        // bedrock to zero. This means that we need to set boundary conditions
        // at both ice/ice-free-ocean and ice/ice-free-bedrock interfaces below
        // to be consistent.
        if (M.ice_free(M_ij)) {
          rhs_uv[i][j].u = 0.0;
          rhs_uv[i][j].v = 0.0;
          continue;
        }

        if (is_marginal(i, j)) {
          PetscInt aMM = 1, aPP = 1, bMM = 1, bPP = 1;
          // direct neighbors
          if (M.ice_free(M_e)) aPP = 0;
          if (M.ice_free(M_w)) aMM = 0;
          if (M.ice_free(M_n)) bPP = 0;
          if (M.ice_free(M_s)) bMM = 0;

          double ice_pressure = ice.rho * standard_gravity * H_ij,
            ocean_pressure;

          const double h_grounded = (*bed)(i,j) + H_ij,
            h_floating = sea_level + (1.0 - ice.rho / ocean_rho) * H_ij,
            H_ij2 = H_ij*H_ij;

          //double h_ij = 0.0, tdx = 0.0, tdy = 0.0;
          double h_ij = 0.0, tdx = taud(i, j).u, tdy = taud(i, j).v;

          if ((*bed)(i,j) < (sea_level - (ice.rho / ocean_rho) * H_ij)) {
            //calving front boundary condition for floating shelf
            ocean_pressure = 0.5 * ice.rho * standard_gravity * (1 - (ice.rho / ocean_rho))*H_ij2;
            // this is not really the ocean_pressure, but the difference between
            // ocean_pressure and isotrop.normal stresses (=pressure) from within
            // the ice
            h_ij = h_floating;
          } else {
            h_ij = h_grounded;
            if( (*bed)(i,j) >= sea_level) {
              // boundary condition for a "cliff" (grounded ice next to
              // ice-free ocean) or grounded margin.
              ocean_pressure = 0.5 * ice.rho * standard_gravity * H_ij2;
              // this is not 'zero' because the isotrop.normal stresses
              // (=pressure) from within the ice figures on RHS
            } else {
              // boundary condition for marine terminating glacier
              ocean_pressure = 0.5 * ice.rho * standard_gravity *
                (H_ij2 - (ocean_rho / ice.rho)*(sea_level - (*bed)(i,j))*(sea_level - (*bed)(i,j)));
            }
          }

          //here we take the direct gradient at the boundary (not centered)

          if (aPP == 0 && aMM == 1) tdx = ice_pressure*h_ij / dx;
          else if (aMM == 0 && aPP == 1) tdx = -ice_pressure*h_ij / dx;
          else if (aPP == 0 && aMM == 0) tdx = 0; //in case of some kind of ice nose, or ice bridge

          if (bPP == 0 && bMM == 1) tdy = ice_pressure*h_ij / dy;
          else if (bMM == 0 && bPP == 1) tdy = -ice_pressure*h_ij / dy;
          else if (bPP == 0 && bMM == 0) tdy = 0;

          // Note that if the current cell is "marginal" but not a CFBC
          // location, the following two lines are equaivalent to the "usual
          // case" below.
          rhs_uv[i][j].u = tdx - (aMM - aPP)*ocean_pressure / dx;
          rhs_uv[i][j].v = tdy - (bMM - bPP)*ocean_pressure / dy;

          continue;
        } // end of "if (is_marginal(i, j))"

        // If we reached this point, then CFBC are enabled, but we are in the
        // interior of a sheet or shelf. See "usual case" below.

      }   // end of "if (use_cfbc)"

      // usual case: use already computed driving stress
      rhs_uv[i][j].u = taud(i, j).u;
      rhs_uv[i][j].v = taud(i, j).v;
    }
  }

  if (use_cfbc) {
    ierr = thickness->end_access(); CHKERRQ(ierr);
    ierr = bed->end_access(); CHKERRQ(ierr);
    ierr = mask->end_access(); CHKERRQ(ierr);
  }

  if (vel_bc && bc_locations) {
    ierr = bc_locations->end_access(); CHKERRQ(ierr);
    ierr = vel_bc->end_access(); CHKERRQ(ierr);
  }

  ierr = taud.end_access(); CHKERRQ(ierr);
  ierr = DAVecRestoreArray(SSADA, rhs, &rhs_uv); CHKERRQ(ierr);

  return 0;
}


PetscErrorCode SSAFD::assemble_matrix(bool include_basal_shear, Mat A) {
  PetscErrorCode  ierr;

  const PetscScalar   dx=grid.dx, dy=grid.dy;
  const PetscScalar   beta_ice_free_bedrock = config.get("beta_ice_free_bedrock");
  const bool use_cfbc = config.get_flag("calving_front_stress_boundary_condition");

  // shortcut:
  IceModelVec2V &vel = velocity;

  ierr = MatZeroEntries(A); CHKERRQ(ierr);

  /* matrix assembly loop */
  
  ierr = nuH.begin_access(); CHKERRQ(ierr);
  ierr = tauc->begin_access(); CHKERRQ(ierr);
  ierr = vel.begin_access(); CHKERRQ(ierr);
  ierr = mask->begin_access(); CHKERRQ(ierr);

  Mask M;

  if (vel_bc && bc_locations) {
    ierr = bc_locations->begin_access(); CHKERRQ(ierr);
  }

  for (PetscInt i=grid.xs; i<grid.xs+grid.xm; ++i) {
    for (PetscInt j=grid.ys; j<grid.ys+grid.ym; ++j) {

      // Handle the easy case: provided Dirichlet boundary conditions
      if (vel_bc && bc_locations && bc_locations->as_int(i,j) == 1) {
        // set diagonal entry to one (scaled); RHS entry will be known velocity;
        ierr = set_diagonal_matrix_entry(A, i, j, scaling); CHKERRQ(ierr); 
        continue;
      }

      /* Provide shorthand for the following staggered coefficients  nu H:
       *      c_n
       *  c_w     c_e
       *      c_s
       */
      const PetscScalar c_w = nuH(i-1,j,0);
      const PetscScalar c_e = nuH(i,j,0);
      const PetscScalar c_s = nuH(i,j-1,1);
      const PetscScalar c_n = nuH(i,j,1);

      // We use DAGetMatrix to obtain the SSA matrix, which means that all 18
      // non-zeros get allocated, even though we use only 13 (or 14). The
      // remaining 5 (or 4) coefficients are zeros, but we set them anyway,
      // because this makes the code easier to understand.
      const PetscInt sten = 18;
      MatStencil row, col[sten];

      PetscInt aMn = 1, aPn = 1, aMM = 1, aPP = 1, aMs = 1, aPs = 1;
      PetscInt bPw = 1, bPP = 1, bPe = 1, bMw = 1, bMM = 1, bMe = 1;

      PetscInt M_ij = mask->as_int(i,j);

      if (use_cfbc) {
        int
          // direct neighbors
          M_e = mask->as_int(i + 1,j),
          M_w = mask->as_int(i - 1,j),
          M_n = mask->as_int(i,j + 1),
          M_s = mask->as_int(i,j - 1),
          // "diagonal" neighbors
          M_ne = mask->as_int(i + 1,j + 1),
          M_se = mask->as_int(i + 1,j - 1),
          M_nw = mask->as_int(i - 1,j + 1),
          M_sw = mask->as_int(i - 1,j - 1);

        // Note: this sets velocities at both ice-free ocean and ice-free
        // bedrock to zero. This means that we need to set boundary conditions
        // at both ice/ice-free-ocean and ice/ice-free-bedrock interfaces below
        // to be consistent.
        if (M.ice_free(M_ij)) {
          ierr = set_diagonal_matrix_entry(A, i, j, scaling); CHKERRQ(ierr);
          continue;
        }

        if (is_marginal(i, j)) {
          // If at least one of the following four conditions is "true", we're
          // at a CFBC location.
          if (M.ice_free(M_e)) aPP = 0;
          if (M.ice_free(M_w)) aMM = 0;
          if (M.ice_free(M_n)) bPP = 0;
          if (M.ice_free(M_s)) bMM = 0;

          // decide whether to use centered or one-sided differences
          if (M.ice_free(M_n) || M.ice_free(M_ne)) aPn = 0;
          if (M.ice_free(M_e) || M.ice_free(M_ne)) bPe = 0;
          if (M.ice_free(M_e) || M.ice_free(M_se)) bMe = 0;
          if (M.ice_free(M_s) || M.ice_free(M_se)) aPs = 0;
          if (M.ice_free(M_s) || M.ice_free(M_sw)) aMs = 0;
          if (M.ice_free(M_w) || M.ice_free(M_sw)) bMw = 0;
          if (M.ice_free(M_w) || M.ice_free(M_nw)) bPw = 0;
          if (M.ice_free(M_n) || M.ice_free(M_nw)) aMn = 0;
        }
      } // end of "if (use_cfbc)"

      /* begin Maxima-generated code */
      const PetscReal dx2 = dx*dx, dy2 = dy*dy, d4 = 4*dx*dy, d2 = 2*dx*dy;

      /* Coefficients of the discretization of the first equation; u first, then v. */
      PetscReal eq1[] = {
        0,  -c_n*bPP/dy2,  0, 
        -4*c_w*aMM/dx2,  (c_n*bPP+c_s*bMM)/dy2+(4*c_e*aPP+4*c_w*aMM)/dx2,  -4*c_e*aPP/dx2, 
        0,  -c_s*bMM/dy2,  0, 
        c_w*aMM*bPw/d2+c_n*aMn*bPP/d4,  (c_n*aPn*bPP-c_n*aMn*bPP)/d4+(c_w*aMM*bPP-c_e*aPP*bPP)/d2,  -c_e*aPP*bPe/d2-c_n*aPn*bPP/d4, 
        (c_w*aMM*bMw-c_w*aMM*bPw)/d2+(c_n*aMM*bPP-c_s*aMM*bMM)/d4,  (c_n*aPP*bPP-c_n*aMM*bPP-c_s*aPP*bMM+c_s*aMM*bMM)/d4+(c_e*aPP*bPP-c_w*aMM*bPP-c_e*aPP*bMM+c_w*aMM*bMM)/d2,  (c_e*aPP*bPe-c_e*aPP*bMe)/d2+(c_s*aPP*bMM-c_n*aPP*bPP)/d4, 
        -c_w*aMM*bMw/d2-c_s*aMs*bMM/d4,  (c_s*aMs*bMM-c_s*aPs*bMM)/d4+(c_e*aPP*bMM-c_w*aMM*bMM)/d2,  c_e*aPP*bMe/d2+c_s*aPs*bMM/d4, 
      };

      /* Coefficients of the discretization of the second equation; u first, then v. */
      PetscReal eq2[] = {
        c_w*aMM*bPw/d4+c_n*aMn*bPP/d2,  (c_n*aPn*bPP-c_n*aMn*bPP)/d2+(c_w*aMM*bPP-c_e*aPP*bPP)/d4,  -c_e*aPP*bPe/d4-c_n*aPn*bPP/d2, 
        (c_w*aMM*bMw-c_w*aMM*bPw)/d4+(c_n*aMM*bPP-c_s*aMM*bMM)/d2,  (c_n*aPP*bPP-c_n*aMM*bPP-c_s*aPP*bMM+c_s*aMM*bMM)/d2+(c_e*aPP*bPP-c_w*aMM*bPP-c_e*aPP*bMM+c_w*aMM*bMM)/d4,  (c_e*aPP*bPe-c_e*aPP*bMe)/d4+(c_s*aPP*bMM-c_n*aPP*bPP)/d2, 
        -c_w*aMM*bMw/d4-c_s*aMs*bMM/d2,  (c_s*aMs*bMM-c_s*aPs*bMM)/d2+(c_e*aPP*bMM-c_w*aMM*bMM)/d4,  c_e*aPP*bMe/d4+c_s*aPs*bMM/d2, 
        0,  -4*c_n*bPP/dy2,  0, 
        -c_w*aMM/dx2,  (4*c_n*bPP+4*c_s*bMM)/dy2+(c_e*aPP+c_w*aMM)/dx2,  -c_e*aPP/dx2, 
        0,  -4*c_s*bMM/dy2,  0, 
      };

      /* i indices */
      const PetscInt I[] = {
        i-1,  i,  i+1, 
        i-1,  i,  i+1, 
        i-1,  i,  i+1, 
        i-1,  i,  i+1, 
        i-1,  i,  i+1, 
        i-1,  i,  i+1, 
      };

      /* j indices */
      const PetscInt J[] = {
        j+1,  j+1,  j+1, 
        j,  j,  j, 
        j-1,  j-1,  j-1, 
        j+1,  j+1,  j+1, 
        j,  j,  j, 
        j-1,  j-1,  j-1, 
      };

      /* component indices */
      const PetscInt C[] = {
        0,  0,  0, 
        0,  0,  0, 
        0,  0,  0, 
        1,  1,  1, 
        1,  1,  1, 
        1,  1,  1, 
      };
      /* end Maxima-generated code */

      /* Dragging ice experiences friction at the bed determined by the
       *    IceBasalResistancePlasticLaw::drag() methods.  These may be a plastic,
       *    pseudo-plastic, or linear friction law.  Dragging is done implicitly
       *    (i.e. on left side of SSA eqns).  */
      PetscReal beta = 0.0;
      if (include_basal_shear) {
        if (M.grounded_ice(M_ij)) {
          beta = basal.drag((*tauc)(i,j), vel(i,j).u, vel(i,j).v);
        } else if (M.ice_free_land(M_ij)) {
          // apply drag even in this case, to help with margins; note ice free
          // areas already have a strength extension
          beta = beta_ice_free_bedrock;
        }
      }

      // add beta to diagonal entries
      eq1[4]  += beta;
      eq2[13] += beta;

      // build equations: NOTE TRANSPOSE
      row.j = i; row.i = j;
      for (PetscInt m = 0; m < sten; m++) {
        col[m].j = I[m]; col[m].i = J[m]; col[m].c = C[m];
      }

      // set coefficients of the first equation:
      row.c = 0;
      ierr = MatSetValuesStencil(A, 1, &row, sten, col, eq1, INSERT_VALUES); CHKERRQ(ierr);

      // set coefficients of the second equation:
      row.c = 1;
      ierr = MatSetValuesStencil(A, 1, &row, sten, col, eq2, INSERT_VALUES); CHKERRQ(ierr);
    }
  }

  if (vel_bc && bc_locations) {
    ierr = bc_locations->end_access(); CHKERRQ(ierr);
  }

  ierr = mask->end_access(); CHKERRQ(ierr);
  ierr = vel.end_access(); CHKERRQ(ierr);
  ierr = tauc->end_access(); CHKERRQ(ierr);
  ierr = nuH.end_access(); CHKERRQ(ierr);

  ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
#ifdef PISM_DEBUG
  ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
#endif

  return 0;
}



PetscErrorCode SSAFD::solve() {
  PetscErrorCode ierr;
  Mat A = SSAStiffnessMatrix; // solve  A SSAX = SSARHS
  PetscReal   norm, normChange;
  PetscInt    ksp_iterations, outer_iterations;
  KSPConvergedReason  reason;

  stdout_ssa = "";

  bool write_ssa_system_to_matlab = config.get_flag("write_ssa_system_to_matlab");
  PetscReal ssaRelativeTolerance = config.get("ssafd_relative_convergence"),
            epsilon              = config.get("epsilon_ssafd");

  PetscInt ssaMaxIterations = static_cast<PetscInt>(config.get("max_iterations_ssafd"));

  ierr = velocity.copy_to(velocity_old); CHKERRQ(ierr);

  // computation of RHS only needs to be done once; does not depend on
  // solution; but matrix changes under nonlinear iteration (loop over k below)
  ierr = assemble_rhs(SSARHS); CHKERRQ(ierr);

  ierr = compute_hardav_staggered(hardness); CHKERRQ(ierr);

  for (PetscInt l=0; ; ++l) { // iterate with increasing regularization parameter
    ierr = compute_nuH_staggered(nuH, epsilon); CHKERRQ(ierr);

    ierr = update_nuH_viewers(); CHKERRQ(ierr);
    // iterate on effective viscosity: "outer nonlinear iteration":
    for (PetscInt k = 0; k < ssaMaxIterations; ++k) {
      if (getVerbosityLevel() > 2) {
        char tempstr[50] = "";  snprintf(tempstr,50, "  %d,%2d:", l, k);
        stdout_ssa += tempstr;
      }

      // in preparation of measuring change of effective viscosity:
      ierr = nuH.copy_to(nuH_old); CHKERRQ(ierr);

      // assemble (or re-assemble) matrix, which depends on updated viscosity
      ierr = assemble_matrix(true, A); CHKERRQ(ierr);
      if (getVerbosityLevel() > 2)
        stdout_ssa += "A:";

      // call PETSc to solve linear system by iterative method; "inner linear iteration"
      ierr = KSPSetOperators(SSAKSP, A, A, SAME_NONZERO_PATTERN); CHKERRQ(ierr);
      ierr = KSPSolve(SSAKSP, SSARHS, SSAX); CHKERRQ(ierr); // SOLVE

      // report to standard out about iteration
      ierr = KSPGetConvergedReason(SSAKSP, &reason); CHKERRQ(ierr);
      if (reason < 0) {
        ierr = verbPrintf(1,grid.com,
            "\n\n\nPISM ERROR:  KSPSolve() reports 'diverged'; reason = %d = '%s';\n"
                  "  see PETSc man page for KSPGetConvergedReason();   ENDING ...\n\n",
            reason,KSPConvergedReasons[reason]); CHKERRQ(ierr);
        PISMEnd();
      }
      ierr = KSPGetIterationNumber(SSAKSP, &ksp_iterations); CHKERRQ(ierr);
      if (getVerbosityLevel() > 2) {
        char tempstr[50] = "";  snprintf(tempstr,50, "S:%d,%d: ", ksp_iterations, reason);
        stdout_ssa += tempstr;
      }

      // Communicate so that we have stencil width for evaluation of effective
      //   viscosity on next "outer" iteration (and geometry etc. if done):
      ierr = velocity.copy_from(SSAX); CHKERRQ(ierr);

      ierr = velocity.beginGhostComm(); CHKERRQ(ierr);
      ierr = velocity.endGhostComm(); CHKERRQ(ierr);

      // update viscosity and check for viscosity convergence
      ierr = compute_nuH_staggered(nuH, epsilon); CHKERRQ(ierr);
      ierr = update_nuH_viewers(); CHKERRQ(ierr);
      ierr = compute_nuH_norm(norm, normChange); CHKERRQ(ierr);
      if (getVerbosityLevel() > 2) {
        char tempstr[100] = "";
        snprintf(tempstr,100, "|nu|_2, |Delta nu|_2/|nu|_2 = %10.3e %10.3e\n",
                         norm, normChange/norm);
        stdout_ssa += tempstr;
      }

      if (getVerbosityLevel() > 2) { // assume that high verbosity shows interest
                                     //   in immediate feedback about SSA iterations
        ierr = verbPrintf(2,grid.com, stdout_ssa.c_str()); CHKERRQ(ierr);
        stdout_ssa = "";
      }

      outer_iterations = k + 1;
      if (norm == 0 || normChange / norm < ssaRelativeTolerance) goto done;

    } // end of the "outer loop" (index: k)

    if (epsilon > 0.0) {
       // this has no units; epsilon goes up by this ratio when previous value failed
       const PetscScalar DEFAULT_EPSILON_MULTIPLIER_SSA = 4.0;
       ierr = verbPrintf(1,grid.com,
                         "WARNING: Effective viscosity not converged after %d iterations\n"
                         "\twith epsilon=%8.2e. Retrying with epsilon * %8.2e.\n",
                         ssaMaxIterations, epsilon, DEFAULT_EPSILON_MULTIPLIER_SSA);
       CHKERRQ(ierr);

       ierr = velocity.copy_from(velocity_old); CHKERRQ(ierr);
       epsilon *= DEFAULT_EPSILON_MULTIPLIER_SSA;
    } else {
       SETERRQ1(1,
         "Effective viscosity not converged after %d iterations; epsilon=0.0.\n"
         "  Stopping.\n",
         ssaMaxIterations);
    }

  } // end of the "outer outer loop" (index: l)

  done:

  if (getVerbosityLevel() > 2) {
    char tempstr[50] = "";
    snprintf(tempstr,50, "... =%5d outer iterations\n", outer_iterations);
    stdout_ssa += tempstr;
  } else if (getVerbosityLevel() == 2) {
    // at default verbosity, just record last normchange and iterations
    char tempstr[50] = "";
    snprintf(tempstr,50, "%5d outer iterations\n", outer_iterations);
    stdout_ssa += tempstr;
  }
  if (getVerbosityLevel() >= 2)
    stdout_ssa = "  SSA: " + stdout_ssa;

  if (write_ssa_system_to_matlab) {
    ierr = writeSSAsystemMatlab(); CHKERRQ(ierr);
  }

  return 0;
}


PetscErrorCode SSAFD::writeSSAsystemMatlab() {
  PetscErrorCode ierr;
  PetscViewer    viewer;
  char           file_name[PETSC_MAX_PATH_LEN], yearappend[PETSC_MAX_PATH_LEN];

  IceModelVec2S component;
  ierr = component.create(grid, "temp_storage", false); CHKERRQ(ierr);

  // FIXME: the file name prefix should be an option
  strcpy(file_name,"pism_SSAFD");
  snprintf(yearappend, PETSC_MAX_PATH_LEN, "_y%.0f.m", grid.year);
  strcat(file_name,yearappend);
  ierr = verbPrintf(2, grid.com,
             "writing Matlab-readable file for SSAFD system A xsoln = rhs to file `%s' ...\n",
             file_name); CHKERRQ(ierr);
  ierr = PetscViewerCreate(grid.com, &viewer);CHKERRQ(ierr);
  ierr = PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);CHKERRQ(ierr);
  ierr = PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);
  ierr = PetscViewerFileSetName(viewer, file_name);CHKERRQ(ierr);

  // get the command which started the run
  string cmdstr = pism_args_string();

  // save linear system; gives system A xsoln = rhs at last (nonlinear) iteration of SSA
  ierr = PetscViewerASCIIPrintf(viewer,
    "%% A PISM linear system report for the finite difference implementation\n"
    "%% of the SSA stress balance from this run:\n"
    "%%   '%s'\n"
    "%% Writes matrix A (sparse), and vectors uv and rhs, for the linear\n"
    "%% system which was solved at the last step of the nonlinear iteration:\n"
    "%%    A * uv = rhs.\n"
    "%% Also writes the year, the coordinates x,y, their gridded versions\n"
    "%% xx,yy, and the thickness (thk) and surface elevation (usurf).\n"
    "%% Also writes i-offsetvalues of vertically-integrated viscosity\n"
    "%% (nuH_0 = nu * H), and j-offset version of same thing (nuH_1 = nu * H);\n"
    "%% these are on the staggered grid.\n",
    cmdstr.c_str());  CHKERRQ(ierr);

  ierr = PetscViewerASCIIPrintf(viewer,"\n\necho off\n");  CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) SSAStiffnessMatrix,"A"); CHKERRQ(ierr);
  ierr = MatView(SSAStiffnessMatrix, viewer);CHKERRQ(ierr);
  ierr = PetscViewerASCIIPrintf(viewer,"clear zzz\n\n");  CHKERRQ(ierr);

  ierr = PetscObjectSetName((PetscObject) SSARHS,"rhs"); CHKERRQ(ierr);
  ierr = VecView(SSARHS, viewer);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject) SSAX,"uv"); CHKERRQ(ierr);
  ierr = VecView(SSAX, viewer);CHKERRQ(ierr);

  // save coordinates (for viewing, primarily)
  ierr = PetscViewerASCIIPrintf(viewer,"\nyear=%10.6f;\n",grid.year);  CHKERRQ(ierr);
  ierr = PetscViewerASCIIPrintf(viewer,
            "x=%12.3f + (0:%d)*%12.3f;\n"
            "y=%12.3f + (0:%d)*%12.3f;\n",
            -grid.Lx,grid.Mx-1,grid.dx,-grid.Ly,grid.My-1,grid.dy); CHKERRQ(ierr);
  ierr = PetscViewerASCIIPrintf(viewer,"[xx,yy]=meshgrid(x,y);\n");  CHKERRQ(ierr);

  // also save thickness and effective viscosity
  ierr = thickness->view_matlab(viewer); CHKERRQ(ierr);
  ierr = surface->view_matlab(viewer); CHKERRQ(ierr);

  ierr = nuH.get_component(0, component); CHKERRQ(ierr);
  ierr = component.set_name("nuH_0"); CHKERRQ(ierr);
  ierr = component.set_attr("long_name",
    "effective viscosity times thickness (i offset) at current time step"); CHKERRQ(ierr);
  ierr = component.view_matlab(viewer); CHKERRQ(ierr);

  ierr = nuH.get_component(0, component); CHKERRQ(ierr);
  ierr = component.set_name("nuH_1"); CHKERRQ(ierr);
  ierr = component.set_attr("long_name",
    "effective viscosity times thickness (j offset) at current time step"); CHKERRQ(ierr);
  ierr = component.view_matlab(viewer); CHKERRQ(ierr);

  ierr = PetscViewerASCIIPrintf(viewer,"echo on\n");  CHKERRQ(ierr);
  ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);

  return 0;
}


PetscErrorCode SSAFD::compute_nuH_norm(PetscReal &norm, PetscReal &norm_change) {
  PetscErrorCode ierr;

  PetscReal nuNorm[2], nuChange[2];

  const PetscScalar area = grid.dx * grid.dy;
#define MY_NORM     NORM_1

  // Test for change in nu
  ierr = nuH_old.add(-1, nuH); CHKERRQ(ierr);

  ierr = nuH_old.norm_all(MY_NORM, nuChange[0], nuChange[1]); CHKERRQ(ierr);
  ierr =     nuH.norm_all(MY_NORM, nuNorm[0],   nuNorm[1]);   CHKERRQ(ierr);

  nuChange[0] *= area;
  nuChange[1] *= area;
  nuNorm[0] *= area;
  nuNorm[1] *= area;

  norm_change = sqrt(PetscSqr(nuChange[0]) + PetscSqr(nuChange[1]));
  norm = sqrt(PetscSqr(nuNorm[0]) + PetscSqr(nuNorm[1]));
  return 0;
}

PetscErrorCode SSAFD::compute_hardav_staggered(IceModelVec2Stag &result) {
  PetscErrorCode ierr;
  PetscScalar *E, *E_ij, *E_offset;

  E = new PetscScalar[grid.Mz];

  ierr = thickness->begin_access(); CHKERRQ(ierr);
  ierr = enthalpy->begin_access(); CHKERRQ(ierr);
  ierr = result.begin_access(); CHKERRQ(ierr);

  for (PetscInt   i = grid.xs; i < grid.xs+grid.xm; ++i) {
    for (PetscInt j = grid.ys; j < grid.ys+grid.ym; ++j) {
      ierr = enthalpy->getInternalColumn(i,j,&E_ij); CHKERRQ(ierr);
      for (PetscInt o=0; o<2; o++) {
        const PetscInt oi = 1-o, oj=o;
        const PetscScalar H = 0.5 * ((*thickness)(i,j) + (*thickness)(i+oi,j+oj));

        if (H == 0) {
          result(i,j,o) = -1e6; // an obviously impossible value
          continue;
        }

        ierr = enthalpy->getInternalColumn(i+oi,j+oj,&E_offset); CHKERRQ(ierr);
        // build a column of enthalpy values a the current location:
        for (int k = 0; k < grid.Mz; ++k) {
          E[k] = 0.5 * (E_ij[k] + E_offset[k]);
        }

        result(i,j,o) = ice.averagedHardness_from_enth(H, grid.kBelowHeight(H),
                                                       grid.zlevels.data(), E); CHKERRQ(ierr);
      } // o
    }   // j
  }     // i

  ierr = result.end_access(); CHKERRQ(ierr);
  ierr = enthalpy->end_access(); CHKERRQ(ierr);
  ierr = thickness->end_access(); CHKERRQ(ierr);

  delete [] E;
  return 0;
}


PetscErrorCode SSAFD::compute_nuH_staggered(IceModelVec2Stag &result, PetscReal epsilon) {
  PetscErrorCode ierr;
  PISMVector2 **uv;

  ierr = result.begin_access(); CHKERRQ(ierr);
  ierr = velocity.get_array(uv); CHKERRQ(ierr);
  ierr = hardness.begin_access(); CHKERRQ(ierr);
  ierr = thickness->begin_access(); CHKERRQ(ierr);

  PetscScalar ssa_enhancement_factor=1.0;

  if (config.get_flag("do_ssa_enhancement"))
    ssa_enhancement_factor=config.get("ssa_enhancement_factor");

  const PetscScalar dx = grid.dx, dy = grid.dy;

  for (PetscInt o=0; o<2; ++o) {
    const PetscInt oi = 1 - o, oj=o;
    for (PetscInt i=grid.xs; i<grid.xs+grid.xm; ++i) {
      for (PetscInt j=grid.ys; j<grid.ys+grid.ym; ++j) {

        const PetscScalar H = 0.5 * ((*thickness)(i,j) + (*thickness)(i+oi,j+oj));

        if (H < strength_extension->get_min_thickness()) {
          // Extends strength of SSA (i.e. nuH coeff) into the ice free region.
          //  Does not add or subtract ice mass.
          result(i,j,o) = strength_extension->get_notional_strength();
          continue;
        }

        PetscScalar u_x, u_y, v_x, v_y;
        // Check the offset to determine how to differentiate velocity
        if (o == 0) {
          u_x = (uv[i+1][j].u - uv[i][j].u) / dx;
          u_y = (uv[i][j+1].u + uv[i+1][j+1].u - uv[i][j-1].u - uv[i+1][j-1].u) / (4*dy);
          v_x = (uv[i+1][j].v - uv[i][j].v) / dx;
          v_y = (uv[i][j+1].v + uv[i+1][j+1].v - uv[i][j-1].v - uv[i+1][j-1].v) / (4*dy);
        } else {
          u_x = (uv[i+1][j].u + uv[i+1][j+1].u - uv[i-1][j].u - uv[i-1][j+1].u) / (4*dx);
          u_y = (uv[i][j+1].u - uv[i][j].u) / dy;
          v_x = (uv[i+1][j].v + uv[i+1][j+1].v - uv[i-1][j].v - uv[i-1][j+1].v) / (4*dx);
          v_y = (uv[i][j+1].v - uv[i][j].v) / dy;
        }

        result(i,j,o) = H * ice.effectiveViscosity(hardness(i,j,o), u_x, u_y, v_x, v_y);

        if (! finite(result(i,j,o)) || false) {
          ierr = PetscPrintf(grid.com, "nuH[%d][%d][%d] = %e\n", o, i, j, result(i,j,o));
          CHKERRQ(ierr);
          ierr = PetscPrintf(grid.com, "  u_x, u_y, v_x, v_y = %e, %e, %e, %e\n",
                             u_x, u_y, v_x, v_y);
          CHKERRQ(ierr);
        }

        // We ensure that nuH is bounded below by a positive constant.
        result(i,j,o) += epsilon;

        // include the SSA enhancement factor; in most cases ssa_enhancement_factor is 1
        result(i,j,o) /= ssa_enhancement_factor;

      } // j
    } // i
  } // o

  ierr = thickness->end_access(); CHKERRQ(ierr);
  ierr = hardness.end_access(); CHKERRQ(ierr);
  ierr = result.end_access(); CHKERRQ(ierr);
  ierr = velocity.end_access(); CHKERRQ(ierr);

  // Some communication
  ierr = result.beginGhostComm(); CHKERRQ(ierr);
  ierr = result.endGhostComm(); CHKERRQ(ierr);
  return 0;
}


PetscErrorCode SSAFD::update_nuH_viewers() {
  PetscErrorCode ierr;
  IceModelVec2S tmp;

  if (!view_nuh) return 0;

  ierr = tmp.create(grid, "nuH", false); CHKERRQ(ierr);
  ierr = tmp.set_attrs("temporary",
                       "log10 of (viscosity * thickness)",
                       "Pa s m", ""); CHKERRQ(ierr);

  ierr = nuH.begin_access(); CHKERRQ(ierr);
  ierr = tmp.begin_access(); CHKERRQ(ierr);

  for (PetscInt   i = grid.xs; i < grid.xs+grid.xm; ++i) {
    for (PetscInt j = grid.ys; j < grid.ys+grid.ym; ++j) {
      PetscReal avg_nuH = 0.5 * (nuH(i,j,0) + nuH(i,j,1));
        if (avg_nuH > 1.0e14) {
          tmp(i,j) = log10(avg_nuH);
        } else {
          tmp(i,j) = 14.0;
        }
    }
  }

  ierr = tmp.end_access(); CHKERRQ(ierr);
  ierr = nuH.end_access(); CHKERRQ(ierr);

  if (nuh_viewer == PETSC_NULL) {
    ierr = grid.create_viewer(nuh_viewer_size, "nuH", nuh_viewer); CHKERRQ(ierr);
  }

  ierr = tmp.view(nuh_viewer, PETSC_NULL); CHKERRQ(ierr);

  return 0;
}

PetscErrorCode SSAFD::set_diagonal_matrix_entry(Mat A, int i, int j,
                                                PetscScalar value) {
  PetscErrorCode ierr;
  MatStencil row, col;
  row.j = i; row.i = j;
  col.j = i; col.i = j;

  row.c = 0; col.c = 0;
  ierr = MatSetValuesStencil(A, 1, &row, 1, &col, &value, INSERT_VALUES); CHKERRQ(ierr);

  row.c = 1; col.c = 1;
  ierr = MatSetValuesStencil(A, 1, &row, 1, &col, &value, INSERT_VALUES); CHKERRQ(ierr);

  return 0;
}


bool SSAFD::is_marginal(int i, int j) {
  const PetscInt M_ij = mask->as_int(i,j),
    // direct neighbors
    M_e = mask->as_int(i + 1,j),
    M_w = mask->as_int(i - 1,j),
    M_n = mask->as_int(i,j + 1),
    M_s = mask->as_int(i,j - 1),
    // "diagonal" neighbors
    M_ne = mask->as_int(i + 1,j + 1),
    M_se = mask->as_int(i + 1,j - 1),
    M_nw = mask->as_int(i - 1,j + 1),
    M_sw = mask->as_int(i - 1,j - 1);

  Mask M;
  
  return (!M.ice_free(M_ij)) &&
    (M.ice_free(M_e) || M.ice_free(M_w) || M.ice_free(M_n) || M.ice_free(M_s) ||
     M.ice_free(M_ne) || M.ice_free(M_se) || M.ice_free(M_nw) || M.ice_free(M_sw));
}