1 #include <petscksp.h> 2 #include <petscpc.h> 3 #include <petscviewer.h> 4 5 typedef struct { 6 PetscInt num_levels; 7 PetscInt *n_per_level; 8 Mat stiff; 9 Mat *ProlongationOps; 10 PetscBT *CFMarkers; 11 KSP kspHypre; 12 } *DataCompression; 13 14 PetscErrorCode Create1dLaplacian(PetscInt, Mat *); 15 PetscErrorCode DataCompExportMats(DataCompression); 16 PetscErrorCode DataCompDestroy(DataCompression); 17 18 int main(int Argc, char **Args) 19 { 20 PetscInt n_nodes = 33; 21 Vec x, b; 22 PC pcHypre; 23 DataCompression data_comp; 24 25 PetscFunctionBeginUser; 26 PetscCall(PetscInitialize(&Argc, &Args, NULL, NULL)); 27 28 PetscCall(PetscNew(&data_comp)); 29 30 // Creating stiffness matrix 31 PetscCall(Create1dLaplacian(n_nodes, &data_comp->stiff)); 32 PetscCall(PetscObjectSetName((PetscObject)data_comp->stiff, "Stiffness")); 33 34 // Set-up BoomerAMG PC to get Prolongation Operators and Coarse/Fine splittings 35 PetscCall(KSPCreate(PETSC_COMM_WORLD, &data_comp->kspHypre)); 36 PetscCall(KSPSetType(data_comp->kspHypre, KSPRICHARDSON)); 37 PetscCall(KSPGetPC(data_comp->kspHypre, &pcHypre)); 38 PetscCall(PCSetType(pcHypre, PCHYPRE)); 39 PetscCall(PCHYPRESetType(pcHypre, "boomeramg")); 40 PetscCall(PCSetFromOptions(pcHypre)); 41 PetscCall(PCSetOperators(pcHypre, data_comp->stiff, data_comp->stiff)); 42 PetscCall(PCSetUp(pcHypre)); 43 44 PetscCall(MatCreateVecs(data_comp->stiff, &x, &b)); 45 PetscCall(PCApply(pcHypre, x, b)); 46 PetscCall(VecDestroy(&x)); 47 PetscCall(VecDestroy(&b)); 48 49 //Viewing the PC and Extracting the Prolongation Operators and CFMarkers 50 PetscCall(PCView(pcHypre, NULL)); 51 PetscCall(PCGetInterpolations(pcHypre, &data_comp->num_levels, &data_comp->ProlongationOps)); 52 PetscCall(PCHYPREGetCFMarkers(pcHypre, &data_comp->n_per_level, &data_comp->CFMarkers)); 53 54 PetscCall(DataCompExportMats(data_comp)); 55 PetscCall(DataCompDestroy(data_comp)); 56 PetscCall(PetscFinalize()); 57 return 0; 58 } 59 60 PetscErrorCode Create1dLaplacian(PetscInt n, Mat *mat) 61 { 62 PetscFunctionBeginUser; 63 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 3, NULL, mat)); 64 PetscCall(MatSetValue(*mat, n - 1, n - 1, 2.0, INSERT_VALUES)); 65 for (PetscInt i = 0; i < n - 1; i++) { 66 PetscCall(MatSetValue(*mat, i, i, 2.0, INSERT_VALUES)); 67 PetscCall(MatSetValue(*mat, i + 1, i, -1.0, INSERT_VALUES)); 68 PetscCall(MatSetValue(*mat, i, i + 1, -1.0, INSERT_VALUES)); 69 } 70 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); 71 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); 72 PetscFunctionReturn(PETSC_SUCCESS); 73 } 74 75 PetscErrorCode DataCompExportMats(DataCompression data_comp) 76 { 77 PetscFunctionBeginUser; 78 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Num levels: %" PetscInt_FMT "\n", data_comp->num_levels)); 79 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " -- Nodes per level --\n")); 80 for (PetscInt i = 0; i < data_comp->num_levels; i++) { PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Level %" PetscInt_FMT ": %" PetscInt_FMT "\n", i, data_comp->n_per_level[i])); } 81 82 for (PetscInt i = 0; i < data_comp->num_levels - 1; i++) { 83 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Prolongation Operator - Level %" PetscInt_FMT "\n", i)); 84 PetscCall(PetscObjectSetName((PetscObject)data_comp->ProlongationOps[i], "P")); 85 PetscCall(MatView(data_comp->ProlongationOps[i], PETSC_VIEWER_STDOUT_WORLD)); 86 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n")); 87 } 88 89 for (PetscInt i = 0; i < data_comp->num_levels - 1; i++) { 90 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Coarse/Fine splitting - Level %" PetscInt_FMT "\n", i + 1)); 91 PetscCall(PetscBTView(data_comp->n_per_level[i + 1], data_comp->CFMarkers[i], PETSC_VIEWER_STDOUT_WORLD)); 92 } 93 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Stiffness matrix, sparse format:\n")); 94 PetscCall(MatViewFromOptions(data_comp->stiff, NULL, "-mat_view_stiff")); 95 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Finished calling the Viewer functions\n")); 96 PetscFunctionReturn(PETSC_SUCCESS); 97 } 98 99 PetscErrorCode DataCompDestroy(DataCompression data_comp) 100 { 101 PetscFunctionBeginUser; 102 if (data_comp == NULL) PetscFunctionReturn(PETSC_SUCCESS); 103 PetscCall(MatDestroy(&data_comp->stiff)); 104 PetscCall(KSPDestroy(&data_comp->kspHypre)); 105 for (PetscInt i = 0; i < data_comp->num_levels - 1; i++) { 106 PetscCall(MatDestroy(&data_comp->ProlongationOps[i])); 107 PetscCall(PetscBTDestroy(&data_comp->CFMarkers[i])); 108 } 109 PetscCall(PetscFree(data_comp->ProlongationOps)); 110 PetscCall(PetscFree(data_comp->n_per_level)); 111 PetscCall(PetscFree(data_comp->CFMarkers)); 112 PetscCall(PetscFree(data_comp)); 113 PetscFunctionReturn(PETSC_SUCCESS); 114 } 115 116 /*TEST 117 118 test: 119 requires: hypre !complex !defined(PETSC_HAVE_HYPRE_DEVICE) 120 args: -pc_hypre_boomeramg_coarsen_type modifiedRuge-Stueben -pc_hypre_boomeramg_interp_type classical -pc_hypre_boomeramg_strong_threshold 0.25 pc_hypre_boomeramg_numfunctions 1 -pc_hypre_boomeramg_max_row_sum 1.0 -mat_view_stiff 121 122 TEST*/ 123