ttimestepping.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
---
ttimestepping.cc (8341B)
---
1 // Copyright (C) 2004-2017, 2019 Jed Brown, Ed Bueler and Constantine Khroulev
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 #include <sstream> // stringstream
20 #include <algorithm> // std::sort
21
22 #include "IceModel.hh"
23 #include "pism/util/IceGrid.hh"
24 #include "pism/util/ConfigInterface.hh"
25 #include "pism/util/Time.hh"
26 #include "pism/util/MaxTimestep.hh"
27 #include "pism/stressbalance/StressBalance.hh"
28 #include "pism/stressbalance/ShallowStressBalance.hh"
29 #include "pism/util/Component.hh" // ...->max_timestep()
30
31 #include "pism/frontretreat/calving/EigenCalving.hh"
32 #include "pism/frontretreat/calving/HayhurstCalving.hh"
33 #include "pism/frontretreat/calving/vonMisesCalving.hh"
34 #include "pism/frontretreat/FrontRetreat.hh"
35
36 #include "pism/energy/EnergyModel.hh"
37 #include "pism/coupler/OceanModel.hh"
38 #include "pism/coupler/FrontalMelt.hh"
39
40 namespace pism {
41
42 //! Compute the maximum time step allowed by the diffusive SIA.
43 /*!
44 If maximum diffusivity is positive (i.e. if there is diffusion going on) then
45 updates dt.
46
47 Note adapt_ratio * 2 is multiplied by dx^2/(2*maxD) so dt <= adapt_ratio *
48 dx^2/maxD (if dx=dy).
49
50 Reference: [\ref MortonMayers] pp 62--63.
51 */
52 MaxTimestep IceModel::max_timestep_diffusivity() {
53 double D_max = m_stress_balance->max_diffusivity();
54
55 if (D_max > 0.0) {
56 const double
57 dx = m_grid->dx(),
58 dy = m_grid->dy(),
59 adaptive_timestepping_ratio = m_config->get_number("time_stepping.adaptive_ratio"),
60 grid_factor = 1.0 / (dx*dx) + 1.0 / (dy*dy);
61
62 return MaxTimestep(adaptive_timestepping_ratio * 2.0 / (D_max * grid_factor),
63 "diffusivity");
64 } else {
65 return MaxTimestep(m_config->get_number("time_stepping.maximum_time_step", "seconds"),
66 "max time step");
67 }
68 }
69
70 /** @brief Compute the skip counter using "long" (usually determined
71 * using the CFL stability criterion) and "short" (typically
72 * determined using the diffusivity-based stability criterion) time
73 * step lengths.
74 *
75 *
76 * @param[in] input_dt long time-step
77 * @param[in] input_dt_diffusivity short time-step
78 *
79 * @return new skip counter
80 */
81 unsigned int IceModel::skip_counter(double input_dt, double input_dt_diffusivity) {
82
83 if (not m_config->get_flag("time_stepping.skip.enabled")) {
84 return 0;
85 }
86
87 const unsigned int skip_max = static_cast<int>(m_config->get_number("time_stepping.skip.max"));
88
89 if (input_dt_diffusivity > 0.0) {
90 const double conservativeFactor = 0.95;
91 const double counter = floor(conservativeFactor * (input_dt / input_dt_diffusivity));
92 const unsigned int result = static_cast<unsigned int>(counter);
93 return std::min(result, skip_max);
94 } else {
95 return skip_max;
96 }
97
98 return 0;
99 }
100
101 //! Use various stability criteria to determine the time step for an evolution run.
102 /*!
103 The main loop in run() approximates many physical processes. Several of these approximations,
104 including the mass continuity and temperature equations in particular, involve stability
105 criteria. This procedure builds the length of the next time step by using these criteria and
106 by incorporating choices made by options (e.g. <c>-max_dt</c>) and by derived classes.
107
108 @param[out] dt_result computed maximum time step
109 @param[in,out] skip_counter_result time-step skipping counter
110 */
111 void IceModel::max_timestep(double &dt_result, unsigned int &skip_counter_result) {
112
113 const double current_time = m_time->current();
114
115 std::vector<MaxTimestep> restrictions;
116
117 // get time-stepping restrictions from sub-models
118 for (auto m : m_submodels) {
119 restrictions.push_back(m.second->max_timestep(current_time));
120 }
121
122 // mechanisms that use a retreat rate
123 if (m_config->get_flag("geometry.front_retreat.use_cfl") and
124 (m_eigen_calving or m_vonmises_calving or m_hayhurst_calving or m_frontal_melt)) {
125 // at least one of front retreat mechanisms is active
126
127 IceModelVec2S &retreat_rate = m_work2d[0];
128 retreat_rate.set(0.0);
129
130 if (m_eigen_calving) {
131 retreat_rate.add(1.0, m_eigen_calving->calving_rate());
132 }
133
134 if (m_hayhurst_calving) {
135 retreat_rate.add(1.0, m_hayhurst_calving->calving_rate());
136 }
137
138 if (m_vonmises_calving) {
139 retreat_rate.add(1.0, m_vonmises_calving->calving_rate());
140 }
141
142 if (m_frontal_melt) {
143 retreat_rate.add(1.0, m_frontal_melt->retreat_rate());
144 }
145
146 assert(m_front_retreat);
147
148 restrictions.push_back(m_front_retreat->max_timestep(m_geometry.cell_type,
149 m_ssa_dirichlet_bc_mask,
150 retreat_rate));
151 }
152
153 // Always consider the maximum allowed time-step length.
154 if (m_config->get_number("time_stepping.maximum_time_step") > 0.0) {
155 restrictions.push_back(MaxTimestep(m_config->get_number("time_stepping.maximum_time_step",
156 "seconds"),
157 "max"));
158 }
159
160 // Never go past the end of a run.
161 const double time_to_end = m_time->end() - current_time;
162 if (time_to_end > 0.0) {
163 restrictions.push_back(MaxTimestep(time_to_end, "end of the run"));
164 }
165
166 // reporting
167 {
168 restrictions.push_back(ts_max_timestep(current_time));
169 restrictions.push_back(extras_max_timestep(current_time));
170 restrictions.push_back(save_max_timestep(current_time));
171 }
172
173 // mass continuity stability criteria
174 if (m_config->get_flag("geometry.update.enabled")) {
175 CFLData cfl = m_stress_balance->max_timestep_cfl_2d();
176
177 restrictions.push_back(MaxTimestep(cfl.dt_max.value(), "2D CFL"));
178 restrictions.push_back(max_timestep_diffusivity());
179 }
180
181 // Hit multiples of X years, if requested.
182 {
183 const int timestep_hit_multiples = static_cast<int>(m_config->get_number("time_stepping.hit_multiples"));
184 if (timestep_hit_multiples > 0) {
185 const double epsilon = 1.0; // 1 second tolerance
186 double
187 next_time = m_timestep_hit_multiples_last_time;
188
189 while (m_time->increment_date(next_time, timestep_hit_multiples) <= current_time + dt_result + epsilon) {
190 next_time = m_time->increment_date(next_time, timestep_hit_multiples);
191 }
192
193 if (next_time > current_time && next_time <= current_time + dt_result + epsilon) {
194 dt_result = next_time - current_time;
195 m_timestep_hit_multiples_last_time = next_time;
196
197 std::stringstream str;
198 str << "hit multiples of " << timestep_hit_multiples << " years";
199
200 restrictions.push_back(MaxTimestep(next_time - current_time, str.str()));
201
202 }
203 }
204 }
205
206 // sort time step restrictions to find the strictest one
207 std::sort(restrictions.begin(), restrictions.end());
208
209 // note that restrictions has at least 2 elements
210 // the first element is the max time step we can take
211 MaxTimestep dt_max = restrictions[0];
212 MaxTimestep dt_other = restrictions[1];
213 dt_result = dt_max.value();
214 m_adaptive_timestep_reason = (dt_max.description() +
215 " (overrides " + dt_other.description() + ")");
216
217 // the "skipping" mechanism
218 {
219 if (dt_max.description() == "diffusivity" and skip_counter_result == 0) {
220 skip_counter_result = skip_counter(dt_other.value(), dt_max.value());
221 }
222
223 // "max" and "end of the run" limit the "big" time-step (in
224 // the context of the "skipping" mechanism), so we might need to
225 // reset the skip_counter_result to 1.
226 if ((m_adaptive_timestep_reason == "max" ||
227 m_adaptive_timestep_reason == "end of the run") &&
228 skip_counter_result > 1) {
229 skip_counter_result = 1;
230 }
231 }
232 }
233
234 } // end of namespace pism