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 (4930B)
---
1 /* Copyright (C) 2016, 2017 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 "timestepping.hh"
21 #include "pism/util/IceGrid.hh"
22 #include "pism/util/iceModelVec.hh"
23 #include "pism/util/IceModelVec2CellType.hh"
24 #include "pism/util/pism_utilities.hh"
25
26 namespace pism {
27
28 CFLData::CFLData() {
29 u_max = 0.0;
30 v_max = 0.0;
31 w_max = 0.0;
32 }
33
34 //! Compute the maximum velocities for time-stepping and reporting to user.
35 /*!
36 Computes the maximum magnitude of the components \f$u,v,w\f$ of the 3D velocity.
37
38 Under BOMBPROOF there is no CFL condition for the vertical advection.
39 The maximum vertical velocity is computed but it does not affect the output.
40 */
41 CFLData max_timestep_cfl_3d(const IceModelVec2S &ice_thickness,
42 const IceModelVec2CellType &cell_type,
43 const IceModelVec3 &u3,
44 const IceModelVec3 &v3,
45 const IceModelVec3 &w3) {
46
47 IceGrid::ConstPtr grid = ice_thickness.grid();
48 Config::ConstPtr config = grid->ctx()->config();
49
50 double dt_max = config->get_number("time_stepping.maximum_time_step", "seconds");
51
52 IceModelVec::AccessList list{&ice_thickness, &u3, &v3, &w3, &cell_type};
53
54 // update global max of abs of velocities for CFL; only velocities under surface
55 const double
56 one_over_dx = 1.0 / grid->dx(),
57 one_over_dy = 1.0 / grid->dy();
58
59 double u_max = 0.0, v_max = 0.0, w_max = 0.0;
60 ParallelSection loop(grid->com);
61 try {
62 for (Points p(*grid); p; p.next()) {
63 const int i = p.i(), j = p.j();
64
65 if (cell_type.icy(i, j)) {
66 const int ks = grid->kBelowHeight(ice_thickness(i, j));
67 const double
68 *u = u3.get_column(i, j),
69 *v = v3.get_column(i, j),
70 *w = w3.get_column(i, j);
71
72 for (int k = 0; k <= ks; ++k) {
73 const double
74 u_abs = fabs(u[k]),
75 v_abs = fabs(v[k]);
76 u_max = std::max(u_max, u_abs);
77 v_max = std::max(v_max, v_abs);
78 const double denom = fabs(u_abs * one_over_dx) + fabs(v_abs * one_over_dy);
79 if (denom > 0.0) {
80 dt_max = std::min(dt_max, 1.0 / denom);
81 }
82 }
83
84 for (int k = 0; k <= ks; ++k) {
85 w_max = std::max(w_max, fabs(w[k]));
86 }
87 }
88 }
89 } catch (...) {
90 loop.failed();
91 }
92 loop.check();
93
94 CFLData result;
95
96 result.u_max = GlobalMax(grid->com, u_max);
97 result.v_max = GlobalMax(grid->com, v_max);
98 result.w_max = GlobalMax(grid->com, w_max);
99 result.dt_max = MaxTimestep(GlobalMin(grid->com, dt_max));
100
101 return result;
102 }
103
104 //! Compute the CFL constant associated to first-order upwinding for the sliding contribution to mass continuity.
105 /*!
106 This procedure computes the maximum horizontal speed in the icy
107 areas. In particular it computes CFL constant for the upwinding, in
108 GeometryEvolution::step(), which applies to the basal component of mass
109 flux.
110
111 That is, because the map-plane mass continuity is advective in the
112 sliding case we have a CFL condition.
113 */
114 CFLData max_timestep_cfl_2d(const IceModelVec2S &ice_thickness,
115 const IceModelVec2CellType &cell_type,
116 const IceModelVec2V &velocity) {
117
118 IceGrid::ConstPtr grid = ice_thickness.grid();
119 Config::ConstPtr config = grid->ctx()->config();
120
121 double dt_max = config->get_number("time_stepping.maximum_time_step", "seconds");
122
123 const double
124 dx = grid->dx(),
125 dy = grid->dy();
126
127 IceModelVec::AccessList list{&velocity, &cell_type};
128
129 double u_max = 0.0, v_max = 0.0;
130 for (Points p(*grid); p; p.next()) {
131 const int i = p.i(), j = p.j();
132
133 if (cell_type.icy(i, j)) {
134 const double
135 u_abs = fabs(velocity(i, j).u),
136 v_abs = fabs(velocity(i, j).v);
137
138 u_max = std::max(u_max, u_abs);
139 v_max = std::max(v_max, v_abs);
140
141 const double denom = u_abs / dx + v_abs / dy;
142 if (denom > 0.0) {
143 dt_max = std::min(dt_max, 1.0 / denom);
144 }
145 }
146 }
147
148 CFLData result;
149
150 result.u_max = GlobalMax(grid->com, u_max);
151 result.v_max = GlobalMax(grid->com, v_max);
152 result.w_max = 0.0;
153 result.dt_max = MaxTimestep(GlobalMin(grid->com, dt_max));
154
155 return result;
156 }
157
158 } // end of namespace pism