tBedSmoother.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tBedSmoother.cc (14502B)
---
1 // Copyright (C) 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 Ed Bueler and Constantine Khroulev
2 //
3 // This file is part of PISM.
4 //
5 // PISM is free software; you can redistribute it and/or modify it under the
6 // terms of the GNU General Public License as published by the Free Software
7 // Foundation; either version 3 of the License, or (at your option) any later
8 // version.
9 //
10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with PISM; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19 #include <cassert>
20
21 #include "BedSmoother.hh"
22 #include "pism/util/Mask.hh"
23 #include "pism/util/IceGrid.hh"
24 #include "pism/util/petscwrappers/Vec.hh"
25 #include "pism/util/IceModelVec2CellType.hh"
26
27 #include "pism/util/error_handling.hh"
28 #include "pism/util/pism_utilities.hh"
29 #include "pism/util/Logger.hh"
30
31 namespace pism {
32 namespace stressbalance {
33
34 BedSmoother::BedSmoother(IceGrid::ConstPtr g, int MAX_GHOSTS)
35 : m_grid(g), m_config(g->ctx()->config()) {
36
37 const Logger &log = *m_grid->ctx()->log();
38
39 {
40 // allocate Vecs that live on all procs; all have to be as "wide" as any of
41 // their prospective uses
42 m_topgsmooth.create(m_grid, "topgsmooth", WITH_GHOSTS, MAX_GHOSTS);
43 m_topgsmooth.set_attrs("bed_smoother_tool",
44 "smoothed bed elevation, in bed roughness parameterization",
45 "m", "m", "", 0);
46 m_maxtl.create(m_grid, "maxtl", WITH_GHOSTS, MAX_GHOSTS);
47 m_maxtl.set_attrs("bed_smoother_tool",
48 "maximum elevation in local topography patch, in bed roughness parameterization",
49 "m", "m", "", 0);
50 m_C2.create(m_grid, "C2bedsmooth", WITH_GHOSTS, MAX_GHOSTS);
51 m_C2.set_attrs("bed_smoother_tool",
52 "polynomial coeff of H^-2, in bed roughness parameterization",
53 "m2", "m2", "", 0);
54 m_C3.create(m_grid, "C3bedsmooth", WITH_GHOSTS, MAX_GHOSTS);
55 m_C3.set_attrs("bed_smoother_tool",
56 "polynomial coeff of H^-3, in bed roughness parameterization",
57 "m3", "m3", "", 0);
58 m_C4.create(m_grid, "C4bedsmooth", WITH_GHOSTS, MAX_GHOSTS);
59 m_C4.set_attrs("bed_smoother_tool",
60 "polynomial coeff of H^-4, in bed roughness parameterization",
61 "m4", "m4", "", 0);
62
63 // allocate Vecs that live on processor 0:
64 m_topgp0 = m_topgsmooth.allocate_proc0_copy();
65 m_topgsmoothp0 = m_topgsmooth.allocate_proc0_copy();
66 m_maxtlp0 = m_maxtl.allocate_proc0_copy();
67 m_C2p0 = m_C2.allocate_proc0_copy();
68 m_C3p0 = m_C3.allocate_proc0_copy();
69 m_C4p0 = m_C4.allocate_proc0_copy();
70 }
71
72 m_Glen_exponent = m_config->get_number("stress_balance.sia.Glen_exponent"); // choice is SIA; see #285
73 m_smoothing_range = m_config->get_number("stress_balance.sia.bed_smoother.range");
74
75 if (m_smoothing_range > 0.0) {
76 log.message(2,
77 "* Initializing bed smoother object with %.3f km half-width ...\n",
78 units::convert(m_grid->ctx()->unit_system(), m_smoothing_range, "m", "km"));
79 }
80
81 // Make sure that Nx and Ny are initialized. In most cases SIAFD::update() will call
82 // preprocess_bed() and set appropriate values, but in a zero-length (-y 0) run IceModel does not
83 // call SIAFD::update()... We may need to re-structure this class so that everything is
84 // initialized right after construction and users don't have to call preprocess_bed() manually.
85 m_Nx = -1;
86 m_Ny = -1;
87 }
88
89
90 BedSmoother::~BedSmoother() {
91 // empty
92 }
93
94 /*!
95 Input lambda gives physical half-width (in m) of square over which to do the
96 average. Only square smoothing domains are allowed with this call, which is the
97 default case.
98 */
99 void BedSmoother::preprocess_bed(const IceModelVec2S &topg) {
100
101 if (m_smoothing_range <= 0.0) {
102 // smoothing completely inactive. we transfer the original bed topg,
103 // including ghosts, to public member topgsmooth ...
104 topg.update_ghosts(m_topgsmooth);
105 // and we tell theta() to return theta=1
106 m_Nx = -1;
107 m_Ny = -1;
108 return;
109 }
110
111 // determine Nx, Ny, which are always at least one if m_smoothing_range > 0
112 m_Nx = static_cast<int>(ceil(m_smoothing_range / m_grid->dx()));
113 m_Ny = static_cast<int>(ceil(m_smoothing_range / m_grid->dy()));
114 if (m_Nx < 1) {
115 m_Nx = 1;
116 }
117 if (m_Ny < 1) {
118 m_Ny = 1;
119 }
120
121 preprocess_bed(topg, m_Nx, m_Ny);
122 }
123
124 const IceModelVec2S& BedSmoother::smoothed_bed() const {
125 return m_topgsmooth;
126 }
127
128 /*!
129 Inputs Nx,Ny gives half-width in number of grid points, over which to do the
130 average.
131 */
132 void BedSmoother::preprocess_bed(const IceModelVec2S &topg,
133 unsigned int Nx, unsigned int Ny) {
134
135 if ((Nx >= m_grid->Mx()) || (Ny >= m_grid->My())) {
136 throw RuntimeError(PISM_ERROR_LOCATION, "input Nx, Ny in bed smoother is too large because\n"
137 "domain of smoothing exceeds IceGrid domain");
138 }
139 m_Nx = Nx;
140 m_Ny = Ny;
141
142 topg.put_on_proc0(*m_topgp0);
143 smooth_the_bed_on_proc0();
144 // next call *does indeed* fill ghosts in topgsmooth
145 m_topgsmooth.get_from_proc0(*m_topgsmoothp0);
146
147 compute_coefficients_on_proc0();
148 // following calls *do* fill the ghosts
149 m_maxtl.get_from_proc0(*m_maxtlp0);
150 m_C2.get_from_proc0(*m_C2p0);
151 m_C3.get_from_proc0(*m_C3p0);
152 m_C4.get_from_proc0(*m_C4p0);
153 }
154
155
156 //! Computes the smoothed bed by a simple average over a rectangle of grid points.
157 void BedSmoother::smooth_the_bed_on_proc0() {
158
159 ParallelSection rank0(m_grid->com);
160 try {
161 if (m_grid->rank() == 0) {
162 const int Mx = (int)m_grid->Mx();
163 const int My = (int)m_grid->My();
164
165 petsc::VecArray2D
166 b0(*m_topgp0, Mx, My),
167 bs(*m_topgsmoothp0, Mx, My);
168
169 for (int j=0; j < My; j++) {
170 for (int i=0; i < Mx; i++) {
171 // average only over those points which are in the grid; do
172 // not wrap periodically
173 double sum = 0.0, count = 0.0;
174 for (int r = -m_Nx; r <= m_Nx; r++) {
175 for (int s = -m_Ny; s <= m_Ny; s++) {
176 if ((i+r >= 0) and (i+r < Mx) and (j+s >= 0) and (j+s < My)) {
177 sum += b0(i+r, j+s);
178 count += 1.0;
179 }
180 }
181 }
182 // unprotected division by count but r=0,s=0 case guarantees count>=1
183 bs(i, j) = sum / count;
184 }
185 }
186 }
187 } catch (...) {
188 rank0.failed();
189 }
190 rank0.check();
191 }
192
193
194 void BedSmoother::compute_coefficients_on_proc0() {
195
196 const unsigned int Mx = m_grid->Mx(), My = m_grid->My();
197
198 ParallelSection rank0(m_grid->com);
199 try {
200 if (m_grid->rank() == 0) {
201 petsc::VecArray2D
202 b0(*m_topgp0, Mx, My),
203 bs(*m_topgsmoothp0, Mx, My),
204 mt(*m_maxtlp0, Mx, My),
205 c2(*m_C2p0, Mx, My),
206 c3(*m_C3p0, Mx, My),
207 c4(*m_C4p0, Mx, My);
208
209 for (int j=0; j < (int)My; j++) {
210 for (int i=0; i < (int)Mx; i++) {
211 // average only over those points which are in the grid
212 // do not wrap periodically
213 double
214 topgs = bs(i, j),
215 maxtltemp = 0.0,
216 sum2 = 0.0,
217 sum3 = 0.0,
218 sum4 = 0.0,
219 count = 0.0;
220
221 for (int r = -m_Nx; r <= m_Nx; r++) {
222 for (int s = -m_Ny; s <= m_Ny; s++) {
223 if ((i+r >= 0) && (i+r < (int)Mx) && (j+s >= 0) && (j+s < (int)My)) {
224 // tl is elevation of local topography at a pt in patch
225 const double tl = b0(i+r, j+s) - topgs;
226 maxtltemp = std::max(maxtltemp, tl);
227 // accumulate 2nd, 3rd, and 4th powers with only 3 multiplications
228 const double tl2 = tl * tl;
229 sum2 += tl2;
230 sum3 += tl2 * tl;
231 sum4 += tl2 * tl2;
232 count += 1.0;
233 }
234 }
235 }
236 mt(i, j) = maxtltemp;
237
238 // unprotected division by count but r=0,s=0 case guarantees count>=1
239 c2(i, j) = sum2 / count;
240 c3(i, j) = sum3 / count;
241 c4(i, j) = sum4 / count;
242 }
243 }
244
245 // scale the coeffs in Taylor series
246 const double
247 n = m_Glen_exponent,
248 k = (n + 2) / n,
249 s2 = k * (2 * n + 2) / (2 * n),
250 s3 = s2 * (3 * n + 2) / (3 * n),
251 s4 = s3 * (4 * n + 2) / (4 * n);
252
253 PetscErrorCode ierr;
254 ierr = VecScale(*m_C2p0,s2);
255 PISM_CHK(ierr, "VecScale");
256
257 ierr = VecScale(*m_C3p0,s3);
258 PISM_CHK(ierr, "VecScale");
259
260 ierr = VecScale(*m_C4p0,s4);
261 PISM_CHK(ierr, "VecScale");
262 }
263 } catch (...) {
264 rank0.failed();
265 }
266 rank0.check();
267 }
268
269
270 //! Computes a smoothed thickness map.
271 /*!
272 The result `thksmooth` is the difference between the given upper surface
273 elevation (`usurf`) and the stored smoothed bed topography (`topgsmooth`),
274 except where the given original thickness (`thk`) is zero. In places where
275 the original thickness is zero, the result `thksmooth` is also set to zero.
276
277 Ghosted values are updated directly and no communication occurs. In fact,
278 we \e assume `usurf`, `thk`, and `thksmooth` all have stencil width at least
279 equal to GHOSTS. We \e check whether `topgsmooth`, which has stencil width
280 maxGHOSTS, has at least GHOSTS stencil width, and throw an error if not.
281
282 Call preprocess_bed() first.
283 */
284 void BedSmoother::smoothed_thk(const IceModelVec2S &usurf,
285 const IceModelVec2S &thk,
286 const IceModelVec2CellType &mask,
287 IceModelVec2S &result) const {
288
289 IceModelVec::AccessList list{&mask, &m_maxtl, &result, &thk, &m_topgsmooth, &usurf};
290
291 unsigned int GHOSTS = result.stencil_width();
292 assert(mask.stencil_width() >= GHOSTS);
293 assert(m_maxtl.stencil_width() >= GHOSTS);
294 assert(thk.stencil_width() >= GHOSTS);
295 assert(m_topgsmooth.stencil_width() >= GHOSTS);
296 assert(usurf.stencil_width() >= GHOSTS);
297
298 ParallelSection loop(m_grid->com);
299 try {
300 for (PointsWithGhosts p(*m_grid, GHOSTS); p; p.next()) {
301 const int i = p.i(), j = p.j();
302
303 if (thk(i, j) < 0.0) {
304 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "BedSmoother detects negative original thickness\n"
305 "at location (i, j) = (%d, %d) ... ending", i, j);
306 } else if (thk(i, j) == 0.0) {
307 result(i, j) = 0.0;
308 } else if (m_maxtl(i, j) >= thk(i, j)) {
309 result(i, j) = thk(i, j);
310 } else {
311 if (mask.grounded(i, j)) {
312 // if grounded, compute smoothed thickness as the difference of ice
313 // surface elevation and smoothed bed elevation
314 const double thks_try = usurf(i, j) - m_topgsmooth(i, j);
315 result(i, j) = (thks_try > 0.0) ? thks_try : 0.0;
316 } else {
317 // if floating, use original thickness (note: surface elevation was
318 // computed using this thickness and the sea level elevation)
319 result(i, j) = thk(i, j);
320 }
321 }
322 }
323 } catch (...) {
324 loop.failed();
325 }
326 loop.check();
327 }
328
329
330 /*!
331 Implements the strategy for computing \f$\theta(h,x,y)\f$ from previously-
332 stored coefficients, described on [Bed roughness parameterization](@ref bedrough) page and in [\ref
333 Schoofbasaltopg2003].
334
335 Specifically, \f$\theta = \omega^{-n}\f$ where \f$\omega\f$ is a local average
336 of a rational function of surface elevation, approximated here by a Taylor polynomial:
337 \f[ \omega = \fint \left(1 - \frac{\tilde b(x_1,x_2,\xi_1,\xi_2)}{H}
338 \right)^{-(n+2)/n}\,d\xi_1\,d\xi_2
339 \approx 1 + C_2 H^{-2} + C_3 H^{-3} + C_4 H^{-4} \f]
340 where \f$h =\f$ usurf, \f$H = h -\f$ topgsmooth and \f$\tilde b\f$ is the local
341 bed topography, a function with mean zero. The coefficients \f$C_2,C_3,C_4\f$,
342 which depend on \f$x,y\f$, are precomputed by `preprocess_bed()`.
343
344 Ghosted values are updated directly and no communication occurs. In fact,
345 we \e assume `usurf` and `theta` have stencil width at least
346 equal to GHOSTS. We \e check whether `topgsmooth`, which has stencil width
347 maxGHOSTS, has at least GHOSTS stencil width, and throw an error if not.
348
349 Call preprocess_bed() first.
350 */
351 void BedSmoother::theta(const IceModelVec2S &usurf, IceModelVec2S &result) const {
352
353 if ((m_Nx < 0) || (m_Ny < 0)) {
354 result.set(1.0);
355 return;
356 }
357
358 IceModelVec::AccessList list{&m_C2, &m_C3, &m_C4, &m_maxtl, &result, &m_topgsmooth, &usurf};
359
360 unsigned int GHOSTS = result.stencil_width();
361 assert(m_C2.stencil_width() >= GHOSTS);
362 assert(m_C3.stencil_width() >= GHOSTS);
363 assert(m_C4.stencil_width() >= GHOSTS);
364 assert(m_maxtl.stencil_width() >= GHOSTS);
365 assert(m_topgsmooth.stencil_width() >= GHOSTS);
366 assert(usurf.stencil_width() >= GHOSTS);
367
368 const double
369 theta_min = m_config->get_number("stress_balance.sia.bed_smoother.theta_min"),
370 theta_max = 1.0;
371
372 ParallelSection loop(m_grid->com);
373 try {
374 for (PointsWithGhosts p(*m_grid, GHOSTS); p; p.next()) {
375 const int i = p.i(), j = p.j();
376
377 const double H = usurf(i, j) - m_topgsmooth(i, j);
378 if (H > m_maxtl(i, j)) {
379 // thickness exceeds maximum variation in patch of local topography,
380 // so ice buries local topography; note maxtl >= 0 always
381 const double Hinv = 1.0 / std::max(H, 1.0);
382 double omega = 1.0 + Hinv*Hinv * (m_C2(i, j) + Hinv * (m_C3(i, j) + Hinv*m_C4(i, j)));
383 if (omega <= 0) { // this check *should not* be necessary: p4(s) > 0
384 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
385 "omega is negative for i=%d, j=%d\n"
386 "in BedSmoother.theta()", i, j);
387 }
388
389 if (omega < 0.001) { // this check *should not* be necessary
390 omega = 0.001;
391 }
392
393 result(i, j) = pow(omega, -m_Glen_exponent);
394 } else {
395 result(i, j) = 0.00;
396 }
397
398 result(i, j) = clip(result(i, j), theta_min, theta_max);
399 }
400 } catch (...) {
401 loop.failed();
402 }
403 loop.check();
404 }
405
406
407 } // end of namespace stressbalance
408 } // end of namespace pism