tssa_testj.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch) (HTM) git clone git://src.adamsgaard.dk/pism (DIR) Log (DIR) Files (DIR) Refs (DIR) LICENSE --- tssa_testj.cc (6281B) --- 1 // Copyright (C) 2010--2018 Ed Bueler, Constantine Khroulev, and David Maxwell 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 "\nSSA_TESTJ\n" 21 " Testing program for the finite element implementation of the SSA.\n" 22 " Does a time-independent calculation. Does not run IceModel or a derived\n" 23 " class thereof. Uses verification test J. Also may be used in a PISM\n" 24 " software (regression) test.\n\n"; 25 26 #include "pism/basalstrength/basal_resistance.hh" // IceBasalResistancePlasticLaw 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/Mask.hh" 31 #include "pism/util/Context.hh" 32 #include "pism/util/VariableMetadata.hh" 33 #include "pism/util/error_handling.hh" 34 #include "pism/util/iceModelVec.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/exactTestsIJ.h" 40 41 namespace pism { 42 namespace stressbalance { 43 44 class SSATestCaseJ: public SSATestCase 45 { 46 public: 47 SSATestCaseJ(Context::Ptr ctx, int Mx, int My, SSAFactory ssafactory) 48 : SSATestCase(ctx, Mx, My, 300e3, 300e3, CELL_CENTER, XY_PERIODIC) { 49 m_config->set_flag("basal_resistance.pseudo_plastic.enabled", false); 50 51 m_enthalpyconverter = EnthalpyConverter::Ptr(new EnthalpyConverter(*m_config)); 52 m_config->set_string("stress_balance.ssa.flow_law", "isothermal_glen"); 53 54 m_ssa = ssafactory(m_grid); 55 } 56 57 protected: 58 virtual void initializeSSACoefficients(); 59 60 virtual void exactSolution(int i, int j, 61 double x, double y, double *u, double *v); 62 }; 63 64 void SSATestCaseJ::initializeSSACoefficients() { 65 m_tauc.set(0.0); // irrelevant for test J 66 m_geometry.bed_elevation.set(-1000.0); // assures shelf is floating (maximum ice thickness is 770 m) 67 m_geometry.cell_type.set(MASK_FLOATING); 68 69 double enth0 = m_enthalpyconverter->enthalpy(273.15, 0.01, 0.0); // 0.01 water fraction 70 m_ice_enthalpy.set(enth0); 71 72 /* use Ritz et al (2001) value of 30 MPa year for typical vertically-averaged viscosity */ 73 double ocean_rho = m_config->get_number("constants.sea_water.density"), 74 ice_rho = m_config->get_number("constants.ice.density"); 75 const double nu0 = units::convert(m_sys, 30.0, "MPa year", "Pa s"); /* = 9.45e14 Pa s */ 76 const double H0 = 500.; /* 500 m typical thickness */ 77 78 // Test J has a viscosity that is independent of velocity. So we force a 79 // constant viscosity by settting the strength_extension 80 // thickness larger than the given ice thickness. (max = 770m). 81 m_ssa->strength_extension->set_notional_strength(nu0 * H0); 82 m_ssa->strength_extension->set_min_thickness(800); 83 84 IceModelVec::AccessList list{&m_geometry.ice_thickness, &m_geometry.ice_surface_elevation, &m_bc_mask, &m_bc_values}; 85 86 for (Points p(*m_grid); p; p.next()) { 87 const int i = p.i(), j = p.j(); 88 89 const double myx = m_grid->x(i), myy = m_grid->y(j); 90 91 // set H,h on regular grid 92 struct TestJParameters J_parameters = exactJ(myx, myy); 93 94 m_geometry.ice_thickness(i,j) = J_parameters.H; 95 m_geometry.ice_surface_elevation(i,j) = (1.0 - ice_rho / ocean_rho) * J_parameters.H; // FIXME issue #15 96 97 // special case at center point: here we set bc_values at (i,j) by 98 // setting bc_mask and bc_values appropriately 99 if ((i == ((int)m_grid->Mx()) / 2) and 100 (j == ((int)m_grid->My()) / 2)) { 101 m_bc_mask(i,j) = 1; 102 m_bc_values(i,j).u = J_parameters.u; 103 m_bc_values(i,j).v = J_parameters.v; 104 } 105 } 106 107 // communicate what we have set 108 m_geometry.ice_surface_elevation.update_ghosts(); 109 m_geometry.ice_thickness.update_ghosts(); 110 m_bc_mask.update_ghosts(); 111 m_bc_values.update_ghosts(); 112 } 113 114 void SSATestCaseJ::exactSolution(int /*i*/, int /*j*/, 115 double x, double y, 116 double *u, double *v) { 117 struct TestJParameters J_parameters = exactJ(x, y); 118 *u = J_parameters.u; 119 *v = J_parameters.v; 120 } 121 122 } // end of namespace stressbalance 123 } // end of namespace pism 124 125 int main(int argc, char *argv[]) { 126 127 using namespace pism; 128 using namespace pism::stressbalance; 129 130 MPI_Comm com = MPI_COMM_WORLD; 131 petsc::Initializer petsc(argc, argv, help); 132 133 com = PETSC_COMM_WORLD; 134 135 /* This explicit scoping forces destructors to be called before PetscFinalize() */ 136 try { 137 Context::Ptr ctx = context_from_options(com, "ssa_testj"); 138 Config::Ptr config = ctx->config(); 139 140 std::string usage = "\n" 141 "usage of SSA_TESTJ:\n" 142 " run ssafe_test -Mx <number> -My <number> -ssa_method <fd|fem>\n" 143 "\n"; 144 145 bool stop = show_usage_check_req_opts(*ctx->log(), "ssa_testj", {}, usage); 146 147 if (stop) { 148 return 0; 149 } 150 151 // Parameters that can be overridden by command line options 152 153 unsigned int Mx = config->get_number("grid.Mx"); 154 unsigned int My = config->get_number("grid.My"); 155 156 auto method = config->get_string("stress_balance.ssa.method"); 157 auto output = config->get_string("output.file_name"); 158 159 // Determine the kind of solver to use. 160 SSAFactory ssafactory = NULL; 161 if (method == "fem") { 162 ssafactory = SSAFEMFactory; 163 } else if (method == "fd") { 164 ssafactory = SSAFDFactory; 165 } else { 166 /* can't happen */ 167 } 168 169 SSATestCaseJ testcase(ctx, Mx, My, ssafactory); 170 testcase.init(); 171 testcase.run(); 172 testcase.report("J"); 173 testcase.write(output); 174 } 175 catch (...) { 176 handle_fatal_errors(com); 177 return 1; 178 } 179 180 return 0; 181 }