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