xref: /libCEED/examples/solids/src/misc.c (revision eccc2849f69a9f016cade2e1c45046d05a5ce45c) 
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, AppCtx appCtx, Vec U,
127                             PetscInt increment,
128                             PetscScalar loadIncrement) {
129   PetscErrorCode ierr;
130   DM dm;
131   PetscViewer viewer;
132   char outputFilename[PETSC_MAX_PATH_LEN];
133   PetscMPIInt rank;
134 
135   PetscFunctionBeginUser;
136 
137   // Create output directory
138   MPI_Comm_rank(comm, &rank);
139   if (!rank) {ierr = PetscMkdir(appCtx->outputdir); CHKERRQ(ierr);}
140 
141   // Build file name
142   ierr = PetscSNPrintf(outputFilename, sizeof outputFilename,
143                        "%s/solution-%03D.vtu", appCtx->outputdir,
144                        increment); CHKERRQ(ierr);
145 
146   // Increment sequence
147   ierr = VecGetDM(U, &dm); CHKERRQ(ierr);
148   ierr = DMSetOutputSequenceNumber(dm, increment, loadIncrement); CHKERRQ(ierr);
149 
150   // Output solution vector
151   ierr = PetscViewerVTKOpen(comm, outputFilename, FILE_MODE_WRITE, &viewer);
152   CHKERRQ(ierr);
153   ierr = VecView(U, viewer); CHKERRQ(ierr);
154   ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
155 
156   PetscFunctionReturn(0);
157 };
158 
159 // -----------------------------------------------------------------------------
160 // Output diagnostic quantities for visualization
161 // -----------------------------------------------------------------------------
162 PetscErrorCode ViewDiagnosticQuantities(MPI_Comm comm, DM dmU,
163                                         UserMult user, AppCtx appCtx, Vec U,
164                                         CeedElemRestriction ErestrictDiagnostic) {
165   PetscErrorCode ierr;
166   Vec Diagnostic, Yloc, MultVec;
167   CeedVector Yceed;
168   CeedScalar *x, *y;
169   PetscMemType xmemtype, ymemtype;
170   PetscInt lsz;
171   PetscViewer viewer;
172   char outputFilename[PETSC_MAX_PATH_LEN];
173 
174   PetscFunctionBeginUser;
175 
176   // ---------------------------------------------------------------------------
177   // PETSc and libCEED vectors
178   // ---------------------------------------------------------------------------
179   ierr = DMCreateGlobalVector(user->dm, &Diagnostic); CHKERRQ(ierr);
180   ierr = PetscObjectSetName((PetscObject)Diagnostic, ""); CHKERRQ(ierr);
181   ierr = DMCreateLocalVector(user->dm, &Yloc); CHKERRQ(ierr);
182   ierr = VecGetSize(Yloc, &lsz); CHKERRQ(ierr);
183   CeedVectorCreate(user->ceed, lsz, &Yceed);
184 
185   // ---------------------------------------------------------------------------
186   // Compute quantities
187   // ---------------------------------------------------------------------------
188   // -- Global-to-local
189   ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr);
190   ierr = DMPlexInsertBoundaryValues(dmU, PETSC_TRUE, user->Xloc,
191                                     user->loadIncrement, NULL, NULL, NULL);
192   CHKERRQ(ierr);
193   ierr = DMGlobalToLocal(dmU, U, INSERT_VALUES, user->Xloc); CHKERRQ(ierr);
194   ierr = VecZeroEntries(Yloc); CHKERRQ(ierr);
195 
196   // -- Setup CEED vectors
197   ierr = VecGetArrayReadAndMemType(user->Xloc, (const PetscScalar **)&x,
198                                    &xmemtype); CHKERRQ(ierr);
199   ierr = VecGetArrayAndMemType(Yloc, &y, &ymemtype); CHKERRQ(ierr);
200   CeedVectorSetArray(user->Xceed, MemTypeP2C(xmemtype), CEED_USE_POINTER, x);
201   CeedVectorSetArray(Yceed, MemTypeP2C(ymemtype), CEED_USE_POINTER, y);
202 
203   // -- Apply CEED operator
204   CeedOperatorApply(user->op, user->Xceed, Yceed, CEED_REQUEST_IMMEDIATE);
205 
206   // -- Restore PETSc vector; keep Yceed viewing memory of Yloc for use below
207   CeedVectorTakeArray(user->Xceed, MemTypeP2C(xmemtype), NULL);
208   ierr = VecRestoreArrayReadAndMemType(user->Xloc, (const PetscScalar **)&x);
209   CHKERRQ(ierr);
210 
211   // -- Local-to-global
212   ierr = VecZeroEntries(Diagnostic); CHKERRQ(ierr);
213   ierr = DMLocalToGlobal(user->dm, Yloc, ADD_VALUES, Diagnostic);
214   CHKERRQ(ierr);
215 
216   // ---------------------------------------------------------------------------
217   // Scale for multiplicity
218   // ---------------------------------------------------------------------------
219   // -- Setup vectors
220   ierr = VecDuplicate(Diagnostic, &MultVec); CHKERRQ(ierr);
221   ierr = VecZeroEntries(Yloc); CHKERRQ(ierr);
222 
223   // -- Compute multiplicity
224   CeedElemRestrictionGetMultiplicity(ErestrictDiagnostic, Yceed);
225 
226   // -- Restore vectors
227   CeedVectorTakeArray(Yceed, MemTypeP2C(ymemtype), NULL);
228   ierr = VecRestoreArrayAndMemType(Yloc, &y); CHKERRQ(ierr);
229 
230   // -- Local-to-global
231   ierr = VecZeroEntries(MultVec); CHKERRQ(ierr);
232   ierr = DMLocalToGlobal(user->dm, Yloc, ADD_VALUES, MultVec);
233   CHKERRQ(ierr);
234 
235   // -- Scale
236   ierr = VecReciprocal(MultVec); CHKERRQ(ierr);
237   ierr = VecPointwiseMult(Diagnostic, Diagnostic, MultVec);
238 
239   // ---------------------------------------------------------------------------
240   // Output solution vector
241   // ---------------------------------------------------------------------------
242   ierr = PetscSNPrintf(outputFilename, sizeof outputFilename,
243                        "%s/diagnostic_quantities.vtu",
244                        appCtx->outputdir); CHKERRQ(ierr);
245   ierr = PetscViewerVTKOpen(comm, outputFilename, FILE_MODE_WRITE, &viewer);
246   CHKERRQ(ierr);
247   ierr = VecView(Diagnostic, viewer); CHKERRQ(ierr);
248   ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
249 
250   // ---------------------------------------------------------------------------
251   // Cleanup
252   // ---------------------------------------------------------------------------
253   ierr = VecDestroy(&Diagnostic); CHKERRQ(ierr);
254   ierr = VecDestroy(&MultVec); CHKERRQ(ierr);
255   ierr = VecDestroy(&Yloc); CHKERRQ(ierr);
256   CeedVectorDestroy(&Yceed);
257 
258   PetscFunctionReturn(0);
259 };
260