PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
btutest.cc
Go to the documentation of this file.
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023 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 #include <petsc.h>
20 
21 static char help[] =
22  "Tests BedThermalUnit using Test K, without IceModel.\n\n";
23 
24 #include "pism/util/pism_options.hh"
25 #include "pism/util/Grid.hh"
26 #include "pism/util/io/File.hh"
27 #include "pism/verification/BTU_Verification.hh"
28 #include "pism/energy/BTU_Minimal.hh"
29 #include "pism/util/Time.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/EnthalpyConverter.hh"
39 #include "pism/util/MaxTimestep.hh"
40 #include "pism/util/Logger.hh"
41 
42 //! Allocate the PISMV (verification) context. Uses ColdEnthalpyConverter.
43 std::shared_ptr<pism::Context> btutest_context(MPI_Comm com, const std::string &prefix) {
44  using namespace pism;
45 
46  // unit system
48 
49  // logger
50  Logger::Ptr logger = logger_from_options(com);
51 
52  // configuration parameters
53  Config::Ptr config = config_from_options(com, *logger, sys);
54 
55  // default vertical grid parameters
56  config->set_number("grid.Mbz", 11);
57  config->set_number("grid.Lbz", 1000);
58 
59  // when Grid constructor is called, these settings are used
60  config->set_string("time.start", "0s");
61  config->set_number("time.run_length", 1.0);
62 
63  set_config_from_options(sys, *config);
64  config->resolve_filenames();
65 
66  print_config(*logger, 3, *config);
67 
68  Time::Ptr time = std::make_shared<Time>(com, config, *logger, sys);
69 
71 
72  return std::shared_ptr<Context>(new Context(com, sys, config, EC, time, logger, prefix));
73 }
74 
75 int main(int argc, char *argv[]) {
76 
77  using namespace pism;
78 
79  MPI_Comm com = MPI_COMM_WORLD;
80  petsc::Initializer petsc(argc, argv, help);
81 
82  com = MPI_COMM_WORLD;
83 
84  try {
85  std::shared_ptr<Context> ctx = btutest_context(com, "btutest");
86  Logger::Ptr log = ctx->log();
87 
88  std::string usage =
89  " btutest -Mbz NN -Lbz 1000.0 [-o OUT.nc -ys A -ye B -dt C -Mz D -Lz E]\n"
90  "where these are required because they are used in BedThermalUnit:\n"
91  " -Mbz number of bedrock thermal layer levels to use\n"
92  " -Lbz 1000.0 depth of bedrock thermal layer (required; Lbz=1000.0 m in Test K)\n"
93  "and these are allowed:\n"
94  " -o output file name; NetCDF format\n"
95  " -ys start year in using Test K\n"
96  " -ye end year in using Test K\n"
97  " -dt time step B (= positive float) in years\n";
98 
99  bool done = show_usage_check_req_opts(*log, "BTUTEST %s (test program for BedThermalUnit)",
100  {"-Mbz"}, usage);
101  if (done) {
102  return 0;
103  }
104 
105  log->message(2,
106  " initializing Grid from options ...\n");
107 
108  Config::Ptr config = ctx->config();
109 
110  grid::Parameters P(*config);
111  P.Mx = 3;
112  P.My = P.Mx;
113  P.Lx = 1500e3;
114  P.Ly = P.Lx;
115 
116  P.vertical_grid_from_options(config);
117  P.ownership_ranges_from_options(ctx->size());
118 
119  // create grid and set defaults
120  std::shared_ptr<Grid> grid(new Grid(ctx, P));
121 
122  auto outname = config->get_string("output.file");
123 
124  options::Real dt_years(ctx->unit_system(),
125  "-dt", "Time-step, in years", "years", 1.0);
126 
127  // allocate tools and Arrays
128  array::Scalar bedtoptemp(grid, "bedtoptemp");
129  bedtoptemp.metadata(0).long_name("temperature at top of bedrock thermal layer").units("K");
130 
131  array::Scalar heat_flux_at_ice_base(grid, "upward_heat_flux_at_ice_base");
132  heat_flux_at_ice_base.metadata(0)
133  .long_name("upward geothermal flux at bedrock thermal layer base")
134  .units("W m-2")
135  .output_units("mW m-2");
136 
137  // initialize BTU object:
138  energy::BTUGrid bedrock_grid = energy::BTUGrid::FromOptions(ctx);
139 
141 
142  if (bedrock_grid.Mbz > 1) {
143  btu.reset(new energy::BTU_Verification(grid, bedrock_grid, 'K', false));
144  } else {
145  btu.reset(new energy::BTU_Minimal(grid));
146  }
147 
148  InputOptions opts = process_input_options(com, config);
149  btu->init(opts);
150 
151  double dt_seconds = units::convert(ctx->unit_system(), dt_years, "years", "seconds");
152 
153  // worry about time step
154  int N = (int)ceil((ctx->time()->end() - ctx->time()->start()) / dt_seconds);
155  dt_seconds = (ctx->time()->end() - ctx->time()->start()) / (double)N;
156  log->message(2,
157  " user set timestep of %.4f years ...\n"
158  " reset to %.4f years to get integer number of steps ... \n",
159  dt_years.value(), units::convert(ctx->unit_system(), dt_seconds, "seconds", "years"));
160  MaxTimestep max_dt = btu->max_timestep(0.0);
161  log->message(2,
162  " BedThermalUnit reports max timestep of %.4f years ...\n",
163  units::convert(ctx->unit_system(), max_dt.value(), "seconds", "years"));
164 
165  // actually do the time-stepping
166  log->message(2," running ...\n");
167  for (int n = 0; n < N; n++) {
168  // time at start of time-step
169  const double time = ctx->time()->start() + dt_seconds * (double)n;
170 
171  // compute exact ice temperature at z=0 at time y
172  {
173  array::AccessScope list(bedtoptemp);
174  for (auto p = grid->points(); p; p.next()) {
175  const int i = p.i(), j = p.j();
176 
177  bedtoptemp(i,j) = exactK(time, 0.0, 0).T;
178  }
179  }
180  // no need to update ghost values
181 
182  // update the temperature inside the thermal layer using bedtoptemp
183  btu->update(bedtoptemp, time, dt_seconds);
184  log->message(2,".");
185  }
186 
187  log->message(2, "\n done ...\n");
188 
189  // compute final output heat flux G_0 at z=0
190  heat_flux_at_ice_base.copy_from(btu->flux_through_top_surface());
191 
192  auto time = ctx->time();
193 
194  // get, and tell stdout, the correct answer from Test K
195  const double FF = exactK(time->end(), 0.0, 0).F;
196  log->message(2,
197  " exact Test K reports upward heat flux at z=0, at end time %s, as G_0 = %.7f W m-2;\n",
198  time->date(time->end()).c_str(), FF);
199 
200  // compute numerical error
201  heat_flux_at_ice_base.shift(-FF);
202  double max_error = heat_flux_at_ice_base.norm(NORM_INFINITY)[0];
203  double avg_error = heat_flux_at_ice_base.norm(NORM_1)[0];
204  heat_flux_at_ice_base.shift(+FF); // shift it back for writing
205  avg_error /= (grid->Mx() * grid->My());
206  log->message(2,
207  "case dt = %.5f:\n", dt_years.value());
208  log->message(1,
209  "NUMERICAL ERRORS in upward heat flux at z=0 relative to exact solution:\n");
210  log->message(1,
211  "bheatflx0 : max prcntmax av\n");
212  log->message(1,
213  " %11.7f %11.7f %11.7f\n",
214  max_error, 100.0*max_error/FF, avg_error);
215  log->message(1, "NUM ERRORS DONE\n");
216 
217  File file(grid->com,
218  outname,
219  string_to_backend(config->get_string("output.format")),
221  ctx->pio_iosys_id());
222 
223  io::define_time(file, *ctx);
224  io::append_time(file, *ctx->config(), ctx->time()->current());
225 
226  btu->write_model_state(file);
227 
228  bedtoptemp.write(file);
229  heat_flux_at_ice_base.write(file);
230 
231  log->message(2, "done.\n");
232  }
233  catch (...) {
234  handle_fatal_errors(com);
235  return 1;
236  }
237 
238  return 0;
239 }
int main(int argc, char *argv[])
Definition: btutest.cc:75
static char help[]
Definition: btutest.cc:21
std::shared_ptr< pism::Context > btutest_context(MPI_Comm com, const std::string &prefix)
Allocate the PISMV (verification) context. Uses ColdEnthalpyConverter.
Definition: btutest.cc:43
An EnthalpyConverter for use in verification tests.
std::shared_ptr< Config > Ptr
std::shared_ptr< EnthalpyConverter > Ptr
High-level PISM I/O class.
Definition: File.hh:56
Describes the PISM grid and the distribution of data across processors.
Definition: Grid.hh:282
std::shared_ptr< Logger > Ptr
Definition: Logger.hh:45
double value() const
Get the value of the maximum time step.
Definition: MaxTimestep.cc:55
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
Definition: MaxTimestep.hh:31
std::shared_ptr< Time > Ptr
Definition: Time.hh:62
VariableMetadata & long_name(const std::string &input)
VariableMetadata & output_units(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:65
void copy_from(const Array2D< T > &source)
Definition: Array2D.hh:73
void shift(double alpha)
Result: v[j] <- v[j] + alpha for all j. Calls VecShift.
Definition: Array.cc:245
void write(const std::string &filename) const
Definition: Array.cc:800
std::vector< double > norm(int n) const
Computes the norm of all the components of an Array.
Definition: Array.cc:746
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: Array.cc:553
std::shared_ptr< BedThermalUnit > Ptr
void vertical_grid_from_options(std::shared_ptr< const Config > config)
Process -Mz and -Lz; set z;.
Definition: Grid.cc:1209
void ownership_ranges_from_options(unsigned int size)
Re-compute ownership ranges. Uses current values of Mx and My.
Definition: Grid.cc:1142
unsigned int Mx
Number of grid points in the X direction.
Definition: Grid.hh:150
double Ly
Domain half-width in the Y direction.
Definition: Grid.hh:144
double Lx
Domain half-width in the X direction.
Definition: Grid.hh:142
unsigned int My
Number of grid points in the Y direction.
Definition: Grid.hh:152
Grid parameters; used to collect defaults before an Grid is allocated.
Definition: Grid.hh:117
std::shared_ptr< System > Ptr
Definition: Units.hh:47
struct TestKParameters exactK(double t, double z, int bedrock_is_ice)
Definition: exactTestK.c:129
#define n
Definition: exactTestM.c:37
void append_time(const File &file, const Config &config, double time_seconds)
Prepare a file for output.
Definition: io_helpers.cc:247
@ PISM_READWRITE_MOVE
create a file for writing, move foo.nc to foo.nc~ if present
Definition: IO_Flags.hh:78
void define_time(const File &file, const Context &ctx)
Prepare a file for output.
Definition: io_helpers.cc:213
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
io::Backend string_to_backend(const std::string &backend)
Definition: File.cc:57
InputOptions process_input_options(MPI_Comm com, Config::ConstPtr config)
Definition: Component.cc:43
Logger::Ptr logger_from_options(MPI_Comm com)
Definition: Logger.cc:107
Config::Ptr config_from_options(MPI_Comm com, const Logger &log, units::System::Ptr unit_system)
Create a configuration database using command-line options.
void handle_fatal_errors(MPI_Comm com)
void set_config_from_options(units::System::Ptr unit_system, Config &config)
Set configuration parameters using command-line options.
bool show_usage_check_req_opts(const Logger &log, const std::string &execname, const std::vector< std::string > &required_options, const std::string &usage)
In a single call a driver program can provide a usage string to the user and check if required option...
Definition: pism_options.cc:51
void print_config(const Logger &log, int verbosity_threshhold, const Config &config)
Report configuration parameters to stdout.
static BTUGrid FromOptions(std::shared_ptr< const Context > ctx)