xref: /libCEED/examples/solids/src/misc.c (revision ce18bed930e8f3bfebcf709a18844aba97fe4630) 
1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED:  http://github.com/ceed
7 
8 /// @file
9 /// Helper functions for solid mechanics example using PETSc
10 
11 #include "../include/misc.h"
12 #include "../include/utils.h"
13 
14 // -----------------------------------------------------------------------------
15 // Create libCEED operator context
16 // -----------------------------------------------------------------------------
17 // Setup context data for Jacobian evaluation
18 PetscErrorCode SetupJacobianCtx(MPI_Comm comm, AppCtx app_ctx, DM dm, Vec V,
19                                 Vec V_loc, CeedData ceed_data, Ceed ceed,
20                                 CeedQFunctionContext ctx_phys,
21                                 CeedQFunctionContext ctx_phys_smoother,
22                                 UserMult jacobian_ctx) {
23   PetscErrorCode ierr;
24 
25   PetscFunctionBeginUser;
26 
27   // PETSc objects
28   jacobian_ctx->comm = comm;
29   jacobian_ctx->dm = dm;
30 
31   // Work vectors
32   jacobian_ctx->X_loc = V_loc;
33   ierr = VecDuplicate(V_loc, &jacobian_ctx->Y_loc); CHKERRQ(ierr);
34   jacobian_ctx->x_ceed = ceed_data->x_ceed;
35   jacobian_ctx->y_ceed = ceed_data->y_ceed;
36 
37   // libCEED operator
38   jacobian_ctx->op = ceed_data->op_jacobian;
39   jacobian_ctx->qf = ceed_data->qf_jacobian;
40 
41   // Ceed
42   jacobian_ctx->ceed = ceed;
43 
44   // Physics
45   jacobian_ctx->ctx_phys = ctx_phys;
46   jacobian_ctx->ctx_phys_smoother = ctx_phys_smoother;
47   PetscFunctionReturn(0);
48 };
49 
50 // Setup context data for prolongation and restriction operators
51 PetscErrorCode SetupProlongRestrictCtx(MPI_Comm comm, AppCtx app_ctx, DM dm_c,
52                                        DM dm_f, Vec V_f, Vec V_loc_c, Vec V_loc_f,
53                                        CeedData ceed_data_c, CeedData ceed_data_f,
54                                        Ceed ceed, UserMultProlongRestr prolong_restr_ctx) {
55   PetscFunctionBeginUser;
56 
57   // PETSc objects
58   prolong_restr_ctx->comm = comm;
59   prolong_restr_ctx->dm_c = dm_c;
60   prolong_restr_ctx->dm_f = dm_f;
61 
62   // Work vectors
63   prolong_restr_ctx->loc_vec_c = V_loc_c;
64   prolong_restr_ctx->loc_vec_f = V_loc_f;
65   prolong_restr_ctx->ceed_vec_c = ceed_data_c->x_ceed;
66   prolong_restr_ctx->ceed_vec_f = ceed_data_f->x_ceed;
67 
68   // libCEED operators
69   prolong_restr_ctx->op_prolong = ceed_data_f->op_prolong;
70   prolong_restr_ctx->op_restrict = ceed_data_f->op_restrict;
71 
72   // Ceed
73   prolong_restr_ctx->ceed = ceed;
74   PetscFunctionReturn(0);
75 };
76 
77 // -----------------------------------------------------------------------------
78 // Jacobian setup
79 // -----------------------------------------------------------------------------
80 PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat J_pre, void *ctx) {
81   PetscErrorCode ierr;
82 
83   PetscFunctionBeginUser;
84 
85   // Context data
86   FormJacobCtx  form_jacob_ctx = (FormJacobCtx)ctx;
87   PetscInt      num_levels = form_jacob_ctx->num_levels;
88   Mat           *jacob_mat = form_jacob_ctx->jacob_mat;
89 
90   // Update Jacobian on each level
91   for (PetscInt level = 0; level < num_levels; level++) {
92     ierr = MatAssemblyBegin(jacob_mat[level], MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
93     ierr = MatAssemblyEnd(jacob_mat[level], MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
94   }
95 
96   // Form coarse assembled matrix
97   CeedOperatorLinearAssemble(form_jacob_ctx->op_coarse,
98                              form_jacob_ctx->coo_values);
99   const CeedScalar *values;
100   CeedVectorGetArrayRead(form_jacob_ctx->coo_values, CEED_MEM_HOST, &values);
101   ierr = MatSetValuesCOO(form_jacob_ctx->jacob_mat_coarse, values, ADD_VALUES);
102   CHKERRQ(ierr);
103   CeedVectorRestoreArrayRead(form_jacob_ctx->coo_values, &values);
104 
105   // J_pre might be AIJ (e.g., when using coloring), so we need to assemble it
106   ierr = MatAssemblyBegin(J_pre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
107   ierr = MatAssemblyEnd(J_pre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
108   if (J != J_pre) {
109     ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
110     ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
111   }
112   PetscFunctionReturn(0);
113 };
114 
115 // -----------------------------------------------------------------------------
116 // Output solution for visualization
117 // -----------------------------------------------------------------------------
118 PetscErrorCode ViewSolution(MPI_Comm comm, AppCtx app_ctx, Vec U,
119                             PetscInt increment, PetscScalar load_increment) {
120   PetscErrorCode ierr;
121   DM dm;
122   PetscViewer viewer;
123   char output_filename[PETSC_MAX_PATH_LEN];
124   PetscMPIInt rank;
125 
126   PetscFunctionBeginUser;
127 
128   // Create output directory
129   MPI_Comm_rank(comm, &rank);
130   if (!rank) {ierr = PetscMkdir(app_ctx->output_dir); CHKERRQ(ierr);}
131 
132   // Build file name
133   ierr = PetscSNPrintf(output_filename, sizeof output_filename,
134                        "%s/solution-%03D.vtu", app_ctx->output_dir,
135                        increment); CHKERRQ(ierr);
136 
137   // Increment sequence
138   ierr = VecGetDM(U, &dm); CHKERRQ(ierr);
139   ierr = DMSetOutputSequenceNumber(dm, increment, load_increment); CHKERRQ(ierr);
140 
141   // Output solution vector
142   ierr = PetscViewerVTKOpen(comm, output_filename, FILE_MODE_WRITE, &viewer);
143   CHKERRQ(ierr);
144   ierr = VecView(U, viewer); CHKERRQ(ierr);
145   ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
146 
147   PetscFunctionReturn(0);
148 };
149 
150 // -----------------------------------------------------------------------------
151 // Output diagnostic quantities for visualization
152 // -----------------------------------------------------------------------------
153 PetscErrorCode ViewDiagnosticQuantities(MPI_Comm comm, DM dmU,
154                                         UserMult user, AppCtx app_ctx, Vec U,
155                                         CeedElemRestriction elem_restr_diagnostic) {
156   PetscErrorCode ierr;
157   Vec Diagnostic, Y_loc, mult_vec;
158   CeedVector y_ceed;
159   CeedScalar *x, *y;
160   PetscMemType x_mem_type, y_mem_type;
161   PetscInt loc_size;
162   PetscViewer viewer;
163   char output_filename[PETSC_MAX_PATH_LEN];
164 
165   PetscFunctionBeginUser;
166 
167   // ---------------------------------------------------------------------------
168   // PETSc and libCEED vectors
169   // ---------------------------------------------------------------------------
170   ierr = DMCreateGlobalVector(user->dm, &Diagnostic); CHKERRQ(ierr);
171   ierr = PetscObjectSetName((PetscObject)Diagnostic, ""); CHKERRQ(ierr);
172   ierr = DMCreateLocalVector(user->dm, &Y_loc); CHKERRQ(ierr);
173   ierr = VecGetSize(Y_loc, &loc_size); CHKERRQ(ierr);
174   CeedVectorCreate(user->ceed, loc_size, &y_ceed);
175 
176   // ---------------------------------------------------------------------------
177   // Compute quantities
178   // ---------------------------------------------------------------------------
179   // -- Global-to-local
180   ierr = VecZeroEntries(user->X_loc); CHKERRQ(ierr);
181   ierr = DMPlexInsertBoundaryValues(dmU, PETSC_TRUE, user->X_loc,
182                                     user->load_increment, NULL, NULL, NULL);
183   CHKERRQ(ierr);
184   ierr = DMGlobalToLocal(dmU, U, INSERT_VALUES, user->X_loc); CHKERRQ(ierr);
185   ierr = VecZeroEntries(Y_loc); CHKERRQ(ierr);
186 
187   // -- Setup CEED vectors
188   ierr = VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x,
189                                    &x_mem_type); CHKERRQ(ierr);
190   ierr = VecGetArrayAndMemType(Y_loc, &y, &y_mem_type); CHKERRQ(ierr);
191   CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x);
192   CeedVectorSetArray(y_ceed, MemTypeP2C(y_mem_type), CEED_USE_POINTER, y);
193 
194   // -- Apply CEED operator
195   CeedOperatorApply(user->op, user->x_ceed, y_ceed, CEED_REQUEST_IMMEDIATE);
196 
197   // -- Restore PETSc vector; keep y_ceed viewing memory of Y_loc for use below
198   CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL);
199   ierr = VecRestoreArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x);
200   CHKERRQ(ierr);
201 
202   // -- Local-to-global
203   ierr = VecZeroEntries(Diagnostic); CHKERRQ(ierr);
204   ierr = DMLocalToGlobal(user->dm, Y_loc, ADD_VALUES, Diagnostic);
205   CHKERRQ(ierr);
206 
207   // ---------------------------------------------------------------------------
208   // Scale for multiplicity
209   // ---------------------------------------------------------------------------
210   // -- Setup vectors
211   ierr = VecDuplicate(Diagnostic, &mult_vec); CHKERRQ(ierr);
212   ierr = VecZeroEntries(Y_loc); CHKERRQ(ierr);
213 
214   // -- Compute multiplicity
215   CeedElemRestrictionGetMultiplicity(elem_restr_diagnostic, y_ceed);
216 
217   // -- Restore vectors
218   CeedVectorTakeArray(y_ceed, MemTypeP2C(y_mem_type), NULL);
219   ierr = VecRestoreArrayAndMemType(Y_loc, &y); CHKERRQ(ierr);
220 
221   // -- Local-to-global
222   ierr = VecZeroEntries(mult_vec); CHKERRQ(ierr);
223   ierr = DMLocalToGlobal(user->dm, Y_loc, ADD_VALUES, mult_vec);
224   CHKERRQ(ierr);
225 
226   // -- Scale
227   ierr = VecReciprocal(mult_vec); CHKERRQ(ierr);
228   ierr = VecPointwiseMult(Diagnostic, Diagnostic, mult_vec);
229 
230   // ---------------------------------------------------------------------------
231   // Output solution vector
232   // ---------------------------------------------------------------------------
233   ierr = PetscSNPrintf(output_filename, sizeof output_filename,
234                        "%s/diagnostic_quantities.vtu",
235                        app_ctx->output_dir); CHKERRQ(ierr);
236   ierr = PetscViewerVTKOpen(comm, output_filename, FILE_MODE_WRITE, &viewer);
237   CHKERRQ(ierr);
238   ierr = VecView(Diagnostic, viewer); CHKERRQ(ierr);
239   ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
240 
241   // ---------------------------------------------------------------------------
242   // Cleanup
243   // ---------------------------------------------------------------------------
244   ierr = VecDestroy(&Diagnostic); CHKERRQ(ierr);
245   ierr = VecDestroy(&mult_vec); CHKERRQ(ierr);
246   ierr = VecDestroy(&Y_loc); CHKERRQ(ierr);
247   CeedVectorDestroy(&y_ceed);
248 
249   PetscFunctionReturn(0);
250 };
251 
252 // -----------------------------------------------------------------------------
253 // Regression testing
254 // -----------------------------------------------------------------------------
255 // test option change. could remove the loading step. Run only with one loading step and compare relatively to ref file
256 // option: expect_final_strain_energy and check against the relative error to ref is within tolerance (10^-5) I.e. one Newton solve then check final energy
257 PetscErrorCode RegressionTests_solids(AppCtx app_ctx, PetscReal energy) {
258   PetscFunctionBegin;
259 
260   if (app_ctx->expect_final_strain >= 0.) {
261     PetscReal energy_ref = app_ctx->expect_final_strain;
262     PetscReal error = PetscAbsReal(energy - energy_ref) / energy_ref;
263 
264     if (error > app_ctx->test_tol) {
265       PetscErrorCode ierr;
266       ierr = PetscPrintf(PETSC_COMM_WORLD,
267                          "Energy %e does not match expected energy %e: relative tolerance %e > %e\n",
268                          (double)energy, (double)energy_ref, (double)error, app_ctx->test_tol);
269       CHKERRQ(ierr);
270     }
271   }
272   PetscFunctionReturn(0);
273 };
274