1 // Copyright (c) 2017-2026, 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 #include "../include/swarmutils.h"
9 #include "../include/matops.h"
10 #include "../qfunctions/swarm/swarmmass.h"
11
12 // ------------------------------------------------------------------------------------------------
13 // Context utilities
14 // ------------------------------------------------------------------------------------------------
DMSwarmCeedContextCreate(DM dm_swarm,const char * ceed_resource,DMSwarmCeedContext * ctx)15 PetscErrorCode DMSwarmCeedContextCreate(DM dm_swarm, const char *ceed_resource, DMSwarmCeedContext *ctx) {
16 DM dm_mesh, dm_coord;
17 CeedElemRestriction elem_restr_u_mesh, elem_restr_x_mesh, elem_restr_x_points, elem_restr_u_points, elem_restr_q_data_points;
18 CeedBasis basis_u, basis_x;
19 CeedVector x_ref_points, q_data_points;
20 CeedInt num_comp;
21
22 PetscFunctionBeginUser;
23 PetscCall(PetscNew(ctx));
24 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh));
25 PetscCall(DMGetCoordinateDM(dm_mesh, &dm_coord));
26
27 CeedInit(ceed_resource, &(*ctx)->ceed);
28 // Background mesh objects
29 {
30 BPData bp_data = {.q_mode = CEED_GAUSS};
31
32 PetscCall(CreateBasisFromPlex((*ctx)->ceed, dm_mesh, NULL, 0, 0, 0, bp_data, &basis_u));
33 PetscCall(CreateBasisFromPlex((*ctx)->ceed, dm_coord, NULL, 0, 0, 0, bp_data, &basis_x));
34 PetscCall(CreateRestrictionFromPlex((*ctx)->ceed, dm_mesh, 0, NULL, 0, &elem_restr_u_mesh));
35 PetscCall(CreateRestrictionFromPlex((*ctx)->ceed, dm_coord, 0, NULL, 0, &elem_restr_x_mesh));
36
37 // -- Mesh vectors
38 CeedElemRestrictionCreateVector(elem_restr_u_mesh, &(*ctx)->u_mesh, NULL);
39 CeedElemRestrictionCreateVector(elem_restr_u_mesh, &(*ctx)->v_mesh, NULL);
40 }
41 // Swarm objects
42 {
43 PetscInt dim;
44 const PetscInt *cell_points;
45 IS is_points;
46 Vec X_ref;
47 CeedInt num_elem;
48
49 PetscCall(DMSwarmCreateReferenceCoordinates(dm_swarm, &is_points, &X_ref));
50 PetscCall(DMGetDimension(dm_mesh, &dim));
51 CeedElemRestrictionGetNumElements(elem_restr_u_mesh, &num_elem);
52 CeedElemRestrictionGetNumComponents(elem_restr_u_mesh, &num_comp);
53
54 PetscCall(ISGetIndices(is_points, &cell_points));
55 PetscInt num_points = cell_points[num_elem + 1] - num_elem - 2;
56 CeedInt offsets[num_elem + 1 + num_points];
57
58 for (PetscInt i = 0; i < num_elem + 1; i++) offsets[i] = cell_points[i + 1] - 1;
59 for (PetscInt i = num_elem + 1; i < num_points + num_elem + 1; i++) offsets[i] = cell_points[i + 1];
60 PetscCall(ISRestoreIndices(is_points, &cell_points));
61
62 // -- Points restrictions
63 CeedElemRestrictionCreateAtPoints((*ctx)->ceed, num_elem, num_points, num_comp, num_points * num_comp, CEED_MEM_HOST, CEED_COPY_VALUES, offsets,
64 &elem_restr_u_points);
65 CeedElemRestrictionCreateAtPoints((*ctx)->ceed, num_elem, num_points, dim, num_points * dim, CEED_MEM_HOST, CEED_COPY_VALUES, offsets,
66 &elem_restr_x_points);
67 CeedElemRestrictionCreateAtPoints((*ctx)->ceed, num_elem, num_points, 1, num_points, CEED_MEM_HOST, CEED_COPY_VALUES, offsets,
68 &elem_restr_q_data_points);
69
70 // -- Points vectors
71 CeedElemRestrictionCreateVector(elem_restr_u_points, &(*ctx)->u_points, NULL);
72 CeedElemRestrictionCreateVector(elem_restr_q_data_points, &q_data_points, NULL);
73
74 // -- Ref coordinates
75 {
76 PetscMemType X_mem_type;
77 const PetscScalar *x;
78
79 CeedVectorCreate((*ctx)->ceed, num_points * dim, &x_ref_points);
80
81 PetscCall(VecGetArrayReadAndMemType(X_ref, (const PetscScalar **)&x, &X_mem_type));
82 CeedVectorSetArray(x_ref_points, MemTypeP2C(X_mem_type), CEED_COPY_VALUES, (CeedScalar *)x);
83 PetscCall(VecRestoreArrayReadAndMemType(X_ref, (const PetscScalar **)&x));
84 }
85
86 // Create Q data
87 {
88 CeedQFunction qf_setup;
89 CeedOperator op_setup;
90 CeedVector x_coord;
91
92 {
93 Vec X_loc;
94 PetscInt len;
95 const PetscScalar *x;
96
97 PetscCall(DMGetCoordinatesLocal(dm_mesh, &X_loc));
98 PetscCall(VecGetLocalSize(X_loc, &len));
99 CeedVectorCreate((*ctx)->ceed, len, &x_coord);
100
101 PetscCall(VecGetArrayRead(X_loc, &x));
102 CeedVectorSetArray(x_coord, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)x);
103 PetscCall(VecRestoreArrayRead(X_loc, &x));
104 }
105
106 // Setup geometric scaling
107 CeedQFunctionCreateInterior((*ctx)->ceed, 1, SetupMass, SetupMass_loc, &qf_setup);
108 CeedQFunctionAddInput(qf_setup, "x", dim * dim, CEED_EVAL_GRAD);
109 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT);
110 CeedQFunctionAddOutput(qf_setup, "rho", 1, CEED_EVAL_NONE);
111
112 CeedOperatorCreateAtPoints((*ctx)->ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup);
113 CeedOperatorSetField(op_setup, "x", elem_restr_x_mesh, basis_x, CEED_VECTOR_ACTIVE);
114 CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
115 CeedOperatorSetField(op_setup, "rho", elem_restr_q_data_points, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
116 CeedOperatorAtPointsSetPoints(op_setup, elem_restr_x_points, x_ref_points);
117
118 CeedOperatorApply(op_setup, x_coord, q_data_points, CEED_REQUEST_IMMEDIATE);
119
120 // Cleanup
121 CeedVectorDestroy(&x_coord);
122 CeedQFunctionDestroy(&qf_setup);
123 CeedOperatorDestroy(&op_setup);
124 }
125
126 // -- Cleanup
127 PetscCall(ISDestroy(&is_points));
128 PetscCall(VecDestroy(&X_ref));
129 }
130
131 PetscCall(DMSetApplicationContext(dm_mesh, (void *)(*ctx)));
132
133 // Create operators
134 // Mesh to points interpolation operator
135 {
136 CeedQFunction qf_mesh_to_points;
137
138 // -- Create operator
139 CeedQFunctionCreateIdentity((*ctx)->ceed, num_comp, CEED_EVAL_INTERP, CEED_EVAL_NONE, &qf_mesh_to_points);
140
141 CeedOperatorCreateAtPoints((*ctx)->ceed, qf_mesh_to_points, NULL, NULL, &(*ctx)->op_mesh_to_points);
142 CeedOperatorSetField((*ctx)->op_mesh_to_points, "input", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE);
143 CeedOperatorSetField((*ctx)->op_mesh_to_points, "output", elem_restr_u_points, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
144 CeedOperatorAtPointsSetPoints((*ctx)->op_mesh_to_points, elem_restr_x_points, x_ref_points);
145
146 // -- Cleanup
147 CeedQFunctionDestroy(&qf_mesh_to_points);
148 }
149
150 // RHS operator
151 {
152 CeedQFunction qf_pts_to_mesh;
153 CeedQFunctionContext qf_ctx;
154
155 // -- Mass QFunction
156 CeedQFunctionCreateInterior((*ctx)->ceed, 1, Mass, Mass_loc, &qf_pts_to_mesh);
157 CeedQFunctionAddInput(qf_pts_to_mesh, "q data", 1, CEED_EVAL_NONE);
158 CeedQFunctionAddInput(qf_pts_to_mesh, "u", num_comp, CEED_EVAL_NONE);
159 CeedQFunctionAddOutput(qf_pts_to_mesh, "v", num_comp, CEED_EVAL_INTERP);
160
161 // -- QFunction context
162 CeedQFunctionContextCreate((*ctx)->ceed, &qf_ctx);
163 CeedQFunctionContextSetData(qf_ctx, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof(num_comp), &num_comp);
164 CeedQFunctionSetContext(qf_pts_to_mesh, qf_ctx);
165
166 // -- Mass Operator
167 CeedOperatorCreateAtPoints((*ctx)->ceed, qf_pts_to_mesh, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &(*ctx)->op_points_to_mesh);
168 CeedOperatorSetField((*ctx)->op_points_to_mesh, "q data", elem_restr_q_data_points, CEED_BASIS_NONE, q_data_points);
169 CeedOperatorSetField((*ctx)->op_points_to_mesh, "u", elem_restr_u_points, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
170 CeedOperatorSetField((*ctx)->op_points_to_mesh, "v", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE);
171 CeedOperatorAtPointsSetPoints((*ctx)->op_points_to_mesh, elem_restr_x_points, x_ref_points);
172
173 // -- Cleanup
174 CeedQFunctionContextDestroy(&qf_ctx);
175 CeedQFunctionDestroy(&qf_pts_to_mesh);
176 }
177
178 // Mass operator
179 {
180 CeedQFunction qf_mass;
181 CeedQFunctionContext ctx_mass;
182
183 // -- Mass QFunction
184 CeedQFunctionCreateInterior((*ctx)->ceed, 1, Mass, Mass_loc, &qf_mass);
185 CeedQFunctionAddInput(qf_mass, "q data", 1, CEED_EVAL_NONE);
186 CeedQFunctionAddInput(qf_mass, "u", num_comp, CEED_EVAL_INTERP);
187 CeedQFunctionAddOutput(qf_mass, "v", num_comp, CEED_EVAL_INTERP);
188
189 // -- QFunction context
190 CeedQFunctionContextCreate((*ctx)->ceed, &ctx_mass);
191 CeedQFunctionContextSetData(ctx_mass, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof(num_comp), &num_comp);
192 CeedQFunctionSetContext(qf_mass, ctx_mass);
193
194 // -- Mass Operator
195 CeedOperatorCreateAtPoints((*ctx)->ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &(*ctx)->op_mass);
196 CeedOperatorSetField((*ctx)->op_mass, "q data", elem_restr_q_data_points, CEED_BASIS_NONE, q_data_points);
197 CeedOperatorSetField((*ctx)->op_mass, "u", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE);
198 CeedOperatorSetField((*ctx)->op_mass, "v", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE);
199 CeedOperatorAtPointsSetPoints((*ctx)->op_mass, elem_restr_x_points, x_ref_points);
200
201 // -- Cleanup
202 CeedQFunctionContextDestroy(&ctx_mass);
203 CeedQFunctionDestroy(&qf_mass);
204 }
205
206 // Cleanup
207 CeedElemRestrictionDestroy(&elem_restr_u_mesh);
208 CeedElemRestrictionDestroy(&elem_restr_x_mesh);
209 CeedElemRestrictionDestroy(&elem_restr_u_points);
210 CeedElemRestrictionDestroy(&elem_restr_x_points);
211 CeedElemRestrictionDestroy(&elem_restr_q_data_points);
212 CeedBasisDestroy(&basis_u);
213 CeedBasisDestroy(&basis_x);
214 CeedVectorDestroy(&x_ref_points);
215 CeedVectorDestroy(&q_data_points);
216 PetscFunctionReturn(PETSC_SUCCESS);
217 }
218
DMSwarmCeedContextDestroy(DMSwarmCeedContext * ctx)219 PetscErrorCode DMSwarmCeedContextDestroy(DMSwarmCeedContext *ctx) {
220 PetscFunctionBeginUser;
221 CeedDestroy(&(*ctx)->ceed);
222 CeedVectorDestroy(&(*ctx)->u_mesh);
223 CeedVectorDestroy(&(*ctx)->v_mesh);
224 CeedVectorDestroy(&(*ctx)->u_points);
225 CeedOperatorDestroy(&(*ctx)->op_mesh_to_points);
226 CeedOperatorDestroy(&(*ctx)->op_points_to_mesh);
227 CeedOperatorDestroy(&(*ctx)->op_mass);
228 PetscCall(PetscFree(*ctx));
229 PetscFunctionReturn(PETSC_SUCCESS);
230 }
231
232 // ------------------------------------------------------------------------------------------------
233 // PETSc-libCEED memory space utilities
234 // ------------------------------------------------------------------------------------------------
DMSwarmPICFieldP2C(DM dm_swarm,const char * field,CeedVector x_ceed)235 PetscErrorCode DMSwarmPICFieldP2C(DM dm_swarm, const char *field, CeedVector x_ceed) {
236 PetscScalar *x;
237
238 PetscFunctionBeginUser;
239 PetscCall(DMSwarmGetField(dm_swarm, field, NULL, NULL, (void **)&x));
240 CeedVectorSetArray(x_ceed, CEED_MEM_HOST, CEED_USE_POINTER, (CeedScalar *)x);
241 PetscFunctionReturn(PETSC_SUCCESS);
242 }
243
DMSwarmPICFieldC2P(DM dm_swarm,const char * field,CeedVector x_ceed)244 PetscErrorCode DMSwarmPICFieldC2P(DM dm_swarm, const char *field, CeedVector x_ceed) {
245 PetscScalar *x;
246
247 PetscFunctionBeginUser;
248 CeedVectorTakeArray(x_ceed, CEED_MEM_HOST, (CeedScalar **)&x);
249 PetscCall(DMSwarmRestoreField(dm_swarm, field, NULL, NULL, (void **)&x));
250 PetscFunctionReturn(PETSC_SUCCESS);
251 }
252
253 // ------------------------------------------------------------------------------------------------
254 // Swarm point location utility
255 // ------------------------------------------------------------------------------------------------
DMSwarmInitalizePointLocations(DM dm_swarm,PointSwarmType point_swarm_type,PetscInt num_points,PetscInt num_points_per_cell)256 PetscErrorCode DMSwarmInitalizePointLocations(DM dm_swarm, PointSwarmType point_swarm_type, PetscInt num_points, PetscInt num_points_per_cell) {
257 PetscFunctionBeginUser;
258 switch (point_swarm_type) {
259 case SWARM_GAUSS:
260 case SWARM_UNIFORM: {
261 // -- Set gauss or uniform point locations in each cell
262 PetscInt num_points_per_cell_1d = round(cbrt(num_points_per_cell * 1.0)), dim = 3;
263 PetscScalar point_coords[num_points_per_cell * 3];
264 CeedScalar points_1d[num_points_per_cell_1d], weights_1d[num_points_per_cell_1d];
265
266 if (point_swarm_type == SWARM_GAUSS) {
267 PetscCall(CeedGaussQuadrature(num_points_per_cell_1d, points_1d, weights_1d));
268 } else {
269 for (PetscInt i = 0; i < num_points_per_cell_1d; i++) points_1d[i] = 2.0 * (PetscReal)(i + 0.5) / (PetscReal)num_points_per_cell_1d - 1;
270 }
271 for (PetscInt i = 0; i < num_points_per_cell_1d; i++) {
272 for (PetscInt j = 0; j < num_points_per_cell_1d; j++) {
273 for (PetscInt k = 0; k < num_points_per_cell_1d; k++) {
274 PetscInt p = (i * num_points_per_cell_1d + j) * num_points_per_cell_1d + k;
275
276 point_coords[p * dim + 0] = points_1d[i];
277 point_coords[p * dim + 1] = points_1d[j];
278 point_coords[p * dim + 2] = points_1d[k];
279 }
280 }
281 }
282 PetscCall(DMSwarmSetPointCoordinatesCellwise(dm_swarm, num_points_per_cell_1d * num_points_per_cell_1d * num_points_per_cell_1d, point_coords));
283 } break;
284 case SWARM_CELL_RANDOM: {
285 DM dm_mesh;
286
287 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh));
288 // DMSwarmSetPointCoordinatesRandom expects local coordinates to be set up to ensure call is non-collective
289 PetscCall(DMGetCoordinatesLocalSetUp(dm_mesh));
290 PetscCall(DMSwarmSetPointCoordinatesRandom(dm_swarm, num_points_per_cell));
291 } break;
292 case SWARM_SINUSOIDAL: {
293 // -- Set points distributed per sinusoidal functions
294 PetscInt dim = 3;
295 PetscScalar *point_coords;
296 DM dm_mesh;
297
298 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh));
299 PetscCall(DMGetDimension(dm_mesh, &dim));
300
301 PetscCall(DMSwarmGetField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&point_coords));
302 for (PetscInt p = 0; p < num_points; p++) {
303 point_coords[p * dim + 0] = -PetscCosReal((PetscReal)(p + 1) / (PetscReal)(num_points + 1) * PETSC_PI);
304 if (dim > 1) point_coords[p * dim + 1] = -PetscSinReal((PetscReal)(p + 1) / (PetscReal)(num_points + 1) * PETSC_PI);
305 if (dim > 2) point_coords[p * dim + 2] = PetscSinReal((PetscReal)(p + 1) / (PetscReal)(num_points + 1) * PETSC_PI);
306 }
307 PetscCall(DMSwarmRestoreField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&point_coords));
308 } break;
309 }
310 PetscCall(DMSwarmMigrate(dm_swarm, PETSC_TRUE));
311 PetscFunctionReturn(PETSC_SUCCESS);
312 }
313
314 /*@C
315 DMSwarmCreateReferenceCoordinates - Compute the cell reference coordinates for local DMSwarm points.
316
317 Collective
318
319 Input Parameter:
320 . dm_swarm - the `DMSwarm`
321
322 Output Parameters:
323 + is_points - The IS object for indexing into points per cell
324 - X_points_ref - Vec holding the cell reference coordinates for local DMSwarm points
325
326 The index set contains ranges of indices for each local cell. This range contains the indices of every point in the cell.
327
328 ```
329 total_num_cells
330 cell_0_start_index
331 cell_1_start_index
332 cell_2_start_index
333 ...
334 cell_n_start_index
335 cell_n_stop_index
336 cell_0_point_0
337 cell_0_point_0
338 ...
339 cell_n_point_m
340 ```
341
342 Level: beginner
343
344 .seealso: `DMSwarm`
345 @*/
DMSwarmCreateReferenceCoordinates(DM dm_swarm,IS * is_points,Vec * X_points_ref)346 PetscErrorCode DMSwarmCreateReferenceCoordinates(DM dm_swarm, IS *is_points, Vec *X_points_ref) {
347 PetscInt cell_start, cell_end, num_cells_local, dim, num_points_local, *cell_points, points_offset;
348 PetscScalar *coords_points_ref;
349 const PetscScalar *coords_points_true;
350 DM dm_mesh;
351
352 PetscFunctionBeginUser;
353 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh));
354
355 // Create vector to hold reference coordinates
356 {
357 Vec X_points_true;
358
359 PetscCall(DMSwarmCreateLocalVectorFromField(dm_swarm, DMSwarmPICField_coor, &X_points_true));
360 PetscCall(VecDuplicate(X_points_true, X_points_ref));
361 PetscCall(DMSwarmDestroyLocalVectorFromField(dm_swarm, DMSwarmPICField_coor, &X_points_true));
362 }
363
364 // Allocate index set array
365 PetscCall(DMPlexGetHeightStratum(dm_mesh, 0, &cell_start, &cell_end));
366 num_cells_local = cell_end - cell_start;
367 points_offset = num_cells_local + 2;
368 PetscCall(VecGetLocalSize(*X_points_ref, &num_points_local));
369 PetscCall(DMGetDimension(dm_mesh, &dim));
370 num_points_local /= dim;
371 PetscCall(PetscMalloc1(num_points_local + num_cells_local + 2, &cell_points));
372 cell_points[0] = num_cells_local;
373
374 // Get reference coordinates for each swarm point wrt the elements in the background mesh
375 PetscCall(DMSwarmSortGetAccess(dm_swarm));
376 PetscCall(DMSwarmGetField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords_points_true));
377 PetscCall(VecGetArray(*X_points_ref, &coords_points_ref));
378 for (PetscInt cell = cell_start, num_points_processed = 0; cell < cell_end; cell++) {
379 PetscInt *points_in_cell, num_points_in_cell, local_cell = cell - cell_start;
380 PetscReal v[3], J[9], invJ[9], detJ, v0_ref[3] = {-1.0, -1.0, -1.0};
381
382 PetscCall(DMSwarmSortGetPointsPerCell(dm_swarm, cell, &num_points_in_cell, &points_in_cell));
383 // -- Reference coordinates for swarm points in background mesh element
384 PetscCall(DMPlexComputeCellGeometryFEM(dm_mesh, cell, NULL, v, J, invJ, &detJ));
385 cell_points[local_cell + 1] = num_points_processed + points_offset;
386 for (PetscInt p = 0; p < num_points_in_cell; p++) {
387 const CeedInt point = points_in_cell[p];
388
389 cell_points[points_offset + (num_points_processed++)] = point;
390 CoordinatesRealToRef(dim, dim, v0_ref, v, invJ, &coords_points_true[point * dim], &coords_points_ref[point * dim]);
391 }
392
393 // -- Cleanup
394 PetscCall(DMSwarmSortRestorePointsPerCell(dm_swarm, cell, &num_points_in_cell, &points_in_cell));
395 }
396 cell_points[points_offset - 1] = num_points_local + points_offset;
397
398 // Cleanup
399 PetscCall(DMSwarmRestoreField(dm_swarm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords_points_true));
400 PetscCall(VecRestoreArray(*X_points_ref, &coords_points_ref));
401 PetscCall(DMSwarmSortRestoreAccess(dm_swarm));
402
403 // Create index set
404 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, num_points_local + points_offset, cell_points, PETSC_OWN_POINTER, is_points));
405 PetscFunctionReturn(PETSC_SUCCESS);
406 }
407
408 // ------------------------------------------------------------------------------------------------
409 // RHS for Swarm to Mesh projection
410 // ------------------------------------------------------------------------------------------------
DMSwarmCreateProjectionRHS(DM dm_swarm,const char * field,Vec U_points,Vec B_mesh)411 PetscErrorCode DMSwarmCreateProjectionRHS(DM dm_swarm, const char *field, Vec U_points, Vec B_mesh) {
412 PetscMemType B_mem_type, U_mem_type;
413 DM dm_mesh;
414 Vec B_mesh_loc;
415 PetscBool has_u_points;
416 DMSwarmCeedContext swarm_ceed_context;
417
418 PetscFunctionBeginUser;
419 // Get mesh DM
420 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh));
421 PetscCall(DMGetApplicationContext(dm_mesh, (void *)&swarm_ceed_context));
422
423 // Get swarm access
424 has_u_points = !!U_points;
425 if (!has_u_points) {
426 PetscCall(DMSwarmSortGetAccess(dm_swarm));
427 PetscCall(DMSwarmCreateLocalVectorFromField(dm_swarm, field, &U_points));
428 }
429
430 // Get mesh values
431 PetscCall(VecReadP2C(U_points, &U_mem_type, swarm_ceed_context->u_points));
432 PetscCall(DMGetLocalVector(dm_mesh, &B_mesh_loc));
433 PetscCall(VecZeroEntries(B_mesh_loc));
434 PetscCall(VecP2C(B_mesh_loc, &B_mem_type, swarm_ceed_context->v_mesh));
435
436 // Interpolate field from swarm points to mesh
437 CeedOperatorApply(swarm_ceed_context->op_points_to_mesh, swarm_ceed_context->u_points, swarm_ceed_context->v_mesh, CEED_REQUEST_IMMEDIATE);
438
439 // Restore PETSc Vecs and Local to Global
440 PetscCall(VecReadC2P(swarm_ceed_context->u_points, U_mem_type, U_points));
441 PetscCall(VecC2P(swarm_ceed_context->v_mesh, B_mem_type, B_mesh_loc));
442 PetscCall(VecZeroEntries(B_mesh));
443 PetscCall(DMLocalToGlobal(dm_mesh, B_mesh_loc, ADD_VALUES, B_mesh));
444
445 // Restore swarm access
446 if (!has_u_points) {
447 PetscCall(DMSwarmDestroyLocalVectorFromField(dm_swarm, field, &U_points));
448 PetscCall(DMSwarmSortRestoreAccess(dm_swarm));
449 }
450
451 // Cleanup
452 PetscCall(DMRestoreLocalVector(dm_mesh, &B_mesh_loc));
453 PetscFunctionReturn(PETSC_SUCCESS);
454 }
455
456 // ------------------------------------------------------------------------------------------------
457 // Swarm "mass matrix"
458 // ------------------------------------------------------------------------------------------------
MatMult_SwarmMass(Mat A,Vec U_mesh,Vec V_mesh)459 PetscErrorCode MatMult_SwarmMass(Mat A, Vec U_mesh, Vec V_mesh) {
460 PetscMemType U_mem_type, V_mem_type;
461 DM dm_mesh;
462 Vec U_mesh_loc, V_mesh_loc;
463 DMSwarmCeedContext swarm_ceed_context;
464
465 PetscFunctionBeginUser;
466 // Get mesh DM
467 PetscCall(MatGetDM(A, &dm_mesh));
468 PetscCall(DMGetApplicationContext(dm_mesh, (void *)&swarm_ceed_context));
469
470 // Global to Local and get PETSc Vec access
471 PetscCall(DMGetLocalVector(dm_mesh, &U_mesh_loc));
472 PetscCall(VecZeroEntries(U_mesh_loc));
473 PetscCall(DMGlobalToLocal(dm_mesh, U_mesh, INSERT_VALUES, U_mesh_loc));
474 PetscCall(VecReadP2C(U_mesh_loc, &U_mem_type, swarm_ceed_context->u_mesh));
475 PetscCall(DMGetLocalVector(dm_mesh, &V_mesh_loc));
476 PetscCall(VecZeroEntries(V_mesh_loc));
477 PetscCall(VecP2C(V_mesh_loc, &V_mem_type, swarm_ceed_context->v_mesh));
478
479 // Apply swarm mass operator
480 CeedOperatorApply(swarm_ceed_context->op_mass, swarm_ceed_context->u_mesh, swarm_ceed_context->v_mesh, CEED_REQUEST_IMMEDIATE);
481
482 // Restore PETSc Vecs and Local to Global
483 PetscCall(VecReadC2P(swarm_ceed_context->u_mesh, U_mem_type, U_mesh_loc));
484 PetscCall(VecC2P(swarm_ceed_context->v_mesh, V_mem_type, V_mesh_loc));
485 PetscCall(VecZeroEntries(V_mesh));
486 PetscCall(DMLocalToGlobal(dm_mesh, V_mesh_loc, ADD_VALUES, V_mesh));
487
488 // Cleanup
489 PetscCall(DMRestoreLocalVector(dm_mesh, &U_mesh_loc));
490 PetscCall(DMRestoreLocalVector(dm_mesh, &V_mesh_loc));
491 PetscFunctionReturn(PETSC_SUCCESS);
492 }
493
494 // ------------------------------------------------------------------------------------------------
495 // Swarm to mesh projection
496 // ------------------------------------------------------------------------------------------------
DMSwarmProjectFromSwarmToCells(DM dm_swarm,const char * field,Vec U_points,Vec U_mesh)497 PetscErrorCode DMSwarmProjectFromSwarmToCells(DM dm_swarm, const char *field, Vec U_points, Vec U_mesh) {
498 PetscBool test_mode;
499 Vec B_mesh;
500 Mat M;
501 KSP ksp;
502 DM dm_mesh;
503 DMSwarmCeedContext swarm_ceed_context;
504 MPI_Comm comm;
505
506 PetscFunctionBeginUser;
507 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Swarm-to-Mesh Projection Options", NULL);
508 PetscCall(PetscOptionsBool("-test", "Testing mode (do not print unless error is large)", NULL, test_mode, &test_mode, NULL));
509 PetscOptionsEnd();
510
511 comm = PetscObjectComm((PetscObject)dm_swarm);
512 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh));
513 PetscCall(DMGetApplicationContext(dm_mesh, (void *)&swarm_ceed_context));
514 PetscCall(VecDuplicate(U_mesh, &B_mesh));
515
516 // Setup "mass matrix"
517 {
518 PetscInt l_size, g_size;
519
520 PetscCall(VecGetLocalSize(U_mesh, &l_size));
521 PetscCall(VecGetSize(U_mesh, &g_size));
522 PetscCall(MatCreateShell(comm, l_size, l_size, g_size, g_size, swarm_ceed_context, &M));
523 PetscCall(MatSetDM(M, dm_mesh));
524 PetscCall(MatShellSetOperation(M, MATOP_MULT, (void (*)(void))MatMult_SwarmMass));
525 }
526
527 // Setup KSP
528 {
529 PC pc;
530
531 PetscCall(KSPCreate(comm, &ksp));
532 PetscCall(KSPGetPC(ksp, &pc));
533 PetscCall(PCSetType(pc, PCJACOBI));
534 PetscCall(PCJacobiSetType(pc, PC_JACOBI_ROWSUM));
535 PetscCall(KSPSetType(ksp, KSPCG));
536 PetscCall(KSPSetNormType(ksp, KSP_NORM_NATURAL));
537 PetscCall(KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
538 PetscCall(KSPSetOperators(ksp, M, M));
539 PetscCall(KSPSetFromOptions(ksp));
540 PetscCall(PetscObjectSetName((PetscObject)ksp, "Swarm-to-Mesh Projection"));
541 PetscCall(KSPViewFromOptions(ksp, NULL, "-ksp_projection_view"));
542 }
543
544 // Setup RHS
545 PetscCall(DMSwarmCreateProjectionRHS(dm_swarm, field, U_points, B_mesh));
546
547 // Solve
548 PetscCall(VecZeroEntries(U_mesh));
549 PetscCall(KSPSolve(ksp, B_mesh, U_mesh));
550
551 // KSP summary
552 {
553 KSPType ksp_type;
554 KSPConvergedReason reason;
555 PetscReal rnorm;
556 PetscInt its;
557 PetscCall(KSPGetType(ksp, &ksp_type));
558 PetscCall(KSPGetConvergedReason(ksp, &reason));
559 PetscCall(KSPGetIterationNumber(ksp, &its));
560 PetscCall(KSPGetResidualNorm(ksp, &rnorm));
561
562 if (!test_mode || reason < 0 || rnorm > 1e-8) {
563 PetscCall(PetscPrintf(comm,
564 "Swarm-to-Mesh Projection KSP Solve:\n"
565 " KSP type: %s\n"
566 " KSP convergence: %s\n"
567 " Total KSP iterations: %" PetscInt_FMT "\n"
568 " Final rnorm: %e\n",
569 ksp_type, KSPConvergedReasons[reason], its, (double)rnorm));
570 }
571 }
572
573 // Optional viewing
574 PetscCall(KSPViewFromOptions(ksp, NULL, "-ksp_view"));
575
576 // Cleanup
577 PetscCall(VecDestroy(&B_mesh));
578 PetscCall(MatDestroy(&M));
579 PetscCall(KSPDestroy(&ksp));
580 PetscFunctionReturn(PETSC_SUCCESS);
581 }
582
583 // ------------------------------------------------------------------------------------------------
584 // BP setup
585 // ------------------------------------------------------------------------------------------------
SetupProblemSwarm(DM dm_swarm,Ceed ceed,BPData bp_data,CeedData data,PetscBool setup_rhs,Vec rhs,Vec target)586 PetscErrorCode SetupProblemSwarm(DM dm_swarm, Ceed ceed, BPData bp_data, CeedData data, PetscBool setup_rhs, Vec rhs, Vec target) {
587 DM dm_mesh, dm_coord;
588 CeedElemRestriction elem_restr_u_mesh, elem_restr_x_mesh, elem_restr_x_points, elem_restr_u_points, elem_restr_q_data_points;
589 CeedBasis basis_u, basis_x;
590 CeedVector x_coord, x_ref_points, q_data_points;
591 PetscInt X_loc_len, dim;
592 CeedInt num_comp, q_data_size = bp_data.q_data_size;
593 CeedScalar R = 1; // radius of the sphere
594 CeedScalar l = 1.0 / PetscSqrtReal(3.0); // half edge of the inscribed cube
595 Vec X_loc;
596 PetscMemType X_mem_type;
597
598 PetscFunctionBeginUser;
599 PetscCall(DMSwarmGetCellDM(dm_swarm, &dm_mesh));
600 PetscCall(DMGetCoordinateDM(dm_mesh, &dm_coord));
601
602 // Get coordinates
603 PetscCall(DMGetCoordinatesLocal(dm_mesh, &X_loc));
604 PetscCall(VecGetLocalSize(X_loc, &X_loc_len));
605 CeedVectorCreate(ceed, X_loc_len, &x_coord);
606 PetscCall(VecReadP2C(X_loc, &X_mem_type, x_coord));
607
608 // Background mesh objects
609 PetscCall(CreateBasisFromPlex(ceed, dm_mesh, NULL, 0, 0, 0, bp_data, &basis_u));
610 PetscCall(CreateBasisFromPlex(ceed, dm_coord, NULL, 0, 0, 0, bp_data, &basis_x));
611 PetscCall(CreateRestrictionFromPlex(ceed, dm_mesh, 0, NULL, 0, &elem_restr_u_mesh));
612 PetscCall(CreateRestrictionFromPlex(ceed, dm_coord, 0, NULL, 0, &elem_restr_x_mesh));
613
614 CeedElemRestrictionCreateVector(elem_restr_u_mesh, &data->x_ceed, NULL);
615 CeedElemRestrictionCreateVector(elem_restr_u_mesh, &data->y_ceed, NULL);
616
617 // Swarm objects
618 {
619 const PetscInt *cell_points;
620 CeedInt *offsets;
621 IS is_points;
622 Vec X_ref;
623 CeedInt num_elem;
624
625 PetscCall(DMSwarmCreateReferenceCoordinates(dm_swarm, &is_points, &X_ref));
626 PetscCall(DMGetDimension(dm_mesh, &dim));
627 CeedElemRestrictionGetNumElements(elem_restr_u_mesh, &num_elem);
628 CeedElemRestrictionGetNumComponents(elem_restr_u_mesh, &num_comp);
629
630 PetscCall(ISGetIndices(is_points, &cell_points));
631 PetscInt num_points = cell_points[num_elem + 1] - num_elem - 2;
632 PetscCall(PetscCalloc1(num_elem + 1 + num_points, &offsets));
633
634 for (PetscInt i = 0; i < num_elem + 1; i++) offsets[i] = cell_points[i + 1] - 1;
635 for (PetscInt i = num_elem + 1; i < num_points + num_elem + 1; i++) offsets[i] = cell_points[i + 1];
636 PetscCall(ISRestoreIndices(is_points, &cell_points));
637
638 // -- Points restrictions
639 CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, num_comp, num_points * num_comp, CEED_MEM_HOST, CEED_COPY_VALUES, offsets,
640 &elem_restr_u_points);
641 CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, dim, num_points * dim, CEED_MEM_HOST, CEED_COPY_VALUES, offsets,
642 &elem_restr_x_points);
643 CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points, q_data_size, num_points * q_data_size, CEED_MEM_HOST, CEED_COPY_VALUES, offsets,
644 &elem_restr_q_data_points);
645
646 // -- Points vectors
647 CeedElemRestrictionCreateVector(elem_restr_q_data_points, &q_data_points, NULL);
648
649 // -- Ref coordinates
650 {
651 PetscMemType X_mem_type;
652 const PetscScalar *x;
653
654 CeedVectorCreate(ceed, num_points * dim, &x_ref_points);
655
656 PetscCall(VecGetArrayReadAndMemType(X_ref, (const PetscScalar **)&x, &X_mem_type));
657 CeedVectorSetArray(x_ref_points, MemTypeP2C(X_mem_type), CEED_COPY_VALUES, (CeedScalar *)x);
658 PetscCall(VecRestoreArrayReadAndMemType(X_ref, (const PetscScalar **)&x));
659 }
660
661 // Create Q data
662 {
663 CeedQFunction qf_setup;
664 CeedOperator op_setup;
665
666 // Setup geometric scaling
667 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_geo, bp_data.setup_geo_loc, &qf_setup);
668 CeedQFunctionAddInput(qf_setup, "x", dim, CEED_EVAL_INTERP);
669 CeedQFunctionAddInput(qf_setup, "dx", dim * dim, CEED_EVAL_GRAD);
670 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT);
671 CeedQFunctionAddOutput(qf_setup, "qdata", q_data_size, CEED_EVAL_NONE);
672
673 CeedOperatorCreateAtPoints(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup);
674 CeedOperatorSetField(op_setup, "x", elem_restr_x_mesh, basis_x, CEED_VECTOR_ACTIVE);
675 CeedOperatorSetField(op_setup, "dx", elem_restr_x_mesh, basis_x, CEED_VECTOR_ACTIVE);
676 CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
677 CeedOperatorSetField(op_setup, "qdata", elem_restr_q_data_points, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
678 CeedOperatorAtPointsSetPoints(op_setup, elem_restr_x_points, x_ref_points);
679
680 CeedOperatorApply(op_setup, x_coord, q_data_points, CEED_REQUEST_IMMEDIATE);
681
682 // Cleanup
683 CeedQFunctionDestroy(&qf_setup);
684 CeedOperatorDestroy(&op_setup);
685 }
686
687 // Cleanup
688 PetscCall(ISDestroy(&is_points));
689 PetscCall(PetscFree(offsets));
690 PetscCall(VecDestroy(&X_ref));
691 }
692
693 // Set up PDE operator
694
695 CeedQFunction qf_apply;
696 CeedOperator op_apply;
697 CeedInt in_scale = bp_data.in_mode == CEED_EVAL_GRAD ? dim : 1;
698 CeedInt out_scale = bp_data.out_mode == CEED_EVAL_GRAD ? dim : 1;
699
700 CeedQFunctionCreateInterior(ceed, 1, bp_data.apply, bp_data.apply_loc, &qf_apply);
701 CeedQFunctionAddInput(qf_apply, "u", num_comp * in_scale, bp_data.in_mode);
702 CeedQFunctionAddInput(qf_apply, "qdata", q_data_size, CEED_EVAL_NONE);
703 CeedQFunctionAddOutput(qf_apply, "v", num_comp * out_scale, bp_data.out_mode);
704
705 // Create the mass or diff operator
706 CeedOperatorCreateAtPoints(ceed, qf_apply, NULL, NULL, &op_apply);
707 CeedOperatorSetField(op_apply, "u", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE);
708 CeedOperatorSetField(op_apply, "qdata", elem_restr_q_data_points, CEED_BASIS_NONE, q_data_points);
709 CeedOperatorSetField(op_apply, "v", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE);
710 CeedOperatorAtPointsSetPoints(op_apply, elem_restr_x_points, x_ref_points);
711
712 // Set up RHS if needed
713 if (setup_rhs) {
714 CeedQFunction qf_setup_rhs;
715 CeedOperator op_setup_rhs;
716 Vec rhs_loc;
717 CeedVector rhs_ceed, target_ceed;
718 PetscMemType rhs_mem_type, target_mem_type;
719
720 // Create RHS vector
721 PetscCall(DMCreateLocalVector(dm_mesh, &rhs_loc));
722
723 CeedElemRestrictionCreateVector(elem_restr_u_points, &target_ceed, NULL);
724 CeedElemRestrictionCreateVector(elem_restr_u_mesh, &rhs_ceed, NULL);
725
726 // Create the q-function that sets up the RHS and true solution
727 CeedQFunctionCreateInterior(ceed, 1, bp_data.setup_rhs, bp_data.setup_rhs_loc, &qf_setup_rhs);
728 CeedQFunctionAddInput(qf_setup_rhs, "x", dim, CEED_EVAL_INTERP);
729 CeedQFunctionAddInput(qf_setup_rhs, "qdata", q_data_size, CEED_EVAL_NONE);
730 CeedQFunctionAddOutput(qf_setup_rhs, "true solution", num_comp, CEED_EVAL_NONE);
731 CeedQFunctionAddOutput(qf_setup_rhs, "rhs", num_comp, CEED_EVAL_INTERP);
732
733 // Create the operator that builds the RHS and true solution
734 CeedOperatorCreateAtPoints(ceed, qf_setup_rhs, NULL, NULL, &op_setup_rhs);
735 CeedOperatorSetField(op_setup_rhs, "x", elem_restr_x_mesh, basis_x, CEED_VECTOR_ACTIVE);
736 CeedOperatorSetField(op_setup_rhs, "qdata", elem_restr_q_data_points, CEED_BASIS_NONE, q_data_points);
737 CeedOperatorSetField(op_setup_rhs, "true solution", elem_restr_u_points, CEED_BASIS_NONE, target_ceed);
738 CeedOperatorSetField(op_setup_rhs, "rhs", elem_restr_u_mesh, basis_u, CEED_VECTOR_ACTIVE);
739 CeedOperatorAtPointsSetPoints(op_setup_rhs, elem_restr_x_points, x_ref_points);
740
741 // Set up the libCEED context
742 CeedQFunctionContext ctx_rhs_setup;
743 CeedQFunctionContextCreate(ceed, &ctx_rhs_setup);
744 CeedScalar rhs_setup_data[2] = {R, l};
745 CeedQFunctionContextSetData(ctx_rhs_setup, CEED_MEM_HOST, CEED_COPY_VALUES, sizeof rhs_setup_data, &rhs_setup_data);
746 CeedQFunctionSetContext(qf_setup_rhs, ctx_rhs_setup);
747 CeedQFunctionContextDestroy(&ctx_rhs_setup);
748
749 // Setup RHS and target
750 PetscCall(VecP2C(rhs_loc, &rhs_mem_type, rhs_ceed));
751 PetscCall(VecP2C(target, &target_mem_type, target_ceed));
752 CeedOperatorApply(op_setup_rhs, x_coord, rhs_ceed, CEED_REQUEST_IMMEDIATE);
753 PetscCall(VecC2P(rhs_ceed, rhs_mem_type, rhs_loc));
754 PetscCall(VecC2P(target_ceed, target_mem_type, target));
755
756 // Local-to-global
757 PetscCall(VecZeroEntries(rhs));
758 PetscCall(DMLocalToGlobal(dm_mesh, rhs_loc, ADD_VALUES, rhs));
759
760 PetscCall(VecViewFromOptions(rhs, NULL, "-rhs_view"));
761
762 // Cleanup
763 PetscCall(DMRestoreLocalVector(dm_mesh, &rhs_loc));
764 CeedVectorDestroy(&rhs_ceed);
765 CeedVectorDestroy(&target_ceed);
766 CeedQFunctionDestroy(&qf_setup_rhs);
767 CeedOperatorDestroy(&op_setup_rhs);
768 }
769
770 // Save libCEED data
771 data->basis_x = basis_x;
772 data->basis_u = basis_u;
773 data->elem_restr_x = elem_restr_x_mesh;
774 data->elem_restr_u = elem_restr_u_mesh;
775 data->elem_restr_u_i = elem_restr_u_points;
776 data->elem_restr_qd_i = elem_restr_q_data_points;
777 data->qf_apply = qf_apply;
778 data->op_apply = op_apply;
779 data->q_data = q_data_points;
780 data->q_data_size = q_data_size;
781
782 // Cleanup
783 PetscCall(VecReadC2P(x_coord, X_mem_type, X_loc));
784 CeedVectorDestroy(&x_coord);
785 CeedVectorDestroy(&x_ref_points);
786 CeedElemRestrictionDestroy(&elem_restr_x_points);
787
788 PetscFunctionReturn(PETSC_SUCCESS);
789 }
790