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 /// Matrix shell operations for solid mechanics example using PETSc
10
11 #include "../include/matops.h"
12
13 #include <ceed.h>
14 #include <petscdmplex.h>
15
16 #include "../include/structs.h"
17 #include "../include/utils.h"
18
19 // -----------------------------------------------------------------------------
20 // libCEED Operators for MatShell
21 // -----------------------------------------------------------------------------
22 // This function uses libCEED to compute the local action of an operator
ApplyLocalCeedOp(Vec X,Vec Y,UserMult user)23 PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, UserMult user) {
24 PetscScalar *x, *y;
25 PetscMemType x_mem_type, y_mem_type;
26
27 PetscFunctionBeginUser;
28
29 // Global-to-local
30 PetscCall(DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->X_loc));
31 PetscCall(VecZeroEntries(user->Y_loc));
32
33 // Setup CEED vectors
34 PetscCall(VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x, &x_mem_type));
35 PetscCall(VecGetArrayAndMemType(user->Y_loc, &y, &y_mem_type));
36 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x);
37 CeedVectorSetArray(user->y_ceed, MemTypeP2C(y_mem_type), CEED_USE_POINTER, y);
38
39 // Apply CEED operator
40 CeedOperatorApply(user->op, user->x_ceed, user->y_ceed, CEED_REQUEST_IMMEDIATE);
41
42 // Restore PETSc vectors
43 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL);
44 CeedVectorTakeArray(user->y_ceed, MemTypeP2C(y_mem_type), NULL);
45 PetscCall(VecRestoreArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x));
46 PetscCall(VecRestoreArrayAndMemType(user->Y_loc, &y));
47
48 // Local-to-global
49 PetscCall(VecZeroEntries(Y));
50 PetscCall(DMLocalToGlobal(user->dm, user->Y_loc, ADD_VALUES, Y));
51
52 PetscFunctionReturn(PETSC_SUCCESS);
53 };
54
55 // This function uses libCEED to compute the non-linear residual
FormResidual_Ceed(SNES snes,Vec X,Vec Y,void * ctx)56 PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx) {
57 UserMult user = (UserMult)ctx;
58
59 PetscFunctionBeginUser;
60
61 // Use computed BCs
62 PetscCall(VecZeroEntries(user->X_loc));
63 PetscCall(DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->X_loc, user->load_increment, NULL, NULL, NULL));
64
65 // libCEED for local action of residual evaluator
66 PetscCall(ApplyLocalCeedOp(X, Y, user));
67
68 // Neumann BCs
69 if (user->neumann_bcs) {
70 PetscCall(VecAXPY(Y, -user->load_increment, user->neumann_bcs));
71 }
72
73 PetscFunctionReturn(PETSC_SUCCESS);
74 };
75
76 // This function uses libCEED to apply the Jacobian for assembly via a SNES
ApplyJacobianCoarse_Ceed(SNES snes,Vec X,Vec Y,void * ctx)77 PetscErrorCode ApplyJacobianCoarse_Ceed(SNES snes, Vec X, Vec Y, void *ctx) {
78 UserMult user = (UserMult)ctx;
79
80 PetscFunctionBeginUser;
81
82 // Zero boundary values
83 PetscCall(VecZeroEntries(user->X_loc));
84
85 // libCEED for local action of residual evaluator
86 PetscCall(ApplyLocalCeedOp(X, Y, user));
87
88 PetscFunctionReturn(PETSC_SUCCESS);
89 };
90
91 // This function uses libCEED to compute the action of the Jacobian
ApplyJacobian_Ceed(Mat A,Vec X,Vec Y)92 PetscErrorCode ApplyJacobian_Ceed(Mat A, Vec X, Vec Y) {
93 UserMult user;
94
95 PetscFunctionBeginUser;
96
97 // Zero boundary values
98 PetscCall(MatShellGetContext(A, &user));
99 PetscCall(VecZeroEntries(user->X_loc));
100
101 // libCEED for local action of Jacobian
102 PetscCall(ApplyLocalCeedOp(X, Y, user));
103
104 PetscFunctionReturn(PETSC_SUCCESS);
105 };
106
107 // This function uses libCEED to compute the action of the prolongation operator
Prolong_Ceed(Mat A,Vec X,Vec Y)108 PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y) {
109 UserMultProlongRestr user;
110 PetscScalar *c, *f;
111 PetscMemType c_mem_type, f_mem_type;
112
113 PetscFunctionBeginUser;
114
115 PetscCall(MatShellGetContext(A, &user));
116
117 // Global-to-local
118 PetscCall(VecZeroEntries(user->loc_vec_c));
119 PetscCall(DMGlobalToLocal(user->dm_c, X, INSERT_VALUES, user->loc_vec_c));
120 PetscCall(VecZeroEntries(user->loc_vec_f));
121
122 // Setup CEED vectors
123 PetscCall(VecGetArrayReadAndMemType(user->loc_vec_c, (const PetscScalar **)&c, &c_mem_type));
124 PetscCall(VecGetArrayAndMemType(user->loc_vec_f, &f, &f_mem_type));
125 CeedVectorSetArray(user->ceed_vec_c, MemTypeP2C(c_mem_type), CEED_USE_POINTER, c);
126 CeedVectorSetArray(user->ceed_vec_f, MemTypeP2C(f_mem_type), CEED_USE_POINTER, f);
127
128 // Apply CEED operator
129 CeedOperatorApply(user->op_prolong, user->ceed_vec_c, user->ceed_vec_f, CEED_REQUEST_IMMEDIATE);
130
131 // Restore PETSc vectors
132 CeedVectorTakeArray(user->ceed_vec_c, MemTypeP2C(c_mem_type), NULL);
133 CeedVectorTakeArray(user->ceed_vec_f, MemTypeP2C(f_mem_type), NULL);
134 PetscCall(VecRestoreArrayReadAndMemType(user->loc_vec_c, (const PetscScalar **)&c));
135 PetscCall(VecRestoreArrayAndMemType(user->loc_vec_f, &f));
136
137 // Local-to-global
138 PetscCall(VecZeroEntries(Y));
139 PetscCall(DMLocalToGlobal(user->dm_f, user->loc_vec_f, ADD_VALUES, Y));
140
141 PetscFunctionReturn(PETSC_SUCCESS);
142 }
143
144 // This function uses libCEED to compute the action of the restriction operator
Restrict_Ceed(Mat A,Vec X,Vec Y)145 PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y) {
146 UserMultProlongRestr user;
147 PetscScalar *c, *f;
148 PetscMemType c_mem_type, f_mem_type;
149
150 PetscFunctionBeginUser;
151
152 PetscCall(MatShellGetContext(A, &user));
153
154 // Global-to-local
155 PetscCall(VecZeroEntries(user->loc_vec_f));
156 PetscCall(DMGlobalToLocal(user->dm_f, X, INSERT_VALUES, user->loc_vec_f));
157 PetscCall(VecZeroEntries(user->loc_vec_c));
158
159 // Setup CEED vectors
160 PetscCall(VecGetArrayReadAndMemType(user->loc_vec_f, (const PetscScalar **)&f, &f_mem_type));
161 PetscCall(VecGetArrayAndMemType(user->loc_vec_c, &c, &c_mem_type));
162 CeedVectorSetArray(user->ceed_vec_f, MemTypeP2C(f_mem_type), CEED_USE_POINTER, f);
163 CeedVectorSetArray(user->ceed_vec_c, MemTypeP2C(c_mem_type), CEED_USE_POINTER, c);
164
165 // Apply CEED operator
166 CeedOperatorApply(user->op_restrict, user->ceed_vec_f, user->ceed_vec_c, CEED_REQUEST_IMMEDIATE);
167
168 // Restore PETSc vectors
169 CeedVectorTakeArray(user->ceed_vec_f, MemTypeP2C(f_mem_type), NULL);
170 CeedVectorTakeArray(user->ceed_vec_c, MemTypeP2C(c_mem_type), NULL);
171 PetscCall(VecRestoreArrayReadAndMemType(user->loc_vec_f, (const PetscScalar **)&f));
172 PetscCall(VecRestoreArrayAndMemType(user->loc_vec_c, &c));
173
174 // Local-to-global
175 PetscCall(VecZeroEntries(Y));
176 PetscCall(DMLocalToGlobal(user->dm_c, user->loc_vec_c, ADD_VALUES, Y));
177
178 PetscFunctionReturn(PETSC_SUCCESS);
179 };
180
181 // This function returns the computed diagonal of the operator
GetDiag_Ceed(Mat A,Vec D)182 PetscErrorCode GetDiag_Ceed(Mat A, Vec D) {
183 UserMult user;
184
185 PetscFunctionBeginUser;
186
187 PetscCall(MatShellGetContext(A, &user));
188
189 // -- Set physics context
190 if (user->ctx_phys_smoother) CeedQFunctionSetContext(user->qf, user->ctx_phys_smoother);
191
192 // Compute Diagonal via libCEED
193 PetscScalar *x;
194 PetscMemType x_mem_type;
195
196 // -- Place PETSc vector in libCEED vector
197 PetscCall(VecGetArrayAndMemType(user->X_loc, &x, &x_mem_type));
198 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x);
199
200 // -- Compute Diagonal
201 CeedOperatorLinearAssembleDiagonal(user->op, user->x_ceed, CEED_REQUEST_IMMEDIATE);
202
203 // -- Reset physics context
204 if (user->ctx_phys_smoother) CeedQFunctionSetContext(user->qf, user->ctx_phys);
205
206 // -- Local-to-Global
207 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL);
208 PetscCall(VecRestoreArrayAndMemType(user->X_loc, &x));
209 PetscCall(VecZeroEntries(D));
210 PetscCall(DMLocalToGlobal(user->dm, user->X_loc, ADD_VALUES, D));
211
212 // Cleanup
213 PetscCall(VecZeroEntries(user->X_loc));
214
215 PetscFunctionReturn(PETSC_SUCCESS);
216 };
217
218 // This function calculates the strain energy in the final solution
ComputeStrainEnergy(DM dmEnergy,UserMult user,CeedOperator op_energy,Vec X,PetscReal * energy)219 PetscErrorCode ComputeStrainEnergy(DM dmEnergy, UserMult user, CeedOperator op_energy, Vec X, PetscReal *energy) {
220 PetscScalar *x;
221 PetscMemType x_mem_type;
222 PetscInt length;
223
224 PetscFunctionBeginUser;
225
226 // Global-to-local
227 PetscCall(VecZeroEntries(user->X_loc));
228 PetscCall(DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->X_loc));
229 PetscCall(DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->X_loc, user->load_increment, NULL, NULL, NULL));
230
231 // Setup libCEED input vector
232 PetscCall(VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x, &x_mem_type));
233 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x);
234
235 // Setup libCEED output vector
236 Vec E_loc;
237 CeedVector e_loc;
238
239 PetscCall(DMCreateLocalVector(dmEnergy, &E_loc));
240 PetscCall(VecGetSize(E_loc, &length));
241 PetscCall(VecDestroy(&E_loc));
242 CeedVectorCreate(user->ceed, length, &e_loc);
243
244 // Apply libCEED operator
245 CeedOperatorApply(op_energy, user->x_ceed, e_loc, CEED_REQUEST_IMMEDIATE);
246
247 // Restore PETSc vector
248 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL);
249 PetscCall(VecRestoreArrayRead(user->X_loc, (const PetscScalar **)&x));
250
251 // Reduce max error
252 const CeedScalar *e;
253 CeedVectorGetArrayRead(e_loc, CEED_MEM_HOST, &e);
254 (*energy) = 0;
255 for (CeedInt i = 0; i < length; i++) (*energy) += e[i];
256 CeedVectorRestoreArrayRead(e_loc, &e);
257 CeedVectorDestroy(&e_loc);
258
259 PetscCall(MPI_Allreduce(MPI_IN_PLACE, energy, 1, MPIU_REAL, MPIU_SUM, user->comm));
260
261 PetscFunctionReturn(PETSC_SUCCESS);
262 };
263