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