tRegionalYieldStress.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tRegionalYieldStress.cc (4832B)
---
1 /* Copyright (C) 2015, 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 "RegionalYieldStress.hh"
21 #include "pism/util/pism_utilities.hh" // pism::combine()
22 #include "pism/util/MaxTimestep.hh"
23
24 namespace pism {
25
26 RegionalYieldStress::RegionalYieldStress(std::shared_ptr<YieldStress> input)
27 : YieldStress(input->grid()), m_input(input) {
28
29 m_high_tauc = m_config->get_number("regional.no_model_yield_stress", "Pa");
30
31 m_name = "regional " + m_input->name();
32 }
33
34 RegionalYieldStress::~RegionalYieldStress() {
35 // empty
36 }
37
38 /*!
39 * Set `basal_yield_stress` to `tauc` in areas indicated using `mask`.
40 */
41 static void set_no_model_yield_stress(double tauc,
42 const IceModelVec2Int &mask,
43 IceModelVec2S &basal_yield_stress) {
44 auto grid = mask.grid();
45
46 IceModelVec::AccessList list{&mask, &basal_yield_stress};
47
48 for (Points p(*grid); p; p.next()) {
49 const int i = p.i(), j = p.j();
50
51 if (mask(i, j) > 0.5) {
52 basal_yield_stress(i, j) = tauc;
53 }
54 }
55 }
56
57 void RegionalYieldStress::restart_impl(const File &input_file, int record) {
58 m_input->restart(input_file, record);
59
60 // m_basal_yield_stress is a part of the model state for all yield stress models, so the
61 // call above should read it in.
62 m_basal_yield_stress.copy_from(m_input->basal_material_yield_stress());
63
64 IceModelVec2Int no_model_mask(m_grid, "no_model_mask", WITHOUT_GHOSTS);
65 no_model_mask.set_attrs("model_state",
66 "mask: zeros (modeling domain) and ones"
67 " (no-model buffer near grid edges)",
68 "", "", "", 0); // no units and no standard name
69 // We are re-starting a simulation, so the input file has to contain no_model_mask.
70 no_model_mask.read(input_file, record);
71 // However, the used can set "-regrid_vars no_model_mask,...", so we have to try this,
72 // too.
73 regrid(name(), no_model_mask);
74
75 set_no_model_yield_stress(m_high_tauc, no_model_mask, m_basal_yield_stress);
76 }
77
78 void RegionalYieldStress::bootstrap_impl(const File &input_file, const YieldStressInputs &inputs) {
79 m_input->bootstrap(input_file, inputs);
80
81 m_basal_yield_stress.copy_from(m_input->basal_material_yield_stress());
82
83 set_no_model_yield_stress(m_high_tauc, *inputs.no_model_mask, m_basal_yield_stress);
84 }
85
86 void RegionalYieldStress::init_impl(const YieldStressInputs &inputs) {
87 m_input->init(inputs);
88
89 m_basal_yield_stress.copy_from(m_input->basal_material_yield_stress());
90
91 set_no_model_yield_stress(m_high_tauc, *inputs.no_model_mask, m_basal_yield_stress);
92 }
93
94 void RegionalYieldStress::update_impl(const YieldStressInputs &inputs,
95 double t, double dt) {
96 m_input->update(inputs, t, dt);
97
98 m_basal_yield_stress.copy_from(m_input->basal_material_yield_stress());
99
100 set_no_model_yield_stress(m_high_tauc, *inputs.no_model_mask, m_basal_yield_stress);
101 }
102
103 void RegionalYieldStress::define_model_state_impl(const File &output) const {
104 m_input->define_model_state(output);
105
106 // define tauc (this is likely to be a no-op because m_input should have defined it by
107 // now)
108 m_basal_yield_stress.define(output);
109 }
110
111 void RegionalYieldStress::write_model_state_impl(const File &output) const {
112 m_input->write_model_state(output);
113 // Write basal yield stress that includes the modification containing high yield stress
114 // in "no model" areas, overwriting the field written by m_input.
115 m_basal_yield_stress.write(output);
116 }
117
118 DiagnosticList RegionalYieldStress::diagnostics_impl() const {
119 // Override the tauc diagnostic with the one that includes the regional modification
120 return combine({{"tauc", Diagnostic::wrap(m_basal_yield_stress)}},
121 m_input->diagnostics());
122 }
123
124 MaxTimestep RegionalYieldStress::max_timestep_impl(double t) const {
125 auto dt = m_input->max_timestep(t);
126
127 if (dt.finite()) {
128 return MaxTimestep(dt.value(), name());
129 }
130
131 return MaxTimestep(name());
132 }
133
134 TSDiagnosticList RegionalYieldStress::ts_diagnostics_impl() const {
135 return m_input->ts_diagnostics();
136 }
137
138 } // end of namespace pism