src/lbmpy/boundaries/boundaryconditions.py

@@ -1012,21 +1012,46 @@ class FixedDensity(LbBoundary):
     Args:
         density: value of the density which should be set.
         name: optional name of the boundary.
+        data_type: data type of the density value. default is double
     """
 
-    def __init__(self, density, name=None):
+    def __init__(self, density, name=None, data_type='double'):
         if name is None:
             name = "Fixed Density " + str(density)
         self.density = density
+        self.density_is_callable = callable(self.density)
+        self.data_type = data_type
 
         super(FixedDensity, self).__init__(name)
 
+    @property
+    def additional_data(self):
+        """ In case of the UBB boundary additional data is a velocity vector. This vector is added to each cell to
+            realize velocity profiles for the inlet."""
+        if self.density_is_callable:
+            return [('rho', create_type(self.data_type))]
+        else:
+            return []
+
+    @property
+    def additional_data_init_callback(self):
+        """Initialise additional data of the boundary. For an example see
+            `tutorial 02 <https://pycodegen.pages.i10git.cs.fau.de/lbmpy/notebooks/02_tutorial_boundary_setup.html>`_
+            or lbmpy.geometry.add_pipe_inflow_boundary"""
+        if self.density_is_callable:
+            return self.density
+
     def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
         def remove_asymmetric_part_of_main_assignments(assignment_collection, degrees_of_freedom):
             new_main_assignments = [Assignment(a.lhs, get_symmetric_part(a.rhs, degrees_of_freedom))
                                     for a in assignment_collection.main_assignments]
             return assignment_collection.copy(new_main_assignments)
 
+        if self.density_is_callable:
+            density = index_field[0]('rho')
+        else:
+            density = self.density
+
         cqc = lb_method.conserved_quantity_computation
         velocity = cqc.velocity_symbols
         symmetric_eq = remove_asymmetric_part_of_main_assignments(lb_method.get_equilibrium(),
@@ -1037,7 +1062,6 @@ class FixedDensity(LbBoundary):
         simplification = create_simplification_strategy(lb_method)
         symmetric_eq = simplification(symmetric_eq)
 
-        density = self.density
         equilibrium_input = cqc.equilibrium_input_equations_from_init_values(density=density)
         equilibrium_input = equilibrium_input.new_without_subexpressions()
         equilibrium_input = equilibrium_input.main_assignments_dict
@@ -1076,6 +1100,7 @@ class DiffusionDirichlet(LbBoundary):
     """
 
     def __init__(self, concentration, velocity_field=None, name=None, data_type='double'):
+        warn("DiffusionDirichlet BC is deprecated. Please use FixedDensity instead", DeprecationWarning)
         if name is None:
             name = "DiffusionDirichlet"
         self.concentration = concentration

src/lbmpy/phasefield_allen_cahn/analytical.py

@@ -15,23 +15,49 @@ def analytic_rising_speed(gravitational_acceleration, bubble_diameter, viscosity
     return result
 
 
-def analytical_solution_microchannel(reference_length, length_x, length_y,
+def analytical_solution_microchannel(reference_length, liquid_top_height, liquid_bottom_height,
                                      kappa_top, kappa_bottom,
                                      t_h, t_c, t_0,
                                      reference_surface_tension, dynamic_viscosity_light_phase,
                                      transpose=True):
     """
-    https://www.sciencedirect.com/science/article/pii/S0021999113005986
+    Analytic solution for thermocapillary-driven convection of superimposed fluids at
+    zero Reynolds and Marangoni numbers https://doi.org/10.1016/j.ijthermalsci.2010.02.003
+    The analytic solution is evaluated at the cell centers.
+
+                        L
+    +----------------------------------------+
+    |                                        |
+    |              Liquid 1                  | liquid_top_height
+    |                                        |
+    |========================================|  reference_length
+    |                                        |
+    |              Liquid 2                  | liquid_bottom_height
+    |                                        |
+    +----------------------------------------+
+
+    reference_length: the height is taken as reference length
+    liquid_top_height: the height of top liquid
+    liquid_bottom_height: the height of the bottom liquid
+    kappa_top: thermal diffusivity of the top liquid
+    kappa_bottom: thermal diffusivity of the bottom liquid
+    t_h: maximal temperature
+    t_0: minimal temperature
+    t_c: base temperature
+    reference_surface_tension: sigma_ref as defined in the original derivation of the analytic solution
+    dynamic_viscosity_light_phase: dynamic viscosity light phase
+    transpose: return solution arrays in transposed form
     """
     l_ref = reference_length
+    length_x = int(2 * l_ref)
     sigma_t = reference_surface_tension
 
     kstar = kappa_top / kappa_bottom
-    mp = (l_ref // 2) - 1
+    mp = l_ref // 2
 
     w = pi / l_ref
-    a = mp * w
-    b = mp * w
+    a = liquid_top_height * w
+    b = liquid_bottom_height * w
 
     f = 1.0 / (kstar * sinh(b) * cosh(a) + sinh(a) * cosh(b))
     g = sinh(a) * f
@@ -50,8 +76,8 @@ def analytical_solution_microchannel(reference_length, length_x, length_y,
 
     umax = -1.0 * (t_0 * sigma_t / dynamic_viscosity_light_phase) * g * h
     jj = 0
-    xx = np.linspace(-l_ref - 0.5, l_ref - 0.5, length_x)
-    yy = np.linspace(-mp, mp, length_y)
+    xx = np.linspace(-l_ref + 0.5, l_ref - 0.5, length_x)
+    yy = np.linspace(-(l_ref // 2) + 0.5, (l_ref // 2) - 0.5, l_ref)
     u_x = np.zeros([len(xx), len(yy)])
     u_y = np.zeros([len(xx), len(yy)])
     t_a = np.zeros([len(xx), len(yy)])