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