Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.1-cd005eec8 committed by Constantine Khrulev on 2025-03-07
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
siafd_test.cc
Go to the documentation of this file.
1// Copyright (C) 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024 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
19static char help[] =
20 "\nSIAFD_TEST\n"
21 " Testing program for SIA, time-independent calculations separate from\n"
22 " IceModel. Uses verification test F. Also may be used in a PISM software"
23 "(regression) test.\n\n";
24
25#include "pism/stressbalance/sia/SIAFD.hh"
26#include "pism/util/EnthalpyConverter.hh"
27#include "pism/stressbalance/StressBalance.hh"
28#include "pism/stressbalance/SSB_Modifier.hh"
29#include "pism/stressbalance/ShallowStressBalance.hh"
30#include "pism/util/Grid.hh"
31#include "pism/util/Mask.hh"
32#include "pism/util/Context.hh"
33#include "pism/util/Time.hh"
34#include "pism/util/error_handling.hh"
35#include "pism/util/io/File.hh"
36#include "pism/util/petscwrappers/PetscInitializer.hh"
37#include "pism/util/pism_utilities.hh"
38#include "pism/util/pism_options.hh"
39#include "pism/verification/tests/exactTestsFG.hh"
40#include "pism/util/io/io_helpers.hh"
41#include "pism/geometry/Geometry.hh"
42
43namespace pism {
44
45static void compute_strain_heating_errors(const array::Array3D &strain_heating,
46 const array::Scalar &thickness,
47 const Grid &grid,
48 double &gmax_strain_heating_err,
49 double &gav_strain_heating_err) {
50 double max_strain_heating_error = 0.0, av_strain_heating_error = 0.0, avcount = 0.0;
51 const int Mz = grid.Mz();
52
53 const double LforFG = 750000; // m
54
55 const double
56 ice_rho = grid.ctx()->config()->get_number("constants.ice.density"),
57 ice_c = grid.ctx()->config()->get_number("constants.ice.specific_heat_capacity");
58
59 array::AccessScope list{&thickness, &strain_heating};
60
61 ParallelSection loop(grid.com);
62 try {
63 for (auto p = grid.points(); p; p.next()) {
64 const int i = p.i(), j = p.j();
65
66 double
67 xx = grid.x(i),
68 yy = grid.y(j),
69 r = sqrt(pow(xx, 2) + pow(yy, 2));
70
71 if ((r >= 1.0) && (r <= LforFG - 1.0)) {
72 // only evaluate error if inside sheet and not at central
73 // singularity
74 TestFGParameters F = exactFG(0.0, r, grid.z(), 0.0);
75
76 for (int k = 0; k < Mz; k++) {
77 F.Sig[k] *= ice_rho * ice_c; // scale exact strain_heating to J/(s m^3)
78 }
79 const int ks = grid.kBelowHeight(thickness(i, j));
80 const double *strain_heating_ij = strain_heating.get_column(i, j);
81 for (int k = 0; k < ks; k++) { // only eval error if below num surface
82 const double _strain_heating_error = fabs(strain_heating_ij[k] - F.Sig[k]);
83 max_strain_heating_error = std::max(max_strain_heating_error, _strain_heating_error);
84 avcount += 1.0;
85 av_strain_heating_error += _strain_heating_error;
86 }
87 }
88 }
89 } catch (...) {
90 loop.failed();
91 }
92 loop.check();
93
94 gmax_strain_heating_err = GlobalMax(grid.com, max_strain_heating_error);
95 gav_strain_heating_err = GlobalSum(grid.com, av_strain_heating_error);
96 double gavcount = GlobalSum(grid.com, avcount);
97 gav_strain_heating_err = gav_strain_heating_err/std::max(gavcount,1.0); // avoid div by zero
98}
99
100
101static void computeSurfaceVelocityErrors(const Grid &grid,
102 const array::Scalar &ice_thickness,
103 const array::Array3D &u3,
104 const array::Array3D &v3,
105 const array::Array3D &w3,
106 double &gmaxUerr,
107 double &gavUerr,
108 double &gmaxWerr,
109 double &gavWerr) {
110 double maxUerr = 0.0, maxWerr = 0.0, avUerr = 0.0, avWerr = 0.0;
111
112 const double LforFG = 750000; // m
113
114 array::AccessScope list{&ice_thickness, &u3, &v3, &w3};
115
116 for (auto p = grid.points(); p; p.next()) {
117 const int i = p.i(), j = p.j();
118
119 double xx = grid.x(i), yy = grid.y(j),
120 r = sqrt(pow(xx, 2) + pow(yy, 2));
121 if ((r >= 1.0) && (r <= LforFG - 1.0)) { // only evaluate error if inside sheet
122 // and not at central singularity
123 const double H = ice_thickness(i, j);
124 std::vector<double> z(1, H);
125 TestFGParameters F = exactFG(0.0, r, z, 0.0);
126
127 const double uex = (xx/r) * F.U[0];
128 const double vex = (yy/r) * F.U[0];
129 // note that because interpolate does linear interpolation and H is not exactly at
130 // a grid point, this causes nonzero errors
131 const double Uerr = sqrt(pow(u3.interpolate(i,j,H) - uex, 2) +
132 pow(v3.interpolate(i,j,H) - vex, 2));
133 maxUerr = std::max(maxUerr,Uerr);
134 avUerr += Uerr;
135 const double Werr = fabs(w3.interpolate(i,j,H) - F.w[0]);
136 maxWerr = std::max(maxWerr,Werr);
137 avWerr += Werr;
138 }
139 }
140
141 gmaxUerr = GlobalMax(grid.com, maxUerr);
142 gmaxWerr = GlobalMax(grid.com, maxWerr);
143 gavUerr = GlobalSum(grid.com, avUerr);
144 gavUerr = gavUerr/(grid.Mx()*grid.My());
145 gavWerr = GlobalSum(grid.com, avWerr);
146 gavWerr = gavWerr/(grid.Mx()*grid.My());
147}
148
149
151 const Grid &grid,
152 const array::Scalar &thickness,
153 const array::Array3D &temperature,
154 array::Array3D &enthalpy) {
155
156 array::AccessScope list{&temperature, &enthalpy, &thickness};
157
158 for (auto p = grid.points(); p; p.next()) {
159 const int i = p.i(), j = p.j();
160
161 const double *T_ij = temperature.get_column(i,j);
162 double *E_ij = enthalpy.get_column(i,j);
163
164 for (unsigned int k=0; k<grid.Mz(); ++k) {
165 double depth = thickness(i,j) - grid.z(k);
166 E_ij[k] = EC.enthalpy_permissive(T_ij[k], 0.0,
167 EC.pressure(depth));
168 }
169
170 }
171
172 enthalpy.update_ghosts();
173}
174
175
176//! \brief Set the test F initial state.
177static void setInitStateF(Grid &grid,
179 array::Scalar &bed,
180 array::Scalar &mask,
181 array::Scalar &surface,
182 array::Scalar &thickness,
183 array::Array3D &enthalpy) {
184
185 double
186 ST = 1.67e-5,
187 Tmin = 223.15, // K
188 LforFG = 750000; // m
189
190 bed.set(0.0);
191 mask.set(MASK_GROUNDED);
192
193 array::AccessScope list{&thickness, &enthalpy};
194
195 for (auto p = grid.points(); p; p.next()) {
196 const int i = p.i(), j = p.j();
197
198 const double
199 r = std::max(grid::radius(grid, i, j), 1.0), // avoid singularity at origin
200 Ts = Tmin + ST * r;
201
202 if (r > LforFG - 1.0) { // if (essentially) outside of sheet
203 thickness(i, j) = 0.0;
204 enthalpy.set_column(i, j, Ts);
205 } else {
206 TestFGParameters F = exactFG(0.0, r, grid.z(), 0.0);
207
208 thickness(i, j) = F.H;
209 enthalpy.set_column(i, j, (F.T).data());
210 }
211 }
212
213
214 thickness.update_ghosts();
215
216 surface.copy_from(thickness);
217
218 enthalpy_from_temperature_cold(EC, grid, thickness, enthalpy,
219 enthalpy);
220}
221
222static void reportErrors(const Grid &grid,
223 units::System::Ptr unit_system,
224 const array::Scalar &thickness,
225 const array::Array3D &u_sia,
226 const array::Array3D &v_sia,
227 const array::Array3D &w_sia,
228 const array::Array3D &strain_heating) {
229
230 Logger::ConstPtr log = grid.ctx()->log();
231
232 // strain_heating errors if appropriate; reported in 10^6 J/(s m^3)
233 double max_strain_heating_error, av_strain_heating_error;
234 compute_strain_heating_errors(strain_heating, thickness,
235 grid,
236 max_strain_heating_error, av_strain_heating_error);
237
238 log->message(1,
239 "Sigma : maxSig avSig\n");
240 log->message(1, " %12.6f%12.6f\n",
241 max_strain_heating_error*1.0e6, av_strain_heating_error*1.0e6);
242
243 // surface velocity errors if exact values are available; reported in m year-1
244 double maxUerr, avUerr, maxWerr, avWerr;
245 computeSurfaceVelocityErrors(grid, thickness,
246 u_sia,
247 v_sia,
248 w_sia,
249 maxUerr, avUerr,
250 maxWerr, avWerr);
251
252 log->message(1,
253 "surf vels : maxUvec avUvec maxW avW\n");
254 log->message(1, " %12.6f%12.6f%12.6f%12.6f\n",
255 units::convert(unit_system, maxUerr, "m second^-1", "m year^-1"),
256 units::convert(unit_system, avUerr, "m second^-1", "m year^-1"),
257 units::convert(unit_system, maxWerr, "m second^-1", "m year^-1"),
258 units::convert(unit_system, avWerr, "m second^-1", "m year^-1"));
259}
260
261} // end of namespace pism
262
263int main(int argc, char *argv[]) {
264
265 using namespace pism;
266 using namespace pism::stressbalance;
267
268 MPI_Comm com = MPI_COMM_WORLD;
269 petsc::Initializer petsc(argc, argv, help);
270
271 com = MPI_COMM_WORLD;
272
273 /* This explicit scoping forces destructors to be called before PetscFinalize() */
274 try {
275 // set default verbosity
276 std::shared_ptr<Context> ctx = context_from_options(com, "siafd_test");
277 Config::Ptr config = ctx->config();
278
279 config->set_flag("stress_balance.sia.grain_size_age_coupling", false);
280 config->set_string("stress_balance.sia.flow_law", "arr");
281
282 set_config_from_options(ctx->unit_system(), *config);
283 config->resolve_filenames();
284
285 std::string usage = "\n"
286 "usage of SIAFD_TEST:\n"
287 " run siafd_test -Mx <number> -My <number> -Mz <number> -o foo.nc\n"
288 "\n";
289
290 bool stop = show_usage_check_req_opts(*ctx->log(), "siafd_test", {}, usage);
291
292 if (stop) {
293 return 0;
294 }
295
296 auto output_file = config->get_string("output.file");
297
298 grid::Parameters P(*config);
299 P.Lx = 900e3;
300 P.Ly = P.Lx;
301
302 double Lz = 4000.0;
303 unsigned int Mz = config->get_number("grid.Mz");
304
305 P.z = grid::compute_vertical_levels(Lz, Mz, grid::EQUAL);
306 P.ownership_ranges_from_options(*ctx->config(), ctx->size());
307
308 // create grid and set defaults
309 auto grid = std::make_shared<Grid>(ctx, P);
310 grid->report_parameters();
311
313
314 const int WIDE_STENCIL = config->get_number("grid.max_stencil_width");
315
317 enthalpy(grid, "enthalpy", array::WITH_GHOSTS, grid->z(), WIDE_STENCIL),
318 age(grid, "age", array::WITHOUT_GHOSTS, grid->z());
319
320 Geometry geometry(grid);
321 geometry.sea_level_elevation.set(0.0);
322
323 // age of the ice; is not used here
324 age.metadata(0).long_name("age of the ice").units("s");
325 age.set(0.0);
326
327 // enthalpy in the ice
328 enthalpy.metadata(0)
329 .long_name("ice enthalpy (includes sensible heat, latent heat, pressure)")
330 .units("J kg^-1");
331 //
332 enthalpy.set(EC->enthalpy(263.15, 0.0, EC->pressure(1000.0)));
333
334 // Create the SIA solver object:
335
336 // We use SIA_Nonsliding and not SIAFD here because we need the z-component
337 // of the ice velocity, which is computed using incompressibility of ice in
338 // StressBalance::compute_vertical_velocity().
339 std::shared_ptr<SIAFD> sia(new SIAFD(grid));
340 std::shared_ptr<ZeroSliding> no_sliding(new ZeroSliding(grid));
341
342 // stress_balance will de-allocate no_sliding and sia.
343 StressBalance stress_balance(grid, no_sliding, sia);
344
345 // fill the fields:
346 setInitStateF(*grid, *EC,
347 geometry.bed_elevation,
348 geometry.cell_type,
349 geometry.ice_surface_elevation,
350 geometry.ice_thickness,
351 enthalpy);
352
353 geometry.ensure_consistency(config->get_number("geometry.ice_free_thickness_standard"));
354
355 // Initialize the SIA solver:
356 stress_balance.init();
357
358 bool full_update = true;
359
361 inputs.geometry = &geometry;
362 inputs.water_column_pressure = nullptr;
363 inputs.enthalpy = &enthalpy;
364 inputs.age = &age;
365
366 stress_balance.update(inputs, full_update);
367
368 // Report errors relative to the exact solution:
369 const array::Array3D
370 &u3 = stress_balance.velocity_u(),
371 &v3 = stress_balance.velocity_v(),
372 &w3 = stress_balance.velocity_w();
373
374 const array::Array3D &sigma = stress_balance.volumetric_strain_heating();
375
376 reportErrors(*grid, ctx->unit_system(),
377 geometry.ice_thickness, u3, v3, w3, sigma);
378
379 // Write results to an output file:
380 File file(grid->com, output_file, io::PISM_NETCDF3, io::PISM_READWRITE_MOVE);
381 io::define_time(file, *ctx);
382 io::append_time(file, *ctx->config(), ctx->time()->current());
383
384 geometry.ice_surface_elevation.write(file);
385 geometry.ice_thickness.write(file);
386 geometry.cell_type.write(file);
387 geometry.bed_elevation.write(file);
388
389 sia->diffusivity().write(file);
390 u3.write(file);
391 v3.write(file);
392 w3.write(file);
393 sigma.write(file);
394 }
395 catch (...) {
396 handle_fatal_errors(com);
397 return 1;
398 }
399
400 return 0;
401}
An EnthalpyConverter for use in verification tests.
std::shared_ptr< Config > Ptr
std::shared_ptr< EnthalpyConverter > Ptr
double enthalpy_permissive(double T, double omega, double P) const
Compute enthalpy more permissively than enthalpy().
double pressure(double depth) const
Get pressure in ice from depth below surface using the hydrostatic assumption.
Converts between specific enthalpy and temperature or liquid content.
High-level PISM I/O class.
Definition File.hh:55
array::Scalar1 sea_level_elevation
Definition Geometry.hh:48
void ensure_consistency(double ice_free_thickness_threshold)
Definition Geometry.cc:114
array::Scalar2 ice_surface_elevation
Definition Geometry.hh:57
array::CellType2 cell_type
Definition Geometry.hh:55
array::Scalar2 ice_thickness
Definition Geometry.hh:51
array::Scalar2 bed_elevation
Definition Geometry.hh:47
Describes the PISM grid and the distribution of data across processors.
Definition Grid.hh:290
std::shared_ptr< const Logger > ConstPtr
Definition Logger.hh:46
void failed()
Indicates a failure of a parallel section.
VariableMetadata & long_name(const std::string &input)
VariableMetadata & units(const std::string &input)
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition Array.hh:64
void copy_from(const Array2D< T > &source)
Definition Array2D.hh:73
double interpolate(int i, int j, double z) const
Return value of scalar quantity at level z (m) above base of ice (by linear interpolation).
Definition Array3D.cc:89
void set_column(int i, int j, double c)
Set all values of scalar quantity to given a single value in a particular column.
Definition Array3D.cc:50
double * get_column(int i, int j)
Definition Array3D.cc:121
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition Array3D.hh:33
void write(const std::string &filename) const
Definition Array.cc:722
void set(double c)
Result: v[j] <- c for all j.
Definition Array.cc:629
void update_ghosts()
Updates ghost points.
Definition Array.cc:615
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition Array.cc:476
std::vector< double > z
Vertical levels.
Definition Grid.hh:165
void ownership_ranges_from_options(const Config &config, unsigned int size)
Re-compute ownership ranges. Uses current values of Mx and My.
Definition Grid.cc:1192
double Ly
Domain half-width in the Y direction.
Definition Grid.hh:151
double Lx
Domain half-width in the X direction.
Definition Grid.hh:149
Grid parameters; used to collect defaults before an Grid is allocated.
Definition Grid.hh:123
const array::Array3D * age
const array::Scalar * water_column_pressure
const array::Array3D * enthalpy
const array::Array3D & velocity_w() const
const array::Array3D & velocity_u() const
Get components of the the 3D velocity field.
void init()
Initialize the StressBalance object.
void update(const Inputs &inputs, bool full_update)
Update all the fields if (full_update), only update diffusive flux and max. diffusivity otherwise.
const array::Array3D & velocity_v() const
const array::Array3D & volumetric_strain_heating() const
The class defining PISM's interface to the shallow stress balance code.
Returns zero velocity field, zero friction heating, and zero for D^2.
std::shared_ptr< System > Ptr
Definition Units.hh:47
const double Ts
Definition exactTestK.c:39
double radius(const Grid &grid, int i, int j)
Returns the distance from the point (i,j) to the origin.
Definition Grid.cc:1377
Stress balance models and related diagnostics.
Definition Isochrones.hh:31
double convert(System::Ptr system, double input, const std::string &spec1, const std::string &spec2)
Convert a quantity from unit1 to unit2.
Definition Units.cc:70
static double F(double SL, double B, double H, double alpha)
static void compute_strain_heating_errors(const array::Array3D &strain_heating, const array::Scalar &thickness, const Grid &grid, double &gmax_strain_heating_err, double &gav_strain_heating_err)
Definition siafd_test.cc:45
static void computeSurfaceVelocityErrors(const Grid &grid, const array::Scalar &ice_thickness, const array::Array3D &u3, const array::Array3D &v3, const array::Array3D &w3, double &gmaxUerr, double &gavUerr, double &gmaxWerr, double &gavWerr)
static void enthalpy_from_temperature_cold(EnthalpyConverter &EC, const Grid &grid, const array::Scalar &thickness, const array::Array3D &temperature, array::Array3D &enthalpy)
void GlobalMax(MPI_Comm comm, double *local, double *result, int count)
static const double k
Definition exactTestP.cc:42
TestFGParameters exactFG(double t, double r, const std::vector< double > &z, double Cp)
@ MASK_GROUNDED
Definition Mask.hh:33
static void reportErrors(const Grid &grid, units::System::Ptr unit_system, const array::Scalar &thickness, const array::Array3D &u_sia, const array::Array3D &v_sia, const array::Array3D &w_sia, const array::Array3D &strain_heating)
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
static void setInitStateF(Grid &grid, EnthalpyConverter &EC, array::Scalar &bed, array::Scalar &mask, array::Scalar &surface, array::Scalar &thickness, array::Array3D &enthalpy)
Set the test F initial state.
int main(int argc, char *argv[])
static char help[]
Definition siafd_test.cc:19