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