22 " Testing program for the finite element implementation of the SSA.\n"
23 " Does a time-independent calculation. Does not run IceModel or a derived\n"
24 " class thereof. Uses verification test I. Also may be used in a PISM\n"
25 " software (regression) test.\n\n";
27#include "pism/stressbalance/ssa/SSAFD.hh"
28#include "pism/stressbalance/ssa/SSAFEM.hh"
29#include "pism/stressbalance/ssa/SSATestCase.hh"
30#include "pism/util/Context.hh"
31#include "pism/util/error_handling.hh"
32#include "pism/util/petscwrappers/PetscInitializer.hh"
33#include "pism/util/pism_options.hh"
34#include "pism/verification/tests/exactTestsIJ.h"
37namespace stressbalance {
47std::shared_ptr<Grid>
ssa_test_i_grid(std::shared_ptr<Context> ctx,
int Mx,
int My) {
49 std::max(60.0e3, ((Mx - 1) / 2) * (2.0 * (3.0 *
L_schoof) / (My - 1))),
65 double x,
double y,
double *u,
double *v);
76 for (
auto p =
m_grid->points(); p; p.next()) {
77 const int i = p.i(), j = p.j();
103 double *u,
double *v) {
112int main(
int argc,
char *argv[]) {
114 using namespace pism;
117 MPI_Comm com = MPI_COMM_WORLD;
120 com = PETSC_COMM_WORLD;
124 std::shared_ptr<Context> ctx = context_from_options(com,
"ssa_testi");
127 std::string usage =
"\n"
128 "usage of SSA_TESTi:\n"
129 " run ssa_testi -Mx <number> -My <number> -ssa_method <fd|fem>\n"
132 bool stop = show_usage_check_req_opts(*ctx->log(),
"ssa_testi", {}, usage);
139 unsigned int Mx = config->get_number(
"grid.Mx");
140 unsigned int My = config->get_number(
"grid.My");
142 auto method = config->get_string(
"stress_balance.ssa.method");
143 auto output_file = config->get_string(
"output.file");
145 bool write_output = config->get_string(
"output.size") !=
"none";
148 config->set_flag(
"basal_resistance.pseudo_plastic.enabled",
false);
150 config->set_string(
"stress_balance.ssa.flow_law",
"isothermal_glen");
151 config->set_number(
"flow_law.isothermal_Glen.ice_softness",
152 pow(B_schoof, -config->get_number(
"stress_balance.ssa.Glen_exponent")));
155 config->set_flag(
"stress_balance.ssa.compute_surface_gradient_inward",
true);
156 config->set_number(
"stress_balance.ssa.epsilon", 0.0);
158 auto grid = ssa_test_i_grid(ctx, Mx, My);
159 SSATestCaseI testcase(SSATestCase::solver(grid, method));
164 testcase.
write(output_file);
168 handle_fatal_errors(com);
std::shared_ptr< Config > Ptr
double enthalpy(double T, double omega, double P) const
Compute enthalpy from absolute temperature, liquid water fraction, and pressure.
Converts between specific enthalpy and temperature or liquid content.
void ensure_consistency(double ice_free_thickness_threshold)
array::Scalar2 ice_surface_elevation
array::Scalar2 ice_thickness
array::Scalar2 bed_elevation
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
void set(double c)
Result: v[j] <- c for all j.
void update_ghosts()
Updates ghost points.
virtual void exactSolution(int i, int j, double x, double y, double *u, double *v)
virtual void initializeSSACoefficients()
Set up the coefficient variables as appropriate for the test case.
SSATestCaseI(std::shared_ptr< SSA > ssa)
std::shared_ptr< const pism::Grid > m_grid
virtual void init()
Initialize the test case at the start of a run.
array::Array3D m_ice_enthalpy
static std::shared_ptr< Grid > grid(std::shared_ptr< Context > ctx, int Mx, int My, double Lx, double Ly, grid::Registration registration, grid::Periodicity periodicity)
array::Vector2 m_bc_values
virtual void report(const std::string &testname)
Report on the generated solution.
virtual void run()
Solve the SSA.
const Config::ConstPtr m_config
virtual void write(const std::string &filename)
Save the computation and data to a file.
struct TestIParameters exactI(const double m, const double x, const double y)
bool domain_edge(const Grid &grid, int i, int j)
const double aspect_schoof
std::shared_ptr< Grid > ssa_test_i_grid(std::shared_ptr< Context > ctx, int Mx, int My)
Stress balance models and related diagnostics.
int main(int argc, char *argv[])