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