Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.2-d6b3a29ca committed by Constantine Khrulev on 2025-03-28
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AgeColumnSystem.cc
Go to the documentation of this file.
1/* Copyright (C) 2016, 2017, 2023 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 "pism/age/AgeColumnSystem.hh"
21#include "pism/util/error_handling.hh"
22
23namespace pism {
24
25AgeColumnSystem::AgeColumnSystem(const std::vector<double>& storage_grid,
26 const std::string &my_prefix,
27 double dx, double dy, double dt,
28 const array::Array3D &age,
29 const array::Array3D &u3,
30 const array::Array3D &v3,
31 const array::Array3D &w3)
32 : columnSystemCtx(storage_grid, my_prefix, dx, dy, dt, u3, v3, w3),
33 m_age3(age) {
34
35 size_t Mz = m_z.size();
36 m_A.resize(Mz);
37 m_A_n.resize(Mz);
38 m_A_e.resize(Mz);
39 m_A_s.resize(Mz);
40 m_A_w.resize(Mz);
41
42 m_nu = m_dt / m_dz; // derived constant
43}
44
45void AgeColumnSystem::init(int i, int j, double thickness) {
46 init_column(i, j, thickness);
47
48 if (m_ks == 0) {
49 return;
50 }
51
52 coarse_to_fine(m_u3, i, j, &m_u[0]);
53 coarse_to_fine(m_v3, i, j, &m_v[0]);
54 coarse_to_fine(m_w3, i, j, &m_w[0]);
55
61}
62
63//! First-order upwind scheme with implicit in the vertical: one column solve.
64/*!
65 The PDE being solved is
66 \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]
67 */
68void AgeColumnSystem::solve(std::vector<double> &x) {
69
71
72 // set up system: 0 <= k < m_ks
73 for (unsigned int k = 0; k < m_ks; k++) {
74 // do lowest-order upwinding, explicitly for horizontal
75 S.RHS(k) = (m_u[k] < 0 ?
76 m_u[k] * (m_A_e[k] - m_A[k]) / m_dx :
77 m_u[k] * (m_A[k] - m_A_w[k]) / m_dx);
78 S.RHS(k) += (m_v[k] < 0 ?
79 m_v[k] * (m_A_n[k] - m_A[k]) / m_dy :
80 m_v[k] * (m_A[k] - m_A_s[k]) / m_dy);
81 // note it is the age eqn: dage/dt = 1.0 and we have moved the hor.
82 // advection terms over to right:
83 S.RHS(k) = m_A[k] + m_dt * (1.0 - S.RHS(k));
84
85 // do lowest-order upwinding, *implicitly* for vertical
86 double AA = m_nu * m_w[k];
87 if (k > 0) {
88 if (AA >= 0) { // upward velocity
89 S.L(k) = - AA;
90 S.D(k) = 1.0 + AA;
91 S.U(k) = 0.0;
92 } else { // downward velocity; note -AA >= 0
93 S.L(k) = 0.0;
94 S.D(k) = 1.0 - AA;
95 S.U(k) = + AA;
96 }
97 } else { // k == 0 case
98 // note L[0] is not used
99 if (AA > 0) { // if strictly upward velocity apply boundary condition:
100 // age = 0 because ice is being added to base
101 S.D(0) = 1.0;
102 S.U(0) = 0.0;
103 S.RHS(0) = 0.0;
104 } else { // downward velocity; note -AA >= 0
105 S.D(0) = 1.0 - AA;
106 S.U(0) = + AA;
107 // keep rhs[0] as is
108 }
109 }
110 } // done "set up system: 0 <= k < m_ks"
111
112 // surface b.c. at m_ks
113 if (m_ks > 0) {
114 S.L(m_ks) = 0;
115 S.D(m_ks) = 1.0; // ignore U[m_ks]
116 S.RHS(m_ks) = 0.0; // age zero at surface
117 }
118
119 // solve it
120 try {
121 S.solve(m_ks + 1, x);
122 }
123 catch (RuntimeError &e) {
124 e.add_context("solving the tri-diagonal system (AgeColumnSystem) at (%d, %d)\n"
125 "saving system to m-file... ", m_i, m_j);
127 throw;
128 }
129
130 // x[k] contains age for k=0,...,ks, but set age of ice above (and
131 // at) surface to zero years
132 for (unsigned int k = m_ks + 1; k < x.size(); k++) {
133 x[k] = 0.0;
134 }
135}
136
137} // end of namespace pism
std::vector< double > m_A_s
std::vector< double > m_A_n
std::vector< double > m_A_w
AgeColumnSystem(const std::vector< double > &storage_grid, const std::string &my_prefix, double dx, double dy, double dt, const array::Array3D &age, const array::Array3D &u3, const array::Array3D &v3, const array::Array3D &w3)
std::vector< double > m_A_e
std::vector< double > m_A
void init(int i, int j, double thickness)
const array::Array3D & m_age3
void solve(std::vector< double > &x)
First-order upwind scheme with implicit in the vertical: one column solve.
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.
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition Array3D.hh:33
const array::Array3D & m_w3
void coarse_to_fine(const array::Array3D &input, int i, int j, double *output) const
void reportColumnZeroPivotErrorMFile(unsigned int M)
Write system matrix and right-hand-side into an Python script. The file name contains ZERO_PIVOT_ERRO...
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
const array::Array3D & m_u3
pointers to 3D velocity components
void init_column(int i, int j, double ice_thickness)
const array::Array3D & 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:42
static double S(unsigned n)
Definition test_cube.c:58