Element Sets in 3-D

An element set is a group of mesh elements that you want to solve with one model, for example all the elements of one material. SELF can split a 3-D mesh into element sets, run a separate discontinuous Galerkin model on each set, and couple the sets where they share faces. The sets can use different model types, so one set can carry an acoustic model while another carries an elastic one. The parent model that does this is MultiSetModel3D, and it extends Model, so the usual time integrators and ForwardStep drive it unchanged.

Partitions and sub-meshes

You describe a partition by giving each element a set id \(s(e) \in \{1, \dots, n_{sets}\}\). Within set \(k\), the elements are renumbered in global order,

\[ g_k(e) = \#\{ e' \le e : s(e') = k \}, \]

so each rank's elements of a set form one contiguous block, and no element changes rank. SELF_ElementSetMesh_3D builds one Mesh3D per set from this numbering. A face between two elements of the same set stays an interior face of that set's sub-mesh. A face between elements of different sets becomes an interface face in both sub-meshes, tagged with the boundary condition id SELF_BC_SET_INTERFACE. An interface face is never a physical boundary: there is always an element of another set on its far side, and the parent model supplies its exterior state. Domain boundary faces keep their boundary condition ids.

A 2:1 or anisotropic mortar must lie entirely inside one set. ValidateSetPartition rejects a partition that splits a mortar, a set id outside \(1 \dots n_{sets}\), and a set with no elements.

Setting up a multi-set model

Model types cannot be constructed generically in Fortran, so setup takes three steps. First, Init validates the partition and builds every sub-mesh and its geometry. Second, you allocate and initialize a model on every set. Third, SetupInterfaces checks the set models and builds the interface map.

type(MultiSetModel3D) :: ms
integer :: k

call ms%Init(mesh,interp,setOfElem,2)

do k = 1,2
  allocate(LinearEuler3D :: ms%sets(k)%model)
  call ms%sets(k)%model%Init(ms%sets(k)%mesh,ms%sets(k)%geometry)
  ! set the initial condition of set k here
enddo

call ms%SetupInterfaces()
call ms%SetTimeIntegrator('rk3')
call ms%ForwardStep(tn=endtime,dt=dt,ioInterval=iointerval)

Here setOfElem(1:mesh%nElem) holds the set id of each element this rank owns. The second step must run on every rank and for every set, including a set that has no elements on that rank, because model Init is collective over the set's communicator. Each model must be initialized on its own set's mesh and geometry, which SetupInterfaces checks. The sets array must not be reassigned or reallocated afterwards, since every model holds pointers into its own set.

What happens in a stage

Every time the integrator asks for a tendency, the parent

1. calls its own PreTendencyHook,
2. interpolates the solution of each set with interface faces on this rank to its element faces,
3. exchanges the face traces across every interface face, so that each face holds its partner's trace in iface%ghost,
4. calls FillInterfaceBoundary(k) for each set, which writes solution%extBoundary on that set's interface faces,
5. calls each set model's CalculateTendency, which runs that set's own PreTendencyHook.

The exchange applies the same face orientation flip as the conforming side exchange, so on an interface face the exterior state is exactly the value the unsplit mesh would have used. Splitting a single model into sets of the same type therefore reproduces the unsplit solution to the last bit, which is what the bitwise tests check.

The exchange runs on its own duplicate of the mesh communicator, with one message per neighbouring rank and a constant tag. Each sub-mesh decomposition also gets its own duplicate communicator, so the halo exchanges of different sets cannot match each other's messages. Free releases all of these.

SetupInterfaces registers a boundary condition for SELF_BC_SET_INTERFACE on every set model that writes nothing, and maps the boundary conditions again. This keeps interface faces out of the unmapped boundary count without anything but the parent writing their exterior state.

Coupling sets with different variables

The default FillInterfaceBoundary copies the partner's trace into the exterior state unchanged, which is only meaningful when both sets carry the same variables. It stops with an error when the two sets have a different nvar. To couple different models, extend MultiSetModel3D and override FillInterfaceBoundary(iset). For each interface face \(f\) of set iset, the face is (iface%faceElem(f), iface%faceSide(f)) in that set's sub-mesh, the partner set is iface%facePartnerSet(f), and iface%ghost(:,:,1:nvar_p,f) holds the partner's trace in the receiving face's orientation. The faces of set iset are iface%setFaceStart(iset) to iface%setFaceStart(iset+1)-1.

On the GPU backend the ghost traces are also on the device in iface%ghost_gpu, and an override must write solution%extBoundary_gpu. Do not call solution%UpdateDevice() inside the step loop: it resets the static variable state of the halo exchange.

Running in parallel

Element ownership never changes, so a rank may own no elements of a set. Such a rank still takes part in everything that is collective for that set: model Init, the boundary condition mapping, the entropy reduction, and the HDF5 output. It skips that set's CalculateTendency and solution updates. For this to be safe, a set model must not call a collective, or post messages on its mesh communicator, inside CalculateTendency or anything it calls, such as PreTendencyHook or a boundary condition handler. Every rank processes the sets in the same order.

Output

WriteModel writes one HDF5 file holding the shared interpolant, the set id of every element in /elementsets/setOfElem, and for each set \(k\) its solution in /set_<k>/controlgrid/solution and node coordinates in /set_<k>/controlgrid/geometry, with <k> written as four digits and elements in set order. ReadModel reads the solutions back into a model built on the same partition and polynomial degree, and stops if the stored partition or the number of nodes per element differs. WriteTecplot writes one rank-local file per set, named set<k>.solution.<rank>.<iter>.tec.

Limitations

Element sets are available in 3-D only. Solution gradients are not supported, so gradient_enabled must be false on every set; the gradient exchange inside each set's tendency has no hook for interface faces. Adaptive mesh refinement, mortars that cross sets, and rebalancing elements between ranks per set are not supported. Sub-meshes cannot be re-extracted after setup. Bitwise equality relies on the split and unsplit runs performing the same floating point operations in the same order; a model whose FillInterfaceBoundary or flux differs between sets will not reproduce the unsplit run to the last bit.

Tests

elementset3d_extract splits a two-material mesh and checks every sub-mesh face against the global mesh. elementset3d_interface exchanges face coordinates through the interface map and checks that each face receives its own coordinates back, which validates the face pairing and the flips. elementset3d_interface_unequal and elementset3d_interface_unequal_mpi exchange traces between a three-variable and a five-variable set, on a mesh with flipped faces and across two ranks, and check every ghost value against a known function of position. multiset_lineareuler3d_bitwise and its two-rank twin multiset_lineareuler3d_bitwise_mpi run linear Euler as one model and as one set, an x-split and a checkerboard, and require identical solutions. multiset_zeroelem_mpi puts one set entirely on rank 0, runs it through ForwardStep, and checks the solution, the entropy and an HDF5 write and read. multiset3d_guard_read_degree checks that a pickup file written at one polynomial degree is refused by a model at another. elementset3d_interface_flips exchanges traces across faces carrying every one of the eight side flips, on 24 pairs of cubes in every relative orientation. multiset_lineareuler3d_integrators repeats the x-split comparison with the euler, rk2 and rk4 integrators and checks the Tecplot output. multiset_lineareuler3d_mortar_mpi runs one set on a 2:1 mortar mesh and on an anisotropic mortar mesh across two ranks and requires identical solutions. multiset3d_read_missing_variable reads a pickup file that lacks one variable of one set and checks that the variable keeps its value.