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 /* Domain and mesh definition */ 14 PetscBool benchmark; 15 PetscInt cells[3]; 16 PetscInt processGrid[3]; 17 PetscInt nodeGrid[3]; 18 } AppCtx; 19 20 21 static PetscErrorCode trig_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 22 { 23 PetscInt d; 24 *u = 0.0; 25 for (d = 0; d < dim; ++d) *u += PetscSinReal(2.0*PETSC_PI*x[d]); 26 return 0; 27 } 28 29 static void f0_trig_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 30 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 31 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 32 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 33 { 34 PetscInt d; 35 for (d = 0; d < dim; ++d) f0[0] += -4.0*PetscSqr(PETSC_PI)*PetscSinReal(2.0*PETSC_PI*x[d]); 36 } 37 38 static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 39 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 40 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 41 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 42 { 43 PetscInt d; 44 for (d = 0; d < dim; ++d) f1[d] = u_x[d]; 45 } 46 47 static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 48 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 49 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 50 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 51 { 52 PetscInt d; 53 for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0; 54 } 55 56 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 57 { 58 PetscErrorCode ierr; 59 PetscInt asz, dim=2; /* should be default of DMPLex (yuck) */ 60 61 PetscFunctionBeginUser; 62 options->benchmark= PETSC_FALSE; 63 for (asz=0;asz<3;asz++) options->processGrid[asz] = options->cells[asz] = options->nodeGrid[asz] = 1; 64 ierr = PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");CHKERRQ(ierr); 65 ierr = PetscOptionsInt("-dm_plex_box_dim","dim in ex13","ex13.c",dim,&dim,NULL);CHKERRQ(ierr); 66 ierr = PetscOptionsBool("-benchmark", "Solve the benchmark problem", "ex13.c", options->benchmark, &options->benchmark, NULL);CHKERRQ(ierr); 67 asz = dim; 68 ierr = PetscOptionsIntArray("-dm_plex_box_faces","Mesh size (cells) for benchmarking ex13","ex13.c",options->cells,&asz,NULL);CHKERRQ(ierr); 69 asz = dim; 70 ierr = PetscOptionsIntArray("-process_grid_size","Number of processors (np) in each dimension (cells[i]%np[i]==0 && prod(np[i]==#procs)","ex13.c",options->processGrid,&asz,NULL);CHKERRQ(ierr); 71 asz = dim; 72 ierr = PetscOptionsIntArray("-node_grid_size","Number of nodes (nnodes) in each dimension (np[i]%nnodes[i]==0)","ex13.c",options->nodeGrid,&asz,NULL);CHKERRQ(ierr); 73 ierr = PetscOptionsEnd(); 74 PetscFunctionReturn(0); 75 } 76 77 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 78 { 79 PetscErrorCode ierr; 80 PetscInt dim; 81 82 PetscFunctionBeginUser; 83 /* Create box mesh */ 84 ierr = DMPlexCreateBoxMesh(comm, 2, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);CHKERRQ(ierr); 85 /* TODO: This should be pulled into the library */ 86 { 87 char convType[256]; 88 PetscBool flg; 89 90 ierr = PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");CHKERRQ(ierr); 91 ierr = PetscOptionsFList("-dm_plex_convert_type","Convert DMPlex to another format","ex12",DMList,DMPLEX,convType,256,&flg);CHKERRQ(ierr); 92 ierr = PetscOptionsEnd(); 93 if (flg) { 94 DM dmConv; 95 96 ierr = DMConvert(*dm,convType,&dmConv);CHKERRQ(ierr); 97 if (dmConv) { 98 ierr = DMDestroy(dm);CHKERRQ(ierr); 99 *dm = dmConv; 100 } 101 } 102 } 103 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr); 104 105 ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr); 106 ierr = DMSetApplicationContext(*dm, user);CHKERRQ(ierr); 107 if (user->benchmark) { 108 PetscPartitioner part; 109 PetscInt cEnd, ii, np,cells_proc[3],procs_node[3]; 110 PetscMPIInt rank, size; 111 PetscInt *sizes = NULL, *points = NULL; 112 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 113 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 114 ierr = DMPlexGetHeightStratum(*dm, 0, NULL, &cEnd);CHKERRQ(ierr); 115 ierr = DMGetDimension(*dm, &dim);CHKERRQ(ierr); 116 for (ii=0,np=1;ii<dim;ii++) np *= user->processGrid[ii]; /* check number of processors */ 117 if (np!=size)SETERRQ2(comm,PETSC_ERR_SUP,"invalid process grid sum = %D, -n %D",np, size); 118 for (ii=0,np=1;ii<dim;ii++) np *= user->cells[ii]; 119 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 120 if (!rank) { 121 if (np!=cEnd)SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP," cell grid %D != num cells = %D",np,cEnd); 122 for (ii=0;ii<dim;ii++) { 123 if (user->processGrid[ii]%user->nodeGrid[ii]) SETERRQ3(comm,PETSC_ERR_SUP,"dir %D, invalid node grid size %D, process grid %D",ii,user->nodeGrid[ii],user->processGrid[ii]); 124 procs_node[ii] = user->processGrid[ii]/user->nodeGrid[ii]; 125 if (user->cells[ii]%user->processGrid[ii]) SETERRQ3(comm,PETSC_ERR_SUP,"dir %D, invalid process grid size %D, cells %D",ii,user->processGrid[ii],user->cells[ii]); 126 cells_proc[ii] = user->cells[ii]/user->processGrid[ii]; 127 ierr = PetscPrintf(comm, "%D) cells_proc=%D procs_node=%D processGrid=%D cells=%D nodeGrid=%D\n",ii,cells_proc[ii],procs_node[ii],user->processGrid[ii],user->cells[ii],user->nodeGrid[ii]);CHKERRQ(ierr); 128 } 129 for (/* */;ii<3;ii++) { 130 procs_node[ii] = cells_proc[ii] = 1; 131 ierr = PetscPrintf(comm, "%D) cells_proc=%D procs_node=%D processGrid=%D cells=%D nodeGrid=%D\n",ii,cells_proc[ii],procs_node[ii],user->processGrid[ii],user->cells[ii],user->nodeGrid[ii]);CHKERRQ(ierr); 132 } 133 PetscInt pi,pj,pk,ni,nj,nk,ci,cj,ck,pid=0; 134 ierr = PetscMalloc2(size, &sizes, cEnd, &points);CHKERRQ(ierr); 135 for (ii=0,np=1;ii<dim;ii++) np *= cells_proc[ii]; 136 for (ii=0;ii<size;ii++) sizes[ii] = np; 137 for (ii=0;ii<cEnd;ii++) points[ii] = -1; 138 for (nk=0;nk<user->nodeGrid[2];nk++) { /* node loop */ 139 PetscInt idx_2 = nk*cells_proc[2]*procs_node[2]*user->cells[0]*user->cells[1]; 140 for (nj=0;nj<user->nodeGrid[1];nj++) { /* node loop */ 141 PetscInt idx_1 = idx_2 + nj*cells_proc[1]*procs_node[1]*user->cells[0]; 142 for (ni=0;ni<user->nodeGrid[0];ni++) { /* node loop */ 143 PetscInt idx_0 = idx_1 + ni*cells_proc[0]*procs_node[0]; 144 for (pk=0;pk<procs_node[2];pk++) { /* process loop */ 145 PetscInt idx_22 = idx_0 + pk*cells_proc[2]*user->cells[0]*user->cells[1]; 146 for (pj=0;pj<procs_node[1];pj++) { /* process loop */ 147 PetscInt idx_11 = idx_22 + pj*cells_proc[1]*user->cells[0]; 148 for (pi=0;pi<procs_node[0];pi++) { /* process loop */ 149 PetscInt idx_00 = idx_11 + pi*cells_proc[0]; 150 for (ck=0;ck<cells_proc[2];ck++) { /* cell loop */ 151 PetscInt idx_222 = idx_00 + ck*user->cells[0]*user->cells[1]; 152 for (cj=0;cj<cells_proc[1];cj++) { /* cell loop */ 153 PetscInt idx_111 = idx_222 + cj*user->cells[0]; 154 for (ci=0;ci<cells_proc[0];ci++) { /* cell loop */ 155 PetscInt idx_000 = idx_111 + ci; 156 points[idx_000] = pid; 157 } 158 } 159 } 160 pid++; 161 } 162 } 163 } 164 } 165 } 166 } 167 if (pid!=size) SETERRQ2(comm,PETSC_ERR_SUP,"pid %D != size %D",pid,size); 168 /* view */ 169 ierr = PetscPrintf(comm, "points:\n");CHKERRQ(ierr); 170 pid=0; 171 for (ck=0;ck<user->cells[2];ck++) { 172 for (cj=0;cj<user->cells[1];cj++) { 173 for (ci=0;ci<user->cells[0];ci++) { 174 ierr = PetscPrintf(comm, "%6D",points[pid++]);CHKERRQ(ierr); 175 } 176 ierr = PetscPrintf(comm, "\n");CHKERRQ(ierr); 177 } 178 ierr = PetscPrintf(comm, "\n");CHKERRQ(ierr); 179 } 180 } 181 ierr = DMPlexGetPartitioner(*dm, &part);CHKERRQ(ierr); 182 ierr = PetscPartitionerSetType(part, PETSCPARTITIONERSHELL);CHKERRQ(ierr); 183 ierr = PetscPartitionerShellSetPartition(part, size, sizes, points);CHKERRQ(ierr); 184 ierr = PetscFree2(sizes, points);CHKERRQ(ierr); 185 } 186 ierr = DMSetFromOptions(*dm);CHKERRQ(ierr); 187 ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr); 188 PetscFunctionReturn(0); 189 } 190 191 static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user) 192 { 193 PetscDS prob; 194 const PetscInt id = 1; 195 PetscErrorCode ierr; 196 197 PetscFunctionBeginUser; 198 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 199 ierr = PetscDSSetResidual(prob, 0, f0_trig_u, f1_u);CHKERRQ(ierr); 200 ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr); 201 ierr = PetscDSSetExactSolution(prob, 0, trig_u, user);CHKERRQ(ierr); 202 ierr = DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) trig_u, NULL, 1, &id, user);CHKERRQ(ierr); 203 PetscFunctionReturn(0); 204 } 205 206 static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user) 207 { 208 DM cdm = dm; 209 PetscFE fe; 210 DMPolytopeType ct; 211 PetscBool simplex; 212 PetscInt dim, cStart; 213 char prefix[PETSC_MAX_PATH_LEN]; 214 PetscErrorCode ierr; 215 216 PetscFunctionBeginUser; 217 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 218 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, NULL);CHKERRQ(ierr); 219 ierr = DMPlexGetCellType(dm, cStart, &ct);CHKERRQ(ierr); 220 simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE; 221 /* Create finite element */ 222 ierr = PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);CHKERRQ(ierr); 223 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, simplex, name ? prefix : NULL, -1, &fe);CHKERRQ(ierr); 224 ierr = PetscObjectSetName((PetscObject) fe, name);CHKERRQ(ierr); 225 /* Set discretization and boundary conditions for each mesh */ 226 ierr = DMSetField(dm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr); 227 ierr = DMCreateDS(dm);CHKERRQ(ierr); 228 ierr = (*setup)(dm, user);CHKERRQ(ierr); 229 while (cdm) { 230 ierr = DMCopyDisc(dm,cdm);CHKERRQ(ierr); 231 /* TODO: Check whether the boundary of coarse meshes is marked */ 232 ierr = DMGetCoarseDM(cdm, &cdm);CHKERRQ(ierr); 233 } 234 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr); 235 PetscFunctionReturn(0); 236 } 237 238 int main(int argc, char **argv) 239 { 240 DM dm; /* Problem specification */ 241 SNES snes; /* Nonlinear solver */ 242 Vec u; /* Solutions */ 243 AppCtx user; /* User-defined work context */ 244 PetscErrorCode ierr; 245 246 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 247 ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr); 248 /* Primal system */ 249 ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr); 250 ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr); 251 ierr = SNESSetDM(snes, dm);CHKERRQ(ierr); 252 ierr = SetupDiscretization(dm, "potential", SetupPrimalProblem, &user);CHKERRQ(ierr); 253 ierr = DMCreateGlobalVector(dm, &u);CHKERRQ(ierr); 254 ierr = VecSet(u, 0.0);CHKERRQ(ierr); 255 ierr = PetscObjectSetName((PetscObject) u, "potential");CHKERRQ(ierr); 256 ierr = DMPlexSetSNESLocalFEM(dm, &user, &user, &user);CHKERRQ(ierr); 257 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 258 ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr); 259 /* Benchmark system */ 260 if (user.benchmark) { 261 #if defined(PETSC_USE_LOG) 262 PetscLogStage stage; 263 #endif 264 KSP ksp; 265 Vec b; 266 ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr); 267 ierr = SNESGetSolution(snes, &u);CHKERRQ(ierr); 268 ierr = VecZeroEntries(u);CHKERRQ(ierr); 269 ierr = SNESGetFunction(snes, &b, NULL, NULL);CHKERRQ(ierr); 270 ierr = SNESComputeFunction(snes, u, b);CHKERRQ(ierr); 271 ierr = PetscLogStageRegister("KSP Solve only", &stage);CHKERRQ(ierr); 272 ierr = PetscLogStagePush(stage);CHKERRQ(ierr); 273 ierr = KSPSolve(ksp, b, u);CHKERRQ(ierr); 274 ierr = PetscLogStagePop();CHKERRQ(ierr); 275 } 276 ierr = SNESGetSolution(snes, &u);CHKERRQ(ierr); 277 ierr = VecViewFromOptions(u, NULL, "-potential_view");CHKERRQ(ierr); 278 /* Cleanup */ 279 ierr = VecDestroy(&u);CHKERRQ(ierr); 280 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 281 ierr = DMDestroy(&dm);CHKERRQ(ierr); 282 ierr = PetscFinalize(); 283 return ierr; 284 } 285 286 /*TEST 287 288 test: 289 suffix: bench 290 nsize: 16 291 args: -dm_plex_box_dim 2 -dm_plex_box_faces 8,8 -ksp_type cg -pc_type gamg -dm_plex_box_simplex 0 -dm_refine 1 \ 292 -potential_petscspace_degree 2 -dm_distribute -petscpartitioner_type simple -dm_view \ 293 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1,1 -process_grid_size 4,4 -node_grid_size 2,2 -benchmark true 294 295 TEST*/ 296