1 static char help[] = "Benchmark Poisson Problem in 2d and 3d with finite elements.\n\ 2 We solve the Poisson problem in a rectangular\n\ 3 domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\ 4 This example supports automatic convergence estimation\n\ 5 and eventually adaptivity.\n\n\n"; 6 7 #include <petscdmplex.h> 8 #include <petscsnes.h> 9 #include <petscds.h> 10 #include <petscconvest.h> 11 12 typedef struct { 13 PetscInt nit; 14 } AppCtx; 15 16 17 static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 18 { 19 PetscInt d; 20 *u = 0.0; 21 for (d = 0; d < dim; ++d) *u += PetscSinReal(2.0*PETSC_PI*x[d]); 22 return 0; 23 } 24 25 static void f0_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 26 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 27 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 28 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 29 { 30 PetscInt d; 31 for (d = 0; d < dim; ++d) f0[0] += -4.0*PetscSqr(PETSC_PI)*PetscSinReal(2.0*PETSC_PI*x[d]); 32 } 33 34 static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 35 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 36 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 37 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 38 { 39 PetscInt d; 40 for (d = 0; d < dim; ++d) f1[d] = u_x[d]; 41 } 42 43 static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 44 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 45 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 46 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 47 { 48 PetscInt d; 49 for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0; 50 } 51 52 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 53 { 54 PetscErrorCode ierr; 55 56 PetscFunctionBeginUser; 57 options->nit = 10; 58 ierr = PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");CHKERRQ(ierr); 59 ierr = PetscOptionsInt("-benchmark_it", "Solve the benchmark problem this many times", "ex13.c", options->nit, &options->nit, NULL);CHKERRQ(ierr); 60 ierr = PetscOptionsEnd(); 61 PetscFunctionReturn(0); 62 } 63 64 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 65 { 66 PetscErrorCode ierr; 67 68 PetscFunctionBeginUser; 69 /* Create box mesh */ 70 ierr = DMPlexCreateBoxMesh(comm, 2, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr); 71 /* TODO: This should be pulled into the library */ 72 { 73 char convType[256]; 74 PetscBool flg; 75 76 ierr = PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");CHKERRQ(ierr); 77 ierr = PetscOptionsFList("-dm_plex_convert_type","Convert DMPlex to another format","ex12",DMList,DMPLEX,convType,256,&flg);CHKERRQ(ierr); 78 ierr = PetscOptionsEnd(); 79 if (flg) { 80 DM dmConv; 81 82 ierr = DMConvert(*dm,convType,&dmConv);CHKERRQ(ierr); 83 if (dmConv) { 84 ierr = DMDestroy(dm);CHKERRQ(ierr); 85 *dm = dmConv; 86 } 87 } 88 } 89 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); 90 91 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr); 92 ierr = DMSetApplicationContext(*dm, user);CHKERRQ(ierr); 93 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 94 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 95 PetscFunctionReturn(0); 96 } 97 98 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 99 { 100 PetscDS prob; 101 const PetscInt id = 1; 102 PetscErrorCode ierr; 103 104 PetscFunctionBeginUser; 105 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 106 ierr = PetscDSSetResidual(prob, 0, f0_trig_u, f1_u);CHKERRQ(ierr); 107 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr); 108 ierr = PetscDSSetExactSolution(prob, 0, trig_u, user);CHKERRQ(ierr); 109 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) trig_u, NULL, 1, &id, user);CHKERRQ(ierr); 110 PetscFunctionReturn(0); 111 } 112 113 static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user) 114 { 115 DM cdm = dm; 116 PetscFE fe; 117 DMPolytopeType ct; 118 PetscBool simplex; 119 PetscInt dim, cStart; 120 char prefix[PETSC_MAX_PATH_LEN]; 121 PetscErrorCode ierr; 122 123 PetscFunctionBeginUser; 124 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 125 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr); 126 ierr = DMPlexGetCellType(dm, cStart, &ct);CHKERRQ(ierr); 127 simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE; 128 /* Create finite element */ 129 ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);CHKERRQ(ierr); 130 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, simplex, name ? prefix : NULL, -1, &fe);CHKERRQ(ierr); 131 ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr); 132 /* Set discretization and boundary conditions for each mesh */ 133 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 134 ierr = DMCreateDS(dm);CHKERRQ(ierr); 135 ierr = (*setup)(dm, user);CHKERRQ(ierr); 136 while (cdm) { 137 ierr = DMCopyDisc(dm,cdm);CHKERRQ(ierr); 138 /* TODO: Check whether the boundary of coarse meshes is marked */ 139 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 140 } 141 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 142 PetscFunctionReturn(0); 143 } 144 145 int main(int argc, char **argv) 146 { 147 DM dm; /* Problem specification */ 148 SNES snes; /* Nonlinear solver */ 149 Vec u; /* Solutions */ 150 AppCtx user; /* User-defined work context */ 151 PetscErrorCode ierr; 152 153 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 154 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 155 /* Primal system */ 156 ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr); 157 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 158 ierr = SNESSetDM(snes, dm);CHKERRQ(ierr); 159 ierr = SetupDiscretization(dm, "potential", SetupPrimalProblem, &user);CHKERRQ(ierr); 160 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr); 161 ierr = VecSet(u, 0.0);CHKERRQ(ierr); 162 ierr = PetscObjectSetName((PetscObject) u, "potential");CHKERRQ(ierr); 163 ierr = DMPlexSetSNESLocalFEM(dm, &user, &user, &user);CHKERRQ(ierr); 164 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 165 ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr); 166 /* Benchmark system */ 167 if (user.nit) { 168 #if defined(PETSC_USE_LOG) 169 PetscLogStage kspstage,pcstage; 170 #endif 171 KSP ksp; 172 PC pc; 173 Mat A,P; 174 Vec b; 175 PetscInt i; 176 ierr = PetscOptionsClearValue(NULL,"-ksp_monitor");CHKERRQ(ierr); 177 ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr); 178 ierr = SNESGetSolution(snes, &u);CHKERRQ(ierr); 179 ierr = SNESGetJacobian(snes, &A, &P, NULL, NULL);CHKERRQ(ierr); 180 ierr = VecSet(u, 0.0);CHKERRQ(ierr); 181 ierr = SNESGetFunction(snes, &b, NULL, NULL);CHKERRQ(ierr); 182 ierr = SNESComputeFunction(snes, u, b);CHKERRQ(ierr); 183 ierr = SNESComputeJacobian(snes, u, A, P);CHKERRQ(ierr); 184 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 185 ierr = PetscLogStageRegister("PCSetUp", &pcstage);CHKERRQ(ierr); 186 ierr = PetscLogStagePush(pcstage);CHKERRQ(ierr); 187 ierr = PCSetUp(pc);CHKERRQ(ierr); 188 ierr = PetscLogStagePop();CHKERRQ(ierr); 189 ierr = PetscLogStageRegister("KSP Solve only", &kspstage);CHKERRQ(ierr); 190 ierr = PetscLogStagePush(kspstage);CHKERRQ(ierr); 191 for (i=0;i<user.nit;i++) { 192 ierr = VecZeroEntries(u);CHKERRQ(ierr); 193 ierr = KSPSolve(ksp, b, u);CHKERRQ(ierr); 194 } 195 ierr = PetscLogStagePop();CHKERRQ(ierr); 196 } 197 ierr = SNESGetSolution(snes, &u);CHKERRQ(ierr); 198 ierr = VecViewFromOptions(u, NULL, "-potential_view");CHKERRQ(ierr); 199 /* Cleanup */ 200 ierr = VecDestroy(&u);CHKERRQ(ierr); 201 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 202 ierr = DMDestroy(&dm);CHKERRQ(ierr); 203 ierr = PetscFinalize(); 204 return ierr; 205 } 206 207 /*TEST 208 209 test: 210 suffix: bench 211 nsize: 4 212 args: -dm_plex_box_dim 3 -dm_plex_box_simplex 0 -dm_plex_box_faces 2,2,8 -dm_refine 1 -dm_distribute \ 213 -petscpartitioner_type simple -petscpartitioner_simple_process_grid 1,1,2 -petscpartitioner_simple_node_grid 1,1,2 \ 214 -potential_petscspace_degree 2 -ksp_type cg -pc_type gamg -benchmark_it 1 -dm_view -snes_rtol 1.e-4 215 216 test: 217 suffix: comparison 218 nsize: 4 219 args: -dm_plex_box_dim 2 -dm_plex_box_faces 4,4 -dm_refine 3 -petscpartitioner_simple_process_grid 2,2 \ 220 -petscpartitioner_simple_node_grid 1,1 -potential_petscspace_degree 2 -dm_distribute -petscpartitioner_type simple \ 221 -dm_plex_box_simplex 0 -snes_monitor_short -snes_type ksponly -dm_view -pc_type gamg -pc_gamg_process_eq_limit 400 -ksp_norm_type unpreconditioned \ 222 -pc_gamg_coarse_eq_limit 10 -snes_converged_reason -ksp_converged_reason -snes_rtol 1.e-4 223 224 test: 225 suffix: cuda 226 nsize: 4 227 requires: cuda 228 output_file: output/ex13_comparison.out 229 args: -dm_plex_box_dim 2 -dm_plex_box_faces 4,4 -dm_refine 3 -petscpartitioner_simple_process_grid 2,2 \ 230 -petscpartitioner_simple_node_grid 1,1 -potential_petscspace_degree 2 -dm_distribute -petscpartitioner_type simple \ 231 -dm_plex_box_simplex 0 -snes_monitor_short -snes_type ksponly -dm_view -pc_type gamg -pc_gamg_process_eq_limit 400 -ksp_norm_type unpreconditioned \ 232 -pc_gamg_coarse_eq_limit 10 -snes_converged_reason -ksp_converged_reason -snes_rtol 1.e-4 -dm_mat_type aijcusparse -dm_vec_type cuda 233 234 test: 235 suffix: kokkos_comp 236 nsize: 4 237 requires: kokkos_kernels 238 output_file: output/ex13_comparison.out 239 args: -dm_plex_box_dim 2 -dm_plex_box_faces 4,4 -dm_refine 3 -petscpartitioner_simple_process_grid 2,2 \ 240 -petscpartitioner_simple_node_grid 1,1 -potential_petscspace_degree 2 -dm_distribute -petscpartitioner_type simple \ 241 -dm_plex_box_simplex 0 -snes_monitor_short -snes_type ksponly -dm_view -pc_type gamg -pc_gamg_process_eq_limit 400 -ksp_norm_type unpreconditioned \ 242 -pc_gamg_coarse_eq_limit 10 -snes_converged_reason -ksp_converged_reason -snes_rtol 1.e-4 -dm_mat_type aijkokkos -dm_vec_type kokkos 243 244 test: 245 nsize: 4 246 requires: kokkos_kernels 247 suffix: kokkos 248 args: -dm_plex_box_dim 2 -dm_plex_box_faces 2,8 -dm_distribute -petscpartitioner_type simple -petscpartitioner_simple_process_grid 2,1 \ 249 -petscpartitioner_simple_node_grid 2,1 -dm_plex_box_simplex 0 -potential_petscspace_degree 1 -dm_refine 1 -ksp_type cg -pc_type gamg -ksp_norm_type unpreconditioned \ 250 -mg_levels_esteig_ksp_type cg -mg_levels_pc_type jacobi -ksp_converged_reason -snes_monitor_short -snes_rtol 1.e-4 -dm_view -dm_mat_type aijkokkos -dm_vec_type kokkos 251 252 test: 253 suffix: aijmkl_comp 254 nsize: 4 255 requires: mkl 256 output_file: output/ex13_comparison.out 257 args: -dm_plex_box_dim 2 -dm_plex_box_faces 4,4 -dm_refine 3 -petscpartitioner_simple_process_grid 2,2 \ 258 -petscpartitioner_simple_node_grid 1,1 -potential_petscspace_degree 2 -dm_distribute -petscpartitioner_type simple \ 259 -dm_plex_box_simplex 0 -snes_monitor_short -snes_type ksponly -dm_view -pc_type gamg -pc_gamg_process_eq_limit 400 -ksp_norm_type unpreconditioned \ 260 -pc_gamg_coarse_eq_limit 10 -snes_converged_reason -ksp_converged_reason -snes_rtol 1.e-4 -dm_mat_type aijmkl 261 262 test: 263 suffix: aijmkl_seq 264 nsize: 1 265 requires: mkl 266 TODO: broken (INDEFINITE PC) 267 args: -dm_plex_box_dim 3 -dm_plex_box_faces 4,4,4 -dm_refine 1 -petscpartitioner_type simple -potential_petscspace_degree 1 -dm_distribute -dm_plex_box_simplex 0 -snes_monitor_short -snes_type ksponly -dm_view -pc_type gamg -pc_gamg_sym_graph 0 -pc_gamg_threshold -1 -pc_gamg_square_graph 10 -pc_gamg_process_eq_limit 400 -pc_gamg_reuse_interpolation -pc_gamg_coarse_eq_limit 10 -mg_levels_esteig_ksp_type cg -mg_levels_pc_type jacobi -ksp_type cg -ksp_norm_type unpreconditioned -snes_converged_reason -ksp_converged_reason -snes_rtol 1.e-4 -dm_mat_type aijmkl -dm_vec_type standard 268 269 TEST*/ 270