ticeModelVec2.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
---
ticeModelVec2.cc (17411B)
---
1 // Copyright (C) 2008--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 <cstring>
20 #include <cstdlib>
21
22 #include <cassert>
23
24 #include <memory>
25 using std::dynamic_pointer_cast;
26
27 #include <petscdraw.h>
28 #include <petscdmda.h>
29
30 #include "pism/util/io/File.hh"
31 #include "iceModelVec.hh"
32 #include "IceGrid.hh"
33 #include "ConfigInterface.hh"
34
35 #include "error_handling.hh"
36 #include "iceModelVec_helpers.hh"
37
38 #include "pism_utilities.hh"
39 #include "io/io_helpers.hh"
40 #include "pism/util/Logger.hh"
41
42 namespace pism {
43
44 // this file contains methods for derived classes IceModelVec2S and IceModelVec2Int
45
46 // methods for base class IceModelVec are in "iceModelVec.cc"
47
48 IceModelVec2::IceModelVec2()
49 : IceModelVec() {
50 // empty
51 }
52
53 IceModelVec2::Ptr IceModelVec2::To2D(IceModelVec::Ptr input) {
54 IceModelVec2::Ptr result = dynamic_pointer_cast<IceModelVec2,IceModelVec>(input);
55 if (not (bool)result) {
56 throw RuntimeError(PISM_ERROR_LOCATION, "dynamic cast failure");
57 }
58 return result;
59 }
60
61 IceModelVec2V::~IceModelVec2V() {
62 // empty
63 }
64
65 IceModelVec2S::Ptr IceModelVec2S::To2DScalar(IceModelVec::Ptr input) {
66 IceModelVec2S::Ptr result = dynamic_pointer_cast<IceModelVec2S,IceModelVec>(input);
67 if (not (bool)result) {
68 throw RuntimeError(PISM_ERROR_LOCATION, "dynamic cast failure");
69 }
70 return result;
71 }
72
73 IceModelVec2S::IceModelVec2S() {
74 m_begin_end_access_use_dof = false;
75 }
76
77 IceModelVec2S::IceModelVec2S(IceGrid::ConstPtr grid, const std::string &name,
78 IceModelVecKind ghostedp, int width) {
79 m_begin_end_access_use_dof = false;
80
81 create(grid, name, ghostedp, width);
82 }
83
84
85 void IceModelVec2S::create(IceGrid::ConstPtr grid, const std::string &name, IceModelVecKind ghostedp, int width) {
86 assert(m_v == NULL);
87 IceModelVec2::create(grid, name, ghostedp, width, m_dof);
88 }
89
90
91 double** IceModelVec2S::get_array() {
92 begin_access();
93 return static_cast<double**>(m_array);
94 }
95
96 //! Sets an IceModelVec2 to the magnitude of a 2D vector field with components `v_x` and `v_y`.
97 /*! Computes the magnitude \b pointwise, so any of v_x, v_y and the IceModelVec
98 this is called on can be the same.
99
100 Does not communicate.
101 */
102 void IceModelVec2S::set_to_magnitude(const IceModelVec2S &v_x,
103 const IceModelVec2S &v_y) {
104 IceModelVec::AccessList list{this, &v_x, &v_y};
105
106 for (Points p(*m_grid); p; p.next()) {
107 const int i = p.i(), j = p.j();
108
109 (*this)(i,j) = sqrt(PetscSqr(v_x(i,j)) + PetscSqr(v_y(i,j)));
110 }
111
112 inc_state_counter(); // mark as modified
113
114 }
115
116 void IceModelVec2S::set_to_magnitude(const IceModelVec2V &input) {
117 IceModelVec::AccessList list{this, &input};
118
119 for (Points p(*m_grid); p; p.next()) {
120 const int i = p.i(), j = p.j();
121
122 (*this)(i, j) = input(i, j).magnitude();
123 }
124 }
125
126 //! Masks out all the areas where \f$ M \le 0 \f$ by setting them to `fill`.
127 void IceModelVec2S::mask_by(const IceModelVec2S &M, double fill) {
128 IceModelVec::AccessList list{this, &M};
129
130 for (Points p(*m_grid); p; p.next()) {
131 const int i = p.i(), j = p.j();
132
133 if (M(i,j) <= 0.0) {
134 (*this)(i,j) = fill;
135 }
136 }
137
138 inc_state_counter(); // mark as modified
139 }
140
141 void IceModelVec2::write_impl(const File &file) const {
142
143 assert(m_v != NULL);
144
145 // The simplest case:
146 if ((m_dof == 1) and (not m_has_ghosts)) {
147 IceModelVec::write_impl(file);
148 return;
149 }
150
151 // Get the dof=1, stencil_width=0 DMDA (components are always scalar
152 // and we just need a global Vec):
153 petsc::DM::Ptr da2 = m_grid->get_dm(1, 0);
154
155 // a temporary one-component vector, distributed across processors
156 // the same way v is
157 petsc::TemporaryGlobalVec tmp(da2);
158
159 Logger::ConstPtr log = m_grid->ctx()->log();
160 log->message(4, " Writing %s...\n", m_name.c_str());
161
162 for (unsigned int j = 0; j < m_dof; ++j) {
163 IceModelVec2::get_dof(da2, tmp, j);
164
165 petsc::VecArray tmp_array(tmp);
166 io::write_spatial_variable(m_metadata[j], *m_grid, file,
167 tmp_array.get());
168 }
169 }
170
171 void IceModelVec2::read_impl(const File &nc, const unsigned int time) {
172
173 if ((m_dof == 1) and (not m_has_ghosts)) {
174 IceModelVec::read_impl(nc, time);
175 return;
176 }
177
178 Logger::ConstPtr log = m_grid->ctx()->log();
179 log->message(4, " Reading %s...\n", m_name.c_str());
180
181 assert(m_v != NULL);
182
183 // Get the dof=1, stencil_width=0 DMDA (components are always scalar
184 // and we just need a global Vec):
185 petsc::DM::Ptr da2 = m_grid->get_dm(1, 0);
186
187 // a temporary one-component vector, distributed across processors
188 // the same way v is
189 petsc::TemporaryGlobalVec tmp(da2);
190
191 for (unsigned int j = 0; j < m_dof; ++j) {
192
193 {
194 petsc::VecArray tmp_array(tmp);
195 io::read_spatial_variable(m_metadata[j], *m_grid, nc, time, tmp_array.get());
196 }
197
198 IceModelVec2::set_dof(da2, tmp, j);
199 }
200
201 // The calls above only set the values owned by a processor, so we need to
202 // communicate if m_has_ghosts == true:
203 if (m_has_ghosts) {
204 update_ghosts();
205 }
206 }
207
208 void IceModelVec2::regrid_impl(const File &file, RegriddingFlag flag,
209 double default_value) {
210 if ((m_dof == 1) and (not m_has_ghosts)) {
211 IceModelVec::regrid_impl(file, flag, default_value);
212 return;
213 }
214
215 // Get the dof=1, stencil_width=0 DMDA (components are always scalar
216 // and we just need a global Vec):
217 petsc::DM::Ptr da2 = m_grid->get_dm(1, 0);
218
219 // a temporary one-component vector, distributed across processors
220 // the same way v is
221 petsc::TemporaryGlobalVec tmp(da2);
222
223 const bool allow_extrapolation = m_grid->ctx()->config()->get_flag("grid.allow_extrapolation");
224
225 for (unsigned int j = 0; j < m_dof; ++j) {
226 {
227 petsc::VecArray tmp_array(tmp);
228 io::regrid_spatial_variable(m_metadata[j], *m_grid, file, flag,
229 m_report_range, allow_extrapolation,
230 default_value, m_interpolation_type, tmp_array.get());
231 }
232
233 IceModelVec2::set_dof(da2, tmp, j);
234 }
235
236 // The calls above only set the values owned by a processor, so we need to
237 // communicate if m_has_ghosts == true:
238 if (m_has_ghosts == true) {
239 update_ghosts();
240 }
241 }
242
243 //! \brief View a 2D field.
244 void IceModelVec2::view(int viewer_size) const {
245 petsc::Viewer::Ptr viewers[2];
246
247 if (m_dof > 2) {
248 throw RuntimeError(PISM_ERROR_LOCATION, "dof > 2 is not supported");
249 }
250
251 for (unsigned int j = 0; j < std::min(m_dof, 2U); ++j) {
252 std::string
253 c_name = m_metadata[j].get_name(),
254 long_name = m_metadata[j].get_string("long_name"),
255 glaciological_units = m_metadata[j].get_string("glaciological_units"),
256 title = long_name + " (" + glaciological_units + ")";
257
258 if (not m_map_viewers[c_name]) {
259 m_map_viewers[c_name].reset(new petsc::Viewer(m_grid->com, title, viewer_size,
260 m_grid->Lx(), m_grid->Ly()));
261 }
262
263 viewers[j] = m_map_viewers[c_name];
264 }
265
266 view(viewers[0], viewers[1]);
267 }
268
269 //! \brief View a 2D vector field using existing PETSc viewers.
270 void IceModelVec2::view(petsc::Viewer::Ptr v1, petsc::Viewer::Ptr v2) const {
271 PetscErrorCode ierr;
272
273 petsc::Viewer::Ptr viewers[2] = {v1, v2};
274
275 // Get the dof=1, stencil_width=0 DMDA (components are always scalar
276 // and we just need a global Vec):
277 petsc::DM::Ptr da2 = m_grid->get_dm(1, 0);
278
279 petsc::TemporaryGlobalVec tmp(da2);
280
281 for (unsigned int i = 0; i < std::min(m_dof, 2U); ++i) {
282 std::string
283 long_name = m_metadata[i].get_string("long_name"),
284 units = m_metadata[i].get_string("units"),
285 glaciological_units = m_metadata[i].get_string("glaciological_units"),
286 title = long_name + " (" + glaciological_units + ")";
287
288 if (not (bool)viewers[i]) {
289 continue;
290 }
291
292 PetscViewer v = *viewers[i].get();
293
294 PetscDraw draw = NULL;
295 ierr = PetscViewerDrawGetDraw(v, 0, &draw);
296 PISM_CHK(ierr, "PetscViewerDrawGetDraw");
297
298 ierr = PetscDrawSetTitle(draw, title.c_str());
299 PISM_CHK(ierr, "PetscDrawSetTitle");
300
301 IceModelVec2::get_dof(da2, tmp, i);
302
303 convert_vec(tmp, m_metadata[i].unit_system(),
304 units, glaciological_units);
305
306 double bounds[2] = {0.0, 0.0};
307 ierr = VecMin(tmp, NULL, &bounds[0]); PISM_CHK(ierr, "VecMin");
308 ierr = VecMax(tmp, NULL, &bounds[1]); PISM_CHK(ierr, "VecMax");
309
310 if (bounds[0] == bounds[1]) {
311 bounds[0] = -1.0;
312 bounds[1] = 1.0;
313 }
314
315 ierr = PetscViewerDrawSetBounds(v, 1, bounds);
316 PISM_CHK(ierr, "PetscViewerDrawSetBounds");
317
318 ierr = VecView(tmp, v);
319 PISM_CHK(ierr, "VecView");
320 }
321 }
322
323 //! \brief Returns the x-derivative at i,j approximated using centered finite
324 //! differences.
325 double IceModelVec2S::diff_x(int i, int j) const {
326 return ((*this)(i + 1,j) - (*this)(i - 1,j)) / (2 * m_grid->dx());
327 }
328
329 //! \brief Returns the y-derivative at i,j approximated using centered finite
330 //! differences.
331 double IceModelVec2S::diff_y(int i, int j) const {
332 return ((*this)(i,j + 1) - (*this)(i,j - 1)) / (2 * m_grid->dy());
333 }
334
335 //! \brief Returns the x-derivative at i,j approximated using centered finite
336 //! differences. Respects grid periodicity and uses one-sided FD at grid edges
337 //! if necessary.
338 double IceModelVec2S::diff_x_p(int i, int j) const {
339 if (m_grid->periodicity() & X_PERIODIC) {
340 return diff_x(i,j);
341 }
342
343 if (i == 0) {
344 return ((*this)(i + 1,j) - (*this)(i,j)) / (m_grid->dx());
345 } else if (i == (int)m_grid->Mx() - 1) {
346 return ((*this)(i,j) - (*this)(i - 1,j)) / (m_grid->dx());
347 } else {
348 return diff_x(i,j);
349 }
350 }
351
352 //! \brief Returns the y-derivative at i,j approximated using centered finite
353 //! differences. Respects grid periodicity and uses one-sided FD at grid edges
354 //! if necessary.
355 double IceModelVec2S::diff_y_p(int i, int j) const {
356 if (m_grid->periodicity() & Y_PERIODIC) {
357 return diff_y(i,j);
358 }
359
360 if (j == 0) {
361 return ((*this)(i,j + 1) - (*this)(i,j)) / (m_grid->dy());
362 } else if (j == (int)m_grid->My() - 1) {
363 return ((*this)(i,j) - (*this)(i,j - 1)) / (m_grid->dy());
364 } else {
365 return diff_y(i,j);
366 }
367 }
368
369 //! Sums up all the values in an IceModelVec2S object. Ignores ghosts.
370 /*! Avoids copying to a "global" vector.
371 */
372 double IceModelVec2S::sum() const {
373 double result = 0;
374
375 IceModelVec::AccessList list(*this);
376
377 // sum up the local part:
378 for (Points p(*m_grid); p; p.next()) {
379 result += (*this)(p.i(), p.j());
380 }
381
382 // find the global sum:
383 return GlobalSum(m_grid->com, result);
384 }
385
386
387 //! Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
388 double IceModelVec2S::max() const {
389 IceModelVec::AccessList list(*this);
390
391 double result = (*this)(m_grid->xs(),m_grid->ys());
392 for (Points p(*m_grid); p; p.next()) {
393 result = std::max(result,(*this)(p.i(), p.j()));
394 }
395
396 return GlobalMax(m_grid->com, result);
397 }
398
399
400 //! Finds maximum over all the absolute values in an IceModelVec2S object. Ignores ghosts.
401 double IceModelVec2S::absmax() const {
402
403 double result = 0.0;
404
405 IceModelVec::AccessList list(*this);
406 for (Points p(*m_grid); p; p.next()) {
407 result = std::max(result, fabs((*this)(p.i(), p.j())));
408 }
409
410 return GlobalMax(m_grid->com, result);
411 }
412
413
414 //! Finds minimum over all the values in an IceModelVec2S object. Ignores ghosts.
415 double IceModelVec2S::min() const {
416 IceModelVec::AccessList list(*this);
417
418 double result = (*this)(m_grid->xs(), m_grid->ys());
419 for (Points p(*m_grid); p; p.next()) {
420 result = std::min(result,(*this)(p.i(), p.j()));
421 }
422
423 return GlobalMin(m_grid->com, result);
424 }
425
426
427 // IceModelVec2
428 IceModelVec2::IceModelVec2(IceGrid::ConstPtr grid, const std::string &short_name,
429 IceModelVecKind ghostedp, unsigned int stencil_width, int dof) {
430 create(grid, short_name, ghostedp, stencil_width, dof);
431 }
432
433 void IceModelVec2::get_component(unsigned int n, IceModelVec2S &result) const {
434
435 IceModelVec2::get_dof(result.dm(), result.m_v, n);
436 }
437
438 void IceModelVec2::set_component(unsigned int n, const IceModelVec2S &source) {
439
440 IceModelVec2::set_dof(source.dm(), source.m_v, n);
441 }
442
443 void IceModelVec2::create(IceGrid::ConstPtr grid, const std::string & name,
444 IceModelVecKind ghostedp,
445 unsigned int stencil_width, int dof) {
446 PetscErrorCode ierr;
447 assert(m_v == NULL);
448
449 m_dof = dof;
450 m_grid = grid;
451
452 if ((m_dof != 1) || (stencil_width > m_grid->ctx()->config()->get_number("grid.max_stencil_width"))) {
453 m_da_stencil_width = stencil_width;
454 } else {
455 m_da_stencil_width = m_grid->ctx()->config()->get_number("grid.max_stencil_width");
456 }
457
458 // initialize the da member:
459 m_da = m_grid->get_dm(this->m_dof, this->m_da_stencil_width);
460
461 if (ghostedp) {
462 ierr = DMCreateLocalVector(*m_da, m_v.rawptr());
463 PISM_CHK(ierr, "DMCreateLocalVector");
464 } else {
465 ierr = DMCreateGlobalVector(*m_da, m_v.rawptr());
466 PISM_CHK(ierr, "DMCreateGlobalVector");
467 }
468
469 m_has_ghosts = (ghostedp == WITH_GHOSTS);
470 m_name = name;
471
472 if (m_dof == 1) {
473 m_metadata.push_back(SpatialVariableMetadata(m_grid->ctx()->unit_system(),
474 name));
475 } else {
476
477 for (unsigned int j = 0; j < m_dof; ++j) {
478 m_metadata.push_back(SpatialVariableMetadata(m_grid->ctx()->unit_system(),
479 pism::printf("%s[%d]", m_name.c_str(), j)));
480 }
481 }
482 }
483
484 void IceModelVec2S::add(double alpha, const IceModelVec &x) {
485 add_2d<IceModelVec2S>(this, alpha, &x, this);
486 }
487
488 void IceModelVec2S::add(double alpha, const IceModelVec &x, IceModelVec &result) const {
489 add_2d<IceModelVec2S>(this, alpha, &x, &result);
490 }
491
492 void IceModelVec2S::copy_from(const IceModelVec &source) {
493 copy_2d<IceModelVec2S>(&source, this);
494 }
495
496 // IceModelVec2Stag
497
498 IceModelVec2Stag::IceModelVec2Stag()
499 : IceModelVec2() {
500 m_dof = 2;
501 m_begin_end_access_use_dof = true;
502 }
503
504 IceModelVec2Stag::Ptr IceModelVec2Stag::ToStaggered(IceModelVec::Ptr input) {
505 IceModelVec2Stag::Ptr result = dynamic_pointer_cast<IceModelVec2Stag,IceModelVec>(input);
506 if (not (bool)result) {
507 throw RuntimeError(PISM_ERROR_LOCATION, "dynamic cast failure");
508 }
509 return result;
510 }
511
512 IceModelVec2Stag::IceModelVec2Stag(IceGrid::ConstPtr grid, const std::string &name,
513 IceModelVecKind ghostedp,
514 unsigned int stencil_width)
515 : IceModelVec2() {
516 m_dof = 2;
517 m_begin_end_access_use_dof = true;
518
519 create(grid, name, ghostedp, stencil_width);
520 }
521
522 void IceModelVec2Stag::create(IceGrid::ConstPtr grid, const std::string &name,
523 IceModelVecKind ghostedp,
524 unsigned int stencil_width) {
525
526 IceModelVec2::create(grid, name, ghostedp, stencil_width, m_dof);
527 }
528
529 //! Averages staggered grid values of a scalar field and puts them on a regular grid.
530 /*!
531 * The current IceModelVec needs to have ghosts.
532 */
533 void IceModelVec2Stag::staggered_to_regular(IceModelVec2S &result) const {
534 IceModelVec::AccessList list{this, &result};
535
536 for (Points p(*m_grid); p; p.next()) {
537 const int i = p.i(), j = p.j();
538
539 result(i,j) = 0.25 * ((*this)(i,j,0) + (*this)(i,j,1)
540 + (*this)(i,j-1,1) + (*this)(i-1,j,0));
541 }
542 }
543
544 //! \brief Averages staggered grid values of a 2D vector field (u on the
545 //! i-offset, v on the j-offset) and puts them on a regular grid.
546 /*!
547 * The current IceModelVec needs to have ghosts.
548 */
549 void IceModelVec2Stag::staggered_to_regular(IceModelVec2V &result) const {
550 IceModelVec::AccessList list{this, &result};
551
552 for (Points p(*m_grid); p; p.next()) {
553 const int i = p.i(), j = p.j();
554
555 result(i,j).u = 0.5 * ((*this)(i-1,j,0) + (*this)(i,j,0));
556 result(i,j).v = 0.5 * ((*this)(i,j-1,1) + (*this)(i,j,1));
557 }
558 }
559
560
561 //! For each component, finds the maximum over all the absolute values. Ignores ghosts.
562 std::vector<double> IceModelVec2Stag::absmaxcomponents() const {
563 std::vector<double> z(2, 0.0);
564
565 IceModelVec::AccessList list(*this);
566 for (Points p(*m_grid); p; p.next()) {
567 const int i = p.i(), j = p.j();
568
569 z[0] = std::max(z[0], fabs((*this)(i, j, 0)));
570 z[1] = std::max(z[1], fabs((*this)(i, j, 1)));
571 }
572
573 z[0] = GlobalMax(m_grid->com, z[0]);
574 z[1] = GlobalMax(m_grid->com, z[1]);
575
576 return z;
577 }
578
579 IceModelVec2Int::IceModelVec2Int()
580 : IceModelVec2S() {
581 m_interpolation_type = NEAREST;
582 }
583
584
585 IceModelVec2Int::IceModelVec2Int(IceGrid::ConstPtr grid, const std::string &name,
586 IceModelVecKind ghostedp, int width)
587 : IceModelVec2S(grid, name, ghostedp, width) {
588 m_interpolation_type = NEAREST;
589 }
590
591
592 } // end of namespace pism