One moment method

lbmweights/mrt_moments.py

@@ -48,96 +48,88 @@ def get_weights(selected_moments, c_s_sq, velocities, order):
     return method.weights
 
 
-def get_v17_moments():
+def get_moments(lattice):
     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)
+    if lattice.q == 17:
+        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()
+    elif lattice.q == 37:
+        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)
+
+        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()
+    else:
+        raise ValueError("Correct moments only found for V17 and 37")
 
     return selected_moments
 
 
-#moments = get_v17_moments()
-moments = get_v37_moments()
+lattice = Lattice.from_name("D2V17")
+moments = get_moments(lattice)
 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)
-
+M = sp.Matrix(get_moments_of_continuous_maxwellian_equilibrium(moments=moments, dim=2, order=4, c_s_sq=lattice.c_s_sq))
+m = make_matrix(moments, lattice.velocities)
+pprint(m.inv() * M)