ticeModelVec3.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
---
ticeModelVec3.cc (9333B)
---
1 // Copyright (C) 2008--2018 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 #include <cassert>
22
23 #include <memory>
24 using std::dynamic_pointer_cast;
25
26 #include <petscdmda.h>
27
28 #include "iceModelVec.hh"
29 #include "IceGrid.hh"
30 #include "ConfigInterface.hh"
31
32 #include "error_handling.hh"
33
34 namespace pism {
35
36 // this file contains method for derived class IceModelVec3
37
38 // methods for base class IceModelVec and derived class IceModelVec2S
39 // are in "iceModelVec.cc"
40
41 IceModelVec3D::IceModelVec3D()
42 : IceModelVec() {
43 m_bsearch_accel = gsl_interp_accel_alloc();
44 if (m_bsearch_accel == NULL) {
45 throw RuntimeError(PISM_ERROR_LOCATION, "Failed to allocate a GSL interpolation accelerator");
46 }
47 }
48
49 IceModelVec3D::~IceModelVec3D() {
50 gsl_interp_accel_free(m_bsearch_accel);
51 }
52
53 IceModelVec3::IceModelVec3() {
54 // empty
55 }
56
57 IceModelVec3::IceModelVec3(IceGrid::ConstPtr grid, const std::string &short_name,
58 IceModelVecKind ghostedp,
59 unsigned int stencil_width) {
60 create(grid, short_name, ghostedp, stencil_width);
61 }
62
63 IceModelVec3::~IceModelVec3() {
64 // empty
65 }
66
67 IceModelVec3::Ptr IceModelVec3::To3DScalar(IceModelVec::Ptr input) {
68 IceModelVec3::Ptr result = dynamic_pointer_cast<IceModelVec3,IceModelVec>(input);
69 if (not (bool)result) {
70 throw RuntimeError(PISM_ERROR_LOCATION, "dynamic cast failure");
71 }
72 return result;
73 }
74
75 //! Allocate a DA and a Vec from information in IceGrid.
76 void IceModelVec3D::allocate(IceGrid::ConstPtr grid, const std::string &name,
77 IceModelVecKind ghostedp, const std::vector<double> &levels,
78 unsigned int stencil_width) {
79 PetscErrorCode ierr;
80 m_grid = grid;
81
82 m_zlevels = levels;
83 m_da_stencil_width = stencil_width;
84
85 m_da = m_grid->get_dm(this->m_zlevels.size(), this->m_da_stencil_width);
86
87 m_has_ghosts = (ghostedp == WITH_GHOSTS);
88
89 if (m_has_ghosts == true) {
90 ierr = DMCreateLocalVector(*m_da, m_v.rawptr());
91 PISM_CHK(ierr, "DMCreateLocalVector");
92 } else {
93 ierr = DMCreateGlobalVector(*m_da, m_v.rawptr());
94 PISM_CHK(ierr, "DMCreateGlobalVector");
95 }
96
97 m_name = name;
98
99 m_metadata.push_back(SpatialVariableMetadata(m_grid->ctx()->unit_system(),
100 name, m_zlevels));
101 }
102
103 bool IceModelVec3D::isLegalLevel(double z) const {
104 double z_min = m_zlevels.front(),
105 z_max = m_zlevels.back();
106 if (z < z_min - 1.0e-6 || z > z_max + 1.0e-6) {
107 return false;
108 }
109 return true;
110 }
111
112 //! Set all values of scalar quantity to given a single value in a particular column.
113 void IceModelVec3D::set_column(int i, int j, double c) {
114 PetscErrorCode ierr;
115 #if (Pism_DEBUG==1)
116 assert(m_array != NULL);
117 check_array_indices(i, j, 0);
118 #endif
119
120 double ***arr = (double***) m_array;
121
122 if (c == 0.0) {
123 ierr = PetscMemzero(arr[j][i], m_zlevels.size() * sizeof(double));
124 PISM_CHK(ierr, "PetscMemzero");
125 } else {
126 unsigned int nlevels = m_zlevels.size();
127 for (unsigned int k=0; k < nlevels; k++) {
128 arr[j][i][k] = c;
129 }
130 }
131 }
132
133 void IceModelVec3D::set_column(int i, int j, const std::vector<double> &data) {
134 double *column = get_column(i, j);
135 for (unsigned int k = 0; k < data.size(); ++k) {
136 column[k] = data[k];
137 }
138
139 }
140
141 const std::vector<double> IceModelVec3D::get_column_vector(int i, int j) const {
142 unsigned int n = m_zlevels.size();
143 std::vector<double> result(n);
144 const double *data = get_column(i, j);
145 for (unsigned int k = 0; k < n; ++k) {
146 result[k] = data[k];
147 }
148 return result;
149 }
150
151
152 //! Return value of scalar quantity at level z (m) above base of ice (by linear interpolation).
153 double IceModelVec3D::getValZ(int i, int j, double z) const {
154 #if (Pism_DEBUG==1)
155 assert(m_array != NULL);
156 check_array_indices(i, j, 0);
157
158 if (not isLegalLevel(z)) {
159 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "IceModelVec3 getValZ(): level %f is not legal; name = %s",
160 z, m_name.c_str());
161 }
162 #endif
163
164 double ***arr = (double***) m_array;
165 if (z >= m_zlevels.back()) {
166 unsigned int nlevels = m_zlevels.size();
167 return arr[j][i][nlevels - 1];
168 } else if (z <= m_zlevels.front()) {
169 return arr[j][i][0];
170 }
171
172 unsigned int mcurr = gsl_interp_accel_find(m_bsearch_accel, &m_zlevels[0], m_zlevels.size(), z);
173
174 const double incr = (z - m_zlevels[mcurr]) / (m_zlevels[mcurr+1] - m_zlevels[mcurr]);
175 const double valm = arr[j][i][mcurr];
176 return valm + incr * (arr[j][i][mcurr+1] - valm);
177 }
178
179 //! Copies a horizontal slice at level z of an IceModelVec3 into a Vec gslice.
180 /*!
181 * FIXME: this method is misnamed: the slice is horizontal in the PISM
182 * coordinate system, not in reality.
183 */
184 void IceModelVec3::getHorSlice(Vec &gslice, double z) const {
185
186 petsc::DM::Ptr da2 = m_grid->get_dm(1, m_grid->ctx()->config()->get_number("grid.max_stencil_width"));
187
188 IceModelVec::AccessList list(*this);
189 petsc::DMDAVecArray slice(da2, gslice);
190 double **slice_val = (double**)slice.get();
191
192 ParallelSection loop(m_grid->com);
193 try {
194 for (Points p(*m_grid); p; p.next()) {
195 const int i = p.i(), j = p.j();
196 slice_val[j][i] = getValZ(i,j,z);
197 }
198 } catch (...) {
199 loop.failed();
200 }
201 loop.check();
202 }
203
204 //! Copies a horizontal slice at level z of an IceModelVec3 into an IceModelVec2S gslice.
205 /*!
206 * FIXME: this method is misnamed: the slice is horizontal in the PISM
207 * coordinate system, not in reality.
208 */
209 void IceModelVec3::getHorSlice(IceModelVec2S &gslice, double z) const {
210 IceModelVec::AccessList list{this, &gslice};
211
212 ParallelSection loop(m_grid->com);
213 try {
214 for (Points p(*m_grid); p; p.next()) {
215 const int i = p.i(), j = p.j();
216 gslice(i, j) = getValZ(i, j, z);
217 }
218 } catch (...) {
219 loop.failed();
220 }
221 loop.check();
222 }
223
224
225 //! Copies the values of an IceModelVec3 at the ice surface (specified by the level myH) to an IceModelVec2S gsurf.
226 void IceModelVec3::getSurfaceValues(IceModelVec2S &surface_values,
227 const IceModelVec2S &H) const {
228 IceModelVec::AccessList list{this, &surface_values, &H};
229
230 ParallelSection loop(m_grid->com);
231 try {
232 for (Points p(*m_grid); p; p.next()) {
233 const int i = p.i(), j = p.j();
234 surface_values(i, j) = getValZ(i, j, H(i, j));
235 }
236 } catch (...) {
237 loop.failed();
238 }
239 loop.check();
240 }
241
242 double* IceModelVec3D::get_column(int i, int j) {
243 #if (Pism_DEBUG==1)
244 check_array_indices(i, j, 0);
245 #endif
246 return ((double***) m_array)[j][i];
247 }
248
249 const double* IceModelVec3D::get_column(int i, int j) const {
250 #if (Pism_DEBUG==1)
251 check_array_indices(i, j, 0);
252 #endif
253 return ((double***) m_array)[j][i];
254 }
255
256 void IceModelVec3D::set_column(int i, int j, const double *valsIN) {
257 #if (Pism_DEBUG==1)
258 check_array_indices(i, j, 0);
259 #endif
260 double ***arr = (double***) m_array;
261 PetscErrorCode ierr = PetscMemcpy(arr[j][i], valsIN, m_zlevels.size()*sizeof(double));
262 PISM_CHK(ierr, "PetscMemcpy");
263 }
264
265 void IceModelVec3::create(IceGrid::ConstPtr grid, const std::string &name,
266 IceModelVecKind ghostedp,
267 unsigned int stencil_width) {
268
269 IceModelVec3D::allocate(grid, name, ghostedp,
270 grid->z(), stencil_width);
271 }
272
273 /** Sum a 3-D vector in the Z direction to create a 2-D vector.
274
275 Note that this sums up all the values in a column, including ones
276 above the ice. This may or may not be what you need. Also, take a look
277 at IceModel::compute_ice_enthalpy(PetscScalar &result) in iMreport.cc.
278
279 As for the difference between IceModelVec2 and IceModelVec2S, the
280 former can store fields with more than 1 "degree of freedom" per grid
281 point (such as 2D fields on the "staggered" grid, with the first
282 degree of freedom corresponding to the i-offset and second to
283 j-offset).
284
285 IceModelVec2S is just IceModelVec2 with "dof == 1", and
286 IceModelVec2V is IceModelVec2 with "dof == 2". (Plus some extra
287 methods, of course.)
288
289 Either one of IceModelVec2 and IceModelVec2S would work in this
290 case.
291
292 Computes output = A*output + B*sum_columns(input) + C
293
294 @see https://github.com/pism/pism/issues/229 */
295 void IceModelVec3::sumColumns(IceModelVec2S &output, double A, double B) const {
296 const unsigned int Mz = m_grid->Mz();
297
298 AccessList access{this, &output};
299
300 for (Points p(*m_grid); p; p.next()) {
301 const int i = p.i(), j = p.j();
302
303 const double *column = this->get_column(i,j);
304
305 double scalar_sum = 0.0;
306 for (unsigned int k = 0; k < Mz; ++k) {
307 scalar_sum += column[k];
308 }
309 output(i,j) = A * output(i,j) + B * scalar_sum;
310 }
311 }
312
313 } // end of namespace pism