1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 // libCEED + PETSc Example: Elasticity 18 // 19 // This example demonstrates a simple usage of libCEED with PETSc to solve 20 // elasticity problems. 21 // 22 // The code uses higher level communication protocols in DMPlex. 23 // 24 // Build with: 25 // 26 // make elasticity [PETSC_DIR=</path/to/petsc>] [CEED_DIR=</path/to/libceed>] 27 // 28 // Sample runs: 29 // 30 // ./elasticity -mesh [.exo file] -degree 2 -E 1 -nu 0.3 -problem linElas -forcing mms 31 // ./elasticity -mesh [.exo file] -degree 2 -E 1 -nu 0.3 -bc_clamp 998,999 -bc_clamp_998_translate 0.1,0.2,0.3 -problem hyperSS -forcing none -ceed /cpu/self 32 // ./elasticity -mesh [.exo file] -degree 2 -E 1 -nu 0.3 -bc_clamp 998,999 -bc_clamp_998 rotate 1,0,0,0.2 -problem hyperFS -forcing none -ceed /gpu/occa 33 // 34 // Sample meshes can be found at https://github.com/jeremylt/ceedSampleMeshes 35 // 36 //TESTARGS -ceed {ceed_resource} -test -degree 2 -nu 0.3 -E 1 -dm_plex_box_faces 3,3,3 37 38 /// @file 39 /// CEED elasticity example using PETSc with DMPlex 40 41 const char help[] = "Solve solid Problems with CEED and PETSc DMPlex\n"; 42 43 #include "elasticity.h" 44 45 int main(int argc, char **argv) { 46 PetscInt ierr; 47 MPI_Comm comm; 48 // Context structs 49 AppCtx appCtx; // Contains problem options 50 Physics phys; // Contains physical constants 51 Units units; // Contains units scaling 52 // PETSc objects 53 PetscLogStage stageDMSetup, stageLibceedSetup, 54 stageSnesSetup, stageSnesSolve; 55 DM dmOrig; // Distributed DM to clone 56 DM *levelDMs; 57 Vec U, *Ug, *Uloc; // U: solution, R: residual, F: forcing 58 Vec R, Rloc, F, Floc; // g: global, loc: local 59 SNES snes, snesCoarse = NULL; 60 Mat *jacobMat, jacobMatCoarse, *prolongRestrMat; 61 // PETSc data 62 UserMult resCtx, jacobCoarseCtx = NULL, *jacobCtx; 63 FormJacobCtx formJacobCtx; 64 UserMultProlongRestr *prolongRestrCtx; 65 PCMGCycleType pcmgCycleType = PC_MG_CYCLE_V; 66 // libCEED objects 67 Ceed ceed, ceedFine = NULL; 68 CeedData *ceedData; 69 CeedQFunction qfRestrict = NULL, qfProlong = NULL; 70 // Parameters 71 PetscInt ncompu = 3; // 3 DoFs in 3D 72 PetscInt numLevels = 1, fineLevel = 0; 73 PetscInt *Ugsz, *Ulsz, *Ulocsz; // sz: size 74 PetscInt snesIts = 0; 75 // Timing 76 double startTime, elapsedTime, minTime, maxTime; 77 78 ierr = PetscInitialize(&argc, &argv, NULL, help); 79 if (ierr) 80 return ierr; 81 82 // --------------------------------------------------------------------------- 83 // Process command line options 84 // --------------------------------------------------------------------------- 85 comm = PETSC_COMM_WORLD; 86 87 // -- Set mesh file, polynomial degree, problem type 88 ierr = PetscCalloc1(1, &appCtx); CHKERRQ(ierr); 89 ierr = ProcessCommandLineOptions(comm, appCtx); CHKERRQ(ierr); 90 numLevels = appCtx->numLevels; 91 fineLevel = numLevels - 1; 92 93 // -- Set Poison's ratio, Young's Modulus 94 ierr = PetscMalloc1(1, &phys); CHKERRQ(ierr); 95 ierr = PetscMalloc1(1, &units); CHKERRQ(ierr); 96 ierr = ProcessPhysics(comm, phys, units); CHKERRQ(ierr); 97 98 // --------------------------------------------------------------------------- 99 // Setup DM 100 // --------------------------------------------------------------------------- 101 // Performance logging 102 ierr = PetscLogStageRegister("DM and Vector Setup Stage", &stageDMSetup); 103 CHKERRQ(ierr); 104 ierr = PetscLogStagePush(stageDMSetup); CHKERRQ(ierr); 105 106 // -- Create distributed DM from mesh file 107 ierr = CreateDistributedDM(comm, appCtx, &dmOrig); CHKERRQ(ierr); 108 109 // -- Setup DM by polynomial degree 110 ierr = PetscMalloc1(numLevels, &levelDMs); CHKERRQ(ierr); 111 for (PetscInt level = 0; level < numLevels; level++) { 112 ierr = DMClone(dmOrig, &levelDMs[level]); CHKERRQ(ierr); 113 ierr = SetupDMByDegree(levelDMs[level], appCtx, appCtx->levelDegrees[level], 114 ncompu); CHKERRQ(ierr); 115 } 116 117 // --------------------------------------------------------------------------- 118 // Setup solution and work vectors 119 // --------------------------------------------------------------------------- 120 // Allocate arrays 121 ierr = PetscMalloc1(numLevels, &Ug); CHKERRQ(ierr); 122 ierr = PetscMalloc1(numLevels, &Uloc); CHKERRQ(ierr); 123 ierr = PetscMalloc1(numLevels, &Ugsz); CHKERRQ(ierr); 124 ierr = PetscMalloc1(numLevels, &Ulsz); CHKERRQ(ierr); 125 ierr = PetscMalloc1(numLevels, &Ulocsz); CHKERRQ(ierr); 126 127 // -- Setup solution vectors for each level 128 for (PetscInt level = 0; level < numLevels; level++) { 129 // -- Create global unknown vector U 130 ierr = DMCreateGlobalVector(levelDMs[level], &Ug[level]); CHKERRQ(ierr); 131 ierr = VecGetSize(Ug[level], &Ugsz[level]); CHKERRQ(ierr); 132 // Note: Local size for matShell 133 ierr = VecGetLocalSize(Ug[level], &Ulsz[level]); CHKERRQ(ierr); 134 135 // -- Create local unknown vector Uloc 136 ierr = DMCreateLocalVector(levelDMs[level], &Uloc[level]); CHKERRQ(ierr); 137 // Note: local size for libCEED 138 ierr = VecGetSize(Uloc[level], &Ulocsz[level]); CHKERRQ(ierr); 139 } 140 141 // -- Create residual and forcing vectors 142 ierr = VecDuplicate(Ug[fineLevel], &U); CHKERRQ(ierr); 143 ierr = VecDuplicate(Ug[fineLevel], &R); CHKERRQ(ierr); 144 ierr = VecDuplicate(Ug[fineLevel], &F); CHKERRQ(ierr); 145 ierr = VecDuplicate(Uloc[fineLevel], &Rloc); CHKERRQ(ierr); 146 ierr = VecDuplicate(Uloc[fineLevel], &Floc); CHKERRQ(ierr); 147 148 // Performance logging 149 ierr = PetscLogStagePop(); 150 151 // --------------------------------------------------------------------------- 152 // Set up libCEED 153 // --------------------------------------------------------------------------- 154 // Performance logging 155 ierr = PetscLogStageRegister("libCEED Setup Stage", &stageLibceedSetup); 156 CHKERRQ(ierr); 157 ierr = PetscLogStagePush(stageLibceedSetup); CHKERRQ(ierr); 158 159 // Initalize backend 160 CeedInit(appCtx->ceedResource, &ceed); 161 if (appCtx->degree > 4 && appCtx->ceedResourceFine[0]) 162 CeedInit(appCtx->ceedResourceFine, &ceedFine); 163 164 // -- Create libCEED local forcing vector 165 CeedVector forceCeed; 166 CeedScalar *f; 167 if (appCtx->forcingChoice != FORCE_NONE) { 168 ierr = VecGetArray(Floc, &f); CHKERRQ(ierr); 169 CeedVectorCreate(ceed, Ulocsz[fineLevel], &forceCeed); 170 CeedVectorSetArray(forceCeed, CEED_MEM_HOST, CEED_USE_POINTER, f); 171 } 172 173 // -- Restriction and prolongation QFunction 174 if (appCtx->multigridChoice != MULTIGRID_NONE) { 175 CeedQFunctionCreateIdentity(ceed, ncompu, CEED_EVAL_NONE, CEED_EVAL_INTERP, 176 &qfRestrict); 177 CeedQFunctionCreateIdentity(ceed, ncompu, CEED_EVAL_INTERP, CEED_EVAL_NONE, 178 &qfProlong); 179 } 180 181 // -- Setup libCEED objects 182 ierr = PetscMalloc1(numLevels, &ceedData); CHKERRQ(ierr); 183 // ---- Setup residual evaluator and geometric information 184 ierr = PetscCalloc1(1, &ceedData[fineLevel]); CHKERRQ(ierr); 185 { 186 bool highOrder = (appCtx->levelDegrees[fineLevel] > 4 && ceedFine); 187 Ceed levelCeed = highOrder ? ceedFine : ceed; 188 ierr = SetupLibceedFineLevel(levelDMs[fineLevel], levelCeed, appCtx, 189 phys, ceedData, fineLevel, ncompu, 190 Ugsz[fineLevel], Ulocsz[fineLevel], forceCeed, 191 qfRestrict, qfProlong); 192 CHKERRQ(ierr); 193 } 194 // ---- Setup Jacobian evaluator and prolongation/restriction 195 for (PetscInt level = 0; level < numLevels; level++) { 196 if (level != fineLevel) { 197 ierr = PetscCalloc1(1, &ceedData[level]); CHKERRQ(ierr); 198 } 199 200 // Note: use high order ceed, if specified and degree > 4 201 bool highOrder = (appCtx->levelDegrees[level] > 4 && ceedFine); 202 Ceed levelCeed = highOrder ? ceedFine : ceed; 203 ierr = SetupLibceedLevel(levelDMs[level], levelCeed, appCtx, phys, 204 ceedData, level, ncompu, Ugsz[level], 205 Ulocsz[level], forceCeed, qfRestrict, 206 qfProlong); CHKERRQ(ierr); 207 } 208 209 // Performance logging 210 ierr = PetscLogStagePop(); 211 212 // --------------------------------------------------------------------------- 213 // Setup global forcing vector 214 // --------------------------------------------------------------------------- 215 ierr = VecZeroEntries(F); CHKERRQ(ierr); 216 217 if (appCtx->forcingChoice != FORCE_NONE) { 218 ierr = VecRestoreArray(Floc, &f); CHKERRQ(ierr); 219 ierr = DMLocalToGlobal(levelDMs[fineLevel], Floc, ADD_VALUES, F); 220 CHKERRQ(ierr); 221 CeedVectorDestroy(&forceCeed); 222 } 223 224 // --------------------------------------------------------------------------- 225 // Print problem summary 226 // --------------------------------------------------------------------------- 227 if (!appCtx->testMode) { 228 const char *usedresource; 229 CeedGetResource(ceed, &usedresource); 230 231 ierr = PetscPrintf(comm, 232 "\n-- Elastisticy Example - libCEED + PETSc --\n" 233 " libCEED:\n" 234 " libCEED Backend : %s\n", 235 usedresource); CHKERRQ(ierr); 236 237 if (ceedFine) { 238 ierr = PetscPrintf(comm, 239 " libCEED Backend - high order : %s\n", 240 appCtx->ceedResourceFine); CHKERRQ(ierr); 241 } 242 243 ierr = PetscPrintf(comm, 244 " Problem:\n" 245 " Problem Name : %s\n" 246 " Forcing Function : %s\n" 247 " Mesh:\n" 248 " File : %s\n" 249 " Number of 1D Basis Nodes (p) : %d\n" 250 " Number of 1D Quadrature Points (q) : %d\n" 251 " Global nodes : %D\n" 252 " Owned nodes : %D\n" 253 " DoF per node : %D\n" 254 " Multigrid:\n" 255 " Type : %s\n" 256 " Number of Levels : %d\n", 257 problemTypesForDisp[appCtx->problemChoice], 258 forcingTypesForDisp[appCtx->forcingChoice], 259 appCtx->meshFile[0] ? appCtx->meshFile : "Box Mesh", 260 appCtx->degree + 1, appCtx->degree + 1, 261 Ugsz[fineLevel]/ncompu, Ulsz[fineLevel]/ncompu, ncompu, 262 (appCtx->degree == 1 && 263 appCtx->multigridChoice != MULTIGRID_NONE) ? 264 "Algebraic multigrid" : 265 multigridTypesForDisp[appCtx->multigridChoice], 266 (appCtx->degree == 1 || 267 appCtx->multigridChoice == MULTIGRID_NONE) ? 268 0 : numLevels); CHKERRQ(ierr); 269 270 if (appCtx->multigridChoice != MULTIGRID_NONE) { 271 for (PetscInt i = 0; i < 2; i++) { 272 CeedInt level = i ? fineLevel : 0; 273 ierr = PetscPrintf(comm, 274 " Level %D (%s):\n" 275 " Number of 1D Basis Nodes (p) : %d\n" 276 " Global Nodes : %D\n" 277 " Owned Nodes : %D\n", 278 level, i ? "fine" : "coarse", 279 appCtx->levelDegrees[level] + 1, 280 Ugsz[level]/ncompu, Ulsz[level]/ncompu); 281 CHKERRQ(ierr); 282 } 283 } 284 } 285 286 // --------------------------------------------------------------------------- 287 // Setup SNES 288 // --------------------------------------------------------------------------- 289 // Performance logging 290 ierr = PetscLogStageRegister("SNES Setup Stage", &stageSnesSetup); 291 CHKERRQ(ierr); 292 ierr = PetscLogStagePush(stageSnesSetup); CHKERRQ(ierr); 293 294 // Create SNES 295 ierr = SNESCreate(comm, &snes); CHKERRQ(ierr); 296 ierr = SNESSetDM(snes, levelDMs[fineLevel]); CHKERRQ(ierr); 297 298 // -- Jacobian evaluators 299 ierr = PetscMalloc1(numLevels, &jacobCtx); CHKERRQ(ierr); 300 ierr = PetscMalloc1(numLevels, &jacobMat); CHKERRQ(ierr); 301 for (PetscInt level = 0; level < numLevels; level++) { 302 // -- Jacobian context for level 303 ierr = PetscMalloc1(1, &jacobCtx[level]); CHKERRQ(ierr); 304 ierr = SetupJacobianCtx(comm, appCtx, levelDMs[level], Ug[level], 305 Uloc[level], ceedData[level], ceed, 306 jacobCtx[level]); CHKERRQ(ierr); 307 308 // -- Form Action of Jacobian on delta_u 309 ierr = MatCreateShell(comm, Ulsz[level], Ulsz[level], Ugsz[level], 310 Ugsz[level], jacobCtx[level], &jacobMat[level]); 311 CHKERRQ(ierr); 312 ierr = MatShellSetOperation(jacobMat[level], MATOP_MULT, 313 (void (*)(void))ApplyJacobian_Ceed); 314 CHKERRQ(ierr); 315 ierr = MatShellSetOperation(jacobMat[level], MATOP_GET_DIAGONAL, 316 (void(*)(void))GetDiag_Ceed); 317 318 } 319 // Note: FormJacobian updates Jacobian matrices on each level 320 // and assembles the Jpre matrix, if needed 321 ierr = PetscMalloc1(1, &formJacobCtx); CHKERRQ(ierr); 322 formJacobCtx->jacobCtx = jacobCtx; 323 formJacobCtx->numLevels = numLevels; 324 formJacobCtx->jacobMat = jacobMat; 325 326 // -- Residual evaluation function 327 ierr = PetscMalloc1(1, &resCtx); CHKERRQ(ierr); 328 ierr = PetscMemcpy(resCtx, jacobCtx[fineLevel], 329 sizeof(*jacobCtx[fineLevel])); CHKERRQ(ierr); 330 resCtx->op = ceedData[fineLevel]->opApply; 331 ierr = SNESSetFunction(snes, R, FormResidual_Ceed, resCtx); CHKERRQ(ierr); 332 333 // -- Prolongation/Restriction evaluation 334 ierr = PetscMalloc1(numLevels, &prolongRestrCtx); CHKERRQ(ierr); 335 ierr = PetscMalloc1(numLevels, &prolongRestrMat); CHKERRQ(ierr); 336 for (PetscInt level = 1; level < numLevels; level++) { 337 // ---- Prolongation/restriction context for level 338 ierr = PetscMalloc1(1, &prolongRestrCtx[level]); CHKERRQ(ierr); 339 ierr = SetupProlongRestrictCtx(comm, levelDMs[level-1], levelDMs[level], 340 Ug[level], Uloc[level-1], Uloc[level], 341 ceedData[level-1], ceedData[level], ceed, 342 prolongRestrCtx[level]); CHKERRQ(ierr); 343 344 // ---- Form Action of Jacobian on delta_u 345 ierr = MatCreateShell(comm, Ulsz[level], Ulsz[level-1], Ugsz[level], 346 Ugsz[level-1], prolongRestrCtx[level], 347 &prolongRestrMat[level]); CHKERRQ(ierr); 348 // Note: In PCMG, restriction is the transpose of prolongation 349 ierr = MatShellSetOperation(prolongRestrMat[level], MATOP_MULT, 350 (void (*)(void))Prolong_Ceed); 351 ierr = MatShellSetOperation(prolongRestrMat[level], MATOP_MULT_TRANSPOSE, 352 (void (*)(void))Restrict_Ceed); 353 CHKERRQ(ierr); 354 } 355 356 // --------------------------------------------------------------------------- 357 // Setup dummy SNES for AMG coarse solve 358 // --------------------------------------------------------------------------- 359 if (appCtx->multigridChoice != MULTIGRID_NONE) { 360 // -- Jacobian Matrix 361 ierr = DMSetMatType(levelDMs[0], MATAIJ); CHKERRQ(ierr); 362 ierr = DMCreateMatrix(levelDMs[0], &jacobMatCoarse); CHKERRQ(ierr); 363 364 if (appCtx->degree > 1) { 365 ierr = SNESCreate(comm, &snesCoarse); CHKERRQ(ierr); 366 ierr = SNESSetDM(snesCoarse, levelDMs[0]); CHKERRQ(ierr); 367 ierr = SNESSetSolution(snesCoarse, Ug[0]); CHKERRQ(ierr); 368 369 // -- Jacobian function 370 ierr = SNESSetJacobian(snesCoarse, jacobMatCoarse, jacobMatCoarse, NULL, 371 NULL); CHKERRQ(ierr); 372 373 // -- Residual evaluation function 374 ierr = PetscMalloc1(1, &jacobCoarseCtx); CHKERRQ(ierr); 375 ierr = PetscMemcpy(jacobCoarseCtx, jacobCtx[0], sizeof(*jacobCtx[0])); 376 CHKERRQ(ierr); 377 ierr = SNESSetFunction(snesCoarse, Ug[0], ApplyJacobianCoarse_Ceed, 378 jacobCoarseCtx); CHKERRQ(ierr); 379 380 // -- Update formJacobCtx 381 formJacobCtx->Ucoarse = Ug[0]; 382 formJacobCtx->snesCoarse = snesCoarse; 383 formJacobCtx->jacobMatCoarse = jacobMatCoarse; 384 } 385 } 386 387 // Set Jacobian function 388 if (appCtx->degree > 1) { 389 ierr = SNESSetJacobian(snes, jacobMat[fineLevel], jacobMat[fineLevel], 390 FormJacobian, formJacobCtx); CHKERRQ(ierr); 391 } else { 392 ierr = SNESSetJacobian(snes, jacobMat[0], jacobMatCoarse, 393 SNESComputeJacobianDefaultColor, NULL); 394 CHKERRQ(ierr); 395 } 396 397 // --------------------------------------------------------------------------- 398 // Setup KSP 399 // --------------------------------------------------------------------------- 400 { 401 PC pc; 402 KSP ksp; 403 404 // -- KSP 405 ierr = SNESGetKSP(snes, &ksp); CHKERRQ(ierr); 406 ierr = KSPSetType(ksp, KSPCG); CHKERRQ(ierr); 407 ierr = KSPSetNormType(ksp, KSP_NORM_NATURAL); CHKERRQ(ierr); 408 ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 409 PETSC_DEFAULT); CHKERRQ(ierr); 410 ierr = KSPSetOptionsPrefix(ksp, "outer_"); CHKERRQ(ierr); 411 412 // -- Preconditioning 413 ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr); 414 ierr = PCSetDM(pc, levelDMs[fineLevel]); CHKERRQ(ierr); 415 ierr = PCSetOptionsPrefix(pc, "outer_"); CHKERRQ(ierr); 416 417 if (appCtx->multigridChoice == MULTIGRID_NONE) { 418 // ---- No Multigrid 419 ierr = PCSetType(pc, PCJACOBI); CHKERRQ(ierr); 420 ierr = PCJacobiSetType(pc, PC_JACOBI_DIAGONAL); CHKERRQ(ierr); 421 } else if (appCtx->degree == 1) { 422 // ---- AMG for degree 1 423 ierr = PCSetType(pc, PCGAMG); CHKERRQ(ierr); 424 } else { 425 // ---- PCMG 426 ierr = PCSetType(pc, PCMG); CHKERRQ(ierr); 427 428 // ------ PCMG levels 429 ierr = PCMGSetLevels(pc, numLevels, NULL); CHKERRQ(ierr); 430 for (PetscInt level = 0; level < numLevels; level++) { 431 // -------- Smoother 432 KSP kspSmoother, kspEst; 433 PC pcSmoother; 434 435 // ---------- Smoother KSP 436 ierr = PCMGGetSmoother(pc, level, &kspSmoother); CHKERRQ(ierr); 437 ierr = KSPSetDM(kspSmoother, levelDMs[level]); CHKERRQ(ierr); 438 ierr = KSPSetDMActive(kspSmoother, PETSC_FALSE); CHKERRQ(ierr); 439 440 // ---------- Chebyshev options 441 ierr = KSPSetType(kspSmoother, KSPCHEBYSHEV); CHKERRQ(ierr); 442 ierr = KSPChebyshevEstEigSet(kspSmoother, 0, 0.1, 0, 1.1); 443 CHKERRQ(ierr); 444 ierr = KSPChebyshevEstEigGetKSP(kspSmoother, &kspEst); CHKERRQ(ierr); 445 ierr = KSPSetType(kspEst, KSPCG); CHKERRQ(ierr); 446 ierr = KSPChebyshevEstEigSetUseNoisy(kspSmoother, PETSC_TRUE); 447 CHKERRQ(ierr); 448 ierr = KSPSetOperators(kspSmoother, jacobMat[level], jacobMat[level]); 449 CHKERRQ(ierr); 450 451 // ---------- Smoother preconditioner 452 ierr = KSPGetPC(kspSmoother, &pcSmoother); CHKERRQ(ierr); 453 ierr = PCSetType(pcSmoother, PCJACOBI); CHKERRQ(ierr); 454 ierr = PCJacobiSetType(pcSmoother, PC_JACOBI_DIAGONAL); CHKERRQ(ierr); 455 456 // -------- Work vector 457 if (level != fineLevel) { 458 ierr = PCMGSetX(pc, level, Ug[level]); CHKERRQ(ierr); 459 } 460 461 // -------- Level prolongation/restriction operator 462 if (level > 0) { 463 ierr = PCMGSetInterpolation(pc, level, prolongRestrMat[level]); 464 CHKERRQ(ierr); 465 ierr = PCMGSetRestriction(pc, level, prolongRestrMat[level]); 466 CHKERRQ(ierr); 467 } 468 } 469 470 // ------ PCMG coarse solve 471 KSP kspCoarse; 472 PC pcCoarse; 473 474 // -------- Coarse KSP 475 ierr = PCMGGetCoarseSolve(pc, &kspCoarse); CHKERRQ(ierr); 476 ierr = KSPSetType(kspCoarse, KSPPREONLY); CHKERRQ(ierr); 477 ierr = KSPSetOperators(kspCoarse, jacobMatCoarse, jacobMatCoarse); 478 CHKERRQ(ierr); 479 ierr = KSPSetOptionsPrefix(kspCoarse, "coarse_"); CHKERRQ(ierr); 480 481 // -------- Coarse preconditioner 482 ierr = KSPGetPC(kspCoarse, &pcCoarse); CHKERRQ(ierr); 483 ierr = PCSetType(pcCoarse, PCGAMG); CHKERRQ(ierr); 484 ierr = PCSetOptionsPrefix(pcCoarse, "coarse_"); CHKERRQ(ierr); 485 486 ierr = KSPSetFromOptions(kspCoarse); CHKERRQ(ierr); 487 ierr = PCSetFromOptions(pcCoarse); CHKERRQ(ierr); 488 489 // ------ PCMG options 490 ierr = PCMGSetType(pc, PC_MG_MULTIPLICATIVE); CHKERRQ(ierr); 491 ierr = PCMGSetNumberSmooth(pc, 3); CHKERRQ(ierr); 492 ierr = PCMGSetCycleType(pc, pcmgCycleType); CHKERRQ(ierr); 493 } 494 ierr = KSPSetFromOptions(ksp); 495 ierr = PCSetFromOptions(pc); 496 } 497 { 498 // Default to critical-point (CP) line search (related to Wolfe's curvature condition) 499 SNESLineSearch linesearch; 500 501 ierr = SNESGetLineSearch(snes, &linesearch); CHKERRQ(ierr); 502 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHCP); CHKERRQ(ierr); 503 } 504 505 ierr = SNESSetFromOptions(snes); CHKERRQ(ierr); 506 507 // Performance logging 508 ierr = PetscLogStagePop(); 509 510 // --------------------------------------------------------------------------- 511 // Set initial guess 512 // --------------------------------------------------------------------------- 513 ierr = VecSet(U, 0.0); CHKERRQ(ierr); 514 515 // View solution 516 if (appCtx->viewSoln) { 517 ierr = ViewSolution(comm, U, 0, 0.0); CHKERRQ(ierr); 518 } 519 520 // --------------------------------------------------------------------------- 521 // Solve SNES 522 // --------------------------------------------------------------------------- 523 PetscBool snesMonitor = PETSC_FALSE; 524 ierr = PetscOptionsHasName(NULL, NULL, "-snes_monitor", &snesMonitor); 525 CHKERRQ(ierr); 526 527 // Performance logging 528 ierr = PetscLogStageRegister("SNES Solve Stage", &stageSnesSolve); 529 CHKERRQ(ierr); 530 ierr = PetscLogStagePush(stageSnesSolve); CHKERRQ(ierr); 531 532 // Timing 533 ierr = PetscBarrier((PetscObject)snes); CHKERRQ(ierr); 534 startTime = MPI_Wtime(); 535 536 // Solve for each load increment 537 PetscInt increment; 538 for (increment = 1; increment <= appCtx->numIncrements; increment++) { 539 // -- Log increment count 540 if (snesMonitor) { 541 ierr = PetscPrintf(comm, "%d Load Increment\n", increment - 1); 542 CHKERRQ(ierr); 543 } 544 545 // -- Scale the problem 546 PetscScalar loadIncrement = 1.0*increment / appCtx->numIncrements, 547 scalingFactor = loadIncrement / 548 (increment == 1 ? 1 : resCtx->loadIncrement); 549 resCtx->loadIncrement = loadIncrement; 550 if (appCtx->numIncrements > 1 && appCtx->forcingChoice != FORCE_NONE) { 551 ierr = VecScale(F, scalingFactor); CHKERRQ(ierr); 552 } 553 554 // -- Solve 555 ierr = SNESSolve(snes, F, U); CHKERRQ(ierr); 556 557 // -- View solution 558 if (appCtx->viewSoln) { 559 ierr = ViewSolution(comm, U, increment, loadIncrement); CHKERRQ(ierr); 560 } 561 562 // -- Update SNES iteration count 563 PetscInt its; 564 ierr = SNESGetIterationNumber(snes, &its); CHKERRQ(ierr); 565 snesIts += its; 566 567 // -- Check for divergence 568 SNESConvergedReason reason; 569 ierr = SNESGetConvergedReason(snes, &reason); CHKERRQ(ierr); 570 if (reason < 0) 571 break; 572 } 573 574 // Timing 575 elapsedTime = MPI_Wtime() - startTime; 576 577 // Performance logging 578 ierr = PetscLogStagePop(); 579 580 // View solution 581 if (appCtx->viewFinalSoln && !appCtx->viewSoln) { 582 ierr = ViewSolution(comm, U, increment - 1, 1.0); CHKERRQ(ierr); 583 } 584 585 // --------------------------------------------------------------------------- 586 // Output summary 587 // --------------------------------------------------------------------------- 588 if (!appCtx->testMode) { 589 // -- SNES 590 SNESType snesType; 591 SNESConvergedReason reason; 592 PetscReal rnorm; 593 ierr = SNESGetType(snes, &snesType); CHKERRQ(ierr); 594 ierr = SNESGetConvergedReason(snes, &reason); CHKERRQ(ierr); 595 ierr = SNESGetFunctionNorm(snes, &rnorm); CHKERRQ(ierr); 596 ierr = PetscPrintf(comm, 597 " SNES:\n" 598 " SNES Type : %s\n" 599 " SNES Convergence : %s\n" 600 " Number of Load Increments : %d\n" 601 " Completed Load Increments : %d\n" 602 " Total SNES Iterations : %D\n" 603 " Final rnorm : %e\n", 604 snesType, SNESConvergedReasons[reason], 605 appCtx->numIncrements, increment - 1, 606 snesIts, (double)rnorm); CHKERRQ(ierr); 607 608 // -- KSP 609 KSP ksp; 610 KSPType kspType; 611 ierr = SNESGetKSP(snes, &ksp); CHKERRQ(ierr); 612 ierr = KSPGetType(ksp, &kspType); CHKERRQ(ierr); 613 ierr = PetscPrintf(comm, 614 " Linear Solver:\n" 615 " KSP Type : %s\n", 616 kspType); CHKERRQ(ierr); 617 618 // -- PC 619 PC pc; 620 PCType pcType; 621 ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr); 622 ierr = PCGetType(pc, &pcType); CHKERRQ(ierr); 623 ierr = PetscPrintf(comm, 624 " PC Type : %s\n", 625 pcType); CHKERRQ(ierr); 626 627 if (!strcmp(pcType, PCMG)) { 628 PCMGType pcmgType; 629 ierr = PCMGGetType(pc, &pcmgType); CHKERRQ(ierr); 630 ierr = PetscPrintf(comm, 631 " P-Multigrid:\n" 632 " PCMG Type : %s\n" 633 " PCMG Cycle Type : %s\n", 634 PCMGTypes[pcmgType], 635 PCMGCycleTypes[pcmgCycleType]); CHKERRQ(ierr); 636 637 // -- Coarse Solve 638 KSP kspCoarse; 639 PC pcCoarse; 640 PCType pcType; 641 642 ierr = PCMGGetCoarseSolve(pc, &kspCoarse); CHKERRQ(ierr); 643 ierr = KSPGetType(kspCoarse, &kspType); CHKERRQ(ierr); 644 ierr = KSPGetPC(kspCoarse, &pcCoarse); CHKERRQ(ierr); 645 ierr = PCGetType(pcCoarse, &pcType); CHKERRQ(ierr); 646 ierr = PetscPrintf(comm, 647 " Coarse Solve:\n" 648 " KSP Type : %s\n" 649 " PC Type : %s\n", 650 kspType, pcType); CHKERRQ(ierr); 651 } 652 } 653 654 // --------------------------------------------------------------------------- 655 // Compute solve time 656 // --------------------------------------------------------------------------- 657 if (!appCtx->testMode) { 658 ierr = MPI_Allreduce(&elapsedTime, &minTime, 1, MPI_DOUBLE, MPI_MIN, comm); 659 CHKERRQ(ierr); 660 ierr = MPI_Allreduce(&elapsedTime, &maxTime, 1, MPI_DOUBLE, MPI_MAX, comm); 661 CHKERRQ(ierr); 662 ierr = PetscPrintf(comm, 663 " Performance:\n" 664 " SNES Solve Time : %g (%g) sec\n", 665 maxTime, minTime); CHKERRQ(ierr); 666 } 667 668 // --------------------------------------------------------------------------- 669 // Compute error 670 // --------------------------------------------------------------------------- 671 if (appCtx->forcingChoice == FORCE_MMS) { 672 CeedScalar l2Error = 1., l2Unorm = 1.; 673 const CeedScalar *truearray; 674 Vec errorVec, trueVec; 675 676 // -- Work vectors 677 ierr = VecDuplicate(U, &errorVec); CHKERRQ(ierr); 678 ierr = VecSet(errorVec, 0.0); CHKERRQ(ierr); 679 ierr = VecDuplicate(U, &trueVec); CHKERRQ(ierr); 680 ierr = VecSet(trueVec, 0.0); CHKERRQ(ierr); 681 682 // -- Assemble global true solution vector 683 CeedVectorGetArrayRead(ceedData[fineLevel]->truesoln, CEED_MEM_HOST, 684 &truearray); 685 ierr = VecPlaceArray(resCtx->Yloc, truearray); CHKERRQ(ierr); 686 ierr = DMLocalToGlobal(resCtx->dm, resCtx->Yloc, INSERT_VALUES, trueVec); 687 CHKERRQ(ierr); 688 ierr = VecResetArray(resCtx->Yloc); CHKERRQ(ierr); 689 CeedVectorRestoreArrayRead(ceedData[fineLevel]->truesoln, &truearray); 690 691 // -- Compute L2 error 692 ierr = VecWAXPY(errorVec, -1.0, U, trueVec); CHKERRQ(ierr); 693 ierr = VecNorm(errorVec, NORM_2, &l2Error); CHKERRQ(ierr); 694 ierr = VecNorm(U, NORM_2, &l2Unorm); CHKERRQ(ierr); 695 l2Error /= l2Unorm; 696 697 // -- Output 698 if (!appCtx->testMode || l2Error > 0.05) { 699 ierr = PetscPrintf(comm, 700 " L2 Error : %e\n", 701 l2Error); CHKERRQ(ierr); 702 } 703 704 // -- Cleanup 705 ierr = VecDestroy(&errorVec); CHKERRQ(ierr); 706 ierr = VecDestroy(&trueVec); CHKERRQ(ierr); 707 } 708 709 // --------------------------------------------------------------------------- 710 // Compute energy 711 // --------------------------------------------------------------------------- 712 if (!appCtx->testMode) { 713 // -- Compute L2 error 714 CeedScalar energy; 715 ierr = ComputeStrainEnergy(resCtx, 716 ceedData[fineLevel]->opEnergy, U, 717 ceedData[fineLevel]->energy, &energy); 718 CHKERRQ(ierr); 719 720 // -- Output 721 ierr = PetscPrintf(comm, 722 " Strain Energy : %e\n", 723 energy); CHKERRQ(ierr); 724 } 725 726 // --------------------------------------------------------------------------- 727 // Free objects 728 // --------------------------------------------------------------------------- 729 // Data in arrays per level 730 for (PetscInt level = 0; level < numLevels; level++) { 731 // Vectors 732 ierr = VecDestroy(&Ug[level]); CHKERRQ(ierr); 733 ierr = VecDestroy(&Uloc[level]); CHKERRQ(ierr); 734 735 // Jacobian matrix and data 736 ierr = VecDestroy(&jacobCtx[level]->Yloc); CHKERRQ(ierr); 737 ierr = MatDestroy(&jacobMat[level]); CHKERRQ(ierr); 738 ierr = PetscFree(jacobCtx[level]); CHKERRQ(ierr); 739 740 // Prolongation/Restriction matrix and data 741 if (level > 0) { 742 ierr = VecDestroy(&prolongRestrCtx[level]->multVec); CHKERRQ(ierr); 743 ierr = PetscFree(prolongRestrCtx[level]); CHKERRQ(ierr); 744 ierr = MatDestroy(&prolongRestrMat[level]); CHKERRQ(ierr); 745 } 746 747 // DM 748 ierr = DMDestroy(&levelDMs[level]); CHKERRQ(ierr); 749 750 // libCEED objects 751 ierr = CeedDataDestroy(level, ceedData[level]); CHKERRQ(ierr); 752 } 753 754 // Arrays 755 ierr = PetscFree(Ug); CHKERRQ(ierr); 756 ierr = PetscFree(Uloc); CHKERRQ(ierr); 757 ierr = PetscFree(Ugsz); CHKERRQ(ierr); 758 ierr = PetscFree(Ulsz); CHKERRQ(ierr); 759 ierr = PetscFree(Ulocsz); CHKERRQ(ierr); 760 ierr = PetscFree(jacobCtx); CHKERRQ(ierr); 761 ierr = PetscFree(jacobMat); CHKERRQ(ierr); 762 ierr = PetscFree(prolongRestrCtx); CHKERRQ(ierr); 763 ierr = PetscFree(prolongRestrMat); CHKERRQ(ierr); 764 ierr = PetscFree(appCtx->levelDegrees); CHKERRQ(ierr); 765 ierr = PetscFree(ceedData); CHKERRQ(ierr); 766 767 // libCEED objects 768 CeedQFunctionDestroy(&qfRestrict); 769 CeedQFunctionDestroy(&qfProlong); 770 CeedDestroy(&ceed); 771 CeedDestroy(&ceedFine); 772 773 // PETSc objects 774 ierr = VecDestroy(&U); CHKERRQ(ierr); 775 ierr = VecDestroy(&R); CHKERRQ(ierr); 776 ierr = VecDestroy(&Rloc); CHKERRQ(ierr); 777 ierr = VecDestroy(&F); CHKERRQ(ierr); 778 ierr = VecDestroy(&Floc); CHKERRQ(ierr); 779 ierr = MatDestroy(&jacobMatCoarse); CHKERRQ(ierr); 780 ierr = SNESDestroy(&snes); CHKERRQ(ierr); 781 ierr = SNESDestroy(&snesCoarse); CHKERRQ(ierr); 782 ierr = DMDestroy(&dmOrig); CHKERRQ(ierr); 783 ierr = PetscFree(levelDMs); CHKERRQ(ierr); 784 785 // Structs 786 ierr = PetscFree(resCtx); CHKERRQ(ierr); 787 ierr = PetscFree(formJacobCtx); CHKERRQ(ierr); 788 ierr = PetscFree(jacobCoarseCtx); CHKERRQ(ierr); 789 ierr = PetscFree(appCtx); CHKERRQ(ierr); 790 ierr = PetscFree(phys); CHKERRQ(ierr); 791 ierr = PetscFree(units); CHKERRQ(ierr); 792 793 return PetscFinalize(); 794 } 795