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