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