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