Commit efb97424 authored by Michael Kuron's avatar Michael Kuron
Browse files

Test the simple force model’s moments

parent 46415911
......@@ -82,7 +82,8 @@ def test_total_momentum(method, force_model, omega):
@pytest.mark.parametrize("stencil", ["D2Q9", "D3Q15", "D3Q19", "D3Q27"])
def test_stress(stencil):
@pytest.mark.parametrize("force_model", ["simple", "schiller"])
def test_modes(stencil, force_model):
"""check Schiller's force term in mode space"""
stencil = get_stencil(stencil)
dim = len(stencil[0])
......@@ -98,31 +99,42 @@ def test_stress(stencil):
stencil=stencil,
relaxation_rates=[omega_s, omega_b, omega_o, omega_e, omega_o, omega_e],
compressible=True,
force_model="schiller",
force_model=force_model,
force=F)
force_moments = sp.simplify(method.moment_matrix * sp.Matrix(method.force_model(method)))
# The momentum modes should contain the force
assert list(force_moments[1:dim+1]) == F
# The mass mode should be zero
assert force_moments[0] == 0
# The stress modes 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 + omega_s) * (traceless(sp.Matrix(u) * sp.Matrix(F).transpose()) + \
traceless(sp.Matrix(F) * sp.Matrix(u).transpose())) + \
sp.Rational(1,3) * (2 + omega_b) * sp.Matrix(u).dot(F) * sp.eye(method.dim)
# The momentum moments should contain the force
assert list(force_moments[1:dim+1]) == F
num_stresses = (dim*dim-dim)//2+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)}
for force_moment, moment in zip(force_moments[dim+1:dim+1+num_stresses],
method.moments[dim+1:dim+1+num_stresses]):
ref = moment.subs(subs)
diff = sp.simplify(ref - force_moment)
if is_bulk_moment(moment, dim):
assert diff == 0 or isinstance(diff, sp.Rational) # difference should be zero or a constant
else:
assert diff == 0 # difference should be zero
if force_model == "schiller":
num_stresses = (dim*dim-dim)//2+dim
# 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 + omega_s) * (traceless(sp.Matrix(u) * sp.Matrix(F).transpose()) + \
traceless(sp.Matrix(F) * sp.Matrix(u).transpose())) + \
sp.Rational(1,3) * (2 + omega_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)}
for force_moment, moment in zip(force_moments[dim+1:dim+1+num_stresses],
method.moments[dim+1:dim+1+num_stresses]):
ref = moment.subs(subs)
diff = sp.simplify(ref - force_moment)
if is_bulk_moment(moment, dim):
assert diff == 0 or isinstance(diff, sp.Rational) # difference should be zero or a constant
else:
assert diff == 0 # difference should be zero
# All other moments should be zero
assert list(force_moments[dim+1+num_stresses:]) == [0] * (len(stencil)-(dim+1+num_stresses))
elif force_model == "simple":
# All other moments should be zero
assert list(force_moments[dim+1:]) == [0] * (len(stencil)-(dim+1))
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