diff --git a/pystencils/fd/finitevolumes.py b/pystencils/fd/finitevolumes.py
index 588a4a45b2bab9b5aeb652d7c632902f667ddc8c..64f9db6900dbd884c7fef5b05b6d9d6e7243e657 100644
--- a/pystencils/fd/finitevolumes.py
+++ b/pystencils/fd/finitevolumes.py
@@ -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}
diff --git a/pystencils_tests/test_fvm.py b/pystencils_tests/test_fvm.py
index 09b3e0be3dfc919e6c73c021009b9e87b65471ec..b33049a52b6775bdbb20b49992561ebc0406d38b 100644
--- a/pystencils_tests/test_fvm.py
+++ b/pystencils_tests/test_fvm.py
@@ -12,7 +12,7 @@ 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]
@@ -26,14 +26,13 @@ def test_advection_diffusion(dim: int):
velocity_field = dh.add_array('v', values_per_cell=dim)
D = 0.0666
- time = 200
- time_advection = 50
+ time = 300
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)
@@ -44,27 +43,27 @@ def test_advection_diffusion(dim: int):
flux.append(ps.Assignment(adv.lhs, adv.rhs + div.rhs))
flux_kernel = ps.create_staggered_kernel(flux).compile()
- adv_kernel = ps.create_staggered_kernel(vof_adv).compile()
-
+
pde_kernel = ps.create_kernel(
fvm_eq.discrete_continuity(j_field)).compile()
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}")
@@ -82,50 +81,24 @@ def test_advection_diffusion(dim: int):
dh.run_kernel(flux_kernel)
dh.run_kernel(pde_kernel)
sync_conc()
-
- calc_density = density(pos - velocity * time, time)
-
- np.testing.assert_allclose(dh.gather_array(
- n_field.name), calc_density, atol=1e-2, rtol=0)
- assert dh.min(n_field.name) >= 0
- for vel in product([0, -0.08, 0.08], repeat=dim):
- run(np.array(vel), 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)
- def run_advection_only(velocity: np.ndarray, time: int):
- print(f"{velocity}, {time}")
- dh.fill(n_field.name, np.nan, ghost_layers=True, inner_ghost_layers=True)
- dh.fill(j_field.name, np.nan, ghost_layers=True, inner_ghost_layers=True)
-
- # set initial values for velocity and density
- 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[domain_center])
+ analytical_maximum = analytic_density[peak_position]
+ relative_maximum_error = np.abs(simulation_density[peak_position] - analytical_maximum) / analytical_maximum
- sync_conc()
- for i in range(time):
- dh.run_kernel(adv_kernel)
- dh.run_kernel(pde_kernel)
- sync_conc()
-
- max_direction = np.unravel_index(np.argmax(dh.gather_array(n_field.name)), domain_size) - np.asarray(domain_center)
- norm = max(abs(max_direction))
- max_direction = max_direction / (norm if norm > 0 else 1)
-
- velocity_direction = velocity
- norm = max(abs(velocity_direction))
- velocity_direction = velocity_direction / (norm if norm > 0 else 1)
-
- np.testing.assert_array_equal(max_direction, velocity_direction)
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.6, 0.6], repeat=dim):
- run_advection_only(np.array(vel), time_advection)
+ for vel in product([0, -0.02, 0.02], repeat=dim):
+ run(np.array(vel), time)
def test_ek():
-
# parameters
L = (40, 40)
@@ -147,7 +120,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)