diff --git a/opencmp/models/tfm.py b/opencmp/models/tfm.py index 1cd6204..9edebb7 100644 --- a/opencmp/models/tfm.py +++ b/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 diff --git a/pytests/models/test_tfm.py b/pytests/models/test_tfm.py index 5bf3d16..d0610db 100644 --- a/pytests/models/test_tfm.py +++ b/pytests/models/test_tfm.py @@ -4,6 +4,7 @@ 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 @@ -11,6 +12,56 @@ 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)