PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
AgeColumnSystem.cc
Go to the documentation of this file.
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 
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 
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);
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
std::vector< double > m_A_s
std::vector< double > m_A_n
std::vector< double > m_A_w
std::vector< double > m_A_e
std::vector< double > m_A
const IceModelVec3 & m_age3
void init(int i, int j, double thickness)
AgeColumnSystem(const std::vector< double > &storage_grid, const std::string &my_prefix, double dx, double dy, double dt, const IceModelVec3 &age, const IceModelVec3 &u3, const IceModelVec3 &v3, const IceModelVec3 &w3)
void solve(std::vector< double > &x)
First-order upwind scheme with implicit in the vertical: one column solve.
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition: iceModelVec.hh:404
void add_context(const std::string &message)
Add a message providing some context. This way we can (sort of) get a stack trace even though C++ exc...
Virtual base class. Abstracts a tridiagonal system to solve in a column of ice and/or bedrock.
Definition: ColumnSystem.hh:82
const IceModelVec3 & m_u3
pointers to 3D velocity components
const IceModelVec3 & m_w3
void reportColumnZeroPivotErrorMFile(unsigned int M)
Write system matrix and right-hand-side into an Python script. The file name contains ZERO_PIVOT_ERRO...
void coarse_to_fine(const IceModelVec3 &coarse, int i, int j, double *fine) const
std::vector< double > m_z
levels of the fine vertical grid
unsigned int m_ks
current system size; corresponds to the highest vertical level within the ice
std::vector< double > m_u
u-component of the ice velocity
void init_column(int i, int j, double ice_thickness)
const IceModelVec3 & m_v3
int m_i
current column indexes
TridiagonalSystem * m_solver
std::vector< double > m_w
w-component of the ice velocity
std::vector< double > m_v
v-component of the ice velocity
Base class for tridiagonal systems in the ice.
static const double k
Definition: exactTestP.cc:45
static double S(unsigned n)
Definition: test_cube.c:58