21 " Testing program for SIA, time-independent calculations separate from\n"
22 " IceModel. Uses verification test F. Also may be used in a PISM software"
23 "(regression) test.\n\n";
25#include "pism/stressbalance/sia/SIAFD.hh"
26#include "pism/util/EnthalpyConverter.hh"
27#include "pism/stressbalance/StressBalance.hh"
28#include "pism/stressbalance/SSB_Modifier.hh"
29#include "pism/stressbalance/ShallowStressBalance.hh"
30#include "pism/util/Grid.hh"
31#include "pism/util/Mask.hh"
32#include "pism/util/Context.hh"
33#include "pism/util/Time.hh"
34#include "pism/util/error_handling.hh"
35#include "pism/util/io/File.hh"
36#include "pism/util/petscwrappers/PetscInitializer.hh"
37#include "pism/util/pism_utilities.hh"
38#include "pism/util/pism_options.hh"
39#include "pism/verification/tests/exactTestsFG.hh"
40#include "pism/util/io/io_helpers.hh"
41#include "pism/geometry/Geometry.hh"
48 double &gmax_strain_heating_err,
49 double &gav_strain_heating_err) {
50 double max_strain_heating_error = 0.0, av_strain_heating_error = 0.0, avcount = 0.0;
51 const int Mz = grid.Mz();
53 const double LforFG = 750000;
56 ice_rho = grid.ctx()->config()->get_number(
"constants.ice.density"),
57 ice_c = grid.ctx()->config()->get_number(
"constants.ice.specific_heat_capacity");
63 for (
auto p = grid.points(); p; p.next()) {
64 const int i = p.i(), j = p.j();
69 r = sqrt(pow(xx, 2) + pow(yy, 2));
71 if ((r >= 1.0) && (r <= LforFG - 1.0)) {
76 for (
int k = 0;
k < Mz;
k++) {
77 F.Sig[
k] *= ice_rho * ice_c;
79 const int ks = grid.kBelowHeight(thickness(i, j));
80 const double *strain_heating_ij = strain_heating.
get_column(i, j);
81 for (
int k = 0;
k < ks;
k++) {
82 const double _strain_heating_error = fabs(strain_heating_ij[
k] -
F.Sig[
k]);
83 max_strain_heating_error = std::max(max_strain_heating_error, _strain_heating_error);
85 av_strain_heating_error += _strain_heating_error;
94 gmax_strain_heating_err =
GlobalMax(grid.com, max_strain_heating_error);
95 gav_strain_heating_err =
GlobalSum(grid.com, av_strain_heating_error);
96 double gavcount =
GlobalSum(grid.com, avcount);
97 gav_strain_heating_err = gav_strain_heating_err/std::max(gavcount,1.0);
110 double maxUerr = 0.0, maxWerr = 0.0, avUerr = 0.0, avWerr = 0.0;
112 const double LforFG = 750000;
116 for (
auto p = grid.points(); p; p.next()) {
117 const int i = p.i(), j = p.j();
119 double xx = grid.x(i), yy = grid.y(j),
120 r = sqrt(pow(xx, 2) + pow(yy, 2));
121 if ((r >= 1.0) && (r <= LforFG - 1.0)) {
123 const double H = ice_thickness(i, j);
124 std::vector<double> z(1, H);
127 const double uex = (xx/r) *
F.U[0];
128 const double vex = (yy/r) *
F.U[0];
131 const double Uerr = sqrt(pow(u3.
interpolate(i,j,H) - uex, 2) +
133 maxUerr = std::max(maxUerr,Uerr);
135 const double Werr = fabs(w3.
interpolate(i,j,H) -
F.w[0]);
136 maxWerr = std::max(maxWerr,Werr);
144 gavUerr = gavUerr/(grid.Mx()*grid.My());
146 gavWerr = gavWerr/(grid.Mx()*grid.My());
158 for (
auto p = grid.points(); p; p.next()) {
159 const int i = p.i(), j = p.j();
161 const double *T_ij = temperature.
get_column(i,j);
164 for (
unsigned int k=0;
k<grid.Mz(); ++
k) {
165 double depth = thickness(i,j) - grid.z(
k);
195 for (
auto p = grid.points(); p; p.next()) {
196 const int i = p.i(), j = p.j();
202 if (r > LforFG - 1.0) {
203 thickness(i, j) = 0.0;
208 thickness(i, j) =
F.H;
233 double max_strain_heating_error, av_strain_heating_error;
236 max_strain_heating_error, av_strain_heating_error);
239 "Sigma : maxSig avSig\n");
240 log->message(1,
" %12.6f%12.6f\n",
241 max_strain_heating_error*1.0e6, av_strain_heating_error*1.0e6);
244 double maxUerr, avUerr, maxWerr, avWerr;
253 "surf vels : maxUvec avUvec maxW avW\n");
254 log->message(1,
" %12.6f%12.6f%12.6f%12.6f\n",
263int main(
int argc,
char *argv[]) {
265 using namespace pism;
268 MPI_Comm com = MPI_COMM_WORLD;
271 com = MPI_COMM_WORLD;
276 std::shared_ptr<Context> ctx = context_from_options(com,
"siafd_test");
279 config->set_flag(
"stress_balance.sia.grain_size_age_coupling",
false);
280 config->set_string(
"stress_balance.sia.flow_law",
"arr");
282 set_config_from_options(ctx->unit_system(), *config);
283 config->resolve_filenames();
285 std::string usage =
"\n"
286 "usage of SIAFD_TEST:\n"
287 " run siafd_test -Mx <number> -My <number> -Mz <number> -o foo.nc\n"
290 bool stop = show_usage_check_req_opts(*ctx->log(),
"siafd_test", {}, usage);
296 auto output_file = config->get_string(
"output.file");
303 unsigned int Mz = config->get_number(
"grid.Mz");
305 P.
z = grid::compute_vertical_levels(Lz, Mz, grid::EQUAL);
309 auto grid = std::make_shared<Grid>(ctx, P);
310 grid->report_parameters();
314 const int WIDE_STENCIL = config->get_number(
"grid.max_stencil_width");
317 enthalpy(grid,
"enthalpy", array::WITH_GHOSTS, grid->z(), WIDE_STENCIL),
318 age(grid,
"age", array::WITHOUT_GHOSTS, grid->z());
329 .long_name(
"ice enthalpy (includes sensible heat, latent heat, pressure)")
332 enthalpy.set(EC->enthalpy(263.15, 0.0, EC->pressure(1000.0)));
339 std::shared_ptr<SIAFD> sia(
new SIAFD(grid));
340 std::shared_ptr<ZeroSliding> no_sliding(
new ZeroSliding(grid));
346 setInitStateF(*grid, *EC,
353 geometry.
ensure_consistency(config->get_number(
"geometry.ice_free_thickness_standard"));
356 stress_balance.
init();
358 bool full_update =
true;
366 stress_balance.
update(inputs, full_update);
376 reportErrors(*grid, ctx->unit_system(),
380 File file(grid->com, output_file, io::PISM_NETCDF3, io::PISM_READWRITE_MOVE);
381 io::define_time(file, *ctx);
382 io::append_time(file, *ctx->config(), ctx->time()->current());
389 sia->diffusivity().write(file);
396 handle_fatal_errors(com);
An EnthalpyConverter for use in verification tests.
std::shared_ptr< Config > Ptr
std::shared_ptr< EnthalpyConverter > Ptr
double enthalpy_permissive(double T, double omega, double P) const
Compute enthalpy more permissively than enthalpy().
double pressure(double depth) const
Get pressure in ice from depth below surface using the hydrostatic assumption.
Converts between specific enthalpy and temperature or liquid content.
High-level PISM I/O class.
array::Scalar1 sea_level_elevation
void ensure_consistency(double ice_free_thickness_threshold)
array::Scalar2 ice_surface_elevation
array::CellType2 cell_type
array::Scalar2 ice_thickness
array::Scalar2 bed_elevation
Describes the PISM grid and the distribution of data across processors.
std::shared_ptr< const Logger > ConstPtr
void failed()
Indicates a failure of a parallel section.
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
void copy_from(const Array2D< T > &source)
double interpolate(int i, int j, double z) const
Return value of scalar quantity at level z (m) above base of ice (by linear interpolation).
void set_column(int i, int j, double c)
Set all values of scalar quantity to given a single value in a particular column.
double * get_column(int i, int j)
A virtual class collecting methods common to ice and bedrock 3D fields.
void write(const std::string &filename) const
void set(double c)
Result: v[j] <- c for all j.
void update_ghosts()
Updates ghost points.
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
std::vector< double > z
Vertical levels.
void ownership_ranges_from_options(const Config &config, unsigned int size)
Re-compute ownership ranges. Uses current values of Mx and My.
double Ly
Domain half-width in the Y direction.
double Lx
Domain half-width in the X direction.
Grid parameters; used to collect defaults before an Grid is allocated.
const array::Array3D & velocity_w() const
const array::Array3D & velocity_u() const
Get components of the the 3D velocity field.
void init()
Initialize the StressBalance object.
void update(const Inputs &inputs, bool full_update)
Update all the fields if (full_update), only update diffusive flux and max. diffusivity otherwise.
const array::Array3D & velocity_v() const
const array::Array3D & volumetric_strain_heating() const
The class defining PISM's interface to the shallow stress balance code.
Returns zero velocity field, zero friction heating, and zero for D^2.
std::shared_ptr< System > Ptr
double radius(const Grid &grid, int i, int j)
Returns the distance from the point (i,j) to the origin.
Stress balance models and related diagnostics.
double convert(System::Ptr system, double input, const std::string &spec1, const std::string &spec2)
Convert a quantity from unit1 to unit2.
static double F(double SL, double B, double H, double alpha)
static void compute_strain_heating_errors(const array::Array3D &strain_heating, const array::Scalar &thickness, const Grid &grid, double &gmax_strain_heating_err, double &gav_strain_heating_err)
static void computeSurfaceVelocityErrors(const Grid &grid, const array::Scalar &ice_thickness, const array::Array3D &u3, const array::Array3D &v3, const array::Array3D &w3, double &gmaxUerr, double &gavUerr, double &gmaxWerr, double &gavWerr)
static void enthalpy_from_temperature_cold(EnthalpyConverter &EC, const Grid &grid, const array::Scalar &thickness, const array::Array3D &temperature, array::Array3D &enthalpy)
void GlobalMax(MPI_Comm comm, double *local, double *result, int count)
TestFGParameters exactFG(double t, double r, const std::vector< double > &z, double Cp)
static void reportErrors(const Grid &grid, units::System::Ptr unit_system, const array::Scalar &thickness, const array::Array3D &u_sia, const array::Array3D &v_sia, const array::Array3D &w_sia, const array::Array3D &strain_heating)
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
static void setInitStateF(Grid &grid, EnthalpyConverter &EC, array::Scalar &bed, array::Scalar &mask, array::Scalar &surface, array::Scalar &thickness, array::Array3D &enthalpy)
Set the test F initial state.
int main(int argc, char *argv[])