PISM, A Parallel Ice Sheet Model: src/base/iMcalving.cc Source File

PISM, A Parallel Ice Sheet Model stable v0.5
00001 // Copyright (C) 2004--2011 Torsten Albrecht and Constantine Khroulev
00002 //
00003 // This file is part of PISM.
00004 //
00005 // PISM is free software; you can redistribute it and/or modify it under the
00006 // terms of the GNU General Public License as published by the Free Software
00007 // Foundation; either version 2 of the License, or (at your option) any later
00008 // version.
00009 //
00010 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
00011 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00012 // FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
00013 // details.
00014 //
00015 // You should have received a copy of the GNU General Public License
00016 // along with PISM; if not, write to the Free Software
00017 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00018 
00019 
00020 #include <cmath>
00021 #include <petscdmda.h>
00022 #include "iceModel.hh"
00023 #include "pism_signal.h"
00024 #include "Mask.hh"
00025 #include "PISMOcean.hh"
00026 
00027 
00029 
00030 
00032 
00035 PetscErrorCode IceModel::eigenCalving() {
00036   const PetscScalar   dx = grid.dx, dy = grid.dy;
00037   PetscErrorCode ierr;
00038   ierr = verbPrintf(4, grid.com, "######### eigenCalving() start \n");    CHKERRQ(ierr);
00039   
00040   PetscScalar
00041     my_discharge_flux = 0,
00042     discharge_flux = 0;
00043 
00044   // is ghost communication really needed here?
00045   ierr = vH.beginGhostComm(); CHKERRQ(ierr);
00046   ierr = vH.endGhostComm(); CHKERRQ(ierr);
00047   ierr = vMask.beginGhostComm(); CHKERRQ(ierr);
00048   ierr = vMask.endGhostComm(); CHKERRQ(ierr);
00049 
00050   double ocean_rho = config.get("sea_water_density");
00051   double ice_rho = config.get("ice_density");
00052 
00053   const PetscScalar eigenCalvFactor = config.get("eigen_calving_K");
00054 
00055   PetscReal sea_level = 0;
00056 
00057   if (ocean != NULL) {
00058     ierr = ocean->sea_level_elevation(sea_level); CHKERRQ(ierr);
00059   } else { SETERRQ(grid.com, 2, "PISM ERROR: ocean == NULL"); }
00060 
00061   IceModelVec2S vHnew = vWork2d[0];
00062   ierr = vH.copy_to(vHnew); CHKERRQ(ierr);
00063 
00064   IceModelVec2S vDiffCalvRate = vWork2d[1];
00065   ierr = vDiffCalvRate.set(0.0); CHKERRQ(ierr);
00066 
00067   if (PetscAbs(dx - dy)/PetscMin(dx,dy) > 1e-2) {
00068     PetscPrintf(grid.com,
00069       "PISMPIK_ERROR: -eigen_calving using a non-square grid cell does not work (yet);\n"
00070       "               since it has no direction!!!\n, dx = %f, dy = %f, rel. diff = %f",
00071       dx,dy,PetscAbs(dx - dy)/PetscMax(dx,dy));
00072     PISMEnd();
00073   }
00074 
00075   // Distance (grid cells) from calving front where strain rate is evaluated
00076   PetscInt offset = 2;
00077 
00078   MaskQuery mask(vMask);
00079 
00080   ierr = vH.begin_access(); CHKERRQ(ierr);
00081   ierr = vHnew.begin_access(); CHKERRQ(ierr);
00082   ierr = vMask.begin_access(); CHKERRQ(ierr);
00083   ierr = vbed.begin_access(); CHKERRQ(ierr);
00084   ierr = vHref.begin_access(); CHKERRQ(ierr);
00085   ierr = vPrinStrain1.begin_access(); CHKERRQ(ierr);
00086   ierr = vPrinStrain2.begin_access(); CHKERRQ(ierr);
00087   ierr = vDiffCalvRate.begin_access(); CHKERRQ(ierr);
00088   for (PetscInt i = grid.xs; i < grid.xs + grid.xm; ++i) {
00089     for (PetscInt j = grid.ys; j < grid.ys + grid.ym; ++j) {
00090       // Average of strain-rate eigenvalues in adjacent floating grid cells to
00091       // be used for eigencalving
00092       PetscScalar eigen1 = 0.0, eigen2 = 0.0;
00093 
00094       // Number of directly adjacent floating boxes
00095       PetscInt N = 0;
00096 
00097       // Neighbor-averaged ice thickness
00098       PetscScalar H_average = 0.0; // is calculated here as average over direct neighbors
00099 
00100       // Counting adjacent floating boxes (with distance "offset")
00101       PetscInt M = 0;
00102 
00103       // find partially filled or empty grid boxes on the icefree ocean, which
00104       // have floating ice neighbors after massContExplicitStep (mask not updated)
00105 
00106       bool floating_e = (vH(i + 1, j) > 0.0 && (vbed(i + 1, j) < (sea_level - ice_rho / ocean_rho*vH(i + 1, j)))), 
00107            floating_w = (vH(i - 1, j) > 0.0 && (vbed(i - 1, j) < (sea_level - ice_rho / ocean_rho*vH(i - 1, j)))),
00108            floating_n = (vH(i, j + 1) > 0.0 && (vbed(i, j + 1) < (sea_level - ice_rho / ocean_rho*vH(i, j + 1)))),
00109            floating_s = (vH(i, j - 1) > 0.0 && (vbed(i, j - 1) < (sea_level - ice_rho / ocean_rho*vH(i, j - 1))));
00110 
00111       bool next_to_floating = floating_e || floating_w || floating_n || floating_s;
00112 
00113       bool ice_free_ocean = ( vH(i, j) == 0.0 && vbed(i, j) < sea_level );
00114 
00115       H_average = 0.0; eigen1 = 0.0; eigen2 = 0.0; M = 0, N = 0;
00116 
00117       if (ice_free_ocean && next_to_floating) {
00118 
00119         if ( floating_e ) { N += 1; H_average += vH(i + 1, j); }
00120         if ( floating_w ) { N += 1; H_average += vH(i - 1, j); }
00121         if ( floating_n ) { N += 1; H_average += vH(i, j + 1); }
00122         if ( floating_s ) { N += 1; H_average += vH(i, j - 1); }
00123 
00124         if (N > 0)
00125           H_average /= N;
00126 
00127         if ( mask.floating_ice(i + offset, j) && !mask.ice_margin(i + offset, j)){
00128           eigen1 += vPrinStrain1(i + offset, j);
00129           eigen2 += vPrinStrain2(i + offset, j);
00130           M += 1;
00131         }
00132 
00133         if ( mask.floating_ice(i - offset, j) && !mask.ice_margin(i - offset, j)){
00134           eigen1 += vPrinStrain1(i - offset, j);
00135           eigen2 += vPrinStrain2(i - offset, j);
00136           M += 1;
00137         }
00138 
00139         if ( mask.floating_ice(i, j + offset) && !mask.ice_margin(i , j + offset)){
00140           eigen1 += vPrinStrain1(i, j + offset);
00141           eigen2 += vPrinStrain2(i, j + offset);
00142           M += 1;
00143         }
00144 
00145         if ( mask.floating_ice(i, j - offset) && !mask.ice_margin(i , j - offset)){
00146           eigen1 += vPrinStrain1(i, j - offset);
00147           eigen2 += vPrinStrain2(i, j - offset);
00148           M += 1;
00149         }
00150 
00151         if (M > 0) {
00152           eigen1 /= M;
00153           eigen2 /= M;
00154         }
00155 
00156         PetscScalar calvrateHorizontal = 0.0,
00157           eigenCalvOffset = 0.0; // if it's not exactly the zero line of
00158                                // transition from compressive to extensive flow regime
00159 
00160         // calving law
00161         if ( eigen2 > eigenCalvOffset && eigen1 > 0.0) { // if spreading in all directions
00162           calvrateHorizontal = eigenCalvFactor * eigen1 * (eigen2 - eigenCalvOffset);
00163           // eigen1 * eigen2 has units [s^ - 2] and calvrateHorizontal [m*s^1]
00164           // hence, eigenCalvFactor has units [m*s]
00165         } else calvrateHorizontal = 0.0;
00166 
00167         // calculate mass loss with respect to the associated ice thickness and the grid size:
00168         PetscScalar calvrate = calvrateHorizontal * H_average / dx; // in m/s
00169 
00170         // apply calving rate at partially filled or empty grid cells
00171         if (calvrate > 0.0) {
00172           my_discharge_flux -= calvrate * dt; // no need to account for diffcalvrate further down, its all in this line
00173           vHref(i, j) -= calvrate * dt; // in m
00174           if(vHref(i, j) < 0.0) { // i.e. partially filled grid cell has completely calved off
00175             vDiffCalvRate(i, j) =  - vHref(i, j) / dt;// in m/s, means additional ice loss
00176             vHref(i, j) = 0.0;
00177             if(N > 0){
00178               vDiffCalvRate(i, j) = vDiffCalvRate(i, j) / N;
00179             }
00180           }
00181         }
00182       }
00183     }
00184   }
00185   ierr = vDiffCalvRate.end_access(); CHKERRQ(ierr);
00186 
00187   ierr = vDiffCalvRate.beginGhostComm(); CHKERRQ(ierr);
00188   ierr = vDiffCalvRate.endGhostComm(); CHKERRQ(ierr);
00189 
00190   ierr = vDiffCalvRate.begin_access(); CHKERRQ(ierr);
00191   for (PetscInt i = grid.xs; i < grid.xs + grid.xm; ++i) {
00192     for (PetscInt j = grid.ys; j < grid.ys + grid.ym; ++j) {
00193 
00194       PetscScalar restCalvRate = 0.0;
00195       bool hereFloating = (vH(i, j) > 0.0 && (vbed(i, j) < (sea_level - ice_rho / ocean_rho*vH(i, j))));
00196 
00197       if (hereFloating &&
00198           (vDiffCalvRate(i + 1, j) > 0.0 || vDiffCalvRate(i - 1, j) > 0.0 ||
00199            vDiffCalvRate(i, j + 1) > 0.0 || vDiffCalvRate(i, j - 1) > 0.0 )) {
00200 
00201         restCalvRate = (vDiffCalvRate(i + 1, j) +
00202                         vDiffCalvRate(i - 1, j) +
00203                         vDiffCalvRate(i, j + 1) +
00204                         vDiffCalvRate(i, j - 1));     // in m/s
00205 
00206         vHref(i, j) = vH(i, j) - (restCalvRate * dt); // in m
00207 
00208         vHnew(i, j) = 0.0;
00209 
00210         if(vHref(i, j) < 0.0) { // i.e. terminal floating ice grid cell has calved off completely.
00211           // We do not account for further calving ice-inwards!
00212           // Alternatively CFL criterion for time stepping could be adjusted to maximum of calving rate.
00213           // ierr = verbPrintf(2, grid.com, "!!!!! calving front would even retreat further at point %d, %d with volume %.2f \n",i,j,-vHref(i, j));    CHKERRQ(ierr);
00214           vHref(i, j) = 0.0;
00215         }
00216       }
00217     }
00218   }
00219   ierr = vHnew.end_access(); CHKERRQ(ierr);
00220   ierr = vH.end_access(); CHKERRQ(ierr);
00221   ierr = vMask.end_access(); CHKERRQ(ierr);
00222   ierr = vbed.end_access(); CHKERRQ(ierr);
00223   ierr = vHref.end_access(); CHKERRQ(ierr);
00224   ierr = vPrinStrain1.end_access(); CHKERRQ(ierr);
00225   ierr = vPrinStrain2.end_access(); CHKERRQ(ierr);
00226   ierr = vDiffCalvRate.end_access(); CHKERRQ(ierr);
00227   
00228   ierr = PISMGlobalSum(&my_discharge_flux,     &discharge_flux,     grid.com); CHKERRQ(ierr);
00229   PetscScalar factor = config.get("ice_density") * (dx * dy);
00230   cumulative_discharge_flux     += discharge_flux     * factor;
00231 
00232   ierr = vHnew.beginGhostComm(vH); CHKERRQ(ierr);
00233   ierr = vHnew.endGhostComm(vH); CHKERRQ(ierr);
00234 
00235   return 0;
00236 }
00237 
00238 
00245 PetscErrorCode IceModel::calvingAtThickness() {
00246   PetscErrorCode ierr;
00247   ierr = verbPrintf(4, grid.com, "######### calvingAtThickness() start\n");    CHKERRQ(ierr);
00248   
00249   PetscScalar
00250     my_discharge_flux = 0,
00251     discharge_flux = 0;
00252   const PetscScalar dx = grid.dx, dy = grid.dy;
00253 
00254   ierr = vH.beginGhostComm(); CHKERRQ(ierr);
00255   ierr = vH.endGhostComm(); CHKERRQ(ierr);
00256 
00257   IceModelVec2S vHnew = vWork2d[0];
00258   ierr = vH.copy_to(vHnew); CHKERRQ(ierr);
00259 
00260   double ocean_rho = config.get("sea_water_density"),
00261          ice_rho   = config.get("ice_density"),
00262          Hcalving  = config.get("calving_at_thickness");
00263 
00264   PetscReal sea_level;
00265   if (ocean != NULL) {
00266     ierr = ocean->sea_level_elevation(sea_level); CHKERRQ(ierr);
00267   } else { SETERRQ(grid.com, 2, "PISM ERROR: ocean == NULL"); }
00268 
00269   ierr = vH.begin_access(); CHKERRQ(ierr);
00270   ierr = vHnew.begin_access(); CHKERRQ(ierr);
00271   ierr = vbed.begin_access(); CHKERRQ(ierr);
00272   for (PetscInt i = grid.xs; i < grid.xs + grid.xm; ++i) {
00273     for (PetscInt j = grid.ys; j < grid.ys + grid.ym; ++j) {
00274       bool hereFloating = (vH(i, j) > 0.0 && (vbed(i, j) < (sea_level - ice_rho / ocean_rho * vH(i, j))));
00275       bool icefreeOceanNeighbor = ( (vH(i + 1, j) == 0.0 && vbed(i + 1, j) < sea_level) ||
00276                                     (vH(i - 1, j) == 0.0 && vbed(i - 1, j) < sea_level) ||
00277                                     (vH(i, j + 1) == 0.0 && vbed(i, j + 1) < sea_level) ||
00278                                     (vH(i, j - 1) == 0.0 && vbed(i, j - 1) < sea_level));
00279       if (hereFloating && vH(i, j) <= Hcalving && icefreeOceanNeighbor) {
00280         my_discharge_flux -= vHnew(i, j);
00281         vHnew(i, j) = 0.0;
00282       }
00283     }
00284   }
00285   ierr = vHnew.end_access(); CHKERRQ(ierr);
00286   ierr = vH.end_access(); CHKERRQ(ierr);
00287   ierr = vbed.end_access(); CHKERRQ(ierr);
00288   
00289   ierr = PISMGlobalSum(&my_discharge_flux,     &discharge_flux,     grid.com); CHKERRQ(ierr);
00290   PetscScalar factor = config.get("ice_density") * (dx * dy);
00291   cumulative_discharge_flux     += discharge_flux     * factor;
00292 
00293   ierr = vHnew.beginGhostComm(vH); CHKERRQ(ierr);
00294   ierr = vHnew.endGhostComm(vH); CHKERRQ(ierr);
00295 
00296   return 0;
00297 }
00298 
00299 
00300 
00302 
00303 PetscErrorCode IceModel::dt_from_eigenCalving() {
00304   const PetscScalar   dx = grid.dx, dy = grid.dy;
00305   PetscErrorCode ierr;
00306   ierr = verbPrintf(4, grid.com, "######### dt_from_eigenCalving() start \n");    CHKERRQ(ierr);
00307 
00308   // is ghost communication really needed here?
00309   ierr = vH.beginGhostComm(); CHKERRQ(ierr);
00310   ierr = vH.endGhostComm(); CHKERRQ(ierr);
00311   ierr = vMask.beginGhostComm(); CHKERRQ(ierr);
00312   ierr = vMask.endGhostComm(); CHKERRQ(ierr);
00313 
00314   double ocean_rho = config.get("sea_water_density");
00315   double ice_rho = config.get("ice_density");
00316 
00317   const PetscScalar eigenCalvFactor = config.get("eigen_calving_K");
00318 
00319   PetscReal sea_level = 0;
00320 
00321   if (ocean != NULL) {
00322     ierr = ocean->sea_level_elevation(sea_level); CHKERRQ(ierr);
00323   } else { SETERRQ(grid.com, 2, "PISM ERROR: ocean == NULL"); }
00324 
00325   if (PetscAbs(dx - dy)/PetscMin(dx,dy) > 1e-2) {
00326     PetscPrintf(grid.com,
00327       "PISMPIK_ERROR: -eigen_calving using a non-square grid cell does not work (yet);\n"
00328       "               since it has no direction!!!\n, dx = %f, dy = %f, rel. diff = %f",
00329       dx,dy,PetscAbs(dx - dy)/PetscMax(dx,dy));
00330     PISMEnd();
00331   }
00332 
00333   PetscScalar dt_from_eigencalving_min = 0.001; //about 9 hours which corresponds to 10000 km/a on a 10 km grid
00334 
00335   // Distance (grid cells) from calving front where strain rate is evaluated
00336   PetscInt offset = 2;
00337   PetscScalar my_maxCalvingRate=0.0, my_meancalvrate=0.0, my_cratecounter=0.0;
00338   PetscInt i0=0, j0=0;
00339 
00340   MaskQuery mask(vMask);
00341 
00342   ierr = vH.begin_access(); CHKERRQ(ierr);
00343   ierr = vMask.begin_access(); CHKERRQ(ierr);
00344   ierr = vbed.begin_access(); CHKERRQ(ierr);
00345   ierr = vPrinStrain1.begin_access(); CHKERRQ(ierr);
00346   ierr = vPrinStrain2.begin_access(); CHKERRQ(ierr);
00347 
00348   for (PetscInt i = grid.xs; i < grid.xs + grid.xm; ++i) {
00349     for (PetscInt j = grid.ys; j < grid.ys + grid.ym; ++j) {
00350       // Average of strain-rate eigenvalues in adjacent floating grid cells to
00351       // be used for eigencalving
00352       PetscScalar eigen1 = 0.0, eigen2 = 0.0;
00353       // Neighbor-averaged ice thickness
00354       PetscInt M = 0;
00355 
00356       // find partially filled or empty grid boxes on the icefree ocean, which
00357       // have floating ice neighbors after massContExplicitStep (mask not updated)
00358       bool next_to_floating =
00359         ((vH(i + 1, j) > 0.0 && (vbed(i + 1, j) < (sea_level - ice_rho / ocean_rho*vH(i + 1, j)))) ||
00360          (vH(i - 1, j) > 0.0 && (vbed(i - 1, j) < (sea_level - ice_rho / ocean_rho*vH(i - 1, j)))) ||
00361          (vH(i, j + 1) > 0.0 && (vbed(i, j + 1) < (sea_level - ice_rho / ocean_rho*vH(i, j + 1)))) ||
00362          (vH(i, j - 1) > 0.0 && (vbed(i, j - 1) < (sea_level - ice_rho / ocean_rho*vH(i, j - 1)))));
00363 
00364       bool ice_free_ocean = ( vH(i, j) == 0.0 && vbed(i, j) < sea_level );
00365 
00366       if (ice_free_ocean && next_to_floating) {
00367 
00368           PetscScalar calvrateHorizontal = 0.0,
00369                       eigenCalvOffset = 0.0; 
00370 
00371           if ( mask.floating_ice(i + offset, j) && !mask.ice_margin(i + offset, j)) {
00372             eigen1 += vPrinStrain1(i + offset, j);
00373             eigen2 += vPrinStrain2(i + offset, j);
00374             M += 1;
00375           }
00376           if ( mask.floating_ice(i - offset, j) && !mask.ice_margin(i - offset, j)){
00377             eigen1 += vPrinStrain1(i - offset, j);
00378             eigen2 += vPrinStrain2(i - offset, j);
00379             M += 1;
00380           }
00381           if ( mask.floating_ice(i, j + offset) && !mask.ice_margin(i , j + offset)){
00382             eigen1 += vPrinStrain1(i, j + offset);
00383             eigen2 += vPrinStrain2(i, j + offset);
00384             M += 1;
00385           }
00386           if ( mask.floating_ice(i, j - offset) && !mask.ice_margin(i , j - offset)){
00387             eigen1 += vPrinStrain1(i, j - offset);
00388             eigen2 += vPrinStrain2(i, j - offset);
00389             M += 1;
00390           }
00391           if (M > 0) {
00392             eigen1 /= M;
00393             eigen2 /= M;
00394           }
00395 
00396           // calving law
00397           if ( eigen2 > eigenCalvOffset && eigen1 > 0.0) { // if spreading in all directions
00398             calvrateHorizontal = eigenCalvFactor * eigen1 * (eigen2 - eigenCalvOffset);
00399             my_cratecounter+=1.0;
00400             my_meancalvrate+=calvrateHorizontal;
00401             if ( my_maxCalvingRate < calvrateHorizontal) {
00402               i0=i;
00403               j0=j;
00404             }
00405             my_maxCalvingRate=PetscMax(my_maxCalvingRate,calvrateHorizontal);
00406           } else calvrateHorizontal = 0.0;
00407        }
00408     }
00409   }
00410 
00411   ierr = vH.end_access(); CHKERRQ(ierr);
00412   ierr = vMask.end_access(); CHKERRQ(ierr);
00413   ierr = vbed.end_access(); CHKERRQ(ierr);
00414   ierr = vPrinStrain1.end_access(); CHKERRQ(ierr);
00415   ierr = vPrinStrain2.end_access(); CHKERRQ(ierr);
00416 
00417   PetscScalar maxCalvingRate=0.0, meancalvrate=0.0, cratecounter=0.0;
00418   ierr = PISMGlobalSum(&my_meancalvrate, &meancalvrate, grid.com); CHKERRQ(ierr);
00419   ierr = PISMGlobalSum(&my_cratecounter, &cratecounter, grid.com); CHKERRQ(ierr);
00420   meancalvrate/=cratecounter;
00421   ierr = PISMGlobalMax(&my_maxCalvingRate, &maxCalvingRate, grid.com); CHKERRQ(ierr);
00422 
00423   PetscScalar denom = maxCalvingRate/dx;
00424   denom += (0.01/secpera)/(dx + dy);  // make sure it's pos.
00425   dt_from_eigencalving = 1.0/denom;
00426  
00427   ierr = verbPrintf(2, grid.com, "!!!!! c_rate = %.0f m/a ( dt=%.5f a ) at point %d, %d with mean_c=%.0f m/a over %.0f cells \n",maxCalvingRate*secpera,dt_from_eigencalving/secpera,i0,j0,meancalvrate*secpera,cratecounter);    CHKERRQ(ierr);
00428          
00429   dt_from_eigencalving = PetscMax(dt_from_eigencalving,dt_from_eigencalving_min);
00430         
00431   return 0;
00432 }
00433