Commit 864ad922 authored by Markus Holzer's avatar Markus Holzer
Browse files

Added diffusion dirichlet

parent 90a303c8
Pipeline #34179 failed with stages
in 5 minutes and 49 seconds
from lbmpy.boundaries.boundaryconditions import (
UBB, FixedDensity, DiffusionDirichlet, ConcentrationBoundary, SimpleExtrapolationOutflow,
UBB, FixedDensity, DiffusionDirichlet, SimpleExtrapolationOutflow,
ExtrapolationOutflow, NeumannByCopy, NoSlip, StreamInConstant)
from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
__all__ = ['NoSlip', 'UBB', 'SimpleExtrapolationOutflow', 'ExtrapolationOutflow',
'FixedDensity', 'DiffusionDirichlet', 'ConcentrationBoundary', 'NeumannByCopy',
'FixedDensity', 'DiffusionDirichlet', 'NeumannByCopy',
'LatticeBoltzmannBoundaryHandling', 'StreamInConstant']
from lbmpy.advanced_streaming.utility import AccessPdfValues, Timestep
from pystencils.simp.assignment_collection import AssignmentCollection
from pystencils import Assignment, Field
from pystencils.sympyextensions import scalar_product
from lbmpy.boundaries.boundaryhandling import LbmWeightInfo
from pystencils.data_types import create_type
from pystencils.sympyextensions import get_symmetric_part
......@@ -461,42 +462,7 @@ class FixedDensity(LbBoundary):
# end class FixedDensity
class DiffusionDirichlet(LbBoundary):
"""Boundary condition for advection-diffusion problems that fixes the concentration at the obstacle.
Args:
concentration: value of the concentration which should be set.
name: optional name of the boundary.
"""
def __init__(self, concentration, name=None):
if name is None:
name = "Diffusion Dirichlet " + str(concentration)
self.concentration = concentration
super(DiffusionDirichlet, self).__init__(name)
def get_additional_code_nodes(self, lb_method):
"""Return a list of code nodes that will be added in the generated code before the index field loop.
Args:
lb_method: Lattice Boltzmann method. See :func:`lbmpy.creationfunctions.create_lb_method`
Returns:
list containing LbmWeightInfo
"""
return [LbmWeightInfo(lb_method)]
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
w_dir = LbmWeightInfo.weight_of_direction(dir_symbol, lb_method)
return [Assignment(f_in(inv_dir[dir_symbol]),
2 * w_dir * self.concentration - f_out(dir_symbol))]
# end class DiffusionDirichlet
class ConcentrationBoundary(LbBoundary):
"""Concentration boundary which is used for concentration or thermal boundary conditions of convection-diffusion
equation Base on https://doi.org/10.1103/PhysRevE.85.016701.
......@@ -504,17 +470,20 @@ class ConcentrationBoundary(LbBoundary):
concentration: can either be a constant, an access into a field, or a callback function.
The callback functions gets a numpy record array with members, 'x','y','z', 'dir' (direction)
and 'concentration' which has to be set to the desired velocity of the corresponding link
velocity_field: if velocity field is given the boundary value is approximated by using the discrete equilibrium.
name: optional name of the boundary.
data_type: data type of the concentration value. default is double
"""
def __init__(self, concentration, name=None, data_type='double'):
def __init__(self, concentration, velocity_field=None, name=None, data_type='double'):
if name is None:
name = "ConcentrationBoundary"
name = "DiffusionDirichlet"
self.concentration = concentration
self.data_type = data_type
self.concentration_is_callable = callable(self.concentration)
self.velocity_field = velocity_field
super(ConcentrationBoundary, self).__init__(name)
super(DiffusionDirichlet, self).__init__(name)
@property
def additional_data(self):
......@@ -542,16 +511,32 @@ class ConcentrationBoundary(LbBoundary):
Returns:
list containing LbmWeightInfo
"""
return [LbmWeightInfo(lb_method)]
if self.velocity_field:
return [LbmWeightInfo(lb_method), NeighbourOffsetArrays(lb_method.stencil)]
else:
return [LbmWeightInfo(lb_method)]
def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
w_dir = LbmWeightInfo.weight_of_direction(dir_symbol, lb_method)
if self.concentration_is_callable:
concentration = index_field[0]('concentration')
else:
concentration = self.concentration
return [Assignment(f_in(inv_dir[dir_symbol]),
2 * w_dir * concentration - f_out(dir_symbol))]
if self.velocity_field:
neighbour_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
u = self.velocity_field
cs = sp.Rational(1, 3)
equilibrium = (1 + scalar_product(neighbour_offset, u.center_vector)**2 / (2 * cs**4)
- scalar_product(u.center_vector, u.center_vector) / (2 * cs**2))
else:
equilibrium = sp.Rational(1, 1)
result = [Assignment(f_in(inv_dir[dir_symbol]),
2.0 * w_dir * concentration * equilibrium - f_out(dir_symbol))]
return result
# end class DiffusionDirichlet
......
import pytest
from lbmpy.creationfunctions import create_lb_function
from lbmpy.scenarios import create_lid_driven_cavity
from pystencils import show_code
from pystencils import get_code_str
def test_creation():
@pytest.mark.parametrize('method', ['srt', 'trt', 'mrt', 'cumulant'])
def test_creation(method):
"""Simple test that makes sure that only float variables are created"""
func = create_lb_function(method='srt', relaxation_rate=1.5,
func = create_lb_function(method=method, relaxation_rates=[1.5]*9,
optimization={'double_precision': False})
code = str(show_code(func.ast))
assert 'double' not in code
code_str = get_code_str(func)
assert 'double' not in code_str
assert 'float' in code_str
def test_scenario():
sc = create_lid_driven_cavity((16, 16, 8), relaxation_rate=1.5,
optimization={'double_precision': False})
sc.run(1)
code_str = str(show_code(sc.ast))
code_str = get_code_str(sc.ast)
assert 'double' not in code_str
assert 'float' in code_str
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment