PISM, A Parallel Ice Sheet Model  stable v2.0.6 committed by Constantine Khrulev on 2023-01-23 15:14:38 -0900
ColumnSystem.hh
Go to the documentation of this file.
1 // Copyright (C) 2009-2011, 2013, 2014, 2015, 2016, 2019, 2021 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 #ifndef PISM_COLUMNSYSTEM_HH
20 #define PISM_COLUMNSYSTEM_HH
21 
22 #include <string>
23 #include <ostream>
24 #include <vector>
25 
26 namespace pism {
27 
28 //! Virtual base class. Abstracts a tridiagonal system to solve in a column of ice and/or bedrock.
29 /*!
30  Because both the age evolution and conservation of energy equations require us to set up
31  and solve a tridiagonal system of equations, this is structure is worth abstracting.
32 
33  This base class just holds the tridiagonal system and the ability to
34  solve it, but does not insert entries into the relevant matrix locations.
35  Derived classes will actually set up instances of the system.
36 
37  The sequence requires setting the column-independent (public) data members,
38  calling the initAllColumns() routine, and then setting up and solving
39  the system in each column.
40 
41  The tridiagonal algorithm here comes from Numerical Recipes in C
42  Section 2.4, page 50. It solves the system:
43 
44 @verbatim
45  [b_1 c_1 0 ... ] [ u_1 ] [ r_1 ]
46  [a_2 b_2 c_2 ... ] [ u_2 ] [ r_2 ]
47  [ ... ].[ ... ] = [ ... ]
48  [ ... a_{N-1} b_{N-1} c_{N-1}] [u_{N-1}] [ r_{N-1} ]
49  [ ... 0 a_N b_N ] [ u_N ] [ r_N ]
50 @endverbatim
51 
52  HOWEVER... the version in this code is different from Numerical
53  Recipes in two ways:
54 
55  - Indexing is zero-based
56  - Variables have been renamed.
57 
58 @verbatim
59  NR PISM
60  ==================
61  a L "Lower Diagonal" (L doesn't use index 0)
62  b D "Diagonal"
63  c U "Upper Diagonal"
64  u x
65  r rhs
66  bet b
67  j k
68  n n
69  gam work
70 @endverbatim
71 
72  Therefore... this version of the code solves the following problem:
73 
74 @verbatim
75  [D_0 U_0 0 ... ] [ x_0 ] [ r_0 ]
76  [L_1 D_1 U_1 ... ] [ x_1 ] [ r_1 ]
77  [ ... ].[ ... ] = [ ... ]
78  [ ... L_{N-2} D_{N-2} U_{N-2}] [x_{N-2}] [ r_{N-2} ]
79  [ ... 0 L_{N-1} D_{N-1}] [x_{N-1}] [ r_{N-1} ]
80 @endverbatim
81 */
83 public:
84  TridiagonalSystem(unsigned int max_size, const std::string &prefix);
85 
86  double norm1(unsigned int system_size) const;
87  double ddratio(unsigned int system_size) const;
88  void reset();
89 
90  // uses an output argument to allow re-using storage and avoid
91  // copying
92  void solve(unsigned int system_size, std::vector<double> &result);
93  void solve(unsigned int system_size, double *result);
94 
95  void save_system_with_solution(const std::string &filename,
96  unsigned int system_size,
97  const std::vector<double> &solution);
98 
99  //! Save the system to a stream using the ASCII MATLAB (Octave)
100  //! format. Virtual to allow saving more info in derived classes.
101  void save_system(std::ostream &output,
102  unsigned int system_size) const;
103 
104  void save_matrix(std::ostream &output,
105  unsigned int system_size,
106  const std::string &variable) const;
107 
108  static void save_vector(std::ostream &output,
109  const std::vector<double> &v,
110  unsigned int system_size,
111  const std::string &variable);
112 
113  std::string prefix() const;
114 
115  double& L(size_t i) {
116  return m_L[i];
117  }
118  double& D(size_t i) {
119  return m_D[i];
120  }
121  double& U(size_t i) {
122  return m_U[i];
123  }
124  double& RHS(size_t i) {
125  return m_rhs[i];
126  }
127 private:
128  unsigned int m_max_system_size; // maximum system size
129  std::vector<double> m_L, m_D, m_U, m_rhs, m_work; // vectors for tridiagonal system
130 
131  std::string m_prefix;
132 };
133 
134 class IceModelVec3;
135 class ColumnInterpolation;
136 
137 //! Base class for tridiagonal systems in the ice.
138 /*! Adds data members used in time-dependent systems with advection
139  (dx, dy, dz, dt, velocity components).
140  */
142 public:
143  columnSystemCtx(const std::vector<double>& storage_grid, const std::string &prefix,
144  double dx, double dy, double dt,
145  const IceModelVec3 &u3, const IceModelVec3 &v3, const IceModelVec3 &w3);
147 
148  void save_to_file(const std::vector<double> &x);
149  void save_to_file(const std::string &filename, const std::vector<double> &x);
150 
151  unsigned int ks() const;
152  double dz() const;
153  const std::vector<double>& z() const;
154  void fine_to_coarse(const std::vector<double> &fine, int i, int j,
155  IceModelVec3& coarse) const;
156 protected:
158 
160 
161  //! current system size; corresponds to the highest vertical level within the ice
162  unsigned int m_ks;
163  //! current column indexes
164  int m_i, m_j;
165 
166  double m_dx, m_dy, m_dz, m_dt;
167 
168  //! u-component of the ice velocity
169  std::vector<double> m_u;
170  //! v-component of the ice velocity
171  std::vector<double> m_v;
172  //! w-component of the ice velocity
173  std::vector<double> m_w;
174  //! levels of the fine vertical grid
175  std::vector<double> m_z;
176 
177  //! pointers to 3D velocity components
179 
180  void init_column(int i, int j, double ice_thickness);
181 
182  void reportColumnZeroPivotErrorMFile(unsigned int M);
183 
184  void init_fine_grid(const std::vector<double>& storage_grid);
185 
186  void coarse_to_fine(const IceModelVec3 &coarse, int i, int j, double* fine) const;
187 };
188 
189 } // end of namespace pism
190 
191 #endif /* PISM_COLUMNSYSTEM_HH */
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition: iceModelVec.hh:404
void save_matrix(std::ostream &output, unsigned int system_size, const std::string &variable) const
View the tridiagonal matrix.
void save_system(std::ostream &output, unsigned int system_size) const
View the tridiagonal system A x = b to an output stream, both A as a full matrix and b as a vector.
std::vector< double > m_work
std::string prefix() const
double & RHS(size_t i)
unsigned int m_max_system_size
TridiagonalSystem(unsigned int max_size, const std::string &prefix)
Allocate a tridiagonal system of maximum size max_system_size.
Definition: ColumnSystem.cc:49
void solve(unsigned int system_size, std::vector< double > &result)
double ddratio(unsigned int system_size) const
Compute diagonal-dominance ratio. If this is less than one then the matrix is strictly diagonally-dom...
double norm1(unsigned int system_size) const
Compute 1-norm, which is max sum of absolute values of columns.
Definition: ColumnSystem.cc:74
double & D(size_t i)
double & U(size_t i)
double & L(size_t i)
static void save_vector(std::ostream &output, const std::vector< double > &v, unsigned int system_size, const std::string &variable)
Utility for simple ascii view of a vector (one-dimensional column) quantity.
std::vector< double > m_D
std::vector< double > m_U
void reset()
Zero all entries.
Definition: ColumnSystem.cc:63
std::vector< double > m_L
void save_system_with_solution(const std::string &filename, unsigned int system_size, const std::vector< double > &solution)
std::vector< double > m_rhs
Virtual base class. Abstracts a tridiagonal system to solve in a column of ice and/or bedrock.
Definition: ColumnSystem.hh:82
void fine_to_coarse(const std::vector< double > &fine, int i, int j, IceModelVec3 &coarse) const
const std::vector< double > & z() const
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 save_to_file(const std::vector< double > &x)
Write system matrix, right-hand-side, and (provided) solution into Python script. Constructs file nam...
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
columnSystemCtx(const std::vector< double > &storage_grid, const std::string &prefix, double dx, double dy, double dt, const IceModelVec3 &u3, const IceModelVec3 &v3, const IceModelVec3 &w3)
A column system is a kind of a tridiagonal system.
unsigned int ks() const
void init_fine_grid(const std::vector< double > &storage_grid)
int m_i
current column indexes
TridiagonalSystem * m_solver
ColumnInterpolation * m_interp
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.