Commit a759b049 authored by Frederik Hennig's avatar Frederik Hennig
Browse files

Patched SRT form of force models

parent f6b0f841
Pipeline #33895 failed with stages
in 14 minutes and 38 seconds
...@@ -240,9 +240,16 @@ class Guo(AbstractForceModel): ...@@ -240,9 +240,16 @@ class Guo(AbstractForceModel):
""" """
def __call__(self, lb_method): def __call__(self, lb_method):
if len(set(lb_method.relaxation_rates)) == 1:
# It's an SRT method!
rr = lb_method.symbolic_relaxation_matrix[0]
force_terms = Luo(self.symbolic_force_vector)(lb_method)
correction_factor = (1 - rr / 2)
result = correction_factor * force_terms
else:
force_terms = self.moment_space_forcing(lb_method) force_terms = self.moment_space_forcing(lb_method)
result = lb_method.moment_matrix.inv() * force_terms result = (lb_method.moment_matrix.inv() * force_terms).expand()
return result.expand() return result
def moment_space_forcing(self, lb_method): def moment_space_forcing(self, lb_method):
luo = Luo(self.symbolic_force_vector) luo = Luo(self.symbolic_force_vector)
...@@ -270,8 +277,7 @@ class He(AbstractForceModel): ...@@ -270,8 +277,7 @@ class He(AbstractForceModel):
Adapts the calculation of the macroscopic velocity as well as the equilibrium velocity Adapts the calculation of the macroscopic velocity as well as the equilibrium velocity
(both shifted by :math:`\frac{\mathbf{F}}{2}`)! (both shifted by :math:`\frac{\mathbf{F}}{2}`)!
Force moments are by default derived from the continuous maxwellian equilibrium Force moments are derived from the continuous maxwellian equilibrium. From the
(to disable this, set ``use_continuous_maxwellian=False`` in constructor). From the
moment integrals of the continuous force term moment integrals of the continuous force term
.. math:: .. math::
...@@ -294,9 +300,8 @@ class He(AbstractForceModel): ...@@ -294,9 +300,8 @@ class He(AbstractForceModel):
\right) \right)
""" """
def __init__(self, force, use_continuous_maxwellian=True): def __init__(self, force):
super(He, self).__init__(force) super(He, self).__init__(force)
self.continuous = use_continuous_maxwellian
def forcing_terms(self, lb_method): def forcing_terms(self, lb_method):
rho = lb_method.zeroth_order_equilibrium_moment_symbol rho = lb_method.zeroth_order_equilibrium_moment_symbol
...@@ -314,7 +319,7 @@ class He(AbstractForceModel): ...@@ -314,7 +319,7 @@ class He(AbstractForceModel):
eu_dot_f = (direction - u).dot(force) eu_dot_f = (direction - u).dot(force)
result.append(eq * eu_dot_f / cs_sq) result.append(eq * eu_dot_f / cs_sq)
return result return sp.Matrix(result)
def continuous_force_raw_moments(self, lb_method): def continuous_force_raw_moments(self, lb_method):
rho = lb_method.zeroth_order_equilibrium_moment_symbol rho = lb_method.zeroth_order_equilibrium_moment_symbol
...@@ -353,35 +358,26 @@ class He(AbstractForceModel): ...@@ -353,35 +358,26 @@ class He(AbstractForceModel):
return (shift_matrix * raw_moments).expand() return (shift_matrix * raw_moments).expand()
def __call__(self, lb_method): def __call__(self, lb_method):
if len(set(lb_method.relaxation_rates)) == 1:
# It's an SRT method!
rr = lb_method.symbolic_relaxation_matrix[0]
force_terms = self.forcing_terms(lb_method)
correction_factor = (1 - rr / 2)
result = correction_factor * force_terms
else:
force_terms = self.moment_space_forcing(lb_method) force_terms = self.moment_space_forcing(lb_method)
result = lb_method.moment_matrix.inv() * force_terms result = (lb_method.moment_matrix.inv() * force_terms).expand()
return result.expand() return result
def moment_space_forcing(self, lb_method): def moment_space_forcing(self, lb_method):
correction_factor = sp.eye(len(lb_method.stencil)) - sp.Rational(1, 2) * lb_method.symbolic_relaxation_matrix correction_factor = sp.eye(len(lb_method.stencil)) - sp.Rational(1, 2) * lb_method.symbolic_relaxation_matrix
if self.continuous:
moments = self.continuous_force_raw_moments(lb_method) moments = self.continuous_force_raw_moments(lb_method)
else:
fq = self.forcing_terms(lb_method)
# u = lb_method.first_order_equilibrium_moment_symbols
moments = lb_method.moment_matrix * sp.Matrix(fq)
# reduced_moments = [remove_higher_order_terms(moment, order=2, symbols=u) for moment in moments]
moments = (correction_factor * moments).expand() moments = (correction_factor * moments).expand()
return moments return moments
def central_moment_space_forcing(self, lb_method): def central_moment_space_forcing(self, lb_method):
correction_factor = sp.eye(len(lb_method.stencil)) - sp.Rational(1, 2) * lb_method.symbolic_relaxation_matrix correction_factor = sp.eye(len(lb_method.stencil)) - sp.Rational(1, 2) * lb_method.symbolic_relaxation_matrix
if self.continuous:
central_moments = self.continuous_force_central_moments(lb_method) central_moments = self.continuous_force_central_moments(lb_method)
else:
fq = self.forcing_terms(lb_method)
# u = lb_method.first_order_equilibrium_moment_symbols
central_moments = lb_method.moment_matrix * sp.Matrix(fq)
# reduced_central_moments = [remove_higher_order_terms(central_moment, order=2, symbols=u)
# for central_moment in central_moments]
central_moments = (correction_factor * central_moments).expand() central_moments = (correction_factor * central_moments).expand()
return central_moments return central_moments
...@@ -410,9 +406,16 @@ class Buick(AbstractForceModel): ...@@ -410,9 +406,16 @@ class Buick(AbstractForceModel):
""" """
def __call__(self, lb_method, **kwargs): def __call__(self, lb_method, **kwargs):
if len(set(lb_method.relaxation_rates)) == 1:
# It's an SRT method!
rr = lb_method.symbolic_relaxation_matrix[0]
force_terms = Simple(self.symbolic_force_vector)(lb_method)
correction_factor = (1 - rr / 2)
result = correction_factor * force_terms
else:
force_terms = self.moment_space_forcing(lb_method) force_terms = self.moment_space_forcing(lb_method)
result = lb_method.moment_matrix.inv() * force_terms result = (lb_method.moment_matrix.inv() * force_terms).expand()
return result.expand() return result
def moment_space_forcing(self, lb_method): def moment_space_forcing(self, lb_method):
simple = Simple(self.symbolic_force_vector) simple = Simple(self.symbolic_force_vector)
......
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