tSSATestCase.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tSSATestCase.cc (13367B)
---
1 // Copyright (C) 2009--2019 Ed Bueler, Constantine Khroulev, and David Maxwell
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 "SSATestCase.hh"
20 #include "SSAFD.hh"
21 #include "SSAFEM.hh"
22 #include "pism/util/Mask.hh"
23 #include "pism/util/Time.hh"
24 #include "pism/util/io/File.hh"
25 #include "pism/util/pism_options.hh"
26 #include "pism/util/io/io_helpers.hh"
27 #include "pism/util/pism_utilities.hh"
28 #include "pism/stressbalance/StressBalance.hh"
29
30 namespace pism {
31 namespace stressbalance {
32
33 SSATestCase::SSATestCase(Context::Ptr ctx, int Mx, int My,
34 double Lx, double Ly,
35 GridRegistration registration,
36 Periodicity periodicity)
37 : m_com(ctx->com()), m_ctx(ctx), m_config(ctx->config()),
38 m_grid(IceGrid::Shallow(m_ctx, Lx, Ly, 0.0, 0.0, Mx, My, registration, periodicity)),
39 m_sys(ctx->unit_system()),
40 m_geometry(m_grid),
41 m_ssa(NULL) {
42
43 const unsigned int WIDE_STENCIL = m_config->get_number("grid.max_stencil_width");
44
45 // yield stress for basal till (plastic or pseudo-plastic model)
46 m_tauc.create(m_grid, "tauc", WITH_GHOSTS, WIDE_STENCIL);
47 m_tauc.set_attrs("diagnostic",
48 "yield stress for basal till (plastic or pseudo-plastic model)",
49 "Pa", "Pa", "", 0);
50
51 // enthalpy
52 m_ice_enthalpy.create(m_grid, "enthalpy", WITH_GHOSTS, WIDE_STENCIL);
53 m_ice_enthalpy.set_attrs("model_state",
54 "ice enthalpy (includes sensible heat, latent heat, pressure)",
55 "J kg-1", "J kg-1", "", 0);
56
57 // dirichlet boundary condition (FIXME: perhaps unused!)
58 m_bc_values.create(m_grid, "_bc", WITH_GHOSTS, WIDE_STENCIL); // u_bc and v_bc
59 m_bc_values.set_attrs("intent",
60 "X-component of the SSA velocity boundary conditions",
61 "m s-1", "m year-1", "", 0);
62 m_bc_values.set_attrs("intent",
63 "Y-component of the SSA velocity boundary conditions",
64 "m s-1", "m year-1", "", 1);
65
66 Config::ConstPtr config = m_grid->ctx()->config();
67 units::System::Ptr sys = m_grid->ctx()->unit_system();
68 double fill_value = units::convert(sys, config->get_number("output.fill_value"), "m year-1", "m second-1");
69
70 auto large_number = units::convert(m_sys, 1e6, "m year-1", "m second-1");
71
72 m_bc_values.metadata(0).set_numbers("valid_range", {-large_number, large_number});
73 m_bc_values.metadata(0).set_number("_FillValue", fill_value);
74
75 m_bc_values.metadata(1).set_numbers("valid_range", {-large_number, large_number});
76 m_bc_values.metadata(1).set_number("_FillValue", fill_value);
77
78 m_bc_values.set(fill_value);
79
80 // Dirichlet B.C. mask
81 m_bc_mask.create(m_grid, "bc_mask", WITH_GHOSTS, WIDE_STENCIL);
82 m_bc_mask.set_attrs("model_state",
83 "grounded_dragging_floating integer mask",
84 "", "", "", 0);
85
86 m_bc_mask.metadata().set_numbers("flag_values", {0.0, 1.0});
87 m_bc_mask.metadata().set_string("flag_meanings",
88 "no_data ssa.dirichlet_bc_location");
89
90 m_melange_back_pressure.create(m_grid, "melange_back_pressure_fraction",
91 WITH_GHOSTS, WIDE_STENCIL);
92 m_melange_back_pressure.set_attrs("boundary_condition",
93 "melange back pressure fraction",
94 "", "", "", 0);
95 m_melange_back_pressure.set(0.0);
96 }
97
98 SSATestCase::~SSATestCase()
99 {
100 delete m_ssa;
101 }
102
103 //! Initialize the test case at the start of a run
104 void SSATestCase::init() {
105
106 m_ssa->init();
107
108 // Allow the subclass to set the coefficients.
109 initializeSSACoefficients();
110 }
111
112 //! Solve the SSA
113 void SSATestCase::run() {
114 m_ctx->log()->message(2, "* Solving the SSA stress balance ...\n");
115
116 m_geometry.ensure_consistency(m_config->get_number("stress_balance.ice_free_thickness_standard"));
117
118 Inputs inputs;
119 inputs.melange_back_pressure = &m_melange_back_pressure;
120 inputs.geometry = &m_geometry;
121 inputs.enthalpy = &m_ice_enthalpy;
122 inputs.basal_yield_stress = &m_tauc;
123 inputs.bc_mask = &m_bc_mask;
124 inputs.bc_values = &m_bc_values;
125
126 bool full_update = true;
127 m_ssa->update(inputs, full_update);
128 }
129
130 //! Report on the generated solution
131 void SSATestCase::report(const std::string &testname) {
132
133 m_ctx->log()->message(3, m_ssa->stdout_report());
134
135 double
136 maxvecerr = 0.0,
137 avvecerr = 0.0,
138 avuerr = 0.0,
139 avverr = 0.0,
140 maxuerr = 0.0,
141 maxverr = 0.0;
142 double
143 gmaxvecerr = 0.0,
144 gavvecerr = 0.0,
145 gavuerr = 0.0,
146 gavverr = 0.0,
147 gmaxuerr = 0.0,
148 gmaxverr = 0.0;
149
150 if (m_config->get_flag("basal_resistance.pseudo_plastic.enabled") &&
151 m_config->get_number("basal_resistance.pseudo_plastic.q") != 1.0) {
152 m_ctx->log()->message(1,
153 "WARNING: numerical errors not valid for pseudo-plastic till\n");
154 }
155 m_ctx->log()->message(1,
156 "NUMERICAL ERRORS in velocity relative to exact solution:\n");
157
158 const IceModelVec2V &vel_ssa = m_ssa->velocity();
159
160 IceModelVec::AccessList list{&vel_ssa};
161
162 double exactvelmax = 0, gexactvelmax = 0;
163 for (Points p(*m_grid); p; p.next()) {
164 const int i = p.i(), j = p.j();
165
166 double uexact, vexact;
167 double myx = m_grid->x(i), myy = m_grid->y(j);
168
169 exactSolution(i,j,myx,myy,&uexact,&vexact);
170
171 double exactnormsq=sqrt(uexact*uexact+vexact*vexact);
172 exactvelmax = std::max(exactnormsq,exactvelmax);
173
174 // compute maximum errors
175 const double uerr = fabs(vel_ssa(i,j).u - uexact);
176 const double verr = fabs(vel_ssa(i,j).v - vexact);
177 avuerr = avuerr + uerr;
178 avverr = avverr + verr;
179 maxuerr = std::max(maxuerr,uerr);
180 maxverr = std::max(maxverr,verr);
181 const double vecerr = sqrt(uerr * uerr + verr * verr);
182 maxvecerr = std::max(maxvecerr,vecerr);
183 avvecerr = avvecerr + vecerr;
184 }
185
186 unsigned int N = (m_grid->Mx()*m_grid->My());
187
188 gexactvelmax = GlobalMax(m_grid->com, exactvelmax);
189 gmaxuerr = GlobalMax(m_grid->com, maxuerr);
190 gmaxverr = GlobalMax(m_grid->com, maxverr);
191 gavuerr = GlobalSum(m_grid->com, avuerr);
192 gavuerr = gavuerr / N;
193 gavverr = GlobalSum(m_grid->com, avverr);
194 gavverr = gavverr / N;
195 gmaxvecerr = GlobalMax(m_grid->com, maxvecerr);
196 gavvecerr = GlobalSum(m_grid->com, avvecerr);
197 gavvecerr = gavvecerr / N;
198
199 using pism::units::convert;
200
201 m_ctx->log()->message(1,
202 "velocity : maxvector prcntavvec maxu maxv avu avv\n");
203 m_ctx->log()->message(1,
204 " %11.4f%13.5f%10.4f%10.4f%10.4f%10.4f\n",
205 convert(m_sys, gmaxvecerr, "m second-1", "m year-1"),
206 (gavvecerr/gexactvelmax)*100.0,
207 convert(m_sys, gmaxuerr, "m second-1", "m year-1"),
208 convert(m_sys, gmaxverr, "m second-1", "m year-1"),
209 convert(m_sys, gavuerr, "m second-1", "m year-1"),
210 convert(m_sys, gavverr, "m second-1", "m year-1"));
211
212 m_ctx->log()->message(1, "NUM ERRORS DONE\n");
213
214 report_netcdf(testname,
215 convert(m_sys, gmaxvecerr, "m second-1", "m year-1"),
216 (gavvecerr/gexactvelmax)*100.0,
217 convert(m_sys, gmaxuerr, "m second-1", "m year-1"),
218 convert(m_sys, gmaxverr, "m second-1", "m year-1"),
219 convert(m_sys, gavuerr, "m second-1", "m year-1"),
220 convert(m_sys, gavverr, "m second-1", "m year-1"));
221 }
222
223 void SSATestCase::report_netcdf(const std::string &testname,
224 double max_vector,
225 double rel_vector,
226 double max_u,
227 double max_v,
228 double avg_u,
229 double avg_v) {
230 TimeseriesMetadata err("N", "N", m_grid->ctx()->unit_system());
231 unsigned int start;
232 VariableMetadata global_attributes("PISM_GLOBAL", m_grid->ctx()->unit_system());
233
234 options::String filename("-report_file", "NetCDF error report file");
235
236 if (not filename.is_set()) {
237 return;
238 }
239
240 err.set_string("units", "1");
241
242 m_ctx->log()->message(2, "Also writing errors to '%s'...\n", filename->c_str());
243
244 bool append = options::Bool("-append", "Append the NetCDF error report");
245
246 IO_Mode mode = PISM_READWRITE;
247 if (not append) {
248 mode = PISM_READWRITE_MOVE;
249 }
250
251 global_attributes.set_string("source", std::string("PISM ") + pism::revision);
252
253 // Find the number of records in this file:
254 File file(m_grid->com, filename, PISM_NETCDF3, mode); // OK to use NetCDF3.
255 start = file.dimension_length("N");
256
257 io::write_attributes(file, global_attributes, PISM_DOUBLE);
258
259 // Write the dimension variable:
260 io::write_timeseries(file, err, (size_t)start, (double)(start + 1), PISM_INT);
261
262 // Always write grid parameters:
263 err.set_name("dx");
264 err.set_string("units", "meters");
265 io::write_timeseries(file, err, (size_t)start, m_grid->dx());
266 err.set_name("dy");
267 io::write_timeseries(file, err, (size_t)start, m_grid->dy());
268
269 // Always write the test name:
270 err.clear_all_strings(); err.clear_all_doubles(); err.set_string("units", "1");
271 err.set_name("test");
272 io::write_timeseries(file, err, (size_t)start, testname[0], PISM_INT);
273
274 err.clear_all_strings(); err.clear_all_doubles(); err.set_string("units", "1");
275 err.set_name("max_velocity");
276 err.set_string("units", "m year-1");
277 err.set_string("long_name", "maximum ice velocity magnitude error");
278 io::write_timeseries(file, err, (size_t)start, max_vector);
279
280 err.clear_all_strings(); err.clear_all_doubles(); err.set_string("units", "1");
281 err.set_name("relative_velocity");
282 err.set_string("units", "percent");
283 err.set_string("long_name", "relative ice velocity magnitude error");
284 io::write_timeseries(file, err, (size_t)start, rel_vector);
285
286 err.clear_all_strings(); err.clear_all_doubles(); err.set_string("units", "1");
287 err.set_name("maximum_u");
288 err.set_string("units", "m year-1");
289 err.set_string("long_name", "maximum error in the X-component of the ice velocity");
290 io::write_timeseries(file, err, (size_t)start, max_u);
291
292 err.clear_all_strings(); err.clear_all_doubles(); err.set_string("units", "1");
293 err.set_name("maximum_v");
294 err.set_string("units", "m year-1");
295 err.set_string("long_name", "maximum error in the Y-component of the ice velocity");
296 io::write_timeseries(file, err, (size_t)start, max_v);
297
298 err.clear_all_strings(); err.clear_all_doubles(); err.set_string("units", "1");
299 err.set_name("average_u");
300 err.set_string("units", "m year-1");
301 err.set_string("long_name", "average error in the X-component of the ice velocity");
302 io::write_timeseries(file, err, (size_t)start, avg_u);
303
304 err.clear_all_strings(); err.clear_all_doubles(); err.set_string("units", "1");
305 err.set_name("average_v");
306 err.set_string("units", "m year-1");
307 err.set_string("long_name", "average error in the Y-component of the ice velocity");
308 io::write_timeseries(file, err, (size_t)start, avg_v);
309
310 file.close();
311 }
312
313 void SSATestCase::exactSolution(int /*i*/, int /*j*/,
314 double /*x*/, double /*y*/,
315 double *u, double *v) {
316 *u=0; *v=0;
317 }
318
319 //! Save the computation and data to a file.
320 void SSATestCase::write(const std::string &filename) {
321
322 // Write results to an output file:
323 File file(m_grid->com, filename, PISM_NETCDF3, PISM_READWRITE_MOVE);
324 io::define_time(file, *m_grid->ctx());
325 io::append_time(file, *m_config, 0.0);
326
327 m_geometry.ice_surface_elevation.write(file);
328 m_geometry.ice_thickness.write(file);
329 m_bc_mask.write(file);
330 m_tauc.write(file);
331 m_geometry.bed_elevation.write(file);
332 m_ice_enthalpy.write(file);
333 m_bc_values.write(file);
334
335 const IceModelVec2V &vel_ssa = m_ssa->velocity();
336 vel_ssa.write(file);
337
338 IceModelVec2V exact;
339 exact.create(m_grid, "_exact", WITHOUT_GHOSTS);
340 exact.set_attrs("diagnostic",
341 "X-component of the SSA exact solution",
342 "m s-1", "m year-1", "", 0);
343 exact.set_attrs("diagnostic",
344 "Y-component of the SSA exact solution",
345 "m s-1", "m year-1", "", 1);
346
347 IceModelVec::AccessList list(exact);
348 for (Points p(*m_grid); p; p.next()) {
349 const int i = p.i(), j = p.j();
350
351 exactSolution(i, j, m_grid->x(i), m_grid->y(j),
352 &(exact(i,j).u), &(exact(i,j).v));
353 }
354 exact.write(file);
355
356 file.close();
357 }
358
359 } // end of namespace stressbalance
360 } // end of namespace pism