1 // Copyright (c) 2017-2024, 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") -ceed {ceed_resource} -test -problem FS-NH -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") -ceed {ceed_resource} -test -problem FS-MR -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 char hostname[PETSC_MAX_PATH_LEN]; 311 const char *used_resource; 312 PetscMPIInt comm_size; 313 314 CeedGetResource(ceed, &used_resource); 315 PetscCall(PetscGetHostName(hostname, sizeof hostname)); 316 PetscCall(MPI_Comm_size(comm, &comm_size)); 317 318 PetscCall(PetscPrintf(comm, 319 "\n-- Elasticity Example - libCEED + PETSc --\n" 320 " MPI:\n" 321 " Hostname : %s\n" 322 " Total ranks : %d\n" 323 " libCEED:\n" 324 " libCEED Backend : %s\n" 325 " libCEED Backend MemType : %s\n", 326 hostname, comm_size, used_resource, CeedMemTypes[mem_type_backend])); 327 328 VecType vecType; 329 PetscCall(VecGetType(U, &vecType)); 330 PetscCall(PetscPrintf(comm, 331 " PETSc:\n" 332 " PETSc Vec Type : %s\n", 333 vecType)); 334 335 PetscCall(PetscPrintf(comm, 336 " Problem:\n" 337 " Problem Name : %s\n" 338 " Forcing Function : %s\n" 339 " Mesh:\n" 340 " File : %s\n" 341 " Number of 1D Basis Nodes (p) : %" PetscInt_FMT "\n" 342 " Number of 1D Quadrature Points (q) : %" PetscInt_FMT "\n" 343 " Global nodes : %" PetscInt_FMT "\n" 344 " Owned nodes : %" PetscInt_FMT "\n" 345 " DoF per node : %" PetscInt_FMT "\n" 346 " Multigrid:\n" 347 " Type : %s\n" 348 " Number of Levels : %" PetscInt_FMT "\n", 349 app_ctx->name_for_disp, forcing_types_for_disp[app_ctx->forcing_choice], 350 app_ctx->mesh_file[0] ? app_ctx->mesh_file : "Box Mesh", app_ctx->degree + 1, app_ctx->degree + 1, 351 U_g_size[fine_level] / num_comp_u, U_l_size[fine_level] / num_comp_u, num_comp_u, 352 (app_ctx->degree == 1 && app_ctx->multigrid_choice != MULTIGRID_NONE) ? "Algebraic multigrid" 353 : multigrid_types_for_disp[app_ctx->multigrid_choice], 354 (app_ctx->degree == 1 || app_ctx->multigrid_choice == MULTIGRID_NONE) ? 0 : num_levels)); 355 356 if (app_ctx->multigrid_choice != MULTIGRID_NONE) { 357 for (PetscInt i = 0; i < 2; i++) { 358 PetscInt level = i ? fine_level : 0; 359 360 PetscCall(PetscPrintf(comm, 361 " Level %" PetscInt_FMT " (%s):\n" 362 " Number of 1D Basis Nodes (p) : %" PetscInt_FMT "\n" 363 " Global Nodes : %" PetscInt_FMT "\n" 364 " Owned Nodes : %" PetscInt_FMT "\n", 365 level, i ? "fine" : "coarse", app_ctx->level_degrees[level] + 1, U_g_size[level] / num_comp_u, 366 U_l_size[level] / num_comp_u)); 367 } 368 } 369 } 370 371 // --------------------------------------------------------------------------- 372 // Setup SNES 373 // --------------------------------------------------------------------------- 374 // Performance logging 375 PetscCall(PetscLogStageRegister("SNES Setup Stage", &stage_snes_setup)); 376 PetscCall(PetscLogStagePush(stage_snes_setup)); 377 378 // Create SNES 379 PetscCall(SNESCreate(comm, &snes)); 380 PetscCall(SNESSetDM(snes, level_dms[fine_level])); 381 382 // -- Jacobian evaluators 383 PetscCall(PetscMalloc1(num_levels, &jacob_ctx)); 384 PetscCall(PetscMalloc1(num_levels, &jacob_mat)); 385 for (PetscInt level = 0; level < num_levels; level++) { 386 // -- Jacobian context for level 387 PetscCall(PetscMalloc1(1, &jacob_ctx[level])); 388 PetscCall(SetupJacobianCtx(comm, app_ctx, level_dms[level], U_g[level], U_loc[level], ceed_data[level], ceed, ctx_phys, ctx_phys_smoother, 389 jacob_ctx[level])); 390 391 // -- Form Action of Jacobian on delta_u 392 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])); 393 PetscCall(MatShellSetOperation(jacob_mat[level], MATOP_MULT, (void (*)(void))ApplyJacobian_Ceed)); 394 PetscCall(MatShellSetOperation(jacob_mat[level], MATOP_GET_DIAGONAL, (void (*)(void))GetDiag_Ceed)); 395 PetscCall(MatShellSetVecType(jacob_mat[level], vectype)); 396 } 397 // Note: FormJacobian updates Jacobian matrices on each level and assembles the Jpre matrix, if needed 398 PetscCall(PetscMalloc1(1, &form_jacob_ctx)); 399 form_jacob_ctx->jacob_ctx = jacob_ctx; 400 form_jacob_ctx->num_levels = num_levels; 401 form_jacob_ctx->jacob_mat = jacob_mat; 402 403 // -- Residual evaluation function 404 PetscCall(PetscCalloc1(1, &res_ctx)); 405 PetscCall(PetscMemcpy(res_ctx, jacob_ctx[fine_level], sizeof(*jacob_ctx[fine_level]))); 406 res_ctx->op = ceed_data[fine_level]->op_residual; 407 res_ctx->qf = ceed_data[fine_level]->qf_residual; 408 if (app_ctx->bc_traction_count > 0) res_ctx->neumann_bcs = neumann_bcs; 409 else res_ctx->neumann_bcs = NULL; 410 PetscCall(SNESSetFunction(snes, R, FormResidual_Ceed, res_ctx)); 411 412 // -- Prolongation/Restriction evaluation 413 PetscCall(PetscMalloc1(num_levels, &prolong_restr_ctx)); 414 PetscCall(PetscMalloc1(num_levels, &prolong_restr_mat)); 415 for (PetscInt level = 1; level < num_levels; level++) { 416 // ---- Prolongation/restriction context for level 417 PetscCall(PetscMalloc1(1, &prolong_restr_ctx[level])); 418 PetscCall(SetupProlongRestrictCtx(comm, app_ctx, level_dms[level - 1], level_dms[level], U_g[level], U_loc[level - 1], U_loc[level], 419 ceed_data[level - 1], ceed_data[level], ceed, prolong_restr_ctx[level])); 420 421 // ---- Form Action of Jacobian on delta_u 422 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], 423 &prolong_restr_mat[level])); 424 // Note: In PCMG, restriction is the transpose of prolongation 425 PetscCall(MatShellSetOperation(prolong_restr_mat[level], MATOP_MULT, (void (*)(void))Prolong_Ceed)); 426 PetscCall(MatShellSetOperation(prolong_restr_mat[level], MATOP_MULT_TRANSPOSE, (void (*)(void))Restrict_Ceed)); 427 PetscCall(MatShellSetVecType(prolong_restr_mat[level], vectype)); 428 } 429 430 // --------------------------------------------------------------------------- 431 // Setup for AMG coarse solve 432 // --------------------------------------------------------------------------- 433 if (app_ctx->multigrid_choice != MULTIGRID_NONE) { 434 // -- Jacobian Matrix 435 PetscCall(DMCreateMatrix(level_dms[0], &jacob_mat_coarse)); 436 437 if (app_ctx->degree > 1) { 438 // -- Assemble sparsity pattern 439 PetscCount num_entries; 440 PetscInt *rows_petsc, *cols_petsc; 441 CeedInt *rows_ceed, *cols_ceed; 442 CeedVector coo_values; 443 ISLocalToGlobalMapping ltog_row, ltog_col; 444 445 CeedOperatorLinearAssembleSymbolic(ceed_data[0]->op_jacobian, &num_entries, &rows_ceed, &cols_ceed); 446 PetscCall(IntArrayCeedToPetsc(num_entries, &rows_ceed, &rows_petsc)); 447 PetscCall(IntArrayCeedToPetsc(num_entries, &cols_ceed, &cols_petsc)); 448 PetscCall(MatGetLocalToGlobalMapping(jacob_mat_coarse, <og_row, <og_col)); 449 PetscCall(ISLocalToGlobalMappingApply(ltog_row, num_entries, rows_petsc, rows_petsc)); 450 PetscCall(ISLocalToGlobalMappingApply(ltog_col, num_entries, cols_petsc, cols_petsc)); 451 PetscCall(MatSetPreallocationCOO(jacob_mat_coarse, num_entries, rows_petsc, cols_petsc)); 452 free(rows_petsc); 453 free(cols_petsc); 454 CeedVectorCreate(ceed, num_entries, &coo_values); 455 456 // -- Update form_jacob_ctx 457 form_jacob_ctx->coo_values = coo_values; 458 form_jacob_ctx->op_coarse = ceed_data[0]->op_jacobian; 459 form_jacob_ctx->jacob_mat_coarse = jacob_mat_coarse; 460 } 461 } 462 463 // Set Jacobian function 464 if (app_ctx->degree > 1) { 465 PetscCall(SNESSetJacobian(snes, jacob_mat[fine_level], jacob_mat[fine_level], FormJacobian, form_jacob_ctx)); 466 } else { 467 PetscCall(SNESSetJacobian(snes, jacob_mat[0], jacob_mat_coarse, SNESComputeJacobianDefaultColor, NULL)); 468 } 469 470 // --------------------------------------------------------------------------- 471 // Setup KSP 472 // --------------------------------------------------------------------------- 473 { 474 PC pc; 475 KSP ksp; 476 477 // -- KSP 478 PetscCall(SNESGetKSP(snes, &ksp)); 479 PetscCall(KSPSetType(ksp, KSPCG)); 480 PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL)); 481 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); 482 PetscCall(KSPSetOptionsPrefix(ksp, "outer_")); 483 484 // -- Preconditioning 485 PetscCall(KSPGetPC(ksp, &pc)); 486 PetscCall(PCSetDM(pc, level_dms[fine_level])); 487 PetscCall(PCSetOptionsPrefix(pc, "outer_")); 488 489 if (app_ctx->multigrid_choice == MULTIGRID_NONE) { 490 // ---- No Multigrid 491 PetscCall(PCSetType(pc, PCJACOBI)); 492 PetscCall(PCJacobiSetType(pc, PC_JACOBI_DIAGONAL)); 493 } else if (app_ctx->degree == 1) { 494 // ---- AMG for degree 1 495 PetscCall(PCSetType(pc, PCGAMG)); 496 } else { 497 // ---- PCMG 498 PetscCall(PCSetType(pc, PCMG)); 499 500 // ------ PCMG levels 501 PetscCall(PCMGSetLevels(pc, num_levels, NULL)); 502 for (PetscInt level = 0; level < num_levels; level++) { 503 // -------- Smoother 504 KSP ksp_smoother, ksp_est; 505 PC pc_smoother; 506 507 // ---------- Smoother KSP 508 PetscCall(PCMGGetSmoother(pc, level, &ksp_smoother)); 509 PetscCall(KSPSetDM(ksp_smoother, level_dms[level])); 510 PetscCall(KSPSetDMActive(ksp_smoother, PETSC_FALSE)); 511 512 // ---------- Chebyshev options 513 PetscCall(KSPSetType(ksp_smoother, KSPCHEBYSHEV)); 514 PetscCall(KSPChebyshevEstEigSet(ksp_smoother, 0, 0.1, 0, 1.1)); 515 PetscCall(KSPChebyshevEstEigGetKSP(ksp_smoother, &ksp_est)); 516 PetscCall(KSPSetType(ksp_est, KSPCG)); 517 PetscCall(KSPChebyshevEstEigSetUseNoisy(ksp_smoother, PETSC_TRUE)); 518 PetscCall(KSPSetOperators(ksp_smoother, jacob_mat[level], jacob_mat[level])); 519 520 // ---------- Smoother preconditioner 521 PetscCall(KSPGetPC(ksp_smoother, &pc_smoother)); 522 PetscCall(PCSetType(pc_smoother, PCJACOBI)); 523 PetscCall(PCJacobiSetType(pc_smoother, PC_JACOBI_DIAGONAL)); 524 525 // -------- Work vector 526 if (level != fine_level) { 527 PetscCall(PCMGSetX(pc, level, U_g[level])); 528 } 529 530 // -------- Level prolongation/restriction operator 531 if (level > 0) { 532 PetscCall(PCMGSetInterpolation(pc, level, prolong_restr_mat[level])); 533 PetscCall(PCMGSetRestriction(pc, level, prolong_restr_mat[level])); 534 } 535 } 536 537 // ------ PCMG coarse solve 538 KSP ksp_coarse; 539 PC pc_coarse; 540 541 // -------- Coarse KSP 542 PetscCall(PCMGGetCoarseSolve(pc, &ksp_coarse)); 543 PetscCall(KSPSetType(ksp_coarse, KSPPREONLY)); 544 PetscCall(KSPSetOperators(ksp_coarse, jacob_mat_coarse, jacob_mat_coarse)); 545 PetscCall(KSPSetOptionsPrefix(ksp_coarse, "coarse_")); 546 547 // -------- Coarse preconditioner 548 PetscCall(KSPGetPC(ksp_coarse, &pc_coarse)); 549 PetscCall(PCSetType(pc_coarse, PCGAMG)); 550 PetscCall(PCSetOptionsPrefix(pc_coarse, "coarse_")); 551 552 PetscCall(KSPSetFromOptions(ksp_coarse)); 553 PetscCall(PCSetFromOptions(pc_coarse)); 554 555 // ------ PCMG options 556 PetscCall(PCMGSetType(pc, PC_MG_MULTIPLICATIVE)); 557 PetscCall(PCMGSetNumberSmooth(pc, 3)); 558 PetscCall(PCMGSetCycleType(pc, pcmg_cycle_type)); 559 } 560 PetscCall(KSPSetFromOptions(ksp)); 561 PetscCall(PCSetFromOptions(pc)); 562 } 563 { 564 // Default to critical-point (CP) line search (related to Wolfe's curvature condition) 565 SNESLineSearch line_search; 566 567 PetscCall(SNESGetLineSearch(snes, &line_search)); 568 PetscCall(SNESLineSearchSetType(line_search, SNESLINESEARCHCP)); 569 } 570 571 PetscCall(SNESSetFromOptions(snes)); 572 573 // Performance logging 574 PetscCall(PetscLogStagePop()); 575 576 // --------------------------------------------------------------------------- 577 // Set initial guess 578 // --------------------------------------------------------------------------- 579 PetscCall(PetscObjectSetName((PetscObject)U, "")); 580 PetscCall(VecSet(U, 0.0)); 581 582 // View solution 583 if (app_ctx->view_soln) { 584 PetscCall(ViewSolution(comm, app_ctx, U, 0, 0.0)); 585 } 586 587 // --------------------------------------------------------------------------- 588 // Solve SNES 589 // --------------------------------------------------------------------------- 590 PetscBool snes_monitor = PETSC_FALSE; 591 PetscCall(PetscOptionsHasName(NULL, NULL, "-snes_monitor", &snes_monitor)); 592 593 // Performance logging 594 PetscCall(PetscLogStageRegister("SNES Solve Stage", &stage_snes_solve)); 595 PetscCall(PetscLogStagePush(stage_snes_solve)); 596 597 // Timing 598 PetscCall(PetscBarrier((PetscObject)snes)); 599 start_time = MPI_Wtime(); 600 601 // Solve for each load increment 602 PetscInt increment; 603 for (increment = 1; increment <= app_ctx->num_increments; increment++) { 604 // -- Log increment count 605 if (snes_monitor) { 606 PetscCall(PetscPrintf(comm, "%" PetscInt_FMT " Load Increment\n", increment - 1)); 607 } 608 609 // -- Scale the problem 610 PetscScalar load_increment = 1.0 * increment / app_ctx->num_increments, 611 scalingFactor = load_increment / (increment == 1 ? 1 : res_ctx->load_increment); 612 res_ctx->load_increment = load_increment; 613 if (app_ctx->num_increments > 1 && app_ctx->forcing_choice != FORCE_NONE) { 614 PetscCall(VecScale(F, scalingFactor)); 615 } 616 617 // -- Solve 618 PetscCall(SNESSolve(snes, F, U)); 619 620 // -- View solution 621 if (app_ctx->view_soln) { 622 PetscCall(ViewSolution(comm, app_ctx, U, increment, load_increment)); 623 } 624 625 // -- Update SNES iteration count 626 PetscInt its; 627 PetscCall(SNESGetIterationNumber(snes, &its)); 628 snes_its += its; 629 PetscCall(SNESGetLinearSolveIterations(snes, &its)); 630 ksp_its += its; 631 632 // -- Check for divergence 633 SNESConvergedReason reason; 634 PetscCall(SNESGetConvergedReason(snes, &reason)); 635 if (reason < 0) break; 636 if (app_ctx->energy_viewer) { 637 // -- Log strain energy for current load increment 638 CeedScalar energy; 639 PetscCall(ComputeStrainEnergy(dm_energy, res_ctx, ceed_data[fine_level]->op_energy, U, &energy)); 640 641 if (!app_ctx->test_mode) { 642 // -- Output 643 PetscCall(PetscPrintf(comm, " Strain Energy : %.12e\n", energy)); 644 } 645 PetscCall(PetscViewerASCIIPrintf(app_ctx->energy_viewer, "%f,%e\n", load_increment, energy)); 646 } 647 } 648 649 // Timing 650 elapsed_time = MPI_Wtime() - start_time; 651 652 // Performance logging 653 PetscCall(PetscLogStagePop()); 654 655 // --------------------------------------------------------------------------- 656 // Output summary 657 // --------------------------------------------------------------------------- 658 if (!app_ctx->test_mode) { 659 // -- SNES 660 SNESType snes_type; 661 SNESConvergedReason reason; 662 PetscReal rnorm; 663 PetscCall(SNESGetType(snes, &snes_type)); 664 PetscCall(SNESGetConvergedReason(snes, &reason)); 665 PetscCall(SNESGetFunctionNorm(snes, &rnorm)); 666 PetscCall(PetscPrintf(comm, 667 " SNES:\n" 668 " SNES Type : %s\n" 669 " SNES Convergence : %s\n" 670 " Number of Load Increments : %" PetscInt_FMT "\n" 671 " Completed Load Increments : %" PetscInt_FMT "\n" 672 " Total SNES Iterations : %" PetscInt_FMT "\n" 673 " Final rnorm : %e\n", 674 snes_type, SNESConvergedReasons[reason], app_ctx->num_increments, increment - 1, snes_its, (double)rnorm)); 675 676 // -- KSP 677 KSP ksp; 678 KSPType ksp_type; 679 PetscCall(SNESGetKSP(snes, &ksp)); 680 PetscCall(KSPGetType(ksp, &ksp_type)); 681 PetscCall(PetscPrintf(comm, 682 " Linear Solver:\n" 683 " KSP Type : %s\n" 684 " Total KSP Iterations : %" PetscInt_FMT "\n", 685 ksp_type, ksp_its)); 686 687 // -- PC 688 PC pc; 689 PCType pc_type; 690 PetscCall(KSPGetPC(ksp, &pc)); 691 PetscCall(PCGetType(pc, &pc_type)); 692 PetscCall(PetscPrintf(comm, " PC Type : %s\n", pc_type)); 693 694 if (!strcmp(pc_type, PCMG)) { 695 PCMGType pcmg_type; 696 PetscCall(PCMGGetType(pc, &pcmg_type)); 697 PetscCall(PetscPrintf(comm, 698 " P-Multigrid:\n" 699 " PCMG Type : %s\n" 700 " PCMG Cycle Type : %s\n", 701 PCMGTypes[pcmg_type], PCMGCycleTypes[pcmg_cycle_type])); 702 703 // -- Coarse Solve 704 KSP ksp_coarse; 705 PC pc_coarse; 706 PCType pc_type; 707 708 PetscCall(PCMGGetCoarseSolve(pc, &ksp_coarse)); 709 PetscCall(KSPGetType(ksp_coarse, &ksp_type)); 710 PetscCall(KSPGetPC(ksp_coarse, &pc_coarse)); 711 PetscCall(PCGetType(pc_coarse, &pc_type)); 712 PetscCall(PetscPrintf(comm, 713 " Coarse Solve:\n" 714 " KSP Type : %s\n" 715 " PC Type : %s\n", 716 ksp_type, pc_type)); 717 } 718 } 719 720 // --------------------------------------------------------------------------- 721 // Compute solve time 722 // --------------------------------------------------------------------------- 723 if (!app_ctx->test_mode) { 724 PetscCall(MPI_Allreduce(&elapsed_time, &min_time, 1, MPI_DOUBLE, MPI_MIN, comm)); 725 PetscCall(MPI_Allreduce(&elapsed_time, &max_time, 1, MPI_DOUBLE, MPI_MAX, comm)); 726 PetscCall(PetscPrintf(comm, 727 " Performance:\n" 728 " SNES Solve Time : %g (%g) sec\n" 729 " DoFs/Sec in SNES : %g (%g) million\n", 730 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)); 731 } 732 733 // --------------------------------------------------------------------------- 734 // Compute error 735 // --------------------------------------------------------------------------- 736 if (app_ctx->forcing_choice == FORCE_MMS) { 737 CeedScalar l2_error = 1., l2_U_norm = 1.; 738 const CeedScalar *true_array; 739 Vec error_vec, true_vec; 740 741 // -- Work vectors 742 PetscCall(VecDuplicate(U, &error_vec)); 743 PetscCall(VecSet(error_vec, 0.0)); 744 PetscCall(VecDuplicate(U, &true_vec)); 745 PetscCall(VecSet(true_vec, 0.0)); 746 747 // -- Assemble global true solution vector 748 CeedVectorGetArrayRead(ceed_data[fine_level]->true_soln, CEED_MEM_HOST, &true_array); 749 PetscCall(VecPlaceArray(res_ctx->Y_loc, (PetscScalar *)true_array)); 750 PetscCall(DMLocalToGlobal(res_ctx->dm, res_ctx->Y_loc, INSERT_VALUES, true_vec)); 751 PetscCall(VecResetArray(res_ctx->Y_loc)); 752 CeedVectorRestoreArrayRead(ceed_data[fine_level]->true_soln, &true_array); 753 754 // -- Compute L2 error 755 PetscCall(VecWAXPY(error_vec, -1.0, U, true_vec)); 756 PetscCall(VecNorm(error_vec, NORM_2, &l2_error)); 757 PetscCall(VecNorm(U, NORM_2, &l2_U_norm)); 758 l2_error /= l2_U_norm; 759 760 // -- Output 761 if (!app_ctx->test_mode || l2_error > 0.05) { 762 PetscCall(PetscPrintf(comm, " L2 Error : %e\n", l2_error)); 763 } 764 765 // -- Cleanup 766 PetscCall(VecDestroy(&error_vec)); 767 PetscCall(VecDestroy(&true_vec)); 768 } 769 770 // --------------------------------------------------------------------------- 771 // Compute energy 772 // --------------------------------------------------------------------------- 773 PetscReal energy; 774 PetscCall(ComputeStrainEnergy(dm_energy, res_ctx, ceed_data[fine_level]->op_energy, U, &energy)); 775 if (!app_ctx->test_mode) { 776 // -- Output 777 PetscCall(PetscPrintf(comm, " Strain Energy : %.12e\n", energy)); 778 } 779 PetscCall(RegressionTests_solids(app_ctx, energy)); 780 781 // --------------------------------------------------------------------------- 782 // Output diagnostic quantities 783 // --------------------------------------------------------------------------- 784 if (app_ctx->view_soln || app_ctx->view_final_soln) { 785 // -- Setup context 786 UserMult diagnostic_ctx; 787 PetscCall(PetscMalloc1(1, &diagnostic_ctx)); 788 PetscCall(PetscMemcpy(diagnostic_ctx, res_ctx, sizeof(*res_ctx))); 789 diagnostic_ctx->dm = dm_diagnostic; 790 diagnostic_ctx->op = ceed_data[fine_level]->op_diagnostic; 791 792 // -- Compute and output 793 PetscCall(ViewDiagnosticQuantities(comm, level_dms[fine_level], diagnostic_ctx, app_ctx, U, ceed_data[fine_level]->elem_restr_diagnostic)); 794 795 // -- Cleanup 796 PetscCall(PetscFree(diagnostic_ctx)); 797 } 798 799 // --------------------------------------------------------------------------- 800 // Free objects 801 // --------------------------------------------------------------------------- 802 // Data in arrays per level 803 for (PetscInt level = 0; level < num_levels; level++) { 804 // Vectors 805 PetscCall(VecDestroy(&U_g[level])); 806 PetscCall(VecDestroy(&U_loc[level])); 807 808 // Jacobian matrix and data 809 PetscCall(VecDestroy(&jacob_ctx[level]->Y_loc)); 810 PetscCall(MatDestroy(&jacob_mat[level])); 811 PetscCall(PetscFree(jacob_ctx[level])); 812 813 // Prolongation/Restriction matrix and data 814 if (level > 0) { 815 PetscCall(PetscFree(prolong_restr_ctx[level])); 816 PetscCall(MatDestroy(&prolong_restr_mat[level])); 817 } 818 819 // DM 820 PetscCall(DMDestroy(&level_dms[level])); 821 822 // libCEED objects 823 PetscCall(CeedDataDestroy(level, ceed_data[level])); 824 } 825 826 PetscCall(PetscViewerDestroy(&app_ctx->energy_viewer)); 827 828 // Arrays 829 PetscCall(PetscFree(U_g)); 830 PetscCall(PetscFree(U_loc)); 831 PetscCall(PetscFree(U_g_size)); 832 PetscCall(PetscFree(U_l_size)); 833 PetscCall(PetscFree(U_loc_size)); 834 PetscCall(PetscFree(jacob_ctx)); 835 PetscCall(PetscFree(jacob_mat)); 836 PetscCall(PetscFree(prolong_restr_ctx)); 837 PetscCall(PetscFree(prolong_restr_mat)); 838 PetscCall(PetscFree(app_ctx->level_degrees)); 839 PetscCall(PetscFree(ceed_data)); 840 841 // libCEED objects 842 CeedVectorDestroy(&form_jacob_ctx->coo_values); 843 CeedQFunctionContextDestroy(&ctx_phys); 844 CeedQFunctionContextDestroy(&ctx_phys_smoother); 845 CeedDestroy(&ceed); 846 847 // PETSc objects 848 PetscCall(VecDestroy(&U)); 849 PetscCall(VecDestroy(&R)); 850 PetscCall(VecDestroy(&R_loc)); 851 PetscCall(VecDestroy(&F)); 852 PetscCall(VecDestroy(&F_loc)); 853 PetscCall(VecDestroy(&neumann_bcs)); 854 PetscCall(VecDestroy(&bcs_loc)); 855 PetscCall(MatDestroy(&jacob_mat_coarse)); 856 PetscCall(SNESDestroy(&snes)); 857 PetscCall(DMDestroy(&dm_orig)); 858 PetscCall(DMDestroy(&dm_energy)); 859 PetscCall(DMDestroy(&dm_diagnostic)); 860 PetscCall(PetscFree(level_dms)); 861 862 // -- Function list 863 PetscCall(PetscFunctionListDestroy(&problem_functions->setupPhysics)); 864 PetscCall(PetscFunctionListDestroy(&problem_functions->setupSmootherPhysics)); 865 PetscCall(PetscFunctionListDestroy(&problem_functions->setupLibceedFineLevel)); 866 PetscCall(PetscFunctionListDestroy(&problem_functions->setupLibceedLevel)); 867 868 // Structs 869 PetscCall(PetscFree(res_ctx)); 870 PetscCall(PetscFree(form_jacob_ctx)); 871 PetscCall(PetscFree(jacob_coarse_ctx)); 872 PetscCall(PetscFree(app_ctx)); 873 PetscCall(PetscFree(problem_functions)); 874 PetscCall(PetscFree(units)); 875 876 return PetscFinalize(); 877 } 878