tenthSystem.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
---
tenthSystem.cc (20110B)
---
1 // Copyright (C) 2009-2018 Andreas Aschwanden and 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 "enthSystem.hh"
20 #include <gsl/gsl_math.h> // GSL_NAN, gsl_isnan()
21 #include "pism/util/ConfigInterface.hh"
22 #include "pism/util/iceModelVec.hh"
23 #include "pism/util/EnthalpyConverter.hh"
24
25 #include "pism/util/error_handling.hh"
26 #include "pism/util/ColumnInterpolation.hh"
27
28 namespace pism {
29 namespace energy {
30
31 enthSystemCtx::enthSystemCtx(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 &Enth3,
36 const IceModelVec3 &u3,
37 const IceModelVec3 &v3,
38 const IceModelVec3 &w3,
39 const IceModelVec3 &strain_heating3,
40 EnthalpyConverter::Ptr EC)
41 : columnSystemCtx(storage_grid, prefix, dx, dy, dt, u3, v3, w3),
42 m_Enth3(Enth3),
43 m_strain_heating3(strain_heating3),
44 m_EC(EC) {
45
46 // set some values so we can check if init was called
47 m_R_cold = -1.0;
48 m_R_temp = -1.0;
49 m_lambda = -1.0;
50 m_D0 = GSL_NAN;
51 m_U0 = GSL_NAN;
52 m_B0 = GSL_NAN;
53 m_L_ks = GSL_NAN;
54 m_D_ks = GSL_NAN;
55 m_U_ks = GSL_NAN;
56 m_B_ks = GSL_NAN;
57
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 m_p_air = config.get_number("surface.pressure");
62
63 m_exclude_horizontal_advection = config.get_flag("energy.margin_exclude_horizontal_advection");
64 m_exclude_vertical_advection = config.get_flag("energy.margin_exclude_vertical_advection");
65 m_exclude_strain_heat = config.get_flag("energy.margin_exclude_strain_heating");
66
67 size_t Mz = m_z.size();
68 m_Enth.resize(Mz);
69 m_Enth_s.resize(Mz);
70 m_strain_heating.resize(Mz);
71 m_R.resize(Mz);
72
73 m_E_n.resize(Mz);
74 m_E_e.resize(Mz);
75 m_E_s.resize(Mz);
76 m_E_w.resize(Mz);
77
78 m_nu = m_dt / m_dz;
79
80 double
81 ratio = config.get_number(prefix + ".temperate_ice_thermal_conductivity_ratio"),
82 K = m_ice_k / m_ice_c,
83 K0 = (ratio * m_ice_k) / m_ice_c;
84
85 m_R_factor = m_dt / (PetscSqr(m_dz) * m_ice_density);
86 m_R_cold = K * m_R_factor;
87 m_R_temp = K0 * m_R_factor;
88
89 if (config.get_flag("energy.temperature_dependent_thermal_conductivity")) {
90 m_k_depends_on_T = true;
91 } else {
92 m_k_depends_on_T = false;
93 }
94 }
95
96
97 enthSystemCtx::~enthSystemCtx() {
98 }
99
100 /*!
101 In this implementation \f$k\f$ does not depend on temperature.
102 */
103 double enthSystemCtx::k_from_T(double T) const {
104
105 if (m_k_depends_on_T) {
106 return 9.828 * exp(-0.0057 * T);
107 }
108
109 return m_ice_k;
110 }
111
112 void enthSystemCtx::init(int i, int j, bool marginal, double ice_thickness) {
113 m_ice_thickness = ice_thickness;
114
115 m_marginal = marginal;
116
117 init_column(i, j, m_ice_thickness);
118
119 if (m_ks == 0) {
120 return;
121 }
122
123 coarse_to_fine(m_u3, m_i, m_j, &m_u[0]);
124 coarse_to_fine(m_v3, m_i, m_j, &m_v[0]);
125
126 if (m_marginal and m_exclude_vertical_advection) {
127 for (unsigned int k = 0; k < m_w.size(); ++k) {
128 m_w[k] = 0.0;
129 }
130 } else {
131 coarse_to_fine(m_w3, m_i, m_j, &m_w[0]);
132 }
133
134 coarse_to_fine(m_strain_heating3, m_i, m_j, &m_strain_heating[0]);
135 coarse_to_fine(m_Enth3, m_i, m_j, &m_Enth[0]);
136
137 coarse_to_fine(m_Enth3, m_i, m_j+1, &m_E_n[0]);
138 coarse_to_fine(m_Enth3, m_i+1, m_j, &m_E_e[0]);
139 coarse_to_fine(m_Enth3, m_i, m_j-1, &m_E_s[0]);
140 coarse_to_fine(m_Enth3, m_i-1, m_j, &m_E_w[0]);
141
142 compute_enthalpy_CTS();
143
144 m_lambda = compute_lambda();
145
146 assemble_R();
147 }
148
149 //! Compute the CTS value of enthalpy in an ice column.
150 /*! m_Enth_s is set to the enthalpy value for the pressure-melting
151 temperature with zero water fraction at the corresponding z level.
152 */
153 void enthSystemCtx::compute_enthalpy_CTS() {
154
155 for (unsigned int k = 0; k <= m_ks; k++) {
156 const double
157 depth = m_ice_thickness - k * m_dz,
158 p = m_EC->pressure(depth); // FIXME issue #15
159 m_Enth_s[k] = m_EC->enthalpy_cts(p);
160 }
161
162 const double Es_air = m_EC->enthalpy_cts(m_p_air);
163 for (unsigned int k = m_ks+1; k < m_Enth_s.size(); k++) {
164 m_Enth_s[k] = Es_air;
165 }
166 }
167
168 //! Compute the lambda for BOMBPROOF.
169 /*!
170 See page \ref bombproofenth.
171 */
172 double enthSystemCtx::compute_lambda() {
173 double result = 1.0; // start with centered implicit for more accuracy
174 const double epsilon = 1e-6 / 3.15569259747e7;
175
176 for (unsigned int k = 0; k <= m_ks; k++) {
177 if (m_Enth[k] > m_Enth_s[k]) { // lambda = 0 if temperate ice present in column
178 result = 0.0;
179 } else {
180 const double denom = (fabs(m_w[k]) + epsilon) * m_ice_density * m_ice_c * m_dz;
181 result = std::min(result, 2.0 * m_ice_k / denom);
182 }
183 }
184 return result;
185 }
186
187
188 void enthSystemCtx::set_surface_dirichlet_bc(double E_surface) {
189 #if (Pism_DEBUG==1)
190 if ((m_nu < 0.0) || (m_R_cold < 0.0) || (m_R_temp < 0.0)) {
191 throw RuntimeError(PISM_ERROR_LOCATION, "setDirichletSurface() should only be called after\n"
192 "initAllColumns() in enthSystemCtx");
193 }
194 #endif
195 m_L_ks = 0.0;
196 m_D_ks = 1.0;
197 m_U_ks = 0.0;
198 m_B_ks = E_surface;
199 }
200
201 static inline double upwind(double u, double E_m, double E, double E_p, double delta_inverse) {
202 return u * delta_inverse * (u < 0 ? (E_p - E) : (E - E_m));
203 }
204
205
206 //! Set the top surface heat flux *into* the ice.
207 /** @param[in] heat_flux prescribed heat flux (positive means flux into the ice)
208 */
209 void enthSystemCtx::set_surface_heat_flux(double heat_flux) {
210 // extract K from R[ks], so this code works even if K=K(T)
211 // recall: R = (ice_K / ice_density) * dt / PetscSqr(dz)
212 const double
213 K = (m_ice_density * PetscSqr(m_dz) * m_R[m_ks]) / m_dt,
214 G = heat_flux / K;
215
216 this->set_surface_neumann_bc(G);
217 }
218
219 //! Set enthalpy flux at the surface.
220 /*! This method should probably be used for debugging only. Its purpose is to allow setting the
221 enthalpy flux even if K == 0, i.e. in a "pure advection" setup.
222 */
223 void enthSystemCtx::set_surface_neumann_bc(double G) {
224 const bool include_horizontal_advection = not (m_marginal and m_exclude_horizontal_advection);
225 const bool include_strain_heating = not (m_marginal and m_exclude_strain_heat);
226
227 const double
228 Rminus = 0.5 * (m_R[m_ks - 1] + m_R[m_ks]), // R_{ks-1/2}
229 Rplus = m_R[m_ks], // R_{ks+1/2}
230 mu_w = 0.5 * m_nu * m_w[m_ks];
231
232 const double A_l = m_w[m_ks] < 0.0 ? 1.0 - m_lambda : m_lambda - 1.0;
233 const double A_d = m_w[m_ks] < 0.0 ? m_lambda - 1.0 : 1.0 - m_lambda;
234 const double A_b = m_w[m_ks] < 0.0 ? m_lambda - 2.0 : -m_lambda;
235
236 // modified lower-diagonal entry:
237 m_L_ks = - Rminus - Rplus + 2.0 * mu_w * A_l;
238 // diagonal entry
239 m_D_ks = 1.0 + Rminus + Rplus + 2.0 * mu_w * A_d;
240 // upper-diagonal entry (not used)
241 m_U_ks = 0.0;
242 // m_Enth[m_ks] (below) is there due to the fully-implicit discretization in time, the second term is
243 // the modification of the right-hand side implementing the Neumann B.C. (similar to
244 // set_basal_heat_flux(); see that method for details)
245 m_B_ks = m_Enth[m_ks] + 2.0 * G * m_dz * (Rplus + mu_w * A_b);
246
247 // treat horizontal velocity using first-order upwinding:
248 double upwind_u = 0.0;
249 double upwind_v = 0.0;
250 if (include_horizontal_advection) {
251 upwind_u = upwind(m_u[m_ks], m_E_w[m_ks], m_Enth[m_ks], m_E_e[m_ks], 1.0 / m_dx);
252 upwind_v = upwind(m_v[m_ks], m_E_s[m_ks], m_Enth[m_ks], m_E_n[m_ks], 1.0 / m_dy);
253 }
254 double Sigma = 0.0;
255 if (include_strain_heating) {
256 Sigma = m_strain_heating[m_ks];
257 }
258
259 m_B_ks += m_dt * ((Sigma / m_ice_density) - upwind_u - upwind_v); // = rhs[m_ks]
260 }
261
262 void enthSystemCtx::checkReadyToSolve() {
263 if (m_nu < 0.0 || m_R_cold < 0.0 || m_R_temp < 0.0) {
264 throw RuntimeError(PISM_ERROR_LOCATION,
265 "not ready to solve: need initAllColumns() in enthSystemCtx");
266 }
267 if (m_lambda < 0.0) {
268 throw RuntimeError(PISM_ERROR_LOCATION,
269 "not ready to solve: need setSchemeParamsThisColumn() in enthSystemCtx");
270 }
271 }
272
273
274 //! Set coefficients in discrete equation for \f$E = Y\f$ at base of ice.
275 /*!
276 This method should only be called if everything but the basal boundary condition
277 is already set.
278 */
279 void enthSystemCtx::set_basal_dirichlet_bc(double Y) {
280 #if (Pism_DEBUG==1)
281 checkReadyToSolve();
282 if (gsl_isnan(m_D0) == 0 || gsl_isnan(m_U0) == 0 || gsl_isnan(m_B0) == 0) {
283 throw RuntimeError(PISM_ERROR_LOCATION, "setting basal boundary conditions twice in enthSystemCtx");
284 }
285 #endif
286 m_D0 = 1.0;
287 m_U0 = 0.0;
288 m_B0 = Y;
289 }
290
291 //! Set coefficients in discrete equation for Neumann condition at base of ice.
292 /*!
293 This method generates the Neumann boundary condition for the linear system.
294
295 The Neumann boundary condition is
296 @f[ \frac{\partial E}{\partial z} = - \frac{\phi}{K} @f]
297 where \f$\phi\f$ is the heat flux. Here \f$K\f$ is allowed to vary, and takes
298 its value from the value computed in assemble_R().
299
300 The boundary condition is combined with the partial differential equation by the
301 technique of introducing an imaginary point at \f$z=-\Delta z\f$ and then
302 eliminating it.
303
304 In other words, we combine the centered finite difference approximation
305 @f[ \frac{ E_{1} - E_{-1} }{2\dz} = -\frac{\phi}{K} @f]
306 with
307
308 @f[ -R_{k-\frac12} E_{k-1} + \left( 1 + R_{k-\frac12} + R_{k+\frac12} \right) E_{k} - R_{k+\frac12} E_{k+1} + \text{advective terms} = \dots @f]
309
310 to get
311
312 @f{align*}{
313 \frac{E_{1}-E_{-1}}{2\,\Delta z} & = -\frac{\phi}{K_{0}}, \\
314 E_{1}-E_{-1} & = -\frac{2\,\Delta z\,\phi}{K_{0}}, \\
315 E_{-1}\,R_{-\frac12}-R_{-\frac12}\,E_{1} & = \frac{2\,R_{-\frac12}\,\Delta z\,\phi}{K_{0}}, \\
316 -R_{\frac12}\,E_{1}+E_{0}\,\left(R_{\frac12}+R_{-\frac12}+1\right)-E_{-1}\,R_{-\frac12} + \text{advective terms} & = \dots, \\
317 \left(-R_{\frac12}-R_{-\frac12}\right)\,E_{1}+E_{0}\,\left(R_{\frac12}+R_{-\frac12}+1\right) + \text{advective terms} & = \frac{2\,R_{-\frac12}\,\Delta z\,\phi}{K_{0}}+\dots.
318 @f}
319
320 The error in the pure conductive and smooth conductivity case is @f$ O(\dz^2) @f$.
321
322 This method should only be called if everything but the basal boundary condition
323 is already set.
324
325 @param[in] heat_flux prescribed heat flux (positive means flux into the ice)
326
327 */
328 void enthSystemCtx::set_basal_heat_flux(double heat_flux) {
329 // extract K from R[0], so this code works even if K=K(T)
330 // recall: R = (ice_K / ice_density) * dt / PetscSqr(dz)
331 const double
332 K = (m_ice_density * PetscSqr(m_dz) * m_R[0]) / m_dt,
333 G = - heat_flux / K;
334
335 this->set_basal_neumann_bc(G);
336 }
337
338 //! Set enthalpy flux at the base.
339 /*! This method should probably be used for debugging only. Its purpose is to allow setting the
340 enthalpy flux even if K == 0, i.e. in a "pure advection" setup.
341 */
342 void enthSystemCtx::set_basal_neumann_bc(double G) {
343 const bool include_horizontal_advection = not (m_marginal and m_exclude_horizontal_advection);
344 const bool include_strain_heating = not (m_marginal and m_exclude_strain_heat);
345
346 const double
347 Rminus = m_R[0], // R_{-1/2}
348 Rplus = 0.5 * (m_R[0] + m_R[1]), // R_{+1/2}
349 mu_w = 0.5 * m_nu * m_w[0];
350
351 const double A_d = m_w[0] < 0.0 ? m_lambda - 1.0 : 1.0 - m_lambda;
352 const double A_u = m_w[0] < 0.0 ? 1.0 - m_lambda : m_lambda - 1.0;
353 const double A_b = m_w[0] < 0.0 ? -m_lambda : m_lambda - 2.0;
354
355 // diagonal entry
356 m_D0 = 1.0 + Rminus + Rplus + 2.0 * mu_w * A_d;
357 // upper-diagonal entry
358 m_U0 = - Rminus - Rplus + 2.0 * mu_w * A_u;
359 // right-hand side, excluding the strain heating term and the horizontal advection
360 m_B0 = m_Enth[0] + 2.0 * G * m_dz * (-Rminus + mu_w * A_b);
361
362 // treat horizontal velocity using first-order upwinding:
363 double upwind_u = 0.0;
364 double upwind_v = 0.0;
365 if (include_horizontal_advection) {
366 upwind_u = upwind(m_u[0], m_E_w[0], m_Enth[0], m_E_e[0], 1.0 / m_dx);
367 upwind_v = upwind(m_v[0], m_E_s[0], m_Enth[0], m_E_n[0], 1.0 / m_dy);
368 }
369 double Sigma = 0.0;
370 if (include_strain_heating) {
371 Sigma = m_strain_heating[0];
372 }
373
374 m_B0 += m_dt * ((Sigma / m_ice_density) - upwind_u - upwind_v); // = rhs[m_ks]
375 }
376
377
378 //! \brief Assemble the R array. The R value switches at the CTS.
379 /*! In a simple abstract diffusion
380 \f[ \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial z^2}, \f]
381 with time steps \f$\Delta t\f$ and spatial steps \f$\Delta z\f$ we define
382 \f[ R = \frac{D \Delta t}{\Delta z^2}. \f]
383 This is used in an implicit method to write each line in the linear system, for
384 example [\ref MortonMayers]:
385 \f[ -R U_{j-1}^{n+1} + (1+2R) U_j^{n+1} - R U_{j+1}^{n+1} = U_j^n. \f]
386
387 In the case of conservation of energy [\ref AschwandenBuelerKhroulevBlatter],
388 \f[ u=E \qquad \text{ and } \qquad D = \frac{K}{\rho} \qquad \text{ and } \qquad K = \frac{k}{c}. \f]
389 Thus
390 \f[ R = \frac{k \Delta t}{\rho c \Delta z^2}. \f]
391 */
392 void enthSystemCtx::assemble_R() {
393 if (not m_k_depends_on_T) {
394
395 for (unsigned int k = 1; k <= m_ks; k++) {
396 m_R[k] = (m_Enth[k] < m_Enth_s[k]) ? m_R_cold : m_R_temp;
397 }
398 //still the cold ice value, if no temperate layer above
399 m_R[0] = (m_Enth[1] < m_Enth_s[1]) ? m_R_cold : m_R_temp;
400
401 } else {
402
403 for (unsigned int k = 1; k <= m_ks; k++) {
404 if (m_Enth[k] < m_Enth_s[k]) {
405 // cold case
406 const double depth = m_ice_thickness - k * m_dz;
407 double T = m_EC->temperature(m_Enth[k],
408 m_EC->pressure(depth)); // FIXME: issue #15
409
410 m_R[k] = ((m_k_depends_on_T ? k_from_T(T) : m_ice_k) / m_EC->c()) * m_R_factor;
411 } else {
412 // temperate case
413 m_R[k] = m_R_temp;
414 }
415 }
416 // still the cold ice value, if no temperate layer above
417 if (m_Enth[1] < m_Enth_s[1]) {
418 double T = m_EC->temperature(m_Enth[0],
419 m_EC->pressure(m_ice_thickness)); // FIXME: issue #15
420 m_R[0] = ((m_k_depends_on_T ? k_from_T(T) : m_ice_k) / m_EC->c()) * m_R_factor;
421 } else {
422 // temperate layer case
423 m_R[0] = m_R_temp;
424 }
425
426 }
427
428 // R[k] for k > m_ks are never used
429 #if (Pism_DEBUG==1)
430 for (unsigned int k = m_ks + 1; k < m_R.size(); ++k) {
431 m_R[k] = GSL_NAN;
432 }
433 #endif
434 }
435
436
437 /*! \brief Solve the tridiagonal system, in a single column, which
438 * determines the new values of the ice enthalpy.
439 *
440 * We are solving a convection-diffusion equation, treating the @f$ z @f$ direction implicitly and
441 * @f$ x, y @f$ directions explicitly. See @ref bombproofenth for the documentation of the treatment
442 * of convection terms. The notes below document the way we treat diffusion using a simplified PDE.
443 *
444 * To discretize
445 * @f[ \diff{}{z} \left( K(E) \diff{E}{z}\right) = \diff{E}{t} @f]
446 *
447 * at a location @f$ k @f$ of the vertical grid, we use centered finite differences in space,
448 * backward differences in time, and evaluate @f$ K(E) @f$ at staggered-grid locations:
449 *
450 * @f[ \frac{K_{k+\frac12}\frac{E_{k+1} - E_{k}}{\dz} - K_{k-\frac12}\frac{E_{k} - E_{k-1}}{\dz}}{\dz} = \frac{E_{k} - E^{n-1}_{k}}{\dt}, @f]
451 *
452 * where @f$ E = E^{n} @f$ for clarity and @f$ K_{k\pm \frac12} = K(E^{n-1}_{k\pm \frac12}) @f$,
453 * %i.e. we linearize this PDE by evaluating @f$ K(E) @f$ at the _previous_ time-step.
454 *
455 * We define
456 *
457 * @f[ R_i = \frac{\dt\, K_i}{\dz^2}, @f]
458 *
459 * and the discretization takes form
460 *
461 * @f[ -R_{k-\frac12} E_{k-1} + \left( 1 + R_{k-\frac12} + R_{k+\frac12} \right) E_{k} - R_{k+\frac12} E_{k+1} = E^{n-1}_{k}. @f]
462 *
463 * In the assembly of the tridiagonal system this corresponds to
464 *
465 * @f{align*}{
466 * L_i &= - \frac12 (R_{i} + R_{i-1}),\\
467 * D_i &= 1 + \frac12 (R_{i} + R_{i-1}) + \frac12 (R_{i} + R_{i+1}),\\
468 * U_i &= - \frac12 (R_{i} + R_{i+1}),\\
469 * b_i &= E^{n-1}_{i},
470 * @f}
471 *
472 * where @f$ L_i, D_i, U_i @f$ are lower-diagonal, diagonal, and upper-diagonal entries
473 * corresponding to an equation @f$ i @f$ and @f$ b_i @f$ is the corresponding right-hand side.
474 * (Staggered-grid values are approximated by interpolating from the regular grid).
475 *
476 * This method is _unconditionally stable_ and has a maximum principle (see [@ref MortonMayers,
477 * section 2.11]).
478 */
479 void enthSystemCtx::solve(std::vector<double> &x) {
480
481 TridiagonalSystem &S = *m_solver;
482
483 #if (Pism_DEBUG==1)
484 checkReadyToSolve();
485 if (gsl_isnan(m_D0) || gsl_isnan(m_U0) || gsl_isnan(m_B0)) {
486 throw RuntimeError(PISM_ERROR_LOCATION,
487 "solveThisColumn() should only be called after\n"
488 " setting basal boundary condition in enthSystemCtx");
489 }
490 #endif
491
492 // k=0 equation is already established
493 // L[0] = 0.0; // not used
494 S.D(0) = m_D0;
495 S.U(0) = m_U0;
496 S.RHS(0) = m_B0;
497
498 const double
499 one_over_rho = 1.0 / m_ice_density,
500 Dx = 1.0 / m_dx,
501 Dy = 1.0 / m_dy;
502
503 const bool include_horizontal_advection = not (m_marginal and m_exclude_horizontal_advection);
504 const bool include_strain_heating = not (m_marginal and m_exclude_strain_heat);
505
506 // generic ice segment in k location (if any; only runs if m_ks >= 2)
507 for (unsigned int k = 1; k < m_ks; k++) {
508 const double
509 Rminus = 0.5 * (m_R[k-1] + m_R[k]), // R_{k-1/2}
510 Rplus = 0.5 * (m_R[k] + m_R[k+1]), // R_{k+1/2}
511 nu_w = m_nu * m_w[k];
512
513 const double
514 A_l = m_w[k] >= 0.0 ? 0.5 * m_lambda - 1.0 : -0.5 * m_lambda,
515 A_d = m_w[k] >= 0.0 ? 1.0 - m_lambda : m_lambda - 1.0,
516 A_u = m_w[k] >= 0.0 ? 0.5 * m_lambda : 1.0 - 0.5 * m_lambda;
517
518 S.L(k) = - Rminus + nu_w * A_l;
519 S.D(k) = 1.0 + Rminus + Rplus + nu_w * A_d;
520 S.U(k) = - Rplus + nu_w * A_u;
521
522 // horizontal velocity and strain heating
523 double upwind_u = 0.0;
524 double upwind_v = 0.0;
525 if (include_horizontal_advection) {
526 upwind_u = upwind(m_u[k], m_E_w[k], m_Enth[k], m_E_e[k], Dx);
527 upwind_v = upwind(m_v[k], m_E_s[k], m_Enth[k], m_E_n[k], Dy);
528 }
529 double Sigma = 0.0;
530 if (include_strain_heating) {
531 Sigma = m_strain_heating[k];
532 }
533
534 S.RHS(k) = m_Enth[k] + m_dt * (one_over_rho * Sigma - upwind_u - upwind_v);
535 }
536
537 // Assemble the top surface equation. Values m_{L,D,U,B}_ks are set using set_surface_dirichlet()
538 // or set_surface_heat_flux().
539 if (m_ks > 0) {
540 S.L(m_ks) = m_L_ks;
541 }
542 S.D(m_ks) = m_D_ks;
543 if (m_ks < m_z.size() - 1) {
544 S.U(m_ks) = m_U_ks;
545 }
546 S.RHS(m_ks) = m_B_ks;
547
548 // Solve it; note drainage is not addressed yet and post-processing may occur
549 try {
550 S.solve(m_ks + 1, x);
551 }
552 catch (RuntimeError &e) {
553 e.add_context("solving the tri-diagonal system (enthSystemCtx) at (%d,%d)\n"
554 "saving system to m-file... ", m_i, m_j);
555 reportColumnZeroPivotErrorMFile(m_ks + 1);
556 throw;
557 }
558
559 // air above
560 for (unsigned int k = m_ks+1; k < x.size(); k++) {
561 x[k] = m_B_ks;
562 }
563
564 #if (Pism_DEBUG==1)
565 // if success, mark column as done by making scheme params and b.c. coeffs invalid
566 m_lambda = -1.0;
567 m_D0 = GSL_NAN;
568 m_U0 = GSL_NAN;
569 m_B0 = GSL_NAN;
570 m_L_ks = GSL_NAN;
571 m_D_ks = GSL_NAN;
572 m_U_ks = GSL_NAN;
573 m_B_ks = GSL_NAN;
574 #endif
575 }
576
577 void enthSystemCtx::save_system(std::ostream &output, unsigned int system_size) const {
578 m_solver->save_system(output, system_size);
579 m_solver->save_vector(output, m_R, system_size, m_solver->prefix() + "_R");
580 }
581
582 } // end of namespace energy
583 } // end of namespace pism