Coverage for dragonfly_sim/core/farfield.py: 96%
154 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-17 03:37 +0000
1"""
2Non-reflecting far-field refinements on top of the characteristic far field
3(CompressibleEulerModel.define_farfield_bc).
5TransverseFarfield ('riemann2') integrates the incoming Riemann invariant's
6transverse terms along the boundary, which turns the plain characteristic
7condition (exact only at normal incidence) into a second-order absorbing one
8(Engquist-Majda / Higdon) for oblique waves. Unsteady runs only.
9"""
10import numpy as np
11from mpi4py import MPI
12import ufl
13import basix
14import dolfinx
15import dolfinx.fem.petsc
17from dragonfly_sim.utils.Euler_utils import pressure
20class TransverseFarfield:
21 """
22 Riemann far-field with transverse terms in the incoming characteristic.
24 At the far field the incoming invariant R- = u.n - 2c/(g-1) obeys
25 (straight boundary, tangent tau = (-n_y, n_x))
27 dR-/dt + (u.n - c) dR-/dn + u_t dR-/dtau - c du_t/dtau = 0.
29 Dropping only the normal (incoming-wave) term and keeping beta times the
30 transverse ones, plus a weak relaxation to the freestream, gives
32 dR-/dt = -beta (u_t dR-/dtau - c du_t/dtau) - K (R- - R-_inf).
34 beta = 0 with R-(0) = R-_inf is exactly riemann_farfield_state. For
35 linear acoustics on a boundary with tangential mean flow the plane-wave
36 reflection coefficient is
38 r(theta) = [a(s-1) + beta r^2] / [a(s+1) - beta r^2],
39 s = cos(theta), r = sin(theta), a = 1 + (1-beta) M r,
41 so beta = 0 gives (s-1)/(s+1), beta = 1 gives (1-s)/(1+s) (no better),
42 and beta = 1/2 gives -((1-s)/(1+s))^2 for any tangential Mach number:
43 the second-order Engquist-Majda/Higdon absorbing condition.
44 K = sigma c_inf (1 - M^2) / L (Poinsot & Lele) only stops long-term drift.
46 The equation is integrated with the solver's own BDF coefficients for the
47 deviation delta = R- - R-_inf (a0 + a1 + a2 = 0 cancels R-_inf), with
48 the transverse terms extrapolated from the two previous levels
49 (IMEX/SBDF2; on a BDF1 step T* = T(U^n)):
51 delta^{n+1} = (-a1 delta^n - a2 delta^{n-1} - dt beta T*) / (a0 + K dt),
52 T* = 2 T(U^n) - T(U^{n-1}).
54 Explicit is safe here: the far-field cells are coarse (h ~ 1.35), so the
55 transverse transport has CFL ~ U dt/h ~ c dt/h ~ 0.04. Keeping T implicit
56 put gradients of the nonlinear primitives inside the Newton Jacobian,
57 whose form took longer to compile than dolfinx's JIT wait on the other
58 ranks; this way the residual only sees delta^{n+1} as a stored field and
59 its Jacobian is the plain Riemann one.
61 T is taken from the interior trace. The levels live in DG fields of the
62 solution degree whose traces on the far-field facets are the boundary
63 data: delta^{n+1} is evaluated at 3 Gauss points per facet and
64 L2-projected (cell-wise least squares; a corner cell's two facets share
65 its corner dof, so there the trace is a compromise). Storing delta rather
66 than R- keeps the freestream exact: delta stays identically zero for
67 beta = 0.
69 Usage: construct it, pass it to model.define_farfield_bc(tag,
70 transverse=...), and register it with the model's TimeIntegrator
71 (attach(integrator)), which then calls prepare() before every physical
72 step's Newton solve and commit() after the step is accepted. Its two
73 history levels are in `functions`, for checkpoints.
75 K defaults to sigma c_inf (1 - M^2) / L with sigma = 0.25 and L the
76 far-field distance `relax_length`.
78 Tangential derivatives at p = 0: the state is constant in each cell, so
79 ufl.grad of it vanishes and T would be identically zero (riemann2 would
80 reduce to riemann). The derivatives are therefore taken from a
81 Green-Gauss gradient of the state at the two history levels
82 (core/reconstruction.py): the model's own reconstruction when it has one
83 (RANS), otherwise a standalone one that is not part of the residual.
84 """
86 def __init__(self, model, farfield_tag, beta=0.5, K=None, relax_length=50.0, relax_sigma=0.25):
87 if model.dimensions != 2:
88 raise NotImplementedError("the transverse far field (riemann2) is 2D only")
89 mesh_obj = model.mesh_obj
90 self.mesh = mesh_obj.mesh
91 self.gamma = gamma = model.gamma
92 self.n = mesh_obj.n
93 self.farfield_tag = farfield_tag
94 ff_bc = model.farfield_conditions
95 V_s = model.functionspaces["V_scalar"]
96 self.delta_n = dolfinx.fem.Function(V_s, name="ff_delta_n")
97 self.delta_nm1 = dolfinx.fem.Function(V_s, name="ff_delta_nm1")
98 self.delta_next = dolfinx.fem.Function(V_s, name="ff_delta_next")
99 self.functions = [self.delta_n, self.delta_nm1]
100 if K is None:
101 K = relax_sigma * (gamma * ff_bc['p'] / ff_bc['rho']) ** 0.5 \
102 * (1.0 - ff_bc['M'] ** 2) / relax_length
103 self.beta, self.K = float(beta), float(K)
104 # this boundary equation's own copies of the step size and BDF
105 # coefficients, set in prepare()
106 self.dt = dolfinx.fem.Constant(self.mesh, 1.0)
107 self.a = tuple(dolfinx.fem.Constant(self.mesh, c) for c in (1.0, -1.0, 0.0))
108 self._coef = (self.dt, self.a, self.beta, self.K)
109 # extrapolation weight dt_n / dt_{n-1}
110 self.ext_w = dolfinx.fem.Constant(self.mesh, 1.0)
111 # p = 0: Green-Gauss gradients of U^n and U^{n-1} (set up in attach)
112 self.model = model
113 self._rec = None
114 self.G_n = self.G_nm1 = None
115 if model.poly_order == 0:
116 from dragonfly_sim.core.reconstruction import GreenGaussReconstruction
117 self._rec = getattr(model, "reconstruction", None) or \
118 GreenGaussReconstruction(model, muscl=False)
119 V_G = self._rec.G.function_space
120 self.G_n = dolfinx.fem.Function(V_G, name="ff_grad_n")
121 self.G_nm1 = dolfinx.fem.Function(V_G, name="ff_grad_nm1")
122 self.functions += [self.G_n, self.G_nm1]
124 c_inf = (gamma * ff_bc['p'] / ff_bc['rho']) ** 0.5
125 self.R_in = ufl.dot(ff_bc['u'], self.n) - 2.0 * c_inf / (gamma - 1.0) \
126 + self.delta_next
127 self._exprs = None
128 self._setup_projection(mesh_obj, V_s)
130 def transverse_term(self, U, G=None):
131 """T = u_t dR-/dtau - c du_t/dtau of state U on the far field (UFL),
132 with grad U from the field G (n_state x 2) when given, else ufl.grad."""
133 gamma, n = self.gamma, self.n
134 n_state = U.ufl_shape[0]
135 Um = ufl.as_vector([U[i] for i in range(4)]) if n_state > 4 else U # mean flow
136 rho = Um[0]
137 vel = ufl.as_vector([Um[1] / rho, Um[2] / rho])
138 c = ufl.sqrt(gamma * pressure(Um, gamma) / rho)
139 tau = ufl.as_vector([-n[1], n[0]])
140 if G is None:
141 grad_vel, grad_c = ufl.grad(vel), ufl.grad(c)
142 else:
143 from dragonfly_sim.utils.NS_utils import primitive_gradients
144 grad_vel, grad_T, _, _ = primitive_gradients(U, G, gamma, 4) # grad u, grad(p/rho)
145 grad_c = gamma / (2.0 * c) * grad_T
146 dvel_dtau = ufl.dot(grad_vel, tau) # grad(vel)[i, j] = d u_i / d x_j
147 # R-_int = u.n - 2c/(g-1); n is constant along the straight far field
148 dR_dtau = ufl.dot(dvel_dtau, n) - 2.0 / (gamma - 1.0) * ufl.dot(grad_c, tau)
149 return ufl.dot(vel, tau) * dR_dtau - c * ufl.dot(dvel_dtau, tau)
151 def bind_history(self, U_n, U_nm1):
152 # Separate BDF1 and BDF2 expressions rather than one with an
153 # extrapolation weight: 0 * T(U^{n-1}) would still be evaluated and
154 # turn NaN if U^{n-1} were not a valid state.
155 dt, (a0, a1, a2), beta, K = self._coef
156 T_n = self.transverse_term(U_n, self.G_n)
157 # linear extrapolation to t^{n+1}; w = dt_n / dt_{n-1} (1 at constant dt)
158 w = self.ext_w
159 T_star = {1: T_n, 2: (1.0 + w) * T_n - w * self.transverse_term(U_nm1, self.G_nm1)}
160 self.delta_new = {k: (-a1 * self.delta_n - a2 * self.delta_nm1 - dt * beta * T)
161 / (a0 + K * dt) for k, T in T_star.items()}
162 self._exprs = {k: dolfinx.fem.Expression(e, self._s[:, None])
163 for k, e in self.delta_new.items()}
165 def _setup_projection(self, mesh_obj, V_s):
166 mesh = self.mesh
167 tdim = mesh.topology.dim
168 n_owned = mesh.topology.index_map(tdim).size_local
169 mesh.topology.create_connectivity(tdim - 1, tdim)
170 mesh.topology.create_connectivity(tdim, tdim - 1)
171 f2c = mesh.topology.connectivity(tdim - 1, tdim)
172 c2f = mesh.topology.connectivity(tdim, tdim - 1)
174 # 3-point Gauss-Legendre on the reference facet [0, 1]
175 g = np.sqrt(3.0 / 5.0)
176 s = np.array([0.5 - 0.5 * g, 0.5, 0.5 + 0.5 * g])
177 w = np.array([5.0, 8.0, 5.0]) / 18.0
178 self._s = s
180 element = V_s.element.basix_element
181 topo = basix.topology(element.cell_type)[tdim - 1]
182 ref = basix.geometry(element.cell_type)
183 tab = [element.tabulate(0, ref[f[0]] + s[:, None] * (ref[f[1]] - ref[f[0]]))[0, :, :, 0]
184 for f in topo]
186 per_cell = {}
187 for f in mesh_obj.meshtags.find(self.farfield_tag):
188 cell = f2c.links(f)[0]
189 if cell >= n_owned:
190 continue
191 per_cell.setdefault(cell, []).append(int(np.where(c2f.links(cell) == f)[0][0]))
193 x_nodes = mesh.geometry.x[:, :tdim]
194 gdofs = mesh.geometry.dofmaps[0]
195 cmap = mesh.geometry.cmaps[0]
196 # Expression.eval on facets applies dolfinx's facet reflection, which
197 # makes both neighbours see a facet's points in the same order. On a
198 # reflected facet the point at s therefore lies at reference
199 # va + (1-s)(vb-va). Detect the orientation per facet from the
200 # physical coordinates instead of relying on the permutation bits.
201 ent_all = np.array([(c, lf) for c, lfs in per_cell.items() for lf in lfs],
202 dtype=np.int32).reshape(-1, 2)
203 # Building an Expression is a collective JIT call: every rank must
204 # build it, including ranks without far-field facets.
205 x_expr = dolfinx.fem.Expression(ufl.SpatialCoordinate(mesh), s[:, None])
206 x_eval = x_expr.eval(mesh, ent_all) if len(ent_all) else np.zeros((0, len(s), tdim))
207 entities, ops, cell_dofs = [], [], []
208 k_ent = 0
209 for cell, lfs in per_cell.items():
210 rows, weights = [], []
211 nodes = mesh.geometry.x[gdofs[cell]]
212 for lf in lfs:
213 va, vb = (x_nodes[gdofs[cell][topo[lf][k]]] for k in (0, 1))
214 ra, rb = ref[topo[lf][0]], ref[topo[lf][1]]
215 fwd = cmap.push_forward(ra + s[:, None] * (rb - ra), nodes)[:, :tdim]
216 got = x_eval[k_ent][:, :tdim]
217 k_ent += 1
218 scale = np.linalg.norm(vb - va)
219 if np.allclose(got, fwd, rtol=0, atol=1e-10 * scale):
220 rows.append(tab[lf])
221 elif np.allclose(got, fwd[::-1], rtol=0, atol=1e-10 * scale):
222 rows.append(tab[lf][::-1])
223 else:
224 raise RuntimeError("cannot match facet quadrature points of cell "
225 "{} facet {}".format(cell, lf))
226 weights.append(w * scale)
227 entities.append((cell, lf))
228 B = np.vstack(rows)
229 sw = np.sqrt(np.concatenate(weights))
230 # weighted least squares, minimum norm for dofs no facet sees
231 ops.append(np.linalg.pinv(sw[:, None] * B) * sw[None, :])
232 cell_dofs.append(V_s.dofmap.cell_dofs(cell))
233 self._entities = np.array(entities, dtype=np.int32).reshape(-1, 2)
234 self._ops, self._cell_dofs = ops, cell_dofs
235 self._facets_per_cell = [len(l) for l in per_cell.values()]
236 self.n_facets = mesh.comm.allreduce(len(self._entities), MPI.SUM)
238 def attach(self, integrator):
239 """Register with a TimeIntegrator: bind the solution history and
240 receive prepare()/commit() calls from it."""
241 self.bind_history(integrator.term.U_n, integrator.term.U_nm1)
242 integrator.boundary_updaters.append(self)
243 if self._rec is not None:
244 # both history gradients from the current state (U^n = U^{n-1} at
245 # the start; a restart overwrites them from the checkpoint)
246 self._gradient_into(self.G_n)
247 self.G_nm1.x.array[:] = self.G_n.x.array
249 def _gradient_into(self, target):
250 self._rec.update_gradient()
251 target.x.array[:] = self._rec.G.x.array
252 target.x.scatter_forward()
254 def prepare(self, order, ratio=None, dt=None):
255 """Before a step's Newton solve: delta^{n+1} into delta_next, for a
256 BDF `order` step of size dt with step ratio dt_n/dt_{n-1}. Collective.
257 Returns max |delta^{n+1}| over the far-field quadrature points."""
258 from dragonfly_sim.core.time_integration import bdf_coefficients
259 for c, val in zip(self.a, bdf_coefficients(order, ratio)):
260 c.value = val
261 if dt is not None:
262 self.dt.value = dt
263 self.ext_w.value = 1.0 if ratio is None else ratio
264 expr = self._exprs[2 if order == 2 else 1]
265 vals = expr.eval(self.mesh, self._entities) if len(self._entities) \
266 else np.zeros((0, 3))
267 new = np.zeros_like(self.delta_next.x.array)
268 i = 0
269 for op, dofs, nf in zip(self._ops, self._cell_dofs, self._facets_per_cell):
270 new[dofs] = op @ vals[i:i + nf].reshape(-1)
271 i += nf
272 self.delta_next.x.array[:] = new
273 self.delta_next.x.scatter_forward()
274 return self.mesh.comm.allreduce(float(np.abs(vals).max(initial=0.0)), MPI.MAX)
276 def commit(self):
277 """After the step is accepted (u_vec = U^{n+1}): shift the boundary
278 history."""
279 self.delta_nm1.x.array[:] = self.delta_n.x.array
280 self.delta_n.x.array[:] = self.delta_next.x.array
281 if self._rec is not None:
282 self.G_nm1.x.array[:] = self.G_n.x.array
283 self._gradient_into(self.G_n)