Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
itischler
lbmpy
Commits
c3dbce3c
Commit
c3dbce3c
authored
Sep 04, 2020
by
Michael Kuron
Browse files
Default to Schiller force model and test against Ladd&Verberg
parent
ec74a5cb
Changes
3
Hide whitespace changes
Inline
Side-by-side
lbmpy/creationfunctions.py
View file @
c3dbce3c
...
...
@@ -41,8 +41,8 @@ General:
compressible.
- ``equilibrium_order=2``: order in velocity, at which the equilibrium moment/cumulant approximation is
truncated. Order 2 is sufficient to approximate Navier-Stokes
- ``force_model=None``: possible values: ``None``, ``'simple'``, ``'luo'``, ``'guo'`` ``'buick'``,
or an instance of a
class implementing the same methods as the classes in :mod:`lbmpy.forcemodels`. For details, see
- ``force_model=None``: possible values: ``None``, ``'simple'``, ``'luo'``, ``'guo'``
,
``'buick'``,
``'schiller'``, or
an instance of a
class implementing the same methods as the classes in :mod:`lbmpy.forcemodels`. For details, see
:mod:`lbmpy.forcemodels`
- ``force=(0,0,0)``: either constant force or a symbolic expression depending on field value
- ``maxwellian_moments=True``: way to compute equilibrium moments/cumulants, if False the standard
...
...
@@ -387,7 +387,7 @@ def create_lb_method(**params):
no_force_model
=
'force_model'
not
in
params
or
params
[
'force_model'
]
==
'none'
or
params
[
'force_model'
]
is
None
if
not
force_is_zero
and
no_force_model
:
params
[
'force_model'
]
=
'
guo
'
params
[
'force_model'
]
=
'
schiller
'
if
'force_model'
in
params
:
force_model
=
force_model_from_string
(
params
[
'force_model'
],
params
[
'force'
][:
dim
])
...
...
lbmpy/forcemodels.py
View file @
c3dbce3c
...
...
@@ -7,7 +7,7 @@ Get started:
------------
This module offers different models to introduce a body force in the lattice Boltzmann scheme.
If you don't know which model to choose, use :class:`lbmpy.forcemodels.
Guo
`.
If you don't know which model to choose, use :class:`lbmpy.forcemodels.
Schiller
`.
For incompressible collision models the :class:`lbmpy.forcemodels.Buick` model can be better.
...
...
@@ -161,7 +161,7 @@ class Guo:
class
Schiller
:
r
"""
Force model by Schiller :cite:`schiller2008thermal`
Force model by Schiller :cite:`schiller2008thermal`
, equation 4.67
Equivalent to Guo but not restricted to SRT.
"""
def
__init__
(
self
,
force
):
...
...
lbmpy_tests/test_force.py
View file @
c3dbce3c
...
...
@@ -78,3 +78,47 @@ def test_total_momentum(method, force_model, omega):
time_loop
(
t
)
total
=
np
.
sum
(
dh
.
gather_array
(
u
.
name
),
axis
=
(
0
,
1
))
assert
np
.
allclose
(
total
/
np
.
prod
(
L
)
/
F
/
t
,
1
)
@
pytest
.
mark
.
parametrize
(
"stencil"
,
[
"D2Q9"
,
"D3Q15"
,
"D3Q19"
,
"D3Q27"
])
def
test_stress
(
stencil
):
"""check Schiller's force term in mode space"""
stencil
=
get_stencil
(
stencil
)
dim
=
len
(
stencil
[
0
])
omega_s
=
sp
.
Symbol
(
"omega_s"
)
omega_b
=
sp
.
Symbol
(
"omega_b"
)
omega_o
=
sp
.
Symbol
(
"omega_o"
)
omega_e
=
sp
.
Symbol
(
"omega_e"
)
F
=
[
sp
.
Symbol
(
f
"F_
{
i
}
"
)
for
i
in
range
(
dim
)]
method
=
create_lb_method
(
method
=
"mrt"
,
weighted
=
True
,
stencil
=
stencil
,
relaxation_rates
=
[
omega_s
,
omega_b
,
omega_o
,
omega_e
,
omega_o
,
omega_e
],
compressible
=
True
,
force_model
=
"schiller"
,
force
=
F
)
force_moments
=
sp
.
simplify
(
method
.
moment_matrix
*
sp
.
Matrix
(
method
.
force_model
(
method
)))
# The momentum modes should contain the force
assert
force_moments
[
1
:
dim
+
1
]
==
F
# The stress modes should match eq. 47 from https://doi.org/10.1023/A:1010414013942
u
=
method
.
first_order_equilibrium_moment_symbols
def
traceless
(
m
):
tr
=
sp
.
simplify
(
sp
.
Trace
(
m
))
return
m
-
tr
/
m
.
shape
[
0
]
*
sp
.
eye
(
m
.
shape
[
0
])
C
=
sp
.
Rational
(
1
,
2
)
*
(
2
+
omega_s
)
*
(
traceless
(
sp
.
Matrix
(
u
)
*
sp
.
Matrix
(
F
).
transpose
())
+
\
traceless
(
sp
.
Matrix
(
F
)
*
sp
.
Matrix
(
u
).
transpose
()))
+
\
sp
.
Rational
(
1
,
3
)
*
(
2
+
omega_b
)
*
sp
.
Matrix
(
u
).
dot
(
F
)
*
sp
.
eye
(
method
.
dim
)
num_stresses
=
(
dim
*
dim
-
dim
)
//
2
+
dim
subs
=
{
sp
.
Symbol
(
chr
(
ord
(
"x"
)
+
i
))
*
sp
.
Symbol
(
chr
(
ord
(
"x"
)
+
j
))
:
C
[
i
,
j
]
for
i
in
range
(
dim
)
for
j
in
range
(
dim
)}
for
force_moment
,
moment
in
zip
(
force_moments
[
dim
+
1
:
dim
+
1
+
num_stresses
],
method
.
moments
[
dim
+
1
:
dim
+
1
+
num_stresses
]):
ref
=
moment
.
subs
(
subs
)
diff
=
sp
.
simplify
(
ref
-
force_moment
)
assert
diff
==
0
or
isinstance
(
diff
,
sp
.
Rational
)
# difference should be zero or a constant
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment