PCTwoLevelsBlockDiagonal.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 #include "PCTwoLevelsBlockDiagonal.hpp"
00030 #include "Exception.hpp"
00031 
00032 #include <iostream>
00033 
00034 PCTwoLevelsBlockDiagonal::PCTwoLevelsBlockDiagonal(KSP& rKspObject, std::vector<PetscInt>& rBathNodes)
00035 {
00036     PCTwoLevelsBlockDiagonalCreate(rKspObject, rBathNodes);
00037     PCTwoLevelsBlockDiagonalSetUp();
00038 }
00039 
00040 PCTwoLevelsBlockDiagonal::~PCTwoLevelsBlockDiagonal()
00041 {
00042     MatDestroy(mPCContext.A11_matrix_subblock);
00043     MatDestroy(mPCContext.A22_B1_matrix_subblock);
00044     MatDestroy(mPCContext.A22_B2_matrix_subblock);
00045 
00046     PCDestroy(mPCContext.PC_amg_A11);
00047     PCDestroy(mPCContext.PC_amg_A22_B1);
00048     PCDestroy(mPCContext.PC_amg_A22_B2);
00049 
00050     VecDestroy(mPCContext.x1_subvector);
00051     VecDestroy(mPCContext.y1_subvector);
00052 
00053     VecDestroy(mPCContext.x21_subvector);
00054     VecDestroy(mPCContext.y21_subvector);
00055 
00056     VecDestroy(mPCContext.x22_subvector);
00057     VecDestroy(mPCContext.y22_subvector);
00058     
00059     VecScatterDestroy(mPCContext.A11_scatter_ctx);
00060     VecScatterDestroy(mPCContext.A22_B1_scatter_ctx);    
00061     VecScatterDestroy(mPCContext.A22_B2_scatter_ctx);    
00062 }
00063 
00064 void PCTwoLevelsBlockDiagonal::PCTwoLevelsBlockDiagonalCreate(KSP& rKspObject, std::vector<PetscInt>& rBathNodes)
00065 {
00066     KSPGetPC(rKspObject, &mPetscPCObject);
00067 
00068     Mat system_matrix, dummy;
00069     MatStructure flag;
00070     KSPGetOperators(rKspObject, &system_matrix, &dummy, &flag);
00071 
00072     PetscInt num_rows, num_columns;
00073     MatGetSize(system_matrix, &num_rows, &num_columns);
00074     assert(num_rows==num_columns);
00075 
00076     PetscInt num_local_rows, num_local_columns;
00077     MatGetLocalSize(system_matrix, &num_local_rows, &num_local_columns);
00078 
00079     // odd number of rows: impossible in Bidomain.
00080     // odd number of local rows: impossible if V_m and phi_e for each node are stored in the same processor.
00081     if ((num_rows%2 != 0) || (num_local_rows%2 != 0))
00082     {
00083         TERMINATE("Wrong matrix parallel layout detected in PCLDUFactorisation.");
00084     }
00085 
00086     // Allocate memory     
00087     unsigned subvector_num_rows = num_rows/2;    
00088     unsigned subvector_local_rows = num_local_rows/2;
00089     
00090     unsigned subvector_num_rows_tissue = subvector_num_rows - rBathNodes.size();    
00091     unsigned subvector_local_rows_tissue = subvector_num_rows_tissue; 
00092 
00093     unsigned subvector_num_rows_bath = rBathNodes.size();    
00094     unsigned subvector_local_rows_bath = subvector_num_rows_bath; 
00095     
00096     assert(PetscTools::IsSequential());
00097     
00098     mPCContext.x1_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
00099     mPCContext.x21_subvector = PetscTools::CreateVec(subvector_num_rows_tissue, subvector_local_rows_tissue);
00100     mPCContext.x22_subvector = PetscTools::CreateVec(subvector_num_rows_bath, subvector_local_rows_bath);
00101     mPCContext.y1_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
00102     mPCContext.y21_subvector = PetscTools::CreateVec(subvector_num_rows_tissue, subvector_local_rows_tissue);
00103     mPCContext.y22_subvector = PetscTools::CreateVec(subvector_num_rows_bath, subvector_local_rows_bath);
00104 
00105     /*
00106      * Define IS objects that will be used throughout the method.
00107      */
00108     IS A11_all_rows;
00109     ISCreateStride(PETSC_COMM_WORLD, num_rows/2, 0, 2, &A11_all_rows);     
00110      
00111     IS A22_all_rows;
00112     PetscInt A22_size;
00113     VecGetSize(mPCContext.x1_subvector, &A22_size);
00115     ISCreateStride(PETSC_COMM_WORLD, A22_size, 1, 2, &A22_all_rows);
00116            
00117     IS A22_bath_rows;
00118     PetscInt* phi_e_bath_rows = new PetscInt[rBathNodes.size()];
00119     for (unsigned index=0; index<rBathNodes.size(); index++)
00120     {
00121         phi_e_bath_rows[index] = 2*rBathNodes[index] + 1;
00122     }
00123     ISCreateGeneralWithArray(PETSC_COMM_WORLD, rBathNodes.size(), phi_e_bath_rows, &A22_bath_rows);        
00124     
00125     IS A22_tissue_rows;
00126     ISDifference(A22_all_rows, A22_bath_rows, &A22_tissue_rows);            
00127 
00128     // Create scatter contexts
00129     {
00130         // Needed by VecScatterCreate in order to find out parallel layout.
00131         Vec dummy_vec = PetscTools::CreateVec(num_rows, num_local_rows);
00132 
00134         IS& A11_rows=A11_all_rows;
00135         IS& A22_B1_rows=A22_tissue_rows;
00136         IS& A22_B2_rows=A22_bath_rows;        
00137         
00138         IS all_vector;    
00139         ISCreateStride(PETSC_COMM_WORLD, num_rows/2, 0, 1, &all_vector);
00140         
00141         IS tissue_vector;    
00142         ISCreateStride(PETSC_COMM_WORLD, (num_rows/2)-rBathNodes.size(), 0, 1, &tissue_vector);
00143 
00144         IS bath_vector;    
00145         ISCreateStride(PETSC_COMM_WORLD, rBathNodes.size(), 0, 1, &bath_vector);
00146         
00147         VecScatterCreate(dummy_vec, A11_rows, mPCContext.x1_subvector, all_vector, &mPCContext.A11_scatter_ctx);    
00148         VecScatterCreate(dummy_vec, A22_B1_rows, mPCContext.x21_subvector, tissue_vector, &mPCContext.A22_B1_scatter_ctx);
00149         VecScatterCreate(dummy_vec, A22_B2_rows, mPCContext.x22_subvector, bath_vector, &mPCContext.A22_B2_scatter_ctx);
00150         
00151         ISDestroy(all_vector);
00152         ISDestroy(tissue_vector);
00153         ISDestroy(bath_vector);
00154         
00155         VecDestroy(dummy_vec);        
00156     }
00157             
00158     // Get matrix sublock A11        
00159     {
00160         // Work out local row range for subblock A11 (same as x1 or y1) 
00161         PetscInt low, high, global_size;
00162         VecGetOwnershipRange(mPCContext.x1_subvector, &low, &high);
00163         VecGetSize(mPCContext.x1_subvector, &global_size);        
00164         assert(global_size == num_rows/2);       
00165         
00166         IS A11_local_rows;
00167         IS& A11_columns=A11_all_rows;
00168         ISCreateStride(PETSC_COMM_WORLD, high-low, 2*low, 2, &A11_local_rows); 
00169     
00170 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR == 1) //PETSc 3.1
00171         MatGetSubMatrix(system_matrix, A11_local_rows, A11_columns,
00172             MAT_INITIAL_MATRIX, &mPCContext.A11_matrix_subblock);
00173 #else
00174         MatGetSubMatrix(system_matrix, A11_local_rows, A11_columns, PETSC_DECIDE, 
00175             MAT_INITIAL_MATRIX, &mPCContext.A11_matrix_subblock);
00176 #endif
00177     
00178         ISDestroy(A11_local_rows);
00179     }
00180 
00181     // Get matrix sublock A22_B1
00182     {
00183 //        // Work out local row range for subblock A22 (same as x2 or y2) 
00184 //        PetscInt low, high, global_size;
00185 //        VecGetOwnershipRange(mPCContext.x21_subvector, &low, &high);
00186 //        VecGetSize(mPCContext.x21_subvector, &global_size);        
00187 //        assert(global_size == (num_rows/2) - (PetscInt) rBathNodes.size());       
00188         
00189         assert(PetscTools::IsSequential());
00190         IS& A22_B1_local_rows = A22_tissue_rows; // wrong in parallel, need to give local rows
00191         IS& A22_B1_columns = A22_tissue_rows;
00192                 
00193 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR == 1) //PETSc 3.1
00194         MatGetSubMatrix(system_matrix, A22_B1_local_rows, A22_B1_columns,
00195             MAT_INITIAL_MATRIX, &mPCContext.A22_B1_matrix_subblock);
00196 #else
00197         MatGetSubMatrix(system_matrix, A22_B1_local_rows, A22_B1_columns, PETSC_DECIDE, 
00198             MAT_INITIAL_MATRIX, &mPCContext.A22_B1_matrix_subblock);
00199 #endif
00200     
00201     }
00202 
00203     // Get matrix sublock A22_B2
00204     {
00205 //        // Work out local row range for subblock A22 (same as x2 or y2) 
00206 //        PetscInt low, high, global_size;
00207 //        VecGetOwnershipRange(mPCContext.x21_subvector, &low, &high);
00208 //        VecGetSize(mPCContext.x21_subvector, &global_size);        
00209 //        assert(global_size == (num_rows/2) - (PetscInt) rBathNodes.size());       
00210         
00211         assert(PetscTools::IsSequential());
00212         IS& A22_B2_local_rows = A22_bath_rows; // wrong in parallel, need to give local rows
00213         IS& A22_B2_columns = A22_bath_rows;
00214                 
00215 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR == 1) //PETSc 3.1
00216         MatGetSubMatrix(system_matrix, A22_B2_local_rows, A22_B2_columns,
00217             MAT_INITIAL_MATRIX, &mPCContext.A22_B2_matrix_subblock);
00218 #else
00219         MatGetSubMatrix(system_matrix, A22_B2_local_rows, A22_B2_columns, PETSC_DECIDE, 
00220             MAT_INITIAL_MATRIX, &mPCContext.A22_B2_matrix_subblock);
00221 #endif
00222     
00223     }
00224     
00225     ISDestroy(A11_all_rows);
00226     ISDestroy(A22_all_rows);
00227     ISDestroy(A22_bath_rows);
00228     delete[] phi_e_bath_rows;
00229     ISDestroy(A22_tissue_rows);
00230 
00231     // Register call-back function and its context
00232     PCSetType(mPetscPCObject, PCSHELL);
00233 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
00234     PCShellSetApply(mPetscPCObject, PCTwoLevelsBlockDiagonalApply, (void*) &mPCContext);
00235 #else
00236     // Register PC context so it gets passed to PCTwoLevelsBlockDiagonalApply
00237     PCShellSetContext(mPetscPCObject, &mPCContext);
00238     // Register call-back function
00239     PCShellSetApply(mPetscPCObject, PCTwoLevelsBlockDiagonalApply);
00240 #endif
00241 }
00242 
00243 void PCTwoLevelsBlockDiagonal::PCTwoLevelsBlockDiagonalSetUp()
00244 {
00245     // These options will get read by PCSetFromOptions
00246     PetscOptionsSetValue("-pc_hypre_boomeramg_max_iter", "1");
00247     PetscOptionsSetValue("-pc_hypre_boomeramg_strong_threshold", "0.0");
00248     PetscOptionsSetValue("-pc_hypre_type", "boomeramg");
00249 
00250     // Set up amg preconditioner for block A11
00251     PCCreate(PETSC_COMM_WORLD, &(mPCContext.PC_amg_A11));
00252     PCSetType(mPCContext.PC_amg_A11, PCBJACOBI);
00253     PCSetOperators(mPCContext.PC_amg_A11, mPCContext.A11_matrix_subblock, mPCContext.A11_matrix_subblock, DIFFERENT_NONZERO_PATTERN);//   SAME_PRECONDITIONER);
00254     PCSetFromOptions(mPCContext.PC_amg_A11);
00255     PCSetUp(mPCContext.PC_amg_A11);
00256 
00257     // Set up amg preconditioner for block A22_B1
00258     PCCreate(PETSC_COMM_WORLD, &(mPCContext.PC_amg_A22_B1));
00259     PCSetType(mPCContext.PC_amg_A22_B1, PCBJACOBI);
00260     PCSetOperators(mPCContext.PC_amg_A22_B1, mPCContext.A22_B1_matrix_subblock, mPCContext.A22_B1_matrix_subblock, DIFFERENT_NONZERO_PATTERN);//   SAME_PRECONDITIONER);
00261     PCSetFromOptions(mPCContext.PC_amg_A22_B1);
00262     PCSetUp(mPCContext.PC_amg_A22_B1);
00263 
00264     // Set up amg preconditioner for block A22_B2
00265     PCCreate(PETSC_COMM_WORLD, &(mPCContext.PC_amg_A22_B2));
00266     PCSetType(mPCContext.PC_amg_A22_B2, PCHYPRE);
00267     //PCHYPRESetType(mPCContext.PC_amg_A22_B2, "boomeramg");
00268     PetscOptionsSetValue("-pc_hypre_type", "boomeramg");
00269 
00270     PetscOptionsSetValue("-pc_hypre_boomeramg_max_iter", "1");
00271     PetscOptionsSetValue("-pc_hypre_boomeramg_strong_threshold", "0.0");
00272     PetscOptionsSetValue("-pc_hypre_boomeramg_coarsen_type", "HMIS");
00273 
00274     PCSetOperators(mPCContext.PC_amg_A22_B2, mPCContext.A22_B2_matrix_subblock, mPCContext.A22_B2_matrix_subblock, DIFFERENT_NONZERO_PATTERN);//   SAME_PRECONDITIONER);
00275     PCSetFromOptions(mPCContext.PC_amg_A22_B2);
00276     PCSetUp(mPCContext.PC_amg_A22_B2);
00277 }
00278 
00279 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR == 1) //PETSc 3.1
00280 PetscErrorCode PCTwoLevelsBlockDiagonalApply(PC pc_object, Vec x, Vec y)
00281 {
00282   void* pc_context;
00283 
00284   PCShellGetContext(pc_object, &pc_context);   
00285 #else
00286 PetscErrorCode PCTwoLevelsBlockDiagonalApply(void* pc_context, Vec x, Vec y)
00287 {
00288 #endif
00289 
00290     // Cast the context pointer to PCTwoLevelsBlockDiagonalContext
00291     PCTwoLevelsBlockDiagonal::PCTwoLevelsBlockDiagonalContext* block_diag_context = (PCTwoLevelsBlockDiagonal::PCTwoLevelsBlockDiagonalContext*) pc_context;
00292     assert(block_diag_context!=NULL);
00293 
00294     /*
00295      * Scatter x = [x1 x21 x22]'
00296      */
00297 //PETSc-3.x.x or PETSc-2.3.3
00298 #if ( (PETSC_VERSION_MAJOR == 3) || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 3 && PETSC_VERSION_SUBMINOR == 3)) //2.3.3 or 3.x.x
00299     VecScatterBegin(block_diag_context->A11_scatter_ctx, x, block_diag_context->x1_subvector, INSERT_VALUES, SCATTER_FORWARD);
00300     VecScatterEnd(block_diag_context->A11_scatter_ctx, x, block_diag_context->x1_subvector, INSERT_VALUES, SCATTER_FORWARD);
00301 #else
00302     VecScatterBegin(x, block_diag_context->x1_subvector, INSERT_VALUES, SCATTER_FORWARD, block_diag_context->A11_scatter_ctx);
00303     VecScatterEnd(x, block_diag_context->x1_subvector, INSERT_VALUES, SCATTER_FORWARD, block_diag_context->A11_scatter_ctx);
00304 #endif
00305 
00306 //PETSc-3.x.x or PETSc-2.3.3
00307 #if ( (PETSC_VERSION_MAJOR == 3) || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 3 && PETSC_VERSION_SUBMINOR == 3)) //2.3.3 or 3.x.x
00308     VecScatterBegin(block_diag_context->A22_B1_scatter_ctx, x, block_diag_context->x21_subvector, INSERT_VALUES, SCATTER_FORWARD);
00309     VecScatterEnd(block_diag_context->A22_B1_scatter_ctx, x, block_diag_context->x21_subvector, INSERT_VALUES, SCATTER_FORWARD);
00310 #else
00311     VecScatterBegin(x, block_diag_context->x21_subvector, INSERT_VALUES, SCATTER_FORWARD, block_diag_context->A22_B1_scatter_ctx);
00312     VecScatterEnd(x, block_diag_context->x21_subvector, INSERT_VALUES, SCATTER_FORWARD, block_diag_context->A22_B1_scatter_ctx);
00313 #endif
00314 
00315 //PETSc-3.x.x or PETSc-2.3.3
00316 #if ( (PETSC_VERSION_MAJOR == 3) || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 3 && PETSC_VERSION_SUBMINOR == 3)) //2.3.3 or 3.x.x
00317     VecScatterBegin(block_diag_context->A22_B2_scatter_ctx, x, block_diag_context->x22_subvector, INSERT_VALUES, SCATTER_FORWARD);
00318     VecScatterEnd(block_diag_context->A22_B2_scatter_ctx, x, block_diag_context->x22_subvector, INSERT_VALUES, SCATTER_FORWARD);
00319 #else
00320     VecScatterBegin(x, block_diag_context->x22_subvector, INSERT_VALUES, SCATTER_FORWARD, block_diag_context->A22_B2_scatter_ctx);
00321     VecScatterEnd(x, block_diag_context->x22_subvector, INSERT_VALUES, SCATTER_FORWARD, block_diag_context->A22_B2_scatter_ctx);
00322 #endif
00323 
00324     /*
00325      *  y1  = ILU(A11)*x1
00326      *  y21 = ILU(A22)*x21
00327      *  y22 = AMG(A22)*x22
00328      */
00329     PCApply(block_diag_context->PC_amg_A11, block_diag_context->x1_subvector, block_diag_context->y1_subvector);
00330     PCApply(block_diag_context->PC_amg_A22_B1, block_diag_context->x21_subvector, block_diag_context->y21_subvector);
00331     PCApply(block_diag_context->PC_amg_A22_B2, block_diag_context->x22_subvector, block_diag_context->y22_subvector);
00332 
00333     /*
00334      * Gather y = [y1 y21 y22]'
00335      */
00336 //PETSc-3.x.x or PETSc-2.3.3
00337 #if ( (PETSC_VERSION_MAJOR == 3) || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 3 && PETSC_VERSION_SUBMINOR == 3)) //2.3.3 or 3.x.x
00338     VecScatterBegin(block_diag_context->A11_scatter_ctx, block_diag_context->y1_subvector, y, INSERT_VALUES, SCATTER_REVERSE);
00339     VecScatterEnd(block_diag_context->A11_scatter_ctx, block_diag_context->y1_subvector, y, INSERT_VALUES, SCATTER_REVERSE);
00340 #else
00341     VecScatterBegin(block_diag_context->y1_subvector, y, INSERT_VALUES, SCATTER_REVERSE, block_diag_context->A11_scatter_ctx);
00342     VecScatterEnd(block_diag_context->y1_subvector, y, INSERT_VALUES, SCATTER_REVERSE, block_diag_context->A11_scatter_ctx);
00343 #endif
00344 
00345 //PETSc-3.x.x or PETSc-2.3.3
00346 #if ( (PETSC_VERSION_MAJOR == 3) || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 3 && PETSC_VERSION_SUBMINOR == 3)) //2.3.3 or 3.x.x
00347     VecScatterBegin(block_diag_context->A22_B1_scatter_ctx, block_diag_context->y21_subvector, y, INSERT_VALUES, SCATTER_REVERSE);
00348     VecScatterEnd(block_diag_context->A22_B1_scatter_ctx, block_diag_context->y21_subvector, y, INSERT_VALUES, SCATTER_REVERSE);
00349 #else
00350     VecScatterBegin(block_diag_context->y21_subvector, y, INSERT_VALUES, SCATTER_REVERSE, block_diag_context->A22_B1_scatter_ctx);
00351     VecScatterEnd(block_diag_context->y21_subvector, y, INSERT_VALUES, SCATTER_REVERSE, block_diag_context->A22_B1_scatter_ctx);
00352 #endif
00353 
00354 //PETSc-3.x.x or PETSc-2.3.3
00355 #if ( (PETSC_VERSION_MAJOR == 3) || (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 3 && PETSC_VERSION_SUBMINOR == 3)) //2.3.3 or 3.x.x
00356     VecScatterBegin(block_diag_context->A22_B2_scatter_ctx, block_diag_context->y22_subvector, y, INSERT_VALUES, SCATTER_REVERSE);
00357     VecScatterEnd(block_diag_context->A22_B2_scatter_ctx, block_diag_context->y22_subvector, y, INSERT_VALUES, SCATTER_REVERSE);
00358 #else
00359     VecScatterBegin(block_diag_context->y22_subvector, y, INSERT_VALUES, SCATTER_REVERSE, block_diag_context->A22_B2_scatter_ctx);
00360     VecScatterEnd(block_diag_context->y22_subvector, y, INSERT_VALUES, SCATTER_REVERSE, block_diag_context->A22_B2_scatter_ctx);
00361 #endif
00362 
00363     return 0;
00364 }

Generated on Tue May 31 14:31:42 2011 for Chaste by  doxygen 1.5.5