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
Alexander Reinauer
pystencils
Commits
1d6d60ad
Commit
1d6d60ad
authored
Mar 31, 2020
by
Alexander Reinauer
Browse files
Fixed VoF-Advection signs and more checks in testcase
parent
573f728f
Pipeline
#22928
passed with stage
in 8 minutes and 51 seconds
Changes
2
Pipelines
1
Hide whitespace changes
Inline
Side-by-side
pystencils/fd/finitevolumes.py
View file @
1d6d60ad
...
...
@@ -15,7 +15,7 @@ class FVM1stOrder:
flux: a list of sympy expressions that specify the flux, one for each cartesian direction
source: a list of sympy expressions that specify the source
"""
def
__init__
(
self
,
field
:
ps
.
field
.
Field
,
flux
=
0
,
source
=
0
):
def
__init__
(
self
,
field
:
ps
.
field
.
Field
,
flux
=
0
,
source
=
0
,
flux_scaling
:
bool
=
False
):
def
normalize
(
f
,
shape
):
shape
=
tuple
(
s
for
s
in
shape
if
s
!=
1
)
if
not
shape
:
...
...
@@ -30,6 +30,7 @@ class FVM1stOrder:
self
.
dim
=
self
.
c
.
spatial_dimensions
self
.
j
=
normalize
(
flux
,
(
self
.
dim
,
)
+
self
.
c
.
index_shape
)
self
.
q
=
normalize
(
source
,
self
.
c
.
index_shape
)
self
.
flux_scaling
=
flux_scaling
def
discrete_flux
(
self
,
flux_field
:
ps
.
field
.
Field
):
"""Return a list of assignments for the discrete fluxes
...
...
@@ -85,14 +86,22 @@ class FVM1stOrder:
for
i
in
range
(
1
,
self
.
dim
):
directional_flux
+=
fluxes
[
i
]
*
int
(
neighbor
[
i
])
discrete_flux
=
discretize
(
directional_flux
,
neighbor
)
if
self
.
flux_scaling
:
discrete_flux
/=
sp
.
Matrix
(
neighbor
).
norm
()
discrete_fluxes
.
append
(
discrete_flux
)
if
self
.
flux_scaling
:
A0
=
sum
([
sp
.
Matrix
(
ps
.
stencil
.
direction_string_to_offset
(
d
)).
norm
()
for
d
in
flux_field
.
staggered_stencil
])
/
self
.
dim
else
:
A0
=
1
if
flux_field
.
index_dimensions
>
1
:
return
[
ps
.
Assignment
(
lhs
,
rhs
)
return
[
ps
.
Assignment
(
lhs
,
rhs
/
A0
)
for
i
,
d
in
enumerate
(
flux_field
.
staggered_stencil
)
if
discrete_fluxes
[
i
]
for
lhs
,
rhs
in
zip
(
flux_field
.
staggered_vector_access
(
d
),
sp
.
simplify
(
discrete_fluxes
[
i
]))]
else
:
return
[
ps
.
Assignment
(
flux_field
.
staggered_access
(
d
),
sp
.
simplify
(
discrete_fluxes
[
i
]))
return
[
ps
.
Assignment
(
flux_field
.
staggered_access
(
d
),
sp
.
simplify
(
discrete_fluxes
[
i
]
/
A0
))
for
i
,
d
in
enumerate
(
flux_field
.
staggered_stencil
)]
def
discrete_source
(
self
):
...
...
@@ -185,13 +194,14 @@ class FVM1stOrder:
neighbors
=
flux_field
.
staggered_stencil
+
[
ps
.
stencil
.
inverse_direction_string
(
d
)
for
d
in
flux_field
.
staggered_stencil
]
divergence
=
flux_field
.
staggered_vector_access
(
neighbors
[
0
])
/
\
sp
.
Matrix
(
ps
.
stencil
.
direction_string_to_offset
(
neighbors
[
0
])).
norm
()
(
sp
.
Matrix
(
ps
.
stencil
.
direction_string_to_offset
(
neighbors
[
0
])).
norm
()
if
not
self
.
flux_scaling
else
1
)
for
d
in
neighbors
[
1
:]:
divergence
+=
flux_field
.
staggered_vector_access
(
d
)
/
\
sp
.
Matrix
(
ps
.
stencil
.
direction_string_to_offset
(
d
)).
norm
()
A0
=
sum
([
sp
.
Matrix
(
ps
.
stencil
.
direction_string_to_offset
(
d
)).
norm
()
for
d
in
flux_field
.
staggered_stencil
])
/
self
.
dim
divergence
/=
A0
(
sp
.
Matrix
(
ps
.
stencil
.
direction_string_to_offset
(
d
)).
norm
()
if
not
self
.
flux_scaling
else
1
)
if
not
self
.
flux_scaling
:
A0
=
sum
([
sp
.
Matrix
(
ps
.
stencil
.
direction_string_to_offset
(
d
)).
norm
()
for
d
in
flux_field
.
staggered_stencil
])
/
self
.
dim
divergence
/=
A0
source
=
self
.
discrete_source
()
source
=
{
s
.
lhs
:
s
.
rhs
for
s
in
source
}
...
...
@@ -217,19 +227,23 @@ def VOF(j: ps.field.Field, v: ps.field.Field, ρ: ps.field.Field):
c
=
ps
.
stencil
.
direction_string_to_offset
(
neighbor
)
v1
=
v
.
neighbor_vector
(
c
)
# sign flips if exactly one directional offset is 1
it
=
iter
(
c
==
1
)
sign
=
-
1
if
any
(
it
)
and
not
any
(
it
)
else
1
# going out
cond
=
sp
.
And
(
*
[
sp
.
Or
(
c
[
i
]
*
v0
[
i
]
>
0
,
c
[
i
]
==
0
)
for
i
in
range
(
len
(
v0
))])
overlap1
=
[
1
-
sp
.
Abs
(
v0
[
i
])
for
i
in
range
(
len
(
v0
))]
overlap2
=
[
c
[
i
]
*
v0
[
i
]
for
i
in
range
(
len
(
v0
))]
overlap2
=
[
-
v0
[
i
]
for
i
in
range
(
len
(
v0
))]
overlap
=
sp
.
Mul
(
*
[(
overlap1
[
i
]
if
c
[
i
]
==
0
else
overlap2
[
i
])
for
i
in
range
(
len
(
v0
))])
fluxes
[
d
].
append
(
ρ
.
center_vector
*
overlap
*
sp
.
Piecewise
((
1
,
cond
),
(
0
,
True
)))
fluxes
[
d
].
append
(
sign
*
ρ
.
center_vector
*
overlap
*
sp
.
Piecewise
((
1
,
cond
),
(
0
,
True
)))
# coming in
cond
=
sp
.
And
(
*
[
sp
.
Or
(
c
[
i
]
*
v1
[
i
]
<
0
,
c
[
i
]
==
0
)
for
i
in
range
(
len
(
v1
))])
overlap1
=
[
1
-
sp
.
Abs
(
v1
[
i
])
for
i
in
range
(
len
(
v1
))]
overlap2
=
[
v1
[
i
]
for
i
in
range
(
len
(
v1
))]
overlap
=
sp
.
Mul
(
*
[(
overlap1
[
i
]
if
c
[
i
]
==
0
else
overlap2
[
i
])
for
i
in
range
(
len
(
v1
))])
fluxes
[
d
].
append
(
ρ
.
neighbor_vector
(
c
)
*
overlap
*
sp
.
Piecewise
((
1
,
cond
),
(
0
,
True
)))
fluxes
[
d
].
append
(
-
sign
*
ρ
.
neighbor_vector
(
c
)
*
overlap
*
sp
.
Piecewise
((
1
,
cond
),
(
0
,
True
)))
for
i
,
ff
in
enumerate
(
fluxes
):
fluxes
[
i
]
=
ff
[
0
]
...
...
pystencils_tests/test_fvm.py
View file @
1d6d60ad
...
...
@@ -12,9 +12,11 @@ def test_advection_diffusion(dim: int):
domain_size
=
(
32
,
32
)
flux_neighbors
=
4
elif
dim
==
3
:
domain_size
=
(
16
,
16
,
16
)
domain_size
=
(
32
,
32
,
32
)
flux_neighbors
=
13
domain_center
=
[
dom
//
2
for
dom
in
domain_size
]
dh
=
ps
.
create_data_handling
(
domain_size
=
domain_size
,
periodicity
=
True
,
default_target
=
'cpu'
)
...
...
@@ -24,13 +26,13 @@ def test_advection_diffusion(dim: int):
velocity_field
=
dh
.
add_array
(
'v'
,
values_per_cell
=
dim
)
D
=
0.0666
time
=
2
00
time
=
3
00
def
grad
(
f
):
return
sp
.
Matrix
([
ps
.
fd
.
diff
(
f
,
i
)
for
i
in
range
(
dim
)])
flux_eq
=
-
D
*
grad
(
n_field
)
fvm_eq
=
ps
.
fd
.
FVM1stOrder
(
n_field
,
flux
=
flux_eq
)
fvm_eq
=
ps
.
fd
.
FVM1stOrder
(
n_field
,
flux
=
flux_eq
,
flux_scaling
=
True
)
vof_adv
=
ps
.
fd
.
VOF
(
j_field
,
velocity_field
,
n_field
)
...
...
@@ -47,20 +49,21 @@ def test_advection_diffusion(dim: int):
sync_conc
=
dh
.
synchronization_function
([
n_field
.
name
])
# analytical density calculation
# analytical density calculation
with first eight neighbors in square domain
def
density
(
pos
:
np
.
ndarray
,
time
:
int
):
return
(
4
*
np
.
pi
*
D
*
time
)
**
(
-
1.5
)
*
\
np
.
exp
(
-
np
.
sum
(
np
.
square
(
pos
),
axis
=
dim
)
/
(
4
*
D
*
time
))
return
sum
((
4
*
np
.
pi
*
D
*
time
)
**
(
-
dim
/
2
)
*
np
.
exp
(
-
np
.
sum
(
np
.
square
(
pos
+
np
.
asarray
(
shift
)),
axis
=
dim
)
/
(
4
*
D
*
time
))
for
shift
in
product
([
0
,
dh
.
shape
[
0
],
-
dh
.
shape
[
0
]],
repeat
=
dim
))
pos
=
np
.
zeros
((
*
domain_size
,
dim
))
xpos
=
np
.
arange
(
-
domain_size
[
0
]
//
2
,
domain_size
[
0
]
//
2
)
ypos
=
np
.
arange
(
-
domain_size
[
1
]
//
2
,
domain_size
[
1
]
//
2
)
if
dim
==
2
:
pos
[...,
1
],
pos
[...,
0
]
=
np
.
meshgrid
(
xpos
,
ypos
)
pos
[...,
0
],
pos
[...,
1
]
=
np
.
meshgrid
(
xpos
,
ypos
,
indexing
=
'ij'
)
elif
dim
==
3
:
zpos
=
np
.
arange
(
-
domain_size
[
2
]
//
2
,
domain_size
[
2
]
//
2
)
pos
[...,
2
],
pos
[...,
1
],
pos
[...,
0
]
=
np
.
meshgrid
(
xpos
,
ypos
,
zpos
)
pos
[...,
0
],
pos
[...,
1
],
pos
[...,
2
]
=
np
.
meshgrid
(
xpos
,
ypos
,
zpos
,
indexing
=
'ij'
)
def
run
(
velocity
:
np
.
ndarray
,
time
:
int
):
print
(
f
"
{
velocity
}
,
{
time
}
"
)
...
...
@@ -71,8 +74,7 @@ def test_advection_diffusion(dim: int):
for
i
in
range
(
dim
):
dh
.
fill
(
velocity_field
.
name
,
velocity
[
i
],
i
,
ghost_layers
=
True
,
inner_ghost_layers
=
True
)
dh
.
fill
(
n_field
.
name
,
0
)
dh
.
fill
(
n_field
.
name
,
1
,
slice_obj
=
ps
.
make_slice
[[
dom
//
2
for
dom
in
domain_size
]])
dh
.
fill
(
n_field
.
name
,
1
,
slice_obj
=
ps
.
make_slice
[
domain_center
])
sync_conc
()
for
i
in
range
(
time
):
...
...
@@ -80,12 +82,19 @@ def test_advection_diffusion(dim: int):
dh
.
run_kernel
(
pde_kernel
)
sync_conc
()
calc_density
=
density
(
pos
-
velocity
*
time
,
time
)
simulation_density
=
dh
.
gather_array
(
n_field
.
name
)
analytic_density
=
density
(
pos
-
velocity
*
time
,
time
)
peak_position
=
tuple
(
int
(
val
)
for
val
in
domain_center
+
velocity
*
time
)
analytical_maximum
=
analytic_density
[
peak_position
]
relative_maximum_error
=
np
.
abs
(
simulation_density
[
peak_position
]
-
analytical_maximum
)
/
analytical_maximum
np
.
testing
.
assert_allclose
(
dh
.
gather_array
(
n_field
.
name
),
calc_density
,
atol
=
1e-2
,
rtol
=
0
)
assert
dh
.
min
(
n_field
.
name
)
>=
0
assert
relative_maximum_error
<=
0.18
np
.
testing
.
assert_array_equal
(
np
.
unravel_index
(
np
.
argmax
(
simulation_density
),
dh
.
shape
),
peak_position
)
np
.
testing
.
assert_allclose
(
simulation_density
,
analytic_density
,
atol
=
5e-4
)
for
vel
in
product
([
0
,
-
0.0
8
,
0.0
8
],
repeat
=
dim
):
for
vel
in
product
([
0
,
-
0.0
2
,
0.0
2
],
repeat
=
dim
):
run
(
np
.
array
(
vel
),
time
)
...
...
@@ -112,7 +121,7 @@ def test_ek():
flux_eq
=
-
D
*
Gradient
(
c
)
+
D
*
z
*
c
.
center
*
Gradient
(
Phi
)
disc
=
ps
.
fd
.
FVM1stOrder
(
c
,
flux_eq
)
disc
=
ps
.
fd
.
FVM1stOrder
(
c
,
flux_eq
,
flux_scaling
=
False
)
flux_assignments
=
disc
.
discrete_flux
(
j
)
advection_assignments
=
ps
.
fd
.
VOF
(
j
,
v
,
c
)
continuity_assignments
=
disc
.
discrete_continuity
(
j
)
...
...
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