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