Coverage for dragonfly_sim/core/form_manager.py: 87%
93 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-17 03:37 +0000
1import numpy as np
2import ufl
3import dolfinx
5# Subdomain IDs of the tagged integration measures
6CELL_TAG = 1
7INTERIOR_FACET_TAG = 1
9def _all_entities_meshtags(mesh, dim, value):
10 # Apply the meshtag `value` to all local + ghost entities of dimension `dim`
11 imap = mesh.topology.index_map(dim)
12 if imap is None:
13 mesh.topology.create_entities(dim)
14 imap = mesh.topology.index_map(dim)
16 n_entities = imap.size_local + imap.num_ghosts
17 indices = np.arange(n_entities, dtype=np.int32)
18 values = np.full(n_entities, value, dtype=np.int32)
19 return dolfinx.mesh.meshtags(mesh, dim, indices, values)
22class FormManager:
23 """
24 Owns the integration measures and the form graph.
26 Instances come in two flavours, because the two things have different
27 lifetimes:
29 - MESH-LEVEL (`model=None`, reached as `mesh_obj.form_manager`): owns the
30 measures. There is one per Mesh, so every model built on that Mesh
31 integrates against the same measures.
33 Mesh itself exposes no dx/dS/ds. It did briefly, as delegating
34 properties, and that was worse than either alternative: a measure
35 reached through the Mesh reads as though the Mesh owned it, which is
36 the exact misreading that let the old construct_* methods overwrite
37 measures in place underneath compiled forms and live postprocessor
38 reads. Naming the owner at every use site is the point.
39 - MODEL-LEVEL (`model=...`, `measure_source=<the mesh-level one>`): owns
40 one model's residual/Jacobian expressions and its compiled forms, and
41 reads the measures through the mesh-level instance.
43 The measures are built incrementally, because the two construction calls
44 are independent and the existing call order varies: set_up_sim
45 builds the boundary measure first (it needs the facet split to tag the
46 wall before the mesh warper is constructed) and dx/dS later, while
47 some tests build dx/dS and never build ds at all.
48 """
50 def __init__(self, mesh_obj, model=None, measure_source=None):
51 self.mesh_obj = mesh_obj
52 self.mesh = mesh_obj.mesh
53 self.model = model
54 self._measure_source = measure_source
56 self._dx = None
57 self._dS = None
58 self._ds = None
60 self._cell_tags = None
61 self._interior_facet_tags = None
63 self.cell_metadata = None
64 self.boundary_metadata = None
66 # {key: ufl.Form} -- the UFL expressions registered on this manager.
67 self._ufl_forms = {}
69 # {key: dolfinx.fem.Form} -- the compiled bindings of those
70 # expressions, cached so each form is compiled at most once.
71 self._compiled = {}
73 # ------------------------------------------------------------------
74 # Measures
75 # ------------------------------------------------------------------
76 def build_cell_and_interior_facet_measures(self, measure_metadata=None):
77 """
78 dx and dS, tagged and backed by all-ones MeshTags.
80 Both come back already restricted to their tag, so a call site writes
81 `* forms.dx` and gets dx(CELL_TAG).
82 """
83 self.cell_metadata = measure_metadata
84 tdim = self.mesh.topology.dim
86 self._cell_tags = _all_entities_meshtags(self.mesh, tdim, CELL_TAG)
87 self._interior_facet_tags = _all_entities_meshtags(
88 self.mesh, tdim - 1, INTERIOR_FACET_TAG)
90 kwargs = {} if measure_metadata is None else {"metadata": measure_metadata}
91 self._dx = ufl.Measure(
92 "dx", domain=self.mesh, subdomain_data=self._cell_tags,
93 **kwargs)(CELL_TAG)
94 self._dS = ufl.Measure(
95 "dS", domain=self.mesh, subdomain_data=self._interior_facet_tags,
96 **kwargs)(INTERIOR_FACET_TAG)
98 def build_boundary_measure(self, measure_metadata=None):
99 """
100 ds, over the boundary facet tagging already on the Mesh.
102 The tagging itself stays Mesh.tag_boundary_facets' job and is NOT
103 built here. It is mesh topology, not a measure: the mesh warper, the
104 FFD wall-node selection and the meshtags XDMF export all read
105 Mesh.meshtags without ever touching a form. Only the measure derived
106 from it belongs to this class.
108 Unlike dx/dS this measure comes back UNrestricted, because every call
109 site applies its own tag: ds(WALL_TAG), ds(FARFIELD_TAG), ...
110 """
111 if self.mesh_obj.meshtags is None:
112 raise RuntimeError(
113 "The boundary facets have not been tagged yet. Call "
114 "Mesh.tag_boundary_facets() before build_boundary_measure().")
116 self.boundary_metadata = measure_metadata
118 kwargs = {} if measure_metadata is None else {"metadata": measure_metadata}
119 self._ds = ufl.Measure(
120 "ds", domain=self.mesh, subdomain_data=self.mesh_obj.meshtags,
121 **kwargs)
123 return self._ds
125 def _measures(self):
126 """The manager that actually holds the measures (self, or the mesh-level one)."""
127 return self if self._measure_source is None else self._measure_source
129 @property
130 def dx(self):
131 src = self._measures()
132 if src is not self:
133 return src.dx
134 if self._dx is None:
135 raise RuntimeError(
136 "dx has not been built yet. Call "
137 "form_manager.build_cell_and_interior_facet_measures() "
138 "before building any form.")
139 return self._dx
141 @property
142 def dS(self):
143 src = self._measures()
144 if src is not self:
145 return src.dS
146 if self._dS is None:
147 raise RuntimeError(
148 "dS has not been built yet. Call "
149 "form_manager.build_cell_and_interior_facet_measures() "
150 "before building any form.")
151 return self._dS
153 @property
154 def ds(self):
155 src = self._measures()
156 if src is not self:
157 return src.ds
158 if self._ds is None:
159 raise RuntimeError(
160 "ds has not been built yet. Call Mesh.tag_boundary_facets() "
161 "and then form_manager.build_boundary_measure() before "
162 "building any boundary form.")
163 return self._ds
165 # ------------------------------------------------------------------
166 # The form graph
167 # ------------------------------------------------------------------
168 def build_residual(self, U_coeff=None):
169 """
170 The physical residual: interior + boundary flux terms.
172 `U_coeff` is accepted for backward compatibility but has no effect.
174 The physics itself stays on CompressibleEulerModel; this composes it.
175 """
176 model = self.model
177 if model is None:
178 raise RuntimeError(
179 "build_residual() needs a model-level FormManager; this one "
180 "owns measures only.")
182 F = model.compute_interior_integral_terms()
183 for F_bc in model.bcs:
184 F += F_bc
186 return F
188 def ufl_form(self, key, factory=None):
189 """The UFL expression registered under `key`, building it on demand."""
190 if key not in self._ufl_forms:
191 if factory is None:
192 raise KeyError(
193 "No UFL form registered under {!r} and no factory given "
194 "to build one.".format(key))
195 self._ufl_forms[key] = factory()
196 return self._ufl_forms[key]
198 def compiled(self, key, factory=None):
199 """
200 The compiled form (dolfinx.fem.form) of the expression registered
201 under `key`, compiled on first request and cached.
202 """
203 if key not in self._compiled:
204 self._compiled[key] = dolfinx.fem.form(self.ufl_form(key, factory))
205 return self._compiled[key]
207 def assemble_matrix_scipy(self, key, factory=None, csr=False):
208 """
209 The compiled form registered under `key`, assembled as a scipy sparse
210 matrix.
212 NOTE the format. MatrixCSR.to_scipy() returns BSR, not CSR, whenever
213 the space is blocked -- which the Euler state space always is
214 (block size 2 + dimensions). scipy's BSR does not implement row
215 indexing at all (it raises NotImplementedError), so pass csr=True
216 for anything that slices rows. The default keeps the block structure,
217 which is the cheaper representation when you only need matvecs.
218 """
219 A = dolfinx.fem.assemble_matrix(self.compiled(key, factory))
220 A.scatter_reverse()
221 sparse = A.to_scipy()
222 return sparse.tocsr() if csr else sparse