diff --git a/lbmpy/phasefield_allen_cahn/derivatives.py b/lbmpy/phasefield_allen_cahn/derivatives.py
deleted file mode 100644
index 2a5f45b981a277489c99e8fa3730f4ca16a71154..0000000000000000000000000000000000000000
--- a/lbmpy/phasefield_allen_cahn/derivatives.py
+++ /dev/null
@@ -1,53 +0,0 @@
-def laplacian(phi_field):
- r"""
- Get a symbolic expression for the laplacian.
- Args:
- phi_field: the phase-field on which the laplacian is applied
- """
- lp_phi = 16.0 * ((phi_field[1, 0, 0]) + (phi_field[-1, 0, 0])
- + (phi_field[0, 1, 0]) + (phi_field[0, -1, 0])
- + (phi_field[0, 0, 1]) + (phi_field[0, 0, -1])) \
- + 1.0 * (
- (phi_field[1, 1, 1]) + (phi_field[-1, 1, 1])
- + (phi_field[1, -1, 1]) + (phi_field[-1, -1, 1])
- + (phi_field[1, 1, -1]) + (phi_field[-1, 1, -1])
- + (phi_field[1, -1, -1]) + (phi_field[-1, -1, -1])) \
- + 4.0 * (
- (phi_field[1, 1, 0]) + (phi_field[-1, 1, 0])
- + (phi_field[1, -1, 0]) + (phi_field[-1, -1, 0])
- + (phi_field[1, 0, 1]) + (phi_field[-1, 0, 1])
- + (phi_field[1, 0, -1]) + (phi_field[-1, 0, -1])
- + (phi_field[0, 1, 1]) + (phi_field[0, -1, 1])
- + (phi_field[0, 1, -1]) + (phi_field[0, -1, -1])) \
- - 152.0 * phi_field[0, 0, 0]
-
- lp_phi = lp_phi / 36
-
- return lp_phi
-
-
-def isotropic_gradient(phi_field):
- r"""
- Get a symbolic expression for the isotropic gradient of the phase-field
- Args:
- phi_field: the phase-field on which the isotropic gradient is applied
- """
- grad_phi_x = 16.00 * (phi_field[1, 0, 0] - phi_field[-1, 0, 0])\
- + (phi_field[1, 1, 1] - phi_field[-1, 1, 1] + phi_field[1, -1, 1] - phi_field[-1, -1, 1]
- + phi_field[1, 1, -1] - phi_field[-1, 1, -1] + phi_field[1, -1, -1] - phi_field[-1, -1, -1])\
- + 4.00 * (phi_field[1, 1, 0] - phi_field[-1, 1, 0] + phi_field[1, -1, 0] - phi_field[-1, -1, 0]
- + phi_field[1, 0, 1] - phi_field[-1, 0, 1] + phi_field[1, 0, -1] - phi_field[-1, 0, -1])
- grad_phi_y = 16.00 * (phi_field[0, 1, 0] - phi_field[0, -1, 0]) \
- + (phi_field[1, 1, 1] + phi_field[-1, 1, 1] - phi_field[1, -1, 1] - phi_field[-1, -1, 1]
- + phi_field[1, 1, -1] + phi_field[-1, 1, -1] - phi_field[1, -1, -1] - phi_field[-1, -1, -1])\
- + 4.00 * (phi_field[1, 1, 0] + phi_field[-1, 1, 0] - phi_field[1, -1, 0] - phi_field[-1, -1, 0]
- + phi_field[0, 1, 1] - phi_field[0, -1, 1] + phi_field[0, 1, -1] - phi_field[0, -1, -1])
- grad_phi_z = 16.00 * (phi_field[0, 0, 1] - phi_field[0, 0, -1]) \
- + (phi_field[1, 1, 1] + phi_field[-1, 1, 1] + phi_field[1, -1, 1] + phi_field[-1, -1, 1]
- - phi_field[1, 1, -1] - phi_field[-1, 1, -1] - phi_field[1, -1, -1] - phi_field[-1, -1, -1])\
- + 4.00 * (phi_field[1, 0, 1] + phi_field[-1, 0, 1] - phi_field[1, 0, -1] - phi_field[-1, 0, -1]
- + phi_field[0, 1, 1] + phi_field[0, -1, 1] - phi_field[0, 1, -1] - phi_field[0, -1, -1])
-
- grad = [grad_phi_x / 72, grad_phi_y / 72, grad_phi_z / 72]
-
- return grad
diff --git a/lbmpy/phasefield_allen_cahn/kernel_equations.py b/lbmpy/phasefield_allen_cahn/kernel_equations.py
index 985702f424aa139040f862429fc3ef73781a5252..b6904aa90b2900d55358786bab1c5f457efa2572 100644
--- a/lbmpy/phasefield_allen_cahn/kernel_equations.py
+++ b/lbmpy/phasefield_allen_cahn/kernel_equations.py
@@ -4,7 +4,6 @@ from lbmpy.fieldaccess import StreamPushTwoFieldsAccessor, CollideOnlyInplaceAcc
from pystencils.fd.derivation import FiniteDifferenceStencilDerivation
from lbmpy.maxwellian_equilibrium import get_weights
from pystencils import Assignment, AssignmentCollection
-from lbmpy.phasefield_allen_cahn.derivatives import laplacian, isotropic_gradient
import sympy as sp
import numpy as np
@@ -15,30 +14,34 @@ def chemical_potential_symbolic(phi_field, stencil, beta, kappa):
Get a symbolic expression for the chemical potential according to equation (5) in PhysRevE.96.053301.
Args:
phi_field: the phase-field on which the chemical potential is applied
- stencil: stencil of the lattice Boltzmann step
+ stencil: stencil to derive the finite difference for the laplacian (2nd order isotropic)
beta: coefficient related to surface tension and interface thickness
kappa: coefficient related to surface tension and interface thickness
"""
dimensions = len(stencil[0])
-
- deriv = FiniteDifferenceStencilDerivation((0, 0), stencil)
- # assume the stencil is symmetric
- deriv.assume_symmetric(0)
- deriv.assume_symmetric(1)
- if dimensions == 3:
- deriv.assume_symmetric(2)
-
- # set weights for missing degrees of freedom in the calculation
- if len(stencil) == 9:
- deriv.set_weight((1, 1), sp.Rational(1, 6))
- deriv.set_weight((0, -1), sp.Rational(2, 3))
- deriv.set_weight((0, 0), sp.Rational(-10, 3))
-
- # assume the stencil is isotropic
- res = deriv.get_stencil(isotropic=True)
- lap = res.apply(phi_field.center)
- else:
- lap = laplacian(phi_field)
+ q = len(stencil)
+ lap = sp.simplify(0)
+ for i in range(dimensions):
+ deriv = FiniteDifferenceStencilDerivation((i, i), stencil)
+ for j in range(dimensions):
+ # assume the stencil is symmetric
+ deriv.assume_symmetric(dim=j, anti_symmetric=False)
+
+ # set weights for missing degrees of freedom in the calculation and assume the stencil is isotropic
+ if q == 9:
+ res = deriv.get_stencil(isotropic=True)
+ lap += res.apply(phi_field.center)
+ elif q == 15:
+ deriv.set_weight((0, 0, 0), sp.Rational(-32, 27))
+ res = deriv.get_stencil(isotropic=True)
+ lap += res.apply(phi_field.center)
+ elif q == 19:
+ res = deriv.get_stencil(isotropic=True)
+ lap += res.apply(phi_field.center)
+ else:
+ deriv.set_weight((0, 0, 0), sp.Rational(-38, 27))
+ res = deriv.get_stencil(isotropic=True)
+ lap += res.apply(phi_field.center)
# get the chemical potential
mu = 4.0 * beta * phi_field.center * (phi_field.center - 1.0) * (phi_field.center - 0.5) - \
@@ -51,9 +54,10 @@ def isotropic_gradient_symbolic(phi_field, stencil):
Get a symbolic expression for the isotropic gradient of the phase-field
Args:
phi_field: the phase-field on which the isotropic gradient is applied
- stencil: stencil of the lattice Boltzmann step
+ stencil: stencil to derive the finite difference for the gradient (2nd order isotropic)
"""
dimensions = len(stencil[0])
+ q = len(stencil)
deriv = FiniteDifferenceStencilDerivation((0,), stencil)
deriv.assume_symmetric(0, anti_symmetric=True)
@@ -61,39 +65,55 @@ def isotropic_gradient_symbolic(phi_field, stencil):
if dimensions == 3:
deriv.assume_symmetric(2, anti_symmetric=False)
- if len(stencil) == 9:
+ # set weights for missing degrees of freedom in the calculation and assume the stencil is isotropic
+ # furthermore the stencils gets rotated to get the y and z components
+ if q == 9:
res = deriv.get_stencil(isotropic=True)
grad = [res.apply(phi_field.center), res.rotate_weights_and_apply(phi_field.center, (0, 1)), 0]
- elif len(stencil) == 19:
- deriv.set_weight((0, 0, 0), sp.Integer(0))
+ elif q == 15:
+ res = deriv.get_stencil(isotropic=True)
+ grad = [res.apply(phi_field.center),
+ res.rotate_weights_and_apply(phi_field.center, (0, 1)),
+ res.rotate_weights_and_apply(phi_field.center, (1, 2))]
+ elif q == 19:
+ deriv.set_weight((0, 0, 0), sp.sympify(0))
deriv.set_weight((1, 0, 0), sp.Rational(1, 6))
- deriv.set_weight((1, 1, 0), sp.Rational(1, 12))
res = deriv.get_stencil(isotropic=True)
grad = [res.apply(phi_field.center),
- res.rotate_weights_and_apply(phi_field.center, (1, 0)),
- res.rotate_weights_and_apply(phi_field.center, (2, 1))]
+ res.rotate_weights_and_apply(phi_field.center, (0, 1)),
+ res.rotate_weights_and_apply(phi_field.center, (1, 2))]
else:
- grad = isotropic_gradient(phi_field)
+ deriv.set_weight((0, 0, 0), sp.sympify(0))
+ deriv.set_weight((1, 0, 0), sp.Rational(2, 9))
+
+ res = deriv.get_stencil(isotropic=True)
+ grad = [res.apply(phi_field.center),
+ res.rotate_weights_and_apply(phi_field.center, (0, 1)),
+ res.rotate_weights_and_apply(phi_field.center, (1, 2))]
return grad
-def normalized_isotropic_gradient_symbolic(phi_field, stencil):
+def normalized_isotropic_gradient_symbolic(phi_field, stencil, fd_stencil=None):
r"""
Get a symbolic expression for the normalized isotropic gradient of the phase-field
Args:
phi_field: the phase-field on which the normalized isotropic gradient is applied
stencil: stencil of the lattice Boltzmann step
+ fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
+ field. If it is not given the stencil of the LB method will be applied.
"""
- iso_grad = isotropic_gradient_symbolic(phi_field, stencil)
- tmp = (sum(map(lambda x: x * x, isotropic_gradient_symbolic(phi_field, stencil))) + 1.e-32) ** 0.5
+ if fd_stencil is None:
+ fd_stencil = stencil
- result = [x / tmp for x in iso_grad]
+ tmp = (sum(map(lambda x: x * x, isotropic_gradient_symbolic(phi_field, fd_stencil))) + 1.e-32) ** 0.5
+
+ result = [x / tmp for x in isotropic_gradient_symbolic(phi_field, fd_stencil)]
return result
-def pressure_force(phi_field, stencil, density_heavy, density_light):
+def pressure_force(phi_field, stencil, density_heavy, density_light, fd_stencil=None):
r"""
Get a symbolic expression for the pressure force
Args:
@@ -101,13 +121,18 @@ def pressure_force(phi_field, stencil, density_heavy, density_light):
stencil: stencil of the lattice Boltzmann step
density_heavy: density of the heavier fluid
density_light: density of the lighter fluid
+ fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
+ field. If it is not given the stencil of the LB method will be applied.
"""
- iso_grad = isotropic_gradient_symbolic(phi_field, stencil)
+ if fd_stencil is None:
+ fd_stencil = stencil
+
+ iso_grad = isotropic_gradient_symbolic(phi_field, fd_stencil)
result = list(map(lambda x: sp.Rational(-1, 3) * sp.symbols("rho") * (density_heavy - density_light) * x, iso_grad))
return result
-def viscous_force(lb_velocity_field, phi_field, mrt_method, tau, density_heavy, density_light):
+def viscous_force(lb_velocity_field, phi_field, mrt_method, tau, density_heavy, density_light, fd_stencil=None):
r"""
Get a symbolic expression for the viscous force
Args:
@@ -117,11 +142,16 @@ def viscous_force(lb_velocity_field, phi_field, mrt_method, tau, density_heavy,
tau: relaxation time of the hydrodynamic lattice boltzmann step
density_heavy: density of the heavier fluid
density_light: density of the lighter fluid
+ fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
+ field. If it is not given the stencil of the LB method will be applied.
"""
stencil = mrt_method.stencil
dimensions = len(stencil[0])
- iso_grad = isotropic_gradient_symbolic(phi_field, stencil)
+ if fd_stencil is None:
+ fd_stencil = stencil
+
+ iso_grad = isotropic_gradient_symbolic(phi_field, fd_stencil)
moment_matrix = mrt_method.moment_matrix
rel = mrt_method.relaxation_rates
@@ -166,7 +196,7 @@ def viscous_force(lb_velocity_field, phi_field, mrt_method, tau, density_heavy,
return [fmx, fmy, fmz]
-def surface_tension_force(phi_field, stencil, beta, kappa):
+def surface_tension_force(phi_field, stencil, beta, kappa, fd_stencil=None):
r"""
Get a symbolic expression for the surface tension force
Args:
@@ -174,14 +204,19 @@ def surface_tension_force(phi_field, stencil, beta, kappa):
stencil: stencil of the lattice Boltzmann step
beta: coefficient related to surface tension and interface thickness
kappa: coefficient related to surface tension and interface thickness
+ fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
+ field. If it is not given the stencil of the LB method will be applied.
"""
- chemical_potential = chemical_potential_symbolic(phi_field, stencil, beta, kappa)
- iso_grad = isotropic_gradient_symbolic(phi_field, stencil)
+ if fd_stencil is None:
+ fd_stencil = stencil
+
+ chemical_potential = chemical_potential_symbolic(phi_field, fd_stencil, beta, kappa)
+ iso_grad = isotropic_gradient_symbolic(phi_field, fd_stencil)
return [chemical_potential * x for x in iso_grad]
def hydrodynamic_force(lb_velocity_field, phi_field, lb_method, tau,
- density_heavy, density_light, kappa, beta, body_force):
+ density_heavy, density_light, kappa, beta, body_force, fd_stencil=None):
r"""
Get a symbolic expression for the hydrodynamic force
Args:
@@ -194,12 +229,18 @@ def hydrodynamic_force(lb_velocity_field, phi_field, lb_method, tau,
beta: coefficient related to surface tension and interface thickness
kappa: coefficient related to surface tension and interface thickness
body_force: force acting on the fluids. Usually the gravity
+ fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
+ field. If it is not given the stencil of the LB method will be applied.
"""
stencil = lb_method.stencil
dimensions = len(stencil[0])
- fp = pressure_force(phi_field, stencil, density_heavy, density_light)
- fm = viscous_force(lb_velocity_field, phi_field, lb_method, tau, density_heavy, density_light)
- fs = surface_tension_force(phi_field, stencil, beta, kappa)
+
+ if fd_stencil is None:
+ fd_stencil = stencil
+
+ fp = pressure_force(phi_field, stencil, density_heavy, density_light, fd_stencil)
+ fm = viscous_force(lb_velocity_field, phi_field, lb_method, tau, density_heavy, density_light, fd_stencil)
+ fs = surface_tension_force(phi_field, stencil, beta, kappa, fd_stencil)
result = []
for i in range(dimensions):
@@ -208,16 +249,21 @@ def hydrodynamic_force(lb_velocity_field, phi_field, lb_method, tau,
return result
-def interface_tracking_force(phi_field, stencil, interface_thickness):
+def interface_tracking_force(phi_field, stencil, interface_thickness, fd_stencil=None):
r"""
Get a symbolic expression for the hydrodynamic force
Args:
phi_field: phase-field
stencil: stencil of the phase-field distribution lattice Boltzmann step
interface_thickness: interface thickness
+ fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
+ field. If it is not given the stencil of the LB method will be applied.
"""
+ if fd_stencil is None:
+ fd_stencil = stencil
+
dimensions = len(stencil[0])
- normal_fd = normalized_isotropic_gradient_symbolic(phi_field, stencil)
+ normal_fd = normalized_isotropic_gradient_symbolic(phi_field, stencil, fd_stencil)
result = []
for i in range(dimensions):
result.append(((1.0 - 4.0 * (phi_field.center - 0.5) ** 2) / interface_thickness) * normal_fd[i])
@@ -264,15 +310,16 @@ def get_update_rules_velocity(src_field, u_in, lb_method, force, density):
return update_u
-def get_collision_assignments_hydro(collision_rule=None, density=1, optimization={}, **kwargs):
+def get_collision_assignments_hydro(density=1, optimization=None, **kwargs):
r"""
Get collision assignments for the hydrodynamic lattice Boltzmann step. Here the force gets applied in the moment
space. Afterwards the transformation back to the pdf space happens.
Args:
- collision_rule: arbitrary collision rule, e.g. created with create_lb_collision_rule
density: the interpolated density of the simulation
optimization: for details see createfunctions.py
"""
+ if optimization is None:
+ optimization = {}
params, opt_params = update_with_default_parameters(kwargs, optimization)
lb_method = params['lb_method']
@@ -288,12 +335,8 @@ def get_collision_assignments_hydro(collision_rule=None, density=1, optimization
if opt_params['symbolic_field'] is not None:
src_field = opt_params['symbolic_field']
- elif opt_params['field_size']:
- field_size = [s + 2 for s in opt_params['field_size']] + [q]
- src_field = Field.create_fixed_size(params['field_name'], field_size, index_dimensions=1,
- layout=opt_params['field_layout'], dtype=field_data_type)
else:
- src_field = Field.create_generic(params['field_name'], spatial_dimensions=collision_rule.method.dim,
+ src_field = Field.create_generic(params['field_name'], spatial_dimensions=lb_method.dim,
index_shape=(q,), layout=opt_params['field_layout'], dtype=field_data_type)
if opt_params['symbolic_temporary_field'] is not None:
@@ -348,26 +391,33 @@ def get_collision_assignments_hydro(collision_rule=None, density=1, optimization
return hydro_lb_update_rule
-def initializer_kernel_phase_field_lb(phi_field_distributions, phi_field, velocity, mrt_method, interface_thickness):
+def initializer_kernel_phase_field_lb(lb_phase_field, phi_field, velocity_field, mrt_method, interface_thickness,
+ fd_stencil=None):
r"""
Returns an assignment list for initializing the phase-field distribution functions
Args:
- phi_field_distributions: source field of phase-field distribution function
+ lb_phase_field: source field of phase-field distribution function
phi_field: phase-field
- velocity: velocity field
+ velocity_field: velocity field
mrt_method: lattice Boltzmann method of the phase-field lattice Boltzmann step
interface_thickness: interface thickness
+ fd_stencil: stencil to derive the finite differences of the isotropic gradient and the laplacian of the phase
+ field. If it is not given the stencil of the LB method will be applied.
"""
stencil = mrt_method.stencil
dimensions = len(stencil[0])
+
+ if fd_stencil is None:
+ fd_stencil = stencil
+
weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
u_symp = sp.symbols("u_:{}".format(dimensions))
- normal_fd = normalized_isotropic_gradient_symbolic(phi_field, stencil)
+ normal_fd = normalized_isotropic_gradient_symbolic(phi_field, stencil, fd_stencil)
gamma = mrt_method.get_equilibrium_terms()
gamma = gamma.subs({sp.symbols("rho"): 1})
- gamma_init = gamma.subs({x: y for x, y in zip(u_symp, velocity.center_vector)})
+ gamma_init = gamma.subs({x: y for x, y in zip(u_symp, velocity_field.center_vector)})
# create the kernels for the initialization of the h field
h_updates = list()
@@ -380,17 +430,17 @@ def initializer_kernel_phase_field_lb(phi_field_distributions, phi_field, veloci
* scalar_product(d, normal_fd[0:dimensions]))
for i, _ in enumerate(stencil):
- h_updates.append(Assignment(phi_field_distributions.center(i), phi_field.center * gamma_init[i] - 0.5 * f[i]))
+ h_updates.append(Assignment(lb_phase_field.center(i), phi_field.center * gamma_init[i] - 0.5 * f[i]))
return h_updates
-def initializer_kernel_hydro_lb(velocity_distributions, velocity, mrt_method):
+def initializer_kernel_hydro_lb(lb_velocity_field, velocity_field, mrt_method):
r"""
Returns an assignment list for initializing the velocity distribution functions
Args:
- velocity_distributions: source field of velocity distribution function
- velocity: velocity field
+ lb_velocity_field: source field of velocity distribution function
+ velocity_field: velocity field
mrt_method: lattice Boltzmann method of the hydrodynamic lattice Boltzmann step
"""
stencil = mrt_method.stencil
@@ -400,10 +450,10 @@ def initializer_kernel_hydro_lb(velocity_distributions, velocity, mrt_method):
gamma = mrt_method.get_equilibrium_terms()
gamma = gamma.subs({sp.symbols("rho"): 1})
- gamma_init = gamma.subs({x: y for x, y in zip(u_symp, velocity.center_vector)})
+ gamma_init = gamma.subs({x: y for x, y in zip(u_symp, velocity_field.center_vector)})
g_updates = list()
for i, _ in enumerate(stencil):
- g_updates.append(Assignment(velocity_distributions.center(i), gamma_init[i] - weights[i]))
+ g_updates.append(Assignment(lb_velocity_field.center(i), gamma_init[i] - weights[i]))
return g_updates
diff --git a/lbmpy_tests/phasefield_allen_cahn/test_codegen_3D.py b/lbmpy_tests/phasefield_allen_cahn/test_codegen_3d.py
similarity index 77%
rename from lbmpy_tests/phasefield_allen_cahn/test_codegen_3D.py
rename to lbmpy_tests/phasefield_allen_cahn/test_codegen_3d.py
index d5b7319c61bc831a699f8907d43a37c9868e2eb3..cf461876169ccc7ac1e5e913b655faae7822a343 100644
--- a/lbmpy_tests/phasefield_allen_cahn/test_codegen_3D.py
+++ b/lbmpy_tests/phasefield_allen_cahn/test_codegen_3d.py
@@ -13,10 +13,9 @@ from lbmpy.phasefield_allen_cahn.parameter_calculation import (
calculate_dimensionless_rising_bubble, calculate_parameters_rti)
from lbmpy.stencils import get_stencil
from pystencils import AssignmentCollection, fields
-from pystencils.simp import sympy_cse
-def test_codegen_3D():
+def test_codegen_3d():
stencil_phase = get_stencil("D3Q15")
stencil_hydro = get_stencil("D3Q27")
assert (len(stencil_phase[0]) == len(stencil_hydro[0]))
@@ -69,7 +68,6 @@ def test_codegen_3D():
# fields
u = fields("vel_field(" + str(dimensions) + "): [" + str(dimensions) + "D]", layout='fzyx')
C = fields("phase_field: [" + str(dimensions) + "D]", layout='fzyx')
- force = fields("force(" + str(dimensions) + "): [" + str(dimensions) + "D]", layout='fzyx')
h = fields("lb_phase_field(" + str(len(stencil_phase)) + "): [" + str(dimensions) + "D]", layout='fzyx')
h_tmp = fields("lb_phase_field_tmp(" + str(len(stencil_phase)) + "): [" + str(dimensions) + "D]", layout='fzyx')
@@ -100,7 +98,8 @@ def test_codegen_3D():
force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W)]
force_model_h = MultiphaseForceModel(force=force_h)
- force_g = hydrodynamic_force(g, C, method_hydro, relaxation_time, density_liquid, density_gas, kappa, beta, body_force)
+ force_g = hydrodynamic_force(g, C, method_hydro,
+ relaxation_time, density_liquid, density_gas, kappa, beta, body_force)
force_model_g = MultiphaseForceModel(force=force_g, rho=density)
h_tmp_symbol_list = [h_tmp.center(i) for i, _ in enumerate(stencil_phase)]
@@ -122,22 +121,28 @@ def test_codegen_3D():
allen_cahn_lb.set_main_assignments_from_dict({**allen_cahn_lb.main_assignments_dict, **{C.center: sum_h}})
allen_cahn_update_rule = AssignmentCollection(main_assignments=allen_cahn_lb.main_assignments,
subexpressions=allen_cahn_lb.subexpressions)
- allen_cahn_update_rule = sympy_cse(allen_cahn_update_rule)
-
# ---------------------------------------------------------------------------------------------------------
method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, force_model=force_model_g)
- hydro_lb_update_rule = get_collision_assignments_hydro(lb_method=method_hydro,
- density=density,
- velocity_input=u,
- force=force_g,
- optimization={"symbolic_field": g,
- "symbolic_temporary_field": g_tmp},
- kernel_type='collide_only')
-
- # streaming of the hydrodynamic distribution
- stream_hydro = create_lb_update_rule(stencil=stencil_hydro,
- optimization={"symbolic_field": g,
- "symbolic_temporary_field": g_tmp},
- kernel_type='stream_pull_only')
+ hydro_lb_update_rule_normal = get_collision_assignments_hydro(lb_method=method_hydro,
+ density=density,
+ velocity_input=u,
+ force=force_g,
+ optimization={"symbolic_field": g,
+ "symbolic_temporary_field": g_tmp},
+ kernel_type='collide_only')
+
+ hydro_lb_update_rule_push = get_collision_assignments_hydro(lb_method=method_hydro,
+ density=density,
+ velocity_input=u,
+ force=force_g,
+ optimization={"symbolic_field": g,
+ "symbolic_temporary_field": g_tmp},
+ kernel_type='collide_stream_push')
+
+ hydro_lb_update_rule_generic_fields = get_collision_assignments_hydro(lb_method=method_hydro,
+ density=density,
+ velocity_input=u,
+ force=force_g,
+ kernel_type='collide_only')
diff --git a/lbmpy_tests/phasefield_allen_cahn/test_phase_field_allen_cahn.ipynb b/lbmpy_tests/phasefield_allen_cahn/test_phase_field_allen_cahn_2d.ipynb
similarity index 100%
rename from lbmpy_tests/phasefield_allen_cahn/test_phase_field_allen_cahn.ipynb
rename to lbmpy_tests/phasefield_allen_cahn/test_phase_field_allen_cahn_2d.ipynb
diff --git a/lbmpy_tests/phasefield_allen_cahn/test_phase_field_derivatives.ipynb b/lbmpy_tests/phasefield_allen_cahn/test_phase_field_derivatives.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..64e03f04c2a69de07547630aefa9240859484bf1
--- /dev/null
+++ b/lbmpy_tests/phasefield_allen_cahn/test_phase_field_derivatives.ipynb
@@ -0,0 +1,323 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from lbmpy.phasefield_allen_cahn.kernel_equations import *\n",
+ "from lbmpy.stencils import get_stencil\n",
+ "from pystencils import fields\n",
+ "\n",
+ "from lbmpy.creationfunctions import create_lb_method"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def apply(field_access: Field.Access, stencil, weights):\n",
+ " f = field_access\n",
+ " return sum(f.get_shifted(*offset) * weight for offset, weight in zip(stencil, weights))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# test chemical potencial"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "stencil = get_stencil(\"D2Q9\")\n",
+ "dimensions = len(stencil[0])\n",
+ "C = fields(\"phase_field: [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "\n",
+ "beta = 0\n",
+ "kappa = 1\n",
+ "\n",
+ "a = chemical_potential_symbolic(C, stencil, beta, kappa)\n",
+ "\n",
+ "expected_result = sp.Array([20, -4, -4, -4, -4, -1, -1, -1, -1]) / 6\n",
+ "b = apply(C.center, stencil, expected_result)\n",
+ "\n",
+ "assert a == b"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "stencil = get_stencil(\"D3Q15\")\n",
+ "dimensions = len(stencil[0])\n",
+ "C = fields(\"phase_field: [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "\n",
+ "beta = 0\n",
+ "kappa = 1\n",
+ "\n",
+ "a = chemical_potential_symbolic(C, stencil, beta, kappa)\n",
+ "\n",
+ "expected_result = sp.Array([256, -28, -28, -28, -28, -28, -28, -11, -11, -11, -11, -11, -11, -11, -11]) / 72\n",
+ "b = apply(C.center, stencil, expected_result)\n",
+ "\n",
+ "assert a == b"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "stencil = get_stencil(\"D3Q19\")\n",
+ "dimensions = len(stencil[0])\n",
+ "C = fields(\"phase_field: [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "\n",
+ "beta = 0\n",
+ "kappa = 1\n",
+ "\n",
+ "a = chemical_potential_symbolic(C, stencil, beta, kappa)\n",
+ "\n",
+ "expected_result = sp.Array([24, -2, -2, -2, -2, -2, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]) / 6\n",
+ "b = apply(C.center, stencil, expected_result)\n",
+ "\n",
+ "assert a == b"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "stencil = get_stencil(\"D3Q27\")\n",
+ "dimensions = len(stencil[0])\n",
+ "C = fields(\"phase_field: [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "\n",
+ "beta = 0\n",
+ "kappa = 1\n",
+ "\n",
+ "a = chemical_potential_symbolic(C, stencil, beta, kappa)\n",
+ "\n",
+ "expected_result = sp.Array([152,\n",
+ " -16, -16, -16, -16, -16, -16,\n",
+ " -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4,\n",
+ " -1, -1, -1, -1, -1, -1, -1, -1]) / 36\n",
+ "b = apply(C.center, stencil, expected_result)\n",
+ "\n",
+ "assert a == b"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# test isotropic gradient"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "stencil = get_stencil(\"D2Q9\")\n",
+ "dimensions = len(stencil[0])\n",
+ "C = fields(\"phase_field: [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "\n",
+ "a = isotropic_gradient_symbolic(C, stencil)\n",
+ "\n",
+ "expected_result = sp.Array([-1, -4, -1, 1, 4, 1]) / 12\n",
+ "expected_grad_stencil = ((-1,-1), (-1,0), (-1,1), (1,-1), (1,0), (1,1))\n",
+ "b = apply(C.center, expected_grad_stencil, expected_result)\n",
+ "\n",
+ "assert a[0] == b"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "stencil = get_stencil(\"D3Q15\")\n",
+ "dimensions = len(stencil[0])\n",
+ "C = fields(\"phase_field: [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "\n",
+ "a = isotropic_gradient_symbolic(C, stencil)\n",
+ "\n",
+ "expected_result = sp.Array([-1, -1, -8, -1, -1, 1, 1, 8, 1, 1]) / 24\n",
+ "expected_grad_stencil = ((-1,-1,-1), (-1,-1,1), (-1,0,0), (-1,1,-1), (-1,1,1), (1,-1,-1), (1,-1,1), (1,0,0), (1,1,-1), (1,1,1))\n",
+ "b = apply(C.center, expected_grad_stencil, expected_result)\n",
+ "\n",
+ "assert a[0] == b"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "stencil = get_stencil(\"D3Q19\")\n",
+ "dimensions = len(stencil[0])\n",
+ "C = fields(\"phase_field: [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "\n",
+ "a = isotropic_gradient_symbolic(C, stencil)\n",
+ "\n",
+ "expected_result = sp.Array([-1, -1, -2, -1, -1, 1, 1, 2, 1, 1]) / 12\n",
+ "expected_grad_stencil = ((-1,-1,0), (-1,0,-1), (-1,0,0), (-1,0,1), (-1,1,0), (1,-1,0), (1,0,-1), (1,0,0), (1,0,1), (1,1,0))\n",
+ "b = apply(C.center, expected_grad_stencil, expected_result)\n",
+ "\n",
+ "assert a[0] == b"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "stencil = get_stencil(\"D3Q27\")\n",
+ "dimensions = len(stencil[0])\n",
+ "C = fields(\"phase_field: [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "\n",
+ "a = isotropic_gradient_symbolic(C, stencil)\n",
+ "\n",
+ "expected_result = sp.Array([-1, -4, -1, -4, -16, -4, -1, -4, -1, 1, 4, 1, 4, 16, 4, 1, 4, 1]) / 72\n",
+ "expected_grad_stencil = ((-1,-1,-1), (-1,-1,0), (-1,-1,1), (-1,0,-1), (-1,0,0), (-1,0,1), (-1,1,-1), (-1,1,0), (-1,1,1),\n",
+ " (1,-1,-1), (1,-1,0), (1,-1,1), (1,0,-1), (1,0,0), (1,0,1), (1,1,-1), (1,1,0), (1,1,1))\n",
+ "b = apply(C.center, expected_grad_stencil, expected_result)\n",
+ "\n",
+ "assert a[0] == b"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# test hydrodynamic force"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "stencil = get_stencil(\"D3Q27\")\n",
+ "dimensions = len(stencil[0])\n",
+ "C = fields(\"phase_field: [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "g = fields(\"lb_velocity_field(\" + str(len(stencil)) + \"): [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "\n",
+ "tau = 0.53\n",
+ "\n",
+ "lb_method = create_lb_method(stencil=stencil, method=\"mrt\", relaxation_rates=[1/tau, 1, 1, 1, 1, 1])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "a = hydrodynamic_force(g, C, lb_method, tau, 1, 0.1, 1, 0, [0, 0, 0] , fd_stencil=None)\n",
+ "b = hydrodynamic_force(g, C, lb_method, tau, 1, 0.1, 1, 0, [0, 0, 0] , fd_stencil=get_stencil(\"D3Q27\"))\n",
+ "c = hydrodynamic_force(g, C, lb_method, tau, 1, 0.1, 1, 0, [0, 0, 0] , fd_stencil=get_stencil(\"D3Q19\"))\n",
+ "d = hydrodynamic_force(g, C, lb_method, tau, 1, 0.1, 1, 0, [0, 0, 0] , fd_stencil=get_stencil(\"D3Q15\"))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "stencil = get_stencil(\"D2Q9\")\n",
+ "dimensions = len(stencil[0])\n",
+ "C = fields(\"phase_field: [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "g = fields(\"lb_velocity_field(\" + str(len(stencil)) + \"): [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "\n",
+ "tau = 0.53\n",
+ "\n",
+ "lb_method = create_lb_method(stencil=stencil, method=\"mrt\", relaxation_rates=[1/tau, 1, 1, 1, 1, 1])\n",
+ "\n",
+ "a = hydrodynamic_force(g, C, lb_method, tau, 1, 0.1, 1, 0, [0, 0, 0] , fd_stencil=None)\n",
+ "b = hydrodynamic_force(g, C, lb_method, tau, 1, 0.1, 1, 0, [0, 0, 0] , fd_stencil=get_stencil(\"D2Q9\"))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "stencil = get_stencil(\"D3Q27\")\n",
+ "dimensions = len(stencil[0])\n",
+ "C = fields(\"phase_field: [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "u = fields(\"vel_field(\" + str(dimensions) + \"): [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "h = fields(\"lb_velocity_field(\" + str(len(stencil)) + \"): [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "\n",
+ "tau = 0.53\n",
+ "\n",
+ "lb_method = create_lb_method(stencil=stencil, method=\"srt\")\n",
+ "\n",
+ "a = initializer_kernel_phase_field_lb(h, C, u, lb_method, 5, fd_stencil=None)\n",
+ "b = initializer_kernel_phase_field_lb(h, C, u, lb_method, 5, fd_stencil=get_stencil(\"D3Q27\"))\n",
+ "c = initializer_kernel_phase_field_lb(h, C, u, lb_method, 5, fd_stencil=get_stencil(\"D3Q19\"))\n",
+ "d = initializer_kernel_phase_field_lb(h, C, u, lb_method, 5, fd_stencil=get_stencil(\"D3Q15\"))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "stencil = get_stencil(\"D2Q9\")\n",
+ "dimensions = len(stencil[0])\n",
+ "C = fields(\"phase_field: [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "u = fields(\"vel_field(\" + str(dimensions) + \"): [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "h = fields(\"lb_velocity_field(\" + str(len(stencil)) + \"): [\" + str(dimensions) + \"D]\", layout='fzyx')\n",
+ "\n",
+ "tau = 0.53\n",
+ "\n",
+ "lb_method = create_lb_method(stencil=stencil, method=\"srt\")\n",
+ "\n",
+ "a = initializer_kernel_phase_field_lb(h, C, u, lb_method, 5, fd_stencil=None)\n",
+ "b = initializer_kernel_phase_field_lb(h, C, u, lb_method, 5, fd_stencil=get_stencil(\"D2Q9\"))"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.7.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}