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