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
47e9a1a9
Commit
47e9a1a9
authored
Sep 07, 2020
by
Michael Kuron
Browse files
Fix sign error in Schiller force model
parent
4c83f435
Changes
2
Hide whitespace changes
Inline
Side-by-side
lbmpy/forcemodels.py
View file @
47e9a1a9
...
...
@@ -175,7 +175,7 @@ class Schiller:
omega
=
get_shear_relaxation_rate
(
lb_method
)
omega_bulk
=
get_bulk_relaxation_rate
(
lb_method
)
G
=
(
u
*
force
.
transpose
()
+
force
*
u
.
transpose
()
-
uf
*
sp
.
Rational
(
2
,
lb_method
.
dim
))
*
sp
.
Rational
(
1
,
2
)
*
\
(
2
+
omega
)
+
uf
*
sp
.
Rational
(
1
,
lb_method
.
dim
)
*
(
2
+
omega_bulk
)
(
2
-
omega
)
+
uf
*
sp
.
Rational
(
1
,
lb_method
.
dim
)
*
(
2
-
omega_bulk
)
result
=
[]
for
direction
,
w_i
in
zip
(
lb_method
.
stencil
,
lb_method
.
weights
):
...
...
lbmpy_tests/test_force.py
View file @
47e9a1a9
...
...
@@ -112,15 +112,16 @@ def test_modes(stencil, force_model):
if
force_model
==
"schiller"
:
num_stresses
=
(
dim
*
dim
-
dim
)
//
2
+
dim
lambda_s
,
lambda_b
=
-
omega_s
,
-
omega_b
# The stress moments 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
+
omeg
a_s
)
*
(
traceless
(
sp
.
Matrix
(
u
)
*
sp
.
Matrix
(
F
).
transpose
())
+
\
C
=
sp
.
Rational
(
1
,
2
)
*
(
2
+
lambd
a_s
)
*
(
traceless
(
sp
.
Matrix
(
u
)
*
sp
.
Matrix
(
F
).
transpose
())
+
\
traceless
(
sp
.
Matrix
(
F
)
*
sp
.
Matrix
(
u
).
transpose
()))
+
\
sp
.
Rational
(
1
,
method
.
dim
)
*
(
2
+
omeg
a_b
)
*
sp
.
Matrix
(
u
).
dot
(
F
)
*
sp
.
eye
(
method
.
dim
)
sp
.
Rational
(
1
,
method
.
dim
)
*
(
2
+
lambd
a_b
)
*
sp
.
Matrix
(
u
).
dot
(
F
)
*
sp
.
eye
(
method
.
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
)}
...
...
@@ -139,7 +140,7 @@ def test_modes(stencil, force_model):
# Check eq. 4.53b from schiller2008thermal
assert
[
sp
.
simplify
(
sum
(
ff
[
i
]
*
stencil
[
i
][
j
]
for
i
in
range
(
len
(
stencil
))))
for
j
in
range
(
dim
)]
==
F
# Check eq. 4.61a from schiller2008thermal
ref
=
(
2
+
omeg
a_s
)
/
2
*
(
traceless
(
sp
.
Matrix
(
u
)
*
sp
.
Matrix
(
F
).
transpose
())
+
\
ref
=
(
2
+
lambd
a_s
)
/
2
*
(
traceless
(
sp
.
Matrix
(
u
)
*
sp
.
Matrix
(
F
).
transpose
())
+
\
traceless
(
sp
.
Matrix
(
F
)
*
sp
.
Matrix
(
u
).
transpose
()))
s
=
sp
.
zeros
(
dim
)
for
i
in
range
(
0
,
len
(
stencil
)):
...
...
@@ -147,7 +148,7 @@ def test_modes(stencil, force_model):
assert
sp
.
simplify
(
s
-
ref
)
==
sp
.
zeros
(
dim
)
# Check eq. 4.61b from schiller2008thermal
assert
sp
.
simplify
(
sum
(
ff
[
i
]
*
stencil
[
i
][
a
]
**
2
for
i
in
range
(
len
(
stencil
))
for
a
in
range
(
dim
))
-
(
2
+
omeg
a_b
)
*
sp
.
Matrix
(
u
).
dot
(
F
))
==
0
-
(
2
+
lambd
a_b
)
*
sp
.
Matrix
(
u
).
dot
(
F
))
==
0
# All other moments should be zero
assert
list
(
force_moments
[
dim
+
1
+
num_stresses
:])
==
[
0
]
*
(
len
(
stencil
)
-
(
dim
+
1
+
num_stresses
))
...
...
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