ttempSystem.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
---
ttempSystem.cc (7744B)
---
1 // Copyright (C) 2004-2011, 2013, 2014, 2015, 2016, 2017, 2018 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 <cassert>
20
21 #include "pism/util/pism_utilities.hh"
22 #include "pism/util/iceModelVec.hh"
23 #include "tempSystem.hh"
24 #include "pism/util/Mask.hh"
25
26 #include "pism/util/error_handling.hh"
27
28 namespace pism {
29 namespace energy {
30
31 tempSystemCtx::tempSystemCtx(const std::vector<double>& storage_grid,
32 const std::string &prefix,
33 double dx, double dy, double dt,
34 const Config &config,
35 const IceModelVec3 &T3,
36 const IceModelVec3 &u3,
37 const IceModelVec3 &v3,
38 const IceModelVec3 &w3,
39 const IceModelVec3 &strain_heating3)
40 : columnSystemCtx(storage_grid, prefix, dx, dy, dt, u3, v3, w3),
41 m_T3(T3),
42 m_strain_heating3(strain_heating3) {
43
44 // set flags to indicate nothing yet set
45 m_surfBCsValid = false;
46 m_basalBCsValid = false;
47
48 size_t Mz = m_z.size();
49 m_T.resize(Mz);
50 m_strain_heating.resize(Mz);
51
52 m_T_n.resize(Mz);
53 m_T_e.resize(Mz);
54 m_T_s.resize(Mz);
55 m_T_w.resize(Mz);
56
57 // set physical constants
58 m_ice_density = config.get_number("constants.ice.density");
59 m_ice_c = config.get_number("constants.ice.specific_heat_capacity");
60 m_ice_k = config.get_number("constants.ice.thermal_conductivity");
61
62 // set derived constants
63 m_nu = m_dt / m_dz;
64 m_rho_c_I = m_ice_density * m_ice_c;
65 m_iceK = m_ice_k / m_rho_c_I;
66 m_iceR = m_iceK * m_dt / (m_dz*m_dz);
67 }
68
69 void tempSystemCtx::initThisColumn(int i, int j, bool is_marginal, MaskValue mask,
70 double ice_thickness) {
71
72 m_is_marginal = is_marginal;
73 m_mask = mask;
74
75 init_column(i, j, ice_thickness);
76
77 if (m_ks == 0) {
78 return;
79 }
80
81 coarse_to_fine(m_u3, m_i, m_j, &m_u[0]);
82 coarse_to_fine(m_v3, m_i, m_j, &m_v[0]);
83 coarse_to_fine(m_w3, m_i, m_j, &m_w[0]);
84 coarse_to_fine(m_strain_heating3, m_i, m_j, &m_strain_heating[0]);
85 coarse_to_fine(m_T3, m_i, m_j, &m_T[0]);
86
87 coarse_to_fine(m_T3, m_i, m_j+1, &m_T_n[0]);
88 coarse_to_fine(m_T3, m_i+1, m_j, &m_T_e[0]);
89 coarse_to_fine(m_T3, m_i, m_j-1, &m_T_s[0]);
90 coarse_to_fine(m_T3, m_i-1, m_j, &m_T_w[0]);
91
92 m_lambda = compute_lambda();
93 }
94
95 void tempSystemCtx::setSurfaceBoundaryValuesThisColumn(double my_Ts) {
96 // allow setting surface BCs only once:
97 assert(not m_surfBCsValid);
98
99 m_Ts = my_Ts;
100 m_surfBCsValid = true;
101 }
102
103
104 void tempSystemCtx::setBasalBoundaryValuesThisColumn(double my_G0,
105 double my_Tshelfbase, double my_Rb) {
106 // allow setting basal BCs only once:
107 assert(not m_basalBCsValid);
108
109 m_G0 = my_G0;
110 m_Tshelfbase = my_Tshelfbase;
111 m_Rb = my_Rb;
112 m_basalBCsValid = true;
113 }
114
115 double tempSystemCtx::compute_lambda() {
116 double result = 1.0; // start with centered implicit for more accuracy
117 const double epsilon = 1e-6 / 3.15569259747e7;
118
119 for (unsigned int k = 0; k <= m_ks; k++) {
120 const double denom = (fabs(m_w[k]) + epsilon) * m_ice_density * m_ice_c * m_dz;
121 result = std::min(result, 2.0 * m_ice_k / denom);
122 }
123 return result;
124 }
125
126 void tempSystemCtx::solveThisColumn(std::vector<double> &x) {
127
128 TridiagonalSystem &S = *m_solver;
129
130 assert(m_surfBCsValid == true);
131 assert(m_basalBCsValid == true);
132
133 // bottom of ice; k=0 eqn
134 if (m_ks == 0) { // no ice; set m_T[0] to surface temp if grounded
135 // note L[0] not allocated
136 S.D(0) = 1.0;
137 S.U(0) = 0.0;
138 // if floating and no ice then worry only about bedrock temps
139 if (mask::ocean(m_mask)) {
140 // essentially no ice but floating ... ask OceanCoupler
141 S.RHS(0) = m_Tshelfbase;
142 } else { // top of bedrock sees atmosphere
143 S.RHS(0) = m_Ts;
144 }
145 } else { // m_ks > 0; there is ice
146 // for w, always difference *up* from base, but make it implicit
147 if (mask::ocean(m_mask)) {
148 // just apply Dirichlet condition to base of column of ice in an ice shelf
149 // note that L[0] is not used
150 S.D(0) = 1.0;
151 S.U(0) = 0.0;
152 S.RHS(0) = m_Tshelfbase; // set by OceanCoupler
153 } else {
154 // there is *grounded* ice; from FV across interface
155 S.RHS(0) = m_T[0] + m_dt * (m_Rb / (m_rho_c_I * m_dz));
156
157 double Sigma = 0.0;
158 if (not m_is_marginal) {
159 Sigma = m_strain_heating[0];
160 }
161 S.RHS(0) += m_dt * 0.5 * Sigma / m_rho_c_I;
162
163 double UpTu = 0.0;
164 double UpTv = 0.0;
165 if (not m_is_marginal) {
166 UpTu = (m_u[0] < 0 ?
167 m_u[0] * (m_T_e[0] - m_T[0]) / m_dx :
168 m_u[0] * (m_T[0] - m_T_w[0]) / m_dx);
169 UpTv = (m_v[0] < 0 ?
170 m_v[0] * (m_T_n[0] - m_T[0]) / m_dy :
171 m_v[0] * (m_T[0] - m_T_s[0]) / m_dy);
172 }
173 S.RHS(0) -= m_dt * (0.5 * (UpTu + UpTv));
174
175 // vertical upwinding
176 // L[0] = 0.0; (is not used)
177 S.D(0) = 1.0 + 2.0 * m_iceR;
178 S.U(0) = - 2.0 * m_iceR;
179 if (m_w[0] < 0.0) { // velocity downward: add velocity contribution
180 const double AA = m_dt * m_w[0] / (2.0 * m_dz);
181 S.D(0) -= AA;
182 S.U(0) += AA;
183 }
184 // apply geothermal flux G0 here
185 S.RHS(0) += 2.0 * m_dt * m_G0 / (m_rho_c_I * m_dz);
186 }
187 }
188
189 // generic ice segment; build 1:m_ks-1 eqns
190 for (unsigned int k = 1; k < m_ks; k++) {
191 const double AA = m_nu * m_w[k];
192 if (m_w[k] >= 0.0) { // velocity upward
193 S.L(k) = - m_iceR - AA * (1.0 - m_lambda/2.0);
194 S.D(k) = 1.0 + 2.0 * m_iceR + AA * (1.0 - m_lambda);
195 S.U(k) = - m_iceR + AA * (m_lambda/2.0);
196 } else { // velocity downward
197 S.L(k) = - m_iceR - AA * (m_lambda/2.0);
198 S.D(k) = 1.0 + 2.0 * m_iceR - AA * (1.0 - m_lambda);
199 S.U(k) = - m_iceR + AA * (1.0 - m_lambda/2.0);
200 }
201 S.RHS(k) = m_T[k];
202
203 double Sigma = 0.0;
204 if (not m_is_marginal) {
205 Sigma = m_strain_heating[k];
206 }
207
208 double UpTu = 0.0;
209 double UpTv = 0.0;
210 if (not m_is_marginal) {
211 UpTu = (m_u[k] < 0 ?
212 m_u[k] * (m_T_e[k] - m_T[k]) / m_dx :
213 m_u[k] * (m_T[k] - m_T_w[k]) / m_dx);
214 UpTv = (m_v[k] < 0 ?
215 m_v[k] * (m_T_n[k] - m_T[k]) / m_dy :
216 m_v[k] * (m_T[k] - m_T_s[k]) / m_dy);
217 }
218
219 S.RHS(k) += m_dt * (Sigma / m_rho_c_I - UpTu - UpTv);
220 }
221
222 // surface b.c.
223 if (m_ks>0) {
224 S.L(m_ks) = 0.0;
225 S.D(m_ks) = 1.0;
226 // ignore U[m_ks]
227 S.RHS(m_ks) = m_Ts;
228 }
229
230 // mark column as done
231 m_surfBCsValid = false;
232 m_basalBCsValid = false;
233
234 // solve it; note melting not addressed yet
235 try {
236 S.solve(m_ks + 1, x);
237 }
238 catch (RuntimeError &e) {
239 e.add_context("solving the tri-diagonal system (tempSystemCtx) at (%d,%d)\n"
240 "saving system to m-file... ", m_i, m_j);
241 reportColumnZeroPivotErrorMFile(m_ks + 1);
242 throw;
243 }
244 }
245
246
247 } // end of namespace energy
248 } // end of namespace pism