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