matops.c (25a77dc8a1b03301bfbe493d597d7d9d199b7a5e) matops.c (b68a8d799acb1d44569fb95028e25f895284bc04)
1// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
2// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
3// reserved. See files LICENSE and NOTICE for details.
4//
5// This file is part of CEED, a collection of benchmarks, miniapps, software
6// libraries and APIs for efficient high-order finite element and spectral
7// element discretizations for exascale applications. For more information and
8// source code availability see http://github.com/ceed.

--- 12 unchanged lines hidden (view full) ---

21
22// -----------------------------------------------------------------------------
23// libCEED Operators for MatShell
24// -----------------------------------------------------------------------------
25// This function uses libCEED to compute the local action of an operator
26PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, UserMult user) {
27 PetscErrorCode ierr;
28 PetscScalar *x, *y;
1// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
2// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
3// reserved. See files LICENSE and NOTICE for details.
4//
5// This file is part of CEED, a collection of benchmarks, miniapps, software
6// libraries and APIs for efficient high-order finite element and spectral
7// element discretizations for exascale applications. For more information and
8// source code availability see http://github.com/ceed.

--- 12 unchanged lines hidden (view full) ---

21
22// -----------------------------------------------------------------------------
23// libCEED Operators for MatShell
24// -----------------------------------------------------------------------------
25// This function uses libCEED to compute the local action of an operator
26PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, UserMult user) {
27 PetscErrorCode ierr;
28 PetscScalar *x, *y;
29 PetscMemType xmemtype, ymemtype;
29
30 PetscFunctionBeginUser;
31
32 // Global-to-local
33 ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->Xloc); CHKERRQ(ierr);
34 ierr = VecZeroEntries(user->Yloc); CHKERRQ(ierr);
35
36 // Setup CEED vectors
30
31 PetscFunctionBeginUser;
32
33 // Global-to-local
34 ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->Xloc); CHKERRQ(ierr);
35 ierr = VecZeroEntries(user->Yloc); CHKERRQ(ierr);
36
37 // Setup CEED vectors
37 ierr = user->VecGetArrayRead(user->Xloc, (const PetscScalar **)&x);
38 CHKERRQ(ierr);
39 ierr = user->VecGetArray(user->Yloc, &y); CHKERRQ(ierr);
40 CeedVectorSetArray(user->Xceed, user->memType, CEED_USE_POINTER, x);
41 CeedVectorSetArray(user->Yceed, user->memType, CEED_USE_POINTER, y);
38 ierr = VecGetArrayReadAndMemType(user->Xloc, (const PetscScalar **)&x,
39 &xmemtype); CHKERRQ(ierr);
40 ierr = VecGetArrayAndMemType(user->Yloc, &y, &ymemtype); CHKERRQ(ierr);
41 CeedVectorSetArray(user->Xceed, MemTypeP2C(xmemtype), CEED_USE_POINTER, x);
42 CeedVectorSetArray(user->Yceed, MemTypeP2C(ymemtype), CEED_USE_POINTER, y);
42
43 // Apply CEED operator
43
44 // Apply CEED operator
44 // Note: We could use VecGetArrayInPlace. Instead, we use SetArray/TakeArray
45 // so we can request host memory for easier debugging.
46 CeedOperatorApply(user->op, user->Xceed, user->Yceed, CEED_REQUEST_IMMEDIATE);
47
48 // Restore PETSc vectors
45 CeedOperatorApply(user->op, user->Xceed, user->Yceed, CEED_REQUEST_IMMEDIATE);
46
47 // Restore PETSc vectors
49 CeedVectorTakeArray(user->Xceed, user->memType, NULL);
50 CeedVectorTakeArray(user->Yceed, user->memType, NULL);
51 ierr = user->VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x);
48 CeedVectorTakeArray(user->Xceed, MemTypeP2C(xmemtype), NULL);
49 CeedVectorTakeArray(user->Yceed, MemTypeP2C(ymemtype), NULL);
50 ierr = VecRestoreArrayReadAndMemType(user->Xloc, (const PetscScalar **)&x);
52 CHKERRQ(ierr);
51 CHKERRQ(ierr);
53 ierr = user->VecRestoreArray(user->Yloc, &y); CHKERRQ(ierr);
52 ierr = VecRestoreArrayAndMemType(user->Yloc, &y); CHKERRQ(ierr);
54
55 // Local-to-global
56 ierr = VecZeroEntries(Y); CHKERRQ(ierr);
57 ierr = DMLocalToGlobal(user->dm, user->Yloc, ADD_VALUES, Y); CHKERRQ(ierr);
58
59 PetscFunctionReturn(0);
60};
61

--- 54 unchanged lines hidden (view full) ---

116 PetscFunctionReturn(0);
117};
118
119// This function uses libCEED to compute the action of the prolongation operator
120PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y) {
121 PetscErrorCode ierr;
122 UserMultProlongRestr user;
123 PetscScalar *c, *f;
53
54 // Local-to-global
55 ierr = VecZeroEntries(Y); CHKERRQ(ierr);
56 ierr = DMLocalToGlobal(user->dm, user->Yloc, ADD_VALUES, Y); CHKERRQ(ierr);
57
58 PetscFunctionReturn(0);
59};
60

--- 54 unchanged lines hidden (view full) ---

115 PetscFunctionReturn(0);
116};
117
118// This function uses libCEED to compute the action of the prolongation operator
119PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y) {
120 PetscErrorCode ierr;
121 UserMultProlongRestr user;
122 PetscScalar *c, *f;
123 PetscMemType cmemtype, fmemtype;
124
125 PetscFunctionBeginUser;
126
127 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr);
128
129 // Global-to-local
130 ierr = VecZeroEntries(user->locVecC); CHKERRQ(ierr);
131 ierr = DMGlobalToLocal(user->dmC, X, INSERT_VALUES, user->locVecC);
132 CHKERRQ(ierr);
133 ierr = VecZeroEntries(user->locVecF); CHKERRQ(ierr);
134
135 // Setup CEED vectors
124
125 PetscFunctionBeginUser;
126
127 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr);
128
129 // Global-to-local
130 ierr = VecZeroEntries(user->locVecC); CHKERRQ(ierr);
131 ierr = DMGlobalToLocal(user->dmC, X, INSERT_VALUES, user->locVecC);
132 CHKERRQ(ierr);
133 ierr = VecZeroEntries(user->locVecF); CHKERRQ(ierr);
134
135 // Setup CEED vectors
136 ierr = user->VecGetArrayRead(user->locVecC, (const PetscScalar **)&c);
137 CHKERRQ(ierr);
138 ierr = user->VecGetArray(user->locVecF, &f); CHKERRQ(ierr);
139 CeedVectorSetArray(user->ceedVecC, user->memType, CEED_USE_POINTER, c);
140 CeedVectorSetArray(user->ceedVecF, user->memType, CEED_USE_POINTER, f);
136 ierr = VecGetArrayReadAndMemType(user->locVecC, (const PetscScalar **)&c,
137 &cmemtype); CHKERRQ(ierr);
138 ierr = VecGetArrayAndMemType(user->locVecF, &f, &fmemtype); CHKERRQ(ierr);
139 CeedVectorSetArray(user->ceedVecC, MemTypeP2C(cmemtype), CEED_USE_POINTER, c);
140 CeedVectorSetArray(user->ceedVecF, MemTypeP2C(fmemtype), CEED_USE_POINTER, f);
141
142 // Apply CEED operator
143 CeedOperatorApply(user->opProlong, user->ceedVecC, user->ceedVecF,
144 CEED_REQUEST_IMMEDIATE);
145
146 // Restore PETSc vectors
141
142 // Apply CEED operator
143 CeedOperatorApply(user->opProlong, user->ceedVecC, user->ceedVecF,
144 CEED_REQUEST_IMMEDIATE);
145
146 // Restore PETSc vectors
147 CeedVectorTakeArray(user->ceedVecC, user->memType, NULL);
148 CeedVectorTakeArray(user->ceedVecF, user->memType, NULL);
149 ierr = user->VecRestoreArrayRead(user->locVecC, (const PetscScalar **)&c);
147 CeedVectorTakeArray(user->ceedVecC, MemTypeP2C(cmemtype), NULL);
148 CeedVectorTakeArray(user->ceedVecF, MemTypeP2C(fmemtype), NULL);
149 ierr = VecRestoreArrayReadAndMemType(user->locVecC, (const PetscScalar **)&c);
150 CHKERRQ(ierr);
150 CHKERRQ(ierr);
151 ierr = user->VecRestoreArray(user->locVecF, &f); CHKERRQ(ierr);
151 ierr = VecRestoreArrayAndMemType(user->locVecF, &f); CHKERRQ(ierr);
152
153 // Local-to-global
154 ierr = VecZeroEntries(Y); CHKERRQ(ierr);
155 ierr = DMLocalToGlobal(user->dmF, user->locVecF, ADD_VALUES, Y);
156 CHKERRQ(ierr);
157
158 PetscFunctionReturn(0);
159}
160
161// This function uses libCEED to compute the action of the restriction operator
162PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y) {
163 PetscErrorCode ierr;
164 UserMultProlongRestr user;
165 PetscScalar *c, *f;
152
153 // Local-to-global
154 ierr = VecZeroEntries(Y); CHKERRQ(ierr);
155 ierr = DMLocalToGlobal(user->dmF, user->locVecF, ADD_VALUES, Y);
156 CHKERRQ(ierr);
157
158 PetscFunctionReturn(0);
159}
160
161// This function uses libCEED to compute the action of the restriction operator
162PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y) {
163 PetscErrorCode ierr;
164 UserMultProlongRestr user;
165 PetscScalar *c, *f;
166 PetscMemType cmemtype, fmemtype;
166
167 PetscFunctionBeginUser;
168
169 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr);
170
171 // Global-to-local
172 ierr = VecZeroEntries(user->locVecF); CHKERRQ(ierr);
173 ierr = DMGlobalToLocal(user->dmF, X, INSERT_VALUES, user->locVecF);
174 CHKERRQ(ierr);
175 ierr = VecZeroEntries(user->locVecC); CHKERRQ(ierr);
176
177 // Setup CEED vectors
167
168 PetscFunctionBeginUser;
169
170 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr);
171
172 // Global-to-local
173 ierr = VecZeroEntries(user->locVecF); CHKERRQ(ierr);
174 ierr = DMGlobalToLocal(user->dmF, X, INSERT_VALUES, user->locVecF);
175 CHKERRQ(ierr);
176 ierr = VecZeroEntries(user->locVecC); CHKERRQ(ierr);
177
178 // Setup CEED vectors
178 ierr = user->VecGetArrayRead(user->locVecF, (const PetscScalar **)&f);
179 CHKERRQ(ierr);
180 ierr = user->VecGetArray(user->locVecC, &c); CHKERRQ(ierr);
181 CeedVectorSetArray(user->ceedVecF, user->memType, CEED_USE_POINTER, f);
182 CeedVectorSetArray(user->ceedVecC, user->memType, CEED_USE_POINTER, c);
179 ierr = VecGetArrayReadAndMemType(user->locVecF, (const PetscScalar **)&f,
180 &fmemtype); CHKERRQ(ierr);
181 ierr = VecGetArrayAndMemType(user->locVecC, &c, &cmemtype); CHKERRQ(ierr);
182 CeedVectorSetArray(user->ceedVecF, MemTypeP2C(fmemtype), CEED_USE_POINTER, f);
183 CeedVectorSetArray(user->ceedVecC, MemTypeP2C(cmemtype), CEED_USE_POINTER, c);
183
184 // Apply CEED operator
185 CeedOperatorApply(user->opRestrict, user->ceedVecF, user->ceedVecC,
186 CEED_REQUEST_IMMEDIATE);
187
188 // Restore PETSc vectors
184
185 // Apply CEED operator
186 CeedOperatorApply(user->opRestrict, user->ceedVecF, user->ceedVecC,
187 CEED_REQUEST_IMMEDIATE);
188
189 // Restore PETSc vectors
189 CeedVectorTakeArray(user->ceedVecF, user->memType, NULL);
190 CeedVectorTakeArray(user->ceedVecC, user->memType, NULL);
191 ierr = user->VecRestoreArrayRead(user->locVecF, (const PetscScalar **)&f);
190 CeedVectorTakeArray(user->ceedVecF, MemTypeP2C(fmemtype), NULL);
191 CeedVectorTakeArray(user->ceedVecC, MemTypeP2C(cmemtype), NULL);
192 ierr = VecRestoreArrayReadAndMemType(user->locVecF, (const PetscScalar **)&f);
192 CHKERRQ(ierr);
193 CHKERRQ(ierr);
193 ierr = user->VecRestoreArray(user->locVecC, &c); CHKERRQ(ierr);
194 ierr = VecRestoreArrayAndMemType(user->locVecC, &c); CHKERRQ(ierr);
194
195 // Local-to-global
196 ierr = VecZeroEntries(Y); CHKERRQ(ierr);
197 ierr = DMLocalToGlobal(user->dmC, user->locVecC, ADD_VALUES, Y);
198 CHKERRQ(ierr);
199
200 PetscFunctionReturn(0);
201};

--- 8 unchanged lines hidden (view full) ---

210 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr);
211
212 // -- Set physics context
213 if (user->ctxPhysSmoother)
214 CeedQFunctionSetContext(user->qf, user->ctxPhysSmoother);
215
216 // Compute Diagonal via libCEED
217 PetscScalar *x;
195
196 // Local-to-global
197 ierr = VecZeroEntries(Y); CHKERRQ(ierr);
198 ierr = DMLocalToGlobal(user->dmC, user->locVecC, ADD_VALUES, Y);
199 CHKERRQ(ierr);
200
201 PetscFunctionReturn(0);
202};

--- 8 unchanged lines hidden (view full) ---

211 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr);
212
213 // -- Set physics context
214 if (user->ctxPhysSmoother)
215 CeedQFunctionSetContext(user->qf, user->ctxPhysSmoother);
216
217 // Compute Diagonal via libCEED
218 PetscScalar *x;
219 PetscMemType xmemtype;
218
219 // -- Place PETSc vector in libCEED vector
220
221 // -- Place PETSc vector in libCEED vector
220 ierr = user->VecGetArray(user->Xloc, &x); CHKERRQ(ierr);
221 CeedVectorSetArray(user->Xceed, user->memType, CEED_USE_POINTER, x);
222 ierr = VecGetArrayAndMemType(user->Xloc, &x, &xmemtype); CHKERRQ(ierr);
223 CeedVectorSetArray(user->Xceed, MemTypeP2C(xmemtype), CEED_USE_POINTER, x);
222
223 // -- Compute Diagonal
224 CeedOperatorLinearAssembleDiagonal(user->op, user->Xceed,
225 CEED_REQUEST_IMMEDIATE);
226
227 // -- Reset physics context
228 if (user->ctxPhysSmoother)
229 CeedQFunctionSetContext(user->qf, user->ctxPhys);
230
231 // -- Local-to-Global
224
225 // -- Compute Diagonal
226 CeedOperatorLinearAssembleDiagonal(user->op, user->Xceed,
227 CEED_REQUEST_IMMEDIATE);
228
229 // -- Reset physics context
230 if (user->ctxPhysSmoother)
231 CeedQFunctionSetContext(user->qf, user->ctxPhys);
232
233 // -- Local-to-Global
232 CeedVectorTakeArray(user->Xceed, user->memType, NULL);
233 ierr = user->VecRestoreArray(user->Xloc, &x); CHKERRQ(ierr);
234 CeedVectorTakeArray(user->Xceed, MemTypeP2C(xmemtype), NULL);
235 ierr = VecRestoreArrayAndMemType(user->Xloc, &x); CHKERRQ(ierr);
234 ierr = VecZeroEntries(D); CHKERRQ(ierr);
235 ierr = DMLocalToGlobal(user->dm, user->Xloc, ADD_VALUES, D); CHKERRQ(ierr);
236
237 // Cleanup
238 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr);
239
240 PetscFunctionReturn(0);
241};
242
243// This function calculates the strain energy in the final solution
244PetscErrorCode ComputeStrainEnergy(DM dmEnergy, UserMult user,
245 CeedOperator opEnergy, Vec X,
246 PetscReal *energy) {
247 PetscErrorCode ierr;
248 PetscScalar *x;
236 ierr = VecZeroEntries(D); CHKERRQ(ierr);
237 ierr = DMLocalToGlobal(user->dm, user->Xloc, ADD_VALUES, D); CHKERRQ(ierr);
238
239 // Cleanup
240 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr);
241
242 PetscFunctionReturn(0);
243};
244
245// This function calculates the strain energy in the final solution
246PetscErrorCode ComputeStrainEnergy(DM dmEnergy, UserMult user,
247 CeedOperator opEnergy, Vec X,
248 PetscReal *energy) {
249 PetscErrorCode ierr;
250 PetscScalar *x;
251 PetscMemType xmemtype;
249 CeedInt length;
250
251 PetscFunctionBeginUser;
252
253 // Global-to-local
254 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr);
255 ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->Xloc); CHKERRQ(ierr);
256 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->Xloc,
257 user->loadIncrement, NULL, NULL, NULL);
258 CHKERRQ(ierr);
259
260 // Setup libCEED input vector
252 CeedInt length;
253
254 PetscFunctionBeginUser;
255
256 // Global-to-local
257 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr);
258 ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->Xloc); CHKERRQ(ierr);
259 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->Xloc,
260 user->loadIncrement, NULL, NULL, NULL);
261 CHKERRQ(ierr);
262
263 // Setup libCEED input vector
261 ierr = user->VecGetArrayRead(user->Xloc, (const PetscScalar **)&x);
262 CHKERRQ(ierr);
263 CeedVectorSetArray(user->Xceed, user->memType, CEED_USE_POINTER, x);
264 ierr = VecGetArrayReadAndMemType(user->Xloc, (const PetscScalar **)&x,
265 &xmemtype); CHKERRQ(ierr);
266 CeedVectorSetArray(user->Xceed, MemTypeP2C(xmemtype), CEED_USE_POINTER, x);
264
265 // Setup libCEED output vector
266 Vec Eloc;
267 CeedVector eloc;
268 ierr = DMCreateLocalVector(dmEnergy, &Eloc); CHKERRQ(ierr);
269 ierr = VecGetSize(Eloc, &length); CHKERRQ(ierr);
270 ierr = VecDestroy(&Eloc); CHKERRQ(ierr);
271 CeedVectorCreate(user->ceed, length, &eloc);
272
273 // Apply libCEED operator
274 CeedOperatorApply(opEnergy, user->Xceed, eloc, CEED_REQUEST_IMMEDIATE);
275
276 // Restore PETSc vector
267
268 // Setup libCEED output vector
269 Vec Eloc;
270 CeedVector eloc;
271 ierr = DMCreateLocalVector(dmEnergy, &Eloc); CHKERRQ(ierr);
272 ierr = VecGetSize(Eloc, &length); CHKERRQ(ierr);
273 ierr = VecDestroy(&Eloc); CHKERRQ(ierr);
274 CeedVectorCreate(user->ceed, length, &eloc);
275
276 // Apply libCEED operator
277 CeedOperatorApply(opEnergy, user->Xceed, eloc, CEED_REQUEST_IMMEDIATE);
278
279 // Restore PETSc vector
277 ierr = user->VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x);
280 CeedVectorTakeArray(user->Xceed, MemTypeP2C(xmemtype), NULL);
281 ierr = VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x);
278 CHKERRQ(ierr);
279
280 // Reduce max error
281 const CeedScalar *e;
282 CeedVectorGetArrayRead(eloc, CEED_MEM_HOST, &e);
283 (*energy) = 0;
284 for (CeedInt i=0; i<length; i++)
285 (*energy) += e[i];
286 CeedVectorRestoreArrayRead(eloc, &e);
287 CeedVectorDestroy(&eloc);
288
289 ierr = MPI_Allreduce(MPI_IN_PLACE, energy, 1, MPIU_REAL, MPIU_SUM,
290 user->comm); CHKERRQ(ierr);
291
292 PetscFunctionReturn(0);
293};
282 CHKERRQ(ierr);
283
284 // Reduce max error
285 const CeedScalar *e;
286 CeedVectorGetArrayRead(eloc, CEED_MEM_HOST, &e);
287 (*energy) = 0;
288 for (CeedInt i=0; i<length; i++)
289 (*energy) += e[i];
290 CeedVectorRestoreArrayRead(eloc, &e);
291 CeedVectorDestroy(&eloc);
292
293 ierr = MPI_Allreduce(MPI_IN_PLACE, energy, 1, MPIU_REAL, MPIU_SUM,
294 user->comm); CHKERRQ(ierr);
295
296 PetscFunctionReturn(0);
297};