Commit e44b4429 authored by Markus Holzer's avatar Markus Holzer
Browse files

Fixed population space forcing

parent 1ce2c1d7
Pipeline #33593 failed with stages
in 9 minutes and 44 seconds
......@@ -537,7 +537,7 @@ def force_model_from_string(force_model_name, force_values):
'edm': forcemodels.EDM,
'schiller': forcemodels.Schiller,
'cumulant': cumulant_force_model.CenteredCumulantForceModel,
'He': forcemodels.He
'he': forcemodels.He
}
if force_model_name.lower() not in force_model_dict:
raise ValueError("Unknown force model %s" % (force_model_name,))
......
......@@ -121,7 +121,7 @@ class AbstractForceModel(abc.ABC):
def __call__(self, lb_method):
r"""
This function returns a vector of size q which is added to the PDFs in the PDF space. It has to be implemented
by all forcing models.
by all forcing models and returns a sympy Matrix containing the q dimensional force vector.
Args:
lb_method: LB method, see lbmpy.creationfunctions.create_lb_method
......@@ -169,11 +169,13 @@ class Simple(AbstractForceModel):
assert len(self._force) == lb_method.dim
cs_sq = sp.Rational(1, 3) # squared speed of sound
return [(w_i / cs_sq) * scalar_product(self._force, direction)
for direction, w_i in zip(lb_method.stencil, lb_method.weights)]
result = [(w_i / cs_sq) * scalar_product(self._force, direction)
for direction, w_i in zip(lb_method.stencil, lb_method.weights)]
return sp.Matrix(result)
def moment_space_forcing(self, lb_method):
return (lb_method.moment_matrix * sp.Matrix(self(lb_method))).expand()
return (lb_method.moment_matrix * self(lb_method)).expand()
def central_moment_space_forcing(self, lb_method):
moments = (lb_method.moment_matrix * sp.Matrix(self(lb_method))).expand()
......@@ -204,13 +206,13 @@ class Luo(AbstractForceModel):
fq = w_i * force.dot(first_summand + second_summand)
result.append(fq.simplify())
return result
return sp.Matrix(result)
def moment_space_forcing(self, lb_method):
return (lb_method.moment_matrix * sp.Matrix(self(lb_method))).expand()
return (lb_method.moment_matrix * self(lb_method)).expand()
def central_moment_space_forcing(self, lb_method):
moments = (lb_method.moment_matrix * sp.Matrix(self(lb_method))).expand()
moments = (lb_method.moment_matrix * self(lb_method)).expand()
return lb_method.shift_matrix * moments
......@@ -225,14 +227,9 @@ class Guo(AbstractForceModel):
super(Guo, self).__init__(force)
def __call__(self, lb_method):
luo = Luo(self._force)
shear_relaxation_rate = get_shear_relaxation_rate(lb_method)
assert len(set(lb_method.relaxation_rates)) == 1, (
"In population space, guo only works for SRT. Use the function moment_space_forcing as call operator "
"instead")
correction_factor = (1 - sp.Rational(1, 2) * shear_relaxation_rate)
return [correction_factor * t for t in luo(lb_method)]
force_terms = self.moment_space_forcing(lb_method)
result = lb_method.moment_matrix.inv() * force_terms
return result.expand()
def moment_space_forcing(self, lb_method):
luo = Luo(self._force)
......@@ -280,14 +277,9 @@ class He(AbstractForceModel):
return result
def __call__(self, lb_method):
fq = self.forcing_terms(lb_method)
shear_relaxation_rate = get_shear_relaxation_rate(lb_method)
assert len(set(lb_method.relaxation_rates)) == 1, (
"In population space, He only works for SRT. Use the function moment_space_forcing as call operator "
"instead")
correction_factor = (1 - sp.Rational(1, 2) * shear_relaxation_rate)
return [correction_factor * t for t in fq]
force_terms = self.moment_space_forcing(lb_method)
result = lb_method.moment_matrix.inv() * force_terms
return result.expand()
def moment_space_forcing(self, lb_method):
fq = self.forcing_terms(lb_method)
......@@ -296,7 +288,8 @@ class He(AbstractForceModel):
correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.relaxation_matrix
moments = correction_factor * (lb_method.moment_matrix * sp.Matrix(fq))
moments = moments.expand()
return [remove_higher_order_terms(moment, order=2, symbols=u) for moment in moments]
reduced_moments = [remove_higher_order_terms(moment, order=2, symbols=u) for moment in moments]
return sp.Matrix(reduced_moments)
def central_moment_space_forcing(self, lb_method):
fq = self.forcing_terms(lb_method)
......@@ -305,7 +298,9 @@ class He(AbstractForceModel):
correction_factor = sp.eye(q) - sp.Rational(1, 2) * lb_method.relaxation_matrix
moments = (lb_method.moment_matrix * sp.Matrix(fq))
central_moments = correction_factor * (lb_method.shift_matrix * moments)
return [remove_higher_order_terms(central_moment, order=2, symbols=u) for central_moment in central_moments]
reduced_central_moments = [remove_higher_order_terms(central_moment, order=2, symbols=u)
for central_moment in central_moments]
return sp.Matrix(reduced_central_moments)
def equilibrium_velocity_shift(self, density):
return default_velocity_shift(density, self._force)
......@@ -348,14 +343,9 @@ class Buick(AbstractForceModel):
super(Buick, self).__init__(force)
def __call__(self, lb_method, **kwargs):
simple = Simple(self._force)
shear_relaxation_rate = get_shear_relaxation_rate(lb_method)
assert len(set(lb_method.relaxation_rates)) == 1, (
"In population space, buick only works for SRT. Use the function moment_space_forcing as call operator "
"instead")
correction_factor = (1 - sp.Rational(1, 2) * shear_relaxation_rate)
return [correction_factor * t for t in simple(lb_method)]
force_terms = self.moment_space_forcing(lb_method)
result = lb_method.moment_matrix.inv() * force_terms
return result.expand()
def moment_space_forcing(self, lb_method):
simple = Simple(self._force)
......
......@@ -9,7 +9,7 @@ from lbmpy.methods.conservedquantitycomputation import AbstractConservedQuantity
from lbmpy.moment_transforms import (FastCentralMomentTransform,
PRE_COLLISION_CENTRAL_MOMENT, POST_COLLISION_CENTRAL_MOMENT)
from lbmpy.moments import (polynomial_to_exponent_representation, MOMENT_SYMBOLS, moment_matrix, set_up_shift_matrix,
from lbmpy.moments import (MOMENT_SYMBOLS, moment_matrix, set_up_shift_matrix,
statistical_quantity_symbol, extract_monomials, exponent_tuple_sort_key)
......@@ -246,7 +246,6 @@ class CentralMomentBasedLbMethod(AbstractLbMethod):
pre_simplification=False,
include_force_terms=False):
stencil = self.stencil
dim = len(self.stencil[0])
f = self.pre_collision_pdf_symbols
density = self.zeroth_order_equilibrium_moment_symbol
velocity = self.first_order_equilibrium_moment_symbols
......
......@@ -5,8 +5,6 @@ import lbmpy.forcemodels
from lbmpy.moments import is_bulk_moment
import pytest
from contextlib import ExitStack as does_not_raise
force_models = [fm.lower() for fm in dir(lbmpy.forcemodels) if fm[0].isupper()]
force_models.remove("abstractforcemodel")
......@@ -29,22 +27,14 @@ def test_total_momentum(method, force_model, omega):
ρ = dh.add_array('rho')
u = dh.add_array('u', values_per_cell=dh.dim)
expectation = does_not_raise()
skip = False
if force_model in ['guo', 'buick'] and method != 'srt':
expectation = pytest.raises(AssertionError)
skip = True
with expectation:
collision = create_lb_update_rule(method=method,
stencil=stencil,
relaxation_rate=omega,
compressible=True,
force_model=force_model,
force=F,
kernel_type='collide_only',
optimization={'symbolic_field': src})
if skip:
return
collision = create_lb_update_rule(method=method,
stencil=stencil,
relaxation_rate=omega,
compressible=True,
force_model=force_model,
force=F,
kernel_type='collide_only',
optimization={'symbolic_field': src})
stream = create_stream_pull_with_output_kernel(collision.method, src, dst,
{'density': ρ, 'velocity': u})
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment