tColumnSystem.hh - 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
---
tColumnSystem.hh (6324B)
---
1 // Copyright (C) 2009-2011, 2013, 2014, 2015, 2016, 2019 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 __columnSystem_hh
20 #define __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 */
82 class TridiagonalSystem {
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
105 void save_matrix(std::ostream &output,
106 unsigned int system_size,
107 const std::string &info) const;
108
109 void save_vector(std::ostream &output,
110 const std::vector<double> &v,
111 unsigned int system_size, const std::string &info) const;
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 */
141 class columnSystemCtx {
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);
146 ~columnSystemCtx();
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:
157 TridiagonalSystem *m_solver;
158
159 ColumnInterpolation *m_interp;
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
178 const IceModelVec3 &m_u3, &m_v3, &m_w3;
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 /* __columnSystem_hh */
192