Add Correction factor to outflow

doc/sphinx/lbmpy.bib

@@ -124,10 +124,22 @@ url = {https://doi.org/10.1016/j.jcp.2017.05.040},
 doi = {10.1016/j.jcp.2017.05.040},
 journal = {J. Comput. Phys.},
 month = {nov},
-pages = {862–888},
+pages = {862-888},
 numpages = {27}
 }
 
+@article{Alihussein2024,
+author = {Alihussein, Hussein and Prasannakumar, Adarsh and Geier, Martin and Krafczyk, Manfred},
+title = {{Direct cumulant lattice Boltzmann simulations of transitional flow in gyroidal structures including experimental validation}},
+journal = {Computers & Mathematics with Applications},
+volume = {157},
+pages = {159-172},
+year = {2024},
+issn = {0898-1221},
+doi = {10.1016/j.camwa.2023.12.029},
+url = {https://www.sciencedirect.com/science/article/pii/S0898122123005904},
+}
+
 @Article{Coreixas2019,
   title = {Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations},
   author = {Coreixas, Christophe and Chopard, Bastien and Latt, Jonas},

src/lbmpy/boundaries/boundaryconditions.py

@@ -900,9 +900,9 @@ class ExtrapolationOutflow(LbBoundary):
     array is used which stores them.
 
     .. math ::
-        f_{\overline{1}jkxyzt} = f_{\overline{1}jk(x - \Delta x)yz(t - \Delta t)} c \theta^{\frac{1}{2}}
+    f_{\overline{1}jkxyzt} = f_{\overline{1}jk(x - \Delta x)yz(t - \Delta t)} c \theta^{\frac{1}{2}}
         \frac{\Delta t}{\Delta x} + \left(1 - c \theta^{\frac{1}{2}} \frac{\Delta t}{\Delta x} \right)
-         f_{\overline{1}jk(x - \Delta x)yzt}
+         f_{\overline{1}jk(x - \Delta x)yzt} + \gamma_\rho w_{ijk} \Delta \rho
 
 
     Args:
@@ -917,11 +917,16 @@ class ExtrapolationOutflow(LbBoundary):
                          positional arguments, specifying the initial density on boundary nodes
         initial_velocity: tuple of floating point constants or callback taking spatial coordinates (x, y [,z]) as
                           positional arguments, specifying the initial velocity on boundary nodes
+        density: The ExtrapolationOutflow does not impose any coniditon on the density (or pressure). A simple
+                 correction term to impose a density condition is proposed by :cite:`Alihussein2024`
+        density_correction_factor: If density condition is added to the boundary (see last term of equation) then this
+                                   value determines how fast the boundary pushes the density to the imposed value.
     """
 
     def __init__(self, normal_direction, lb_method, dt=1, dx=1, name=None,
                  streaming_pattern='pull', zeroth_timestep=Timestep.BOTH,
-                 initial_density=None, initial_velocity=None, data_type='double'):
+                 initial_density=None, initial_velocity=None, data_type='double',
+                 density=None, density_correction_factor=0.01):
 
         self.lb_method = lb_method
         self.stencil = lb_method.stencil
@@ -947,6 +952,8 @@ class ExtrapolationOutflow(LbBoundary):
         self.equilibrium_calculation = None
 
         self.data_type = data_type
+        self.density = density
+        self.density_correction_factor = density_correction_factor
 
         if initial_density and initial_velocity:
             equilibrium = lb_method.get_equilibrium(conserved_quantity_equations=AssignmentCollection([]))
@@ -989,12 +996,17 @@ class ExtrapolationOutflow(LbBoundary):
             #   Initial fluid cell PDF values
             entry['pdf'] = pdf_acc.read_pdf(boundary_data.pdf_array, domain_cell, inv_dir)
             entry['pdf_nd'] = get_boundary_cell_pdfs(domain_cell, outflow_cell, inv_dir)
+            if self.density:
+                entry['rho'] = self.density - 1.0 if self.lb_method.conserved_quantity_computation.zero_centered_pdfs else self.density
+
 
     @property
     def additional_data(self):
         """Used internally only. For the ExtrapolationOutflow information of the previous PDF values is needed. This
         information is stored in the index vector."""
         data = [('pdf', create_type(self.data_type)), ('pdf_nd', create_type(self.data_type))]
+        if self.density:
+            data.append(('rho', create_type(self.data_type)))
         return data
 
     @property
@@ -1011,7 +1023,11 @@ class ExtrapolationOutflow(LbBoundary):
             list containing NeighbourOffsetArrays
 
         """
-        return [NeighbourOffsetArrays(lb_method.stencil)]
+        code_nodes = [NeighbourOffsetArrays(lb_method.stencil)]
+        if self.density:
+            code_nodes.append(LbmWeightInfo(lb_method, self.data_type))
+        return code_nodes
+    
 
     def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field, force_vector):
         subexpressions = []
@@ -1021,9 +1037,23 @@ class ExtrapolationOutflow(LbBoundary):
         neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
         tangential_offset = tuple(offset - normal for offset, normal in zip(neighbor_offset, self.normal_direction))
 
+        correction_term = 0
+        if self.density:
+            cqc = lb_method.conserved_quantity_computation
+            density_symbol = sp.Symbol("rho")
+            pdf_field_accesses = [f_out[tangential_offset](i) for i in range(len(lb_method.stencil))]
+            density_equations = cqc.output_equations_from_pdfs(pdf_field_accesses, {'density': density_symbol})
+            subexpressions.append(density_equations.all_assignments)
+
+            weight_info = LbmWeightInfo(lb_method, data_type=self.data_type)
+            weight_of_direction = weight_info.weight_of_direction
+            weight = weight_of_direction(dir_symbol, lb_method)
+            correction_term = self.density_correction_factor * weight * index_field[0]('rho')
+
+
         interpolated_pdf_sym = sp.Symbol('pdf_inter')
         interpolated_pdf_asm = Assignment(interpolated_pdf_sym, (index_field[0]('pdf') * (self.c * dtdx))
-                                          + ((sp.Number(1) - self.c * dtdx) * index_field[0]('pdf_nd')))
+                                          + ((sp.Number(1) - self.c * dtdx) * index_field[0]('pdf_nd')) + correction_term)
         subexpressions.append(interpolated_pdf_asm)
 
         asm = Assignment(f_in.center(inv_dir[dir_symbol]), interpolated_pdf_sym)
@@ -1035,6 +1065,10 @@ class ExtrapolationOutflow(LbBoundary):
         asm = Assignment(index_field[0]('pdf_nd'), interpolated_pdf_sym)
         boundary_assignments.append(asm)
 
+        if self.density:
+            asm = Assignment(index_field[0]('rho'), density_symbol)
+            boundary_assignments.append(asm)
+
         return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)