tAgeModel.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 --- tAgeModel.cc (8050B) --- 1 /* Copyright (C) 2016, 2017, 2019 PISM Authors 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 20 #include "AgeModel.hh" 21 22 #include "pism/age/AgeColumnSystem.hh" 23 #include "pism/util/error_handling.hh" 24 #include "pism/util/Vars.hh" 25 #include "pism/util/io/File.hh" 26 27 namespace pism { 28 29 AgeModelInputs::AgeModelInputs(const IceModelVec2S *thickness, 30 const IceModelVec3 *u, 31 const IceModelVec3 *v, 32 const IceModelVec3 *w) 33 : ice_thickness(thickness), u3(u), v3(v), w3(w) { 34 // empty 35 } 36 37 AgeModelInputs::AgeModelInputs() { 38 ice_thickness = NULL; 39 u3 = NULL; 40 v3 = NULL; 41 w3 = NULL; 42 } 43 44 static void check_input(const IceModelVec *ptr, const char *name) { 45 if (ptr == NULL) { 46 throw RuntimeError::formatted(PISM_ERROR_LOCATION, 47 "ice age model input %s was not provided", name); 48 } 49 } 50 51 void AgeModelInputs::check() const { 52 check_input(ice_thickness, "ice_thickness"); 53 check_input(u3, "u3"); 54 check_input(v3, "v3"); 55 check_input(w3, "w3"); 56 } 57 58 AgeModel::AgeModel(IceGrid::ConstPtr grid, stressbalance::StressBalance *stress_balance) 59 : Component(grid), 60 // FIXME: should be able to use width=1... 61 m_ice_age(m_grid, "age", WITH_GHOSTS, m_config->get_number("grid.max_stencil_width")), 62 m_work(m_grid, "work_vector", WITHOUT_GHOSTS), 63 m_stress_balance(stress_balance) { 64 65 m_ice_age.set_attrs("model_state", "age of ice", 66 "s", "years", "" /* no standard name*/, 0); 67 68 m_ice_age.metadata().set_number("valid_min", 0.0); 69 70 m_work.set_attrs("internal", "new values of age during time step", 71 "s", "s", "", 0); 72 } 73 74 /*! 75 Let \f$\tau(t,x,y,z)\f$ be the age of the ice. Denote the three-dimensional 76 velocity field within the ice fluid as \f$(u,v,w)\f$. The age equation 77 is \f$d\tau/dt = 1\f$, that is, ice may move but it gets one year older in one 78 year. Thus 79 \f[ \frac{\partial \tau}{\partial t} + u \frac{\partial \tau}{\partial x} 80 + v \frac{\partial \tau}{\partial y} + w \frac{\partial \tau}{\partial z} = 1 \f] 81 This equation is purely advective and hyperbolic. The right-hand side is "1" as 82 long as age \f$\tau\f$ and time \f$t\f$ are measured in the same units. 83 Because the velocity field is incompressible, \f$\nabla \cdot (u,v,w) = 0\f$, 84 we can rewrite the equation as 85 \f[ \frac{\partial \tau}{\partial t} + \nabla \left( (u,v,w) \tau \right) = 1 \f] 86 There is a conservative first-order numerical method; see AgeColumnSystem::solveThisColumn(). 87 88 The boundary condition is that when the ice falls as snow it has age zero. 89 That is, \f$\tau(t,x,y,h(t,x,y)) = 0\f$ in accumulation areas. There is no 90 boundary condition elsewhere on the ice upper surface, as the characteristics 91 go outward in the ablation zone. If the velocity in the bottom cell of ice 92 is upward (\f$w>0\f$) then we also apply a zero age boundary condition, 93 \f$\tau(t,x,y,0) = 0\f$. This is the case where ice freezes on at the base, 94 either grounded basal ice freezing on stored water in till, or marine basal ice. 95 (Note that the water that is frozen-on as ice might be quite "old" in the sense 96 that its most recent time in the atmosphere was long ago; this comment is 97 relevant to any analysis which relates isotope ratios to modeled age.) 98 99 The numerical method is a conservative form of first-order upwinding, but the 100 vertical advection term is computed implicitly. Thus there is no CFL-type 101 stability condition from the vertical velocity; CFL is only for the horizontal 102 velocity. We use a finely-spaced, equally-spaced vertical grid in the 103 calculation. Note that the columnSystemCtx methods coarse_to_fine() and 104 fine_to_coarse() interpolate back and forth between this fine grid and 105 the storage grid. The storage grid may or may not be equally-spaced. See 106 AgeColumnSystem::solve() for the actual method. 107 */ 108 void AgeModel::update(double t, double dt, const AgeModelInputs &inputs) { 109 110 // fix a compiler warning 111 (void) t; 112 113 inputs.check(); 114 115 const IceModelVec2S &ice_thickness = *inputs.ice_thickness; 116 117 const IceModelVec3 118 &u3 = *inputs.u3, 119 &v3 = *inputs.v3, 120 &w3 = *inputs.w3; 121 122 AgeColumnSystem system(m_grid->z(), "age", 123 m_grid->dx(), m_grid->dy(), dt, 124 m_ice_age, u3, v3, w3); // linear system to solve in each column 125 126 size_t Mz_fine = system.z().size(); 127 std::vector<double> x(Mz_fine); // space for solution 128 129 IceModelVec::AccessList list{&ice_thickness, &u3, &v3, &w3, &m_ice_age, &m_work}; 130 131 unsigned int Mz = m_grid->Mz(); 132 133 ParallelSection loop(m_grid->com); 134 try { 135 for (Points p(*m_grid); p; p.next()) { 136 const int i = p.i(), j = p.j(); 137 138 system.init(i, j, ice_thickness(i, j)); 139 140 if (system.ks() == 0) { 141 // if no ice, set the entire column to zero age 142 m_work.set_column(i, j, 0.0); 143 } else { 144 // general case: solve advection PDE 145 146 // solve the system for this column; call checks that params set 147 system.solve(x); 148 149 // put solution in IceModelVec3 150 system.fine_to_coarse(x, i, j, m_work); 151 152 // Ensure that the age of the ice is non-negative. 153 // 154 // FIXME: this is a kludge. We need to ensure that our numerical method has the maximum 155 // principle instead. (We may still need this for correctness, though.) 156 double *column = m_work.get_column(i, j); 157 for (unsigned int k = 0; k < Mz; ++k) { 158 if (column[k] < 0.0) { 159 column[k] = 0.0; 160 } 161 } 162 } 163 } 164 } catch (...) { 165 loop.failed(); 166 } 167 loop.check(); 168 169 m_work.update_ghosts(m_ice_age); 170 } 171 172 const IceModelVec3 & AgeModel::age() const { 173 return m_ice_age; 174 } 175 176 MaxTimestep AgeModel::max_timestep_impl(double t) const { 177 // fix a compiler warning 178 (void) t; 179 180 if (m_stress_balance == NULL) { 181 throw RuntimeError::formatted(PISM_ERROR_LOCATION, 182 "AgeModel: no stress balance provided." 183 " Cannot compute max. time step."); 184 } 185 186 return MaxTimestep(m_stress_balance->max_timestep_cfl_3d().dt_max.value(), "age model"); 187 } 188 189 void AgeModel::init(const InputOptions &opts) { 190 191 m_log->message(2, "* Initializing the age model...\n"); 192 193 194 double initial_age_years = m_config->get_number("age.initial_value", "years"); 195 196 if (opts.type == INIT_RESTART) { 197 File input_file(m_grid->com, opts.filename, PISM_GUESS, PISM_READONLY); 198 199 if (input_file.find_variable("age")) { 200 m_ice_age.read(input_file, opts.record); 201 } else { 202 m_log->message(2, 203 "PISM WARNING: input file '%s' does not have the 'age' variable.\n" 204 " Setting it to %f years...\n", 205 opts.filename.c_str(), initial_age_years); 206 m_ice_age.set(m_config->get_number("age.initial_value", "seconds")); 207 } 208 } else { 209 m_log->message(2, " - setting initial age to %.4f years\n", initial_age_years); 210 m_ice_age.set(m_config->get_number("age.initial_value", "seconds")); 211 } 212 213 regrid("Age Model", m_ice_age, REGRID_WITHOUT_REGRID_VARS); 214 } 215 216 void AgeModel::define_model_state_impl(const File &output) const { 217 m_ice_age.define(output); 218 } 219 220 void AgeModel::write_model_state_impl(const File &output) const { 221 m_ice_age.write(output); 222 } 223 224 } // end of namespace pism