tAgeColumnSystem.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
---
tAgeColumnSystem.cc (4490B)
---
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 "AgeColumnSystem.hh"
21
22 #include "pism/util/error_handling.hh"
23
24 namespace pism {
25
26 AgeColumnSystem::AgeColumnSystem(const std::vector<double>& storage_grid,
27 const std::string &my_prefix,
28 double dx, double dy, double dt,
29 const IceModelVec3 &age,
30 const IceModelVec3 &u3,
31 const IceModelVec3 &v3,
32 const IceModelVec3 &w3)
33 : columnSystemCtx(storage_grid, my_prefix, dx, dy, dt, u3, v3, w3),
34 m_age3(age) {
35
36 size_t Mz = m_z.size();
37 m_A.resize(Mz);
38 m_A_n.resize(Mz);
39 m_A_e.resize(Mz);
40 m_A_s.resize(Mz);
41 m_A_w.resize(Mz);
42
43 m_nu = m_dt / m_dz; // derived constant
44 }
45
46 void AgeColumnSystem::init(int i, int j, double thickness) {
47 init_column(i, j, thickness);
48
49 if (m_ks == 0) {
50 return;
51 }
52
53 coarse_to_fine(m_u3, i, j, &m_u[0]);
54 coarse_to_fine(m_v3, i, j, &m_v[0]);
55 coarse_to_fine(m_w3, i, j, &m_w[0]);
56
57 coarse_to_fine(m_age3, m_i, m_j, &m_A[0]);
58 coarse_to_fine(m_age3, m_i, m_j+1, &m_A_n[0]);
59 coarse_to_fine(m_age3, m_i+1, m_j, &m_A_e[0]);
60 coarse_to_fine(m_age3, m_i, m_j-1, &m_A_s[0]);
61 coarse_to_fine(m_age3, m_i-1, m_j, &m_A_w[0]);
62 }
63
64 //! First-order upwind scheme with implicit in the vertical: one column solve.
65 /*!
66 The PDE being solved is
67 \f[ \frac{\partial \tau}{\partial t} + \frac{\partial}{\partial x}\left(u \tau\right) + \frac{\partial}{\partial y}\left(v \tau\right) + \frac{\partial}{\partial z}\left(w \tau\right) = 1. \f]
68 */
69 void AgeColumnSystem::solve(std::vector<double> &x) {
70
71 TridiagonalSystem &S = *m_solver;
72
73 // set up system: 0 <= k < m_ks
74 for (unsigned int k = 0; k < m_ks; k++) {
75 // do lowest-order upwinding, explicitly for horizontal
76 S.RHS(k) = (m_u[k] < 0 ?
77 m_u[k] * (m_A_e[k] - m_A[k]) / m_dx :
78 m_u[k] * (m_A[k] - m_A_w[k]) / m_dx);
79 S.RHS(k) += (m_v[k] < 0 ?
80 m_v[k] * (m_A_n[k] - m_A[k]) / m_dy :
81 m_v[k] * (m_A[k] - m_A_s[k]) / m_dy);
82 // note it is the age eqn: dage/dt = 1.0 and we have moved the hor.
83 // advection terms over to right:
84 S.RHS(k) = m_A[k] + m_dt * (1.0 - S.RHS(k));
85
86 // do lowest-order upwinding, *implicitly* for vertical
87 double AA = m_nu * m_w[k];
88 if (k > 0) {
89 if (AA >= 0) { // upward velocity
90 S.L(k) = - AA;
91 S.D(k) = 1.0 + AA;
92 S.U(k) = 0.0;
93 } else { // downward velocity; note -AA >= 0
94 S.L(k) = 0.0;
95 S.D(k) = 1.0 - AA;
96 S.U(k) = + AA;
97 }
98 } else { // k == 0 case
99 // note L[0] is not used
100 if (AA > 0) { // if strictly upward velocity apply boundary condition:
101 // age = 0 because ice is being added to base
102 S.D(0) = 1.0;
103 S.U(0) = 0.0;
104 S.RHS(0) = 0.0;
105 } else { // downward velocity; note -AA >= 0
106 S.D(0) = 1.0 - AA;
107 S.U(0) = + AA;
108 // keep rhs[0] as is
109 }
110 }
111 } // done "set up system: 0 <= k < m_ks"
112
113 // surface b.c. at m_ks
114 if (m_ks > 0) {
115 S.L(m_ks) = 0;
116 S.D(m_ks) = 1.0; // ignore U[m_ks]
117 S.RHS(m_ks) = 0.0; // age zero at surface
118 }
119
120 // solve it
121 try {
122 S.solve(m_ks + 1, x);
123 }
124 catch (RuntimeError &e) {
125 e.add_context("solving the tri-diagonal system (AgeColumnSystem) at (%d, %d)\n"
126 "saving system to m-file... ", m_i, m_j);
127 reportColumnZeroPivotErrorMFile(m_ks + 1);
128 throw;
129 }
130
131 // x[k] contains age for k=0,...,ks, but set age of ice above (and
132 // at) surface to zero years
133 for (unsigned int k = m_ks + 1; k < x.size(); k++) {
134 x[k] = 0.0;
135 }
136 }
137
138 } // end of namespace pism