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