Moments for lbmweights

lbmweights/mrt_moments.py

@@ -133,3 +133,89 @@ M = sp.Matrix(get_moments_of_continuous_maxwellian_equilibrium(moments=moments,
 m = make_matrix(moments, lattice.velocities)
 pprint(m.inv() * M)
 
+tau = 0.55
+lb_method = create_with_continuous_maxwellian_eq_moments(stencil=lattice.velocities,
+                                                         moment_to_relaxation_rate_dict={m: 1/tau for m in moments},
+                                                         compressible=True,
+                                                         equilibrium_order=None,
+                                                         c_s_sq=lattice.c_s_sq,
+                                                         )
+
+def l2norm(approximation, solution):
+    """
+    Calculating the l2norm of two matrices.
+    The matrices have to have the same shape
+
+    :return Value that indicates how close mat2 i
+    """
+    assert (
+        approximation.shape == solution.shape
+    ), "The shape of the matrices have to be the same"
+
+    diff = ((approximation - solution) ** 2).sum()
+    abs = (solution ** 2).sum()
+    norm = diff / abs
+    return np.sqrt(norm)
+
+
+def taylor_green(data_shape, t, tau, c_s_sq, test=False):
+    N_x = data_shape[0]
+    N_y = data_shape[1]
+    u_0 = 0.03
+    c_s_sq = float(c_s_sq)
+    p_0 = c_s_sq
+    roh_avg = 0.1
+
+    # Reynolds configs
+    viscocity = c_s_sq * (tau - 0.5)
+    c_s_sq = float(c_s_sq)
+
+    kx = 2 * np.pi
+    ky = 2 * np.pi
+
+    assert abs(c_s_sq * (tau - 0.5) - viscocity) <= 1e-10
+
+    X, Y = np.meshgrid(
+        np.linspace(0, 1, num=N_x, endpoint=True),
+        np.linspace(0, 1, num=N_y, endpoint=True),
+        indexing='ij'
+    )
+
+    t_d = float(1 / (viscocity * ((2 * np.pi/N_x) ** 2 + (2 * np.pi/N_y) ** 2)))
+    u_x = u_0 * (-1 * np.cos(kx * X) * np.sin(ky * Y)) * np.exp(- t / t_d)
+    u_y = u_0 * (np.sin(kx * X) * np.cos(ky * Y)) * np.exp(- t / t_d)
+    rho = (p_0 - roh_avg * u_0 ** 2 / 4 * (np.cos(2 * kx * X) + np.cos(2 * ky * Y)) * np.exp(-2*t/t_d))/c_s_sq
+    return rho, u_x, u_y
+
+def get_error(sc):
+    rho_analytical, ux_analytical, uy_analytical = taylor_green(sc.domain_size, sc.time_steps_run, tau, 1/3)
+    uy_analytical = - uy_analytical
+    simulated = sc.velocity[:, :, 0] ** 2 + sc.velocity[:, :, 1] ** 2
+    analytical = ux_analytical**2 + uy_analytical**2
+    return l2norm(simulated, analytical)
+
+
+print(lattice.velocities)
+sc = LatticeBoltzmannStep(stencil="D2V37", domain_size=(100, 100), periodicity=True, lb_method=lb_method)
+rho, u_x, u_y = taylor_green((100, 100), 0, tau, lattice.c_s_sq)
+init_u = np.stack((u_x, u_y), axis=2)
+np.copyto(sc.data_handling.cpu_arrays[sc.density_data_name][1:-1, 1:-1],
+          rho)
+np.copyto(sc.data_handling.cpu_arrays[sc.velocity_data_name][1:-1, 1:-1, 0],
+          u_x)
+np.copyto(sc.data_handling.cpu_arrays[sc.velocity_data_name][1:-1, 1:-1, 1],
+          u_y)
+sc.set_pdf_fields_from_macroscopic_values()
+
+sc.run_iterative_initialization(check_residuum_after=5)
+
+time = []
+errors = []
+for x in range(1000):
+    sc.run(1)
+    time.append(sc.time_steps_run)
+    errors.append(get_error(sc))
+
+plt.plot(time, errors)
+plt.show()
+