Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • ravi.k.ayyala/lbmpy
  • brendan-waters/lbmpy
  • anirudh.jonnalagadda/lbmpy
  • jbadwaik/lbmpy
  • alexander.reinauer/lbmpy
  • itischler/lbmpy
  • he66coqe/lbmpy
  • ev81oxyl/lbmpy
  • Bindgen/lbmpy
  • da15siwa/lbmpy
  • holzer/lbmpy
  • RudolfWeeber/lbmpy
  • pycodegen/lbmpy
13 results
Select Git revision
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 4500 deletions
lbmpy/methods/momentbased/momentbasedmethod.py deleted 100644 → 0
View file @ cf17a9b9
from collections import OrderedDict
import sympy as sp
from lbmpy.maxwellian_equilibrium import get_weights
from lbmpy.methods.abstractlbmethod import AbstractLbMethod, LbmCollisionRule, RelaxationInfo
from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantityComputation
from lbmpy.moments import MOMENT_SYMBOLS, moment_matrix
from pystencils import Assignment
from pystencils.sympyextensions import subs_additive
class MomentBasedLbMethod(AbstractLbMethod):
def __init__(self, stencil, moment_to_relaxation_info_dict, conserved_quantity_computation=None, force_model=None):
"""
Moment based LBM is a class to represent the single (SRT), two (TRT) and multi relaxation time (MRT) methods.
These methods work by transforming the pdfs into moment space using a linear transformation. In the moment
space each component (moment) is relaxed to an equilibrium moment by a certain relaxation rate. These
equilibrium moments can e.g. be determined by taking the equilibrium moments of the continuous Maxwellian.
Args:
stencil: see :func:`lbmpy.stencils.get_stencil`
moment_to_relaxation_info_dict: a dictionary mapping moments in either tuple or polynomial formulation
to a RelaxationInfo, which consists of the corresponding equilibrium moment
and a relaxation rate
conserved_quantity_computation: instance of :class:`lbmpy.methods.AbstractConservedQuantityComputation`.
This determines how conserved quantities are computed, and defines
the symbols used in the equilibrium moments like e.g. density and velocity
force_model: force model instance, or None if no forcing terms are required
"""
assert isinstance(conserved_quantity_computation, AbstractConservedQuantityComputation)
super(MomentBasedLbMethod, self).__init__(stencil)
self._forceModel = force_model
self._momentToRelaxationInfoDict = OrderedDict(moment_to_relaxation_info_dict.items())
self._conservedQuantityComputation = conserved_quantity_computation
self._weights = None
@property
def force_model(self):
return self._forceModel
@property
def relaxation_info_dict(self):
return self._momentToRelaxationInfoDict
@property
def conserved_quantity_computation(self):
return self._conservedQuantityComputation
@property
def moments(self):
return tuple(self._momentToRelaxationInfoDict.keys())
@property
def moment_equilibrium_values(self):
return tuple([e.equilibrium_value for e in self._momentToRelaxationInfoDict.values()])
@property
def relaxation_rates(self):
return tuple([e.relaxation_rate for e in self._momentToRelaxationInfoDict.values()])
@property
def zeroth_order_equilibrium_moment_symbol(self, ):
return self._conservedQuantityComputation.zeroth_order_moment_symbol
@property
def first_order_equilibrium_moment_symbols(self, ):
return self._conservedQuantityComputation.first_order_moment_symbols
@property
def weights(self):
if self._weights is None:
self._weights = self._compute_weights()
return self._weights
def override_weights(self, weights):
assert len(weights) == len(self.stencil)
self._weights = weights
def get_equilibrium(self, conserved_quantity_equations=None, include_force_terms=False):
relaxation_matrix = sp.eye(len(self.relaxation_rates))
return self._collision_rule_with_relaxation_matrix(relaxation_matrix,
conserved_quantity_equations=conserved_quantity_equations,
include_force_terms=include_force_terms)
def get_equilibrium_terms(self):
equilibrium = self.get_equilibrium()
return sp.Matrix([eq.rhs for eq in equilibrium.main_assignments])
def get_collision_rule(self, conserved_quantity_equations=None, pre_simplification=True):
d = self.relaxation_matrix
relaxation_rate_sub_expressions, d = self._generate_relaxation_matrix(d, pre_simplification)
ac = self._collision_rule_with_relaxation_matrix(d, relaxation_rate_sub_expressions,
True, conserved_quantity_equations)
return ac
def set_zeroth_moment_relaxation_rate(self, relaxation_rate):
e = sp.Rational(1, 1)
prev_entry = self._momentToRelaxationInfoDict[e]
new_entry = RelaxationInfo(prev_entry[0], relaxation_rate)
self._momentToRelaxationInfoDict[e] = new_entry
def set_first_moment_relaxation_rate(self, relaxation_rate):
for e in MOMENT_SYMBOLS[:self.dim]:
assert e in self._momentToRelaxationInfoDict, "First moments are not relaxed separately by this method"
for e in MOMENT_SYMBOLS[:self.dim]:
prev_entry = self._momentToRelaxationInfoDict[e]
new_entry = RelaxationInfo(prev_entry[0], relaxation_rate)
self._momentToRelaxationInfoDict[e] = new_entry
def set_conserved_moments_relaxation_rate(self, relaxation_rate):
self.set_zeroth_moment_relaxation_rate(relaxation_rate)
self.set_first_moment_relaxation_rate(relaxation_rate)
def set_force_model(self, force_model):
self._forceModel = force_model
@property
def collision_matrix(self):
pdfs_to_moments = self.moment_matrix
d = self.relaxation_matrix
return pdfs_to_moments.inv() * d * pdfs_to_moments
@property
def inverse_collision_matrix(self):
pdfs_to_moments = self.moment_matrix
inverse_relaxation_matrix = self.relaxation_matrix.inv()
return pdfs_to_moments.inv() * inverse_relaxation_matrix * pdfs_to_moments
@property
def moment_matrix(self):
return moment_matrix(self.moments, self.stencil)
@property
def relaxation_matrix(self):
d = sp.zeros(len(self.relaxation_rates))
for i in range(0, len(self.relaxation_rates)):
d[i, i] = self.relaxation_rates[i]
return d
@property
def is_orthogonal(self):
return (self.moment_matrix * self.moment_matrix.T).is_diagonal()
@property
def is_weighted_orthogonal(self):
w = get_weights(self.stencil, sp.Rational(1, 3))
return (sp.matrix_multiply_elementwise(self.moment_matrix, sp.Matrix([w] * len(w))) * self.moment_matrix.T
).is_diagonal()
def __getstate__(self):
# Workaround for a bug in joblib
self._momentToRelaxationInfoDictToPickle = [i for i in self._momentToRelaxationInfoDict.items()]
return self.__dict__
def _repr_html_(self):
table = """
<table style="border:none; width: 100%">
<tr {nb}>
<th {nb} >Moment</th>
<th {nb} >Eq. Value </th>
<th {nb} >Relaxation Rate</th>
</tr>
{content}
</table>
"""
content = ""
for moment, (eq_value, rr) in self._momentToRelaxationInfoDict.items():
vals = {
'rr': sp.latex(rr),
'moment': sp.latex(moment),
'eq_value': sp.latex(eq_value),
'nb': 'style="border:none"',
}
content += """<tr {nb}>
<td {nb}>${moment}$</td>
<td {nb}>${eq_value}$</td>
<td {nb}>${rr}$</td>
</tr>\n""".format(**vals)
return table.format(content=content, nb='style="border:none"')
def _compute_weights(self):
replacements = self._conservedQuantityComputation.default_values
ac = self.get_equilibrium(include_force_terms=False)
ac = ac.new_with_substitutions(replacements, substitute_on_lhs=False).new_without_subexpressions()
new_assignments = [Assignment(e.lhs,
subs_additive(e.rhs, sp.sympify(1), sum(self.pre_collision_pdf_symbols),
required_match_replacement=1.0))
for e in ac.main_assignments]
ac = ac.copy(new_assignments)
weights = []
for eq in ac.main_assignments:
value = eq.rhs.expand()
assert len(value.atoms(sp.Symbol)) == 0, "Failed to compute weights " + str(value)
weights.append(value)
return weights
def _collision_rule_with_relaxation_matrix(self, d, additional_subexpressions=(), include_force_terms=True,
conserved_quantity_equations=None):
f = sp.Matrix(self.pre_collision_pdf_symbols)
pdf_to_moment = self.moment_matrix
m_eq = sp.Matrix(self.moment_equilibrium_values)
collision_rule = f + pdf_to_moment.inv() * d * (m_eq - pdf_to_moment * f)
collision_eqs = [Assignment(lhs, rhs) for lhs, rhs in zip(self.post_collision_pdf_symbols, collision_rule)]
if conserved_quantity_equations is None:
conserved_quantity_equations = self._conservedQuantityComputation.equilibrium_input_equations_from_pdfs(f)
simplification_hints = conserved_quantity_equations.simplification_hints.copy()
simplification_hints.update(self._conservedQuantityComputation.defined_symbols())
simplification_hints['relaxation_rates'] = [d[i, i] for i in range(d.rows)]
all_subexpressions = list(additional_subexpressions) + conserved_quantity_equations.all_assignments
if self._forceModel is not None and include_force_terms:
force_model_terms = self._forceModel(self)
force_term_symbols = sp.symbols("forceTerm_:%d" % (len(force_model_terms,)))
force_subexpressions = [Assignment(sym, force_model_term)
for sym, force_model_term in zip(force_term_symbols, force_model_terms)]
all_subexpressions += force_subexpressions
collision_eqs = [Assignment(eq.lhs, eq.rhs + force_term_symbol)
for eq, force_term_symbol in zip(collision_eqs, force_term_symbols)]
simplification_hints['force_terms'] = force_term_symbols
return LbmCollisionRule(self, collision_eqs, all_subexpressions,
simplification_hints)
@staticmethod
def _generate_relaxation_matrix(relaxation_matrix, keep_rr_symbolic):
"""
For SRT and TRT the equations can be easier simplified if the relaxation times are symbols, not numbers.
This function replaces the numbers in the relaxation matrix with symbols in this case, and returns also
the subexpressions, that assign the number to the newly introduced symbol
"""
rr = [relaxation_matrix[i, i] for i in range(relaxation_matrix.rows)]
if keep_rr_symbolic <= 2:
unique_relaxation_rates = set(rr)
subexpressions = {}
for rt in unique_relaxation_rates:
rt = sp.sympify(rt)
if not isinstance(rt, sp.Symbol):
rt_symbol = sp.Symbol("rr_%d" % (len(subexpressions),))
subexpressions[rt] = rt_symbol
new_rr = [subexpressions[sp.sympify(e)] if sp.sympify(e) in subexpressions else e
for e in rr]
substitutions = [Assignment(e[1], e[0]) for e in subexpressions.items()]
return substitutions, sp.diag(*new_rr)
else:
return [], relaxation_matrix
lbmpy/phasefield_allen_cahn/analytical.py deleted 100644 → 0
View file @ cf17a9b9
def analytic_rising_speed(gravitational_acceleration, bubble_diameter, viscosity_gas):
r"""
Calculated the analytical rising speed of a bubble. This is the expected end rising speed.
Args:
gravitational_acceleration: the gravitational acceleration acting in the simulation scenario. Usually it gets
calculated based on dimensionless parameters which describe the scenario
bubble_diameter: the diameter of the bubble at the beginning of the simulation
viscosity_gas: the viscosity of the fluid inside the bubble
"""
result = -(gravitational_acceleration * bubble_diameter * bubble_diameter) / (12.0 * viscosity_gas)
return result
lbmpy/phasefield_allen_cahn/force_model.py deleted 100644 → 0
View file @ cf17a9b9
import sympy as sp
import numpy as np
from lbmpy.forcemodels import Simple
class MultiphaseForceModel:
r"""
A force model based on PhysRevE.96.053301. This model realises the modified equilibrium distributions meaning the
force gets shifted by minus one half multiplied with the collision operator
"""
def __init__(self, force, rho=1):
self._force = force
self._rho = rho
def __call__(self, lb_method):
stencil = lb_method.stencil
force_symp = sp.symbols("force_:{}".format(lb_method.dim))
simple = Simple(force_symp)
force = [f / self._rho for f in simple(lb_method)]
moment_matrix = lb_method.moment_matrix
relaxation_rates = sp.Matrix(np.diag(lb_method.relaxation_rates))
mrt_collision_op = moment_matrix.inv() * relaxation_rates * moment_matrix
result = -0.5 * mrt_collision_op * sp.Matrix(force) + sp.Matrix(force)
for i in range(0, len(stencil)):
result[i] = result[i].simplify()
subs_dict = dict(zip(force_symp, self._force))
for i in range(0, len(stencil)):
result[i] = result[i].subs(subs_dict)
return result
lbmpy/phasefield_allen_cahn/kernel_equations.py deleted 100644 → 0
View file @ cf17a9b9
from pystencils.fd.derivation import FiniteDifferenceStencilDerivation
from pystencils import Assignment, AssignmentCollection
from lbmpy.maxwellian_equilibrium import get_weights
from lbmpy.fieldaccess import StreamPushTwoFieldsAccessor, CollideOnlyInplaceAccessor
import sympy as sp
import numpy as np
def chemical_potential_symbolic(phi_field, stencil, beta, kappa):
r"""
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 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])
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) - \
kappa * lap
return mu
def isotropic_gradient_symbolic(phi_field, stencil):
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
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)
deriv.assume_symmetric(1, anti_symmetric=False)
if dimensions == 3:
deriv.assume_symmetric(2, anti_symmetric=False)
# 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 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))
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))]
else:
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, 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.
"""
if fd_stencil is None:
fd_stencil = stencil
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, fd_stencil=None):
r"""
Get a symbolic expression for the pressure force
Args:
phi_field: phase-field
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.
"""
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, fd_stencil=None):
r"""
Get a symbolic expression for the viscous force
Args:
lb_velocity_field: hydrodynamic distribution function
phi_field: phase-field
mrt_method: mrt lattice boltzmann method used for hydrodynamics
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])
if fd_stencil is None:
fd_stencil = stencil
iso_grad = isotropic_gradient_symbolic(phi_field, fd_stencil)
non_equilibrium = lb_velocity_field.center_vector - mrt_method.get_equilibrium_terms()
stress_tensor = [0] * 6
# Calculate Stress Tensor MRT
for i, d in enumerate(stencil):
stress_tensor[0] = sp.Add(stress_tensor[0], non_equilibrium[i] * (d[0] * d[0]))
stress_tensor[1] = sp.Add(stress_tensor[1], non_equilibrium[i] * (d[1] * d[1]))
if dimensions == 3:
stress_tensor[2] = sp.Add(stress_tensor[2], non_equilibrium[i] * (d[2] * d[2]))
stress_tensor[3] = sp.Add(stress_tensor[3], non_equilibrium[i] * (d[1] * d[2]))
stress_tensor[4] = sp.Add(stress_tensor[4], non_equilibrium[i] * (d[0] * d[2]))
stress_tensor[5] = sp.Add(stress_tensor[5], non_equilibrium[i] * (d[0] * d[1]))
density_difference = density_heavy - density_light
# Calculate Viscous Force MRT
fmx = (0.5 - tau) * (stress_tensor[0] * iso_grad[0]
+ stress_tensor[5] * iso_grad[1]
+ stress_tensor[4] * iso_grad[2]) * density_difference
fmy = (0.5 - tau) * (stress_tensor[5] * iso_grad[0]
+ stress_tensor[1] * iso_grad[1]
+ stress_tensor[3] * iso_grad[2]) * density_difference
fmz = (0.5 - tau) * (stress_tensor[4] * iso_grad[0]
+ stress_tensor[3] * iso_grad[1]
+ stress_tensor[2] * iso_grad[2]) * density_difference
return [fmx, fmy, fmz]
def surface_tension_force(phi_field, stencil, beta, kappa, fd_stencil=None):
r"""
Get a symbolic expression for the surface tension force
Args:
phi_field: the phase-field on which the chemical potential is applied
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.
"""
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, fd_stencil=None):
r"""
Get a symbolic expression for the hydrodynamic force
Args:
lb_velocity_field: hydrodynamic distribution function
phi_field: phase-field
lb_method: Lattice boltzmann method used for hydrodynamics
tau: relaxation time of the hydrodynamic lattice boltzmann step
density_heavy: density of the heavier fluid
density_light: density of the lighter fluid
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])
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):
result.append(fs[i] + fp[i] + fm[i] + body_force[i])
return result
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, 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])
return result
def get_update_rules_velocity(src_field, u_in, lb_method, force, density, sub_iterations=2):
r"""
Get assignments to update the velocity with a force shift
Args:
src_field: the source field of the hydrodynamic distribution function
u_in: velocity field
lb_method: mrt lattice boltzmann method used for hydrodynamics
force: force acting on the hydrodynamic lb step
density: the interpolated density of the simulation
sub_iterations: number of updates of the velocity field
"""
stencil = lb_method.stencil
dimensions = len(stencil[0])
moment_matrix = lb_method.moment_matrix
eq = lb_method.moment_equilibrium_values
first_eqs = lb_method.first_order_equilibrium_moment_symbols
indices = list()
for i in range(dimensions):
indices.append(eq.index(first_eqs[i]))
src = [src_field.center(i) for i, _ in enumerate(stencil)]
m0 = np.dot(moment_matrix.tolist(), src)
update_u = list()
update_u.append(Assignment(sp.symbols("rho"), m0[0]))
index = 0
u_symp = sp.symbols("u_:{}".format(dimensions))
aleph = sp.symbols("aleph_:{}".format(dimensions * sub_iterations))
for i in range(dimensions):
update_u.append(Assignment(aleph[i], u_in.center_vector[i]))
index += 1
for k in range(sub_iterations - 1):
subs_dict = dict(zip(u_symp, aleph[k * dimensions:index]))
for i in range(dimensions):
update_u.append(Assignment(aleph[index], m0[indices[i]] + force[i].subs(subs_dict) / density / 2))
index += 1
subs_dict = dict(zip(u_symp, aleph[index - dimensions:index]))
for i in range(dimensions):
update_u.append(Assignment(u_symp[i], m0[indices[i]] + force[i].subs(subs_dict) / density / 2))
# update_u.append(Assignment(u_in.center_vector[i], m0[indices[i]] + force[i].subs(subs_dict) / density / 2))
return update_u
def get_collision_assignments_hydro(lb_method, density, velocity_input, force, sub_iterations, symbolic_fields,
kernel_type):
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:
lb_method: moment based lattice Boltzmann method
density: the interpolated density of the simulation
velocity_input: velocity field for the hydrodynamic and Allen-Chan LB step
force: force vector containing a summation of the surface tension-, pressure-, viscous- and bodyforce vector
sub_iterations: number of updates of the velocity field
symbolic_fields: PDF fields for source and destination
kernel_type: collide_stream_push or collide_only
"""
stencil = lb_method.stencil
dimensions = len(stencil[0])
src_field = symbolic_fields['symbolic_field']
dst_field = symbolic_fields['symbolic_temporary_field']
moment_matrix = lb_method.moment_matrix
rel = lb_method.relaxation_rates
eq = lb_method.moment_equilibrium_values
first_eqs = lb_method.first_order_equilibrium_moment_symbols
indices = list()
for i in range(dimensions):
indices.append(eq.index(first_eqs[i]))
eq = np.array(eq)
g_vals = [src_field.center(i) for i, _ in enumerate(stencil)]
m0 = np.dot(moment_matrix.tolist(), g_vals)
mf = np.zeros(len(stencil), dtype=object)
for i in range(dimensions):
mf[indices[i]] = force[i] / density
m = sp.symbols("m_:{}".format(len(stencil)))
update_m = get_update_rules_velocity(src_field, velocity_input, lb_method, force,
density, sub_iterations=sub_iterations)
u_symp = sp.symbols("u_:{}".format(dimensions))
for i in range(0, len(stencil)):
update_m.append(Assignment(m[i], m0[i] - (m0[i] - eq[i] + mf[i] / 2) * rel[i] + mf[i]))
update_g = list()
var = np.dot(moment_matrix.inv().tolist(), m)
if kernel_type == 'collide_stream_push':
push_accessor = StreamPushTwoFieldsAccessor()
post_collision_accesses = push_accessor.write(dst_field, stencil)
else:
collide_accessor = CollideOnlyInplaceAccessor()
post_collision_accesses = collide_accessor.write(src_field, stencil)
for i in range(0, len(stencil)):
update_g.append(Assignment(post_collision_accesses[i], var[i]))
for i in range(dimensions):
update_g.append(Assignment(velocity_input.center_vector[i], u_symp[i]))
hydro_lb_update_rule = AssignmentCollection(main_assignments=update_g,
subexpressions=update_m)
return hydro_lb_update_rule
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:
lb_phase_field: source field of phase-field distribution function
phi_field: phase-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, 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_field.center_vector)})
# create the kernels for the initialization of the h field
h_updates = list()
def scalar_product(a, b):
return sum(a_i * b_i for a_i, b_i in zip(a, b))
f = []
for i, d in enumerate(stencil):
f.append(weights[i] * ((1.0 - 4.0 * (phi_field.center - 0.5) ** 2) / interface_thickness)
* scalar_product(d, normal_fd[0:dimensions]))
for i, _ in enumerate(stencil):
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(lb_velocity_field, velocity_field, mrt_method):
r"""
Returns an assignment list for initializing the velocity distribution functions
Args:
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
dimensions = len(stencil[0])
weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
u_symp = sp.symbols("u_:{}".format(dimensions))
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_field.center_vector)})
g_updates = list()
for i, _ in enumerate(stencil):
g_updates.append(Assignment(lb_velocity_field.center(i), gamma_init[i] - weights[i]))
return g_updates
lbmpy/phasefield_allen_cahn/parameter_calculation.py deleted 100644 → 0
View file @ cf17a9b9
import math
def calculate_parameters_rti(reference_length=256,
reference_time=16000,
density_heavy=1.0,
capillary_number=0.26,
reynolds_number=3000,
atwood_number=0.5,
peclet_number=1000,
density_ratio=3,
viscosity_ratio=1):
r"""
Calculate the simulation parameters for the Rayleigh-Taylor instability.
The calculation is done according to the description in part B of PhysRevE.96.053301.
Args:
reference_length: reference length of the RTI
reference_time: chosen reference time
density_heavy: density of the heavier fluid
capillary_number: capillary number of the simulation represents the relative effect of viscous drag
forces versus surface tension forces
reynolds_number: reynolds number of the simulation is the ratio between the viscous forces in a fluid
and the inertial forces
atwood_number: atwood number of the simulation is a dimensionless density ratio
peclet_number: peclet number of the simulation is the ratio of the rate of advection
of a physical quantity by the flow to the rate of diffusion of the same quantity
driven by an appropriate gradient
density_ratio: density ratio of the heavier and the lighter fluid
viscosity_ratio: viscosity ratio of the heavier and the lighter fluid
"""
# calculate the gravitational acceleration and the reference velocity
g = reference_length / (reference_time ** 2 * atwood_number)
reference_velocity = math.sqrt(g * reference_length)
dynamic_viscosity_heavy = (density_heavy * reference_velocity * reference_length) / reynolds_number
dynamic_viscosity_light = dynamic_viscosity_heavy / viscosity_ratio
density_light = density_heavy / density_ratio
kinematic_viscosity_heavy = dynamic_viscosity_heavy / density_heavy
kinematic_viscosity_light = dynamic_viscosity_light / density_light
relaxation_time_heavy = 3.0 * kinematic_viscosity_heavy
relaxation_time_light = 3.0 * kinematic_viscosity_light
surface_tension = (dynamic_viscosity_heavy * reference_velocity) / capillary_number
mobility = (reference_velocity * reference_length) / peclet_number
parameters = {
"density_light": density_light,
"dynamic_viscosity_heavy": dynamic_viscosity_heavy,
"dynamic_viscosity_light": dynamic_viscosity_light,
"relaxation_time_heavy": relaxation_time_heavy,
"relaxation_time_light": relaxation_time_light,
"gravitational_acceleration": -g,
"mobility": mobility,
"surface_tension": surface_tension
}
return parameters
def calculate_dimensionless_rising_bubble(reference_time=18000,
density_heavy=1.0,
bubble_radius=16,
bond_number=1,
reynolds_number=40,
density_ratio=1000,
viscosity_ratio=100):
r"""
Calculate the simulation parameters for a rising bubble. The parameter calculation follows the ideas of Weber and
is based on the Reynolds number. This means the rising velocity of the bubble is implicitly stated with the
Reynolds number
Args:
reference_time: chosen reference time
density_heavy: density of the heavier fluid
bubble_radius: initial radius of the rising bubble
bond_number: also called eötvös number is measuring the importance of gravitational forces compared to
surface tension forces
reynolds_number: reynolds number of the simulation is the ratio between the viscous forces in a fluid
and the inertial forces
density_ratio: density ratio of the heavier and the lighter fluid
viscosity_ratio: viscosity ratio of the heavier and the lighter fluid
"""
bubble_diameter = bubble_radius * 2
g = bubble_diameter / (reference_time ** 2)
density_light = density_heavy / density_ratio
dynamic_viscosity_heavy = (density_heavy * math.sqrt(g * bubble_diameter ** 3)) / reynolds_number
dynamic_viscosity_light = dynamic_viscosity_heavy / viscosity_ratio
kinematic_viscosity_heavy = dynamic_viscosity_heavy / density_heavy
kinematic_viscosity_light = dynamic_viscosity_light / density_light
relaxation_time_heavy = 3 * kinematic_viscosity_heavy
relaxation_time_light = 3 * kinematic_viscosity_light
surface_tension = (density_heavy - density_light) * g * bubble_diameter ** 2 / bond_number
# calculation of the Morton number
# Mo = gravitational_acceleration * dynamic_viscosity_heavy / (density_heavy * surface_tension ** 3)
parameters = {
"density_light": density_light,
"dynamic_viscosity_heavy": dynamic_viscosity_heavy,
"dynamic_viscosity_light": dynamic_viscosity_light,
"relaxation_time_heavy": relaxation_time_heavy,
"relaxation_time_light": relaxation_time_light,
"gravitational_acceleration": -g,
"surface_tension": surface_tension
}
return parameters
lbmpy/simplificationfactory.py deleted 100644 → 0
View file @ cf17a9b9
import sympy as sp
from lbmpy.innerloopsplit import create_lbm_split_groups
from lbmpy.methods.momentbased.momentbasedmethod import MomentBasedLbMethod
from lbmpy.methods.centeredcumulant import CenteredCumulantBasedLbMethod
from lbmpy.methods.centeredcumulant.simplification import insert_aliases, insert_zeros, insert_constants
from lbmpy.methods.momentbased.momentbasedsimplifications import (
factor_density_after_factoring_relaxation_times, factor_relaxation_rates,
replace_common_quadratic_and_constant_term, replace_density_and_velocity, replace_second_order_velocity_products)
from pystencils.simp import (
SimplificationStrategy, add_subexpressions_for_divisions, apply_to_all_assignments,
subexpression_substitution_in_main_assignments)
def create_simplification_strategy(lb_method, split_inner_loop=False):
s = SimplificationStrategy()
expand = apply_to_all_assignments(sp.expand)
if isinstance(lb_method, MomentBasedLbMethod):
if len(set(lb_method.relaxation_rates)) <= 2:
s.add(expand)
s.add(replace_second_order_velocity_products)
s.add(expand)
s.add(factor_relaxation_rates)
s.add(replace_density_and_velocity)
s.add(replace_common_quadratic_and_constant_term)
s.add(factor_density_after_factoring_relaxation_times)
s.add(subexpression_substitution_in_main_assignments)
if split_inner_loop:
s.add(create_lbm_split_groups)
s.add(add_subexpressions_for_divisions)
else:
s.add(subexpression_substitution_in_main_assignments)
if split_inner_loop:
s.add(create_lbm_split_groups)
elif isinstance(lb_method, CenteredCumulantBasedLbMethod):
s.add(insert_zeros)
s.add(insert_aliases)
s.add(insert_constants)
s.add(lambda ac: ac.new_without_unused_subexpressions())
return s
lbmpy/stencils.py deleted 100644 → 0
View file @ cf17a9b9
def get_stencil(name, ordering='walberla'):
"""
Stencils are tuples of discrete velocities. They are commonly labeled in the 'DxQy' notation, where d is the
dimension (length of the velocity tuples) and y is number of discrete velocities.
Args:
name: DxQy notation
ordering: the LBM literature does not use a common order of the discrete velocities, therefore here
different common orderings are available. All orderings lead to the same method, it just has
to be used consistently. Here more orderings are available to compare intermediate results with
the literature.
"""
try:
return get_stencil.data[name.upper()][ordering.lower()]
except KeyError:
err_msg = ""
for stencil, ordering_names in get_stencil.data.items():
err_msg += " %s: %s\n" % (stencil, ", ".join(ordering_names.keys()))
raise ValueError("No such stencil available. "
"Available stencils: <stencil_name>( <ordering_names> )\n" + err_msg)
get_stencil.data = {
'D2Q9': {
'walberla': ((0, 0),
(0, 1), (0, -1), (-1, 0), (1, 0),
(-1, 1), (1, 1), (-1, -1), (1, -1),),
'counterclockwise': ((0, 0),
(1, 0), (0, 1), (-1, 0), (0, -1),
(1, 1), (-1, 1), (-1, -1), (1, -1)),
'braunschweig': ((0, 0),
(-1, 1), (-1, 0), (-1, -1), (0, -1),
(1, -1), (1, 0), (1, 1), (0, 1)),
'uk': ((0, 0),
(1, 0), (-1, 0), (0, 1), (0, -1),
(1, 1), (-1, -1), (-1, 1), (1, -1),
)
},
'D2V17': {
'walberla': (
(0, 0), (0, -1), (-1, 0), (1, 0), (0, 1), (-1, -1), (1, -1), (-1, 1), (1, 1), (-2, -2), (2, -2), (-2, 2),
(2, 2), (0, -3), (-3, 0), (3, 0), (0, 3)),
},
'D2V37': {
'walberla': (
(0, 0), (0, -1), (-1, 0), (1, 0), (0, 1), (-1, -1), (1, -1), (-1, 1), (1, 1), (0, -2), (-2, 0), (2, 0),
(0, 2), (-1, -2), (1, -2), (-2, -1), (2, -1), (-2, 1), (2, 1), (-1, 2), (1, 2), (-2, -2), (2, -2), (-2, 2),
(2, 2), (0, -3), (-3, 0), (3, 0), (0, 3), (-1, -3), (1, -3), (-3, -1), (3, -1), (-3, 1), (3, 1), (-1, 3),
(1, 3))
},
'D3Q15': {
'walberla':
((0, 0, 0),
(0, 1, 0), (0, -1, 0), (-1, 0, 0), (1, 0, 0), (0, 0, 1), (0, 0, -1),
(1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1), (1, 1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, -1)),
'premnath': ((0, 0, 0),
(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
(1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1),
(1, 1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, -1)),
},
'D3Q19': {
'walberla': ((0, 0, 0),
(0, 1, 0), (0, -1, 0), (-1, 0, 0), (1, 0, 0), (0, 0, 1), (0, 0, -1),
(-1, 1, 0), (1, 1, 0), (-1, -1, 0), (1, -1, 0),
(0, 1, 1), (0, -1, 1), (-1, 0, 1), (1, 0, 1),
(0, 1, -1), (0, -1, -1), (-1, 0, -1), (1, 0, -1)),
'braunschweig': ((0, 0, 0),
(1, 0, 0), (-1, 0, 0),
(0, 1, 0), (0, -1, 0),
(0, 0, 1), (0, 0, -1),
(1, 1, 0), (-1, -1, 0),
(1, -1, 0), (-1, 1, 0),
(1, 0, 1), (-1, 0, -1),
(1, 0, -1), (-1, 0, 1),
(0, 1, 1), (0, -1, -1),
(0, 1, -1), (0, -1, 1)),
'premnath': ((0, 0, 0),
(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
(1, 1, 0), (-1, 1, 0), (1, -1, 0), (-1, -1, 0),
(1, 0, 1), (-1, 0, 1), (1, 0, -1), (-1, 0, -1),
(0, 1, 1), (0, -1, 1), (0, 1, -1), (0, -1, -1)),
},
'D3Q27': {
'walberla': ((0, 0, 0),
(0, 1, 0), (0, -1, 0), (-1, 0, 0), (1, 0, 0), (0, 0, 1), (0, 0, -1),
(-1, 1, 0), (1, 1, 0), (-1, -1, 0), (1, -1, 0),
(0, 1, 1), (0, -1, 1), (-1, 0, 1), (1, 0, 1),
(0, 1, -1), (0, -1, -1), (-1, 0, -1), (1, 0, -1),
(1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1), (1, 1, -1), (-1, 1, -1), (1, -1, -1),
(-1, -1, -1)),
'premnath': ((0, 0, 0),
(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
(1, 1, 0), (-1, 1, 0), (1, -1, 0), (-1, -1, 0),
(1, 0, 1), (-1, 0, 1), (1, 0, -1), (-1, 0, -1),
(0, 1, 1), (0, -1, 1), (0, 1, -1), (0, -1, -1),
(1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1),
(1, 1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, -1)),
'fakhari': ((0, 0, 0),
(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1),
(1, 1, 1), (-1, 1, 1), (1, -1, 1), (-1, -1, 1),
(1, 1, -1), (-1, 1, -1), (1, -1, -1), (-1, -1, -1),
(1, 1, 0), (-1, 1, 0), (1, -1, 0), (-1, -1, 0),
(1, 0, 1), (-1, 0, 1), (1, 0, -1), (-1, 0, -1), (0, 1, 1), (0, -1, 1), (0, 1, -1), (0, -1, -1)),
}
}
lbmpy/turbulence_models.py deleted 100644 → 0
View file @ cf17a9b9
import sympy as sp
from lbmpy.relaxationrates import get_shear_relaxation_rate
from pystencils import Assignment
def second_order_moment_tensor(function_values, stencil):
"""Returns (D x D) Matrix of second order moments of the given function where D is the dimension"""
assert len(function_values) == len(stencil)
dim = len(stencil[0])
return sp.Matrix(dim, dim, lambda i, j: sum(c[i] * c[j] * f for f, c in zip(function_values, stencil)))
def frobenius_norm(matrix, factor=1):
"""Computes the Frobenius norm of a matrix defined as the square root of the sum of squared matrix elements
The optional factor is added inside the square root"""
return sp.sqrt(sum(i * i for i in matrix) * factor)
def add_smagorinsky_model(collision_rule, smagorinsky_constant, omega_output_field=None):
r""" Adds a smagorinsky model to a lattice Boltzmann collision rule. To add the Smagorinsky model to a LB scheme
one has to first compute the strain rate tensor $S_{ij}$ in each cell, and compute the turbulent
viscosity :math:`nu_t` from it. Then the local relaxation rate has to be adapted to match the total viscosity
:math `\nu_{total}` instead of the standard viscosity :math `\nu_0`.
A fortunate property of LB methods is, that the strain rate tensor can be computed locally from the
non-equilibrium part of the distribution function. This is somewhat surprising, since the strain rate tensor
contains first order derivatives. The strain rate tensor can be obtained by
.. math ::
S_{ij} = - \frac{3 \omega_s}{2 \rho_{(0)}} \Pi_{ij}^{(neq)}
where :math `\omega_s` is the relaxation rate that determines the viscosity, :math `\rho_{(0)}` is :math `\rho`
in compressible models and :math `1` for incompressible schemes.
:math `\Pi_{ij}^{(neq)}` is the second order moment tensor of the non-equilibrium part of
the distribution functions
:math `f^{(neq)} = f - f^{(eq)}` and can be computed as
.. math ::
\Pi_{ij}^{(neq)} = \sum_q c_{qi} c_{qj} \; f_q^{(neq)}
"""
method = collision_rule.method
omega_s = get_shear_relaxation_rate(method)
if isinstance(omega_s, float) or isinstance(omega_s, int):
raise ValueError("For the smagorinsky model the shear relaxation rate has to be a symbol")
f_neq = sp.Matrix(method.pre_collision_pdf_symbols) - method.get_equilibrium_terms()
tau_0 = sp.Symbol("tau_0_")
second_order_neq_moments = sp.Symbol("Pi")
rho = method.zeroth_order_equilibrium_moment_symbol if method.conserved_quantity_computation.compressible else 1
adapted_omega = sp.Symbol("smagorinsky_omega")
collision_rule = collision_rule.new_with_substitutions({omega_s: adapted_omega})
# for derivation see notebook demo_custom_LES_model.ipynb
eqs = [Assignment(tau_0, 1 / omega_s),
Assignment(second_order_neq_moments,
frobenius_norm(second_order_moment_tensor(f_neq, method.stencil), factor=2) / rho),
Assignment(adapted_omega,
2 / (tau_0 + sp.sqrt(18 * smagorinsky_constant ** 2 * second_order_neq_moments + tau_0 ** 2)))]
collision_rule.subexpressions += eqs
collision_rule.topological_sort(sort_subexpressions=True, sort_main_assignments=False)
if omega_output_field:
collision_rule.main_assignments.append(Assignment(omega_output_field.center, adapted_omega))
return collision_rule
lbmpy_tests/advanced_setups/Lees_Edwards_boundary_conditions.ipynb deleted 100644 → 0
View file @ cf17a9b9
Source diff could not be displayed: it is too large. Options to address this: view the blob.
lbmpy_tests/advanced_streaming/test_communication.py deleted 100644 → 0
View file @ cf17a9b9
import pystencils as ps
import numpy as np
from lbmpy.stencils import get_stencil
from pystencils.slicing import get_slice_before_ghost_layer, get_ghost_region_slice
from lbmpy.creationfunctions import create_lb_update_rule
from lbmpy.advanced_streaming.communication import get_communication_slices, _fix_length_one_slices, \
LBMPeriodicityHandling
from lbmpy.advanced_streaming.utility import streaming_patterns, Timestep
import pytest
@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q15', 'D3Q19', 'D3Q27'])
@pytest.mark.parametrize('streaming_pattern', streaming_patterns)
@pytest.mark.parametrize('timestep', [Timestep.EVEN, Timestep.ODD])
def test_slices_not_empty(stencil, streaming_pattern, timestep):
stencil = get_stencil(stencil)
dim = len(stencil[0])
q = len(stencil)
arr = np.zeros((4,) * dim + (q,))
slices = get_communication_slices(stencil, streaming_pattern=streaming_pattern, prev_timestep=timestep,
ghost_layers=1)
for _, slices_list in slices.items():
for src, dst in slices_list:
assert all(s != 0 for s in arr[src].shape)
assert all(s != 0 for s in arr[dst].shape)
@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q15', 'D3Q19', 'D3Q27'])
@pytest.mark.parametrize('streaming_pattern', streaming_patterns)
@pytest.mark.parametrize('timestep', [Timestep.EVEN, Timestep.ODD])
def test_src_dst_same_shape(stencil, streaming_pattern, timestep):
stencil = get_stencil(stencil)
dim = len(stencil[0])
q = len(stencil)
arr = np.zeros((4,) * dim + (q,))
slices = get_communication_slices(stencil, streaming_pattern=streaming_pattern, prev_timestep=timestep,
ghost_layers=1)
for _, slices_list in slices.items():
for src, dst in slices_list:
src_shape = arr[src].shape
dst_shape = arr[dst].shape
assert src_shape == dst_shape
@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q15', 'D3Q19', 'D3Q27'])
def test_pull_communication_slices(stencil):
stencil = get_stencil(stencil)
slices = get_communication_slices(
stencil, streaming_pattern='pull', prev_timestep=Timestep.BOTH, ghost_layers=1)
for i, d in enumerate(stencil):
if i == 0:
continue
for s in slices[d]:
if s[0][-1] == i:
src = s[0][:-1]
dst = s[1][:-1]
break
inner_slice = _fix_length_one_slices(get_slice_before_ghost_layer(d, ghost_layers=1))
inv_dir = (-e for e in d)
gl_slice = _fix_length_one_slices(get_ghost_region_slice(inv_dir, ghost_layers=1))
assert src == inner_slice
assert dst == gl_slice
@pytest.mark.parametrize('stencil_name', ['D2Q9', 'D3Q15', 'D3Q19', 'D3Q27'])
def test_optimised_and_full_communication_equivalence(stencil_name):
target = 'cpu'
stencil = get_stencil(stencil_name)
dim = len(stencil[0])
domain_size = (4, ) * dim
dh = ps.create_data_handling(domain_size, periodicity=(True, ) * dim,
parallel=False, default_target=target)
pdf = dh.add_array("pdf", values_per_cell=len(stencil), dtype=np.int64)
dh.fill("pdf", 0, ghost_layers=True)
pdf_tmp = dh.add_array("pdf_tmp", values_per_cell=len(stencil), dtype=np.int64)
dh.fill("pdf_tmp", 0, ghost_layers=True)
gl = dh.ghost_layers_of_field("pdf")
num = 0
for idx, x in np.ndenumerate(dh.cpu_arrays['pdf']):
dh.cpu_arrays['pdf'][idx] = num
dh.cpu_arrays['pdf_tmp'][idx] = num
num += 1
ac = create_lb_update_rule(stencil=stencil,
optimization={"symbolic_field": pdf,
"symbolic_temporary_field": pdf_tmp},
kernel_type='stream_pull_only')
ast = ps.create_kernel(ac, target=dh.default_target, cpu_openmp=True)
stream = ast.compile()
full_communication = dh.synchronization_function(pdf.name, target=dh.default_target, optimization={"openmp": True})
full_communication()
dh.run_kernel(stream)
dh.swap("pdf", "pdf_tmp")
pdf_full_communication = np.copy(dh.cpu_arrays['pdf'])
num = 0
for idx, x in np.ndenumerate(dh.cpu_arrays['pdf']):
dh.cpu_arrays['pdf'][idx] = num
dh.cpu_arrays['pdf_tmp'][idx] = num
num += 1
optimised_communication = LBMPeriodicityHandling(stencil=stencil, data_handling=dh, pdf_field_name=pdf.name,
streaming_pattern='pull')
optimised_communication()
dh.run_kernel(stream)
dh.swap("pdf", "pdf_tmp")
if dim == 3:
for i in range(gl, domain_size[0]):
for j in range(gl, domain_size[1]):
for k in range(gl, domain_size[2]):
for f in range(len(stencil)):
assert dh.cpu_arrays['pdf'][i, j, k, f] == pdf_full_communication[i, j, k, f], print(f)
else:
for i in range(gl, domain_size[0]):
for j in range(gl, domain_size[1]):
for f in range(len(stencil)):
assert dh.cpu_arrays['pdf'][i, j, f] == pdf_full_communication[i, j, f]
lbmpy_tests/centeredcumulant/test_central_moment_transform.py deleted 100644 → 0
View file @ cf17a9b9
import sympy as sp
from lbmpy.moments import get_default_moment_set_for_stencil, extract_monomials
import pytest
from lbmpy.stencils import get_stencil
from lbmpy.methods.centeredcumulant.centered_cumulants import statistical_quantity_symbol
from lbmpy.methods.momentbased.moment_transforms import (
PdfsToCentralMomentsByMatrix, FastCentralMomentTransform, PdfsToCentralMomentsByShiftMatrix,
PRE_COLLISION_CENTRAL_MOMENT, POST_COLLISION_CENTRAL_MOMENT
)
@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q19', 'D3Q27'])
def test_forward_transform(stencil):
stencil = get_stencil(stencil)
dim = len(stencil[0])
q = len(stencil)
moment_exponents = list(extract_monomials(get_default_moment_set_for_stencil(stencil)))
moment_exponents = sorted(moment_exponents, key=sum)
pdfs = sp.symbols(f"f_:{q}")
rho = sp.Symbol('rho')
u = sp.symbols(f"u_:{dim}")
matrix_transform = PdfsToCentralMomentsByMatrix(stencil, moment_exponents, rho, u)
fast_transform = FastCentralMomentTransform(stencil, moment_exponents, rho, u)
shift_transform = PdfsToCentralMomentsByShiftMatrix(stencil, moment_exponents, rho, u)
f_to_k_matrix = matrix_transform.forward_transform(pdfs, moment_symbol_base=PRE_COLLISION_CENTRAL_MOMENT)
f_to_k_matrix = f_to_k_matrix.new_without_subexpressions().main_assignments_dict
f_to_k_fast = fast_transform.forward_transform(pdfs, moment_symbol_base=PRE_COLLISION_CENTRAL_MOMENT)
f_to_k_fast = f_to_k_fast.new_without_subexpressions().main_assignments_dict
f_to_k_shift = shift_transform.forward_transform(pdfs, moment_symbol_base=PRE_COLLISION_CENTRAL_MOMENT,
simplification=False)
f_to_k_shift = f_to_k_shift.new_without_subexpressions().main_assignments_dict
for e in moment_exponents:
moment_symbol = statistical_quantity_symbol(PRE_COLLISION_CENTRAL_MOMENT, e)
rhs_matrix = f_to_k_matrix[moment_symbol].expand()
rhs_fast = f_to_k_fast[moment_symbol].expand()
rhs_shift = f_to_k_shift[moment_symbol].expand()
assert (rhs_matrix - rhs_fast) == 0, f"Mismatch between matrix and fast transform at {moment_symbol}."
assert (rhs_matrix - rhs_shift) == 0, f"Mismatch between matrix and shift-matrix transform at {moment_symbol}."
@pytest.mark.parametrize('stencil', ['D2Q9', 'D3Q19', 'D3Q27'])
def test_backward_transform(stencil):
stencil = get_stencil(stencil)
dim = len(stencil[0])
q = len(stencil)
moment_exponents = list(extract_monomials(get_default_moment_set_for_stencil(stencil)))
moment_exponents = sorted(moment_exponents, key=sum)
pdfs = sp.symbols(f"f_:{q}")
rho = sp.Symbol('rho')
u = sp.symbols(f"u_:{dim}")
matrix_transform = PdfsToCentralMomentsByMatrix(stencil, moment_exponents, rho, u)
fast_transform = FastCentralMomentTransform(stencil, moment_exponents, rho, u)
shift_transform = PdfsToCentralMomentsByShiftMatrix(stencil, moment_exponents, rho, u)
k_to_f_matrix = matrix_transform.backward_transform(pdfs, moment_symbol_base=POST_COLLISION_CENTRAL_MOMENT)
k_to_f_matrix = k_to_f_matrix.new_without_subexpressions().main_assignments_dict
k_to_f_fast = fast_transform.backward_transform(pdfs, moment_symbol_base=POST_COLLISION_CENTRAL_MOMENT)
k_to_f_fast = k_to_f_fast.new_without_subexpressions().main_assignments_dict
k_to_f_shift = shift_transform.backward_transform(pdfs, moment_symbol_base=POST_COLLISION_CENTRAL_MOMENT)
k_to_f_shift = k_to_f_shift.new_without_subexpressions().main_assignments_dict
for f in pdfs:
rhs_matrix = k_to_f_matrix[f].expand()
rhs_fast = k_to_f_fast[f].expand()
rhs_shift = k_to_f_shift[f].expand()
assert (rhs_matrix - rhs_fast) == 0, f"Mismatch between matrix and fast transform at {f}."
assert (rhs_matrix - rhs_shift) == 0, f"Mismatch between matrix and shift-matrix transform at {f}."
lbmpy_tests/centeredcumulant/test_periodic_pipe_flow.ipynb deleted 100644 → 0
View file @ cf17a9b9
%% Cell type:code id: tags:
``` python
import pytest
pytest.importorskip('pycuda')
```
%% Output
<module 'pycuda' from '/home/markus/miniconda3/envs/pystencils/lib/python3.8/site-packages/pycuda/__init__.py'>
%% Cell type:code id: tags:
``` python
%load_ext autoreload
%autoreload 2
```
%% Cell type:code id: tags:
``` python
import sympy as sp
import numpy as np
from numpy.testing import assert_allclose
from pystencils.session import *
from lbmpy.session import *
from lbmpy.stencils import get_stencil
from lbmpy.methods.centeredcumulant import CenteredCumulantForceModel
from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
```
%% Cell type:markdown id: tags:
## Pipe Flow Scenario
%% Cell type:code id: tags:
``` python
class PeriodicPipeFlow:
def __init__(self, method_params, optimization=dict()):
wall_boundary = NoSlip()
self.stencil = method_params.get('stencil', get_stencil("D2Q9"))
self.streaming_pattern = method_params.get('streaming_pattern', 'pull')
self.target = optimization.get('target', 'cpu')
self.gpu = self.target in ['gpu', 'opencl']
# Stencil
self.q = len(self.stencil)
self.dim = len(self.stencil[0])
# Streaming
self.inplace = is_inplace(self.streaming_pattern)
self.timesteps = get_timesteps(self.streaming_pattern)
self.zeroth_timestep = self.timesteps[0]
# Domain, Data Handling and PDF fields
self.pipe_length = 60
self.pipe_radius = 15
self.domain_size = (self.pipe_length, ) + (2 * self.pipe_radius,) * (self.dim - 1)
self.periodicity = (True, ) + (False, ) * (self.dim - 1)
force = (0.0001, ) + (0.0,) * (self.dim - 1)
self.force = method_params.get('force', force)
self.dh = create_data_handling(domain_size=self.domain_size,
periodicity=self.periodicity, default_target=self.target)
self.pdfs = self.dh.add_array('pdfs', self.q)
if not self.inplace:
self.pdfs_tmp = self.dh.add_array_like('pdfs_tmp', self.pdfs.name)
method_params['force'] = self.force
optimization['symbolic_field'] = self.pdfs
if not self.inplace:
optimization['symbolic_temporary_field'] = self.pdfs_tmp
self.lb_collision = create_lb_collision_rule(optimization=optimization, **method_params)
self.lb_method = self.lb_collision.method
self.lb_kernels = []
for t in self.timesteps:
self.lb_kernels.append(create_lb_function(collision_rule=self.lb_collision,
optimization=optimization,
timestep=t,
**method_params))
# Macroscopic Values
self.density = 1.0
self.density_field = self.dh.add_array('rho', 1)
u_x = 0.0
self.velocity = (u_x,) * self.dim
self.velocity_field = self.dh.add_array('u', self.dim)
setter = macroscopic_values_setter(
self.lb_method, self.density, self.velocity, self.pdfs,
streaming_pattern=self.streaming_pattern, previous_timestep=self.zeroth_timestep)
self.init_kernel = create_kernel(setter, ghost_layers=1, target=self.target).compile()
self.getter_kernels = []
for t in self.timesteps:
getter = macroscopic_values_getter(
self.lb_method, self.density_field, self.velocity_field, self.pdfs,
streaming_pattern=self.streaming_pattern, previous_timestep=t)
self.getter_kernels.append(create_kernel(getter, ghost_layers=1, target=self.target).compile())
# Periodicity
self.periodicity_handler = LBMPeriodicityHandling(
self.stencil, self.dh, self.pdfs.name, streaming_pattern=self.streaming_pattern)
# Boundary Handling
self.wall = wall_boundary
self.bh = LatticeBoltzmannBoundaryHandling(
self.lb_method, self.dh, self.pdfs.name,
streaming_pattern=self.streaming_pattern, target=self.target)
self.bh.set_boundary(boundary_obj=self.wall, mask_callback=self.mask_callback)
self.current_timestep = self.zeroth_timestep
def mask_callback(self, x, y, z=None):
y = y - self.pipe_radius
z = z - self.pipe_radius if z is not None else 0
return np.sqrt(y**2 + z**2) >= self.pipe_radius
def init(self):
self.current_timestep = self.zeroth_timestep
self.dh.run_kernel(self.init_kernel)
def step(self):
# Order matters! First communicate, then boundaries, otherwise
# periodicity handling overwrites reflected populations
# Periodicty
self.periodicity_handler(self.current_timestep)
# Boundaries
self.bh(prev_timestep=self.current_timestep)
# Here, the next time step begins
self.current_timestep = self.current_timestep.next()
# LBM Step
self.dh.run_kernel(self.lb_kernels[self.current_timestep.idx])
# Field Swaps
if not self.inplace:
self.dh.swap(self.pdfs.name, self.pdfs_tmp.name)
# Macroscopic Values
self.dh.run_kernel(self.getter_kernels[self.current_timestep.idx])
def run(self, iterations):
for _ in range(iterations):
self.step()
@property
def velocity_array(self):
if self.gpu:
self.dh.to_cpu(self.velocity_field.name)
return self.dh.gather_array(self.velocity_field.name)
def get_trimmed_velocity_array(self):
if self.gpu:
self.dh.to_cpu(self.velocity_field.name)
u = np.copy(self.dh.gather_array(self.velocity_field.name))
mask = self.bh.get_mask(None, self.wall)
for idx in np.ndindex(u.shape[:-1]):
if mask[idx] != 0:
u[idx] = np.full((self.dim, ), np.nan)
return u
```
%% Cell type:markdown id: tags:
## General Setup
%% Cell type:code id: tags:
``` python
stencil = get_stencil('D3Q19')
dim = len(stencil[0])
target = 'gpu'
streaming_pattern = 'aa'
optimization = { 'target': target }
viscous_rr = 1.1
force = (0.0001, ) + (0.0,) * (dim - 1)
```
%% Cell type:markdown id: tags:
## 1. Reference: SRT Method
%% Cell type:code id: tags:
``` python
srt_params = {
'stencil': stencil,
'method': 'srt',
'relaxation_rate': viscous_rr,
'force_model': 'guo',
'force' : force,
'streaming_pattern': streaming_pattern
}
srt_flow = PeriodicPipeFlow(srt_params, optimization)
srt_flow.init()
srt_flow.run(400)
```
%% Cell type:code id: tags:
``` python
srt_u = srt_flow.get_trimmed_velocity_array()
ps.plot.vector_field_magnitude(srt_u[30,:,:,:])
ps.plot.colorbar()
```
%% Output
<matplotlib.colorbar.Colorbar at 0x7f20b2f75bb0>
%% Cell type:markdown id: tags:
## 2. Centered Cumulant Method with Implicit Forcing
The default setup of the centered cumulant method does not use a force model to compute the force contributions on all populations, but instead applies the force in cumulant space simply by setting the momentum relaxation rate $\omega_1 = 2$. Due to the half-force shift of the equilibrium input velocity, the first central moments (which correspond to momentum relative to the moving frame of reference) are not zero, but equal to $-F/2$. The relaxation process causes the first central moments to simply change sign:
$$
\kappa_{100}^* = - \kappa_{100}.
$$
Thus, $\kappa_{100}^* = F_x /2$. In total, the entire acceleration given by the force is added onto the momentum.
%% Cell type:code id: tags:
``` python
cm_method_params = {
'method' : 'monomial_cumulant',
'stencil': stencil,
'relaxation_rate': viscous_rr, # Specify viscous relaxation rate only
'force_model': CenteredCumulantForceModel(force),
'force' : force,
'streaming_pattern' : streaming_pattern
}
optimization = {
'target': target,
'pre_simplification' : True
}
cm_impl_f_flow = PeriodicPipeFlow(cm_method_params, optimization)
```
%% Cell type:code id: tags:
``` python
lb_method = cm_impl_f_flow.lb_method
assert all(rr == 2 for rr in lb_method.relaxation_rates[1:4])
assert all(rr == viscous_rr for rr in lb_method.relaxation_rates[4:9])
```
%% Cell type:code id: tags:
``` python
cm_impl_f_flow.init()
cm_impl_f_flow.run(400)
```
%% Cell type:code id: tags:
``` python
cm_impl_f_u = cm_impl_f_flow.get_trimmed_velocity_array()
ps.plot.vector_field_magnitude(cm_impl_f_u[30,:,:,:])
ps.plot.colorbar()
```
%% Output
<matplotlib.colorbar.Colorbar at 0x7f20aa846670>
%% Cell type:code id: tags:
``` python
assert_allclose(cm_impl_f_u, srt_u, rtol=1, atol=1e-4)
```
%% Cell type:markdown id: tags:
## 3. Centered Cumulant Method with Explicit Forcing
%% Cell type:code id: tags:
``` python
from lbmpy.forcemodels import Schiller
cm_method_params_expl_force = {
'method' : 'cumulant',
'stencil': stencil,
'relaxation_rates': [viscous_rr], # Specify viscous relaxation rate only
'force_model': Schiller(force),
'force' : force,
'streaming_pattern' : streaming_pattern
}
optimization = {
'target': target,
'pre_simplification' : True
}
cm_expl_f_flow = PeriodicPipeFlow(cm_method_params_expl_force, optimization)
```
%% Cell type:code id: tags:
``` python
lb_method = cm_expl_f_flow.lb_method
assert all(rr == 0 for rr in lb_method.relaxation_rates[1:4])
assert all(rr == viscous_rr for rr in lb_method.relaxation_rates[4:9])
```
%% Cell type:code id: tags:
``` python
cm_expl_f_flow.init()
cm_expl_f_flow.run(400)
```
%% Cell type:code id: tags:
``` python
cm_expl_f_u = cm_expl_f_flow.get_trimmed_velocity_array()
ps.plot.vector_field_magnitude(cm_expl_f_u[30,:,:,:])
ps.plot.colorbar()
```
%% Output
<matplotlib.colorbar.Colorbar at 0x7f20aa641070>
%% Cell type:code id: tags:
``` python
assert_allclose(cm_expl_f_u, srt_u, rtol=1, atol=1e-4)
assert_allclose(cm_expl_f_u, cm_impl_f_u, rtol=1, atol=1e-4)
```
lbmpy_tests/phasefield/test_analytical_3phase_model.ipynb deleted 100644 → 0
View file @ cf17a9b9
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.phasefield.analytical import *
from pystencils.fd import evaluate_diffs
```
%% Cell type:markdown id: tags:
# Analytical checks for 3-Phase model
Here you can inspect the components of the free energy. View only bulk or interface contributions, before and after transformation from $U \rightarrow (\rho, \phi,\psi)$:
%% Cell type:code id: tags:
``` python
order_parameters = sp.symbols("rho phi psi")
rho, phi, psi = order_parameters
F, _ = free_energy_functional_3_phases(include_bulk=True,
include_interface=True,
transformed=True,
expand_derivatives=False)
F
```
%% Output
$$\frac{\alpha^{2} \kappa_{0}}{2} {\partial (\frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}) }^{2} + \frac{\alpha^{2} \kappa_{1}}{2} {\partial (- \frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}) }^{2} + \frac{\alpha^{2} \kappa_{2}}{2} {\partial \psi}^{2} + \frac{\kappa_{0}}{2} \left(\frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}\right)^{2} \left(- \frac{\phi}{2} + \frac{\psi}{2} - \frac{\rho}{2} + 1\right)^{2} + \frac{\kappa_{1}}{2} \left(- \frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}\right)^{2} \left(\frac{\phi}{2} + \frac{\psi}{2} - \frac{\rho}{2} + 1\right)^{2} + \frac{\kappa_{2} \psi^{2}}{2} \left(- \psi + 1\right)^{2}$$
2 2 2 2 2
α ⋅κ₀⋅D(phi/2 - psi/2 + rho/2) α ⋅κ₁⋅D(-phi/2 - psi/2 + rho/2) α ⋅κ₂⋅D(p
─────────────────────────────── + ──────────────────────────────── + ─────────
2 2 2
2 2 2 2
⎛φ ψ ρ⎞ ⎛ φ ψ ρ ⎞ ⎛ φ ψ ρ⎞ ⎛φ ψ ρ ⎞
2 κ₀⋅⎜─ - ─ + ─⎟ ⋅⎜- ─ + ─ - ─ + 1⎟ κ₁⋅⎜- ─ - ─ + ─⎟ ⋅⎜─ + ─ - ─ + 1⎟
si) ⎝2 2 2⎠ ⎝ 2 2 2 ⎠ ⎝ 2 2 2⎠ ⎝2 2 2 ⎠
──── + ────────────────────────────────── + ──────────────────────────────────
2 2
2 2
κ₂⋅ψ ⋅(-ψ + 1)
+ ───────────────
2
%% Cell type:markdown id: tags:
### Analytically checking the phase transition profile
%% Cell type:markdown id: tags:
Automatically deriving chemical potential as functional derivative of free energy
%% Cell type:code id: tags:
``` python
F, _ = free_energy_functional_3_phases(order_parameters)
mu_diff_eq = chemical_potentials_from_free_energy(F, order_parameters)
mu_diff_eq[0]
```
%% Output
$$- \frac{\alpha^{2} \kappa_{0}}{4} {\partial {\partial \phi}} + \frac{\alpha^{2} \kappa_{0}}{4} {\partial {\partial \psi}} - \frac{\alpha^{2} \kappa_{0}}{4} {\partial {\partial \rho}} + \frac{\alpha^{2} \kappa_{1}}{4} {\partial {\partial \phi}} + \frac{\alpha^{2} \kappa_{1}}{4} {\partial {\partial \psi}} - \frac{\alpha^{2} \kappa_{1}}{4} {\partial {\partial \rho}} + \frac{\kappa_{0} \phi^{3}}{8} - \frac{3 \kappa_{0}}{8} \phi^{2} \psi + \frac{3 \kappa_{0}}{8} \phi^{2} \rho - \frac{3 \kappa_{0}}{8} \phi^{2} + \frac{3 \kappa_{0}}{8} \phi \psi^{2} - \frac{3 \kappa_{0}}{4} \phi \psi \rho + \frac{3 \kappa_{0}}{4} \phi \psi + \frac{3 \kappa_{0}}{8} \phi \rho^{2} - \frac{3 \kappa_{0}}{4} \phi \rho + \frac{\kappa_{0} \phi}{4} - \frac{\kappa_{0} \psi^{3}}{8} + \frac{3 \kappa_{0}}{8} \psi^{2} \rho - \frac{3 \kappa_{0}}{8} \psi^{2} - \frac{3 \kappa_{0}}{8} \psi \rho^{2} + \frac{3 \kappa_{0}}{4} \psi \rho - \frac{\kappa_{0} \psi}{4} + \frac{\kappa_{0} \rho^{3}}{8} - \frac{3 \kappa_{0}}{8} \rho^{2} + \frac{\kappa_{0} \rho}{4} - \frac{\kappa_{1} \phi^{3}}{8} - \frac{3 \kappa_{1}}{8} \phi^{2} \psi + \frac{3 \kappa_{1}}{8} \phi^{2} \rho - \frac{3 \kappa_{1}}{8} \phi^{2} - \frac{3 \kappa_{1}}{8} \phi \psi^{2} + \frac{3 \kappa_{1}}{4} \phi \psi \rho - \frac{3 \kappa_{1}}{4} \phi \psi - \frac{3 \kappa_{1}}{8} \phi \rho^{2} + \frac{3 \kappa_{1}}{4} \phi \rho - \frac{\kappa_{1} \phi}{4} - \frac{\kappa_{1} \psi^{3}}{8} + \frac{3 \kappa_{1}}{8} \psi^{2} \rho - \frac{3 \kappa_{1}}{8} \psi^{2} - \frac{3 \kappa_{1}}{8} \psi \rho^{2} + \frac{3 \kappa_{1}}{4} \psi \rho - \frac{\kappa_{1} \psi}{4} + \frac{\kappa_{1} \rho^{3}}{8} - \frac{3 \kappa_{1}}{8} \rho^{2} + \frac{\kappa_{1} \rho}{4}$$
2 2 2 2 2
α ⋅κ₀⋅D(D(phi)) α ⋅κ₀⋅D(D(psi)) α ⋅κ₀⋅D(D(rho)) α ⋅κ₁⋅D(D(phi)) α ⋅κ
- ─────────────── + ─────────────── - ─────────────── + ─────────────── + ────
4 4 4 4
2 3 2 2 2
₁⋅D(D(psi)) α ⋅κ₁⋅D(D(rho)) κ₀⋅φ 3⋅κ₀⋅φ ⋅ψ 3⋅κ₀⋅φ ⋅ρ 3⋅κ₀⋅φ 3⋅κ₀
─────────── - ─────────────── + ───── - ───────── + ───────── - ─────── + ────
4 4 8 8 8 8
2 2 3 2
⋅φ⋅ψ 3⋅κ₀⋅φ⋅ψ⋅ρ 3⋅κ₀⋅φ⋅ψ 3⋅κ₀⋅φ⋅ρ 3⋅κ₀⋅φ⋅ρ κ₀⋅φ κ₀⋅ψ 3⋅κ₀⋅ψ ⋅
───── - ────────── + ──────── + ───────── - ──────── + ──── - ───── + ────────
8 4 4 8 4 4 8 8
2 2 3 2 3
ρ 3⋅κ₀⋅ψ 3⋅κ₀⋅ψ⋅ρ 3⋅κ₀⋅ψ⋅ρ κ₀⋅ψ κ₀⋅ρ 3⋅κ₀⋅ρ κ₀⋅ρ κ₁⋅φ 3
─ - ─────── - ───────── + ──────── - ──── + ───── - ─────── + ──── - ───── - ─
8 8 4 4 8 8 4 8
2 2 2 2 2
⋅κ₁⋅φ ⋅ψ 3⋅κ₁⋅φ ⋅ρ 3⋅κ₁⋅φ 3⋅κ₁⋅φ⋅ψ 3⋅κ₁⋅φ⋅ψ⋅ρ 3⋅κ₁⋅φ⋅ψ 3⋅κ₁⋅φ⋅ρ
──────── + ───────── - ─────── - ───────── + ────────── - ──────── - ─────────
8 8 8 8 4 4 8
3 2 2 2
3⋅κ₁⋅φ⋅ρ κ₁⋅φ κ₁⋅ψ 3⋅κ₁⋅ψ ⋅ρ 3⋅κ₁⋅ψ 3⋅κ₁⋅ψ⋅ρ 3⋅κ₁⋅ψ⋅ρ κ₁⋅ψ
+ ──────── - ──── - ───── + ───────── - ─────── - ───────── + ──────── - ────
4 4 8 8 8 8 4 4
3 2
κ₁⋅ρ 3⋅κ₁⋅ρ κ₁⋅ρ
+ ───── - ─────── + ────
8 8 4
%% Cell type:markdown id: tags:
Checking if expected profile is a solution of the differential equation:
%% Cell type:code id: tags:
``` python
x = sp.symbols("x")
expectedProfile = analytic_interface_profile(x)
expectedProfile
```
%% Output
$$\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}$$
⎛ x ⎞
tanh⎜───⎟
⎝2⋅α⎠ 1
───────── + ─
2 2
%% Cell type:markdown id: tags:
Checking a phase transition from $C_0$ to $C_2$. This means that $\rho=1$ while $phi$ and $psi$ are the analytical profile or 1-analytical profile
%% Cell type:code id: tags:
``` python
for eq in mu_diff_eq:
eq = eq.subs({rho: 1,
phi: 1 - expectedProfile,
psi: expectedProfile})
eq = evaluate_diffs(eq, x).expand()
assert eq == 0
```
%% Cell type:markdown id: tags:
### Checking the surface tensions parameters
%% Cell type:code id: tags:
``` python
F, _ = free_energy_functional_3_phases(order_parameters)
F = expand_diff_linear(F, functions=order_parameters) # expand derivatives using product rule
two_phase_free_energy = F.subs({rho: 1,
phi: 1 - expectedProfile,
psi: expectedProfile})
two_phase_free_energy = sp.simplify(evaluate_diffs(two_phase_free_energy, x))
two_phase_free_energy
```
%% Output
$$\frac{\kappa_{0} + \kappa_{2}}{16 \cosh^{4}{\left (\frac{x}{2 \alpha} \right )}}$$
κ₀ + κ₂
─────────────
4⎛ x ⎞
16⋅cosh ⎜───⎟
⎝2⋅α⎠
%% Cell type:code id: tags:
``` python
gamma = cosh_integral(two_phase_free_energy, x)
gamma
```
%% Output
$$\frac{\alpha}{6} \left(\kappa_{0} + \kappa_{2}\right)$$
α⋅(κ₀ + κ₂)
───────────
6
%% Cell type:code id: tags:
``` python
alpha, k0, k2 = sp.symbols("alpha, kappa_0, kappa_2")
assert gamma == alpha/6 * (k0 + k2)
```
lbmpy_tests/phasefield/test_analytical_n_phase.ipynb deleted 100644 → 0
View file @ cf17a9b9
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.phasefield.analytical import *
from functools import partial
```
%% Cell type:markdown id: tags:
# Analytical checks for N-Phase model
%% Cell type:markdown id: tags:
### Formulation of free energy
%% Cell type:markdown id: tags:
In the next cell you can inspect the formulation of the free energy functional.
Bulk and interface terms can be viewed separately.
%% Cell type:code id: tags:
``` python
num_phases = 3
order_params = symbolic_order_parameters(num_symbols=num_phases-1)
f2 = partial(n_phases_correction_function, beta=1)
F = free_energy_functional_n_phases(order_parameters=order_params,
include_interface=False, f2=f2,
include_bulk=True )
F
```
%% Output
$$\frac{3}{2 \alpha} \left(\tau_{0 1} \left(\phi_{0}^{2} \left(- \phi_{0} + 1\right)^{2} + \phi_{1}^{2} \left(- \phi_{1} + 1\right)^{2} - \begin{cases} - \left(\phi_{0} + \phi_{1}\right)^{2} & \text{for}\: \phi_{0} + \phi_{1} < 0 \\- \left(- \phi_{0} - \phi_{1} + 1\right)^{2} & \text{for}\: \phi_{0} + \phi_{1} > 1 \\\left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} & \text{otherwise} \end{cases}\right) + \tau_{0 2} \left(\phi_{0}^{2} \left(- \phi_{0} + 1\right)^{2} + \left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} - \begin{cases} - \left(- \phi_{1} + 1\right)^{2} & \text{for}\: - \phi_{1} + 1 < 0 \\- \phi_{1}^{2} & \text{for}\: - \phi_{1} + 1 > 1 \\\phi_{1}^{2} \left(- \phi_{1} + 1\right)^{2} & \text{otherwise} \end{cases}\right) + \tau_{1 2} \left(\phi_{1}^{2} \left(- \phi_{1} + 1\right)^{2} + \left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} - \begin{cases} - \left(- \phi_{0} + 1\right)^{2} & \text{for}\: - \phi_{0} + 1 < 0 \\- \phi_{0}^{2} & \text{for}\: - \phi_{0} + 1 > 1 \\\phi_{0}^{2} \left(- \phi_{0} + 1\right)^{2} & \text{otherwise} \end{cases}\right)\right)$$
⎛ ⎛ ⎛⎧ 2
⎜ ⎜ ⎜⎪ -(φ₀ + φ₁) for φ
⎜ ⎜ ⎜⎪
⎜ ⎜ 2 2 2 2 ⎜⎪ 2
3⋅⎜τ₀ ₁⋅⎜φ₀ ⋅(-φ₀ + 1) + φ₁ ⋅(-φ₁ + 1) - ⎜⎨ -(-φ₀ - φ₁ + 1) for φ
⎜ ⎜ ⎜⎪
⎜ ⎜ ⎜⎪ 2 2
⎜ ⎜ ⎜⎪(φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1) ot
⎝ ⎝ ⎝⎩
──────────────────────────────────────────────────────────────────────────────
⎞⎞ ⎛ ⎛⎧
₀ + φ₁ < 0⎟⎟ ⎜ ⎜⎪ -(-φ₁ +
⎟⎟ ⎜ ⎜⎪
⎟⎟ ⎜ 2 2 2 2 ⎜⎪ 2
₀ + φ₁ > 1⎟⎟ + τ₀ ₂⋅⎜φ₀ ⋅(-φ₀ + 1) + (φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1) - ⎜⎨ -φ₁
⎟⎟ ⎜ ⎜⎪
⎟⎟ ⎜ ⎜⎪ 2
herwise ⎟⎟ ⎜ ⎜⎪φ₁ ⋅(-φ₁
⎠⎠ ⎝ ⎝⎩
──────────────────────────────────────────────────────────────────────────────
2⋅α
2 ⎞⎞ ⎛
1) for -φ₁ + 1 < 0⎟⎟ ⎜
⎟⎟ ⎜
⎟⎟ ⎜ 2 2 2 2
for -φ₁ + 1 > 1⎟⎟ + τ₁ ₂⋅⎜φ₁ ⋅(-φ₁ + 1) + (φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1) -
⎟⎟ ⎜
2 ⎟⎟ ⎜
+ 1) otherwise ⎟⎟ ⎜
⎠⎠ ⎝
──────────────────────────────────────────────────────────────────────────────
⎛⎧ 2 ⎞⎞⎞
⎜⎪ -(-φ₀ + 1) for -φ₀ + 1 < 0⎟⎟⎟
⎜⎪ ⎟⎟⎟
⎜⎪ 2 ⎟⎟⎟
⎜⎨ -φ₀ for -φ₀ + 1 > 1⎟⎟⎟
⎜⎪ ⎟⎟⎟
⎜⎪ 2 2 ⎟⎟⎟
⎜⎪φ₀ ⋅(-φ₀ + 1) otherwise ⎟⎟⎟
⎝⎩ ⎠⎠⎠
─────────────────────────────────────
%% Cell type:markdown id: tags:
### Analytically checking the phase transition profile
First we define the order parameters and free energy, including bulk and interface terms:
%% Cell type:code id: tags:
``` python
F = free_energy_functional_n_phases(order_parameters=symbolic_order_parameters(num_symbols=num_phases-1))
```
%% Cell type:markdown id: tags:
Then we automatically derive the differential equations for the chemial potential $\mu$
%% Cell type:code id: tags:
``` python
mu_diff_eq = chemical_potentials_from_free_energy(F, order_params)
# there is one equation less than phases
assert len(mu_diff_eq) == num_phases - 1
# show the first one
mu_diff_eq[0]
```
%% Output
$$3 \alpha \tau_{0 1} {\partial {\partial \phi_{1}}} - 6 \alpha \tau_{0 2} {\partial {\partial \phi_{0}}} - 3 \alpha \tau_{0 2} {\partial {\partial \phi_{1}}} - 3 \alpha \tau_{1 2} {\partial {\partial \phi_{1}}} + \frac{12 \phi_{0}^{3}}{\alpha} \tau_{0 2} - \frac{18 \phi_{1}}{\alpha} \phi_{0}^{2} \tau_{0 1} + \frac{18 \phi_{1}}{\alpha} \phi_{0}^{2} \tau_{0 2} + \frac{18 \phi_{1}}{\alpha} \phi_{0}^{2} \tau_{1 2} - \frac{18 \phi_{0}^{2}}{\alpha} \tau_{0 2} - \frac{18 \phi_{0}}{\alpha} \phi_{1}^{2} \tau_{0 1} + \frac{18 \phi_{0}}{\alpha} \phi_{1}^{2} \tau_{0 2} + \frac{18 \phi_{0}}{\alpha} \phi_{1}^{2} \tau_{1 2} + \frac{18 \phi_{0}}{\alpha} \phi_{1} \tau_{0 1} - \frac{18 \phi_{0}}{\alpha} \phi_{1} \tau_{0 2} - \frac{18 \phi_{0}}{\alpha} \phi_{1} \tau_{1 2} + \frac{6 \phi_{0}}{\alpha} \tau_{0 2} - \frac{6 \phi_{1}^{3}}{\alpha} \tau_{0 1} + \frac{6 \phi_{1}^{3}}{\alpha} \tau_{0 2} + \frac{6 \phi_{1}^{3}}{\alpha} \tau_{1 2} + \frac{9 \phi_{1}^{2}}{\alpha} \tau_{0 1} - \frac{9 \phi_{1}^{2}}{\alpha} \tau_{0 2} - \frac{9 \phi_{1}^{2}}{\alpha} \tau_{1 2} - \frac{3 \phi_{1}}{\alpha} \tau_{0 1} + \frac{3 \phi_{1}}{\alpha} \tau_{0 2} + \frac{3 \phi_{1}}{\alpha} \tau_{1 2}$$
3⋅α⋅τ₀ ₁⋅D(D(phi_1)) - 6⋅α⋅τ₀ ₂⋅D(D(phi_0)) - 3⋅α⋅τ₀ ₂⋅D(D(phi_1)) - 3⋅α⋅τ₁ ₂⋅
3 2 2 2
12⋅φ₀ ⋅τ₀ ₂ 18⋅φ₀ ⋅φ₁⋅τ₀ ₁ 18⋅φ₀ ⋅φ₁⋅τ₀ ₂ 18⋅φ₀ ⋅φ₁⋅τ₁ ₂
D(D(phi_1)) + ─────────── - ────────────── + ────────────── + ────────────── -
α α α α
2 2 2 2
18⋅φ₀ ⋅τ₀ ₂ 18⋅φ₀⋅φ₁ ⋅τ₀ ₁ 18⋅φ₀⋅φ₁ ⋅τ₀ ₂ 18⋅φ₀⋅φ₁ ⋅τ₁ ₂ 18⋅φ₀⋅φ₁⋅τ₀
─────────── - ────────────── + ────────────── + ────────────── + ────────────
α α α α α
3 3
₁ 18⋅φ₀⋅φ₁⋅τ₀ ₂ 18⋅φ₀⋅φ₁⋅τ₁ ₂ 6⋅φ₀⋅τ₀ ₂ 6⋅φ₁ ⋅τ₀ ₁ 6⋅φ₁ ⋅τ₀ ₂ 6⋅φ₁
─ - ───────────── - ───────────── + ───────── - ────────── + ────────── + ────
α α α α α
3 2 2 2
⋅τ₁ ₂ 9⋅φ₁ ⋅τ₀ ₁ 9⋅φ₁ ⋅τ₀ ₂ 9⋅φ₁ ⋅τ₁ ₂ 3⋅φ₁⋅τ₀ ₁ 3⋅φ₁⋅τ₀ ₂ 3⋅φ₁⋅τ
────── + ────────── - ────────── - ────────── - ───────── + ───────── + ──────
α α α α α α α
₁ ₂
───
%% Cell type:markdown id: tags:
Next we check, that the $tanh$ is indeed a solution of this equation...
%% Cell type:code id: tags:
``` python
x = sp.symbols("x")
expected_profile = analytic_interface_profile(x)
expected_profile
```
%% Output
$$\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}$$
⎛ x ⎞
tanh⎜───⎟
⎝2⋅α⎠ 1
───────── + ─
2 2
%% Cell type:markdown id: tags:
... by inserting it for $\phi_0$, and setting all other order parameters to zero
%% Cell type:code id: tags:
``` python
# zero other parameters
diff_eq = mu_diff_eq[0].subs({p: 0 for p in order_params[1:]})
# insert analytical solution
diff_eq = diff_eq.subs(order_params[0], expected_profile)
diff_eq.factor()
```
%% Output
$$- \frac{3 \tau_{0 2}}{2 \alpha} \left(4 \alpha^{2} {\partial {\partial (\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}) }} - \tanh^{3}{\left (\frac{x}{2 \alpha} \right )} + \tanh{\left (\frac{x}{2 \alpha} \right )}\right)$$
⎛ 2 3⎛ x ⎞ ⎛ x ⎞⎞
-3⋅τ₀ ₂⋅⎜4⋅α ⋅D(D(tanh(x/(2*alpha))/2 + 1/2)) - tanh ⎜───⎟ + tanh⎜───⎟⎟
⎝ ⎝2⋅α⎠ ⎝2⋅α⎠⎠
────────────────────────────────────────────────────────────────────────
2⋅α
%% Cell type:markdown id: tags:
finally the differentials have to be evaluated...
%% Cell type:code id: tags:
``` python
from pystencils.fd import evaluate_diffs
diff_eq = evaluate_diffs(diff_eq, x)
diff_eq
```
%% Output
$$\frac{12}{\alpha} \tau_{0 2} \left(\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}\right)^{3} - \frac{18}{\alpha} \tau_{0 2} \left(\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}\right)^{2} + \frac{6}{\alpha} \tau_{0 2} \left(\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}\right) + \frac{3 \tau_{0 2}}{2 \alpha} \left(- \tanh^{2}{\left (\frac{x}{2 \alpha} \right )} + 1\right) \tanh{\left (\frac{x}{2 \alpha} \right )}$$
3 2
⎛ ⎛ x ⎞ ⎞ ⎛ ⎛ x ⎞ ⎞ ⎛ ⎛ x ⎞ ⎞
⎜tanh⎜───⎟ ⎟ ⎜tanh⎜───⎟ ⎟ ⎜tanh⎜───⎟ ⎟
⎜ ⎝2⋅α⎠ 1⎟ ⎜ ⎝2⋅α⎠ 1⎟ ⎜ ⎝2⋅α⎠ 1⎟
12⋅τ₀ ₂⋅⎜───────── + ─⎟ 18⋅τ₀ ₂⋅⎜───────── + ─⎟ 6⋅τ₀ ₂⋅⎜───────── + ─⎟
⎝ 2 2⎠ ⎝ 2 2⎠ ⎝ 2 2⎠
──────────────────────── - ──────────────────────── + ────────────────────── +
α α α
⎛ 2⎛ x ⎞ ⎞ ⎛ x ⎞
3⋅τ₀ ₂⋅⎜- tanh ⎜───⎟ + 1⎟⋅tanh⎜───⎟
⎝ ⎝2⋅α⎠ ⎠ ⎝2⋅α⎠
───────────────────────────────────
2⋅α
%% Cell type:markdown id: tags:
.. and the result simplified...
%% Cell type:code id: tags:
``` python
assert diff_eq.expand() == 0
diff_eq.expand()
```
%% Output
$$0$$
0
%% Cell type:markdown id: tags:
...and indeed the expected tanh profile satisfies this differential equation.
Next lets check for the interface between phase 0 and phase 1:
%% Cell type:code id: tags:
``` python
for diff_eq in mu_diff_eq:
eq = diff_eq.subs({order_params[0]: expected_profile,
order_params[1]: 1 - expected_profile})
assert evaluate_diffs(eq, x).expand() == 0
```
%% Cell type:markdown id: tags:
### Checking the surface tensions parameters
Computing the excess free energy per unit area of an interface between two phases.
This should be exactly the surface tension parameter.
%% Cell type:code id: tags:
``` python
order_params = symbolic_order_parameters(num_symbols=num_phases-1)
F = free_energy_functional_n_phases(order_parameters=order_params)
```
%% Cell type:code id: tags:
``` python
two_phase_free_energy = F.subs({order_params[0]: expected_profile,
order_params[1]: 1 - expected_profile})
# Evaluate differentials and simplification
two_phase_free_energy = sp.simplify(evaluate_diffs(two_phase_free_energy, x))
```
%% Cell type:code id: tags:
``` python
excess_free_energy = sp.Integral(two_phase_free_energy, x)
excess_free_energy
```
%% Output
$$\int \frac{3 \tau_{0 1}}{8 \alpha \cosh^{4}{\left (\frac{x}{2 \alpha} \right )}}\, dx$$
⌠
⎮ 3⋅τ₀ ₁
⎮ ────────────── dx
⎮ 4⎛ x ⎞
⎮ 8⋅α⋅cosh ⎜───⎟
⎮ ⎝2⋅α⎠
⌡
%% Cell type:markdown id: tags:
Sympy cannot integrate this automatically - help with a manual substitution $\frac{x}{2\alpha} \rightarrow u$
%% Cell type:code id: tags:
``` python
coshTerm = list(excess_free_energy.atoms(sp.cosh))[0]
transformed_int = excess_free_energy.transform(coshTerm.args[0], sp.Symbol("u", real=True))
transformed_int
```
%% Output
$$\int \frac{3 \tau_{0 1}}{4 \cosh^{4}{\left (u \right )}}\, du$$
⌠
⎮ 3⋅τ₀ ₁
⎮ ────────── du
⎮ 4
⎮ 4⋅cosh (u)
⌡
%% Cell type:markdown id: tags:
Now the integral can be done:
%% Cell type:code id: tags:
``` python
result = sp.integrate(transformed_int.args[0], (transformed_int.args[1][0], -sp.oo, sp.oo))
assert result == symmetric_symbolic_surface_tension(0,1)
result
```
%% Output
$$\tau_{0 1}$$
τ₀ ₁
lbmpy_tests/phasefield/test_liquid_lens_angles.ipynb deleted 100644 → 0
View file @ cf17a9b9
Source diff could not be displayed: it is too large. Options to address this: view the blob.
lbmpy_tests/phasefield/test_n_phase_boyer_noncoupled.ipynb deleted 100644 → 0
View file @ cf17a9b9
%% Cell type:code id: tags:
``` python
import pytest
pytest.importorskip('pycuda')
```
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.phasefield.n_phase_boyer import *
from lbmpy.phasefield.kerneleqs import *
from lbmpy.phasefield.contact_angle_circle_fitting import *
from scipy.ndimage.filters import gaussian_filter
from pystencils.simp import sympy_cse_on_assignment_list
one = sp.sympify(1)
import pyximport
pyximport.install(language_level=3)
from lbmpy.phasefield.simplex_projection import simplex_projection_2d # NOQA
```
%% Cell type:markdown id: tags:
# Simulation arbitrary surface tension case
%% Cell type:code id: tags:
``` python
n = 4
dx, dt = 1, 1
mobility = 2e-3
domain_size = (150, 150)
ε = one * 4
penalty_factor = 0
stabilization_factor = 10
κ = (one, one/2, one/3, one/4)
sigma_factor = one / 15
σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
#σ
```
%% Cell type:code id: tags:
``` python
dh = create_data_handling(domain_size, periodicity=True, default_target='gpu')
c = dh.add_array('c', values_per_cell=n)
c_tmp = dh.add_array_like('c_tmp', 'c')
μ = dh.add_array('mu', values_per_cell=n)
cvec = c.center_vector
μvec = μ.center_vector
```
%% Cell type:code id: tags:
``` python
α, _ = diffusion_coefficients(σ)
f = lambda c: c**2 * ( 1 - c ) **2
a, b = compute_ab(f)
capital_f = capital_f0(cvec, σ) + correction_g(cvec, σ) + stabilization_factor * stabilization_term(cvec, α)
f_bulk = free_energy_bulk(capital_f, b, ε) + penalty_factor * (one - sum(cvec))
f_if = free_energy_interfacial(cvec, σ, a, ε)
f = f_bulk + f_if
```
%% Cell type:code id: tags:
``` python
#f_bulk
```
%% Cell type:code id: tags:
``` python
μ_assignments = mu_kernel(f, cvec, c, μ)
μ_assignments = [Assignment(a.lhs, a.rhs.doit()) for a in μ_assignments]
μ_assignments = sympy_cse_on_assignment_list(μ_assignments)
```
%% Cell type:code id: tags:
``` python
discretize = fd.Discretization2ndOrder(dx=dx, dt=dt)
```
%% Cell type:code id: tags:
``` python
def lapl(e):
return sum(ps.fd.diff(e, d, d) for d in range(dh.dim))
```
%% Cell type:code id: tags:
``` python
rhs = α * μvec
discretized_rhs = [discretize(fd.expand_diff_full( lapl(mobility * rhs_i) + fd.transient(cvec[i], idx=i), functions=μvec))
for i, rhs_i in enumerate(rhs)]
c_assignments = [Assignment(lhs, rhs)
for lhs, rhs in zip(c_tmp.center_vector, discretized_rhs)]
```
%% Cell type:code id: tags:
``` python
#c_assignments
```
%% Cell type:code id: tags:
``` python
μ_sync = dh.synchronization_function(μ.name)
c_sync = dh.synchronization_function(c.name)
optimization = {'cpu_openmp': 4, 'cpu_vectorize_info': None}
μ_kernel = create_kernel(μ_assignments, target=dh.default_target, **optimization).compile()
c_kernel = create_kernel(c_assignments, target=dh.default_target, **optimization).compile()
def set_c(slice_obj, values):
for block in dh.iterate(slice_obj):
arr = block[c.name]
arr[..., : ] = values
def smooth():
for block in dh.iterate(ghost_layers=True):
c_arr = block[c.name]
for i in range(n):
gaussian_filter(c_arr[..., i], sigma=2, output=c_arr[..., i])
def time_loop(steps):
dh.all_to_gpu()
for t in range(steps):
c_sync()
dh.run_kernel(μ_kernel)
μ_sync()
dh.run_kernel(c_kernel)
dh.swap(c.name, c_tmp.name)
#simplex_projection_2d(dh.cpu_arrays[c.name])
dh.all_to_cpu()
```
%% Cell type:code id: tags:
``` python
set_c(make_slice[:, :], [0, 0, 0, 0])
set_c(make_slice[:, 0.5:], [1, 0, 0, 0])
set_c(make_slice[:, :0.5], [0, 1, 0, 0])
set_c(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1, 0])
smooth()
```
%% Cell type:code id: tags:
``` python
#dh.load_all('n_phases_state_size200_stab10.npz')
```
%% Cell type:code id: tags:
``` python
plt.phase_plot(dh.gather_array(c.name))
```
%% Output
%% Cell type:code id: tags:
``` python
neumann_angles_from_surface_tensions(lambda i, j: float(σ[i, j]))
```
%% Output
$\displaystyle \left[ 146.44269023807928, \ 117.81813928465394, \ 95.73917047726678\right]$
[146.44269023807928, 117.81813928465394, 95.73917047726678]
%% Cell type:code id: tags:
``` python
import time
for i in range(10):
start = time.perf_counter()
time_loop(1_000)
end = time.perf_counter()
try:
print(i, end - start, liquid_lens_neumann_angles(dh.gather_array(c.name)))
except Exception:
print(i, end - start, "none found")
```
%% Output
0 0.30607624700132874 [83.97888710273061, 100.48794346625529, 175.5331694310141]
1 0.2600655169990205 [83.73094685534376, 100.65854574856168, 175.6105073960945]
2 0.2601136189987301 [83.49914818603683, 100.82173327601079, 175.67911853795226]
3 0.25987518599868054 [83.31519592224448, 100.94468140501989, 175.74012267273554]
4 0.2651959220020217 [83.14239972296966, 101.06100094405181, 175.79659933297853]
5 0.25910847799968906 [82.984481834461, 101.16731750637399, 175.8482006591651]
6 0.259863024999504 [82.84781128433397, 101.2570276449976, 175.89516107066834]
7 0.2606479199966998 [82.7456965110742, 101.31687551766585, 175.93742797125986]
8 0.25991897900166805 [82.67010885583116, 101.35099855297112, 175.97889259119805]
9 0.2590353729974595 [75.9000280154447, 108.9652166787719, 175.1347553057833]
%% Cell type:code id: tags:
``` python
plt.subplot(1,3,1)
t = dh.gather_array(c.name, make_slice[25, :]).squeeze()
plt.plot(t);
plt.subplot(1,3,2)
plt.phase_plot(dh.gather_array(c.name), linewidth=1)
plt.subplot(1,3,3)
plt.scalar_field(dh.gather_array(μ.name)[:, :, 2])
plt.colorbar();
```
%% Output
%% Cell type:code id: tags:
``` python
assert not np.isnan(dh.max(c.name))
```
%% Cell type:code id: tags:
``` python
t = dh.gather_array(c.name, make_slice[25, 55:90]).squeeze()
plt.hlines(0.5, 0, 30)
plt.plot(t);
```
%% Output
lbmpy_tests/phasefield/test_n_phase_paper_boyer_lbm.ipynb deleted 100644 → 0
View file @ cf17a9b9
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.phasefield.analytical import *
from lbmpy.phasefield.n_phase_boyer import *
from lbmpy.phasefield.kerneleqs import *
from pystencils.fd.spatial import discretize_spatial
from pystencils.simp import sympy_cse_on_assignment_list
from lbmpy.phasefield.cahn_hilliard_lbm import *
from pystencils.fd.derivation import *
one = sp.sympify(1)
```
%% Cell type:markdown id: tags:
# Chemical potential
Current state:
- not stable (yet)
- without LB coupling the model is stable
%% Cell type:code id: tags:
``` python
domain_size = (100, 100)
n = 4
c = sp.symbols(f"c_:{n}")
simple_potential = False
omega_h = 1.3
if simple_potential:
f = free_energy_functional_n_phases_penalty_term(c, interface_width=1, kappa=(0.05, 0.05/2, 0.05/4))
else:
ε = one * 4
mobility = one * 2 / 1000
κ = (one, one/2, one/3, one/4)
sigma_factor = one / 15
σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
α, _ = diffusion_coefficients(σ)
f_b, f_if, mu_b, mu_if = chemical_potential_n_phase_boyer(c, ε, σ, one,
assume_nonnegative=True, zero_threshold=1e-10)
```
%% Cell type:markdown id: tags:
# Data Setup
%% Cell type:code id: tags:
``` python
dh = create_data_handling(domain_size, periodicity=True, default_ghost_layers=2)
p_linearization = symmetric_tensor_linearization(dh.dim)
c_field = dh.add_array("c", values_per_cell=n)
rho_field = dh.add_array("rho")
mu_field = dh.add_array("mu", values_per_cell=n, latex_name="\\mu")
pressure_field = dh.add_array("P", values_per_cell=len(p_linearization))
force_field = dh.add_array("F", values_per_cell=dh.dim)
u_field = dh.add_array("u", values_per_cell=dh.dim)
# Distribution functions for each order parameter
pdf_field = []
pdf_dst_field = []
for i in range(n):
pdf_field_local = dh.add_array("f%d" % i, values_per_cell=9) # 9 for D2Q9
pdf_dst_field_local = dh.add_array("f%d_dst"%i, values_per_cell=9, latex_name="f%d_{dst}" % i)
pdf_field.append(pdf_field_local)
pdf_dst_field.append(pdf_dst_field_local)
# Distribution functions for the hydrodynamics
pdf_hydro_field = dh.add_array("fh", values_per_cell=9)
pdf_hydro_dst_field = dh.add_array("fh_dst", values_per_cell=9, latex_name="fh_{dst}")
```
%% Cell type:markdown id: tags:
### μ-Kernel
%% Cell type:code id: tags:
``` python
if simple_potential:
mu_assignments = mu_kernel(f, c, c_field, mu_field, discretization='isotropic')
else:
mu_subs = {a: b for a, b in zip(c, c_field.center_vector)}
mu_if_discrete = [discretize_spatial(e.subs(mu_subs), dx=1, stencil='isotropic') for e in mu_if]
mu_b_discrete = [e.subs(mu_subs) for e in mu_b]
mu_assignments = [Assignment(mu_field(i),
mu_if_discrete[i] + mu_b_discrete[i]) for i in range(n)]
mu_assignments = sympy_cse_on_assignment_list(mu_assignments)
μ_kernel = create_kernel(mu_assignments).compile()
```
%% Cell type:code id: tags:
``` python
second_neighbor_stencil = [(i, j)
for i in (-2, -1, 0, 1, 2)
for j in (-2, -1, 0, 1, 2)
]
x_diff = FiniteDifferenceStencilDerivation((0,), second_neighbor_stencil)
x_diff.set_weight((2, 0), sp.Rational(1, 10))
x_diff.assume_symmetric(0, anti_symmetric=True)
x_diff.assume_symmetric(1)
x_diff_stencil = x_diff.get_stencil(isotropic=True)
y_diff = FiniteDifferenceStencilDerivation((1,), second_neighbor_stencil)
y_diff.set_weight((0, 2), sp.Rational(1, 10))
y_diff.assume_symmetric(1, anti_symmetric=True)
y_diff.assume_symmetric(0)
y_diff_stencil = y_diff.get_stencil(isotropic=True)
```
%% Cell type:code id: tags:
``` python
μ = mu_field
μ_discretization = {}
for i in range(n):
μ_discretization.update({Diff(μ(i), 0): x_diff_stencil.apply(μ(i)),
Diff(μ(i), 1): y_diff_stencil.apply(μ(i))})
force_rhs = force_from_phi_and_mu(order_parameters=c_field.center_vector, dim=dh.dim, mu=mu_field.center_vector)
force_rhs = force_rhs.subs(μ_discretization)
force_assignments = [Assignment(force_field(i), force_rhs[i]) for i in range(dh.dim)]
force_kernel = create_kernel(force_assignments).compile()
```
%% Cell type:markdown id: tags:
## Lattice Boltzmann kernels
For each order parameter a Cahn-Hilliard LB method computes the time evolution.
%% Cell type:code id: tags:
``` python
if simple_potential:
mu_alpha = mu_field.center_vector
else:
mu_alpha = mobility * α * mu_field.center_vector
```
%% Cell type:code id: tags:
``` python
# Defining the Cahn-Hilliard Collision assignments
ch_kernels = []
for i in range(n):
ch_method = cahn_hilliard_lb_method(get_stencil("D2Q9"), mu_alpha[i],
relaxation_rate=1.0, gamma=1.0)
kernel = create_lb_function(lb_method=ch_method,
velocity_input=u_field.center_vector,
compressible=True,
output={'density': c_field(i)},
optimization={"symbolic_field":pdf_field[i],
"symbolic_temporary_field": pdf_dst_field[i]})
ch_kernels.append(kernel)
```
%% Cell type:code id: tags:
``` python
hydro_lbm = create_lb_function(relaxation_rate=omega_h, force=force_field,
compressible=True,
optimization={"symbolic_field": pdf_hydro_field,
"symbolic_temporary_field": pdf_hydro_dst_field},
output={'velocity': u_field}
)
#hydro_lbm.update_rule
```
%% Cell type:markdown id: tags:
# Initialization
%% Cell type:code id: tags:
``` python
def set_c(slice_obj, values):
for block in dh.iterate(slice_obj):
arr = block[c_field.name]
arr[..., : ] = values
def init():
dh.fill(u_field.name, 0)
dh.fill(mu_field.name, 0)
dh.fill(force_field.name, 0)
set_c(make_slice[:, :], [0, 0, 0, 0])
set_c(make_slice[:, 0.5:], [1, 0, 0, 0])
set_c(make_slice[:, :0.5], [0, 1, 0, 0])
set_c(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1, 0])
pdf_sync_fns = dh.synchronization_function([f.name for f in pdf_field])
hydro_sync_fn=dh.synchronization_function([pdf_hydro_field.name])
c_sync_fn = dh.synchronization_function([c_field.name])
mu_sync_fn = dh.synchronization_function([mu_field.name])
def time_loop(steps):
for i in range(steps):
c_sync_fn()
dh.run_kernel(μ_kernel)
mu_sync_fn()
dh.run_kernel(force_kernel)
hydro_sync_fn()
dh.run_kernel(hydro_lbm)
dh.swap(pdf_hydro_field.name, pdf_hydro_dst_field.name)
pdf_sync_fns()
for i in range(n):
dh.run_kernel(ch_kernels[i])
dh.swap(pdf_field[i].name,pdf_dst_field[i].name)
```
%% Cell type:code id: tags:
``` python
init()
plt.phase_plot(dh.gather_array(c_field.name))
```
%% Output
/local/bauer/miniconda3/envs/pystencils_dev/lib/python3.7/site-packages/matplotlib/contour.py:1243: UserWarning: No contour levels were found within the data range.
warnings.warn("No contour levels were found"
%% Cell type:code id: tags:
``` python
time_loop(10)
```
%% Cell type:code id: tags:
``` python
#plt.scalar_field(dh.cpu_arrays[mu_field.name][..., 0])
plt.scalar_field(dh.cpu_arrays[u_field.name][..., 0])
plt.colorbar();
```
%% Output
lbmpy_tests/phasefield/test_n_phase_penaltyterm_model.ipynb deleted 100644 → 0
View file @ cf17a9b9
Source diff could not be displayed: it is too large. Options to address this: view the blob.
lbmpy_tests/phasefield/test_numerical_1D_3phase_model.ipynb deleted 100644 → 0
View file @ cf17a9b9
Source diff could not be displayed: it is too large. Options to address this: view the blob.
lbmpy_tests/phasefield/test_numerical_1D_nphase_model.ipynb deleted 100644 → 0
View file @ cf17a9b9
Source diff could not be displayed: it is too large. Options to address this: view the blob.