Fixed some bugs related to the TFM solver

opencmp/models/tfm.py

@@ -16,6 +16,7 @@
 ########################################################################################################################
 
 import logging
+import math
 from typing import Dict, List, Optional, Union
 
 import ngsolve as ngs
@@ -504,6 +505,8 @@ def construct_bilinear_time_ODE(self,
         a += (ngs.div(ud) * q) * ngs.dx
         a += (-Ac * ngs.div(ud) * q) * ngs.dx
         a += (-ngs.grad(Ac) * ud * q) * ngs.dx
+        # the following term ccount for discontinuous alpha_c in the conservative mixture flux.
+        a += (-avg(q) * jump(Ac) * ((uc - ud) * n)) * ngs.dx(skeleton=True)
         if self.mean_zero_pressure:
             pressure_mean, pressure_mean_test = U[-1], V[-1]
             a += (pressure_mean * q + p * pressure_mean_test) * ngs.dx
@@ -608,6 +611,20 @@ def construct_bilinear_time_ODE(self,
             a += (dt * drag_factor_d * (ud - uc) * vd) * ngs.dx
             a += (dt * drag_factor_c * (ud - uc) * vc) * ngs.dx
 
+        if self.Disp_switch and self.drag_switch:
+            Cd = self._get_drag_coeff(wd, wc, Ad, ts)
+            H = 1 - 1.166*Ad + 0.5*Ad**2
+            Cdis = self.Cdis[ts]
+            relative_wind = wd - wc
+            relative_trial = ud - uc
+            slip_product = ngs.InnerProduct(relative_wind, relative_trial)
+            disp_factor_d = 0.75 * Cd * Cdis * rho_c/rho_d * H
+            disp_factor_c = -0.75 * Cd * Cdis * Ad/(Ac + 1e-30) * H
+            # The physical forces were formerly RHS terms, so their
+            # velocity-dependent linearizations enter the left with minus signs.
+            a += (-dt * disp_factor_d * slip_product * ngs.grad(Ac) * vd) * ngs.dx
+            a += (-dt * disp_factor_c * slip_product * ngs.grad(Ac) * vc) * ngs.dx
+
         # — Virtual mass: a_dt (NOT dt-scaled) + a (dt-scaled)
         if self.VM_switch:
             C_VM = self.C_VM[ts]
@@ -683,6 +700,11 @@ def construct_bilinear_time_ODE(self,
             a += (dt * lift_mom_d) * ngs.dx
             a += (dt * lift_mom_c) * ngs.dx
 
+        if self.injection_switch and self.inj_mass_flowrate:
+            dx_inj = ngs.dx(definedon=self.mesh.Materials(self.injection_region))
+            Ad_safe = Max(Ad, ngs.CoefficientFunction(1e-5))
+            a += (dt * self.inj_mass_flowrate / Ad_safe * ud * vd) * dx_inj
+
         return [a]
 
     def construct_bilinear_time_coefficient(self,
@@ -825,30 +847,19 @@ def construct_linear(self,
                      * self._ds(marker)
 
         # ============================================================
-        # 5. Laminar dispersion (linear source term)
-        # ============================================================
-        if self.Disp_switch and self.drag_switch:
-            Cd   = self._get_drag_coeff(wd, wc, Ad, ts)
-            H    = 1 - 1.166*(1 - Ac) + 0.5*(1 - Ac)**2
-            Cdis = self.Cdis[ts]
-            disp_d = 0.75*Cd*Cdis*rho_c/rho_d * H * ngs.Norm(wd - wc)**2
-            disp_c = -0.75*Cd*Cdis*Ad/(Ac + 1e-30) * H * ngs.Norm(wd - wc)**2
-            L += (dt * disp_d * ngs.grad(Ac) * vd) * ngs.dx
-            L += (dt * disp_c * ngs.grad(Ac) * vc) * ngs.dx
-
-        # ============================================================
-        # 6. Interior injection source (on the 'injection' material region)
+        # 5. Interior injection source (on the 'injection' material region)
         # ============================================================
         if self.injection_switch and self.inj_mass_flowrate:
             mdot = self.inj_mass_flowrate
             dx_inj = ngs.dx(definedon=self.mesh.Materials(self.injection_region))
             L += (dt * -mdot * r) * dx_inj          # dispersed-phase mass
-            L += (dt *  mdot * q) * dx_inj          # mixture mass
+            L += (mdot * q) * dx_inj                 # mixture mass (not dt-scaled)
             if self.mesh.dim == 2:
                 inj_vec = ngs.CoefficientFunction((0, self.inj_velocity))
             else:
                 inj_vec = ngs.CoefficientFunction((0, self.inj_velocity, 0))
-            L += (dt * mdot * inj_vec * vc) * dx_inj  # continuous-phase momentum
+            Ad_safe = Max(Ad, ngs.CoefficientFunction(1e-5))
+            L += (dt * mdot / Ad_safe * inj_vec * vd) * dx_inj
 
         return [L]
 
@@ -911,6 +922,10 @@ def solve_single_step(self,
                 ac = gfu.components[comp['alpha_c']]
                 self._lim.bezier_bound(ac, ac.space, ord_s, (0.0, 1.0))
 
+            if not math.isfinite(gfu.vec.Norm()):
+                logging.warning('TFM Picard produced a non-finite solution.')
+                return False
+
             # Match INS convergence testing: evaluate each local-error component
             # independently and exclude pressure.
             converged = True
@@ -924,7 +939,8 @@ def solve_single_step(self,
                 solution_norm = mean(gfu.components[i], self.mesh)
                 tolerance = (self.abs_nonlinear_tolerance
                              + self.rel_nonlinear_tolerance * solution_norm)
-                if err >= tolerance:
+                if not (math.isfinite(err) and math.isfinite(tolerance)
+                        and err < tolerance):
                     converged = False
                     break
 

pytests/models/test_tfm.py

@@ -4,13 +4,64 @@
 
 import ngsolve as ngs
 import pytest
+from netgen.geom2d import unit_square
 
 from opencmp.config_functions import ConfigParser
 from opencmp.config_functions.boundary_conditions import BCFunctions
 from opencmp.models.tfm import TwoFluidModel
 import opencmp.models.tfm as tfm_module
 
 
+def _injection_test_model() -> TwoFluidModel:
+    model = TwoFluidModel.__new__(TwoFluidModel)
+    model.mesh = ngs.Mesh(unit_square.GenerateMesh(maxh=0.4))
+    model.fes = ngs.FESpace([
+        ngs.HDiv(model.mesh, order=1, dgjumps=True),
+        ngs.HDiv(model.mesh, order=1, dgjumps=True),
+        ngs.L2(model.mesh, order=0, dgjumps=True),
+        ngs.L2(model.mesh, order=0, dgjumps=True),
+    ], dgjumps=True)
+    model.model_components = {'u_c': 0, 'u_d': 1, 'p': 2, 'alpha_c': 3}
+    model.UIter = ngs.GridFunction(model.fes)
+    model.UIter.components[3].Set(0.8)
+    model.n = ngs.specialcf.normal(model.mesh.dim)
+    model.IM = ngs.Id(model.mesh.dim)
+    model.penalty_interior = 10.0
+    model.rho_c = [1000.0]
+    model.rho_d = [1.0]
+    model.nu_c = [0.0]
+    model.nu_d = [0.0]
+    model.dp = [0.003]
+    model.BC = {}
+    model.dirichlet_names = {}
+    model.f = {}
+    model.diffusion_switch = False
+    model.canonical_form = 'B-TFM'
+    model.gravity = ngs.CoefficientFunction((0.0, 0.0))
+    model.drag_switch = False
+    model.VM_switch = False
+    model.Disp_switch = False
+    model.Lift_switch = False
+    model.mean_zero_pressure = False
+    model.injection_switch = True
+    model.inj_mass_flowrate = 0.3
+    model.injection_region = '.*'
+    model.inj_velocity = 2.0
+    return model
+
+
+def _boundary_trace_integral(mesh, expression) -> float:
+    """Integrate the adjacent volume trace of an L2 field on the boundary."""
+    space = ngs.L2(mesh, order=0)
+    test = space.TestFunction()
+    one = ngs.GridFunction(space)
+    one.Set(1.0)
+    form = ngs.LinearForm(space)
+    form += expression * test * ngs.ds(skeleton=True, bonus_intorder=10)
+    form.Assemble()
+    return ngs.InnerProduct(form.vec, one.vec)
+
+
 def _tfm_from_config(tmp_path: Path, tfm: dict[str, str], other=None) -> TwoFluidModel:
     config_path = tmp_path / 'config'
     config_path.write_text('')
@@ -187,3 +238,85 @@ def test_ime_closures_are_finite_and_zero_without_a_driver() -> None:
             magnitude_squared = ngs.Integrate(closure * closure, model.mesh)
             assert math.isfinite(magnitude_squared)
             assert magnitude_squared == pytest.approx(0.0, abs=1e-24)
+
+
+def test_injection_sources_use_consistent_equation_scaling_and_phase() -> None:
+    model = _injection_test_model()
+    dt = ngs.Parameter(2e-3)
+    U, V = model.fes.TrialFunction(), model.fes.TestFunction()
+
+    a = ngs.BilinearForm(model.fes)
+    a += model.construct_bilinear_time_ODE(U, V, dt, 0)[0]
+    L = ngs.LinearForm(model.fes)
+    L += model.construct_linear(V, None, dt, 0)[0]
+    a.Assemble()
+    L.Assemble()
+
+    probes = []
+    for component in range(4):
+        probe = ngs.GridFunction(model.fes)
+        if component < 2:
+            probe.components[component].Set(ngs.CoefficientFunction((0.0, 1.0)))
+        else:
+            probe.components[component].Set(1.0)
+        probes.append(probe)
+
+    continuous_momentum = ngs.InnerProduct(L.vec, probes[0].vec)
+    dispersed_momentum = ngs.InnerProduct(L.vec, probes[1].vec)
+    mixture_mass = ngs.InnerProduct(L.vec, probes[2].vec)
+    dispersed_mass = ngs.InnerProduct(L.vec, probes[3].vec)
+
+    assert continuous_momentum == pytest.approx(0.0, abs=1e-14)
+    assert dispersed_momentum == pytest.approx(6e-3)
+    assert mixture_mass == pytest.approx(0.3)
+    assert dispersed_mass == pytest.approx(-6e-4)
+
+    injected_velocity = ngs.GridFunction(model.fes)
+    injected_velocity.components[1].Set(ngs.CoefficientFunction((0.0, 2.0)))
+    dilution = ngs.InnerProduct(a.mat * injected_velocity.vec, probes[1].vec)
+    assert dilution == pytest.approx(dispersed_momentum)
+
+
+def test_mixture_equation_includes_discontinuous_volume_fraction_flux() -> None:
+    model = _injection_test_model()
+    model.injection_switch = False
+    model.UIter.components[3].Set(ngs.IfPos(ngs.x - 0.5, 0.8, 0.2))
+    U, V = model.fes.TrialFunction(), model.fes.TestFunction()
+    a = ngs.BilinearForm(model.fes)
+    a += model.construct_bilinear_time_ODE(U, V, ngs.Parameter(2e-3), 0)[0]
+    a.Assemble()
+
+    state = ngs.GridFunction(model.fes)
+    state.components[0].Set(ngs.CoefficientFunction((1.0, 0.0)))
+    constant_q = ngs.GridFunction(model.fes)
+    constant_q.components[2].Set(1.0)
+
+    weak_mixture_flux = ngs.InnerProduct(a.mat * state.vec, constant_q.vec)
+    physical_boundary_flux = _boundary_trace_integral(
+        model.mesh, model.UIter.components[3] * state.components[0] * model.n)
+
+    assert physical_boundary_flux > 0.1
+    assert weak_mixture_flux == pytest.approx(physical_boundary_flux, abs=1e-12)
+
+
+def test_picard_rejects_nonfinite_solution() -> None:
+    model = _injection_test_model()
+    model.IC = ngs.GridFunction(model.fes)
+    model.IC.components[3].Set(0.8)
+    model.nonlinear_max_iters = 1
+    model.relax_factors = [1.0] * 4
+    model.slope_limiter = False
+    model.apply_dirichlet_bcs_to = lambda *_args: None
+
+    def nonfinite_solve(_a, _L, _preconditioner, gfu):
+        gfu.vec.FV().NumPy()[:] = float('nan')
+
+    model.linear_solve = nonfinite_solve
+
+    class Form:
+        def Assemble(self):
+            pass
+
+    gfu = ngs.GridFunction(model.fes)
+    gfu.vec.data = model.IC.vec
+    assert not model.solve_single_step([Form()], [Form()], [None], gfu)