xref: /libCEED/examples/solids/src/misc.c (revision 9d488d037550e1c84877377029f08e264b1941e6) 
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   CeedVectorTakeArray(ceedDataF->xceed, appCtx->memTypeRequested, NULL);
127   ierr = prolongRestrCtx->VecRestoreArray(VlocF, &multArray); CHKERRQ(ierr);
128 
129   // -- Local-to-global
130   ierr = VecZeroEntries(VF); CHKERRQ(ierr);
131   ierr = DMLocalToGlobal(dmF, VlocF, ADD_VALUES, VF); CHKERRQ(ierr);
132 
133   // -- Global-to-local
134   ierr = VecDuplicate(VlocF, &prolongRestrCtx->multVec); CHKERRQ(ierr);
135   ierr = DMGlobalToLocal(dmF, VF, INSERT_VALUES, prolongRestrCtx->multVec);
136   CHKERRQ(ierr);
137 
138   // -- Reciprocal
139   ierr = VecReciprocal(prolongRestrCtx->multVec); CHKERRQ(ierr);
140 
141   // -- Reset work arrays
142   ierr = VecZeroEntries(VF); CHKERRQ(ierr);
143   ierr = VecZeroEntries(VlocF); CHKERRQ(ierr);
144 
145   PetscFunctionReturn(0);
146 };
147 
148 // -----------------------------------------------------------------------------
149 // Jacobian setup
150 // -----------------------------------------------------------------------------
151 PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat Jpre, void *ctx) {
152   PetscErrorCode ierr;
153 
154   PetscFunctionBeginUser;
155 
156   // Context data
157   FormJacobCtx  formJacobCtx = (FormJacobCtx)ctx;
158   PetscInt      numLevels = formJacobCtx->numLevels;
159   Mat           *jacobMat = formJacobCtx->jacobMat;
160 
161   // Update Jacobian on each level
162   for (PetscInt level = 0; level < numLevels; level++) {
163     ierr = MatAssemblyBegin(jacobMat[level], MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
164     ierr = MatAssemblyEnd(jacobMat[level], MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
165   }
166 
167   // Form coarse assembled matrix
168   ierr = VecZeroEntries(formJacobCtx->Ucoarse); CHKERRQ(ierr);
169   ierr = SNESComputeJacobianDefaultColor(formJacobCtx->snesCoarse,
170                                          formJacobCtx->Ucoarse,
171                                          formJacobCtx->jacobMat[0],
172                                          formJacobCtx->jacobMatCoarse, NULL);
173   CHKERRQ(ierr);
174 
175   // Jpre might be AIJ (e.g., when using coloring), so we need to assemble it
176   ierr = MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
177   ierr = MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
178   if (J != Jpre) {
179     ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
180     ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
181   }
182   PetscFunctionReturn(0);
183 };
184 
185 // -----------------------------------------------------------------------------
186 // Output solution for visualization
187 // -----------------------------------------------------------------------------
188 PetscErrorCode ViewSolution(MPI_Comm comm, Vec U, PetscInt increment,
189                             PetscScalar loadIncrement) {
190   PetscErrorCode ierr;
191   DM dm;
192   PetscViewer viewer;
193   char outputFilename[PETSC_MAX_PATH_LEN];
194 
195   PetscFunctionBeginUser;
196 
197   // Build file name
198   ierr = PetscSNPrintf(outputFilename, sizeof outputFilename,
199                        "solution-%03D.vtu", increment); CHKERRQ(ierr);
200 
201   // Increment sequence
202   ierr = VecGetDM(U, &dm); CHKERRQ(ierr);
203   ierr = DMSetOutputSequenceNumber(dm, increment, loadIncrement); CHKERRQ(ierr);
204 
205   // Output solution vector
206   ierr = PetscViewerVTKOpen(comm, outputFilename, FILE_MODE_WRITE, &viewer);
207   CHKERRQ(ierr);
208   ierr = VecView(U, viewer); CHKERRQ(ierr);
209   ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
210 
211   PetscFunctionReturn(0);
212 };
213 
214 // -----------------------------------------------------------------------------
215 // Output diagnostic quantities for visualization
216 // -----------------------------------------------------------------------------
217 PetscErrorCode ViewDiagnosticQuantities(MPI_Comm comm, DM dmU,
218                                         UserMult user, Vec U,
219                                         CeedElemRestriction ErestrictDiagnostic) {
220   PetscErrorCode ierr;
221   Vec Diagnostic, Yloc, MultVec;
222   CeedVector Yceed;
223   CeedScalar *x, *y;
224   PetscInt lsz;
225   PetscViewer viewer;
226   const char *outputFilename = "diagnostic_quantities.vtu";
227 
228   PetscFunctionBeginUser;
229 
230   // ---------------------------------------------------------------------------
231   // PETSc and libCEED vectors
232   // ---------------------------------------------------------------------------
233   ierr = DMCreateGlobalVector(user->dm, &Diagnostic); CHKERRQ(ierr);
234   ierr = PetscObjectSetName((PetscObject)Diagnostic, ""); CHKERRQ(ierr);
235   ierr = DMCreateLocalVector(user->dm, &Yloc); CHKERRQ(ierr);
236   ierr = VecGetSize(Yloc, &lsz); CHKERRQ(ierr);
237   CeedVectorCreate(user->ceed, lsz, &Yceed);
238 
239   // ---------------------------------------------------------------------------
240   // Compute quantities
241   // ---------------------------------------------------------------------------
242   // -- Global-to-local
243   ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr);
244   ierr = DMPlexInsertBoundaryValues(dmU, PETSC_TRUE, user->Xloc,
245                                     user->loadIncrement, NULL, NULL, NULL);
246   CHKERRQ(ierr);
247   ierr = DMGlobalToLocal(dmU, U, INSERT_VALUES, user->Xloc); CHKERRQ(ierr);
248   ierr = VecZeroEntries(Yloc); CHKERRQ(ierr);
249 
250   // -- Setup CEED vectors
251   ierr = user->VecGetArrayRead(user->Xloc, (const PetscScalar **)&x);
252   CHKERRQ(ierr);
253   ierr = user->VecGetArray(Yloc, &y); CHKERRQ(ierr);
254   CeedVectorSetArray(user->Xceed, user->memType, CEED_USE_POINTER, x);
255   CeedVectorSetArray(Yceed, user->memType, CEED_USE_POINTER, y);
256 
257   // -- Apply CEED operator
258   CeedOperatorApply(user->op, user->Xceed, Yceed, CEED_REQUEST_IMMEDIATE);
259 
260   // -- Restore PETSc vector
261   CeedVectorTakeArray(user->Xceed, user->memType, NULL);
262   ierr = user->VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x);
263   CHKERRQ(ierr);
264 
265   // -- Local-to-global
266   ierr = VecZeroEntries(Diagnostic); CHKERRQ(ierr);
267   ierr = DMLocalToGlobal(user->dm, Yloc, ADD_VALUES, Diagnostic);
268   CHKERRQ(ierr);
269 
270   // ---------------------------------------------------------------------------
271   // Scale for multiplicity
272   // ---------------------------------------------------------------------------
273   // -- Setup vectors
274   ierr = VecDuplicate(Diagnostic, &MultVec); CHKERRQ(ierr);
275   ierr = VecZeroEntries(Yloc); CHKERRQ(ierr);
276 
277   // -- Compute multiplicity
278   CeedElemRestrictionGetMultiplicity(ErestrictDiagnostic, Yceed);
279 
280   // -- Restore vectors
281   CeedVectorTakeArray(Yceed, user->memType, NULL);
282   ierr = user->VecRestoreArray(Yloc, &y); CHKERRQ(ierr);
283 
284   // -- Local-to-global
285   ierr = VecZeroEntries(MultVec); CHKERRQ(ierr);
286   ierr = DMLocalToGlobal(user->dm, Yloc, ADD_VALUES, MultVec);
287   CHKERRQ(ierr);
288 
289   // -- Scale
290   ierr = VecReciprocal(MultVec); CHKERRQ(ierr);
291   ierr = VecPointwiseMult(Diagnostic, Diagnostic, MultVec);
292 
293 
294   // ---------------------------------------------------------------------------
295   // Output solution vector
296   // ---------------------------------------------------------------------------
297   ierr = PetscViewerVTKOpen(comm, outputFilename, FILE_MODE_WRITE, &viewer);
298   CHKERRQ(ierr);
299   ierr = VecView(Diagnostic, viewer); CHKERRQ(ierr);
300   ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
301 
302   // ---------------------------------------------------------------------------
303   // Cleanup
304   // ---------------------------------------------------------------------------
305   ierr = VecDestroy(&Diagnostic); CHKERRQ(ierr);
306   ierr = VecDestroy(&MultVec); CHKERRQ(ierr);
307   ierr = VecDestroy(&Yloc); CHKERRQ(ierr);
308   CeedVectorDestroy(&Yceed);
309 
310   PetscFunctionReturn(0);
311 };
312