Commit 93a1473a authored by mischa's avatar mischa
Browse files

Mrt Moments

parent efb48fda
from lbmpy.session import *
from lbmpy.moments import *
from collections import OrderedDict
from lbmpy.methods.creationfunctions import create_with_continuous_maxwellian_eq_moments
from lbmpy.maxwellian_equilibrium import *
from lbmweights import Lattice
from pprint import pprint
import numpy as np
from numpy.linalg import svd
x, y, z = MOMENT_SYMBOLS
def get_duplicate_moment(selected_moments, velocities):
duplicate = []
m = moment_matrix(selected_moments, velocities)
rows =[m.row(i) for i in range(m.rows)]
for i, r1 in enumerate(rows):
for j, r2 in enumerate(rows):
if (r1 == r2) and i != j and i > j:
print("{}, {}".format(i, j))
try:
duplicate.append(j)
#selected_moments.remove(selected_moments_repr[j])
#selected_moments.remove(selected_moments_repr[i])
except ValueError:
pass
return duplicate
def make_matrix(selected_moments, stencil, p=True):
m = moment_matrix(selected_moments, stencil)
if p:
print("Shape: {}".format(m.shape))
print("Rank: {}".format(m.rank()))
return m
def moments_as_polynomials(selected_moments):
return list(exponents_to_polynomial_representations(selected_moments))
def get_weights(selected_moments, c_s_sq, velocities, order):
mom_to_rr = OrderedDict((k, sp.Symbol("omega")) for k in selected_moments)
method = create_with_continuous_maxwellian_eq_moments(velocities, mom_to_rr, compressible=True,
equilibrium_order=order,
c_s_sq=c_s_sq)
return method.weights
def get_v17_moments():
lattice = Lattice.from_name("D2V17")
d2v17 = lattice.velocities
selected_moments_repr = moments_up_to_order(6, dim=2)
selected_moments_repr = list(exponents_to_polynomial_representations(selected_moments_repr))
print(len(selected_moments_repr))
selected_moments = moments_up_to_order(6, dim=2)
selected_moments = list(exponents_to_polynomial_representations(selected_moments))
duplicate = get_duplicate_moment(selected_moments, lattice.velocities)
duplicate
#m = moment_matrix(selected_moments, d2v17)
selected_moments.remove(selected_moments_repr[13])
selected_moments.remove(selected_moments_repr[16])
selected_moments.remove(selected_moments_repr[17])
selected_moments.remove(selected_moments_repr[22])
selected_moments.remove(selected_moments_repr[23])
selected_moments.remove(selected_moments_repr[26])
m = make_matrix(selected_moments, d2v17)
#moments_as_polynomials(selected_moments)
selected_moments
selected_moments.remove(x*y**3)
selected_moments.remove(x**3*y**2)
selected_moments.remove(x**4*y)
selected_moments.remove(x**4*y**2)
selected_moments.append(x**6 + y**6)
selected_moments.remove(y**6)
selected_moments.remove(x**6)
m = make_matrix(selected_moments, d2v17)
selected_moments
assert (np.array(get_weights(selected_moments, lattice.c_s_sq, d2v17, 6)) == lattice.weights).all()
return selected_moments
def get_v37_moments():
lattice = Lattice.from_name("D2V37")
d2v37 = lattice.velocities
selected_moments_repr = moments_up_to_order(8, dim=2)
selected_moments_repr = list(exponents_to_polynomial_representations(selected_moments_repr))
selected_moments = moments_up_to_order(8, dim=2)
selected_moments = list(exponents_to_polynomial_representations(selected_moments))
duplicate = get_duplicate_moment(selected_moments, lattice.velocities)
m = make_matrix(selected_moments=selected_moments, stencil=lattice.velocities, p=False)
print("Amount of moment: " + str(len(selected_moments)))
pprint(m)
selected_moments.remove(x**8)
selected_moments.remove(y**8)
selected_moments.remove(x**7*y)
selected_moments.remove(x**1*y**7)
selected_moments.remove(x**5*y**3)
selected_moments.remove(x**3*y**5)
selected_moments.remove(x**7)
selected_moments.remove(y**7)
#selected_moments.remove(x**6*y**2)
#selected_moments.remove(x**2*y**6)
#selected_moments.append(x**2*y**6 + x**6*y**2)
pprint(selected_moments)
m = make_matrix(selected_moments, lattice.velocities)
print(m.shape)
calc_weights = get_weights(selected_moments, lattice.c_s_sq, d2v37, 8)
actual_weights = [float(w) for w in lattice.weights]
print(calc_weights)
print(actual_weights)
assert (np.array(calc_weights) - actual_weights <= 1e-15).all()
selected_moments = gram_schmidt(selected_moments, lattice.velocities, lattice.weights)
M = sp.Matrix(get_moments_of_continuous_maxwellian_equilibrium(moments=selected_moments, dim=2, order=4, c_s_sq=lattice.c_s_sq))
pprint(m.inv() * M)
return selected_moments
#moments = get_v17_moments()
moments = get_v37_moments()
pprint(moments)
#M = sp.Matrix(get_moments_of_continuous_maxwellian_equilibrium(moments=new_moments, dim=2, order=3, c_s_sq=lattice.c_s_sq))
#pprint(m.inv() * M)
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