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