Chaste  Release::2018.1
PCTwoLevelsBlockDiagonal.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 "PCTwoLevelsBlockDiagonal.hpp"
37 #include "Exception.hpp"
38 
39 #include <iostream>
40 
41 PCTwoLevelsBlockDiagonal::PCTwoLevelsBlockDiagonal(KSP& rKspObject, std::vector<PetscInt>& rBathNodes)
42 {
43  PCTwoLevelsBlockDiagonalCreate(rKspObject, rBathNodes);
45 }
46 
48 {
52 
56 
59 
62 
65 
69 }
70 
71 void 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
305 PetscErrorCode PCTwoLevelsBlockDiagonalApply(PC pc_object, Vec x, Vec y)
306 {
307  void* pc_context;
308 
309  PCShellGetContext(pc_object, &pc_context);
310 #else
311 PetscErrorCode 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 }
static void SetOption(const char *pOptionName, const char *pOptionValue)
Definition: PetscTools.hpp:384
#define PETSC_DESTROY_PARAM(x)
Definition: PetscTools.hpp:70
static Vec CreateVec(int size, int localSize=PETSC_DECIDE, bool ignoreOffProcEntries=true)
Definition: PetscTools.cpp:214
#define TERMINATE(message)
Definition: Exception.hpp:168
PCTwoLevelsBlockDiagonalContext mPCContext
static bool IsSequential()
Definition: PetscTools.cpp:91
void PCTwoLevelsBlockDiagonalCreate(KSP &rKspObject, std::vector< PetscInt > &rBathNodes)
static void Destroy(Vec &rVec)
Definition: PetscTools.hpp:352
PCTwoLevelsBlockDiagonal(KSP &rKspObject, std::vector< PetscInt > &rBathNodes)