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