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
Markus Holzer
lbmpy
Commits
4dab0828
Commit
4dab0828
authored
Mar 31, 2018
by
Martin Bauer
Browse files
Code Cleanup
- assignment collection - sympyextensions
parent
54b343b9
Changes
26
Hide whitespace changes
Inline
Side-by-side
boundaries/boundaryconditions.py
View file @
4dab0828
import
sympy
as
sp
from
pystencils
import
Field
,
Assignment
from
pystencils.astnodes
import
SympyAssignment
from
pystencils.sympyextensions
import
get
S
ymmetric
P
art
from
pystencils.sympyextensions
import
get
_s
ymmetric
_p
art
from
pystencils.data_types
import
createType
from
lbmpy.simplificationfactory
import
createSimplificationStrategy
from
lbmpy.boundaries.boundaryhandling
import
BoundaryOffsetInfo
,
LbmWeightInfo
...
...
@@ -126,7 +126,7 @@ class UBB(Boundary):
if
self
.
_adaptVelocityToForce
:
cqc
=
lbMethod
.
conservedQuantityComputation
shiftedVelEqs
=
cqc
.
equilibriumInputEquationsFromInitValues
(
velocity
=
velocity
)
velocity
=
[
eq
.
rhs
for
eq
in
shiftedVelEqs
.
extract
(
cqc
.
firstOrderMomentSymbols
).
main
A
ssignments
]
velocity
=
[
eq
.
rhs
for
eq
in
shiftedVelEqs
.
new_filtered
(
cqc
.
firstOrderMomentSymbols
).
main
_a
ssignments
]
c_s_sq
=
sp
.
Rational
(
1
,
3
)
weightOfDirection
=
LbmWeightInfo
.
weightOfDirection
...
...
@@ -140,8 +140,8 @@ class UBB(Boundary):
densitySymbol
=
sp
.
Symbol
(
"rho"
)
pdfFieldAccesses
=
[
pdfField
(
i
)
for
i
in
range
(
len
(
lbMethod
.
stencil
))]
densityEquations
=
cqc
.
outputEquationsFromPdfs
(
pdfFieldAccesses
,
{
'density'
:
densitySymbol
})
densitySymbol
=
lbMethod
.
conservedQuantityComputation
.
defined
S
ymbols
()[
'density'
]
result
=
densityEquations
.
all
Equation
s
densitySymbol
=
lbMethod
.
conservedQuantityComputation
.
defined
_s
ymbols
()[
'density'
]
result
=
densityEquations
.
all
_assignment
s
result
+=
[
Assignment
(
pdfField
[
neighbor
](
inverseDir
),
pdfField
(
direction
)
-
velTerm
*
densitySymbol
)]
return
result
...
...
@@ -159,31 +159,31 @@ class FixedDensity(Boundary):
def
__call__
(
self
,
pdfField
,
directionSymbol
,
lbMethod
,
**
kwargs
):
"""Boundary condition that fixes the density/pressure at the obstacle"""
def
removeAsymmetricPartOfmain
A
ssignments
(
eqColl
,
dofs
):
newmain
A
ssignments
=
[
Assignment
(
e
.
lhs
,
get
S
ymmetric
P
art
(
e
.
rhs
,
dofs
))
for
e
in
eqColl
.
main
A
ssignments
]
return
eqColl
.
copy
(
newmain
A
ssignments
)
def
removeAsymmetricPartOfmain
_a
ssignments
(
eqColl
,
dofs
):
newmain
_a
ssignments
=
[
Assignment
(
e
.
lhs
,
get
_s
ymmetric
_p
art
(
e
.
rhs
,
dofs
))
for
e
in
eqColl
.
main
_a
ssignments
]
return
eqColl
.
copy
(
newmain
_a
ssignments
)
neighbor
=
BoundaryOffsetInfo
.
offsetFromDir
(
directionSymbol
,
lbMethod
.
dim
)
inverseDir
=
BoundaryOffsetInfo
.
invDir
(
directionSymbol
)
cqc
=
lbMethod
.
conservedQuantityComputation
velocity
=
cqc
.
defined
S
ymbols
()[
'velocity'
]
symmetricEq
=
removeAsymmetricPartOfmain
A
ssignments
(
lbMethod
.
getEquilibrium
(),
dofs
=
velocity
)
velocity
=
cqc
.
defined
_s
ymbols
()[
'velocity'
]
symmetricEq
=
removeAsymmetricPartOfmain
_a
ssignments
(
lbMethod
.
getEquilibrium
(),
dofs
=
velocity
)
substitutions
=
{
sym
:
pdfField
(
i
)
for
i
,
sym
in
enumerate
(
lbMethod
.
preCollisionPdfSymbols
)}
symmetricEq
=
symmetricEq
.
copyW
ith
S
ubstitutions
Applied
(
substitutions
)
symmetricEq
=
symmetricEq
.
new_w
ith
_s
ubstitutions
(
substitutions
)
simplification
=
createSimplificationStrategy
(
lbMethod
)
symmetricEq
=
simplification
(
symmetricEq
)
densitySymbol
=
cqc
.
defined
S
ymbols
()[
'density'
]
densitySymbol
=
cqc
.
defined
_s
ymbols
()[
'density'
]
density
=
self
.
_density
densityEq
=
cqc
.
equilibriumInputEquationsFromInitValues
(
density
=
density
).
insertS
ubexpressions
().
main
A
ssignments
[
0
]
densityEq
=
cqc
.
equilibriumInputEquationsFromInitValues
(
density
=
density
).
new_without_s
ubexpressions
().
main
_a
ssignments
[
0
]
assert
densityEq
.
lhs
==
densitySymbol
transformedDensity
=
densityEq
.
rhs
conditions
=
[(
eq_i
.
rhs
,
sp
.
Equality
(
directionSymbol
,
i
))
for
i
,
eq_i
in
enumerate
(
symmetricEq
.
main
A
ssignments
)]
+
[(
0
,
True
)]
for
i
,
eq_i
in
enumerate
(
symmetricEq
.
main
_a
ssignments
)]
+
[(
0
,
True
)]
eq_component
=
sp
.
Piecewise
(
*
conditions
)
subExprs
=
[
Assignment
(
eq
.
lhs
,
transformedDensity
if
eq
.
lhs
==
densitySymbol
else
eq
.
rhs
)
...
...
chapman_enskog/chapman_enskog.py
View file @
4dab0828
...
...
@@ -11,8 +11,8 @@ from lbmpy.chapman_enskog.derivative import collectDerivatives, createNestedDiff
from
lbmpy.moments
import
discreteMoment
,
momentMatrix
,
polynomialToExponentRepresentation
,
getMomentIndices
,
\
nonAliasedMoment
from
pystencils.cache
import
diskcache
from
pystencils.sympyextensions
import
normalize
P
roduct
,
multidimensional
Summation
,
kronecker
D
elta
from
pystencils.sympyextensions
import
productSymmetric
from
pystencils.sympyextensions
import
normalize
_p
roduct
,
multidimensional
_sum
,
kronecker
_d
elta
from
pystencils.sympyextensions
import
symmetric_product
# --------------------------------------------- Helper Functions -------------------------------------------------------
...
...
@@ -77,7 +77,7 @@ class CeMoment(sp.Symbol):
class
LbMethodEqMoments
:
def
__init__
(
self
,
lbMethod
):
self
.
_eq
=
tuple
(
e
.
rhs
for
e
in
lbMethod
.
getEquilibrium
().
main
A
ssignments
)
self
.
_eq
=
tuple
(
e
.
rhs
for
e
in
lbMethod
.
getEquilibrium
().
main
_a
ssignments
)
self
.
_momentCache
=
dict
()
self
.
_postCollisionMomentCache
=
dict
()
self
.
_stencil
=
lbMethod
.
stencil
...
...
@@ -189,7 +189,7 @@ def takeMoments(eqn, pdfToMomentName=(('f', '\Pi'), ('\Omega f', '\\Upsilon')),
derivativeTerm
=
None
cIndices
=
[]
rest
=
1
for
factor
in
normalize
P
roduct
(
productTerm
):
for
factor
in
normalize
_p
roduct
(
productTerm
):
if
isinstance
(
factor
,
Diff
):
assert
fIndex
is
None
fIndex
=
determineFIndex
(
factor
.
getArgRecursive
())
...
...
@@ -269,13 +269,13 @@ def chainSolveAndSubstitute(eqSequence, unknownSelector, normalizingFunc=diffExp
def
countVars
(
expr
,
variables
):
factorList
=
normalize
P
roduct
(
expr
)
factorList
=
normalize
_p
roduct
(
expr
)
diffsToUnpack
=
[
e
for
e
in
factorList
if
isinstance
(
e
,
Diff
)]
factorList
=
[
e
for
e
in
factorList
if
not
isinstance
(
e
,
Diff
)]
while
diffsToUnpack
:
d
=
diffsToUnpack
.
pop
()
args
=
normalize
P
roduct
(
d
.
arg
)
args
=
normalize
_p
roduct
(
d
.
arg
)
for
a
in
args
:
if
isinstance
(
a
,
Diff
):
diffsToUnpack
.
append
(
a
)
...
...
@@ -402,7 +402,7 @@ def matchEquationsToNavierStokes(conservationEquations, rho=sp.Symbol("rho"), u=
ref
=
diffSimplify
(
Diff
(
factor
*
u
[
i
],
t
)
+
sum
(
Diff
(
factor
*
u
[
i
]
*
u
[
j
],
j
)
for
j
in
range
(
dim
)))
shearAndPressureEqs
.
append
(
diffSimplify
(
momentEqs
[
i
])
-
ref
)
#
extract
pressure term
#
new_filtered
pressure term
coefficentArgSets
=
[]
for
i
,
eq
in
enumerate
(
shearAndPressureEqs
):
coefficentArgSets
.
append
(
set
())
...
...
@@ -460,8 +460,8 @@ class ChapmanEnskogAnalysis(object):
cqc
=
method
.
conservedQuantityComputation
self
.
_method
=
method
self
.
_momentCache
=
LbMethodEqMoments
(
method
)
self
.
rho
=
cqc
.
defined
S
ymbols
(
order
=
0
)[
1
]
self
.
u
=
cqc
.
defined
S
ymbols
(
order
=
1
)[
1
]
self
.
rho
=
cqc
.
defined
_s
ymbols
(
order
=
0
)[
1
]
self
.
u
=
cqc
.
defined
_s
ymbols
(
order
=
1
)[
1
]
self
.
t
=
sp
.
Symbol
(
"t"
)
self
.
epsilon
=
sp
.
Symbol
(
"epsilon"
)
...
...
@@ -471,7 +471,7 @@ class ChapmanEnskogAnalysis(object):
# Taking moments
c
=
sp
.
Matrix
([
expandedSymbol
(
"c"
,
subscript
=
i
)
for
i
in
range
(
self
.
_method
.
dim
)])
momentsUntilOrder1
=
[
1
]
+
list
(
c
)
momentsOrder2
=
[
c_i
*
c_j
for
c_i
,
c_j
in
productSymmetric
(
c
,
c
)]
momentsOrder2
=
[
c_i
*
c_j
for
c_i
,
c_j
in
symmetric_product
(
c
,
c
)]
symbolicRelaxationRates
=
[
rr
for
rr
in
method
.
relaxationRates
if
isinstance
(
rr
,
sp
.
Symbol
)]
if
constants
is
None
:
...
...
@@ -541,7 +541,7 @@ class ChapmanEnskogAnalysis(object):
def
getShearViscosityCandidates
(
self
):
result
=
set
()
dim
=
self
.
_method
.
dim
for
i
,
j
in
productSymmetric
(
range
(
dim
),
range
(
dim
),
with
D
iagonal
=
False
):
for
i
,
j
in
symmetric_product
(
range
(
dim
),
range
(
dim
),
with
_d
iagonal
=
False
):
result
.
add
(
-
sp
.
cancel
(
self
.
_sigmaWithoutErrorTerms
[
i
,
j
]
/
(
Diff
(
self
.
u
[
i
],
j
)
+
Diff
(
self
.
u
[
j
],
i
))))
return
result
...
...
@@ -554,7 +554,7 @@ class ChapmanEnskogAnalysis(object):
# Check that shear viscosity does not depend on any u derivatives - create conditions (equations) that
# have to be fulfilled for this to be the case
viscosityReference
=
self
.
_sigmaWithoutErrorTerms
[
0
,
1
].
expand
().
coeff
(
Diff
(
self
.
u
[
0
],
1
))
for
i
,
j
in
productSymmetric
(
range
(
dim
),
range
(
dim
),
with
D
iagonal
=
False
):
for
i
,
j
in
symmetric_product
(
range
(
dim
),
range
(
dim
),
with
_d
iagonal
=
False
):
term
=
self
.
_sigmaWithoutErrorTerms
[
i
,
j
]
equalCrossTermCondition
=
sp
.
expand
(
term
.
coeff
(
Diff
(
self
.
u
[
i
],
j
))
-
viscosityReference
)
term
=
term
.
subs
({
Diff
(
self
.
u
[
i
],
j
):
0
,
...
...
@@ -656,7 +656,7 @@ class SteadyStateChapmanEnskogAnalysis(object):
moment
=
summand
.
atoms
(
CeMoment
)
moment
=
moment
.
pop
()
collisionOperatorExponent
=
normalize
P
roduct
(
summand
).
count
(
B
)
collisionOperatorExponent
=
normalize
_p
roduct
(
summand
).
count
(
B
)
if
collisionOperatorExponent
==
0
:
result
+=
summand
else
:
...
...
@@ -716,8 +716,8 @@ class SteadyStateChapmanEnskogAnalysis(object):
dim
=
self
.
method
.
dim
D
=
createNestedDiff
s
=
functools
.
partial
(
multidimensional
Summation
,
dim
=
dim
)
kd
=
kronecker
D
elta
s
=
functools
.
partial
(
multidimensional
_sum
,
dim
=
dim
)
kd
=
kronecker
_d
elta
eta
,
eta_b
=
sp
.
symbols
(
"nu nu_B"
)
u
=
sp
.
symbols
(
"u_:3"
)[:
dim
]
...
...
chapman_enskog/steady_state.py
View file @
4dab0828
import
sympy
as
sp
from
pystencils.sympyextensions
import
normalize
P
roduct
from
pystencils.sympyextensions
import
normalize
_p
roduct
from
lbmpy.chapman_enskog.derivative
import
Diff
,
DiffOperator
,
expandUsingLinearity
,
normalizeDiffOrder
from
lbmpy.chapman_enskog.chapman_enskog
import
expandedSymbol
,
useChapmanEnskogAnsatz
...
...
@@ -133,7 +133,7 @@ class SteadyStateChapmanEnskogAnalysis(object):
derivative
=
None
newProd
=
1
for
arg
in
reversed
(
normalize
P
roduct
(
product
)):
for
arg
in
reversed
(
normalize
_p
roduct
(
product
)):
if
isinstance
(
arg
,
Diff
):
assert
derivative
is
None
,
"More than one derivative term in the product"
derivative
=
arg
...
...
continuous_distribution_measures.py
View file @
4dab0828
...
...
@@ -6,18 +6,18 @@ import sympy as sp
from
lbmpy.moments
import
polynomialToExponentRepresentation
from
pystencils.cache
import
diskcache
,
memorycache
from
pystencils.sympyextensions
import
makeExponentialFuncArgumentSquares
from
pystencils.sympyextensions
import
complete_the_squares_in_exp
@
memorycache
()
def
momentGeneratingFunction
(
function
,
symbols
,
symbolsInResult
):
def
momentGeneratingFunction
(
generating_
function
,
symbols
,
symbolsInResult
):
"""
Computes the moment generating function of a probability distribution. It is defined as:
.. math ::
F[f(\mathbf{x})](\mathbf{t}) = \int e^{<\mathbf{x}, \mathbf{t}>} f(x)\; dx
:param function: sympy expression
:param
generating_
function: sympy expression
:param symbols: a sequence of symbols forming the vector x
:param symbolsInResult: a sequence forming the vector t
:return: transformation result F: an expression that depends now on symbolsInResult
...
...
@@ -33,7 +33,7 @@ def momentGeneratingFunction(function, symbols, symbolsInResult):
"""
assert
len
(
symbols
)
==
len
(
symbolsInResult
)
for
t_i
,
v_i
in
zip
(
symbolsInResult
,
symbols
):
function
*=
sp
.
exp
(
t_i
*
v_i
)
generating_
function
*=
sp
.
exp
(
t_i
*
v_i
)
# This is a custom transformation that speeds up the integrating process
# of a MaxwellBoltzmann distribution
...
...
@@ -43,11 +43,11 @@ def momentGeneratingFunction(function, symbols, symbolsInResult):
# Without this transformation the following assumptions are required for the u and v variables of Maxwell Boltzmann
# 2D: real=True ( without assumption it will not work)
# 3D: no assumption ( with assumptions it will not work )
function
=
makeExponentialFuncArgumentSquares
(
function
,
symbols
)
function
=
function
.
collect
(
symbols
)
generating_function
=
complete_the_squares_in_exp
(
generating_function
.
simplify
()
,
symbols
)
generating_function
=
generating_
function
.
collect
(
symbols
)
bounds
=
[(
s_i
,
-
sp
.
oo
,
sp
.
oo
)
for
s_i
in
symbols
]
result
=
sp
.
integrate
(
function
,
*
bounds
)
result
=
sp
.
integrate
(
generating_
function
,
*
bounds
)
return
sp
.
simplify
(
result
)
...
...
creationfunctions.py
View file @
4dab0828
...
...
@@ -192,7 +192,7 @@ def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs):
params
[
'optimizationParams'
]
=
optimizationParams
updateRule
=
createLatticeBoltzmannUpdateRule
(
**
params
)
fieldTypes
=
set
(
fa
.
field
.
dtype
for
fa
in
updateRule
.
defined
S
ymbols
if
isinstance
(
fa
,
Field
.
Access
))
fieldTypes
=
set
(
fa
.
field
.
dtype
for
fa
in
updateRule
.
defined
_s
ymbols
if
isinstance
(
fa
,
Field
.
Access
))
res
=
createKernel
(
updateRule
,
target
=
optParams
[
'target'
],
dataType
=
collateTypes
(
fieldTypes
),
cpuOpenMP
=
optParams
[
'openMP'
],
cpuVectorizeInfo
=
optParams
[
'vectorization'
],
gpuIndexing
=
optParams
[
'gpuIndexing'
],
gpuIndexingParams
=
optParams
[
'gpuIndexingParams'
],
...
...
@@ -205,17 +205,17 @@ def createLatticeBoltzmannAst(updateRule=None, optimizationParams={}, **kwargs):
@
diskcacheNoFallback
def
createLatticeBoltzmannUpdateRule
(
collisionRule
=
None
,
optimizationParams
=
{},
**
kwargs
):
params
,
opt
P
arams
=
updateWithDefaultParameters
(
kwargs
,
optimizationParams
)
params
,
opt
_p
arams
=
updateWithDefaultParameters
(
kwargs
,
optimizationParams
)
if
collisionRule
is
None
:
collisionRule
=
createLatticeBoltzmannCollisionRule
(
**
params
,
optimizationParams
=
opt
P
arams
)
collisionRule
=
createLatticeBoltzmannCollisionRule
(
**
params
,
optimizationParams
=
opt
_p
arams
)
lb
M
ethod
=
collisionRule
.
method
lb
_m
ethod
=
collisionRule
.
method
if
params
[
'output'
]
and
params
[
'kernelType'
]
==
'streamPullCollide'
:
cqc
=
lb
M
ethod
.
conservedQuantityComputation
output
E
qs
=
cqc
.
outputEquationsFromPdfs
(
lb
M
ethod
.
preCollisionPdfSymbols
,
params
[
'output'
])
collisionRule
=
collisionRule
.
merge
(
output
E
qs
)
cqc
=
lb
_m
ethod
.
conservedQuantityComputation
output
_e
qs
=
cqc
.
outputEquationsFromPdfs
(
lb
_m
ethod
.
preCollisionPdfSymbols
,
params
[
'output'
])
collisionRule
=
collisionRule
.
new_
merge
d
(
output
_e
qs
)
if
params
[
'entropic'
]:
if
params
[
'entropicNewtonIterations'
]:
...
...
@@ -228,33 +228,33 @@ def createLatticeBoltzmannUpdateRule(collisionRule=None, optimizationParams={},
else
:
collisionRule
=
addEntropyCondition
(
collisionRule
,
omegaOutputField
=
params
[
'omegaOutputField'
])
elif
params
[
'smagorinsky'
]:
smagorinsky
C
onstant
=
0.12
if
params
[
'smagorinsky'
]
is
True
else
params
[
'smagorinsky'
]
collisionRule
=
addSmagorinskyModel
(
collisionRule
,
smagorinsky
C
onstant
,
smagorinsky
_c
onstant
=
0.12
if
params
[
'smagorinsky'
]
is
True
else
params
[
'smagorinsky'
]
collisionRule
=
addSmagorinskyModel
(
collisionRule
,
smagorinsky
_c
onstant
,
omegaOutputField
=
params
[
'omegaOutputField'
])
field
D
type
=
'float64'
if
opt
P
arams
[
'doublePrecision'
]
else
'float32'
field
_data_
type
=
'float64'
if
opt
_p
arams
[
'doublePrecision'
]
else
'float32'
if
opt
P
arams
[
'symbolicField'
]
is
not
None
:
src
F
ield
=
opt
P
arams
[
'symbolicField'
]
elif
opt
P
arams
[
'fieldSize'
]:
field
S
ize
=
[
s
+
2
for
s
in
opt
P
arams
[
'fieldSize'
]]
+
[
len
(
collisionRule
.
stencil
)]
src
F
ield
=
Field
.
createFixedSize
(
params
[
'fieldName'
],
field
S
ize
,
indexDimensions
=
1
,
layout
=
opt
P
arams
[
'fieldLayout'
],
dtype
=
field
D
type
)
if
opt
_p
arams
[
'symbolicField'
]
is
not
None
:
src
_f
ield
=
opt
_p
arams
[
'symbolicField'
]
elif
opt
_p
arams
[
'fieldSize'
]:
field
_s
ize
=
[
s
+
2
for
s
in
opt
_p
arams
[
'fieldSize'
]]
+
[
len
(
collisionRule
.
stencil
)]
src
_f
ield
=
Field
.
createFixedSize
(
params
[
'fieldName'
],
field
_s
ize
,
indexDimensions
=
1
,
layout
=
opt
_p
arams
[
'fieldLayout'
],
dtype
=
field
_data_
type
)
else
:
src
F
ield
=
Field
.
createGeneric
(
params
[
'fieldName'
],
spatialDimensions
=
collisionRule
.
method
.
dim
,
indexDimensions
=
1
,
layout
=
opt
P
arams
[
'fieldLayout'
],
dtype
=
field
D
type
)
src
_f
ield
=
Field
.
createGeneric
(
params
[
'fieldName'
],
spatialDimensions
=
collisionRule
.
method
.
dim
,
indexDimensions
=
1
,
layout
=
opt
_p
arams
[
'fieldLayout'
],
dtype
=
field
_data_
type
)
dst
F
ield
=
src
F
ield
.
newFieldWithDifferentName
(
params
[
'secondFieldName'
])
dst
_f
ield
=
src
_f
ield
.
newFieldWithDifferentName
(
params
[
'secondFieldName'
])
if
params
[
'kernelType'
]
==
'streamPullCollide'
:
accessor
=
StreamPullTwoFieldsAccessor
if
any
(
opt
P
arams
[
'builtinPeriodicity'
]):
accessor
=
PeriodicTwoFieldsAccessor
(
opt
P
arams
[
'builtinPeriodicity'
],
ghostLayers
=
1
)
return
createLBMKernel
(
collisionRule
,
src
F
ield
,
dst
F
ield
,
accessor
)
if
any
(
opt
_p
arams
[
'builtinPeriodicity'
]):
accessor
=
PeriodicTwoFieldsAccessor
(
opt
_p
arams
[
'builtinPeriodicity'
],
ghostLayers
=
1
)
return
createLBMKernel
(
collisionRule
,
src
_f
ield
,
dst
_f
ield
,
accessor
)
elif
params
[
'kernelType'
]
==
'collideOnly'
:
return
createLBMKernel
(
collisionRule
,
src
F
ield
,
src
F
ield
,
CollideOnlyInplaceAccessor
)
return
createLBMKernel
(
collisionRule
,
src
_f
ield
,
src
_f
ield
,
CollideOnlyInplaceAccessor
)
elif
params
[
'kernelType'
]
==
'streamPullOnly'
:
return
createStreamPullWithOutputKernel
(
lb
M
ethod
,
src
F
ield
,
dst
F
ield
,
params
[
'output'
])
return
createStreamPullWithOutputKernel
(
lb
_m
ethod
,
src
_f
ield
,
dst
_f
ield
,
params
[
'output'
])
else
:
raise
ValueError
(
"Invalid value of parameter 'kernelType'"
,
params
[
'kernelType'
])
...
...
cumulants.py
View file @
4dab0828
...
...
@@ -9,8 +9,8 @@ from lbmpy.continuous_distribution_measures import multiDifferentiation
from
lbmpy.moments
import
momentsUpToComponentOrder
from
pystencils.cache
import
memorycache
# ------------------------------------------- Internal Functions -------------------------------------------------------
from
pystencils.sympyextensions
import
fast
S
ubs
from
pystencils.sympyextensions
import
scalar
P
roduct
from
pystencils.sympyextensions
import
fast
_s
ubs
from
pystencils.sympyextensions
import
scalar
_p
roduct
def
__getIndexedSymbols
(
passedSymbols
,
prefix
,
indices
):
...
...
@@ -79,7 +79,7 @@ def __cumulantRawMomentTransform(index, dependentVarDict, outerFunction, default
partitionList
.
append
(
i
)
if
len
(
partitionList
)
==
0
:
# special case for zero index
return
fast
S
ubs
(
outerFunction
(
zerothMoment
),
subsDict
)
return
fast
_s
ubs
(
outerFunction
(
zerothMoment
),
subsDict
)
# implementation of Faa di Bruno's formula:
result
=
0
...
...
@@ -98,14 +98,14 @@ def __cumulantRawMomentTransform(index, dependentVarDict, outerFunction, default
index
[
i
]
=
1
result
=
result
.
subs
(
createMomentSymbol
(
index
),
0
)
return
fast
S
ubs
(
result
,
subsDict
)
return
fast
_s
ubs
(
result
,
subsDict
)
@
memorycache
(
maxsize
=
16
)
def
__getDiscreteCumulantGeneratingFunction
(
function
,
stencil
,
waveNumbers
):
assert
len
(
stencil
)
==
len
(
function
)
laplaceTrafo
=
sum
([
factor
*
sp
.
exp
(
scalar
P
roduct
(
waveNumbers
,
e
))
for
factor
,
e
in
zip
(
function
,
stencil
)])
laplaceTrafo
=
sum
([
factor
*
sp
.
exp
(
scalar
_p
roduct
(
waveNumbers
,
e
))
for
factor
,
e
in
zip
(
function
,
stencil
)])
return
sp
.
ln
(
laplaceTrafo
)
...
...
innerloopsplit.py
View file @
4dab0828
...
...
@@ -16,7 +16,7 @@ def createLbmSplitGroups(lbmCollisionEqs):
Required simplification hints:
- velocity: sequence of velocity symbols
"""
sh
=
lbmCollisionEqs
.
simplification
H
ints
sh
=
lbmCollisionEqs
.
simplification
_h
ints
assert
'velocity'
in
sh
,
"Needs simplification hint 'velocity': Sequence of velocity symbols"
preCollisionSymbols
=
set
(
lbmCollisionEqs
.
method
.
preCollisionPdfSymbols
)
...
...
@@ -26,10 +26,10 @@ def createLbmSplitGroups(lbmCollisionEqs):
stencil
=
lbmCollisionEqs
.
method
.
stencil
importantSubExpressions
=
{
e
.
lhs
for
e
in
lbmCollisionEqs
.
subexpressions
if
preCollisionSymbols
.
intersection
(
lbmCollisionEqs
.
getD
ependent
S
ymbols
([
e
.
lhs
]))}
if
preCollisionSymbols
.
intersection
(
lbmCollisionEqs
.
d
ependent
_s
ymbols
([
e
.
lhs
]))}
otherWrittenFields
=
[]
for
eq
in
lbmCollisionEqs
.
main
A
ssignments
:
for
eq
in
lbmCollisionEqs
.
main
_a
ssignments
:
if
eq
.
lhs
not
in
postCollisionSymbols
and
isinstance
(
eq
.
lhs
,
Field
.
Access
):
otherWrittenFields
.
append
(
eq
.
lhs
)
if
eq
.
lhs
not
in
nonCenterPostCollisionSymbols
:
...
...
@@ -44,7 +44,7 @@ def createLbmSplitGroups(lbmCollisionEqs):
directionGroups
=
defaultdict
(
list
)
dim
=
len
(
stencil
[
0
])
for
direction
,
eq
in
zip
(
stencil
,
lbmCollisionEqs
.
main
A
ssignments
):
for
direction
,
eq
in
zip
(
stencil
,
lbmCollisionEqs
.
main
_a
ssignments
):
if
direction
==
tuple
([
0
]
*
dim
):
splitGroups
[
0
].
append
(
eq
.
lhs
)
continue
...
...
@@ -57,5 +57,5 @@ def createLbmSplitGroups(lbmCollisionEqs):
directionGroups
[
direction
].
append
(
eq
.
lhs
)
splitGroups
+=
directionGroups
.
values
()
lbmCollisionEqs
.
simplification
H
ints
[
'splitGroups'
]
=
splitGroups
lbmCollisionEqs
.
simplification
_h
ints
[
'splitGroups'
]
=
splitGroups
return
lbmCollisionEqs
lbstep.py
View file @
4dab0828
...
...
@@ -251,8 +251,8 @@ class LatticeBoltzmannStep:
inpEqs
=
cqc
.
equilibriumInputEquationsFromInitValues
(
rhoField
,
[
velField
(
i
)
for
i
in
range
(
D
)])
setterEqs
=
lbMethod
.
getEquilibrium
(
conservedQuantityEquations
=
inpEqs
)
setterEqs
=
setterEqs
.
copyW
ith
S
ubstitutions
Applied
({
sym
:
pdfField
(
i
)
for
i
,
sym
in
enumerate
(
lbMethod
.
postCollisionPdfSymbols
)})
setterEqs
=
setterEqs
.
new_w
ith
_s
ubstitutions
({
sym
:
pdfField
(
i
)
for
i
,
sym
in
enumerate
(
lbMethod
.
postCollisionPdfSymbols
)})
setterEqs
=
createSimplificationStrategy
(
lbMethod
)(
setterEqs
)
setterKernel
=
createKernel
(
setterEqs
,
target
=
'cpu'
).
compile
()
...
...
macroscopic_value_kernels.py
View file @
4dab0828
...
...
@@ -53,7 +53,7 @@ def compileMacroscopicValuesGetter(lbMethod, outputQuantities, pdfArr=None, fiel
stencil
=
lbMethod
.
stencil
pdfSymbols
=
[
pdfField
(
i
)
for
i
in
range
(
len
(
stencil
))]
eqs
=
cqc
.
outputEquationsFromPdfs
(
pdfSymbols
,
outputMapping
).
all
Equation
s
eqs
=
cqc
.
outputEquationsFromPdfs
(
pdfSymbols
,
outputMapping
).
all
_assignment
s
if
target
==
'cpu'
:
import
pystencils.cpu
as
cpu
...
...
@@ -118,10 +118,10 @@ def compileMacroscopicValuesSetter(lbMethod, quantitiesToSet, pdfArr=None, field
simplification
=
createSimplificationStrategy
(
lbMethod
)
eq
=
simplification
(
eq
)
else
:
eq
=
eq
.
insertS
ubexpressions
()
eq
=
eq
.
new_without_s
ubexpressions
()
substitutions
=
{
sym
:
pdfField
(
i
)
for
i
,
sym
in
enumerate
(
lbMethod
.
postCollisionPdfSymbols
)}
eq
=
eq
.
copyW
ith
S
ubstitutions
Applied
(
substitutions
).
all
Equation
s
eq
=
eq
.
new_w
ith
_s
ubstitutions
(
substitutions
).
all
_assignment
s
if
target
==
'cpu'
:
import
pystencils.cpu
as
cpu
...
...
@@ -146,18 +146,18 @@ def createAdvancedVelocitySetterCollisionRule(lbMethod, velocityArray, velocityR
velocityField
=
Field
.
createFromNumpyArray
(
'velInput'
,
velocityArray
,
indexDimensions
=
1
)
cqc
=
lbMethod
.
conservedQuantityComputation
densitySymbol
=
cqc
.
defined
S
ymbols
(
order
=
0
)[
1
]
velocitySymbols
=
cqc
.
defined
S
ymbols
(
order
=
1
)[
1
]
densitySymbol
=
cqc
.
defined
_s
ymbols
(
order
=
0
)[
1
]
velocitySymbols
=
cqc
.
defined
_s
ymbols
(
order
=
1
)[
1
]
# density is computed from pdfs
eqInputFromPdfs
=
cqc
.
equilibriumInputEquationsFromPdfs
(
lbMethod
.
preCollisionPdfSymbols
)
eqInputFromPdfs
=
eqInputFromPdfs
.
extract
([
densitySymbol
])
eqInputFromPdfs
=
eqInputFromPdfs
.
new_filtered
([
densitySymbol
])
# velocity is read from input field
velSymbols
=
[
velocityField
(
i
)
for
i
in
range
(
lbMethod
.
dim
)]
eqInputFromField
=
cqc
.
equilibriumInputEquationsFromInitValues
(
velocity
=
velSymbols
)
eqInputFromField
=
eqInputFromField
.
extract
(
velocitySymbols
)
eqInputFromField
=
eqInputFromField
.
new_filtered
(
velocitySymbols
)
# then both are merged together
eqInput
=
eqInputFromPdfs
.
merge
(
eqInputFromField
)
eqInput
=
eqInputFromPdfs
.
new_
merge
d
(
eqInputFromField
)
# set first order relaxation rate
lbMethod
=
deepcopy
(
lbMethod
)
...
...
maxwellian_equilibrium.py
View file @
4dab0828
...
...
@@ -157,7 +157,7 @@ def getMomentsOfContinuousMaxwellianEquilibrium(moments, dim, rho=sp.Symbol("rho
>>> getMomentsOfContinuousMaxwellianEquilibrium( ( (0,0,0), (1,0,0), (0,1,0), (0,0,1), (2,0,0) ), dim=3 )
[rho, rho*u_0, rho*u_1, rho*u_2, rho*(c_s**2 + u_0**2)]
"""
from
pystencils.sympyextensions
import
remove
H
igher
O
rder
T
erms
from
pystencils.sympyextensions
import
remove
_h
igher
_o
rder
_t
erms
from
lbmpy.moments
import
MOMENT_SYMBOLS
from
lbmpy.continuous_distribution_measures
import
continuousMoment
...
...
@@ -167,7 +167,7 @@ def getMomentsOfContinuousMaxwellianEquilibrium(moments, dim, rho=sp.Symbol("rho
mb
=
continuousMaxwellianEquilibrium
(
dim
,
rho
,
u
,
MOMENT_SYMBOLS
[:
dim
],
c_s_sq_helper
)
result
=
[
continuousMoment
(
mb
,
moment
,
MOMENT_SYMBOLS
[:
dim
]).
subs
(
c_s_sq_helper
,
c_s_sq
)
for
moment
in
moments
]
if
order
is
not
None
:
result
=
[
remove
H
igher
O
rder
T
erms
(
r
,
order
,
u
)
for
r
in
result
]
result
=
[
remove
_h
igher
_o
rder
_t
erms
(
r
,
order
=
order
,
symbols
=
u
)
for
r
in
result
]
return
result
...
...
@@ -219,7 +219,7 @@ def getCumulantsOfContinuousMaxwellianEquilibrium(cumulants, dim, rho=sp.Symbol(
c_s_sq
=
sp
.
Symbol
(
"c_s"
)
**
2
,
order
=
None
):
from
lbmpy.moments
import
MOMENT_SYMBOLS
from
lbmpy.continuous_distribution_measures
import
continuousCumulant
from
pystencils.sympyextensions
import
remove
H
igher
O
rder
T
erms
from
pystencils.sympyextensions
import
remove
_h
igher
_o
rder
_t
erms
# trick to speed up sympy integration (otherwise it takes multiple minutes, or aborts):
# use a positive, real symbol to represent c_s_sq -> then replace this symbol afterwards with the real c_s_sq
...
...
@@ -227,7 +227,7 @@ def getCumulantsOfContinuousMaxwellianEquilibrium(cumulants, dim, rho=sp.Symbol(
mb
=
continuousMaxwellianEquilibrium
(
dim
,
rho
,
u
,
MOMENT_SYMBOLS
[:
dim
],
c_s_sq_helper
)
result
=
[
continuousCumulant
(
mb
,
cumulant
,
MOMENT_SYMBOLS
[:
dim
]).
subs
(
c_s_sq_helper
,
c_s_sq
)
for
cumulant
in
cumulants
]
if
order
is
not
None
:
result
=
[
remove
H
igher
O
rder
T
erms
(
r
,
order
,
u
)
for
r
in
result
]
result
=
[
remove
_h
igher
_o
rder
_t
erms
(
r
,
order
=
order
,
symbols
=
u
)
for
r
in
result
]
return
result
...
...
methods/abstractlbmethod.py
View file @
4dab0828
...
...
@@ -8,9 +8,9 @@ RelaxationInfo = namedtuple('RelaxationInfo', ['equilibriumValue', 'relaxationRa
class
LbmCollisionRule
(
AssignmentCollection
):
def
__init__
(
self
,
lbm
M
ethod
,
*
args
,
**
kwargs
):
def
__init__
(
self
,
lb
_
method
,
*
args
,
**
kwargs
):
super
(
LbmCollisionRule
,
self
).
__init__
(
*
args
,
**
kwargs
)
self
.
method
=
lbm
M
ethod
self
.
method
=
lb
_
method
class
AbstractLbMethod
(
abc
.
ABCMeta
(
'ABC'
,
(
object
,),
{})):
...
...
methods/conservedquantitycomputation.py
View file @
4dab0828
...
...
@@ -32,7 +32,7 @@ class AbstractConservedQuantityComputation(abc.ABCMeta('ABC', (object,), {})):
and :func:`equilibriumInputEquationsFromInitValues`
"""
def
defined
S
ymbols
(
self
,
order
=
'all'
):
def
defined
_s
ymbols
(
self
,
order
=
'all'
):
"""
Returns a dict, mapping names of conserved quantities to their symbols
"""
...
...
@@ -98,7 +98,7 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
def
compressible
(
self
):
return
self
.
_compressible
def
defined
S
ymbols
(
self
,
order
=
'all'
):
def
defined
_s
ymbols
(
self
,
order
=
'all'
):
if
order
==
'all'
:
return
{
'density'
:
self
.
_symbolOrder0
,
'velocity'
:
self
.
_symbolsOrder1
}
...
...
@@ -171,7 +171,7 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
ac
=
applyForceModelShift
(
'macroscopicVelocityShift'
,
dim
,
ac
,
self
.
_forceModel
,
self
.
_compressible
)
main_assignments
=
[]
eqs
=
OrderedDict
([(
eq
.
lhs
,
eq
.
rhs
)
for
eq
in
ac
.
all
Equation
s
])
eqs
=
OrderedDict
([(
eq
.
lhs
,
eq
.
rhs
)
for
eq
in
ac
.
all
_assignment
s
])
if
'density'
in
outputQuantityNamesToSymbols
:
density_output_symbol
=
outputQuantityNamesToSymbols
[
'density'
]
...
...
@@ -203,11 +203,11 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
self
.
_symbolOrder0
,
self
.
_symbolsOrder1
)
mom_density_eq_coll
=
applyForceModelShift
(
'macroscopicVelocityShift'
,
dim
,
mom_density_eq_coll
,
self
.
_forceModel
,
self
.
_compressible
)
for
sym
,
val
in
zip
(
momentum_density_output_symbols
,
mom_density_eq_coll
.
main
A
ssignments
[
1
:]):
for
sym
,
val
in
zip
(
momentum_density_output_symbols
,
mom_density_eq_coll
.
main
_a
ssignments
[
1
:]):
main_assignments
.
append
(
Assignment
(
sym
,
val
.
rhs
))
ac
=
ac
.
copy
(
main_assignments
,
[
Assignment
(
a
,
b
)
for
a
,
b
in
eqs
.
items
()])
return
ac
.
new
W
ithout
U
nused
S
ubexpressions
()
return
ac
.
new
_w
ithout
_u
nused
_s
ubexpressions
()
def
__repr__
(
self
):
return
"ConservedValueComputation for %s"
%
(
", "
.
join
(
self
.
conservedQuantities
.
keys
()),)
...
...
@@ -278,10 +278,10 @@ def divideFirstOrderMomentsByRho(assignment_collection, dim):
Returns a new equation collection where the u terms (first order moments) are divided by rho.
The dim parameter specifies the number of first order moments. All subsequent equations are just copied over.
"""
old
E
qs
=
assignment_collection
.
main
A
ssignments
rho
=
old
E
qs
[
0
].
lhs
new_first_order_moment_eq
=
[
Assignment
(
eq
.
lhs
,
eq
.
rhs
/
rho
)
for
eq
in
old
E
qs
[
1
:
dim
+
1
]]
new_eqs
=
[
old
E
qs
[
0
]]
+
new_first_order_moment_eq
+
old
E
qs
[
dim
+
1
:]
old
_e
qs
=
assignment_collection
.
main
_a
ssignments
rho
=
old
_e
qs
[
0
].
lhs
new_first_order_moment_eq
=
[
Assignment
(
eq
.
lhs
,
eq
.
rhs
/
rho
)
for
eq
in
old
_e
qs
[
1
:
dim
+
1
]]
new_eqs
=
[
old
_e
qs
[
0
]]
+
new_first_order_moment_eq
+
old
_e
qs
[
dim
+
1
:]
return
assignment_collection
.
copy
(
new_eqs
)
...
...
@@ -289,9 +289,9 @@ def addDensityOffset(assignment_collection, offset=sp.Rational(1, 1)):
"""
Assumes that first equation is the density (zeroth moment). Changes the density equations by adding offset to it.
"""
old
E
qs
=
assignment_collection
.
main
A
ssignments
new
D
ensity
=
Assignment
(
old
E
qs
[
0
].
lhs
,
old
E
qs
[
0
].
rhs
+
offset
)
return
assignment_collection
.
copy
([
new
D
ensity
]
+
old
E
qs
[
1
:])
old
_e
qs
=
assignment_collection
.
main
_a
ssignments
new
_d
ensity
=
Assignment
(
old
_e
qs
[
0
].
lhs
,
old
_e
qs
[
0
].
rhs
+
offset
)
return
assignment_collection
.
copy
([
new
_d
ensity
]
+
old
_e
qs
[
1
:])
def
applyForceModelShift
(
shiftMemberName
,
dim
,
assignment_collection
,
forceModel
,
compressible
,
reverse
=
False
):
...
...
@@ -301,7 +301,7 @@ def applyForceModelShift(shiftMemberName, dim, assignment_collection, forceModel
equation collection are assumed to be the velocity equations.
"""
if
forceModel
is
not
None
and
hasattr
(
forceModel
,
shiftMemberName
):
old_eqs
=
assignment_collection
.
main
A
ssignments
old_eqs
=
assignment_collection
.
main
_a
ssignments
density
=
old_eqs
[
0
].
lhs
if
compressible
else
sp
.
Rational
(
1
,
1
)
old_vel_eqs
=
old_eqs
[
1
:
dim
+
1
]
shift_func
=
getattr
(
forceModel
,
shiftMemberName
)
...
...
methods/creationfunctions.py
View file @
4dab0828
...
...
@@ -11,7 +11,7 @@ from lbmpy.stencils import stencilsHaveSameEntries, getStencil
from
lbmpy.moments
import
isEven
,
gramSchmidt
,
getDefaultMomentSetForStencil
,
MOMENT_SYMBOLS
,
\
exponentsToPolynomialRepresentations
,
momentsOfOrder
,
momentsUpToComponentOrder
,
sortMomentsIntoGroupsOfSameOrder
,
\
getOrder
,
discreteMoment