1 // Copyright (c) 2017-2026, 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
main(int argc,char ** argv)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