xref: /libCEED/examples/solids/src/misc.c (revision 91dfd1cde504855b06db4052ed3242027f67af19) 
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   PetscFunctionBeginUser;
79 
80   // PETSc objects
81   prolongRestrCtx->comm = comm;
82   prolongRestrCtx->dmC = dmC;
83   prolongRestrCtx->dmF = dmF;
84 
85   // Work vectors
86   prolongRestrCtx->locVecC = VlocC;
87   prolongRestrCtx->locVecF = VlocF;
88   prolongRestrCtx->ceedVecC = ceedDataC->xceed;
89   prolongRestrCtx->ceedVecF = ceedDataF->xceed;
90 
91   // libCEED operators
92   prolongRestrCtx->opProlong = ceedDataF->opProlong;
93   prolongRestrCtx->opRestrict = ceedDataF->opRestrict;
94 
95   // Ceed
96   prolongRestrCtx->ceed = ceed;
97 
98   // Get/Restore Array
99   prolongRestrCtx->memType = appCtx->memTypeRequested;
100   if (appCtx->memTypeRequested == CEED_MEM_HOST) {
101     prolongRestrCtx->VecGetArray = VecGetArray;
102     prolongRestrCtx->VecGetArrayRead = VecGetArrayRead;
103     prolongRestrCtx->VecRestoreArray = VecRestoreArray;
104     prolongRestrCtx->VecRestoreArrayRead = VecRestoreArrayRead;
105   } else {
106     prolongRestrCtx->VecGetArray = VecCUDAGetArray;
107     prolongRestrCtx->VecGetArrayRead = VecCUDAGetArrayRead;
108     prolongRestrCtx->VecRestoreArray = VecCUDARestoreArray;
109     prolongRestrCtx->VecRestoreArrayRead = VecCUDARestoreArrayRead;
110   }
111 
112   PetscFunctionReturn(0);
113 };
114 
115 // -----------------------------------------------------------------------------
116 // Jacobian setup
117 // -----------------------------------------------------------------------------
118 PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat Jpre, void *ctx) {
119   PetscErrorCode ierr;
120 
121   PetscFunctionBeginUser;
122 
123   // Context data
124   FormJacobCtx  formJacobCtx = (FormJacobCtx)ctx;
125   PetscInt      numLevels = formJacobCtx->numLevels;
126   Mat           *jacobMat = formJacobCtx->jacobMat;
127 
128   // Update Jacobian on each level
129   for (PetscInt level = 0; level < numLevels; level++) {
130     ierr = MatAssemblyBegin(jacobMat[level], MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
131     ierr = MatAssemblyEnd(jacobMat[level], MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
132   }
133 
134   // Form coarse assembled matrix
135   ierr = VecZeroEntries(formJacobCtx->Ucoarse); CHKERRQ(ierr);
136   ierr = SNESComputeJacobianDefaultColor(formJacobCtx->snesCoarse,
137                                          formJacobCtx->Ucoarse,
138                                          formJacobCtx->jacobMat[0],
139                                          formJacobCtx->jacobMatCoarse, NULL);
140   CHKERRQ(ierr);
141 
142   // Jpre might be AIJ (e.g., when using coloring), so we need to assemble it
143   ierr = MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
144   ierr = MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
145   if (J != Jpre) {
146     ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
147     ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
148   }
149   PetscFunctionReturn(0);
150 };
151 
152 // -----------------------------------------------------------------------------
153 // Output solution for visualization
154 // -----------------------------------------------------------------------------
155 PetscErrorCode ViewSolution(MPI_Comm comm, Vec U, PetscInt increment,
156                             PetscScalar loadIncrement) {
157   PetscErrorCode ierr;
158   DM dm;
159   PetscViewer viewer;
160   char outputFilename[PETSC_MAX_PATH_LEN];
161 
162   PetscFunctionBeginUser;
163 
164   // Build file name
165   ierr = PetscSNPrintf(outputFilename, sizeof outputFilename,
166                        "solution-%03D.vtu", increment); CHKERRQ(ierr);
167 
168   // Increment sequence
169   ierr = VecGetDM(U, &dm); CHKERRQ(ierr);
170   ierr = DMSetOutputSequenceNumber(dm, increment, loadIncrement); CHKERRQ(ierr);
171 
172   // Output solution vector
173   ierr = PetscViewerVTKOpen(comm, outputFilename, FILE_MODE_WRITE, &viewer);
174   CHKERRQ(ierr);
175   ierr = VecView(U, viewer); CHKERRQ(ierr);
176   ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
177 
178   PetscFunctionReturn(0);
179 };
180 
181 // -----------------------------------------------------------------------------
182 // Output diagnostic quantities for visualization
183 // -----------------------------------------------------------------------------
184 PetscErrorCode ViewDiagnosticQuantities(MPI_Comm comm, DM dmU,
185                                         UserMult user, Vec U,
186                                         CeedElemRestriction ErestrictDiagnostic) {
187   PetscErrorCode ierr;
188   Vec Diagnostic, Yloc, MultVec;
189   CeedVector Yceed;
190   CeedScalar *x, *y;
191   PetscInt lsz;
192   PetscViewer viewer;
193   const char *outputFilename = "diagnostic_quantities.vtu";
194 
195   PetscFunctionBeginUser;
196 
197   // ---------------------------------------------------------------------------
198   // PETSc and libCEED vectors
199   // ---------------------------------------------------------------------------
200   ierr = DMCreateGlobalVector(user->dm, &Diagnostic); CHKERRQ(ierr);
201   ierr = PetscObjectSetName((PetscObject)Diagnostic, ""); CHKERRQ(ierr);
202   ierr = DMCreateLocalVector(user->dm, &Yloc); CHKERRQ(ierr);
203   ierr = VecGetSize(Yloc, &lsz); CHKERRQ(ierr);
204   CeedVectorCreate(user->ceed, lsz, &Yceed);
205 
206   // ---------------------------------------------------------------------------
207   // Compute quantities
208   // ---------------------------------------------------------------------------
209   // -- Global-to-local
210   ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr);
211   ierr = DMPlexInsertBoundaryValues(dmU, PETSC_TRUE, user->Xloc,
212                                     user->loadIncrement, NULL, NULL, NULL);
213   CHKERRQ(ierr);
214   ierr = DMGlobalToLocal(dmU, U, INSERT_VALUES, user->Xloc); CHKERRQ(ierr);
215   ierr = VecZeroEntries(Yloc); CHKERRQ(ierr);
216 
217   // -- Setup CEED vectors
218   ierr = user->VecGetArrayRead(user->Xloc, (const PetscScalar **)&x);
219   CHKERRQ(ierr);
220   ierr = user->VecGetArray(Yloc, &y); CHKERRQ(ierr);
221   CeedVectorSetArray(user->Xceed, user->memType, CEED_USE_POINTER, x);
222   CeedVectorSetArray(Yceed, user->memType, CEED_USE_POINTER, y);
223 
224   // -- Apply CEED operator
225   CeedOperatorApply(user->op, user->Xceed, Yceed, CEED_REQUEST_IMMEDIATE);
226 
227   // -- Restore PETSc vector
228   CeedVectorTakeArray(user->Xceed, user->memType, NULL);
229   ierr = user->VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x);
230   CHKERRQ(ierr);
231 
232   // -- Local-to-global
233   ierr = VecZeroEntries(Diagnostic); CHKERRQ(ierr);
234   ierr = DMLocalToGlobal(user->dm, Yloc, ADD_VALUES, Diagnostic);
235   CHKERRQ(ierr);
236 
237   // ---------------------------------------------------------------------------
238   // Scale for multiplicity
239   // ---------------------------------------------------------------------------
240   // -- Setup vectors
241   ierr = VecDuplicate(Diagnostic, &MultVec); CHKERRQ(ierr);
242   ierr = VecZeroEntries(Yloc); CHKERRQ(ierr);
243 
244   // -- Compute multiplicity
245   CeedElemRestrictionGetMultiplicity(ErestrictDiagnostic, Yceed);
246 
247   // -- Restore vectors
248   CeedVectorTakeArray(Yceed, user->memType, NULL);
249   ierr = user->VecRestoreArray(Yloc, &y); CHKERRQ(ierr);
250 
251   // -- Local-to-global
252   ierr = VecZeroEntries(MultVec); CHKERRQ(ierr);
253   ierr = DMLocalToGlobal(user->dm, Yloc, ADD_VALUES, MultVec);
254   CHKERRQ(ierr);
255 
256   // -- Scale
257   ierr = VecReciprocal(MultVec); CHKERRQ(ierr);
258   ierr = VecPointwiseMult(Diagnostic, Diagnostic, MultVec);
259 
260   // ---------------------------------------------------------------------------
261   // Output solution vector
262   // ---------------------------------------------------------------------------
263   ierr = PetscViewerVTKOpen(comm, outputFilename, FILE_MODE_WRITE, &viewer);
264   CHKERRQ(ierr);
265   ierr = VecView(Diagnostic, viewer); CHKERRQ(ierr);
266   ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
267 
268   // ---------------------------------------------------------------------------
269   // Cleanup
270   // ---------------------------------------------------------------------------
271   ierr = VecDestroy(&Diagnostic); CHKERRQ(ierr);
272   ierr = VecDestroy(&MultVec); CHKERRQ(ierr);
273   ierr = VecDestroy(&Yloc); CHKERRQ(ierr);
274   CeedVectorDestroy(&Yceed);
275 
276   PetscFunctionReturn(0);
277 };
278