|
PISM, A Parallel Ice Sheet Model stable 0.4.1779
|
src/base/stressbalance/SIA_Sliding.cc
Go to the documentation of this file.
00001 // Copyright (C) 2004--2011 Jed Brown, Ed Bueler 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 #include "SIA_Sliding.hh" 00020 #include "Mask.hh" 00021 00022 PetscErrorCode SIA_Sliding::allocate() { 00023 PetscErrorCode ierr; 00024 PetscInt WIDE_STENCIL = grid.max_stencil_width; 00025 00026 for (int i = 0; i < 2; ++i) { 00027 char namestr[30]; 00028 00029 ierr = work_2d_stag[i].create(grid, "work_vector", true); CHKERRQ(ierr); 00030 snprintf(namestr, sizeof(namestr), "work_vector_2d_stag_%d", i); 00031 ierr = work_2d_stag[i].set_name(namestr); CHKERRQ(ierr); 00032 00033 } 00034 00035 ierr = work_2d.create(grid, "work_vector_2d", true, WIDE_STENCIL); CHKERRQ(ierr); 00036 00037 return 0; 00038 } 00039 00040 PetscErrorCode SIA_Sliding::init(PISMVars &vars) { 00041 PetscErrorCode ierr; 00042 00043 ierr = ShallowStressBalance::init(vars); CHKERRQ(ierr); 00044 00045 standard_gravity = config.get("standard_gravity"); 00046 verification_mode = config.get_flag("verification_mode"); 00047 00048 if (config.has("EISMINT_II_experiment")) 00049 eisII_experiment = config.get_string("EISMINT_II_experiment"); 00050 00051 thickness = dynamic_cast<IceModelVec2S*>(vars.get("land_ice_thickness")); 00052 if (thickness == NULL) SETERRQ(1, "land_ice_thickness is not available"); 00053 00054 mask = dynamic_cast<IceModelVec2Int*>(vars.get("mask")); 00055 if (mask == NULL) SETERRQ(1, "mask is not available"); 00056 00057 surface = dynamic_cast<IceModelVec2S*>(vars.get("surface_altitude")); 00058 if (surface == NULL) SETERRQ(1, "surface_altitude is not available"); 00059 00060 bed = dynamic_cast<IceModelVec2S*>(vars.get("bedrock_altitude")); 00061 if (bed == NULL) SETERRQ(1, "bedrock_altitude is not available"); 00062 00063 enthalpy = dynamic_cast<IceModelVec3*>(vars.get("enthalpy")); 00064 if (enthalpy == NULL) SETERRQ(1, "enthalpy is not available"); 00065 00066 return 0; 00067 } 00068 00070 00082 PetscErrorCode SIA_Sliding::update(bool /*fast*/) { 00083 PetscErrorCode ierr; 00084 IceModelVec2Stag h_x = work_2d_stag[0], h_y = work_2d_stag[1]; 00085 00086 ierr = compute_surface_gradient(h_x, h_y); CHKERRQ(ierr); 00087 00088 ierr = D2.set(0.0); CHKERRQ(ierr); 00089 00090 double mu_sliding = config.get("mu_sliding"); 00091 double minimum_temperature_for_sliding = config.get("minimum_temperature_for_sliding"); 00092 00093 MaskQuery m(*mask); 00094 00095 ierr = h_x.begin_access(); CHKERRQ(ierr); 00096 ierr = h_y.begin_access(); CHKERRQ(ierr); 00097 00098 ierr = mask->begin_access(); CHKERRQ(ierr); 00099 ierr = surface->begin_access(); CHKERRQ(ierr); 00100 ierr = bed->begin_access(); CHKERRQ(ierr); 00101 ierr = enthalpy->begin_access(); CHKERRQ(ierr); 00102 00103 ierr = velocity.begin_access(); CHKERRQ(ierr); 00104 ierr = basal_frictional_heating.begin_access(); CHKERRQ(ierr); 00105 for (PetscInt i=grid.xs; i<grid.xs+grid.xm; i++) { 00106 for (PetscInt j=grid.ys; j<grid.ys+grid.ym; j++) { 00107 if (m.ocean(i,j)) { 00108 velocity(i,j).u = 0.0; 00109 velocity(i,j).v = 0.0; 00110 basal_frictional_heating(i,j) = 0.0; 00111 } else { 00112 // basal velocity from SIA-type sliding law: not recommended! 00113 const PetscScalar 00114 myx = grid.x[i], 00115 myy = grid.y[j], 00116 myhx = 0.25 * ( h_x(i,j,0) + h_x(i-1,j,0) 00117 + h_x(i,j,1) + h_x(i,j-1,1)), 00118 myhy = 0.25 * ( h_y(i,j,0) + h_y(i-1,j,0) 00119 + h_y(i,j,1) + h_y(i,j-1,1)), 00120 alpha = sqrt(PetscSqr(myhx) + PetscSqr(myhy)); 00121 PetscScalar T, basalC; 00122 00123 // change r1200: new meaning of H 00124 const PetscScalar H = (*surface)(i,j) - (*bed)(i,j); 00125 00126 ierr = EC.getAbsTemp(enthalpy->getValZ(i,j,0.0), 00127 EC.getPressureFromDepth(H), T); CHKERRQ(ierr); 00128 00129 basalC = basalVelocitySIA(myx, myy, H, T, 00130 alpha, mu_sliding, 00131 minimum_temperature_for_sliding); 00132 velocity(i,j).u = - basalC * myhx; 00133 velocity(i,j).v = - basalC * myhy; 00134 // basal frictional heating; note P * dh/dx is x comp. of basal shear stress 00135 // in ice streams this result will be *overwritten* by 00136 // correctBasalFrictionalHeating() if useSSAVelocities==TRUE 00137 const PetscScalar P = ice.rho * standard_gravity * H; 00138 basal_frictional_heating(i,j) = - (P * myhx) * velocity(i,j).u - (P * myhy) * velocity(i,j).v; 00139 } 00140 } 00141 } 00142 00143 ierr = velocity.end_access(); CHKERRQ(ierr); 00144 ierr = basal_frictional_heating.end_access(); CHKERRQ(ierr); 00145 00146 ierr = h_y.end_access(); CHKERRQ(ierr); 00147 ierr = h_x.end_access(); CHKERRQ(ierr); 00148 00149 ierr = surface->end_access(); CHKERRQ(ierr); 00150 ierr = bed->end_access(); CHKERRQ(ierr); 00151 ierr = mask->end_access(); CHKERRQ(ierr); 00152 ierr = enthalpy->end_access(); CHKERRQ(ierr); 00153 00154 ierr = velocity.beginGhostComm(); CHKERRQ(ierr); 00155 ierr = velocity.endGhostComm(); CHKERRQ(ierr); 00156 00157 return 0; 00158 } 00159 00162 00182 PetscScalar SIA_Sliding::basalVelocitySIA(PetscScalar xIN, PetscScalar yIN, 00183 PetscScalar H, PetscScalar T, 00184 PetscScalar /*alpha*/, PetscScalar mu, 00185 PetscScalar min_T) const { 00186 00187 if (verification_mode) { 00188 // test 'E' mode 00189 const PetscScalar r1 = 200e3, r2 = 700e3, /* define region of sliding */ 00190 theta1 = 10 * (pi/180), theta2 = 40 * (pi/180); 00191 const PetscScalar x = fabs(xIN), y = fabs(yIN); 00192 const PetscScalar r = sqrt(x * x + y * y); 00193 PetscScalar theta; 00194 if (x < 1.0) 00195 theta = pi / 2.0; 00196 else 00197 theta = atan(y / x); 00198 00199 if ((r > r1) && (r < r2) && (theta > theta1) && (theta < theta2)) { 00200 // now INSIDE sliding region 00201 const PetscScalar rbot = (r2 - r1) * (r2 - r1), 00202 thetabot = (theta2 - theta1) * (theta2 - theta1); 00203 const PetscScalar mu_max = 2.5e-11; /* Pa^-1 m s^-1; max sliding coeff */ 00204 PetscScalar muE = mu_max * (4.0 * (r - r1) * (r2 - r) / rbot) 00205 * (4.0 * (theta - theta1) * (theta2 - theta) / thetabot); 00206 return muE * ice.rho * standard_gravity * H; 00207 } else 00208 return 0.0; 00209 } 00210 00211 if ((eisII_experiment == "G") || (eisII_experiment == "H")) { 00212 const PetscScalar Bfactor = 1e-3 / secpera; // m s^-1 Pa^-1 00213 PetscReal pressure = EC.getPressureFromDepth(H), E; 00214 EC.getEnthPermissive(T, 0.0, pressure, E); 00215 00216 if (eisII_experiment == "G") { 00217 return Bfactor * ice.rho * standard_gravity * H; 00218 } else if (eisII_experiment == "H") { 00219 if (EC.isTemperate(E, pressure)) { 00220 return Bfactor * ice.rho * standard_gravity * H; // ditto case G 00221 } else { 00222 return 0.0; 00223 } 00224 } 00225 return 0.0; // zero sliding for other tests 00226 } 00227 00228 // the "usual" case: 00229 if (T + ice.beta_CC_grad * H > min_T) { 00230 const PetscScalar p_over = ice.rho * standard_gravity * H; 00231 return mu * p_over; 00232 } else { 00233 return 0; 00234 } 00235 } 00236 00237 PetscErrorCode SIA_Sliding::compute_surface_gradient(IceModelVec2Stag &h_x, IceModelVec2Stag &h_y) { 00238 PetscErrorCode ierr; 00239 00240 const string method = config.get_string("surface_gradient_method"); 00241 00242 if ((method != "eta") && (method != "mahaffy") && (method != "haseloff")) { 00243 verbPrintf(1, grid.com, 00244 "PISM ERROR: value of surface_gradient_method, option -gradient, not valid ... ending\n"); 00245 PISMEnd(); 00246 } 00247 00248 if (method == "eta") { 00249 00250 ierr = surface_gradient_eta(h_x, h_y); CHKERRQ(ierr); 00251 00252 } else if (method == "haseloff") { 00253 00254 ierr = surface_gradient_haseloff(h_x, h_y); CHKERRQ(ierr); 00255 00256 } else if (method == "mahaffy") { 00257 00258 ierr = surface_gradient_mahaffy(h_x, h_y); CHKERRQ(ierr); 00259 00260 } else { 00261 SETERRQ(1, "can't happen"); 00262 } 00263 00264 return 0; 00265 } 00266 00268 PetscErrorCode SIA_Sliding::surface_gradient_eta(IceModelVec2Stag &h_x, IceModelVec2Stag &h_y) { 00269 PetscErrorCode ierr; 00270 00271 const PetscScalar n = ice.exponent(), // presumably 3.0 00272 etapow = (2.0 * n + 2.0)/n, // = 8/3 if n = 3 00273 invpow = 1.0 / etapow, 00274 dinvpow = (- n - 2.0) / (2.0 * n + 2.0); 00275 const PetscScalar dx = grid.dx, dy = grid.dy; // convenience 00276 IceModelVec2S eta = work_2d; 00277 00278 // compute eta = H^{8/3}, which is more regular, on reg grid 00279 ierr = thickness->begin_access(); CHKERRQ(ierr); 00280 ierr = eta.begin_access(); CHKERRQ(ierr); 00281 PetscInt GHOSTS = 2; 00282 for (PetscInt i = grid.xs - GHOSTS; i < grid.xs+grid.xm + GHOSTS; ++i) { 00283 for (PetscInt j = grid.ys - GHOSTS; j < grid.ys+grid.ym + GHOSTS; ++j) { 00284 eta(i,j) = pow((*thickness)(i,j), etapow); 00285 } 00286 } 00287 ierr = eta.end_access(); CHKERRQ(ierr); 00288 ierr = thickness->end_access(); CHKERRQ(ierr); 00289 00290 ierr = h_x.begin_access(); CHKERRQ(ierr); 00291 ierr = h_y.begin_access(); CHKERRQ(ierr); 00292 00293 // now use Mahaffy on eta to get grad h on staggered; 00294 // note grad h = (3/8) eta^{-5/8} grad eta + grad b because h = H + b 00295 ierr = bed->begin_access(); CHKERRQ(ierr); 00296 ierr = eta.begin_access(); CHKERRQ(ierr); 00297 for (PetscInt o=0; o<2; o++) { 00298 00299 GHOSTS = 1; 00300 for (PetscInt i = grid.xs - GHOSTS; i < grid.xs+grid.xm + GHOSTS; ++i) { 00301 for (PetscInt j = grid.ys - GHOSTS; j < grid.ys+grid.ym + GHOSTS; ++j) { 00302 if (o==0) { // If I-offset 00303 const PetscScalar mean_eta = 0.5 * (eta(i+1,j) + eta(i,j)); 00304 if (mean_eta > 0.0) { 00305 const PetscScalar factor = invpow * pow(mean_eta, dinvpow); 00306 h_x(i,j,o) = factor * (eta(i+1,j) - eta(i,j)) / dx; 00307 h_y(i,j,o) = factor * (+ eta(i+1,j+1) + eta(i,j+1) 00308 - eta(i+1,j-1) - eta(i,j-1)) / (4.0*dy); 00309 } else { 00310 h_x(i,j,o) = 0.0; 00311 h_y(i,j,o) = 0.0; 00312 } 00313 // now add bed slope to get actual h_x,h_y 00314 h_x(i,j,o) += bed->diff_x_stagE(i,j); 00315 h_y(i,j,o) += bed->diff_y_stagE(i,j); 00316 } else { // J-offset 00317 const PetscScalar mean_eta = 0.5 * (eta(i,j+1) + eta(i,j)); 00318 if (mean_eta > 0.0) { 00319 const PetscScalar factor = invpow * pow(mean_eta, dinvpow); 00320 h_y(i,j,o) = factor * (eta(i,j+1) - eta(i,j)) / dy; 00321 h_x(i,j,o) = factor * (+ eta(i+1,j+1) + eta(i+1,j) 00322 - eta(i-1,j+1) - eta(i-1,j)) / (4.0*dx); 00323 } else { 00324 h_y(i,j,o) = 0.0; 00325 h_x(i,j,o) = 0.0; 00326 } 00327 // now add bed slope to get actual h_x,h_y 00328 h_y(i,j,o) += bed->diff_y_stagN(i,j); 00329 h_x(i,j,o) += bed->diff_x_stagN(i,j); 00330 } 00331 } 00332 } 00333 } 00334 ierr = eta.end_access(); CHKERRQ(ierr); 00335 ierr = bed->end_access(); CHKERRQ(ierr); 00336 00337 ierr = h_y.end_access(); CHKERRQ(ierr); 00338 ierr = h_x.end_access(); CHKERRQ(ierr); 00339 00340 return 0; 00341 } 00342 00344 00348 PetscErrorCode SIA_Sliding::surface_gradient_haseloff(IceModelVec2Stag &h_x, IceModelVec2Stag &h_y) { 00349 PetscErrorCode ierr; 00350 00351 const PetscScalar Hicefree = 0.0; // standard for ice-free, in Haseloff 00352 const PetscScalar dx = grid.dx, dy = grid.dy; // convenience 00353 00354 ierr = h_x.begin_access(); CHKERRQ(ierr); 00355 ierr = h_y.begin_access(); CHKERRQ(ierr); 00356 00357 PetscScalar **h, **b, **H; 00358 ierr = bed->get_array(b); CHKERRQ(ierr); 00359 ierr = thickness->get_array(H); CHKERRQ(ierr); 00360 ierr = surface->get_array(h); CHKERRQ(ierr); 00361 for (PetscInt o=0; o<2; o++) { 00362 00363 PetscInt GHOSTS = 1; 00364 for (PetscInt i = grid.xs - GHOSTS; i < grid.xs+grid.xm + GHOSTS; ++i) { 00365 for (PetscInt j = grid.ys - GHOSTS; j < grid.ys+grid.ym + GHOSTS; ++j) { 00366 if (o==0) { // If I-offset 00367 const bool icefreeP = (H[i][j] <= Hicefree), 00368 icefreeE = (H[i+1][j] <= Hicefree), 00369 icefreeN = (H[i][j+1] <= Hicefree), 00370 icefreeS = (H[i][j-1] <= Hicefree), 00371 icefreeNE = (H[i+1][j+1] <= Hicefree), 00372 icefreeSE = (H[i+1][j-1] <= Hicefree); 00373 00374 PetscScalar hhE = h[i+1][j]; // east pseudo-surface elevation 00375 if (icefreeE && (b[i+1][j] > h[i][j] )) hhE = h[i][j]; 00376 if (icefreeP && (b[i][j] > h[i+1][j] )) hhE = h[i][j]; 00377 h_x(i,j,o) = (hhE - h[i][j]) / dx; 00378 00379 PetscScalar hhN = h[i][j+1]; // north pseudo-surface elevation 00380 if (icefreeN && (b[i][j+1] > h[i][j] )) hhN = h[i][j]; 00381 if (icefreeP && (b[i][j] > h[i][j+1] )) hhN = h[i][j]; 00382 PetscScalar hhS = h[i][j-1]; // south pseudo-surface elevation 00383 if (icefreeS && (b[i][j-1] > h[i][j] )) hhS = h[i][j]; 00384 if (icefreeP && (b[i][j] > h[i][j-1] )) hhS = h[i][j]; 00385 PetscScalar hhNE = h[i+1][j+1];// northeast pseudo-surface elevation 00386 if (icefreeNE && (b[i+1][j+1] > h[i+1][j] )) hhNE = h[i+1][j]; 00387 if (icefreeE && (b[i+1][j] > h[i+1][j+1])) hhNE = h[i+1][j]; 00388 PetscScalar hhSE = h[i+1][j-1];// southeast pseudo-surface elevation 00389 if (icefreeSE && (b[i+1][j-1] > h[i+1][j] )) hhSE = h[i+1][j]; 00390 if (icefreeE && (b[i+1][j] > h[i+1][j-1])) hhSE = h[i+1][j]; 00391 h_y(i,j,o) = (hhNE + hhN - hhSE - hhS) / (4.0 * dy); 00392 } else { // J-offset 00393 const bool icefreeP = (H[i][j] <= Hicefree), 00394 icefreeN = (H[i][j+1] <= Hicefree), 00395 icefreeE = (H[i+1][j] <= Hicefree), 00396 icefreeW = (H[i-1][j] <= Hicefree), 00397 icefreeNE = (H[i+1][j+1] <= Hicefree), 00398 icefreeNW = (H[i-1][j+1] <= Hicefree); 00399 00400 PetscScalar hhN = h[i][j+1]; // north pseudo-surface elevation 00401 if (icefreeN && (b[i][j+1] > h[i][j] )) hhN = h[i][j]; 00402 if (icefreeP && (b[i][j] > h[i][j+1] )) hhN = h[i][j]; 00403 h_y(i,j,o) = (hhN - h[i][j]) / dy; 00404 00405 PetscScalar hhE = h[i+1][j]; // east pseudo-surface elevation 00406 if (icefreeE && (b[i+1][j] > h[i][j] )) hhE = h[i][j]; 00407 if (icefreeP && (b[i][j] > h[i+1][j] )) hhE = h[i][j]; 00408 PetscScalar hhW = h[i-1][j]; // west pseudo-surface elevation 00409 if (icefreeW && (b[i-1][j] > h[i][j] )) hhW = h[i][j]; 00410 if (icefreeP && (b[i][j] > h[i-1][j] )) hhW = h[i][j]; 00411 PetscScalar hhNE = h[i+1][j+1];// northeast pseudo-surface elevation 00412 if (icefreeNE && (b[i+1][j+1] > h[i][j+1] )) hhNE = h[i][j+1]; 00413 if (icefreeN && (b[i][j+1] > h[i+1][j+1])) hhNE = h[i][j+1]; 00414 PetscScalar hhNW = h[i-1][j+1];// northwest pseudo-surface elevation 00415 if (icefreeNW && (b[i-1][j+1] > h[i][j+1] )) hhNW = h[i][j+1]; 00416 if (icefreeN && (b[i][j+1] > h[i-1][j+1])) hhNW = h[i][j+1]; 00417 h_x(i,j,o) = (hhNE + hhE - hhNW - hhW) / (4.0 * dx); 00418 } 00419 00420 } // j 00421 } // i 00422 } // o 00423 ierr = thickness->end_access(); CHKERRQ(ierr); 00424 ierr = bed->end_access(); CHKERRQ(ierr); 00425 ierr = surface->end_access(); CHKERRQ(ierr); 00426 00427 ierr = h_y.end_access(); CHKERRQ(ierr); 00428 ierr = h_x.end_access(); CHKERRQ(ierr); 00429 00430 return 0; 00431 } 00432 00435 PetscErrorCode SIA_Sliding::surface_gradient_mahaffy(IceModelVec2Stag &h_x, IceModelVec2Stag &h_y) { 00436 PetscErrorCode ierr; 00437 const PetscScalar dx = grid.dx, dy = grid.dy; // convenience 00438 00439 PetscScalar **h; 00440 ierr = h_x.begin_access(); CHKERRQ(ierr); 00441 ierr = h_y.begin_access(); CHKERRQ(ierr); 00442 ierr = surface->get_array(h); CHKERRQ(ierr); 00443 00444 for (PetscInt o=0; o<2; o++) { 00445 PetscInt GHOSTS = 1; 00446 for (PetscInt i = grid.xs - GHOSTS; i < grid.xs+grid.xm + GHOSTS; ++i) { 00447 for (PetscInt j = grid.ys - GHOSTS; j < grid.ys+grid.ym + GHOSTS; ++j) { 00448 if (o==0) { // If I-offset 00449 h_x(i,j,o) = (h[i+1][j] - h[i][j]) / dx; 00450 h_y(i,j,o) = (+ h[i+1][j+1] + h[i][j+1] 00451 - h[i+1][j-1] - h[i][j-1]) / (4.0*dy); 00452 } else { // J-offset 00453 h_y(i,j,o) = (h[i][j+1] - h[i][j]) / dy; 00454 h_x(i,j,o) = (+ h[i+1][j+1] + h[i+1][j] 00455 - h[i-1][j+1] - h[i-1][j]) / (4.0*dx); 00456 } 00457 } 00458 } 00459 } 00460 00461 ierr = surface->end_access(); CHKERRQ(ierr); 00462 ierr = h_y.end_access(); CHKERRQ(ierr); 00463 ierr = h_x.end_access(); CHKERRQ(ierr); 00464 00465 return 0; 00466 }
