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