tpismv.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
---
tpismv.cc (6556B)
---
1 // Copyright (C) 2004-2017, 2019 Jed Brown, 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 "Ice sheet driver for PISM (SIA and SSA) verification. Uses exact solutions\n"
21 " to various coupled subsystems. Computes difference between exact solution\n"
22 " and numerical solution.\n"
23 " Currently implements tests A, B, C, D, E, F, G, H, K, L.\n\n";
24
25 #include <string>
26
27 #include "pism/util/IceGrid.hh"
28 #include "pism/util/Config.hh"
29 #include "pism/util/error_handling.hh"
30 #include "pism/util/petscwrappers/PetscInitializer.hh"
31 #include "pism/util/pism_options.hh"
32 #include "pism/verification/iceCompModel.hh"
33 #include "pism/util/Context.hh"
34 #include "pism/util/Logger.hh"
35 #include "pism/util/Time.hh"
36 #include "pism/util/EnthalpyConverter.hh"
37
38 using namespace pism;
39
40 //! Allocate the PISMV (verification) context. Uses ColdEnthalpyConverter.
41 Context::Ptr pismv_context(MPI_Comm com, const std::string &prefix) {
42 // unit system
43 units::System::Ptr sys(new units::System);
44
45 // logger
46 Logger::Ptr logger = logger_from_options(com);
47
48 // configuration parameters
49 Config::Ptr config = config_from_options(com, *logger, sys);
50
51 config->set_string("time.calendar", "none");
52 config->set_string("grid.periodicity", "none");
53 config->set_string("grid.registration", "corner");
54
55 set_config_from_options(*config);
56
57 print_config(*logger, 3, *config);
58
59 Time::Ptr time = time_from_options(com, config, sys);
60
61 EnthalpyConverter::Ptr EC = EnthalpyConverter::Ptr(new ColdEnthalpyConverter(*config));
62
63 return Context::Ptr(new Context(com, sys, config, EC, time, logger, prefix));
64 }
65
66 GridParameters pismv_grid_defaults(Config::Ptr config,
67 char testname) {
68 // This sets the defaults for each test; command-line options can override this.
69
70 GridParameters P;
71
72 // use the cell corner grid registration
73 P.registration = CELL_CORNER;
74 // use the non-periodic grid:
75 P.periodicity = NOT_PERIODIC;
76 // equal spacing is the default for all the tests except K
77 P.Lx = config->get_number("grid.Lx");
78 P.Ly = config->get_number("grid.Ly");
79
80 P.Mx = config->get_number("grid.Mx");
81 P.My = config->get_number("grid.My");
82
83 SpacingType spacing = EQUAL;
84 double Lz = config->get_number("grid.Lz");
85 unsigned int Mz = config->get_number("grid.Mz");
86
87 switch (testname) {
88 case 'A':
89 case 'B':
90 case 'H':
91 // use 2400km by 2400km by 4000m rectangular domain
92 P.Lx = 1200e3;
93 P.Ly = P.Lx;
94 Lz = 4000;
95 break;
96 case 'C':
97 case 'D':
98 // use 2000km by 2000km by 4000m rectangular domain
99 P.Lx = 1000e3;
100 P.Ly = P.Lx;
101 Lz = 4000;
102 break;
103 case 'F':
104 case 'G':
105 case 'L':
106 // use 1800km by 1800km by 4000m rectangular domain
107 P.Lx = 900e3;
108 P.Ly = P.Lx;
109 Lz = 4000;
110 break;
111 case 'K':
112 case 'O':
113 // use 2000km by 2000km by 4000m rectangular domain, but make truely periodic
114 config->set_number("grid.Mbz", 2);
115 config->set_number("grid.Lbz", 1000);
116 P.Lx = 1000e3;
117 P.Ly = P.Lx;
118 Lz = 4000;
119 P.periodicity = XY_PERIODIC;
120 spacing = QUADRATIC;
121 break;
122 case 'V':
123 P.My = 3; // it's a flow-line setup
124 P.Lx = 500e3; // 500 km long
125 P.periodicity = Y_PERIODIC;
126 break;
127 default:
128 throw RuntimeError(PISM_ERROR_LOCATION, "desired test not implemented\n");
129 }
130
131 P.z = IceGrid::compute_vertical_levels(Lz, Mz, spacing,
132 config->get_number("grid.lambda"));
133 return P;
134 }
135
136 IceGrid::Ptr pismv_grid(Context::Ptr ctx, char testname) {
137 auto config = ctx->config();
138
139 auto input_file = config->get_string("input.file");
140
141 if (config->get_flag("input.bootstrap")) {
142 throw RuntimeError(PISM_ERROR_LOCATION, "pisms does not support bootstrapping");
143 }
144
145 if (not input_file.empty()) {
146 GridRegistration r = string_to_registration(ctx->config()->get_string("grid.registration"));
147
148 // get grid from a PISM input file
149 return IceGrid::FromFile(ctx, input_file, {"enthalpy", "temp"}, r);
150 } else {
151 // use defaults set by pismv_grid_defaults()
152 GridParameters P = pismv_grid_defaults(ctx->config(), testname);
153 P.horizontal_size_from_options();
154 P.horizontal_extent_from_options();
155 P.vertical_grid_from_options(ctx->config());
156 P.ownership_ranges_from_options(ctx->size());
157
158 return IceGrid::Ptr(new IceGrid(ctx, P));
159 }
160 }
161
162 int main(int argc, char *argv[]) {
163 MPI_Comm com = MPI_COMM_WORLD;
164
165 petsc::Initializer petsc(argc, argv, help);
166 com = PETSC_COMM_WORLD;
167
168 /* This explicit scoping forces destructors to be called before PetscFinalize() */
169 try {
170 Context::Ptr ctx = pismv_context(com, "pismv");
171 Logger::Ptr log = ctx->log();
172
173 std::string usage =
174 " pismv -test x [-no_report] [OTHER PISM & PETSc OPTIONS]\n"
175 "where:\n"
176 " -test x SIA-type verification test (x = A|B|C|D|F|G|H|K|L)\n"
177 " -no_report do not give error report at end of run\n"
178 "(see User's Manual for tests I and J).\n";
179
180 std::vector<std::string> required(1, "-test");
181
182 bool done = show_usage_check_req_opts(*log, "PISMV (verification mode)", required, usage);
183 if (done) {
184 return 0;
185 }
186
187 Config::Ptr config = ctx->config();
188
189 // determine test (and whether to report error)
190 std::string testname = options::Keyword("-test", "Specifies PISM verification test",
191 "A,B,C,D,F,G,H,K,L,V", "A");
192
193 IceGrid::Ptr g = pismv_grid(ctx, testname[0]);
194
195 IceCompModel m(g, ctx, testname[0]);
196
197 m.init();
198
199 m.run();
200 log->message(2, "done with run\n");
201
202 m.reportErrors();
203
204 // provide a default output file name if no -o option is given.
205 m.save_results();
206
207 print_unused_parameters(*log, 3, *config);
208 }
209 catch (...) {
210 handle_fatal_errors(com);
211 return 1;
212 }
213
214 return 0;
215 }
216