stable0.3/src/base/iMgeometry.cc

Go to the documentation of this file.

// Copyright (C) 2004-2010 Jed Brown, Ed Bueler and Constantine Khroulev
//
// 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 <cmath>
#include <cstring>
#include <petscda.h>
#include "iceModel.hh"



PetscErrorCode IceModel::computeDrivingStress(IceModelVec2S &vtaudx, IceModelVec2S &vtaudy) {
  PetscErrorCode ierr;

  const PetscScalar n       = ice->exponent(), // frequently n = 3
                    etapow  = (2.0 * n + 2.0)/n,  // = 8/3 if n = 3
                    invpow  = 1.0 / etapow,  // = 3/8
                    dinvpow = (- n - 2.0) / (2.0 * n + 2.0); // = -5/8
  const PetscScalar minThickEtaTransform = 5.0; // m
  const PetscScalar dx=grid.dx, dy=grid.dy;

  bool compute_surf_grad_inward_ssa = config.get_flag("compute_surf_grad_inward_ssa");
  bool use_eta = config.get_flag("use_eta_transformation");

  ierr =    vh.begin_access();    CHKERRQ(ierr);
  ierr =    vH.begin_access();  CHKERRQ(ierr);
  ierr =  vbed.begin_access();  CHKERRQ(ierr);
  ierr = vMask.begin_access();  CHKERRQ(ierr);

  ierr = vtaudx.begin_access(); CHKERRQ(ierr);
  ierr = vtaudy.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) {
      const PetscScalar pressure = ice->rho * standard_gravity * vH(i,j);
      if (pressure <= 0.0) {
        vtaudx(i,j) = 0.0;
        vtaudy(i,j) = 0.0;
      } else {
        PetscScalar h_x = 0.0, h_y = 0.0;
        // FIXME: we need to handle grid periodicity correctly.
        if (vMask.is_grounded(i,j) && (use_eta == true)) {
          // in grounded case, differentiate eta = H^{8/3} by chain rule
          if (vH(i,j) > 0.0) {
            const PetscScalar myH = (vH(i,j) < minThickEtaTransform) ?
              minThickEtaTransform : vH(i,j);
            const PetscScalar eta = pow(myH, etapow), factor = invpow * pow(eta, dinvpow);
            h_x = factor * (pow(vH(i+1,j),etapow) - pow(vH(i-1,j),etapow)) / (2*dx);
            h_y = factor * (pow(vH(i,j+1),etapow) - pow(vH(i,j-1),etapow)) / (2*dy);
          }
          // now add bed slope to get actual h_x,h_y
          // FIXME: there is no reason to assume user's bed is periodized; see vertical
          //   velocity computation
          h_x += vbed.diff_x(i,j);
          h_y += vbed.diff_y(i,j);
        } else {  // floating or eta transformation is not used
          if (compute_surf_grad_inward_ssa) {
            h_x = vh.diff_x_p(i,j);
            h_y = vh.diff_y_p(i,j);
          } else {
            h_x = vh.diff_x(i,j);
            h_y = vh.diff_y(i,j);
          }
        }

        vtaudx(i,j) = - pressure * h_x;
        vtaudy(i,j) = - pressure * h_y;
      }
    }
  }

  ierr =   vbed.end_access(); CHKERRQ(ierr);
  ierr =     vh.end_access(); CHKERRQ(ierr);
  ierr =     vH.end_access(); CHKERRQ(ierr);
  ierr =  vMask.end_access(); CHKERRQ(ierr);
  ierr = vtaudx.end_access(); CHKERRQ(ierr);
  ierr = vtaudy.end_access(); CHKERRQ(ierr);

  return 0;
}



PetscErrorCode IceModel::updateSurfaceElevationAndMask() {
  PetscErrorCode ierr;

  ierr = update_mask(); CHKERRQ(ierr);
  ierr = update_surface_elevation(); CHKERRQ(ierr);

  return 0;
}

PetscErrorCode IceModel::update_mask() {
  PetscErrorCode ierr;
  PetscScalar **bed, **H, **mask;

  if (ocean == PETSC_NULL) {  SETERRQ(1,"PISM ERROR: ocean == PETSC_NULL");  }
  PetscReal currentSeaLevel;
  ierr = ocean->sea_level_elevation(grid.year, dt / secpera, currentSeaLevel); CHKERRQ(ierr);

  bool use_ssa_when_grounded = config.get_flag("use_ssa_when_grounded"),
    use_ssa_velocity = config.get_flag("use_ssa_velocity");

  double ocean_rho = config.get("sea_water_density");

  ierr =    vH.get_array(H);    CHKERRQ(ierr);
  ierr =  vbed.get_array(bed);  CHKERRQ(ierr);
  ierr = vMask.get_array(mask); CHKERRQ(ierr);

#ifdef LOCAL_GHOST_UPDATE
  PetscInt GHOSTS = 2;
#else
  PetscInt GHOSTS = 0;
#endif
  for (PetscInt   i = grid.xs - GHOSTS; i < grid.xs+grid.xm + GHOSTS; ++i) {
    for (PetscInt j = grid.ys - GHOSTS; j < grid.ys+grid.ym + GHOSTS; ++j) {

      const PetscScalar hgrounded = bed[i][j] + vH(i,j),
        hfloating = currentSeaLevel + (1.0 - ice->rho/ocean_rho) * vH(i,j);

      // points marked as "ocean at time zero" are not updated
      if (vMask.value(i,j) == MASK_OCEAN_AT_TIME_0)
        continue;

      if (vMask.is_floating(i,j)) { // floating

        if (hgrounded > hfloating+1.0) { // flotation criterion says it is grounded
          if (use_ssa_velocity) {
            if (use_ssa_when_grounded) {
              mask[i][j] = MASK_DRAGGING_SHEET;
            } else {
              mask[i][j] = MASK_SHEET;
            }
          } else {
            // we do not have any ice handled by SSA, so it must be SHEET
            mask[i][j] = MASK_SHEET;
          }
        }

      } else {   // grounded

        // apply the flotation criterion:
        if (hgrounded > hfloating-1.0) { // flotation criterion says it is grounded

          // we are using SSA-as-a-sliding-law, so grounded points become DRAGGING
          if (use_ssa_velocity && use_ssa_when_grounded)
            mask[i][j] = MASK_DRAGGING_SHEET;

        } else {
          mask[i][j] = MASK_FLOATING;
        }

      }
    } // inner for loop (j)
  } // outer for loop (i)

  ierr =         vH.end_access(); CHKERRQ(ierr);
  ierr =       vbed.end_access(); CHKERRQ(ierr);
  ierr =      vMask.end_access(); CHKERRQ(ierr);

#ifndef LOCAL_GHOST_UPDATE
  ierr = vMask.beginGhostComm(); CHKERRQ(ierr);
  ierr = vMask.endGhostComm(); CHKERRQ(ierr);
#endif

  return 0;
}

PetscErrorCode IceModel::update_surface_elevation() {
  PetscErrorCode ierr;
  PetscScalar **h, **bed, **H, **mask;

  if (ocean == PETSC_NULL) {  SETERRQ(1,"PISM ERROR: ocean == PETSC_NULL");  }
  PetscReal currentSeaLevel;
  ierr = ocean->sea_level_elevation(grid.year, dt / secpera, currentSeaLevel); CHKERRQ(ierr);

  bool is_dry_simulation = config.get_flag("is_dry_simulation");

  double ocean_rho = config.get("sea_water_density");

  ierr =    vh.get_array(h);    CHKERRQ(ierr);
  ierr =    vH.get_array(H);    CHKERRQ(ierr);
  ierr =  vbed.get_array(bed);  CHKERRQ(ierr);
  ierr = vMask.get_array(mask); CHKERRQ(ierr);

#ifdef LOCAL_GHOST_UPDATE
  PetscInt GHOSTS = 2;
#else
  PetscInt GHOSTS = 0;
#endif
  for (PetscInt   i = grid.xs - GHOSTS; i < grid.xs+grid.xm + GHOSTS; ++i) {
    for (PetscInt j = grid.ys - GHOSTS; j < grid.ys+grid.ym + GHOSTS; ++j) {
      // take this opportunity to check that vH(i,j) >= 0
      if (vH(i,j) < 0) {
        SETERRQ2(1,"Thickness negative at point i=%d, j=%d",i,j);
      }

      const PetscScalar hgrounded = bed[i][j] + vH(i,j),
        hfloating = currentSeaLevel + (1.0 - ice->rho/ocean_rho) * vH(i,j);

      if (is_dry_simulation) {
        // Don't update mask; potentially one would want to do SSA dragging ice
        //   shelf in dry case and/or ignore mean sea level elevation.
        h[i][j] = hgrounded;
        continue;               // go to the next grid point
      }

      if (vMask.value(i,j) == MASK_OCEAN_AT_TIME_0) {
        // mask takes priority over bed in this case (note sea level may change).
        // Example Greenland case: if mask say Ellesmere is OCEAN0,
        //   then never want ice on Ellesmere.
        // If mask says OCEAN0 then don't change the mask and also don't change
        // the thickness; massContExplicitStep() is in charge of that.
        // Almost always the next line is equivalent to h[i][j] = 0.
        h[i][j] = hfloating;  // ignore bed and treat it like deep ocean
        continue;             // go to the next grid point
      }

      if (vMask.is_floating(i,j)) {
        h[i][j] = hfloating; // actually floating so update h
      } else { 
        h[i][j] = hgrounded; // actually grounded so update h
      }
    }
        
  }

  ierr =         vh.end_access(); CHKERRQ(ierr);
  ierr =         vH.end_access(); CHKERRQ(ierr);
  ierr =       vbed.end_access(); CHKERRQ(ierr);
  ierr =      vMask.end_access(); CHKERRQ(ierr);

#ifndef LOCAL_GHOST_UPDATE
  ierr = vh.beginGhostComm(); CHKERRQ(ierr);
  ierr = vh.endGhostComm(); CHKERRQ(ierr);
#endif

  return 0;
}


PetscErrorCode IceModel::massContExplicitStep() {
  PetscErrorCode ierr;

  const PetscScalar   dx = grid.dx, dy = grid.dy;
  bool do_ocean_kill = config.get_flag("ocean_kill"),
    floating_ice_killed = config.get_flag("floating_ice_killed"),
    include_bmr_in_continuity = config.get_flag("include_bmr_in_continuity"),
    do_superpose = config.get_flag("do_superpose");

  if (surface != NULL) {
    ierr = surface->ice_surface_mass_flux(grid.year, dt / secpera, acab); CHKERRQ(ierr);
  } else { SETERRQ(1,"PISM ERROR: surface == NULL"); }

  if (ocean != NULL) {
    ierr = ocean->shelf_base_mass_flux(grid.year, dt / secpera, shelfbmassflux); CHKERRQ(ierr);
  } else { SETERRQ(2,"PISM ERROR: ocean == NULL"); }

  const PetscScalar inC_fofv = 1.0e-4 * PetscSqr(secpera),
                    outC_fofv = 2.0 / pi;

  IceModelVec2S vHnew = vWork2d[0];
  ierr = vH.copy_to(vHnew); CHKERRQ(ierr);

  PetscScalar **bmr_gnded;

  ierr = vH.begin_access(); CHKERRQ(ierr);
  ierr = vbmr.get_array(bmr_gnded); CHKERRQ(ierr);
  ierr = uvbar.begin_access(); CHKERRQ(ierr);
  ierr = vel_basal.begin_access(); CHKERRQ(ierr);
  ierr = acab.begin_access(); CHKERRQ(ierr);
  ierr = shelfbmassflux.begin_access(); CHKERRQ(ierr);
  ierr = vMask.begin_access();  CHKERRQ(ierr);
  ierr = vHnew.begin_access(); CHKERRQ(ierr);

  ierr = vel_ssa.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) {

      // get thickness averaged onto staggered grid;
      //    note Div Q = Div (- f(v) D grad h + (1-f(v)) U_b H) 
      //    in  -ssa -super case; note f(v) is on regular grid;
      //    compare broadcastSSAVelocity(); note uvbar[o] is SIA result:
      //    uvbar[0] H = - D h_x
      PetscScalar He, Hw, Hn, Hs;
      if ( do_superpose && (vMask.value(i,j) == MASK_DRAGGING_SHEET) ) {
        const PetscScalar
          fv  = 1.0 - outC_fofv * atan( inC_fofv *
                      ( PetscSqr(vel_ssa(i,  j).u) + PetscSqr(vel_ssa(i,  j).v) ) ),
          fve = 1.0 - outC_fofv * atan( inC_fofv *
                      ( PetscSqr(vel_ssa(i+1,j).u) + PetscSqr(vel_ssa(i+1,j).v) ) ),
          fvw = 1.0 - outC_fofv * atan( inC_fofv *
                      ( PetscSqr(vel_ssa(i-1,j).u) + PetscSqr(vel_ssa(i-1,j).v) ) ),
          fvn = 1.0 - outC_fofv * atan( inC_fofv *
                      ( PetscSqr(vel_ssa(i,j+1).u) + PetscSqr(vel_ssa(i,j+1).v) ) ),
          fvs = 1.0 - outC_fofv * atan( inC_fofv *
                      ( PetscSqr(vel_ssa(i,j-1).u) + PetscSqr(vel_ssa(i,j-1).v) ) );
        const PetscScalar fvH = fv * vH(i,j);
        He = 0.5 * (fvH + fve * vH(i+1,j));
        Hw = 0.5 * (fvw * vH(i-1,j) + fvH);
        Hn = 0.5 * (fvH + fvn * vH(i,j+1));
        Hs = 0.5 * (fvs * vH(i,j-1) + fvH);
      } else {
        He = 0.5 * (vH(i,j) + vH(i+1,j));
        Hw = 0.5 * (vH(i-1,j) + vH(i,j));
        Hn = 0.5 * (vH(i,j) + vH(i,j+1));
        Hs = 0.5 * (vH(i,j-1) + vH(i,j));
      }

      // staggered grid Div(Q) for SIA (non-sliding) deformation part;
      //    Q = - D grad h = Ubar H    in non-sliding case
      PetscScalar divQ = 0.0;
      if (computeSIAVelocities == PETSC_TRUE) {
        divQ =  (uvbar(i,j,0) * He - uvbar(i-1,j,0) * Hw) / dx
              + (uvbar(i,j,1) * Hn - uvbar(i,j-1,1) * Hs) / dy;
      }

      // basal sliding part: split  Div(v H)  by product rule into  v . grad H
      //    and  (Div v) H; use upwinding on first and centered on second
      divQ +=  vel_basal(i,j).u * ( vel_basal(i,j).u < 0 ? vH(i+1,j)-vH(i,j)
                                 : vH(i,j)-vH(i-1,j) ) / dx
             + vel_basal(i,j).v * ( vel_basal(i,j).v < 0 ? vH(i,j+1)-vH(i,j)
                                 : vH(i,j)-vH(i,j-1) ) / dy;

      divQ += vH(i,j) * ( (vel_basal(i+1,j).u - vel_basal(i-1,j).u) / (2.0*dx)
                          + (vel_basal(i,j+1).v - vel_basal(i,j-1).v) / (2.0*dy) );

      vHnew(i,j) += (acab(i,j) - divQ) * dt; // include M

      if (include_bmr_in_continuity) { // include S
        if (vMask.is_floating(i,j)) {
          vHnew(i,j) -= shelfbmassflux(i,j) * dt;
        } else {
           vHnew(i,j) -= bmr_gnded[i][j] * dt;
        }
      }

      // apply free boundary rule: negative thickness becomes zero
      if (vHnew(i,j) < 0)
        vHnew(i,j) = 0.0;

      // the following conditionals, both -ocean_kill and -float_kill, are also applied in 
      //   IceModel::computeMax2DSlidingSpeed() when determining CFL
      
      // force zero thickness at points which were originally ocean (if "-ocean_kill");
      //   this is calving at original calving front location
      if ( do_ocean_kill && (vMask.value(i,j) == MASK_OCEAN_AT_TIME_0) )
        vHnew(i,j) = 0.0;

      // force zero thickness at points which are floating (if "-float_kill");
      //   this is calving at grounding line
      if ( floating_ice_killed && vMask.is_floating(i,j) )
        vHnew(i,j) = 0.0;

    } // end of the inner for loop
  } // end of the outer for loop

  ierr = vbmr.end_access(); CHKERRQ(ierr);
  ierr = vMask.end_access(); CHKERRQ(ierr);
  ierr = uvbar.end_access(); CHKERRQ(ierr);
  ierr = vel_basal.end_access(); CHKERRQ(ierr);
  ierr = vel_ssa.end_access(); CHKERRQ(ierr);
  ierr = acab.end_access(); CHKERRQ(ierr);
  ierr = shelfbmassflux.end_access(); CHKERRQ(ierr);
  ierr = vH.end_access(); CHKERRQ(ierr);
  ierr = vHnew.end_access(); CHKERRQ(ierr);

  // compute dH/dt (thickening rate) for viewing and for saving at end; only diagnostic
  ierr = vHnew.add(-1.0, vH, vdHdt); CHKERRQ(ierr); // vdHdt = vHnew - vH
  ierr = vdHdt.scale(1.0/dt); CHKERRQ(ierr);        // vdHdt = vdHdt / dt

  // d(volume)/dt
  {
    PetscScalar dvol=0.0;
  
    ierr = vdHdt.begin_access(); CHKERRQ(ierr);
    
    if (cell_area.was_created()) {
      ierr = cell_area.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) {
          dvol += vdHdt(i,j) * cell_area(i,j);
        }
      }  
      ierr = cell_area.end_access(); CHKERRQ(ierr);
    } else {
      const PetscScalar a = grid.dx * grid.dy; // cell area
      for (PetscInt i=grid.xs; i<grid.xs+grid.xm; ++i) {
        for (PetscInt j=grid.ys; j<grid.ys+grid.ym; ++j) {
          dvol += vdHdt(i,j) * a;
        }
      }  
    }

    ierr = vdHdt.end_access(); CHKERRQ(ierr);
    ierr = PetscGlobalSum(&dvol, &dvoldt, grid.com); CHKERRQ(ierr);
  }

  // average value of dH/dt; 
  PetscScalar ice_area;
  ierr = compute_ice_area(ice_area); CHKERRQ(ierr);

  ierr = vdHdt.sum(gdHdtav); CHKERRQ(ierr);
  gdHdtav = gdHdtav / ice_area; // m/s


  // finally copy vHnew into vH and communicate ghosted values
  ierr = vHnew.beginGhostComm(vH); CHKERRQ(ierr);
  ierr = vHnew.endGhostComm(vH); CHKERRQ(ierr);

  // Check if the ice thickness is exceeded the height of the computational box
  // and extend the grid if necessary:
  ierr = check_maximum_thickness(); CHKERRQ(ierr);

  return 0;
}

---