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