tbtutest.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
---
tbtutest.cc (8464B)
---
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 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 static char help[] =
20 "Tests BedThermalUnit using Test K, without IceModel.\n\n";
21
22 #include "pism/util/pism_options.hh"
23 #include "pism/util/IceGrid.hh"
24 #include "pism/util/io/File.hh"
25 #include "pism/util/VariableMetadata.hh"
26 #include "pism/verification/BTU_Verification.hh"
27 #include "pism/energy/BTU_Minimal.hh"
28 #include "pism/util/Time.hh"
29 #include "pism/util/Vars.hh"
30 #include "pism/util/ConfigInterface.hh"
31
32 #include "pism/verification/tests/exactTestK.h"
33
34 #include "pism/util/petscwrappers/PetscInitializer.hh"
35 #include "pism/util/error_handling.hh"
36 #include "pism/util/io/io_helpers.hh"
37 #include "pism/util/Context.hh"
38 #include "pism/util/Config.hh"
39 #include "pism/util/EnthalpyConverter.hh"
40 #include "pism/util/MaxTimestep.hh"
41 #include "pism/util/Logger.hh"
42
43 //! Allocate the PISMV (verification) context. Uses ColdEnthalpyConverter.
44 pism::Context::Ptr btutest_context(MPI_Comm com, const std::string &prefix) {
45 using namespace pism;
46
47 // unit system
48 units::System::Ptr sys(new units::System);
49
50 // logger
51 Logger::Ptr logger = logger_from_options(com);
52
53 // configuration parameters
54 Config::Ptr config = config_from_options(com, *logger, sys);
55
56 // default vertical grid parameters
57 config->set_number("grid.Mbz", 11);
58 config->set_number("grid.Lbz", 1000);
59
60 config->set_string("time.calendar", "none");
61 // when IceGrid constructor is called, these settings are used
62 config->set_number("time.start_year", 0.0);
63 config->set_number("time.run_length", 1.0);
64
65 set_config_from_options(*config);
66
67 print_config(*logger, 3, *config);
68
69 Time::Ptr time = time_from_options(com, config, sys);
70
71 EnthalpyConverter::Ptr EC = EnthalpyConverter::Ptr(new ColdEnthalpyConverter(*config));
72
73 return Context::Ptr(new Context(com, sys, config, EC, time, logger, prefix));
74 }
75
76 int main(int argc, char *argv[]) {
77
78 using namespace pism;
79
80 MPI_Comm com = MPI_COMM_WORLD;
81 petsc::Initializer petsc(argc, argv, help);
82
83 com = PETSC_COMM_WORLD;
84
85 try {
86 Context::Ptr ctx = btutest_context(com, "btutest");
87 Logger::Ptr log = ctx->log();
88
89 std::string usage =
90 " btutest -Mbz NN -Lbz 1000.0 [-o OUT.nc -ys A -ye B -dt C -Mz D -Lz E]\n"
91 "where these are required because they are used in BedThermalUnit:\n"
92 " -Mbz number of bedrock thermal layer levels to use\n"
93 " -Lbz 1000.0 depth of bedrock thermal layer (required; Lbz=1000.0 m in Test K)\n"
94 "and these are allowed:\n"
95 " -o output file name; NetCDF format\n"
96 " -ys start year in using Test K\n"
97 " -ye end year in using Test K\n"
98 " -dt time step B (= positive float) in years\n";
99
100 bool done = show_usage_check_req_opts(*log, "BTUTEST %s (test program for BedThermalUnit)",
101 {"-Mbz"}, usage);
102 if (done) {
103 return 0;
104 }
105
106 log->message(2,
107 " initializing IceGrid from options ...\n");
108
109 Config::Ptr config = ctx->config();
110
111 GridParameters P(config);
112 P.Mx = 3;
113 P.My = P.Mx;
114 P.Lx = 1500e3;
115 P.Ly = P.Lx;
116
117 P.vertical_grid_from_options(config);
118 P.ownership_ranges_from_options(ctx->size());
119
120 // create grid and set defaults
121 IceGrid::Ptr grid(new IceGrid(ctx, P));
122
123 ctx->time()->init(*log);
124
125 auto outname = config->get_string("output.file_name");
126
127 options::Real dt_years("-dt", "Time-step, in years", 1.0);
128
129 // allocate tools and IceModelVecs
130 IceModelVec2S bedtoptemp, heat_flux_at_ice_base;
131 {
132 heat_flux_at_ice_base.create(grid, "upward_heat_flux_at_ice_base", WITHOUT_GHOSTS);
133 heat_flux_at_ice_base.set_attrs("",
134 "upward geothermal flux at bedrock thermal layer base",
135 "W m-2", "mW m-2", "", 0);
136
137 bedtoptemp.create(grid, "bedtoptemp", WITHOUT_GHOSTS);
138 bedtoptemp.set_attrs("",
139 "temperature at top of bedrock thermal layer",
140 "K", "K", "", 0);
141 }
142
143 // initialize BTU object:
144 energy::BTUGrid bedrock_grid = energy::BTUGrid::FromOptions(ctx);
145
146 energy::BedThermalUnit::Ptr btu;
147
148 if (bedrock_grid.Mbz > 1) {
149 btu.reset(new energy::BTU_Verification(grid, bedrock_grid, 'K', false));
150 } else {
151 btu.reset(new energy::BTU_Minimal(grid));
152 }
153
154 InputOptions opts = process_input_options(com, config);
155 btu->init(opts);
156
157 double dt_seconds = units::convert(ctx->unit_system(), dt_years, "years", "seconds");
158
159 // worry about time step
160 int N = (int)ceil((ctx->time()->end() - ctx->time()->start()) / dt_seconds);
161 dt_seconds = (ctx->time()->end() - ctx->time()->start()) / (double)N;
162 log->message(2,
163 " user set timestep of %.4f years ...\n"
164 " reset to %.4f years to get integer number of steps ... \n",
165 dt_years.value(), units::convert(ctx->unit_system(), dt_seconds, "seconds", "years"));
166 MaxTimestep max_dt = btu->max_timestep(0.0);
167 log->message(2,
168 " BedThermalUnit reports max timestep of %.4f years ...\n",
169 units::convert(ctx->unit_system(), max_dt.value(), "seconds", "years"));
170
171 // actually do the time-stepping
172 log->message(2," running ...\n");
173 for (int n = 0; n < N; n++) {
174 // time at start of time-step
175 const double time = ctx->time()->start() + dt_seconds * (double)n;
176
177 // compute exact ice temperature at z=0 at time y
178 {
179 IceModelVec::AccessList list(bedtoptemp);
180 for (Points p(*grid); p; p.next()) {
181 const int i = p.i(), j = p.j();
182
183 bedtoptemp(i,j) = exactK(time, 0.0, 0).T;
184 }
185 }
186 // no need to update ghost values
187
188 // update the temperature inside the thermal layer using bedtoptemp
189 btu->update(bedtoptemp, time, dt_seconds);
190 log->message(2,".");
191 }
192
193 log->message(2, "\n done ...\n");
194
195 // compute final output heat flux G_0 at z=0
196 heat_flux_at_ice_base.copy_from(btu->flux_through_top_surface());
197
198 // get, and tell stdout, the correct answer from Test K
199 const double FF = exactK(ctx->time()->end(), 0.0, 0).F;
200 log->message(2,
201 " exact Test K reports upward heat flux at z=0, at end time %s, as G_0 = %.7f W m-2;\n",
202 ctx->time()->end_date().c_str(), FF);
203
204 // compute numerical error
205 heat_flux_at_ice_base.shift(-FF);
206 double max_error = heat_flux_at_ice_base.norm(NORM_INFINITY);
207 double avg_error = heat_flux_at_ice_base.norm(NORM_1);
208 heat_flux_at_ice_base.shift(+FF); // shift it back for writing
209 avg_error /= (grid->Mx() * grid->My());
210 log->message(2,
211 "case dt = %.5f:\n", dt_years.value());
212 log->message(1,
213 "NUMERICAL ERRORS in upward heat flux at z=0 relative to exact solution:\n");
214 log->message(1,
215 "bheatflx0 : max prcntmax av\n");
216 log->message(1,
217 " %11.7f %11.7f %11.7f\n",
218 max_error, 100.0*max_error/FF, avg_error);
219 log->message(1, "NUM ERRORS DONE\n");
220
221 File file(grid->com,
222 outname,
223 string_to_backend(config->get_string("output.format")),
224 PISM_READWRITE_MOVE,
225 ctx->pio_iosys_id());
226
227 io::define_time(file, *ctx);
228 io::append_time(file, *ctx->config(), ctx->time()->current());
229
230 btu->write_model_state(file);
231
232 bedtoptemp.write(file);
233 heat_flux_at_ice_base.write(file);
234
235 log->message(2, "done.\n");
236 }
237 catch (...) {
238 handle_fatal_errors(com);
239 return 1;
240 }
241
242 return 0;
243 }