xref: /libCEED/examples/solids/src/misc.c (revision 1dfeea39f8876cb55d53d6091f649d2e2329df3d) 
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.
9 //
10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
11 // a collaborative effort of two U.S. Department of Energy organizations (Office
12 // of Science and the National Nuclear Security Administration) responsible for
13 // the planning and preparation of a capable exascale ecosystem, including
14 // software, applications, hardware, advanced system engineering and early
15 // testbed platforms, in support of the nation's exascale computing imperative.
16 
17 /// @file
18 /// Helper functions for solid mechanics example using PETSc
19 
20 #include "../elasticity.h"
21 
22 // -----------------------------------------------------------------------------
23 // Create libCEED operator context
24 // -----------------------------------------------------------------------------
25 // Setup context data for Jacobian evaluation
26 PetscErrorCode SetupJacobianCtx(MPI_Comm comm, AppCtx appCtx, DM dm, Vec V,
27                                 Vec Vloc, CeedData ceedData, Ceed ceed,
28                                 Physics phys, Physics physSmoother,
29                                 UserMult jacobianCtx) {
30   PetscErrorCode ierr;
31 
32   PetscFunctionBeginUser;
33 
34   // PETSc objects
35   jacobianCtx->comm = comm;
36   jacobianCtx->dm = dm;
37 
38   // Work vectors
39   jacobianCtx->Xloc = Vloc;
40   ierr = VecDuplicate(Vloc, &jacobianCtx->Yloc); CHKERRQ(ierr);
41   jacobianCtx->Xceed = ceedData->xceed;
42   jacobianCtx->Yceed = ceedData->yceed;
43 
44   // libCEED operator
45   jacobianCtx->op = ceedData->opJacob;
46   jacobianCtx->qf = ceedData->qfJacob;
47 
48   // Ceed
49   jacobianCtx->ceed = ceed;
50 
51   // Physics
52   jacobianCtx->phys = phys;
53   jacobianCtx->physSmoother = physSmoother;
54 
55   // Get/Restore Array
56   jacobianCtx->memType = appCtx->memTypeRequested;
57   if (appCtx->memTypeRequested == CEED_MEM_HOST) {
58     jacobianCtx->VecGetArray = VecGetArray;
59     jacobianCtx->VecGetArrayRead = VecGetArrayRead;
60     jacobianCtx->VecRestoreArray = VecRestoreArray;
61     jacobianCtx->VecRestoreArrayRead = VecRestoreArrayRead;
62   } else {
63     jacobianCtx->VecGetArray = VecCUDAGetArray;
64     jacobianCtx->VecGetArrayRead = VecCUDAGetArrayRead;
65     jacobianCtx->VecRestoreArray = VecCUDARestoreArray;
66     jacobianCtx->VecRestoreArrayRead = VecCUDARestoreArrayRead;
67   }
68 
69   PetscFunctionReturn(0);
70 };
71 
72 // Setup context data for prolongation and restriction operators
73 PetscErrorCode SetupProlongRestrictCtx(MPI_Comm comm, AppCtx appCtx, DM dmC,
74                                        DM dmF, Vec VF, Vec VlocC, Vec VlocF,
75                                        CeedData ceedDataC, CeedData ceedDataF,
76                                        Ceed ceed,
77                                        UserMultProlongRestr prolongRestrCtx) {
78   PetscErrorCode ierr;
79   PetscScalar *multArray;
80 
81   PetscFunctionBeginUser;
82 
83   // PETSc objects
84   prolongRestrCtx->comm = comm;
85   prolongRestrCtx->dmC = dmC;
86   prolongRestrCtx->dmF = dmF;
87 
88   // Work vectors
89   prolongRestrCtx->locVecC = VlocC;
90   prolongRestrCtx->locVecF = VlocF;
91   prolongRestrCtx->ceedVecC = ceedDataC->xceed;
92   prolongRestrCtx->ceedVecF = ceedDataF->xceed;
93 
94   // libCEED operators
95   prolongRestrCtx->opProlong = ceedDataF->opProlong;
96   prolongRestrCtx->opRestrict = ceedDataF->opRestrict;
97 
98   // Ceed
99   prolongRestrCtx->ceed = ceed;
100 
101   // Get/Restore Array
102   prolongRestrCtx->memType = appCtx->memTypeRequested;
103   if (appCtx->memTypeRequested == CEED_MEM_HOST) {
104     prolongRestrCtx->VecGetArray = VecGetArray;
105     prolongRestrCtx->VecGetArrayRead = VecGetArrayRead;
106     prolongRestrCtx->VecRestoreArray = VecRestoreArray;
107     prolongRestrCtx->VecRestoreArrayRead = VecRestoreArrayRead;
108   } else {
109     prolongRestrCtx->VecGetArray = VecCUDAGetArray;
110     prolongRestrCtx->VecGetArrayRead = VecCUDAGetArrayRead;
111     prolongRestrCtx->VecRestoreArray = VecCUDARestoreArray;
112     prolongRestrCtx->VecRestoreArrayRead = VecCUDARestoreArrayRead;
113   }
114 
115   // Multiplicity vector
116   // -- Set libCEED vector
117   ierr = VecZeroEntries(VlocF);
118   ierr = prolongRestrCtx->VecGetArray(VlocF, &multArray); CHKERRQ(ierr);
119   CeedVectorSetArray(ceedDataF->xceed, appCtx->memTypeRequested,
120                      CEED_USE_POINTER, multArray);
121 
122   // -- Compute multiplicity
123   CeedElemRestrictionGetMultiplicity(ceedDataF->Erestrictu, ceedDataF->xceed);
124 
125   // -- Restore PETSc vector
126   ierr = prolongRestrCtx->VecRestoreArray(VlocF, &multArray); CHKERRQ(ierr);
127 
128   // -- Local-to-global
129   ierr = VecZeroEntries(VF); CHKERRQ(ierr);
130   ierr = DMLocalToGlobal(dmF, VlocF, ADD_VALUES, VF); CHKERRQ(ierr);
131 
132   // -- Global-to-local
133   ierr = VecDuplicate(VlocF, &prolongRestrCtx->multVec); CHKERRQ(ierr);
134   ierr = DMGlobalToLocal(dmF, VF, INSERT_VALUES, prolongRestrCtx->multVec);
135   CHKERRQ(ierr);
136 
137   // -- Reciprocal
138   ierr = VecReciprocal(prolongRestrCtx->multVec); CHKERRQ(ierr);
139 
140   // -- Reset work arrays
141   ierr = VecZeroEntries(VF); CHKERRQ(ierr);
142   ierr = VecZeroEntries(VlocF); CHKERRQ(ierr);
143 
144   PetscFunctionReturn(0);
145 };
146 
147 // -----------------------------------------------------------------------------
148 // Jacobian setup
149 // -----------------------------------------------------------------------------
150 PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat Jpre, void *ctx) {
151   PetscErrorCode ierr;
152 
153   PetscFunctionBeginUser;
154 
155   // Context data
156   FormJacobCtx  formJacobCtx = (FormJacobCtx)ctx;
157   PetscInt      numLevels = formJacobCtx->numLevels;
158   Mat           *jacobMat = formJacobCtx->jacobMat;
159 
160   // Update Jacobian on each level
161   for (PetscInt level = 0; level < numLevels; level++) {
162     ierr = MatAssemblyBegin(jacobMat[level], MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
163     ierr = MatAssemblyEnd(jacobMat[level], MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
164   }
165 
166   // Form coarse assembled matrix
167   ierr = VecZeroEntries(formJacobCtx->Ucoarse); CHKERRQ(ierr);
168   ierr = SNESComputeJacobianDefaultColor(formJacobCtx->snesCoarse,
169                                          formJacobCtx->Ucoarse,
170                                          formJacobCtx->jacobMat[0],
171                                          formJacobCtx->jacobMatCoarse, NULL);
172   CHKERRQ(ierr);
173 
174   // Jpre might be AIJ (e.g., when using coloring), so we need to assemble it
175   ierr = MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
176   ierr = MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
177   if (J != Jpre) {
178     ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
179     ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
180   }
181   PetscFunctionReturn(0);
182 };
183 
184 // -----------------------------------------------------------------------------
185 // Output solution for visualization
186 // -----------------------------------------------------------------------------
187 PetscErrorCode ViewSolution(MPI_Comm comm, Vec U, PetscInt increment,
188                             PetscScalar loadIncrement) {
189   PetscErrorCode ierr;
190   DM dm;
191   PetscViewer viewer;
192   char outputFilename[PETSC_MAX_PATH_LEN];
193 
194   PetscFunctionBeginUser;
195 
196   // Build file name
197   ierr = PetscSNPrintf(outputFilename, sizeof outputFilename,
198                        "solution-%03D.vtu", increment); CHKERRQ(ierr);
199 
200   // Increment sequence
201   ierr = VecGetDM(U, &dm); CHKERRQ(ierr);
202   ierr = DMSetOutputSequenceNumber(dm, increment, loadIncrement); CHKERRQ(ierr);
203 
204   // Output solution vector
205   ierr = PetscViewerVTKOpen(comm, outputFilename, FILE_MODE_WRITE, &viewer);
206   CHKERRQ(ierr);
207   ierr = VecView(U, viewer); CHKERRQ(ierr);
208   ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
209 
210   PetscFunctionReturn(0);
211 };
212 
213 // -----------------------------------------------------------------------------
214 // Output diagnostic quantities for visualization
215 // -----------------------------------------------------------------------------
216 PetscErrorCode ViewDiagnosticQuantities(MPI_Comm comm, DM dmU,
217                                         UserMult user, Vec U,
218                                         CeedElemRestriction ErestrictDiagnostic) {
219   PetscErrorCode ierr;
220   Vec Diagnostic, Yloc, MultVec;
221   CeedVector Yceed;
222   CeedScalar *x, *y;
223   PetscInt lsz;
224   PetscViewer viewer;
225   const char *outputFilename = "diagnostic_quantities.vtu";
226 
227   PetscFunctionBeginUser;
228 
229   // ---------------------------------------------------------------------------
230   // PETSc and libCEED vectors
231   // ---------------------------------------------------------------------------
232   ierr = DMCreateGlobalVector(user->dm, &Diagnostic); CHKERRQ(ierr);
233   ierr = PetscObjectSetName((PetscObject)Diagnostic, ""); CHKERRQ(ierr);
234   ierr = DMCreateLocalVector(user->dm, &Yloc); CHKERRQ(ierr);
235   ierr = VecGetSize(Yloc, &lsz); CHKERRQ(ierr);
236   CeedVectorCreate(user->ceed, lsz, &Yceed);
237 
238   // ---------------------------------------------------------------------------
239   // Compute quantities
240   // ---------------------------------------------------------------------------
241   // -- Global-to-local
242   ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr);
243   ierr = DMPlexInsertBoundaryValues(dmU, PETSC_TRUE, user->Xloc,
244                                     user->loadIncrement, NULL, NULL, NULL);
245   CHKERRQ(ierr);
246   ierr = DMGlobalToLocal(dmU, U, INSERT_VALUES, user->Xloc); CHKERRQ(ierr);
247   ierr = VecZeroEntries(Yloc); CHKERRQ(ierr);
248 
249   // -- Setup CEED vectors
250   ierr = user->VecGetArrayRead(user->Xloc, (const PetscScalar **)&x);
251   CHKERRQ(ierr);
252   ierr = user->VecGetArray(Yloc, &y); CHKERRQ(ierr);
253   CeedVectorSetArray(user->Xceed, user->memType, CEED_USE_POINTER, x);
254   CeedVectorSetArray(Yceed, user->memType, CEED_USE_POINTER, y);
255 
256   // -- Apply CEED operator
257   CeedOperatorApply(user->op, user->Xceed, Yceed, CEED_REQUEST_IMMEDIATE);
258   CeedVectorSyncArray(Yceed, user->memType);
259 
260   // -- Restore PETSc vectors
261   ierr = user->VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x);
262   CHKERRQ(ierr);
263 
264   // -- Local-to-global
265   ierr = VecZeroEntries(Diagnostic); CHKERRQ(ierr);
266   ierr = DMLocalToGlobal(user->dm, Yloc, ADD_VALUES, Diagnostic);
267   CHKERRQ(ierr);
268 
269   // ---------------------------------------------------------------------------
270   // Scale for multiplicity
271   // ---------------------------------------------------------------------------
272   // -- Setup vectors
273   ierr = VecDuplicate(Diagnostic, &MultVec); CHKERRQ(ierr);
274   ierr = VecZeroEntries(Yloc); CHKERRQ(ierr);
275 
276   // -- Compute multiplicity
277   CeedElemRestrictionGetMultiplicity(ErestrictDiagnostic, Yceed);
278 
279   // -- Restore vectors
280   ierr = user->VecRestoreArray(Yloc, &y); CHKERRQ(ierr);
281 
282   // -- Local-to-global
283   ierr = VecZeroEntries(MultVec); CHKERRQ(ierr);
284   ierr = DMLocalToGlobal(user->dm, Yloc, ADD_VALUES, MultVec);
285   CHKERRQ(ierr);
286 
287   // -- Scale
288   ierr = VecReciprocal(MultVec); CHKERRQ(ierr);
289   ierr = VecPointwiseMult(Diagnostic, Diagnostic, MultVec);
290 
291 
292   // ---------------------------------------------------------------------------
293   // Output solution vector
294   // ---------------------------------------------------------------------------
295   ierr = PetscViewerVTKOpen(comm, outputFilename, FILE_MODE_WRITE, &viewer);
296   CHKERRQ(ierr);
297   ierr = VecView(Diagnostic, viewer); CHKERRQ(ierr);
298   ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
299 
300   // ---------------------------------------------------------------------------
301   // Cleanup
302   // ---------------------------------------------------------------------------
303   ierr = VecDestroy(&Diagnostic); CHKERRQ(ierr);
304   ierr = VecDestroy(&MultVec); CHKERRQ(ierr);
305   ierr = VecDestroy(&Yloc); CHKERRQ(ierr);
306   CeedVectorDestroy(&Yceed);
307 
308   PetscFunctionReturn(0);
309 };
310