Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
PCBlockDiagonal.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 <iostream>
37
38#include "PetscVecTools.hpp" // Includes Ublas so must come first
39#include "PCBlockDiagonal.hpp"
40#include "Exception.hpp"
41#include "Warnings.hpp"
42
43#ifdef TRACE_KSP
44#include "Timer.hpp"
45#endif
46
48{
49#ifdef TRACE_KSP
50 mPCContext.mScatterTime = 0.0;
51 mPCContext.mA1PreconditionerTime = 0.0;
52 mPCContext.mA2PreconditionerTime = 0.0;
53 mPCContext.mGatherTime = 0.0;;
54#endif
55
56 PCBlockDiagonalCreate(rKspObject);
58}
59
60PCBlockDiagonal::~PCBlockDiagonal()
61{
62#ifdef TRACE_KSP
64 {
65 std::cout << " -- Block diagonal preconditioner profile information: " << std::endl;
66 std::cout << "\t mScatterTime: " << mPCContext.mScatterTime << std::endl;
67 std::cout << "\t mA1PreconditionerTime: " << mPCContext.mA1PreconditionerTime << std::endl;
68 std::cout << "\t mA2PreconditionerTime: " << mPCContext.mA2PreconditionerTime << std::endl;
69 std::cout << "\t mGatherTime: " << mPCContext.mGatherTime << std::endl;
70 }
71#endif
72
75
78
81
84
87}
88
90{
91 KSPGetPC(rKspObject, &mPetscPCObject);
92
93 Mat system_matrix, dummy;
94#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=5)
95 KSPGetOperators(rKspObject, &system_matrix, &dummy);
96#else
97 MatStructure flag;
98 KSPGetOperators(rKspObject, &system_matrix, &dummy, &flag);
99#endif
100
101 PetscInt num_rows, num_columns;
102 MatGetSize(system_matrix, &num_rows, &num_columns);
103 assert(num_rows==num_columns);
104
105 PetscInt num_local_rows, num_local_columns;
106 MatGetLocalSize(system_matrix, &num_local_rows, &num_local_columns);
107
108 // Odd number of rows: impossible in Bidomain.
109 // Odd number of local rows: impossible if V_m and phi_e for each node are stored in the same processor.
110 if ((num_rows%2 != 0) || (num_local_rows%2 != 0))
111 {
112 TERMINATE("Wrong matrix parallel layout detected in PCBlockDiagonal."); // LCOV_EXCL_LINE
113 }
114
115 // Allocate memory
116 unsigned subvector_num_rows = num_rows/2;
117 unsigned subvector_local_rows = num_local_rows/2;
118 mPCContext.x1_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
119 mPCContext.x2_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
120 mPCContext.y1_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
121 mPCContext.y2_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
122
123 // Create scatter contexts
124 {
125 // Needed by SetupInterleavedVectorScatterGather in order to find out parallel layout.
126 Vec dummy_vec = PetscTools::CreateVec(num_rows, num_local_rows);
127
129
130 PetscTools::Destroy(dummy_vec);
131 }
132
133 // Get matrix sublock A11
134 {
135 // Work out local row range for subblock A11 (same as x1 or y1)
136 PetscInt low, high, global_size;
137 VecGetOwnershipRange(mPCContext.x1_subvector, &low, &high);
138 VecGetSize(mPCContext.x1_subvector, &global_size);
139 assert(global_size == num_rows/2);
140
141 IS A11_local_rows;
142 IS A11_columns;
143 ISCreateStride(PETSC_COMM_WORLD, high-low, 2*low, 2, &A11_local_rows);
144 ISCreateStride(PETSC_COMM_WORLD, global_size, 0, 2, &A11_columns);
145
146#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 8) //PETSc 3.8 or later
147 MatCreateSubMatrix(system_matrix, A11_local_rows, A11_local_rows,
148 MAT_INITIAL_MATRIX, &mPCContext.A11_matrix_subblock);
149#else
150 MatGetSubMatrix(system_matrix, A11_local_rows, A11_local_rows,
151 MAT_INITIAL_MATRIX, &mPCContext.A11_matrix_subblock);
152#endif
153
154
155 ISDestroy(PETSC_DESTROY_PARAM(A11_local_rows));
156 ISDestroy(PETSC_DESTROY_PARAM(A11_columns));
157 }
158
159 // Get matrix sublock A22
160 {
161 // Work out local row range for subblock A22 (same as x2 or y2)
162 PetscInt low, high, global_size;
163 VecGetOwnershipRange(mPCContext.x2_subvector, &low, &high);
164 VecGetSize(mPCContext.x2_subvector, &global_size);
165 assert(global_size == num_rows/2);
166
167 IS A22_local_rows;
168 IS A22_columns;
169 ISCreateStride(PETSC_COMM_WORLD, high-low, 2*low+1, 2, &A22_local_rows);
170 ISCreateStride(PETSC_COMM_WORLD, global_size, 1, 2, &A22_columns);
171
172#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 8) //PETSc 3.8 or later
173 MatCreateSubMatrix(system_matrix, A22_local_rows, A22_local_rows,
174 MAT_INITIAL_MATRIX, &mPCContext.A22_matrix_subblock);
175#else
176 MatGetSubMatrix(system_matrix, A22_local_rows, A22_local_rows,
177 MAT_INITIAL_MATRIX, &mPCContext.A22_matrix_subblock);
178#endif
179
180 ISDestroy(PETSC_DESTROY_PARAM(A22_local_rows));
181 ISDestroy(PETSC_DESTROY_PARAM(A22_columns));
182 }
183
184 // Register call-back function and its context
185 PCSetType(mPetscPCObject, PCSHELL);
186#if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
187 PCShellSetApply(mPetscPCObject, PCBlockDiagonalApply, (void*) &mPCContext);
188#else
189 // Register PC context so it gets passed to PCBlockDiagonalApply
190 PCShellSetContext(mPetscPCObject, &mPCContext);
191
192 // Register call-back function
193 PCShellSetApply(mPetscPCObject, PCBlockDiagonalApply);
194#endif
195
196}
197
199{
200 // These options will get read by PCSetFromOptions
201// PetscTools::SetOption("-pc_hypre_boomeramg_max_iter", "1");
202// PetscTools::SetOption("-pc_hypre_boomeramg_strong_threshold", "0.0");
203// PetscTools::SetOption("-pc_hypre_type", "boomeramg");
204
205 // Set up amg preconditioner for block A11
206 PCCreate(PETSC_COMM_WORLD, &(mPCContext.PC_amg_A11));
207
208 // PCSetType(mPCContext.PC_amg_A11, PCBJACOBI);
209
210 // We are expecting an error from PETSC on systems that don't have the hypre library, so suppress it
211 // in case it aborts
212 PetscPushErrorHandler(PetscIgnoreErrorHandler, nullptr);
213 PetscErrorCode pc_set_error = PCSetType(mPCContext.PC_amg_A11, PCHYPRE);
214 if (pc_set_error != 0)
215 {
216 WARNING("PETSc hypre preconditioning library is not installed"); // LCOV_EXCL_LINE
217 }
218 // Stop supressing error
219 PetscPopErrorHandler();
220
221#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<=5)
222 PetscTools::SetOption("-pc_hypre_type", "euclid");
223 PetscTools::SetOption("-pc_hypre_euclid_levels", "0");
224#else
225/* euclid was removed in 3.6. Use the previously-commented-out alternative */
226 PetscTools::SetOption("-pc_hypre_type", "boomeramg");
227 PetscTools::SetOption("-pc_hypre_boomeramg_max_iter", "1");
228 PetscTools::SetOption("-pc_hypre_boomeramg_strong_threshold", "0.0");
229 PetscTools::SetOption("-pc_hypre_boomeramg_coarsen_type", "HMIS");
230#endif
231
232#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=5)
233 // Attempt to emulate SAME_PRECONDITIONER below
234 PCSetReusePreconditioner(mPCContext.PC_amg_A11, PETSC_TRUE);
236#else
238#endif
239 PCSetFromOptions(mPCContext.PC_amg_A11);
240 PCSetUp(mPCContext.PC_amg_A11);
241
242 // Set up amg preconditioner for block A22
243 PCCreate(PETSC_COMM_WORLD, &(mPCContext.PC_amg_A22));
244
245 /* Full AMG in the block */
246 PetscPushErrorHandler(PetscIgnoreErrorHandler, nullptr);
247 PCSetType(mPCContext.PC_amg_A22, PCHYPRE);
248 // Stop supressing error
249 PetscPopErrorHandler();
250
251 //PCHYPRESetType(mPCContext.PC_amg_A22, "boomeramg");
252 PetscTools::SetOption("-pc_hypre_type", "boomeramg");
253 PetscTools::SetOption("-pc_hypre_boomeramg_max_iter", "1");
254 PetscTools::SetOption("-pc_hypre_boomeramg_strong_threshold", "0.0");
255 PetscTools::SetOption("-pc_hypre_boomeramg_coarsen_type", "HMIS");
256 //PetscTools::SetOption("-pc_hypre_boomeramg_relax_type_all","Jacobi");
257 //PetscTools::SetOption("-pc_hypre_boomeramg_max_levels","10");
258 //PetscTools::SetOption("-pc_hypre_boomeramg_agg_nl", "1");
259 //PetscTools::SetOption("-pc_hypre_boomeramg_print_statistics","");
260 //PetscTools::SetOption("-pc_hypre_boomeramg_interp_type","standard");
261 //PetscTools::SetOption("-pc_hypre_boomeramg_interp_type","ext+i");
262
263 // PCHYPRESetType(mPCContext.PC_amg_A22, "euclid");
264
265 /* Block Jacobi with AMG at each block */
266 // PCSetType(mPCContext.PC_amg_A22, PCBJACOBI);
267
268 // PetscTools::SetOption("-sub_pc_type", "hypre");
269
270 // PetscTools::SetOption("-sub_pc_hypre_type", "boomeramg");
271 // PetscTools::SetOption("-sub_pc_hypre_boomeramg_max_iter", "1");
272 // PetscTools::SetOption("-sub_pc_hypre_boomeramg_strong_threshold", "0.0");
273
274 // PetscTools::SetOption("-pc_hypre_type", "boomeramg");
275 // PetscTools::SetOption("-pc_hypre_boomeramg_max_iter", "1");
276 // PetscTools::SetOption("-pc_hypre_boomeramg_strong_threshold", "0.0");
277
278 /* Additive Schwarz with AMG at each block */
279// PCSetType(mPCContext.PC_amg_A22, PCASM);
280
281// PetscTools::SetOption("-pc_asm_type", "basic");
282// PetscTools::SetOption("-pc_asm_overlap", "1");
283
284// PetscTools::SetOption("-sub_ksp_type", "preonly");
285
286// PetscTools::SetOption("-sub_pc_type", "hypre");
287
288// PetscTools::SetOption("-sub_pc_hypre_type", "boomeramg");
289// PetscTools::SetOption("-sub_pc_hypre_boomeramg_max_iter", "1");
290// PetscTools::SetOption("-sub_pc_hypre_boomeramg_strong_threshold", "0.0");
291
292#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=5)
293 // Attempt to emulate SAME_PRECONDITIONER below
294 PCSetReusePreconditioner(mPCContext.PC_amg_A22, PETSC_TRUE);
296#else
298#endif
299 PCSetFromOptions(mPCContext.PC_amg_A22);
300 PCSetUp(mPCContext.PC_amg_A22);
301}
302
303#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 1) //PETSc 3.1 or later
304PetscErrorCode PCBlockDiagonalApply(PC pc_object, Vec x, Vec y)
305{
306 void* pc_context;
307
308 PCShellGetContext(pc_object, &pc_context);
309#else
310PetscErrorCode PCBlockDiagonalApply(void* pc_context, Vec x, Vec y)
311{
312#endif
313
314 // Cast the context pointer to PCBlockDiagonalContext
316 assert(block_diag_context!=nullptr);
317
318 /*
319 * Scatter x = [x1 x2]'
320 */
321#ifdef TRACE_KSP
322 Timer::Reset();
323#endif
324
325 PetscVecTools::DoInterleavedVecScatter(x, block_diag_context->A11_scatter_ctx, block_diag_context->x1_subvector, block_diag_context->A22_scatter_ctx, block_diag_context->x2_subvector);
326
327#ifdef TRACE_KSP
328 block_diag_context->mScatterTime += Timer::GetElapsedTime();
329#endif
330
331 /*
332 * y1 = AMG(A11)*x1
333 * y2 = AMG(A22)*x2
334 */
335#ifdef TRACE_KSP
336 Timer::Reset();
337#endif
338 PCApply(block_diag_context->PC_amg_A11, block_diag_context->x1_subvector, block_diag_context->y1_subvector);
339#ifdef TRACE_KSP
340 block_diag_context->mA1PreconditionerTime += Timer::GetElapsedTime();
341#endif
342
343#ifdef TRACE_KSP
344 Timer::Reset();
345#endif
346 PCApply(block_diag_context->PC_amg_A22, block_diag_context->x2_subvector, block_diag_context->y2_subvector);
347#ifdef TRACE_KSP
348 block_diag_context->mA2PreconditionerTime += Timer::GetElapsedTime();
349#endif
350
351 /*
352 * Gather y = [y1 y2]'
353 */
354//PETSc-3.x.x or PETSc-2.3.3
355#ifdef TRACE_KSP
356 Timer::Reset();
357#endif
358
359 PetscVecTools::DoInterleavedVecGather(y, block_diag_context->A11_scatter_ctx, block_diag_context->y1_subvector, block_diag_context->A22_scatter_ctx, block_diag_context->y2_subvector);
360
361#ifdef TRACE_KSP
362 block_diag_context->mGatherTime += Timer::GetElapsedTime();
363#endif
364 return 0;
365}
#define TERMINATE(message)
#define PETSC_DESTROY_PARAM(x)
void PCBlockDiagonalCreate(KSP &rKspObject)
PCBlockDiagonal(KSP &rKspObject)
PCBlockDiagonalContext mPCContext
static void Destroy(Vec &rVec)
static bool AmMaster()
static Vec CreateVec(int size, int localSize=PETSC_DECIDE, bool ignoreOffProcEntries=true)
static void SetOption(const char *pOptionName, const char *pOptionValue)
static void DoInterleavedVecScatter(Vec interleavedVec, VecScatter firstVariableScatterContext, Vec firstVariableVec, VecScatter secondVariableScatterContext, Vec secondVariableVec)
static void DoInterleavedVecGather(Vec interleavedVec, VecScatter firstVariableScatterContext, Vec firstVariableVec, VecScatter secondVariableScatterContext, Vec secondVariableVec)
static void SetupInterleavedVectorScatterGather(Vec interleavedVec, VecScatter &rFirstVariableScatterContext, VecScatter &rSecondVariableScatterContext)
static double GetElapsedTime()
Definition Timer.cpp:54
static void Reset()
Definition Timer.cpp:44