Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
PCLDUFactorisation.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 "PCLDUFactorisation.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 PCLDUFactorisationCreate(rKspObject);
58}
59
60PCLDUFactorisation::~PCLDUFactorisation()
61{
62#ifdef TRACE_KSP
64 {
65 std::cout << " -- LDU factorisation 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 mExtraLAOperations: " << mPCContext.mExtraLAOperations << std::endl;
70 std::cout << "\t mGatherTime: " << mPCContext.mGatherTime << std::endl;
71 }
72#endif
73
77
80
87
90}
91
93{
94 KSPGetPC(rKspObject, &mPetscPCObject);
95
96 Mat system_matrix, dummy;
97#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=5)
98 KSPGetOperators(rKspObject, &system_matrix, &dummy);
99#else
100 MatStructure flag;
101 KSPGetOperators(rKspObject, &system_matrix, &dummy, &flag);
102#endif
103
104 PetscInt num_rows, num_columns;
105 MatGetSize(system_matrix, &num_rows, &num_columns);
106
107 PetscInt num_local_rows, num_local_columns;
108 MatGetLocalSize(system_matrix, &num_local_rows, &num_local_columns);
109
110 // Odd number of rows: impossible in Bidomain.
111 // Odd number of local rows: impossible if V_m and phi_e for each node are stored in the same processor.
112 if ((num_rows%2 != 0) || (num_local_rows%2 != 0))
113 {
114 TERMINATE("Wrong matrix parallel layout detected in PCLDUFactorisation."); // LCOV_EXCL_LINE
115 }
116
117 // Allocate memory
118 unsigned subvector_num_rows = num_rows/2;
119 unsigned subvector_local_rows = num_local_rows/2;
120 mPCContext.x1_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
121 mPCContext.x2_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
122 mPCContext.y1_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
123 mPCContext.y2_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
124 mPCContext.z = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
125 mPCContext.temp = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
126
127 // Create scatter contexts
128 {
129 // Needed by VecScatterCreate in order to find out parallel layout.
130 Vec dummy_vec = PetscTools::CreateVec(num_rows, num_local_rows);
131
133
134 PetscTools::Destroy(dummy_vec);
135 }
136
137 // Get matrix sublock A11
138 {
139 // Work out local row range for subblock A11 (same as x1 or y1)
140 PetscInt low, high, global_size;
141 VecGetOwnershipRange(mPCContext.x1_subvector, &low, &high);
142 VecGetSize(mPCContext.x1_subvector, &global_size);
143 assert(global_size == num_rows/2);
144
145 IS A11_local_rows;
146 IS A11_columns;
147 ISCreateStride(PETSC_COMM_WORLD, high-low, 2*low, 2, &A11_local_rows);
148 ISCreateStride(PETSC_COMM_WORLD, global_size, 0, 2, &A11_columns);
149
150#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 8) //PETSc 3.8 or later
151 MatCreateSubMatrix(system_matrix, A11_local_rows, A11_local_rows,
152 MAT_INITIAL_MATRIX, &mPCContext.A11_matrix_subblock);
153#else
154 MatGetSubMatrix(system_matrix, A11_local_rows, A11_local_rows,
155 MAT_INITIAL_MATRIX, &mPCContext.A11_matrix_subblock);
156#endif
157 ISDestroy(PETSC_DESTROY_PARAM(A11_local_rows));
158 ISDestroy(PETSC_DESTROY_PARAM(A11_columns));
159 }
160
161 // Get matrix sublock A22
162 {
163 // Work out local row range for subblock A22 (same as x2 or y2)
164 PetscInt low, high, global_size;
165 VecGetOwnershipRange(mPCContext.x2_subvector, &low, &high);
166 VecGetSize(mPCContext.x2_subvector, &global_size);
167 assert(global_size == num_rows/2);
168
169 IS A22_local_rows;
170 IS A22_columns;
171 ISCreateStride(PETSC_COMM_WORLD, high-low, 2*low+1, 2, &A22_local_rows);
172 ISCreateStride(PETSC_COMM_WORLD, global_size, 1, 2, &A22_columns);
173
174#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 8) //PETSc 3.8 or later
175 MatCreateSubMatrix(system_matrix, A22_local_rows, A22_local_rows,
176 MAT_INITIAL_MATRIX, &mPCContext.A22_matrix_subblock);
177#else
178 MatGetSubMatrix(system_matrix, A22_local_rows, A22_local_rows,
179 MAT_INITIAL_MATRIX, &mPCContext.A22_matrix_subblock);
180#endif
181
182 ISDestroy(PETSC_DESTROY_PARAM(A22_local_rows));
183 ISDestroy(PETSC_DESTROY_PARAM(A22_columns));
184 }
185
186 // Get matrix sublock B (the upper triangular one)
187 {
188 // Work out local row range for subblock B (same as A11, x1 or y1)
189 PetscInt low, high, global_size;
190 VecGetOwnershipRange(mPCContext.x1_subvector, &low, &high);
191 VecGetSize(mPCContext.x1_subvector, &global_size);
192 assert(global_size == num_rows/2);
193
194 IS B_local_rows;
195 IS B_columns;
196 ISCreateStride(PETSC_COMM_WORLD, high-low, 2*low, 2, &B_local_rows);
197
198 ISCreateStride(PETSC_COMM_WORLD, high-low, 2*low+1, 2, &B_columns);
199
200#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 8) //PETSc 3.8 or later
201 MatCreateSubMatrix(system_matrix, B_local_rows, B_columns,
202 MAT_INITIAL_MATRIX, &mPCContext.B_matrix_subblock);
203#else
204 MatGetSubMatrix(system_matrix, B_local_rows, B_columns,
205 MAT_INITIAL_MATRIX, &mPCContext.B_matrix_subblock);
206#endif
207
208 ISDestroy(PETSC_DESTROY_PARAM(B_local_rows));
209 ISDestroy(PETSC_DESTROY_PARAM(B_columns));
210 }
211
212 /*
213 * Experimental (#1082): in PP removing the mass matrix from the A22 block seems to work better.
214 * This is equivalent to do A22 = A22 + B in this implementation.
215 */
216//#if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
217// PetscScalar petsc_one = 1.0;
218// MatAXPY(&petsc_one, mPCContext.B_matrix_subblock, mPCContext.A22_matrix_subblock, DIFFERENT_NONZERO_PATTERN);
219// #else
220// MatAXPY(mPCContext.A22_matrix_subblock, 1.0, mPCContext.B_matrix_subblock, DIFFERENT_NONZERO_PATTERN);
221//#endif
222
223// // Shift the block
224//#if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
225// PetscScalar shift = -1e-8;
226// MatShift(&shift, mPCContext.A22_matrix_subblock);
227// #else
228// PetscScalar shift = -1e-8;
229// MatShift(mPCContext.A22_matrix_subblock, shift);
230//#endif
231
232 PCSetType(mPetscPCObject, PCSHELL);
233#if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
234 // Register PC context and call-back function
235 PCShellSetApply(mPetscPCObject, PCLDUFactorisationApply, (void*) &mPCContext);
236#else
237 // Register PC context so it gets passed to PCBlockDiagonalApply
238 PCShellSetContext(mPetscPCObject, &mPCContext);
239 // Register call-back function
240 PCShellSetApply(mPetscPCObject, PCLDUFactorisationApply);
241#endif
242}
243
245{
246 // These options will get read by PCSetFromOptions
247// PetscTools::SetOption("-pc_hypre_boomeramg_max_iter", "1");
248// PetscTools::SetOption("-pc_hypre_boomeramg_strong_threshold", "0.0");
249// PetscTools::SetOption("-pc_hypre_type", "boomeramg");
250
251 /*
252 * Set up preconditioner for block A11
253 */
254 PCCreate(PETSC_COMM_WORLD, &(mPCContext.PC_amg_A11));
255#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=5)
256 // Attempt to emulate SAME_PRECONDITIONER below
257 PCSetReusePreconditioner(mPCContext.PC_amg_A11, PETSC_TRUE);
259#else
261#endif
262
263 // Choose between the two following blocks in order to approximate inv(A11) with one AMG cycle
264 // or with an CG solve with high tolerance
266// PCSetType(mPCContext.PC_amg_A11, PCBJACOBI);
267
268 // We are expecting an error from PETSC on systems that don't have the hypre library, so suppress it
269 // in case it aborts
270 PetscPushErrorHandler(PetscIgnoreErrorHandler, nullptr);
271 PetscErrorCode pc_set_error = PCSetType(mPCContext.PC_amg_A11, PCHYPRE);
272 if (pc_set_error != 0)
273 {
274 WARNING("PETSc hypre preconditioning library is not installed"); // LCOV_EXCL_LINE
275 }
276 // Stop supressing error
277 PetscPopErrorHandler();
278
279#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<=5)
280 PetscTools::SetOption("-pc_hypre_type", "euclid");
281 PetscTools::SetOption("-pc_hypre_euclid_levels", "0");
282#else
283 /* euclid was removed in 3.6. Use the previously-commented-out alternative */
284 PetscTools::SetOption("-pc_hypre_type", "boomeramg");
285 PetscTools::SetOption("-pc_hypre_boomeramg_max_iter", "1");
286 PetscTools::SetOption("-pc_hypre_boomeramg_strong_threshold", "0.0");
287 PetscTools::SetOption("-pc_hypre_boomeramg_coarsen_type", "HMIS");
288#endif
289
291// PCSetType(mPCContext.PC_amg_A11, PCKSP);
292// KSP ksp1;
293// PCKSPGetKSP(mPCContext.PC_amg_A11,&ksp1);
294// KSPSetType(ksp1, KSPCG);
295// KSPSetTolerances(ksp1, 0.1, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
296//
297// PC prec1;
298// KSPGetPC(ksp1, &prec1);
299// PCSetType(prec1, PCBJACOBI);
300// PCSetFromOptions(prec1);
301// PCSetOperators(prec1, mPCContext.A11_matrix_subblock, mPCContext.A11_matrix_subblock, DIFFERENT_NONZERO_PATTERN);// SAME_PRECONDITIONER);
302// PCSetUp(prec1);
303//
304// KSPSetFromOptions(ksp1);
305// KSPSetUp(ksp1);
307
308 PCSetFromOptions(mPCContext.PC_amg_A11);
309 PCSetUp(mPCContext.PC_amg_A11);
310
311 /*
312 * Set up amg preconditioner for block A22
313 */
314 PCCreate(PETSC_COMM_WORLD, &(mPCContext.PC_amg_A22));
315#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=5)
316 // Attempt to emulate SAME_PRECONDITIONER below
317 PCSetReusePreconditioner(mPCContext.PC_amg_A22, PETSC_TRUE);
319#else
321#endif
322
323 // Choose between the two following blocks in order to approximate inv(A11) with one AMG cycle
324 // or with an CG solve with high tolerance
326 // We are expecting an error from PETSC on systems that don't have the hypre library, so suppress it
327 // in case it aborts
328 PetscPushErrorHandler(PetscIgnoreErrorHandler, nullptr);
329 PCSetType(mPCContext.PC_amg_A22, PCHYPRE);
330 // Stop supressing error
331 PetscPopErrorHandler();
332
333 PetscTools::SetOption("-pc_hypre_type", "boomeramg");
334 PetscTools::SetOption("-pc_hypre_boomeramg_max_iter", "1");
335 PetscTools::SetOption("-pc_hypre_boomeramg_strong_threshold", "0.0");
336 PetscTools::SetOption("-pc_hypre_boomeramg_coarsen_type", "HMIS");
337 // PetscTools::SetOption("-pc_hypre_boomeramg_interp_type","ext+i");
338
340// PCSetType(mPCContext.PC_amg_A22, PCKSP);
341// KSP ksp2;
342// PCKSPGetKSP(mPCContext.PC_amg_A22,&ksp2);
343// KSPSetType(ksp2, KSPCG);
344// KSPSetTolerances(ksp2, 0.1, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
345//
346// PC prec2;
347// KSPGetPC(ksp2, &prec2);
348// PCSetType(prec2, PCBJACOBI);
349// PCSetFromOptions(prec2);
350// PCSetOperators(prec2, mPCContext.A22_matrix_subblock, mPCContext.A22_matrix_subblock, DIFFERENT_NONZERO_PATTERN);// SAME_PRECONDITIONER);
351// PCSetUp(prec2);
352//
353// KSPSetFromOptions(ksp2);
354// KSPSetUp(ksp2);
356
357 PCSetFromOptions(mPCContext.PC_amg_A22);
358 PCSetUp(mPCContext.PC_amg_A22);
359}
360#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 1) //PETSc 3.1 or later
361PetscErrorCode PCLDUFactorisationApply(PC pc_object, Vec x, Vec y)
362{
363 void* pc_context;
364
365 PCShellGetContext(pc_object, &pc_context);
366#else
367PetscErrorCode PCLDUFactorisationApply(void* pc_context, Vec x, Vec y)
368{
369#endif
371
372 // Cast the pointer to a PC context to our defined type
374 assert(block_diag_context!=nullptr);
375
376 /*
377 * Split vector x into two. x = [x1 x2]'
378 */
379#ifdef TRACE_KSP
380 Timer::Reset();
381#endif
382
383 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);
384
385#ifdef TRACE_KSP
386 block_diag_context->mScatterTime += Timer::GetElapsedTime();
387#endif
388
389 /*
390 * Apply preconditioner: [y1 y2]' = inv(P)[x1 x2]'
391 *
392 * z = inv(A11)*x1
393 * y2 = inv(A22)*(x2 - B*z)
394 * y1 = z - inv(A11)(B*y2)
395 */
396#ifdef TRACE_KSP
397 Timer::Reset();
398#endif
399 //z = inv(A11)*x1
400 PCApply(block_diag_context->PC_amg_A11, block_diag_context->x1_subvector, block_diag_context->z);
401#ifdef TRACE_KSP
402 block_diag_context->mA1PreconditionerTime += Timer::GetElapsedTime();
403#endif
404
405 //y2 = inv(A22)*(x2 - B*z)
406#ifdef TRACE_KSP
407 Timer::Reset();
408#endif
409 MatMult(block_diag_context->B_matrix_subblock,block_diag_context->z,block_diag_context->temp); //temp = B*z
410 double minus_one = -1.0;
411#if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
412 VecAYPX(&minus_one, block_diag_context->x2_subvector, block_diag_context->temp); // temp <-- x2 - temp
413#else
414 VecAYPX(block_diag_context->temp, minus_one, block_diag_context->x2_subvector); // temp <-- x2 - temp
415#endif
416#ifdef TRACE_KSP
417 block_diag_context->mExtraLAOperations += Timer::GetElapsedTime();
418#endif
419
420
421#ifdef TRACE_KSP
422 Timer::Reset();
423#endif
424 PCApply(block_diag_context->PC_amg_A22, block_diag_context->temp, block_diag_context->y2_subvector); // y2 = inv(A22)*temp
425#ifdef TRACE_KSP
426 block_diag_context->mA2PreconditionerTime += Timer::GetElapsedTime();
427#endif
428
429 // y1 = z - inv(A11)(B*y2)
430#ifdef TRACE_KSP
431 Timer::Reset();
432#endif
433 MatMult(block_diag_context->B_matrix_subblock,block_diag_context->y2_subvector,block_diag_context->temp); //temp = B*y2
434#ifdef TRACE_KSP
435 block_diag_context->mExtraLAOperations += Timer::GetElapsedTime();
436#endif
437#ifdef TRACE_KSP
438 Timer::Reset();
439#endif
440 PCApply(block_diag_context->PC_amg_A11, block_diag_context->temp, block_diag_context->y1_subvector); // y1 = inv(A11)*temp
441#ifdef TRACE_KSP
442 block_diag_context->mA1PreconditionerTime += Timer::GetElapsedTime();
443#endif
444
445#ifdef TRACE_KSP
446 Timer::Reset();
447#endif
448#if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
449 VecAYPX(&minus_one, block_diag_context->z, block_diag_context->y1_subvector); // y1 <-- z - y1
450#else
451 VecAYPX(block_diag_context->y1_subvector, minus_one, block_diag_context->z); // y1 <-- z - y1
452#endif
453#ifdef TRACE_KSP
454 block_diag_context->mExtraLAOperations += Timer::GetElapsedTime();
455#endif
456
457 /*
458 * Gather vectors y1 and y2. y = [y1 y2]'
459 */
460#ifdef TRACE_KSP
461 Timer::Reset();
462#endif
463
464 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);
465
466#ifdef TRACE_KSP
467 block_diag_context->mGatherTime += Timer::GetElapsedTime();
468#endif
469
470 return 0;
471}
#define TERMINATE(message)
#define PETSC_DESTROY_PARAM(x)
PCLDUFactorisationContext mPCContext
void PCLDUFactorisationCreate(KSP &rKspObject)
PCLDUFactorisation(KSP &rKspObject)
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