Loading [MathJax]/jax/output/HTML-CSS/config.js
PISM, A Parallel Ice Sheet Model 2.2.1-cd005eec8 committed by Constantine Khrulev on 2025-03-07
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
timestepping.cc
Go to the documentation of this file.
1// Copyright (C) 2004-2017, 2019, 2020, 2021, 2022, 2023 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 <algorithm> // std::sort
20#include <cmath> // std::floor
21#include <cassert>
22
23#include "pism/icemodel/IceModel.hh"
24#include "pism/util/Grid.hh"
25#include "pism/util/ConfigInterface.hh"
26#include "pism/util/Time.hh"
27#include "pism/util/MaxTimestep.hh"
28#include "pism/stressbalance/StressBalance.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/coupler/FrontalMelt.hh"
37
38namespace pism {
39
40//! Compute the maximum time step allowed by the diffusive SIA.
41/*!
42If maximum diffusivity is positive (i.e. if there is diffusion going on) then
43updates dt.
44
45Note adapt_ratio * 2 is multiplied by dx^2/(2*maxD) so dt <= adapt_ratio *
46dx^2/maxD (if dx=dy).
47
48Reference: [\ref MortonMayers] pp 62--63.
49 */
51 double D_max = m_stress_balance->max_diffusivity();
52
53 double dx = m_grid->dx(), dy = m_grid->dy(),
54 adaptive_timestepping_ratio = m_config->get_number("time_stepping.adaptive_ratio");
55
56 auto dt_diffusivity = ::pism::max_timestep_diffusivity(D_max, dx, dy, adaptive_timestepping_ratio);
57
58 MaxTimestep dt_max(m_config->get_number("time_stepping.maximum_time_step", "seconds"),
59 "max time step");
60
61 return std::min(dt_diffusivity, dt_max);
62}
63
64/** @brief Compute the skip counter using "long" (usually determined
65 * using the CFL stability criterion) and "short" (typically
66 * determined using the diffusivity-based stability criterion) time
67 * step lengths.
68 *
69 *
70 * @param[in] dt "long" time-step
71 * @param[in] dt_diffusivity "short" time-step
72 *
73 * @return new skip counter
74 */
75unsigned int IceModel::skip_counter(double dt, double dt_diffusivity) {
76
77 if (not m_config->get_flag("time_stepping.skip.enabled")) {
78 return 0;
79 }
80
81 const int skip_max = static_cast<int>(m_config->get_number("time_stepping.skip.max"));
82
83 if (dt_diffusivity > 0.0) {
84 const double conservative_factor = 0.95;
85 const double counter = floor(conservative_factor * (dt / dt_diffusivity));
86 return std::min(static_cast<int>(counter), skip_max);
87 }
88
89 return skip_max;
90}
91
92//! Use various stability criteria to determine the time step for an evolution run.
93/*!
94The main loop in run() approximates many physical processes. Several of these approximations,
95including the mass continuity and temperature equations in particular, involve stability
96criteria. This procedure builds the length of the next time step by using these criteria and
97by incorporating choices made by options (e.g. `-max_dt`) and by derived classes.
98
99@param[in] counter current time-step skipping counter
100 */
102
103 const double current_time = m_time->current();
104
105 std::vector<MaxTimestep> restrictions;
106
107 // get time-stepping restrictions from sub-models
108 for (auto m : m_submodels) {
109 restrictions.push_back(m.second->max_timestep(current_time));
110 }
111
112 // mechanisms that use a retreat rate
113 bool front_retreat =
115 if (front_retreat and m_config->get_flag("geometry.front_retreat.use_cfl")) {
116 // at least one of front retreat mechanisms is active *and* PISM is told to use a CFL
117 // restriction
118
119 array::Scalar &retreat_rate = *m_work2d[0];
120 retreat_rate.set(0.0);
121
122 if (m_eigen_calving) {
123 retreat_rate.add(1.0, m_eigen_calving->calving_rate());
124 }
125
126 if (m_hayhurst_calving) {
127 retreat_rate.add(1.0, m_hayhurst_calving->calving_rate());
128 }
129
130 if (m_vonmises_calving) {
131 retreat_rate.add(1.0, m_vonmises_calving->calving_rate());
132 }
133
134 if (m_frontal_melt) {
135 retreat_rate.add(1.0, m_frontal_melt->retreat_rate());
136 }
137
138 assert(m_front_retreat);
139
140 restrictions.push_back(
141 m_front_retreat->max_timestep(m_geometry.cell_type, m_ice_thickness_bc_mask, retreat_rate));
142 }
143
144 const char *end = "end of the run";
145 const char *max = "max";
146
147 // Always consider the maximum allowed time-step length.
148 double max_timestep = m_config->get_number("time_stepping.maximum_time_step", "seconds");
149 if (max_timestep > 0.0) {
150 restrictions.push_back(MaxTimestep(max_timestep, max));
151 }
152
153 // Never go past the end of a run.
154 const double time_to_end = m_time->end() - current_time;
155 if (time_to_end > 0.0) {
156 restrictions.push_back(MaxTimestep(time_to_end, end));
157 }
158
159 // reporting
160 {
161 restrictions.push_back(ts_max_timestep(current_time));
162 restrictions.push_back(extras_max_timestep(current_time));
163 restrictions.push_back(save_max_timestep(current_time));
164 }
165
166 // mass continuity stability criteria
167 if (m_config->get_flag("geometry.update.enabled")) {
168 auto cfl = m_stress_balance->max_timestep_cfl_2d();
169
170 restrictions.push_back(MaxTimestep(cfl.dt_max.value(), "2D CFL"));
171 restrictions.push_back(max_timestep_diffusivity());
172 }
173
174 // sort time step restrictions to find the strictest one
175 std::sort(restrictions.begin(), restrictions.end());
176
177 // note that restrictions has at least 2 elements
178 // the first element is the max time step we can take
179 assert(restrictions.size() >= 2);
180 auto dt_max = restrictions[0];
181 auto dt_other = restrictions[1];
182
183 TimesteppingInfo result;
184 result.dt = dt_max.value();
185 result.reason = (dt_max.description() + " (overrides " + dt_other.description() + ")");
186 result.skip_counter = 0;
187
188 double resolution = m_config->get_number("time_stepping.resolution", "seconds");
189
190 // Hit all multiples of X years, if requested.
191 {
192 int year_increment = static_cast<int>(m_config->get_number("time_stepping.hit_multiples"));
193 if (year_increment > 0) {
194 auto next_time = m_time->increment_date(m_timestep_hit_multiples_last_time,
195 year_increment);
196
197 if (std::fabs(current_time - next_time) < resolution) {
198 // the current time is a multiple of year_increment
200 next_time = m_time->increment_date(current_time, year_increment);
201 }
202
203 auto dt = next_time - current_time;
204 assert(dt > resolution);
205
206 if (dt < result.dt) {
207 result.dt = dt;
208 result.reason = pism::printf("hit multiples of %d years (overrides %s)",
209 year_increment, dt_max.description().c_str());
210 }
211 }
212 }
213
214 // the "skipping" mechanism
215 {
216 if (dt_max.description() == "diffusivity" and counter == 0) {
217 result.skip_counter = skip_counter(dt_other.value(), dt_max.value());
218 } else {
219 result.skip_counter = counter;
220 }
221
222 // "max" and "end of the run" limit the "big" time-step (in
223 // the context of the "skipping" mechanism), so we might need to
224 // reset the skip_counter_result to 1.
225 if (member(dt_max.description(), {max, end}) and counter > 1) {
226 result.skip_counter = 1;
227 }
228 }
229
230 if (resolution > 0.0) {
231 double dt = std::floor(result.dt * resolution) / resolution;
232
233 // Ensure that the resulting time step is never zero. This may happen if the length of
234 // the run is not an integer multiple of "resolution".
235 if (dt >= resolution) {
236 result.dt = dt;
237 }
238 }
239
240 return result;
241}
242
243} // end of namespace pism
array::CellType2 cell_type
Definition Geometry.hh:55
std::map< std::string, const Component * > m_submodels
the list of sub-models, for writing model states and obtaining diagnostics
Definition IceModel.hh:252
MaxTimestep extras_max_timestep(double my_t)
Computes the maximum time-step we can take and still hit all -extra_times.
std::shared_ptr< stressbalance::StressBalance > m_stress_balance
Definition IceModel.hh:395
const Time::Ptr m_time
Time manager.
Definition IceModel.hh:246
Geometry m_geometry
Definition IceModel.hh:288
double m_timestep_hit_multiples_last_time
Definition IceModel.hh:465
std::shared_ptr< frontalmelt::FrontalMelt > m_frontal_melt
Definition IceModel.hh:281
virtual unsigned int skip_counter(double input_dt, double input_dt_diffusivity)
Compute the skip counter using "long" (usually determined using the CFL stability criterion) and "sho...
virtual MaxTimestep max_timestep_diffusivity()
Compute the maximum time step allowed by the diffusive SIA.
MaxTimestep ts_max_timestep(double my_t)
Computes the maximum time-step we can take and still hit all -ts_times.
Definition output_ts.cc:110
std::shared_ptr< calving::EigenCalving > m_eigen_calving
Definition IceModel.hh:268
std::vector< std::shared_ptr< array::Scalar2 > > m_work2d
Definition IceModel.hh:393
double dt() const
Definition IceModel.cc:126
std::shared_ptr< calving::HayhurstCalving > m_hayhurst_calving
Definition IceModel.hh:269
Config::Ptr m_config
Configuration flags and parameters.
Definition IceModel.hh:238
array::Scalar1 m_ice_thickness_bc_mask
Mask prescribing locations where ice thickness is held constant.
Definition IceModel.hh:307
MaxTimestep save_max_timestep(double my_t)
Computes the maximum time-step we can take and still hit all -save_times.
std::shared_ptr< calving::vonMisesCalving > m_vonmises_calving
Definition IceModel.hh:270
const std::shared_ptr< Grid > m_grid
Computational grid.
Definition IceModel.hh:236
virtual TimesteppingInfo max_timestep(unsigned int counter)
Use various stability criteria to determine the time step for an evolution run.
std::shared_ptr< FrontRetreat > m_front_retreat
Definition IceModel.hh:277
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
void add(double alpha, const Array2D< T > &x)
Definition Array2D.hh:65
void set(double c)
Result: v[j] <- c for all j.
Definition Array.cc:629
std::string printf(const char *format,...)
bool member(const std::string &string, const std::set< std::string > &set)
MaxTimestep max_timestep_diffusivity(double D_max, double dx, double dy, double adaptive_timestepping_ratio)