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