tPoisson.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
---
tPoisson.cc (7879B)
---
1 /* Copyright (C) 2019, 2020 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 "Poisson.hh"
21 #include "pism/util/error_handling.hh"
22
23 namespace pism {
24
25 Poisson::Poisson(IceGrid::ConstPtr grid)
26 : m_grid(grid),
27 m_log(grid->ctx()->log()),
28 m_b(grid, "poisson_rhs", WITHOUT_GHOSTS),
29 m_x(grid, "poisson_x", WITHOUT_GHOSTS),
30 m_mask(grid, "poisson_mask", WITH_GHOSTS){
31
32 m_da = m_x.dm();
33
34 // PETSc objects and settings
35 {
36 PetscErrorCode ierr;
37 ierr = DMSetMatType(*m_da, MATAIJ);
38 PISM_CHK(ierr, "DMSetMatType");
39
40 ierr = DMCreateMatrix(*m_da, m_A.rawptr());
41 PISM_CHK(ierr, "DMCreateMatrix");
42
43 ierr = KSPCreate(m_grid->com, m_KSP.rawptr());
44 PISM_CHK(ierr, "KSPCreate");
45
46 ierr = KSPSetOptionsPrefix(m_KSP, "poisson_");
47 PISM_CHK(ierr, "KSPSetOptionsPrefix");
48
49 // Process options:
50 ierr = KSPSetFromOptions(m_KSP);
51 PISM_CHK(ierr, "KSPSetFromOptions");
52 }
53 }
54
55 /*!
56 * Solve the Poisson equation on the domain defined by `mask == 1` with Dirichlet BC
57 * provided in `bc` (used only where `mask == 0`, possibly redundant away from the domain)
58 * with the constant right hand side `rhs`.
59 *
60 * Set the mask to 2 to use zero Neumann BC.
61 */
62 int Poisson::solve(const IceModelVec2Int& mask, const IceModelVec2S& bc, double rhs,
63 bool reuse_matrix) {
64
65 PetscErrorCode ierr;
66
67 // make a ghosted copy of the mask
68 m_mask.copy_from(mask);
69
70 if (reuse_matrix) {
71 // Use non-zero initial guess. I assume that re-using the matrix means that the BC and
72 // RHS provided are close to the ones used in the previous call and the solution
73 // stored in m_x is a good initial guess for the current problem.
74 ierr = KSPSetInitialGuessNonzero(m_KSP, PETSC_TRUE);
75 PISM_CHK(ierr, "KSPSetInitialGuessNonzero");
76 } else {
77 ierr = KSPSetInitialGuessNonzero(m_KSP, PETSC_FALSE);
78 PISM_CHK(ierr, "KSPSetInitialGuessNonzero");
79
80 assemble_matrix(m_mask, m_A);
81 }
82
83 assemble_rhs(rhs, m_mask, bc, m_b);
84
85 // Call PETSc to solve linear system by iterative method.
86 ierr = KSPSetOperators(m_KSP, m_A, m_A);
87 PISM_CHK(ierr, "KSPSetOperator");
88
89 ierr = KSPSolve(m_KSP, m_b.vec(), m_x.vec());
90 PISM_CHK(ierr, "KSPSolve");
91
92 // Check if diverged
93 KSPConvergedReason reason;
94 ierr = KSPGetConvergedReason(m_KSP, &reason);
95 PISM_CHK(ierr, "KSPGetConvergedReason");
96
97 if (reason < 0) {
98 // KSP diverged
99 m_log->message(1,
100 "PISM ERROR: KSP iteration failed while solving the Poisson equation\n"
101 " reason = %d = '%s'\n",
102 reason, KSPConvergedReasons[reason]);
103
104 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "KSP iteration failed: %s",
105 KSPConvergedReasons[reason]);
106 }
107
108 // report on KSP success
109 PetscInt ksp_iterations = 0;
110 ierr = KSPGetIterationNumber(m_KSP, &ksp_iterations);
111 PISM_CHK(ierr, "KSPGetIterationNumber");
112
113 return ksp_iterations;
114 }
115
116 const IceModelVec2S& Poisson::solution() const {
117 return m_x;
118 }
119
120 // Maxima code deriving the discretization
121 //
122 // /* Shift in x and y directions. */
123 // shift(expr, dx, dy) := op(expr)[args(expr)[1] + dx, args(expr)[2] + dy]$
124 //
125 // d_px(var) := (shift(var, 1, 0) - var)$
126 // d_mx(var) := (var - shift(var, -1, 0))$
127 //
128 // d_py(var) := (shift(var, 0, 1) - var)$
129 // d_my(var) := (var - shift(var, 0, -1))$
130 //
131 // constants : [C_x = 1 / (dx^2), C_y = 1 / (dy^2)]$
132 //
133 // /* discretization of -\nabla \dot (D \nabla W) */
134 // L: (E * d_px(W[i, j]) - W * d_mx(W[i, j])) / dx^2 +
135 // (N * d_py(W[i, j]) - S * d_my(W[i, j])) / dy^2$
136 //
137 // /* discretization of the Poisson equation */
138 // eq: - L = f;
139 //
140 // /* cleaned up equation */
141 // s : map(lambda([x,y], solve(x, y)[1]), constants, [dx^2, dy^2]);
142 // eq2 : at(eq, s);
143 //
144 // /* compute matrix coefficients */
145 // for m: -1 thru 1 do (for n: -1 thru 1 do
146 // (c[2 - n, m + 2] : factor(ratcoef(lhs(eq2), W[i + m, j + n]))))$
147 //
148 // /* print results */
149 // A : genmatrix(c, 3, 3);
150 //
151 // b : rhs(eq2);
152 //
153 // print(''out)$
154 void Poisson::assemble_matrix(const IceModelVec2Int &mask, Mat A) {
155 PetscErrorCode ierr = 0;
156
157 const double
158 dx = m_grid->dx(),
159 dy = m_grid->dy(),
160 C_x = 1.0 / (dx * dx),
161 C_y = 1.0 / (dy * dy);
162
163 const int
164 nrow = 1,
165 ncol = 5,
166 Mx = m_grid->Mx(),
167 My = m_grid->My();
168
169 ierr = MatZeroEntries(A); PISM_CHK(ierr, "MatZeroEntries");
170
171 IceModelVec::AccessList list{&mask};
172
173 /* matrix assembly loop */
174 ParallelSection loop(m_grid->com);
175 try {
176 MatStencil row, col[ncol];
177 row.c = 0;
178
179 for (int m = 0; m < ncol; m++) {
180 col[m].c = 0;
181 }
182
183 for (Points p(*m_grid); p; p.next()) {
184 const int i = p.i(), j = p.j();
185
186 /* i indices */
187 const int I[] = {i, i - 1, i, i + 1, i};
188
189 /* j indices */
190 const int J[] = {j + 1, j, j, j, j - 1};
191
192 row.i = i;
193 row.j = j;
194
195 for (int m = 0; m < ncol; m++) {
196 col[m].i = I[m];
197 col[m].j = J[m];
198 }
199
200 auto M = mask.int_star(i, j);
201
202 if (M.ij == 1) {
203 // Regular location: use coefficients of the discretization of the Laplacian
204
205 // Use zero Neumann BC if a neighbor is marked as a Neumann boundary
206 double
207 N = M.n == 2 ? 0.0 : 1.0,
208 E = M.e == 2 ? 0.0 : 1.0,
209 W = M.w == 2 ? 0.0 : 1.0,
210 S = M.s == 2 ? 0.0 : 1.0;
211
212 // Use zero Neumann BC at edges of the computational domain
213 {
214 N = j == My - 1 ? 0.0 : 1.0;
215 E = i == Mx - 1 ? 0.0 : 1.0;
216 W = i == 0 ? 0.0 : 1.0;
217 S = j == 0 ? 0.0 : 1.0;
218 }
219
220 // discretization of the Laplacian
221 double L[ncol] = {- N * C_y,
222 - W * C_x, (W + E) * C_x + (N + S) * C_y, - E * C_x,
223 - S * C_y};
224
225 ierr = MatSetValuesStencil(A, nrow, &row, ncol, col, L, INSERT_VALUES);
226 PISM_CHK(ierr, "MatSetValuesStencil");
227 } else {
228 // Boundary or outside the domain: assemble trivial equations (1 on the diagonal,
229 // 0 elsewhere)
230 double D[ncol] = {0.0,
231 0.0, 1.0, 0.0,
232 0.0};
233
234 ierr = MatSetValuesStencil(A, nrow, &row, ncol, col, D, INSERT_VALUES);
235 PISM_CHK(ierr, "MatSetValuesStencil");
236 }
237
238 } // i,j-loop
239 } catch (...) {
240 loop.failed();
241 }
242 loop.check();
243
244 ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); PISM_CHK(ierr, "MatAssemblyBegin");
245 ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); PISM_CHK(ierr, "MatAssemblyif");
246
247 #if (Pism_DEBUG==1)
248 ierr = MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
249 PISM_CHK(ierr, "MatSetOption");
250 #endif
251
252 }
253
254 void Poisson::assemble_rhs(double rhs,
255 const IceModelVec2Int &mask,
256 const IceModelVec2S &bc,
257 IceModelVec2S &b) {
258 IceModelVec::AccessList list{&mask, &bc, &b};
259
260 for (Points p(*m_grid); p; p.next()) {
261 const int i = p.i(), j = p.j();
262
263 if (mask.as_int(i, j) == 1) {
264 // inside the domain
265 b(i, j) = rhs;
266 } else {
267 // at the boundary or outside the domain
268 b(i, j) = bc(i, j);
269 }
270 }
271 }
272
273 } // end of namespace pism