pism

[fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
git clone git://src.adamsgaard.dk/pism # fast
git clone https://src.adamsgaard.dk/pism.git # slow
Log | Files | Refs | README | LICENSE Back to index

EnergyModel.cc (13772B)


      1 /* Copyright (C) 2016, 2017, 2018, 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 "EnergyModel.hh"
     21 #include "pism/util/MaxTimestep.hh"
     22 #include "pism/stressbalance/StressBalance.hh"
     23 #include "pism/util/io/File.hh"
     24 #include "pism/util/Vars.hh"
     25 #include "utilities.hh"
     26 #include "pism/util/EnthalpyConverter.hh"
     27 #include "bootstrapping.hh"
     28 #include "pism/util/pism_utilities.hh"
     29 #include "pism/util/error_handling.hh"
     30 #include "pism/util/IceModelVec2CellType.hh"
     31 #include "pism/util/pism_options.hh"
     32 #include "pism/util/Profiling.hh"
     33 
     34 namespace pism {
     35 namespace energy {
     36 
     37 static void check_input(const IceModelVec *ptr, const char *name) {
     38   if (ptr == NULL) {
     39     throw RuntimeError::formatted(PISM_ERROR_LOCATION, "energy balance model input %s was not provided", name);
     40   }
     41 }
     42 
     43 Inputs::Inputs() {
     44   basal_frictional_heating = NULL;
     45   basal_heat_flux          = NULL;
     46   cell_type                = NULL;
     47   ice_thickness            = NULL;
     48   surface_liquid_fraction  = NULL;
     49   shelf_base_temp          = NULL;
     50   surface_temp             = NULL;
     51   till_water_thickness     = NULL;
     52 
     53   volumetric_heating_rate  = NULL;
     54   u3                       = NULL;
     55   v3                       = NULL;
     56   w3                       = NULL;
     57 
     58   no_model_mask = NULL;
     59 }
     60 
     61 void Inputs::check() const {
     62   check_input(cell_type,                "cell_type");
     63   check_input(basal_frictional_heating, "basal_frictional_heating");
     64   check_input(basal_heat_flux,          "basal_heat_flux");
     65   check_input(ice_thickness,            "ice_thickness");
     66   check_input(surface_liquid_fraction,  "surface_liquid_fraction");
     67   check_input(shelf_base_temp,          "shelf_base_temp");
     68   check_input(surface_temp,             "surface_temp");
     69   check_input(till_water_thickness,     "till_water_thickness");
     70 
     71   check_input(volumetric_heating_rate, "volumetric_heating_rate");
     72   check_input(u3, "u3");
     73   check_input(v3, "v3");
     74   check_input(w3, "w3");
     75 }
     76 
     77 EnergyModelStats::EnergyModelStats() {
     78   bulge_counter            = 0;
     79   reduced_accuracy_counter = 0;
     80   low_temperature_counter  = 0;
     81   liquified_ice_volume     = 0.0;
     82 }
     83 
     84 EnergyModelStats& EnergyModelStats::operator+=(const EnergyModelStats &other) {
     85   bulge_counter            += other.bulge_counter;
     86   reduced_accuracy_counter += other.reduced_accuracy_counter;
     87   low_temperature_counter  += other.low_temperature_counter;
     88   liquified_ice_volume     += other.liquified_ice_volume;
     89   return *this;
     90 }
     91 
     92 
     93 bool marginal(const IceModelVec2S &thickness, int i, int j, double threshold) {
     94   int
     95     n = j + 1,
     96     e = i + 1,
     97     s = j - 1,
     98     w = i - 1;
     99 
    100   const double
    101     N  = thickness(i, n),
    102     E  = thickness(e, j),
    103     S  = thickness(i, s),
    104     W  = thickness(w, j),
    105     NW = thickness(w, n),
    106     SW = thickness(w, s),
    107     NE = thickness(e, n),
    108     SE = thickness(e, s);
    109 
    110   return ((E  < threshold) or
    111           (NE < threshold) or
    112           (N  < threshold) or
    113           (NW < threshold) or
    114           (W  < threshold) or
    115           (SW < threshold) or
    116           (S  < threshold) or
    117           (SE < threshold));
    118 }
    119 
    120 
    121 void EnergyModelStats::sum(MPI_Comm com) {
    122   bulge_counter            = GlobalSum(com, bulge_counter);
    123   reduced_accuracy_counter = GlobalSum(com, reduced_accuracy_counter);
    124   low_temperature_counter  = GlobalSum(com, low_temperature_counter);
    125   liquified_ice_volume     = GlobalSum(com, liquified_ice_volume);
    126 }
    127 
    128 
    129 EnergyModel::EnergyModel(IceGrid::ConstPtr grid,
    130                          stressbalance::StressBalance *stress_balance)
    131   : Component(grid), m_stress_balance(stress_balance) {
    132 
    133   const unsigned int WIDE_STENCIL = m_config->get_number("grid.max_stencil_width");
    134 
    135   {
    136     m_ice_enthalpy.create(m_grid, "enthalpy", WITH_GHOSTS, WIDE_STENCIL);
    137     // POSSIBLE standard name = land_ice_enthalpy
    138     m_ice_enthalpy.set_attrs("model_state",
    139                              "ice enthalpy (includes sensible heat, latent heat, pressure)",
    140                              "J kg-1", "J kg-1", "", 0);
    141   }
    142 
    143   {
    144     m_basal_melt_rate.create(m_grid, "basal_melt_rate_grounded", WITHOUT_GHOSTS);
    145     // ghosted to allow the "redundant" computation of tauc
    146     m_basal_melt_rate.set_attrs("model_state",
    147                                 "ice basal melt rate from energy conservation, in ice thickness per time (valid in grounded areas)",
    148                                 "m s-1", "m year-1", "", 0);
    149     // We could use land_ice_basal_melt_rate, but that way both basal_melt_rate_grounded and bmelt
    150     // have this standard name.
    151     m_basal_melt_rate.metadata().set_string("comment", "positive basal melt rate corresponds to ice loss");
    152   }
    153 
    154   // a 3d work vector
    155   {
    156     m_work.create(m_grid, "work_vector", WITHOUT_GHOSTS);
    157     m_work.set_attrs("internal",
    158                      "usually new values of temperature or enthalpy during time step",
    159                      "", "", "", 0);
    160   }
    161 }
    162 
    163 void EnergyModel::init_enthalpy(const File &input_file, bool do_regrid, int record) {
    164 
    165   if (input_file.find_variable("enthalpy")) {
    166     if (do_regrid) {
    167       m_ice_enthalpy.regrid(input_file, CRITICAL);
    168     } else {
    169       m_ice_enthalpy.read(input_file, record);
    170     }
    171   } else if (input_file.find_variable("temp")) {
    172     IceModelVec3
    173       &temp    = m_work,
    174       &liqfrac = m_ice_enthalpy;
    175 
    176     {
    177       temp.set_name("temp");
    178       temp.metadata(0).set_name("temp");
    179       temp.set_attrs("temporary", "ice temperature",
    180                      "Kelvin", "Kelvin", "land_ice_temperature", 0);
    181 
    182       if (do_regrid) {
    183         temp.regrid(input_file, CRITICAL);
    184       } else {
    185         temp.read(input_file, record);
    186       }
    187     }
    188 
    189     const IceModelVec2S & ice_thickness = *m_grid->variables().get_2d_scalar("land_ice_thickness");
    190 
    191     if (input_file.find_variable("liqfrac")) {
    192       SpatialVariableMetadata enthalpy_metadata = m_ice_enthalpy.metadata();
    193 
    194       liqfrac.set_name("liqfrac");
    195       liqfrac.metadata(0).set_name("liqfrac");
    196       liqfrac.set_attrs("temporary", "ice liquid water fraction",
    197                         "1", "1", "", 0);
    198 
    199       if (do_regrid) {
    200         liqfrac.regrid(input_file, CRITICAL);
    201       } else {
    202         liqfrac.read(input_file, record);
    203       }
    204 
    205       m_ice_enthalpy.set_name(enthalpy_metadata.get_name());
    206       m_ice_enthalpy.metadata() = enthalpy_metadata;
    207 
    208       m_log->message(2,
    209                      " - Computing enthalpy using ice temperature,"
    210                      "  liquid water fraction and thickness...\n");
    211       compute_enthalpy(temp, liqfrac, ice_thickness, m_ice_enthalpy);
    212     } else {
    213       m_log->message(2,
    214                      " - Computing enthalpy using ice temperature and thickness "
    215                      "and assuming zero liquid water fraction...\n");
    216       compute_enthalpy_cold(temp, ice_thickness, m_ice_enthalpy);
    217     }
    218   } else {
    219     throw RuntimeError::formatted(PISM_ERROR_LOCATION,
    220                                   "neither enthalpy nor temperature was found in '%s'.\n",
    221                                   input_file.filename().c_str());
    222   }
    223 }
    224 
    225 /*!
    226  * The `-regrid_file` may contain enthalpy, temperature, or *both* temperature and liquid water
    227  * fraction.
    228  */
    229 void EnergyModel::regrid_enthalpy() {
    230 
    231   auto regrid_filename = m_config->get_string("input.regrid.file");
    232   auto regrid_vars     = set_split(m_config->get_string("input.regrid.vars"), ',');
    233 
    234 
    235   if (regrid_filename.empty()) {
    236     return;
    237   }
    238 
    239   std::string enthalpy_name = m_ice_enthalpy.metadata().get_name();
    240 
    241   if (regrid_vars.empty() or member(enthalpy_name, regrid_vars)) {
    242     File regrid_file(m_grid->com, regrid_filename, PISM_GUESS, PISM_READONLY);
    243     init_enthalpy(regrid_file, true, 0);
    244   }
    245 }
    246 
    247 
    248 void EnergyModel::restart(const File &input_file, int record) {
    249   this->restart_impl(input_file, record);
    250 }
    251 
    252 void EnergyModel::bootstrap(const File &input_file,
    253                             const IceModelVec2S &ice_thickness,
    254                             const IceModelVec2S &surface_temperature,
    255                             const IceModelVec2S &climatic_mass_balance,
    256                             const IceModelVec2S &basal_heat_flux) {
    257   this->bootstrap_impl(input_file,
    258                        ice_thickness, surface_temperature,
    259                        climatic_mass_balance, basal_heat_flux);
    260 }
    261 
    262 void EnergyModel::initialize(const IceModelVec2S &basal_melt_rate,
    263                              const IceModelVec2S &ice_thickness,
    264                              const IceModelVec2S &surface_temperature,
    265                              const IceModelVec2S &climatic_mass_balance,
    266                              const IceModelVec2S &basal_heat_flux) {
    267   this->initialize_impl(basal_melt_rate,
    268                         ice_thickness,
    269                         surface_temperature,
    270                         climatic_mass_balance,
    271                         basal_heat_flux);
    272 }
    273 
    274 void EnergyModel::update(double t, double dt, const Inputs &inputs) {
    275   // reset standard out flags at the beginning of every time step
    276   m_stdout_flags = "";
    277   m_stats = EnergyModelStats();
    278 
    279   const Profiling &profiling = m_grid->ctx()->profiling();
    280 
    281   profiling.begin("ice_energy");
    282   {
    283     // this call should fill m_work with new values of enthalpy
    284     this->update_impl(t, dt, inputs);
    285 
    286     m_work.update_ghosts(m_ice_enthalpy);
    287   }
    288   profiling.end("ice_energy");
    289 
    290   // globalize m_stats and update m_stdout_flags
    291   {
    292     char buffer[50] = "";
    293 
    294     m_stats.sum(m_grid->com);
    295 
    296     if (m_stats.reduced_accuracy_counter > 0.0) { // count of when BOMBPROOF switches to lower accuracy
    297       const double reduced_accuracy_percentage = 100.0 * (m_stats.reduced_accuracy_counter / (m_grid->Mx() * m_grid->My()));
    298       const double reporting_threshold = 5.0; // only report if above 5%
    299 
    300       if (reduced_accuracy_percentage > reporting_threshold and m_log->get_threshold() > 2) {
    301         snprintf(buffer, 50, "  [BPsacr=%.4f%%] ", reduced_accuracy_percentage);
    302         m_stdout_flags = buffer + m_stdout_flags;
    303       }
    304     }
    305 
    306     if (m_stats.bulge_counter > 0) {
    307       // count of when advection bulges are limited; frequently it is identically zero
    308       snprintf(buffer, 50, " BULGE=%d ", m_stats.bulge_counter);
    309       m_stdout_flags = buffer + m_stdout_flags;
    310     }
    311   }
    312 }
    313 
    314 MaxTimestep EnergyModel::max_timestep_impl(double t) const {
    315   // silence a compiler warning
    316   (void) t;
    317 
    318   if (m_stress_balance == NULL) {
    319     throw RuntimeError::formatted(PISM_ERROR_LOCATION,
    320                                   "EnergyModel: no stress balance provided."
    321                                   " Cannot compute max. time step.");
    322   }
    323 
    324   return MaxTimestep(m_stress_balance->max_timestep_cfl_3d().dt_max.value(), "energy");
    325 }
    326 
    327 const std::string& EnergyModel::stdout_flags() const {
    328   return m_stdout_flags;
    329 }
    330 
    331 const EnergyModelStats& EnergyModel::stats() const {
    332   return m_stats;
    333 }
    334 
    335 const IceModelVec3 & EnergyModel::enthalpy() const {
    336   return m_ice_enthalpy;
    337 }
    338 
    339 /*! @brief Basal melt rate in grounded areas. (It is set to zero elsewhere.) */
    340 const IceModelVec2S & EnergyModel::basal_melt_rate() const {
    341   return m_basal_melt_rate;
    342 }
    343 
    344 /*! @brief Ice loss "flux" due to ice liquefaction. */
    345 class LiquifiedIceFlux : public TSDiag<TSFluxDiagnostic,EnergyModel> {
    346 public:
    347   LiquifiedIceFlux(const EnergyModel *m)
    348     : TSDiag<TSFluxDiagnostic, EnergyModel>(m, "liquified_ice_flux") {
    349 
    350     set_units("m3 / second", "m3 / year");
    351     m_ts.variable().set_string("long_name",
    352                                "rate of ice loss due to liquefaction,"
    353                                " averaged over the reporting interval");
    354     m_ts.variable().set_string("comment", "positive means ice loss");
    355     m_ts.variable().set_string("cell_methods", "time: mean");
    356   }
    357 protected:
    358   double compute() {
    359     // liquified ice volume during the last time step
    360     return model->stats().liquified_ice_volume;
    361   }
    362 };
    363 
    364 namespace diagnostics {
    365 /*! @brief Report ice enthalpy. */
    366 class Enthalpy : public Diag<EnergyModel>
    367 {
    368 public:
    369   Enthalpy(const EnergyModel *m)
    370     : Diag<EnergyModel>(m) {
    371     m_vars = {model->enthalpy().metadata()};
    372   }
    373 
    374 protected:
    375   IceModelVec::Ptr compute_impl() const {
    376 
    377     IceModelVec3::Ptr result(new IceModelVec3(m_grid, "enthalpy", WITHOUT_GHOSTS));
    378     result->metadata(0) = m_vars[0];
    379 
    380     const IceModelVec3 &input = model->enthalpy();
    381 
    382     // FIXME: implement IceModelVec3::copy_from()
    383 
    384     IceModelVec::AccessList list {result.get(), &input};
    385     ParallelSection loop(m_grid->com);
    386     try {
    387       for (Points p(*m_grid); p; p.next()) {
    388         const int i = p.i(), j = p.j();
    389 
    390         result->set_column(i, j, input.get_column(i, j));
    391       }
    392     } catch (...) {
    393       loop.failed();
    394     }
    395     loop.check();
    396 
    397 
    398     return result;
    399   }
    400 };
    401 
    402 } // end of namespace diagnostics
    403 
    404 DiagnosticList EnergyModel::diagnostics_impl() const {
    405   DiagnosticList result;
    406   result = {
    407     {"enthalpy",                 Diagnostic::Ptr(new diagnostics::Enthalpy(this))},
    408     {"basal_melt_rate_grounded", Diagnostic::wrap(m_basal_melt_rate)}
    409   };
    410   return result;
    411 }
    412 
    413 TSDiagnosticList EnergyModel::ts_diagnostics_impl() const {
    414   return {
    415     {"liquified_ice_flux", TSDiagnostic::Ptr(new LiquifiedIceFlux(this))}
    416   };
    417 }
    418 
    419 } // end of namespace energy
    420 
    421 } // end of namespace pism