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 /// @file
9 /// Helper functions for solid mechanics example using PETSc
10
11 #include "../include/misc.h"
12
13 #include <ceed.h>
14 #include <petscdmplex.h>
15 #include <petscmat.h>
16
17 #include "../include/utils.h"
18
19 // -----------------------------------------------------------------------------
20 // Create libCEED operator context
21 // -----------------------------------------------------------------------------
22 // Setup context data for Jacobian evaluation
SetupJacobianCtx(MPI_Comm comm,AppCtx app_ctx,DM dm,Vec V,Vec V_loc,CeedData ceed_data,Ceed ceed,CeedQFunctionContext ctx_phys,CeedQFunctionContext ctx_phys_smoother,UserMult jacobian_ctx)23 PetscErrorCode SetupJacobianCtx(MPI_Comm comm, AppCtx app_ctx, DM dm, Vec V, Vec V_loc, CeedData ceed_data, Ceed ceed, CeedQFunctionContext ctx_phys,
24 CeedQFunctionContext ctx_phys_smoother, UserMult jacobian_ctx) {
25 PetscFunctionBeginUser;
26
27 // PETSc objects
28 jacobian_ctx->comm = comm;
29 jacobian_ctx->dm = dm;
30
31 // Work vectors
32 jacobian_ctx->X_loc = V_loc;
33 PetscCall(VecDuplicate(V_loc, &jacobian_ctx->Y_loc));
34 jacobian_ctx->x_ceed = ceed_data->x_ceed;
35 jacobian_ctx->y_ceed = ceed_data->y_ceed;
36
37 // libCEED operator
38 jacobian_ctx->op = ceed_data->op_jacobian;
39 jacobian_ctx->qf = ceed_data->qf_jacobian;
40
41 // Ceed
42 jacobian_ctx->ceed = ceed;
43
44 // Physics
45 jacobian_ctx->ctx_phys = ctx_phys;
46 jacobian_ctx->ctx_phys_smoother = ctx_phys_smoother;
47 PetscFunctionReturn(PETSC_SUCCESS);
48 };
49
50 // Setup context data for prolongation and restriction operators
SetupProlongRestrictCtx(MPI_Comm comm,AppCtx app_ctx,DM dm_c,DM dm_f,Vec V_f,Vec V_loc_c,Vec V_loc_f,CeedData ceed_data_c,CeedData ceed_data_f,Ceed ceed,UserMultProlongRestr prolong_restr_ctx)51 PetscErrorCode SetupProlongRestrictCtx(MPI_Comm comm, AppCtx app_ctx, DM dm_c, DM dm_f, Vec V_f, Vec V_loc_c, Vec V_loc_f, CeedData ceed_data_c,
52 CeedData ceed_data_f, Ceed ceed, UserMultProlongRestr prolong_restr_ctx) {
53 PetscFunctionBeginUser;
54
55 // PETSc objects
56 prolong_restr_ctx->comm = comm;
57 prolong_restr_ctx->dm_c = dm_c;
58 prolong_restr_ctx->dm_f = dm_f;
59
60 // Work vectors
61 prolong_restr_ctx->loc_vec_c = V_loc_c;
62 prolong_restr_ctx->loc_vec_f = V_loc_f;
63 prolong_restr_ctx->ceed_vec_c = ceed_data_c->x_ceed;
64 prolong_restr_ctx->ceed_vec_f = ceed_data_f->x_ceed;
65
66 // libCEED operators
67 prolong_restr_ctx->op_prolong = ceed_data_f->op_prolong;
68 prolong_restr_ctx->op_restrict = ceed_data_f->op_restrict;
69
70 // Ceed
71 prolong_restr_ctx->ceed = ceed;
72 PetscFunctionReturn(PETSC_SUCCESS);
73 };
74
75 // -----------------------------------------------------------------------------
76 // Jacobian setup
77 // -----------------------------------------------------------------------------
FormJacobian(SNES snes,Vec U,Mat J,Mat J_pre,void * ctx)78 PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat J_pre, void *ctx) {
79 PetscFunctionBeginUser;
80
81 // Context data
82 FormJacobCtx form_jacob_ctx = (FormJacobCtx)ctx;
83 PetscInt num_levels = form_jacob_ctx->num_levels;
84 Mat *jacob_mat = form_jacob_ctx->jacob_mat;
85
86 // Update Jacobian on each level
87 for (PetscInt level = 0; level < num_levels; level++) {
88 PetscCall(MatAssemblyBegin(jacob_mat[level], MAT_FINAL_ASSEMBLY));
89 PetscCall(MatAssemblyEnd(jacob_mat[level], MAT_FINAL_ASSEMBLY));
90 }
91
92 // Form coarse assembled matrix
93 CeedOperatorLinearAssemble(form_jacob_ctx->op_coarse, form_jacob_ctx->coo_values);
94 const CeedScalar *values;
95 CeedVectorGetArrayRead(form_jacob_ctx->coo_values, CEED_MEM_HOST, &values);
96 PetscCall(MatSetValuesCOO(form_jacob_ctx->jacob_mat_coarse, values, ADD_VALUES));
97 CeedVectorRestoreArrayRead(form_jacob_ctx->coo_values, &values);
98
99 // J_pre might be AIJ (e.g., when using coloring), so we need to assemble it
100 PetscCall(MatAssemblyBegin(J_pre, MAT_FINAL_ASSEMBLY));
101 PetscCall(MatAssemblyEnd(J_pre, MAT_FINAL_ASSEMBLY));
102 if (J != J_pre) {
103 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
104 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
105 }
106 PetscFunctionReturn(PETSC_SUCCESS);
107 };
108
109 // -----------------------------------------------------------------------------
110 // Output solution for visualization
111 // -----------------------------------------------------------------------------
ViewSolution(MPI_Comm comm,AppCtx app_ctx,Vec U,PetscInt increment,PetscScalar load_increment)112 PetscErrorCode ViewSolution(MPI_Comm comm, AppCtx app_ctx, Vec U, PetscInt increment, PetscScalar load_increment) {
113 DM dm;
114 PetscViewer viewer;
115 char output_filename[PETSC_MAX_PATH_LEN];
116 PetscMPIInt rank;
117
118 PetscFunctionBeginUser;
119
120 // Create output directory
121 MPI_Comm_rank(comm, &rank);
122 if (!rank) {
123 PetscCall(PetscMkdir(app_ctx->output_dir));
124 }
125
126 // Build file name
127 PetscCall(PetscSNPrintf(output_filename, sizeof output_filename, "%s/solution-%03" PetscInt_FMT ".vtu", app_ctx->output_dir, increment));
128
129 // Increment sequence
130 PetscCall(VecGetDM(U, &dm));
131 PetscCall(DMSetOutputSequenceNumber(dm, increment, load_increment));
132
133 // Output solution vector
134 PetscCall(PetscViewerVTKOpen(comm, output_filename, FILE_MODE_WRITE, &viewer));
135 PetscCall(VecView(U, viewer));
136 PetscCall(PetscViewerDestroy(&viewer));
137
138 PetscFunctionReturn(PETSC_SUCCESS);
139 };
140
141 // -----------------------------------------------------------------------------
142 // Output diagnostic quantities for visualization
143 // -----------------------------------------------------------------------------
ViewDiagnosticQuantities(MPI_Comm comm,DM dmU,UserMult user,AppCtx app_ctx,Vec U,CeedElemRestriction elem_restr_diagnostic)144 PetscErrorCode ViewDiagnosticQuantities(MPI_Comm comm, DM dmU, UserMult user, AppCtx app_ctx, Vec U, CeedElemRestriction elem_restr_diagnostic) {
145 Vec Diagnostic, Y_loc, mult_vec;
146 CeedVector y_ceed;
147 CeedScalar *x, *y;
148 PetscMemType x_mem_type, y_mem_type;
149 PetscInt loc_size;
150 PetscViewer viewer;
151 char output_filename[PETSC_MAX_PATH_LEN];
152
153 PetscFunctionBeginUser;
154
155 // ---------------------------------------------------------------------------
156 // PETSc and libCEED vectors
157 // ---------------------------------------------------------------------------
158 PetscCall(DMCreateGlobalVector(user->dm, &Diagnostic));
159 PetscCall(PetscObjectSetName((PetscObject)Diagnostic, ""));
160 PetscCall(DMCreateLocalVector(user->dm, &Y_loc));
161 PetscCall(VecGetSize(Y_loc, &loc_size));
162 CeedVectorCreate(user->ceed, loc_size, &y_ceed);
163
164 // ---------------------------------------------------------------------------
165 // Compute quantities
166 // ---------------------------------------------------------------------------
167 // -- Global-to-local
168 PetscCall(VecZeroEntries(user->X_loc));
169 PetscCall(DMPlexInsertBoundaryValues(dmU, PETSC_TRUE, user->X_loc, user->load_increment, NULL, NULL, NULL));
170 PetscCall(DMGlobalToLocal(dmU, U, INSERT_VALUES, user->X_loc));
171 PetscCall(VecZeroEntries(Y_loc));
172
173 // -- Setup CEED vectors
174 PetscCall(VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x, &x_mem_type));
175 PetscCall(VecGetArrayAndMemType(Y_loc, &y, &y_mem_type));
176 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x);
177 CeedVectorSetArray(y_ceed, MemTypeP2C(y_mem_type), CEED_USE_POINTER, y);
178
179 // -- Apply CEED operator
180 CeedOperatorApply(user->op, user->x_ceed, y_ceed, CEED_REQUEST_IMMEDIATE);
181
182 // -- Restore PETSc vector; keep y_ceed viewing memory of Y_loc for use below
183 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL);
184 PetscCall(VecRestoreArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x));
185
186 // -- Local-to-global
187 PetscCall(VecZeroEntries(Diagnostic));
188 PetscCall(DMLocalToGlobal(user->dm, Y_loc, ADD_VALUES, Diagnostic));
189
190 // ---------------------------------------------------------------------------
191 // Scale for multiplicity
192 // ---------------------------------------------------------------------------
193 // -- Setup vectors
194 PetscCall(VecDuplicate(Diagnostic, &mult_vec));
195 PetscCall(VecZeroEntries(Y_loc));
196
197 // -- Compute multiplicity
198 CeedElemRestrictionGetMultiplicity(elem_restr_diagnostic, y_ceed);
199
200 // -- Restore vectors
201 CeedVectorTakeArray(y_ceed, MemTypeP2C(y_mem_type), NULL);
202 PetscCall(VecRestoreArrayAndMemType(Y_loc, &y));
203
204 // -- Local-to-global
205 PetscCall(VecZeroEntries(mult_vec));
206 PetscCall(DMLocalToGlobal(user->dm, Y_loc, ADD_VALUES, mult_vec));
207
208 // -- Scale
209 PetscCall(VecReciprocal(mult_vec));
210 PetscCall(VecPointwiseMult(Diagnostic, Diagnostic, mult_vec));
211
212 // ---------------------------------------------------------------------------
213 // Output solution vector
214 // ---------------------------------------------------------------------------
215 PetscCall(PetscSNPrintf(output_filename, sizeof output_filename, "%s/diagnostic_quantities.vtu", app_ctx->output_dir));
216 PetscCall(PetscViewerVTKOpen(comm, output_filename, FILE_MODE_WRITE, &viewer));
217 PetscCall(VecView(Diagnostic, viewer));
218 PetscCall(PetscViewerDestroy(&viewer));
219
220 // ---------------------------------------------------------------------------
221 // Cleanup
222 // ---------------------------------------------------------------------------
223 PetscCall(VecDestroy(&Diagnostic));
224 PetscCall(VecDestroy(&mult_vec));
225 PetscCall(VecDestroy(&Y_loc));
226 CeedVectorDestroy(&y_ceed);
227
228 PetscFunctionReturn(PETSC_SUCCESS);
229 };
230
231 // -----------------------------------------------------------------------------
232 // Regression testing
233 // -----------------------------------------------------------------------------
234 // test option change. could remove the loading step. Run only with one loading step and compare relatively to ref file
235 // option: expect_final_strain_energy and check against the relative error to ref is within tolerance (10^-5) I.e. one Newton solve then check final
236 // energy
RegressionTests_solids(AppCtx app_ctx,PetscReal energy)237 PetscErrorCode RegressionTests_solids(AppCtx app_ctx, PetscReal energy) {
238 PetscFunctionBegin;
239
240 if (app_ctx->expect_final_strain >= 0.) {
241 PetscReal energy_ref = app_ctx->expect_final_strain;
242 PetscReal error = PetscAbsReal(energy - energy_ref) / energy_ref;
243
244 if (error > app_ctx->test_tol) {
245 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Energy %e does not match expected energy %e: relative tolerance %e > %e\n", (double)energy,
246 (double)energy_ref, (double)error, app_ctx->test_tol));
247 }
248 }
249 PetscFunctionReturn(PETSC_SUCCESS);
250 };
251