Commit 53f369e9 authored by Alexander Reinauer's avatar Alexander Reinauer
Browse files

Fixed VoF-Advection signs and more checks in testcase

parent 573f728f
Pipeline #22927 canceled with stage
in 1 minute and 35 seconds
...@@ -15,7 +15,7 @@ class FVM1stOrder: ...@@ -15,7 +15,7 @@ class FVM1stOrder:
flux: a list of sympy expressions that specify the flux, one for each cartesian direction 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 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): def normalize(f, shape):
shape = tuple(s for s in shape if s != 1) shape = tuple(s for s in shape if s != 1)
if not shape: if not shape:
...@@ -30,6 +30,7 @@ class FVM1stOrder: ...@@ -30,6 +30,7 @@ class FVM1stOrder:
self.dim = self.c.spatial_dimensions self.dim = self.c.spatial_dimensions
self.j = normalize(flux, (self.dim, ) + self.c.index_shape) self.j = normalize(flux, (self.dim, ) + self.c.index_shape)
self.q = normalize(source, 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): def discrete_flux(self, flux_field: ps.field.Field):
"""Return a list of assignments for the discrete fluxes """Return a list of assignments for the discrete fluxes
...@@ -85,14 +86,22 @@ class FVM1stOrder: ...@@ -85,14 +86,22 @@ class FVM1stOrder:
for i in range(1, self.dim): for i in range(1, self.dim):
directional_flux += fluxes[i] * int(neighbor[i]) directional_flux += fluxes[i] * int(neighbor[i])
discrete_flux = discretize(directional_flux, neighbor) discrete_flux = discretize(directional_flux, neighbor)
if self.flux_scaling:
discrete_flux /= sp.Matrix(neighbor).norm()
discrete_fluxes.append(discrete_flux) 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: 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 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]))] for lhs, rhs in zip(flux_field.staggered_vector_access(d), sp.simplify(discrete_fluxes[i]))]
else: 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)] for i, d in enumerate(flux_field.staggered_stencil)]
def discrete_source(self): def discrete_source(self):
...@@ -185,13 +194,14 @@ class FVM1stOrder: ...@@ -185,13 +194,14 @@ class FVM1stOrder:
neighbors = flux_field.staggered_stencil + [ps.stencil.inverse_direction_string(d) neighbors = flux_field.staggered_stencil + [ps.stencil.inverse_direction_string(d)
for d in flux_field.staggered_stencil] for d in flux_field.staggered_stencil]
divergence = flux_field.staggered_vector_access(neighbors[0]) / \ 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:]: for d in neighbors[1:]:
divergence += flux_field.staggered_vector_access(d) / \ divergence += flux_field.staggered_vector_access(d) / \
sp.Matrix(ps.stencil.direction_string_to_offset(d)).norm() (sp.Matrix(ps.stencil.direction_string_to_offset(d)).norm() if not self.flux_scaling else 1)
A0 = sum([sp.Matrix(ps.stencil.direction_string_to_offset(d)).norm() if not self.flux_scaling:
for d in flux_field.staggered_stencil]) / self.dim A0 = sum([sp.Matrix(ps.stencil.direction_string_to_offset(d)).norm()
divergence /= A0 for d in flux_field.staggered_stencil]) / self.dim
divergence /= A0
source = self.discrete_source() source = self.discrete_source()
source = {s.lhs: s.rhs for s in 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): ...@@ -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) c = ps.stencil.direction_string_to_offset(neighbor)
v1 = v.neighbor_vector(c) 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 # going out
cond = sp.And(*[sp.Or(c[i] * v0[i] > 0, c[i] == 0) for i in range(len(v0))]) 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))] 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))]) 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 # coming in
cond = sp.And(*[sp.Or(c[i] * v1[i] < 0, c[i] == 0) for i in range(len(v1))]) 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))] overlap1 = [1 - sp.Abs(v1[i]) for i in range(len(v1))]
overlap2 = [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))]) 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): for i, ff in enumerate(fluxes):
fluxes[i] = ff[0] fluxes[i] = ff[0]
......
...@@ -12,9 +12,11 @@ def test_advection_diffusion(dim: int): ...@@ -12,9 +12,11 @@ def test_advection_diffusion(dim: int):
domain_size = (32, 32) domain_size = (32, 32)
flux_neighbors = 4 flux_neighbors = 4
elif dim == 3: elif dim == 3:
domain_size = (16, 16, 16) domain_size = (32, 32, 32)
flux_neighbors = 13 flux_neighbors = 13
domain_center = [dom // 2 for dom in domain_size]
dh = ps.create_data_handling( dh = ps.create_data_handling(
domain_size=domain_size, periodicity=True, default_target='cpu') domain_size=domain_size, periodicity=True, default_target='cpu')
...@@ -24,13 +26,13 @@ def test_advection_diffusion(dim: int): ...@@ -24,13 +26,13 @@ def test_advection_diffusion(dim: int):
velocity_field = dh.add_array('v', values_per_cell=dim) velocity_field = dh.add_array('v', values_per_cell=dim)
D = 0.0666 D = 0.0666
time = 200 time = 300
def grad(f): def grad(f):
return sp.Matrix([ps.fd.diff(f, i) for i in range(dim)]) return sp.Matrix([ps.fd.diff(f, i) for i in range(dim)])
flux_eq = - D * grad(n_field) 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) vof_adv = ps.fd.VOF(j_field, velocity_field, n_field)
...@@ -47,20 +49,21 @@ def test_advection_diffusion(dim: int): ...@@ -47,20 +49,21 @@ def test_advection_diffusion(dim: int):
sync_conc = dh.synchronization_function([n_field.name]) 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): def density(pos: np.ndarray, time: int):
return (4 * np.pi * D * time)**(-1.5) * \ return sum((4 * np.pi * D * time) ** (-dim / 2) * np.exp(
np.exp(-np.sum(np.square(pos), axis=dim) / (4 * D * time)) -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)) pos = np.zeros((*domain_size, dim))
xpos = np.arange(-domain_size[0] // 2, domain_size[0] // 2) xpos = np.arange(-domain_size[0] // 2, domain_size[0] // 2)
ypos = np.arange(-domain_size[1] // 2, domain_size[1] // 2) ypos = np.arange(-domain_size[1] // 2, domain_size[1] // 2)
if dim == 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: elif dim == 3:
zpos = np.arange(-domain_size[2] // 2, domain_size[2] // 2) 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): def run(velocity: np.ndarray, time: int):
print(f"{velocity}, {time}") print(f"{velocity}, {time}")
...@@ -71,8 +74,7 @@ def test_advection_diffusion(dim: int): ...@@ -71,8 +74,7 @@ def test_advection_diffusion(dim: int):
for i in range(dim): for i in range(dim):
dh.fill(velocity_field.name, velocity[i], i, ghost_layers=True, inner_ghost_layers=True) 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, 0)
dh.fill(n_field.name, 1, slice_obj=ps.make_slice[[ dh.fill(n_field.name, 1, slice_obj=ps.make_slice[domain_center])
dom // 2 for dom in domain_size]])
sync_conc() sync_conc()
for i in range(time): for i in range(time):
...@@ -80,17 +82,23 @@ def test_advection_diffusion(dim: int): ...@@ -80,17 +82,23 @@ def test_advection_diffusion(dim: int):
dh.run_kernel(pde_kernel) dh.run_kernel(pde_kernel)
sync_conc() 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( assert dh.min(n_field.name) >= 0
n_field.name), calc_density, atol=1e-2, rtol=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.08, 0.08], repeat=dim): for vel in product([0, -0.02, 0.02], repeat=dim):
run(np.array(vel), time) run(np.array(vel), time)
def test_ek(): def test_ek():
# parameters # parameters
L = (40, 40) L = (40, 40)
...@@ -112,7 +120,7 @@ def test_ek(): ...@@ -112,7 +120,7 @@ def test_ek():
flux_eq = -D * Gradient(c) + D * z * c.center * Gradient(Phi) 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) flux_assignments = disc.discrete_flux(j)
advection_assignments = ps.fd.VOF(j, v, c) advection_assignments = ps.fd.VOF(j, v, c)
continuity_assignments = disc.discrete_continuity(j) continuity_assignments = disc.discrete_continuity(j)
......
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