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/petscutils.h"
9
10 // -----------------------------------------------------------------------------
11 // Convert PETSc MemType to libCEED MemType
12 // -----------------------------------------------------------------------------
MemTypeP2C(PetscMemType mem_type)13 CeedMemType MemTypeP2C(PetscMemType mem_type) { return PetscMemTypeDevice(mem_type) ? CEED_MEM_DEVICE : CEED_MEM_HOST; }
14
15 // ------------------------------------------------------------------------------------------------
16 // PETSc-libCEED memory space utilities
17 // ------------------------------------------------------------------------------------------------
VecP2C(Vec X_petsc,PetscMemType * mem_type,CeedVector x_ceed)18 PetscErrorCode VecP2C(Vec X_petsc, PetscMemType *mem_type, CeedVector x_ceed) {
19 PetscScalar *x;
20
21 PetscFunctionBeginUser;
22 PetscCall(VecGetArrayAndMemType(X_petsc, &x, mem_type));
23 CeedVectorSetArray(x_ceed, MemTypeP2C(*mem_type), CEED_USE_POINTER, x);
24 PetscFunctionReturn(PETSC_SUCCESS);
25 }
26
VecC2P(CeedVector x_ceed,PetscMemType mem_type,Vec X_petsc)27 PetscErrorCode VecC2P(CeedVector x_ceed, PetscMemType mem_type, Vec X_petsc) {
28 PetscScalar *x;
29
30 PetscFunctionBeginUser;
31 CeedVectorTakeArray(x_ceed, MemTypeP2C(mem_type), &x);
32 PetscCall(VecRestoreArrayAndMemType(X_petsc, &x));
33 PetscFunctionReturn(PETSC_SUCCESS);
34 }
35
VecReadP2C(Vec X_petsc,PetscMemType * mem_type,CeedVector x_ceed)36 PetscErrorCode VecReadP2C(Vec X_petsc, PetscMemType *mem_type, CeedVector x_ceed) {
37 PetscScalar *x;
38
39 PetscFunctionBeginUser;
40 PetscCall(VecGetArrayReadAndMemType(X_petsc, (const PetscScalar **)&x, mem_type));
41 CeedVectorSetArray(x_ceed, MemTypeP2C(*mem_type), CEED_USE_POINTER, x);
42 PetscFunctionReturn(PETSC_SUCCESS);
43 }
44
VecReadC2P(CeedVector x_ceed,PetscMemType mem_type,Vec X_petsc)45 PetscErrorCode VecReadC2P(CeedVector x_ceed, PetscMemType mem_type, Vec X_petsc) {
46 PetscScalar *x;
47
48 PetscFunctionBeginUser;
49 CeedVectorTakeArray(x_ceed, MemTypeP2C(mem_type), &x);
50 PetscCall(VecRestoreArrayReadAndMemType(X_petsc, (const PetscScalar **)&x));
51 PetscFunctionReturn(PETSC_SUCCESS);
52 }
53
54 // -----------------------------------------------------------------------------
55 // Apply 3D Kershaw mesh transformation
56 // -----------------------------------------------------------------------------
57 // Transition from a value of "a" for x=0, to a value of "b" for x=1. Optionally
58 // smooth -- see the commented versions at the end.
step(const double a,const double b,double x)59 static double step(const double a, const double b, double x) {
60 if (x <= 0) return a;
61 if (x >= 1) return b;
62 return a + (b - a) * (x);
63 }
64
65 // 1D transformation at the right boundary
right(const double eps,const double x)66 static double right(const double eps, const double x) { return (x <= 0.5) ? (2 - eps) * x : 1 + eps * (x - 1); }
67
68 // 1D transformation at the left boundary
left(const double eps,const double x)69 static double left(const double eps, const double x) { return 1 - right(eps, 1 - x); }
70
71 // Apply 3D Kershaw mesh transformation
72 // The eps parameters are in (0, 1]
73 // Uniform mesh is recovered for eps=1
Kershaw(DM dm_orig,PetscScalar eps)74 PetscErrorCode Kershaw(DM dm_orig, PetscScalar eps) {
75 Vec coord;
76 PetscInt ncoord;
77 PetscScalar *c;
78
79 PetscFunctionBeginUser;
80 PetscCall(DMGetCoordinatesLocal(dm_orig, &coord));
81 PetscCall(VecGetLocalSize(coord, &ncoord));
82 PetscCall(VecGetArray(coord, &c));
83
84 for (PetscInt i = 0; i < ncoord; i += 3) {
85 PetscScalar x = c[i], y = c[i + 1], z = c[i + 2];
86 PetscInt layer = x * 6;
87 PetscScalar lambda = (x - layer / 6.0) * 6;
88 c[i] = x;
89
90 switch (layer) {
91 case 0:
92 c[i + 1] = left(eps, y);
93 c[i + 2] = left(eps, z);
94 break;
95 case 1:
96 case 4:
97 c[i + 1] = step(left(eps, y), right(eps, y), lambda);
98 c[i + 2] = step(left(eps, z), right(eps, z), lambda);
99 break;
100 case 2:
101 c[i + 1] = step(right(eps, y), left(eps, y), lambda / 2);
102 c[i + 2] = step(right(eps, z), left(eps, z), lambda / 2);
103 break;
104 case 3:
105 c[i + 1] = step(right(eps, y), left(eps, y), (1 + lambda) / 2);
106 c[i + 2] = step(right(eps, z), left(eps, z), (1 + lambda) / 2);
107 break;
108 default:
109 c[i + 1] = right(eps, y);
110 c[i + 2] = right(eps, z);
111 }
112 }
113 PetscCall(VecRestoreArray(coord, &c));
114 PetscFunctionReturn(PETSC_SUCCESS);
115 }
116
117 // -----------------------------------------------------------------------------
118 // Create BC label
119 // -----------------------------------------------------------------------------
CreateBCLabel(DM dm,const char name[])120 static PetscErrorCode CreateBCLabel(DM dm, const char name[]) {
121 DMLabel label;
122
123 PetscFunctionBeginUser;
124 PetscCall(DMCreateLabel(dm, name));
125 PetscCall(DMGetLabel(dm, name, &label));
126 PetscCall(DMPlexMarkBoundaryFaces(dm, PETSC_DETERMINE, label));
127 PetscCall(DMPlexLabelComplete(dm, label));
128 PetscFunctionReturn(PETSC_SUCCESS);
129 }
130
131 // -----------------------------------------------------------------------------
132 // This function sets up a DM for a given degree
133 // -----------------------------------------------------------------------------
SetupDMByDegree(DM dm,PetscInt p_degree,PetscInt q_extra,PetscInt num_comp_u,PetscInt dim,bool enforce_bc)134 PetscErrorCode SetupDMByDegree(DM dm, PetscInt p_degree, PetscInt q_extra, PetscInt num_comp_u, PetscInt dim, bool enforce_bc) {
135 PetscInt marker_ids[1] = {1};
136 PetscInt q_degree = p_degree + q_extra;
137 PetscFE fe;
138 MPI_Comm comm;
139 PetscBool is_simplex = PETSC_TRUE;
140
141 PetscFunctionBeginUser;
142 // Check if simplex or tensor-product mesh
143 PetscCall(DMPlexIsSimplex(dm, &is_simplex));
144 // Setup FE
145 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
146 PetscCall(PetscFECreateLagrange(comm, dim, num_comp_u, is_simplex, p_degree, q_degree, &fe));
147 PetscCall(DMAddField(dm, NULL, (PetscObject)fe));
148 PetscCall(DMCreateDS(dm));
149
150 {
151 // create FE field for coordinates
152 PetscFE fe_coords;
153 PetscInt num_comp_coord;
154 PetscCall(DMGetCoordinateDim(dm, &num_comp_coord));
155 PetscCall(PetscFECreateLagrange(comm, dim, num_comp_coord, is_simplex, 1, q_degree, &fe_coords));
156 PetscCall(DMSetCoordinateDisc(dm, fe_coords, PETSC_TRUE));
157 PetscCall(PetscFEDestroy(&fe_coords));
158 }
159
160 // Setup Dirichlet BC
161 // Note bp1, bp2 are projection and we don't need to apply BC
162 // For bp3,bp4, the target function is zero on the boundaries
163 // So we pass bcFunc = NULL in DMAddBoundary function
164 if (enforce_bc) {
165 PetscBool has_label;
166 PetscCall(DMHasLabel(dm, "marker", &has_label));
167 if (!has_label) {
168 PetscCall(CreateBCLabel(dm, "marker"));
169 }
170 DMLabel label;
171 PetscCall(DMGetLabel(dm, "marker", &label));
172 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, marker_ids, 0, 0, NULL, NULL, NULL, NULL, NULL));
173 PetscCall(DMSetOptionsPrefix(dm, "final_"));
174 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
175 }
176
177 if (!is_simplex) {
178 DM dm_coord;
179 PetscCall(DMGetCoordinateDM(dm, &dm_coord));
180 PetscCall(DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL));
181 PetscCall(DMPlexSetClosurePermutationTensor(dm_coord, PETSC_DETERMINE, NULL));
182 }
183 PetscCall(PetscFEDestroy(&fe));
184 PetscFunctionReturn(PETSC_SUCCESS);
185 }
186
187 // -----------------------------------------------------------------------------
188 // Get CEED restriction data from DMPlex
189 // -----------------------------------------------------------------------------
CreateRestrictionFromPlex(Ceed ceed,DM dm,CeedInt height,DMLabel domain_label,CeedInt value,CeedElemRestriction * elem_restr)190 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr) {
191 PetscInt num_elem, elem_size, num_dof, num_comp, *elem_restr_offsets_petsc;
192 CeedInt *elem_restr_offsets_ceed;
193
194 PetscFunctionBeginUser;
195 PetscCall(DMPlexGetLocalOffsets(dm, domain_label, value, height, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets_petsc));
196
197 PetscCall(IntArrayPetscToCeed(num_elem * elem_size, &elem_restr_offsets_petsc, &elem_restr_offsets_ceed));
198 CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp, 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, elem_restr_offsets_ceed, elem_restr);
199 PetscCall(PetscFree(elem_restr_offsets_ceed));
200 PetscFunctionReturn(PETSC_SUCCESS);
201 }
202
203 // -----------------------------------------------------------------------------
204 // Utility function - convert from DMPolytopeType to CeedElemTopology
205 // -----------------------------------------------------------------------------
ElemTopologyP2C(DMPolytopeType cell_type)206 CeedElemTopology ElemTopologyP2C(DMPolytopeType cell_type) {
207 switch (cell_type) {
208 case DM_POLYTOPE_TRIANGLE:
209 return CEED_TOPOLOGY_TRIANGLE;
210 case DM_POLYTOPE_QUADRILATERAL:
211 return CEED_TOPOLOGY_QUAD;
212 case DM_POLYTOPE_TETRAHEDRON:
213 return CEED_TOPOLOGY_TET;
214 case DM_POLYTOPE_HEXAHEDRON:
215 return CEED_TOPOLOGY_HEX;
216 default:
217 return 0;
218 }
219 }
220
221 // -----------------------------------------------------------------------------
222 // Convert DM field to DS field
223 // -----------------------------------------------------------------------------
DMFieldToDSField(DM dm,DMLabel domain_label,PetscInt dm_field,PetscInt * ds_field)224 PetscErrorCode DMFieldToDSField(DM dm, DMLabel domain_label, PetscInt dm_field, PetscInt *ds_field) {
225 PetscDS ds;
226 IS field_is;
227 const PetscInt *fields;
228 PetscInt num_fields;
229
230 PetscFunctionBeginUser;
231 // Translate dm_field to ds_field
232 PetscCall(DMGetRegionDS(dm, domain_label, &field_is, &ds, NULL));
233 PetscCall(ISGetIndices(field_is, &fields));
234 PetscCall(ISGetSize(field_is, &num_fields));
235 for (PetscInt i = 0; i < num_fields; i++) {
236 if (dm_field == fields[i]) {
237 *ds_field = i;
238 break;
239 }
240 }
241 PetscCall(ISRestoreIndices(field_is, &fields));
242
243 if (*ds_field == -1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Could not find dm_field %" PetscInt_FMT " in DS", dm_field);
244 PetscFunctionReturn(PETSC_SUCCESS);
245 }
246
247 // -----------------------------------------------------------------------------
248 // Create libCEED Basis from PetscTabulation
249 // -----------------------------------------------------------------------------
BasisCreateFromTabulation(Ceed ceed,DM dm,DMLabel domain_label,PetscInt label_value,PetscInt height,PetscInt face,PetscFE fe,PetscTabulation basis_tabulation,PetscQuadrature quadrature,CeedBasis * basis)250 PetscErrorCode BasisCreateFromTabulation(Ceed ceed, DM dm, DMLabel domain_label, PetscInt label_value, PetscInt height, PetscInt face, PetscFE fe,
251 PetscTabulation basis_tabulation, PetscQuadrature quadrature, CeedBasis *basis) {
252 PetscInt first_point;
253 PetscInt ids[1] = {label_value};
254 DMLabel depth_label;
255 DMPolytopeType cell_type;
256 CeedElemTopology elem_topo;
257 PetscScalar *q_points, *interp, *grad;
258 const PetscScalar *q_weights;
259 PetscDualSpace dual_space;
260 PetscInt num_dual_basis_vectors;
261 PetscInt dim, num_comp, P, Q;
262
263 PetscFunctionBeginUser;
264 // General basis information
265 PetscCall(PetscFEGetSpatialDimension(fe, &dim));
266 PetscCall(PetscFEGetNumComponents(fe, &num_comp));
267 PetscCall(PetscFEGetDualSpace(fe, &dual_space));
268 PetscCall(PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors));
269 P = num_dual_basis_vectors / num_comp;
270
271 // Use depth label if no domain label present
272 if (!domain_label) {
273 PetscInt depth;
274
275 PetscCall(DMPlexGetDepth(dm, &depth));
276 PetscCall(DMPlexGetDepthLabel(dm, &depth_label));
277 ids[0] = depth - height;
278 }
279
280 // Get cell interp, grad, and quadrature data
281 PetscCall(DMGetFirstLabeledPoint(dm, dm, domain_label ? domain_label : depth_label, 1, ids, height, &first_point, NULL));
282 PetscCall(DMPlexGetCellType(dm, first_point, &cell_type));
283 elem_topo = ElemTopologyP2C(cell_type);
284 if (!elem_topo) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DMPlex topology not supported");
285 {
286 size_t q_points_size;
287 const PetscScalar *q_points_petsc;
288 PetscInt q_dim;
289
290 PetscCall(PetscQuadratureGetData(quadrature, &q_dim, NULL, &Q, &q_points_petsc, &q_weights));
291 q_points_size = Q * dim * sizeof(CeedScalar);
292 PetscCall(PetscCalloc(q_points_size, &q_points));
293 for (PetscInt q = 0; q < Q; q++) {
294 for (PetscInt d = 0; d < q_dim; d++) q_points[q * dim + d] = q_points_petsc[q * q_dim + d];
295 }
296 }
297
298 // Convert to libCEED orientation
299 {
300 PetscBool is_simplex = PETSC_FALSE;
301 IS permutation = NULL;
302 const PetscInt *permutation_indices;
303
304 PetscCall(DMPlexIsSimplex(dm, &is_simplex));
305 if (!is_simplex) {
306 PetscSection section;
307
308 // -- Get permutation
309 PetscCall(DMGetLocalSection(dm, §ion));
310 PetscCall(PetscSectionGetClosurePermutation(section, (PetscObject)dm, dim, num_comp * P, &permutation));
311 PetscCall(ISGetIndices(permutation, &permutation_indices));
312 }
313
314 // -- Copy interp, grad matrices
315 PetscCall(PetscCalloc(P * Q * sizeof(CeedScalar), &interp));
316 PetscCall(PetscCalloc(P * Q * dim * sizeof(CeedScalar), &grad));
317 const CeedInt c = 0;
318 for (CeedInt q = 0; q < Q; q++) {
319 for (CeedInt p_ceed = 0; p_ceed < P; p_ceed++) {
320 CeedInt p_petsc = is_simplex ? (p_ceed * num_comp) : permutation_indices[p_ceed * num_comp];
321
322 interp[q * P + p_ceed] = basis_tabulation->T[0][((face * Q + q) * P * num_comp + p_petsc) * num_comp + c];
323 for (CeedInt d = 0; d < dim; d++) {
324 grad[(d * Q + q) * P + p_ceed] = basis_tabulation->T[1][(((face * Q + q) * P * num_comp + p_petsc) * num_comp + c) * dim + d];
325 }
326 }
327 }
328
329 // -- Cleanup
330 if (permutation) PetscCall(ISRestoreIndices(permutation, &permutation_indices));
331 PetscCall(ISDestroy(&permutation));
332 }
333
334 // Finally, create libCEED basis
335 CeedBasisCreateH1(ceed, elem_topo, num_comp, P, Q, interp, grad, q_points, q_weights, basis);
336 PetscCall(PetscFree(q_points));
337 PetscCall(PetscFree(interp));
338 PetscCall(PetscFree(grad));
339 PetscFunctionReturn(PETSC_SUCCESS);
340 }
341
342 // -----------------------------------------------------------------------------
343 // Get CEED Basis from DMPlex
344 // -----------------------------------------------------------------------------
CreateBasisFromPlex(Ceed ceed,DM dm,DMLabel domain_label,CeedInt label_value,CeedInt height,CeedInt dm_field,BPData bp_data,CeedBasis * basis)345 PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedInt label_value, CeedInt height, CeedInt dm_field, BPData bp_data,
346 CeedBasis *basis) {
347 PetscDS ds;
348 PetscFE fe;
349 PetscQuadrature quadrature;
350 PetscBool is_simplex = PETSC_TRUE;
351 PetscInt ds_field = -1;
352
353 PetscFunctionBeginUser;
354 // Get element information
355 PetscCall(DMGetRegionDS(dm, domain_label, NULL, &ds, NULL));
356 PetscCall(DMFieldToDSField(dm, domain_label, dm_field, &ds_field));
357 PetscCall(PetscDSGetDiscretization(ds, ds_field, (PetscObject *)&fe));
358 PetscCall(PetscFEGetHeightSubspace(fe, height, &fe));
359 PetscCall(PetscFEGetQuadrature(fe, &quadrature));
360
361 // Check if simplex or tensor-product mesh
362 PetscCall(DMPlexIsSimplex(dm, &is_simplex));
363
364 // Build libCEED basis
365 if (is_simplex) {
366 PetscTabulation basis_tabulation;
367 PetscInt num_derivatives = 1, face = 0;
368
369 PetscCall(PetscFEGetCellTabulation(fe, num_derivatives, &basis_tabulation));
370 PetscCall(BasisCreateFromTabulation(ceed, dm, domain_label, label_value, height, face, fe, basis_tabulation, quadrature, basis));
371 } else {
372 PetscDualSpace dual_space;
373 PetscInt num_dual_basis_vectors;
374 PetscInt dim, num_comp, P, Q;
375
376 PetscCall(PetscFEGetSpatialDimension(fe, &dim));
377 PetscCall(PetscFEGetNumComponents(fe, &num_comp));
378 PetscCall(PetscFEGetDualSpace(fe, &dual_space));
379 PetscCall(PetscDualSpaceGetDimension(dual_space, &num_dual_basis_vectors));
380 P = num_dual_basis_vectors / num_comp;
381 PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &Q, NULL, NULL));
382
383 CeedInt P_1d = (CeedInt)round(pow(P, 1.0 / dim));
384 CeedInt Q_1d = (CeedInt)round(pow(Q, 1.0 / dim));
385
386 CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, P_1d, Q_1d, bp_data.q_mode, basis);
387 }
388 PetscFunctionReturn(PETSC_SUCCESS);
389 }
390
391 // -----------------------------------------------------------------------------
392 // Utilities
393 // -----------------------------------------------------------------------------
394
395 // Utility function, compute three factors of an integer
Split3(PetscInt size,PetscInt m[3],bool reverse)396 static void Split3(PetscInt size, PetscInt m[3], bool reverse) {
397 for (PetscInt d = 0, size_left = size; d < 3; d++) {
398 PetscInt try = (PetscInt)PetscCeilReal(PetscPowReal(size_left, 1. / (3 - d)));
399 while (try * (size_left / try) != size_left) try++;
400 m[reverse ? 2 - d : d] = try;
401 size_left /= try;
402 }
403 }
404
Max3(const PetscInt a[3])405 static int Max3(const PetscInt a[3]) { return PetscMax(a[0], PetscMax(a[1], a[2])); }
406
Min3(const PetscInt a[3])407 static int Min3(const PetscInt a[3]) { return PetscMin(a[0], PetscMin(a[1], a[2])); }
408
409 // -----------------------------------------------------------------------------
410 // Create distribute dm
411 // -----------------------------------------------------------------------------
CreateDistributedDM(RunParams rp,DM * dm)412 PetscErrorCode CreateDistributedDM(RunParams rp, DM *dm) {
413 PetscFunctionBeginUser;
414 // Setup DM
415 if (rp->read_mesh) {
416 PetscCall(DMPlexCreateFromFile(PETSC_COMM_WORLD, rp->filename, NULL, PETSC_TRUE, dm));
417 } else {
418 if (rp->user_l_nodes) {
419 // Find a nicely composite number of elements no less than global nodes
420 PetscMPIInt size;
421 PetscCall(MPI_Comm_size(rp->comm, &size));
422 for (PetscInt g_elem = PetscMax(1, size * rp->local_nodes / PetscPowInt(rp->degree, rp->dim));; g_elem++) {
423 Split3(g_elem, rp->mesh_elem, true);
424 if (Max3(rp->mesh_elem) / Min3(rp->mesh_elem) <= 2) break;
425 }
426 }
427
428 PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, rp->dim, rp->simplex, rp->mesh_elem, NULL, NULL, NULL, PETSC_TRUE, 0, PETSC_FALSE, dm));
429 }
430
431 PetscCall(DMSetFromOptions(*dm));
432 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
433 PetscFunctionReturn(PETSC_SUCCESS);
434 }
435
436 // -----------------------------------------------------------------------------
437