tprintout.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
---
tprintout.cc (9169B)
---
1 // Copyright (C) 2004-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 <cstring>
20 #include <cstdlib>
21
22 #include <petscsys.h>
23
24 #include "IceModel.hh"
25
26 #include "pism/stressbalance/StressBalance.hh"
27
28 #include "pism/util/IceGrid.hh"
29 #include "pism/util/ConfigInterface.hh"
30 #include "pism/util/Time.hh"
31
32 #include "pism/util/pism_utilities.hh"
33
34 namespace pism {
35
36
37 //! Because of the -skip mechanism it is still possible that we can have CFL violations: count them.
38 /*! This applies to the horizontal part of the 3D advection problem solved by AgeModel and the
39 horizontal part of the 3D convection-diffusion problems solved by EnthalpyModel and
40 TemperatureModel.
41 */
42 unsigned int count_CFL_violations(const IceModelVec3 &u3,
43 const IceModelVec3 &v3,
44 const IceModelVec2S &ice_thickness,
45 double dt) {
46
47 if (dt == 0.0) {
48 return 0;
49 }
50
51 IceGrid::ConstPtr grid = u3.grid();
52
53 const double
54 CFL_x = grid->dx() / dt,
55 CFL_y = grid->dy() / dt;
56
57 IceModelVec::AccessList list{&ice_thickness, &u3, &v3};
58
59 unsigned int CFL_violation_count = 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 const int ks = grid->kBelowHeight(ice_thickness(i,j));
66
67 const double
68 *u = u3.get_column(i, j),
69 *v = v3.get_column(i, j);
70
71 // check horizontal CFL conditions at each point
72 for (int k = 0; k <= ks; k++) {
73 if (fabs(u[k]) > CFL_x) {
74 CFL_violation_count += 1;
75 }
76 if (fabs(v[k]) > CFL_y) {
77 CFL_violation_count += 1;
78 }
79 }
80 }
81 } catch (...) {
82 loop.failed();
83 }
84 loop.check();
85
86 return (unsigned int)GlobalMax(grid->com, CFL_violation_count);
87 }
88
89 void IceModel::print_summary(bool tempAndAge) {
90
91 const IceModelVec3
92 &u3 = m_stress_balance->velocity_u(),
93 &v3 = m_stress_balance->velocity_v();
94
95 unsigned int n_CFL_violations = count_CFL_violations(u3, v3, m_geometry.ice_thickness,
96 tempAndAge ? dt_TempAge : m_dt);
97
98 // report CFL violations
99 if (n_CFL_violations > 0.0) {
100 const double CFLviolpercent = 100.0 * n_CFL_violations / (m_grid->Mx() * m_grid->My() * m_grid->Mz());
101 // at default verbosity level, only report CFL viols if above:
102 const double CFLVIOL_REPORT_VERB2_PERCENT = 0.1;
103 if (CFLviolpercent > CFLVIOL_REPORT_VERB2_PERCENT ||
104 m_log->get_threshold() > 2) {
105 char tempstr[90] = "";
106 snprintf(tempstr,90,
107 " [!CFL#=%d (=%5.2f%% of 3D grid)] ",
108 n_CFL_violations,CFLviolpercent);
109 m_stdout_flags = tempstr + m_stdout_flags;
110 }
111 }
112
113 // get maximum diffusivity
114 double max_diffusivity = m_stress_balance->max_diffusivity();
115 // get volumes in m^3 and areas in m^2
116 double volume = ice_volume(m_geometry, 0.0);
117 double area = ice_area(m_geometry, 0.0);
118
119 double meltfrac = 0.0;
120 if (tempAndAge or m_log->get_threshold() >= 3) {
121 meltfrac = compute_temperate_base_fraction(area);
122 }
123
124 // main report: 'S' line
125 print_summary_line(false, tempAndAge, m_dt,
126 volume, area, meltfrac, max_diffusivity);
127 }
128
129
130 //! Print a line to stdout which summarizes the state of the modeled ice sheet at the end of the time step.
131 /*!
132 This method is for casual inspection of model behavior, and to provide the user
133 with some indication of the state of the run.
134
135 Generally, two lines are printed to stdout, the first starting with a space
136 and the second starting with the character 'S' in the left-most column (column 1).
137
138 The first line shows flags for which processes executed, and the length of the
139 time step (and/or substeps under option -skip). See IceModel::run()
140 for meaning of these flags.
141
142 If printPrototype is TRUE then the first line does not appear and
143 the second line has alternate appearance. Specifically, different column 1
144 characters are printed:
145 - 'P' line gives names of the quantities reported in the 'S' line, the
146 "prototype", while
147 - 'U' line gives units of these quantities.
148 This column 1 convention allows automatic tools to read PISM stdout
149 and produce time-series. The 'P' and 'U' lines are intended to appear once at
150 the beginning of the run, while an 'S' line appears at every time step.
151
152 These quantities are reported in this base class version:
153 - `time` is the current model time
154 - `ivol` is the total ice sheet volume
155 - `iarea` is the total area occupied by positive thickness ice
156 - `max_diffusivity` is the maximum diffusivity
157 - `max_hor_vel` is the maximum diffusivity
158
159 Configuration parameters `output.runtime.time_unit_name`, `output.runtime.volume_scale_factor_log10`,
160 and `output.runtime.area_scale_factor_log10` control the appearance and units.
161
162 For more description and examples, see the PISM User's Manual.
163 Derived classes of IceModel may redefine this method and print alternate
164 information.
165 */
166 void IceModel::print_summary_line(bool printPrototype, bool tempAndAge,
167 double delta_t,
168 double volume, double area,
169 double /* meltfrac */, double max_diffusivity) {
170 const bool do_energy = m_config->get_flag("energy.enabled");
171 const int log10scalevol = static_cast<int>(m_config->get_number("output.runtime.volume_scale_factor_log10")),
172 log10scalearea = static_cast<int>(m_config->get_number("output.runtime.area_scale_factor_log10"));
173 const std::string time_units = m_config->get_string("output.runtime.time_unit_name");
174 const bool use_calendar = m_config->get_flag("output.runtime.time_use_calendar");
175
176 const double scalevol = pow(10.0, static_cast<double>(log10scalevol)),
177 scalearea = pow(10.0, static_cast<double>(log10scalearea));
178 char volscalestr[10] = " ", areascalestr[10] = " "; // blank when 10^0 = 1 scaling
179 if (log10scalevol != 0) {
180 snprintf(volscalestr, sizeof(volscalestr), "10^%1d_", log10scalevol);
181 }
182 if (log10scalearea != 0) {
183 snprintf(areascalestr, sizeof(areascalestr), "10^%1d_", log10scalearea);
184 }
185
186 if (printPrototype == true) {
187 m_log->message(2,
188 "P time: ivol iarea max_diffusivity max_hor_vel\n");
189 m_log->message(2,
190 "U %s %skm^3 %skm^2 m^2 s^-1 m/%s\n",
191 time_units.c_str(),volscalestr,areascalestr,time_units.c_str());
192 return;
193 }
194
195 // this version keeps track of what has been done so as to minimize stdout:
196 // FIXME: turn these static variables into class members.
197 static std::string stdout_flags_count0;
198 static int mass_cont_sub_counter = 0;
199 static double mass_cont_sub_dtsum = 0.0;
200 if (mass_cont_sub_counter == 0) {
201 stdout_flags_count0 = m_stdout_flags;
202 }
203 if (delta_t > 0.0) {
204 mass_cont_sub_counter++;
205 mass_cont_sub_dtsum += delta_t;
206 }
207
208 if ((tempAndAge == true) || (!do_energy) || (m_log->get_threshold() > 2)) {
209 char tempstr[90] = "";
210
211 const double major_dt = m_time->convert_time_interval(mass_cont_sub_dtsum, time_units);
212 if (mass_cont_sub_counter <= 1) {
213 snprintf(tempstr,90, " (dt=%.5f)", major_dt);
214 } else {
215 snprintf(tempstr,90, " (dt=%.5f in %d substeps; av dt_sub_mass_cont=%.5f)",
216 major_dt, mass_cont_sub_counter, major_dt / mass_cont_sub_counter);
217 }
218 stdout_flags_count0 += tempstr;
219
220 if (delta_t > 0.0) { // avoids printing an empty line if we have not done anything
221 stdout_flags_count0 += "\n";
222 m_log->message(2, stdout_flags_count0);
223 }
224
225 if (use_calendar) {
226 snprintf(tempstr,90, "%12s", m_time->date().c_str());
227 } else {
228 snprintf(tempstr,90, "%.3f", m_time->convert_time_interval(m_time->current(), time_units));
229 }
230
231
232
233 const CFLData cfl = m_stress_balance->max_timestep_cfl_2d();
234 std::string velocity_units = "meters / (" + time_units + ")";
235 const double maxvel = units::convert(m_sys, std::max(cfl.u_max, cfl.v_max),
236 "m second-1", velocity_units);
237
238 m_log->message(2,
239 "S %s: %8.5f %9.5f %12.5f %12.5f\n",
240 tempstr,
241 volume/(scalevol*1.0e9), area/(scalearea*1.0e6),
242 max_diffusivity, maxvel);
243
244 mass_cont_sub_counter = 0;
245 mass_cont_sub_dtsum = 0.0;
246 }
247 }
248
249
250 } // end of namespace pism