PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
RegionalYieldStress.cc
Go to the documentation of this file.
1 /* Copyright (C) 2015, 2016, 2017, 2018, 2019, 2022, 2023 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 "pism/regional/RegionalYieldStress.hh"
21 #include "pism/util/pism_utilities.hh" // pism::combine()
22 #include "pism/util/MaxTimestep.hh"
23 #include "pism/util/array/Scalar.hh"
24 #include "pism/util/interpolation.hh"
25 
26 namespace pism {
27 
28 RegionalYieldStress::RegionalYieldStress(std::shared_ptr<YieldStress> input)
29  : YieldStress(input->grid()), m_input(input) {
30 
31  m_high_tauc = m_config->get_number("regional.no_model_yield_stress", "Pa");
32 
33  m_name = "regional " + m_input->name();
34 }
35 
36 /*!
37  * Set `basal_yield_stress` to `tauc` in areas indicated using `mask`.
38  */
39 static void set_no_model_yield_stress(double tauc,
40  const array::Scalar &mask,
41  array::Scalar &basal_yield_stress) {
42  auto grid = mask.grid();
43 
44  array::AccessScope list{&mask, &basal_yield_stress};
45 
46  for (auto p = grid->points(); p; p.next()) {
47  const int i = p.i(), j = p.j();
48 
49  if (mask(i, j) > 0.5) {
50  basal_yield_stress(i, j) = tauc;
51  }
52  }
53 }
54 
55 void RegionalYieldStress::restart_impl(const File &input_file, int record) {
56  m_input->restart(input_file, record);
57 
58  // m_basal_yield_stress is a part of the model state for all yield stress models, so the
59  // call above should read it in.
60  m_basal_yield_stress.copy_from(m_input->basal_material_yield_stress());
61 
62  array::Scalar no_model_mask(m_grid, "no_model_mask");
63  no_model_mask.set_interpolation_type(NEAREST);
64  no_model_mask.metadata(0)
65  .long_name(
66  "mask: zeros (modeling domain) and ones (no-model buffer near grid edges)"); // no units and no standard name
67  // We are re-starting a simulation, so the input file has to contain no_model_mask.
68  no_model_mask.read(input_file, record);
69  // However, the used can set "-regrid_vars no_model_mask,...", so we have to try this,
70  // too.
71  regrid(name(), no_model_mask);
72 
74 }
75 
76 void RegionalYieldStress::bootstrap_impl(const File &input_file, const YieldStressInputs &inputs) {
77  m_input->bootstrap(input_file, inputs);
78 
79  m_basal_yield_stress.copy_from(m_input->basal_material_yield_stress());
80 
82 }
83 
85  m_input->init(inputs);
86 
87  m_basal_yield_stress.copy_from(m_input->basal_material_yield_stress());
88 
90 }
91 
93  double t, double dt) {
94  m_input->update(inputs, t, dt);
95 
96  m_basal_yield_stress.copy_from(m_input->basal_material_yield_stress());
97 
99 }
100 
102  m_input->define_model_state(output);
103 
104  // define tauc (this is likely to be a no-op because m_input should have defined it by
105  // now)
107 }
108 
110  m_input->write_model_state(output);
111  // Write basal yield stress that includes the modification containing high yield stress
112  // in "no model" areas, overwriting the field written by m_input.
113  m_basal_yield_stress.write(output);
114 }
115 
117  // Override the tauc diagnostic with the one that includes the regional modification
118  return combine({{"tauc", Diagnostic::wrap(m_basal_yield_stress)}},
119  m_input->diagnostics());
120 }
121 
123  auto dt = m_input->max_timestep(t);
124 
125  if (dt.finite()) {
126  return MaxTimestep(dt.value(), name());
127  }
128 
129  return MaxTimestep(name());
130 }
131 
133  return m_input->ts_diagnostics();
134 }
135 
136 } // end of namespace pism
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:158
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition: Component.hh:156
void regrid(const std::string &module_name, array::Array &variable, RegriddingFlag flag=NO_REGRID_WITHOUT_REGRID_VARS)
Definition: Component.cc:159
static Ptr wrap(const T &input)
Definition: Diagnostic.hh:160
High-level PISM I/O class.
Definition: File.hh:56
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
Definition: MaxTimestep.hh:31
std::shared_ptr< YieldStress > m_input
TSDiagnosticList ts_diagnostics_impl() const
DiagnosticList diagnostics_impl() const
void write_model_state_impl(const File &output) const
The default (empty implementation).
MaxTimestep max_timestep_impl(double t) const
void update_impl(const YieldStressInputs &inputs, double t, double dt)
void restart_impl(const File &input_file, int record)
RegionalYieldStress(std::shared_ptr< YieldStress > input)
void bootstrap_impl(const File &input_file, const YieldStressInputs &inputs)
void init_impl(const YieldStressInputs &inputs)
void define_model_state_impl(const File &output) const
VariableMetadata & long_name(const std::string &input)
const array::Scalar * no_model_mask
Definition: YieldStress.hh:39
std::string m_name
Definition: YieldStress.hh:76
std::string name() const
Definition: YieldStress.cc:108
array::Scalar2 m_basal_yield_stress
Definition: YieldStress.hh:74
The PISM basal yield stress model interface (virtual base class)
Definition: YieldStress.hh:43
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
void copy_from(const Array2D< T > &source)
Definition: Array2D.hh:73
void read(const std::string &filename, unsigned int time)
Definition: Array.cc:809
void set_interpolation_type(InterpolationType type)
Definition: Array.cc:179
void define(const File &file, io::Type default_type) const
Define variables corresponding to an Array in a file opened using file.
Definition: Array.cc:540
void write(const std::string &filename) const
Definition: Array.cc:800
std::shared_ptr< const Grid > grid() const
Definition: Array.cc:132
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: Array.cc:553
@ PISM_DOUBLE
Definition: IO_Flags.hh:52
static void set_no_model_yield_stress(double tauc, const array::Scalar &mask, array::Scalar &basal_yield_stress)
std::map< std::string, TSDiagnostic::Ptr > TSDiagnosticList
Definition: Diagnostic.hh:343
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:125
T combine(const T &a, const T &b)