Chaste  Release::2018.1
PCBlockDiagonal.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, 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 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF 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 
47 PCBlockDiagonal::PCBlockDiagonal(KSP& rKspObject)
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);
57  PCBlockDiagonalSetUp();
58 }
59 
60 PCBlockDiagonal::~PCBlockDiagonal()
61 {
62 #ifdef TRACE_KSP
63  if (PetscTools::AmMaster())
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 
73  PetscTools::Destroy(mPCContext.A11_matrix_subblock);
74  PetscTools::Destroy(mPCContext.A22_matrix_subblock);
75 
76  PCDestroy(PETSC_DESTROY_PARAM(mPCContext.PC_amg_A11));
77  PCDestroy(PETSC_DESTROY_PARAM(mPCContext.PC_amg_A22));
78 
79  PetscTools::Destroy(mPCContext.x1_subvector);
80  PetscTools::Destroy(mPCContext.y1_subvector);
81 
82  PetscTools::Destroy(mPCContext.x2_subvector);
83  PetscTools::Destroy(mPCContext.y2_subvector);
84 
85  VecScatterDestroy(PETSC_DESTROY_PARAM(mPCContext.A11_scatter_ctx));
86  VecScatterDestroy(PETSC_DESTROY_PARAM(mPCContext.A22_scatter_ctx));
87 }
88 
89 void PCBlockDiagonal::PCBlockDiagonalCreate(KSP& rKspObject)
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 
128  PetscVecTools::SetupInterleavedVectorScatterGather(dummy_vec, mPCContext.A11_scatter_ctx, mPCContext.A22_scatter_ctx);
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 
198 void PCBlockDiagonal::PCBlockDiagonalSetUp()
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");
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);
235  PCSetOperators(mPCContext.PC_amg_A11, mPCContext.A11_matrix_subblock, mPCContext.A11_matrix_subblock);
236 #else
237  PCSetOperators(mPCContext.PC_amg_A11, mPCContext.A11_matrix_subblock, mPCContext.A11_matrix_subblock, SAME_PRECONDITIONER);
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);
295  PCSetOperators(mPCContext.PC_amg_A22, mPCContext.A22_matrix_subblock, mPCContext.A22_matrix_subblock);
296 #else
297  PCSetOperators(mPCContext.PC_amg_A22, mPCContext.A22_matrix_subblock, mPCContext.A22_matrix_subblock, SAME_PRECONDITIONER);
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
304 PetscErrorCode PCBlockDiagonalApply(PC pc_object, Vec x, Vec y)
305 {
306  void* pc_context;
307 
308  PCShellGetContext(pc_object, &pc_context);
309 #else
310 PetscErrorCode PCBlockDiagonalApply(void* pc_context, Vec x, Vec y)
311 {
312 #endif
313 
314  // Cast the context pointer to PCBlockDiagonalContext
315  PCBlockDiagonal::PCBlockDiagonalContext* block_diag_context = (PCBlockDiagonal::PCBlockDiagonalContext*) pc_context;
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 }
PetscVecTools::SetupInterleavedVectorScatterGather
static void SetupInterleavedVectorScatterGather(Vec interleavedVec, VecScatter &rFirstVariableScatterContext, VecScatter &rSecondVariableScatterContext)
Definition: PetscVecTools.cpp:154
PCBlockDiagonal::PCBlockDiagonalCreate
void PCBlockDiagonalCreate(KSP &rKspObject)
Definition: PCBlockDiagonal.cpp:89
PetscTools::SetOption
static void SetOption(const char *pOptionName, const char *pOptionValue)
Definition: PetscTools.hpp:384
PETSC_DESTROY_PARAM
#define PETSC_DESTROY_PARAM(x)
Definition: PetscTools.hpp:70
PetscTools::AmMaster
static bool AmMaster()
Definition: PetscTools.cpp:120
PetscTools::CreateVec
static Vec CreateVec(int size, int localSize=PETSC_DECIDE, bool ignoreOffProcEntries=true)
Definition: PetscTools.cpp:214
TERMINATE
#define TERMINATE(message)
Definition: Exception.hpp:168
Timer::GetElapsedTime
static double GetElapsedTime()
Definition: Timer.cpp:54
PetscVecTools::DoInterleavedVecScatter
static void DoInterleavedVecScatter(Vec interleavedVec, VecScatter firstVariableScatterContext, Vec firstVariableVec, VecScatter secondVariableScatterContext, Vec secondVariableVec)
Definition: PetscVecTools.cpp:184
PCBlockDiagonal::mPCContext
PCBlockDiagonalContext mPCContext
Definition: PCBlockDiagonal.hpp:112
PetscTools::Destroy
static void Destroy(Vec &rVec)
Definition: PetscTools.hpp:352
PetscVecTools::DoInterleavedVecGather
static void DoInterleavedVecGather(Vec interleavedVec, VecScatter firstVariableScatterContext, Vec firstVariableVec, VecScatter secondVariableScatterContext, Vec secondVariableVec)
Definition: PetscVecTools.cpp:254
Timer::Reset
static void Reset()
Definition: Timer.cpp:44
PCBlockDiagonal::PCBlockDiagonal
PCBlockDiagonal(KSP &rKspObject)
Definition: PCBlockDiagonal.cpp:47