xref: /libCEED/examples/solids/src/misc.c (revision d83b856e16e1a05ebd491d7f65187d3c08a4ff17) 
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   ierr = VecZeroEntries(form_jacob_ctx->u_coarse); CHKERRQ(ierr);
91   ierr = SNESComputeJacobianDefaultColor(form_jacob_ctx->snes_coarse,
92                                          form_jacob_ctx->u_coarse,
93                                          form_jacob_ctx->jacob_mat[0],
94                                          form_jacob_ctx->jacob_mat_coarse, NULL);
95   CHKERRQ(ierr);
96 
97   // J_pre might be AIJ (e.g., when using coloring), so we need to assemble it
98   ierr = MatAssemblyBegin(J_pre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
99   ierr = MatAssemblyEnd(J_pre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
100   if (J != J_pre) {
101     ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
102     ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
103   }
104   PetscFunctionReturn(0);
105 };
106 
107 // -----------------------------------------------------------------------------
108 // Output solution for visualization
109 // -----------------------------------------------------------------------------
110 PetscErrorCode ViewSolution(MPI_Comm comm, AppCtx app_ctx, Vec U,
111                             PetscInt increment, PetscScalar load_increment) {
112   PetscErrorCode ierr;
113   DM dm;
114   PetscViewer viewer;
115   char output_filename[PETSC_MAX_PATH_LEN];
116   PetscMPIInt rank;
117 
118   PetscFunctionBeginUser;
119 
120   // Create output directory
121   MPI_Comm_rank(comm, &rank);
122   if (!rank) {ierr = PetscMkdir(app_ctx->output_dir); CHKERRQ(ierr);}
123 
124   // Build file name
125   ierr = PetscSNPrintf(output_filename, sizeof output_filename,
126                        "%s/solution-%03D.vtu", app_ctx->output_dir,
127                        increment); CHKERRQ(ierr);
128 
129   // Increment sequence
130   ierr = VecGetDM(U, &dm); CHKERRQ(ierr);
131   ierr = DMSetOutputSequenceNumber(dm, increment, load_increment); CHKERRQ(ierr);
132 
133   // Output solution vector
134   ierr = PetscViewerVTKOpen(comm, output_filename, FILE_MODE_WRITE, &viewer);
135   CHKERRQ(ierr);
136   ierr = VecView(U, viewer); CHKERRQ(ierr);
137   ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
138 
139   PetscFunctionReturn(0);
140 };
141 
142 // -----------------------------------------------------------------------------
143 // Output diagnostic quantities for visualization
144 // -----------------------------------------------------------------------------
145 PetscErrorCode ViewDiagnosticQuantities(MPI_Comm comm, DM dmU,
146                                         UserMult user, AppCtx app_ctx, Vec U,
147                                         CeedElemRestriction elem_restr_diagnostic) {
148   PetscErrorCode ierr;
149   Vec Diagnostic, Y_loc, mult_vec;
150   CeedVector y_ceed;
151   CeedScalar *x, *y;
152   PetscMemType x_mem_type, y_mem_type;
153   PetscInt loc_size;
154   PetscViewer viewer;
155   char output_filename[PETSC_MAX_PATH_LEN];
156 
157   PetscFunctionBeginUser;
158 
159   // ---------------------------------------------------------------------------
160   // PETSc and libCEED vectors
161   // ---------------------------------------------------------------------------
162   ierr = DMCreateGlobalVector(user->dm, &Diagnostic); CHKERRQ(ierr);
163   ierr = PetscObjectSetName((PetscObject)Diagnostic, ""); CHKERRQ(ierr);
164   ierr = DMCreateLocalVector(user->dm, &Y_loc); CHKERRQ(ierr);
165   ierr = VecGetSize(Y_loc, &loc_size); CHKERRQ(ierr);
166   CeedVectorCreate(user->ceed, loc_size, &y_ceed);
167 
168   // ---------------------------------------------------------------------------
169   // Compute quantities
170   // ---------------------------------------------------------------------------
171   // -- Global-to-local
172   ierr = VecZeroEntries(user->X_loc); CHKERRQ(ierr);
173   ierr = DMPlexInsertBoundaryValues(dmU, PETSC_TRUE, user->X_loc,
174                                     user->load_increment, NULL, NULL, NULL);
175   CHKERRQ(ierr);
176   ierr = DMGlobalToLocal(dmU, U, INSERT_VALUES, user->X_loc); CHKERRQ(ierr);
177   ierr = VecZeroEntries(Y_loc); CHKERRQ(ierr);
178 
179   // -- Setup CEED vectors
180   ierr = VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x,
181                                    &x_mem_type); CHKERRQ(ierr);
182   ierr = VecGetArrayAndMemType(Y_loc, &y, &y_mem_type); CHKERRQ(ierr);
183   CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x);
184   CeedVectorSetArray(y_ceed, MemTypeP2C(y_mem_type), CEED_USE_POINTER, y);
185 
186   // -- Apply CEED operator
187   CeedOperatorApply(user->op, user->x_ceed, y_ceed, CEED_REQUEST_IMMEDIATE);
188 
189   // -- Restore PETSc vector; keep y_ceed viewing memory of Y_loc for use below
190   CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL);
191   ierr = VecRestoreArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x);
192   CHKERRQ(ierr);
193 
194   // -- Local-to-global
195   ierr = VecZeroEntries(Diagnostic); CHKERRQ(ierr);
196   ierr = DMLocalToGlobal(user->dm, Y_loc, ADD_VALUES, Diagnostic);
197   CHKERRQ(ierr);
198 
199   // ---------------------------------------------------------------------------
200   // Scale for multiplicity
201   // ---------------------------------------------------------------------------
202   // -- Setup vectors
203   ierr = VecDuplicate(Diagnostic, &mult_vec); CHKERRQ(ierr);
204   ierr = VecZeroEntries(Y_loc); CHKERRQ(ierr);
205 
206   // -- Compute multiplicity
207   CeedElemRestrictionGetMultiplicity(elem_restr_diagnostic, y_ceed);
208 
209   // -- Restore vectors
210   CeedVectorTakeArray(y_ceed, MemTypeP2C(y_mem_type), NULL);
211   ierr = VecRestoreArrayAndMemType(Y_loc, &y); CHKERRQ(ierr);
212 
213   // -- Local-to-global
214   ierr = VecZeroEntries(mult_vec); CHKERRQ(ierr);
215   ierr = DMLocalToGlobal(user->dm, Y_loc, ADD_VALUES, mult_vec);
216   CHKERRQ(ierr);
217 
218   // -- Scale
219   ierr = VecReciprocal(mult_vec); CHKERRQ(ierr);
220   ierr = VecPointwiseMult(Diagnostic, Diagnostic, mult_vec);
221 
222   // ---------------------------------------------------------------------------
223   // Output solution vector
224   // ---------------------------------------------------------------------------
225   ierr = PetscSNPrintf(output_filename, sizeof output_filename,
226                        "%s/diagnostic_quantities.vtu",
227                        app_ctx->output_dir); CHKERRQ(ierr);
228   ierr = PetscViewerVTKOpen(comm, output_filename, FILE_MODE_WRITE, &viewer);
229   CHKERRQ(ierr);
230   ierr = VecView(Diagnostic, viewer); CHKERRQ(ierr);
231   ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
232 
233   // ---------------------------------------------------------------------------
234   // Cleanup
235   // ---------------------------------------------------------------------------
236   ierr = VecDestroy(&Diagnostic); CHKERRQ(ierr);
237   ierr = VecDestroy(&mult_vec); CHKERRQ(ierr);
238   ierr = VecDestroy(&Y_loc); CHKERRQ(ierr);
239   CeedVectorDestroy(&y_ceed);
240 
241   PetscFunctionReturn(0);
242 };
243 
244 // -----------------------------------------------------------------------------
245 // Regression testing
246 // -----------------------------------------------------------------------------
247 // test option change. could remove the loading step. Run only with one loading step and compare relatively to ref file
248 // 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
249 PetscErrorCode RegressionTests_solids(AppCtx app_ctx, PetscReal energy) {
250   PetscFunctionBegin;
251 
252   if (app_ctx->expect_final_strain >= 0.) {
253     PetscReal energy_ref = app_ctx->expect_final_strain;
254     PetscReal error = PetscAbsReal(energy - energy_ref) / energy_ref;
255 
256     if (error > app_ctx->test_tol) {
257       PetscErrorCode ierr;
258       ierr = PetscPrintf(PETSC_COMM_WORLD,
259                          "Energy %e does not match expected energy %e: relative tolerance %e > %e\n",
260                          (double)energy, (double)energy_ref, (double)error, app_ctx->test_tol);
261       CHKERRQ(ierr);
262     }
263   }
264   PetscFunctionReturn(0);
265 };
266