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