Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
PCTwoLevelsBlockDiagonal.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#include "PCTwoLevelsBlockDiagonal.hpp"
37#include "Exception.hpp"
38
39#include <iostream>
40
41PCTwoLevelsBlockDiagonal::PCTwoLevelsBlockDiagonal(KSP& rKspObject, std::vector<PetscInt>& rBathNodes)
42{
43 PCTwoLevelsBlockDiagonalCreate(rKspObject, rBathNodes);
45}
46
70
71void PCTwoLevelsBlockDiagonal::PCTwoLevelsBlockDiagonalCreate(KSP& rKspObject, std::vector<PetscInt>& rBathNodes)
72{
73 KSPGetPC(rKspObject, &mPetscPCObject);
74
75 Mat system_matrix, dummy;
76#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=5)
77 KSPGetOperators(rKspObject, &system_matrix, &dummy);
78#else
79 MatStructure flag;
80 KSPGetOperators(rKspObject, &system_matrix, &dummy, &flag);
81#endif
82
83 PetscInt num_rows, num_columns;
84 MatGetSize(system_matrix, &num_rows, &num_columns);
85 assert(num_rows==num_columns);
86
87 PetscInt num_local_rows, num_local_columns;
88 MatGetLocalSize(system_matrix, &num_local_rows, &num_local_columns);
89
90 // Odd number of rows: impossible in Bidomain.
91 // Odd number of local rows: impossible if V_m and phi_e for each node are stored in the same processor.
92 if ((num_rows%2 != 0) || (num_local_rows%2 != 0))
93 {
94 TERMINATE("Wrong matrix parallel layout detected in PCLDUFactorisation."); // LCOV_EXCL_LINE
95 }
96
97 // Allocate memory
98 unsigned subvector_num_rows = num_rows/2;
99 unsigned subvector_local_rows = num_local_rows/2;
100
101 unsigned subvector_num_rows_tissue = subvector_num_rows - rBathNodes.size();
102 unsigned subvector_local_rows_tissue = subvector_num_rows_tissue;
103
104 unsigned subvector_num_rows_bath = rBathNodes.size();
105 unsigned subvector_local_rows_bath = subvector_num_rows_bath;
106
107 assert(PetscTools::IsSequential());
108
109 mPCContext.x1_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
110 mPCContext.x21_subvector = PetscTools::CreateVec(subvector_num_rows_tissue, subvector_local_rows_tissue);
111 mPCContext.x22_subvector = PetscTools::CreateVec(subvector_num_rows_bath, subvector_local_rows_bath);
112 mPCContext.y1_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
113 mPCContext.y21_subvector = PetscTools::CreateVec(subvector_num_rows_tissue, subvector_local_rows_tissue);
114 mPCContext.y22_subvector = PetscTools::CreateVec(subvector_num_rows_bath, subvector_local_rows_bath);
115
116 // Define IS objects that will be used throughout the method
117 IS A11_all_rows;
118 ISCreateStride(PETSC_COMM_WORLD, num_rows/2, 0, 2, &A11_all_rows);
119
120 IS A22_all_rows;
121 PetscInt A22_size;
122 VecGetSize(mPCContext.x1_subvector, &A22_size);
124 ISCreateStride(PETSC_COMM_WORLD, A22_size, 1, 2, &A22_all_rows);
125
126 IS A22_bath_rows;
127 PetscInt* phi_e_bath_rows = new PetscInt[rBathNodes.size()];
128 for (unsigned index=0; index<rBathNodes.size(); index++)
129 {
130 phi_e_bath_rows[index] = 2*rBathNodes[index] + 1;
131 }
132#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 2) //PETSc 3.2 or later
133 ISCreateGeneral(PETSC_COMM_WORLD, rBathNodes.size(), phi_e_bath_rows, PETSC_USE_POINTER, &A22_bath_rows);
134 #else
135 ISCreateGeneralWithArray(PETSC_COMM_WORLD, rBathNodes.size(), phi_e_bath_rows, &A22_bath_rows);
136#endif
137
138 IS A22_tissue_rows;
139 ISDifference(A22_all_rows, A22_bath_rows, &A22_tissue_rows);
140
141 // Create scatter contexts
142 {
143 // Needed by VecScatterCreate in order to find out parallel layout.
144 Vec dummy_vec = PetscTools::CreateVec(num_rows, num_local_rows);
145
147 IS& A11_rows=A11_all_rows;
148 IS& A22_B1_rows=A22_tissue_rows;
149 IS& A22_B2_rows=A22_bath_rows;
150
151 IS all_vector;
152 ISCreateStride(PETSC_COMM_WORLD, num_rows/2, 0, 1, &all_vector);
153
154 IS tissue_vector;
155 ISCreateStride(PETSC_COMM_WORLD, (num_rows/2)-rBathNodes.size(), 0, 1, &tissue_vector);
156
157 IS bath_vector;
158 ISCreateStride(PETSC_COMM_WORLD, rBathNodes.size(), 0, 1, &bath_vector);
159
160 VecScatterCreate(dummy_vec, A11_rows, mPCContext.x1_subvector, all_vector, &mPCContext.A11_scatter_ctx);
161 VecScatterCreate(dummy_vec, A22_B1_rows, mPCContext.x21_subvector, tissue_vector, &mPCContext.A22_B1_scatter_ctx);
162 VecScatterCreate(dummy_vec, A22_B2_rows, mPCContext.x22_subvector, bath_vector, &mPCContext.A22_B2_scatter_ctx);
163
164 ISDestroy(PETSC_DESTROY_PARAM(all_vector));
165 ISDestroy(PETSC_DESTROY_PARAM(tissue_vector));
166 ISDestroy(PETSC_DESTROY_PARAM(bath_vector));
167
168 PetscTools::Destroy(dummy_vec);
169 }
170
171 // Get matrix sublock A11
172 {
173 // Work out local row range for subblock A11 (same as x1 or y1)
174 PetscInt low, high, global_size;
175 VecGetOwnershipRange(mPCContext.x1_subvector, &low, &high);
176 VecGetSize(mPCContext.x1_subvector, &global_size);
177 assert(global_size == num_rows/2);
178
179 IS A11_local_rows;
180 IS& A11_columns=A11_all_rows;
181 ISCreateStride(PETSC_COMM_WORLD, high-low, 2*low, 2, &A11_local_rows);
182
183#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 8) //PETSc 3.8 or later
184 MatCreateSubMatrix(system_matrix, A11_local_rows, A11_columns,
185 MAT_INITIAL_MATRIX, &mPCContext.A11_matrix_subblock);
186#else
187 MatGetSubMatrix(system_matrix, A11_local_rows, A11_columns,
188 MAT_INITIAL_MATRIX, &mPCContext.A11_matrix_subblock);
189#endif
190
191 ISDestroy(PETSC_DESTROY_PARAM(A11_local_rows));
192 }
193
194 // Get matrix sublock A22_B1
195 {
196// // Work out local row range for subblock A22 (same as x2 or y2)
197// PetscInt low, high, global_size;
198// VecGetOwnershipRange(mPCContext.x21_subvector, &low, &high);
199// VecGetSize(mPCContext.x21_subvector, &global_size);
200// assert(global_size == (num_rows/2) - (PetscInt) rBathNodes.size());
201
202 assert(PetscTools::IsSequential());
203 IS& A22_B1_local_rows = A22_tissue_rows; // wrong in parallel, need to give local rows
204 IS& A22_B1_columns = A22_tissue_rows;
205
206#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 8) //PETSc 3.8 or later
207 MatCreateSubMatrix(system_matrix, A22_B1_local_rows, A22_B1_columns,
208 MAT_INITIAL_MATRIX, &mPCContext.A22_B1_matrix_subblock);
209#else
210 MatGetSubMatrix(system_matrix, A22_B1_local_rows, A22_B1_columns,
211 MAT_INITIAL_MATRIX, &mPCContext.A22_B1_matrix_subblock);
212#endif
213
214 }
215
216 // Get matrix sublock A22_B2
217 {
218// // Work out local row range for subblock A22 (same as x2 or y2)
219// PetscInt low, high, global_size;
220// VecGetOwnershipRange(mPCContext.x21_subvector, &low, &high);
221// VecGetSize(mPCContext.x21_subvector, &global_size);
222// assert(global_size == (num_rows/2) - (PetscInt) rBathNodes.size());
223
224 assert(PetscTools::IsSequential());
225 IS& A22_B2_local_rows = A22_bath_rows; // wrong in parallel, need to give local rows
226 IS& A22_B2_columns = A22_bath_rows;
227
228#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 8) //PETSc 3.8 or later
229 MatCreateSubMatrix(system_matrix, A22_B2_local_rows, A22_B2_columns,
230 MAT_INITIAL_MATRIX, &mPCContext.A22_B2_matrix_subblock);
231#else
232 MatGetSubMatrix(system_matrix, A22_B2_local_rows, A22_B2_columns,
233 MAT_INITIAL_MATRIX, &mPCContext.A22_B2_matrix_subblock);
234#endif
235 }
236
237 ISDestroy(PETSC_DESTROY_PARAM(A11_all_rows));
238 ISDestroy(PETSC_DESTROY_PARAM(A22_all_rows));
239 ISDestroy(PETSC_DESTROY_PARAM(A22_bath_rows));
240 delete[] phi_e_bath_rows;
241 ISDestroy(PETSC_DESTROY_PARAM(A22_tissue_rows));
242
243 // Register call-back function and its context
244 PCSetType(mPetscPCObject, PCSHELL);
245#if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
246 PCShellSetApply(mPetscPCObject, PCTwoLevelsBlockDiagonalApply, (void*) &mPCContext);
247#else
248 // Register PC context so it gets passed to PCTwoLevelsBlockDiagonalApply
249 PCShellSetContext(mPetscPCObject, &mPCContext);
250
251 // Register call-back function
252 PCShellSetApply(mPetscPCObject, PCTwoLevelsBlockDiagonalApply);
253#endif
254}
255
257{
258 // These options will get read by PCSetFromOptions
259 PetscTools::SetOption("-pc_hypre_boomeramg_max_iter", "1");
260 PetscTools::SetOption("-pc_hypre_boomeramg_strong_threshold", "0.0");
261 PetscTools::SetOption("-pc_hypre_type", "boomeramg");
262
263 // Set up amg preconditioner for block A11
264 PCCreate(PETSC_COMM_WORLD, &(mPCContext.PC_amg_A11));
265 PCSetType(mPCContext.PC_amg_A11, PCBJACOBI);
266#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=5)
268#else
269 PCSetOperators(mPCContext.PC_amg_A11, mPCContext.A11_matrix_subblock, mPCContext.A11_matrix_subblock, DIFFERENT_NONZERO_PATTERN);// SAME_PRECONDITIONER);
270#endif
271 PCSetFromOptions(mPCContext.PC_amg_A11);
272 PCSetUp(mPCContext.PC_amg_A11);
273
274 // Set up amg preconditioner for block A22_B1
275 PCCreate(PETSC_COMM_WORLD, &(mPCContext.PC_amg_A22_B1));
276 PCSetType(mPCContext.PC_amg_A22_B1, PCBJACOBI);
277#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=5)
279#else
280 PCSetOperators(mPCContext.PC_amg_A22_B1, mPCContext.A22_B1_matrix_subblock, mPCContext.A22_B1_matrix_subblock, DIFFERENT_NONZERO_PATTERN);// SAME_PRECONDITIONER);
281#endif
282 PCSetFromOptions(mPCContext.PC_amg_A22_B1);
283 PCSetUp(mPCContext.PC_amg_A22_B1);
284
285 // Set up amg preconditioner for block A22_B2
286 PCCreate(PETSC_COMM_WORLD, &(mPCContext.PC_amg_A22_B2));
287 PCSetType(mPCContext.PC_amg_A22_B2, PCHYPRE);
288 //PCHYPRESetType(mPCContext.PC_amg_A22_B2, "boomeramg");
289 PetscTools::SetOption("-pc_hypre_type", "boomeramg");
290
291 PetscTools::SetOption("-pc_hypre_boomeramg_max_iter", "1");
292 PetscTools::SetOption("-pc_hypre_boomeramg_strong_threshold", "0.0");
293 PetscTools::SetOption("-pc_hypre_boomeramg_coarsen_type", "HMIS");
294
295#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=5)
297#else
298 PCSetOperators(mPCContext.PC_amg_A22_B2, mPCContext.A22_B2_matrix_subblock, mPCContext.A22_B2_matrix_subblock, DIFFERENT_NONZERO_PATTERN);// SAME_PRECONDITIONER);
299#endif
300 PCSetFromOptions(mPCContext.PC_amg_A22_B2);
301 PCSetUp(mPCContext.PC_amg_A22_B2);
302}
303
304#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 1) //PETSc 3.1 or later
305PetscErrorCode PCTwoLevelsBlockDiagonalApply(PC pc_object, Vec x, Vec y)
306{
307 void* pc_context;
308
309 PCShellGetContext(pc_object, &pc_context);
310#else
311PetscErrorCode PCTwoLevelsBlockDiagonalApply(void* pc_context, Vec x, Vec y)
312{
313#endif
314
315 // Cast the context pointer to PCTwoLevelsBlockDiagonalContext
317 assert(block_diag_context!=nullptr);
318
319 /*
320 * Scatter x = [x1 x21 x22]'
321 */
322//PETSc-3.x.x or PETSc-2.3.3
323#if ((PETSC_VERSION_MAJOR == 3) || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 3 && PETSC_VERSION_SUBMINOR == 3)) //2.3.3 or 3.x.x
324 VecScatterBegin(block_diag_context->A11_scatter_ctx, x, block_diag_context->x1_subvector, INSERT_VALUES, SCATTER_FORWARD);
325 VecScatterEnd(block_diag_context->A11_scatter_ctx, x, block_diag_context->x1_subvector, INSERT_VALUES, SCATTER_FORWARD);
326#else
327 VecScatterBegin(x, block_diag_context->x1_subvector, INSERT_VALUES, SCATTER_FORWARD, block_diag_context->A11_scatter_ctx);
328 VecScatterEnd(x, block_diag_context->x1_subvector, INSERT_VALUES, SCATTER_FORWARD, block_diag_context->A11_scatter_ctx);
329#endif
330
331//PETSc-3.x.x or PETSc-2.3.3
332#if ((PETSC_VERSION_MAJOR == 3) || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 3 && PETSC_VERSION_SUBMINOR == 3)) //2.3.3 or 3.x.x
333 VecScatterBegin(block_diag_context->A22_B1_scatter_ctx, x, block_diag_context->x21_subvector, INSERT_VALUES, SCATTER_FORWARD);
334 VecScatterEnd(block_diag_context->A22_B1_scatter_ctx, x, block_diag_context->x21_subvector, INSERT_VALUES, SCATTER_FORWARD);
335#else
336 VecScatterBegin(x, block_diag_context->x21_subvector, INSERT_VALUES, SCATTER_FORWARD, block_diag_context->A22_B1_scatter_ctx);
337 VecScatterEnd(x, block_diag_context->x21_subvector, INSERT_VALUES, SCATTER_FORWARD, block_diag_context->A22_B1_scatter_ctx);
338#endif
339
340//PETSc-3.x.x or PETSc-2.3.3
341#if ((PETSC_VERSION_MAJOR == 3) || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 3 && PETSC_VERSION_SUBMINOR == 3)) //2.3.3 or 3.x.x
342 VecScatterBegin(block_diag_context->A22_B2_scatter_ctx, x, block_diag_context->x22_subvector, INSERT_VALUES, SCATTER_FORWARD);
343 VecScatterEnd(block_diag_context->A22_B2_scatter_ctx, x, block_diag_context->x22_subvector, INSERT_VALUES, SCATTER_FORWARD);
344#else
345 VecScatterBegin(x, block_diag_context->x22_subvector, INSERT_VALUES, SCATTER_FORWARD, block_diag_context->A22_B2_scatter_ctx);
346 VecScatterEnd(x, block_diag_context->x22_subvector, INSERT_VALUES, SCATTER_FORWARD, block_diag_context->A22_B2_scatter_ctx);
347#endif
348
349 /*
350 * y1 = ILU(A11)*x1
351 * y21 = ILU(A22)*x21
352 * y22 = AMG(A22)*x22
353 */
354 PCApply(block_diag_context->PC_amg_A11, block_diag_context->x1_subvector, block_diag_context->y1_subvector);
355 PCApply(block_diag_context->PC_amg_A22_B1, block_diag_context->x21_subvector, block_diag_context->y21_subvector);
356 PCApply(block_diag_context->PC_amg_A22_B2, block_diag_context->x22_subvector, block_diag_context->y22_subvector);
357
358 /*
359 * Gather y = [y1 y21 y22]'
360 */
361//PETSc-3.x.x or PETSc-2.3.3
362#if ((PETSC_VERSION_MAJOR == 3) || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 3 && PETSC_VERSION_SUBMINOR == 3)) //2.3.3 or 3.x.x
363 VecScatterBegin(block_diag_context->A11_scatter_ctx, block_diag_context->y1_subvector, y, INSERT_VALUES, SCATTER_REVERSE);
364 VecScatterEnd(block_diag_context->A11_scatter_ctx, block_diag_context->y1_subvector, y, INSERT_VALUES, SCATTER_REVERSE);
365#else
366 VecScatterBegin(block_diag_context->y1_subvector, y, INSERT_VALUES, SCATTER_REVERSE, block_diag_context->A11_scatter_ctx);
367 VecScatterEnd(block_diag_context->y1_subvector, y, INSERT_VALUES, SCATTER_REVERSE, block_diag_context->A11_scatter_ctx);
368#endif
369
370//PETSc-3.x.x or PETSc-2.3.3
371#if ((PETSC_VERSION_MAJOR == 3) || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 3 && PETSC_VERSION_SUBMINOR == 3)) //2.3.3 or 3.x.x
372 VecScatterBegin(block_diag_context->A22_B1_scatter_ctx, block_diag_context->y21_subvector, y, INSERT_VALUES, SCATTER_REVERSE);
373 VecScatterEnd(block_diag_context->A22_B1_scatter_ctx, block_diag_context->y21_subvector, y, INSERT_VALUES, SCATTER_REVERSE);
374#else
375 VecScatterBegin(block_diag_context->y21_subvector, y, INSERT_VALUES, SCATTER_REVERSE, block_diag_context->A22_B1_scatter_ctx);
376 VecScatterEnd(block_diag_context->y21_subvector, y, INSERT_VALUES, SCATTER_REVERSE, block_diag_context->A22_B1_scatter_ctx);
377#endif
378
379//PETSc-3.x.x or PETSc-2.3.3
380#if ((PETSC_VERSION_MAJOR == 3) || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 3 && PETSC_VERSION_SUBMINOR == 3)) //2.3.3 or 3.x.x
381 VecScatterBegin(block_diag_context->A22_B2_scatter_ctx, block_diag_context->y22_subvector, y, INSERT_VALUES, SCATTER_REVERSE);
382 VecScatterEnd(block_diag_context->A22_B2_scatter_ctx, block_diag_context->y22_subvector, y, INSERT_VALUES, SCATTER_REVERSE);
383#else
384 VecScatterBegin(block_diag_context->y22_subvector, y, INSERT_VALUES, SCATTER_REVERSE, block_diag_context->A22_B2_scatter_ctx);
385 VecScatterEnd(block_diag_context->y22_subvector, y, INSERT_VALUES, SCATTER_REVERSE, block_diag_context->A22_B2_scatter_ctx);
386#endif
387
388 return 0;
389}
#define TERMINATE(message)
#define PETSC_DESTROY_PARAM(x)
void PCTwoLevelsBlockDiagonalCreate(KSP &rKspObject, std::vector< PetscInt > &rBathNodes)
PCTwoLevelsBlockDiagonalContext mPCContext
PCTwoLevelsBlockDiagonal(KSP &rKspObject, std::vector< PetscInt > &rBathNodes)
static void Destroy(Vec &rVec)
static bool IsSequential()
static Vec CreateVec(int size, int localSize=PETSC_DECIDE, bool ignoreOffProcEntries=true)
static void SetOption(const char *pOptionName, const char *pOptionValue)