# -*- coding: utf-8 -*- """ Created on Mon Feb 1 11:42:59 2021 @author: Daniel Powell """ import numpy as np from numba import njit, prange from matplotlib import pyplot as plt import matplotlib import matplotlib.cm as cm from scipy.stats import skew, kurtosis import fvm_2d #matplotlib.rc('text', usetex=True) #matplotlib.rcParams['text.latex.preamble'] = r"\usepackage{amsmath}" @njit(nogil=True) def get_RK_ks(f, tn, vn, dt, params): k1 = f(tn, vn, params) k2 = f(tn + ((1 / 5) * dt), vn + ((1 / 5) * dt * k1), params) k3 = f(tn + ((3 / 10) * dt), vn + ((dt / 40) * ((3 * k1) + (9 * k2))), params) k4 = f(tn + ((3 / 5) * dt), vn + (dt * (((3 / 10) * k1) - ((9 / 10) * k2) + ((6 / 5) * k3))), params) k5 = f(tn + dt, vn + (dt * (((-11 / 54) * k1) + ((5 / 2) * k2) - ((70 / 27) * k3) + ((35 / 27) * k4))), params) k6 = f(tn + ((7 / 8) * dt), vn + (dt * (((1631 / 55296) * k1) + ((175 / 512) * k2) + ((575 / 13824) * k3) + ((44275 / 110592) * k4) + ((253 / 4096) * k5))), params) return (k1, k2, k3, k4, k5, k6) @njit(nogil=True) def RK5_step(f, tn, vn, dt, params): (k1, k2, k3, k4, k5, k6) = get_RK_ks(f, tn, vn, dt, params) return vn + (dt * (((37 / 378) * k1) + ((250 / 621) * k3) + ((125 / 594) * k4) + ((512 / 1771) * k6))) @njit(nogil=True) def fluid_vel(X, params): return params[0] * np.array([np.cos(params[1] * X[0]) * np.sin(params[1] * X[1]), -np.sin(params[1] * X[0]) * np.cos(params[1] * X[1])]) @njit(nogil=True) def drag(t, Z, params): u = fluid_vel(Z[:2], params) a = params[2] * (u - Z[2:]) return np.array([Z[2], Z[3], a[0], a[1]]) @njit(nogil=True) def boundary_condition(Z): if Z[0] > l: Z[0] -= l if Z[0] <= -l: Z[0] += l if Z[1] > l: Z[1] -= l if Z[1] <= -l: Z[1] += l return Z @njit(nogil=True, parallel=True) def boundary_condition_all(Z): for i in prange(len(Z)): Z[i, :] = boundary_condition(Z[i, :]) return Z @njit(nogil=True) def single_loop(Z_array, dt, params): t = 0 for i in range(np.shape(Z_array)[0] - 1): Z_array[i + 1, :] = RK5_step(drag, t, Z_array[i, :], dt, params) t += dt Z_array[i + 1, :] = boundary_condition(Z_array[i + 1, :]) return Z_array @njit(nogil=True, parallel=True) def many_loop(Z_array, dt, params): for i in prange(np.shape(Z_array)[0]): Z_array[i, :, :] = single_loop(Z_array[i, :, :], dt, params) return Z_array # taylor_green_scales2021.py a = 35221 # /m umag = 0.29283 # m/s (eddy speed) l = np.pi / a # m (eddy diameter) L = 2 * l # m (domain size) tau_f = l / umag # cool trajectory def single_orbit(Stk, N, filename='1em3'): tau_p = Stk * tau_f inv_tau_p = 1.0 / tau_p # completes a full loop after 65568 iterations because # theta = np.arctan2(Z[:, 1], Z[:, 0]) # np.where(np.abs(theta - theta[0]) < 1e-5) params = np.array([umag, a, inv_tau_p]) dt = tau_p / 10 Z = np.zeros((N, 4)) Z[0, :] = np.array([l * 0.45, l * 0.45, 0, 0]) Z = single_loop(Z, dt, params) (X, Y, XF, YF, dx, Area, V) = fvm_2d.create_uniform_grid(256, l) u, v, w, neighbours = fvm_2d.velocities_and_neighbours(X, Y, 256, umag, a, tau_p) # examine particle velocity error between Lagrangian and EE along trajectory def ee_tg(x, y, tau_p, umag, a): v_x = umag * ((np.cos(a * x) * np.sin(a * y)) + (0.5 * tau_p * umag * a * np.sin(2 * a * x))) v_y = umag * ((-np.sin(a * x) * np.cos(a * y)) + (0.5 * tau_p * umag * a * np.sin(2 * a * y))) return np.array([v_x, v_y]) v_ee = ee_tg(Z[:, 0], Z[:, 1], tau_p, umag, a) v_ee = v_ee.reshape(N, 2) v_error = np.linalg.norm(Z[:, :2] - v_ee, axis=1) angle_error = np.zeros_like(v_error) for i in range(N): dot = np.dot(Z[i, 2:], v_ee[i]) abs_a = np.sqrt((Z[i, 2] * Z[i, 2]) + (Z[i, 3] * Z[i, 3])) abs_b = np.sqrt((v_ee[i, 0] * v_ee[i, 0]) + (v_ee[i, 1] * v_ee[i, 1])) angle_error[i] = np.arccos(dot / (abs_a * abs_b)) # maxima diff = np.diff(v_error) increasing = np.where(diff > 0)[0] array_increase = np.where(np.diff(increasing) > 1)[0] maxima = increasing[array_increase] maxima = maxima[v_error[maxima] > 0.3] t = np.linspace(0, dt * N, N) plt.figure() plt.plot(t / tau_p, v_error) plt.xlabel(r"$\frac{t}{\tau_p}$", fontsize=16) plt.ylabel("Normalised speed error", fontsize=16) plt.ylim((0, 0.45)) #plt.scatter(maxima, v_error_stk3[maxima]) plt.tight_layout() plt.savefig("tg_single_orbit_error_speed_stk" + filename + ".eps") fig, ax = plt.subplots() plt.contourf(X / (2 * l), Y / (2 * l), np.linalg.norm(u / umag, axis=2)) plt.plot(Z[:, 0] / (2 * l), Z[:, 1] / (2 * l), 'k') plt.scatter(Z[maxima, 0] / (2 * l), Z[maxima, 1] / (2 * l), color='r') plt.xlabel(r"$\frac{x}{2l}$", fontsize=16) plt.ylabel(r"$\frac{y}{2l}$", fontsize=16) plt.axis('square') plt.tight_layout() plt.savefig("tg_single_orbit_error_stk" + filename + ".eps") plt.figure() plt.plot(t / tau_p, angle_error) plt.xlabel(r"$\frac{t}{\tau_p}$", fontsize=16) plt.ylabel("Angle error", fontsize=16) plt.ylim((0, 3.5)) #plt.scatter(maxima, v_error_stk3[maxima]) plt.tight_layout() plt.savefig("tg_single_orbit_error_angle_stk" + filename + ".eps") return Z, t, v_error # for Stk = 1e-3, 1e-2 # N = 65568, 7205 Z, t, v_error = single_orbit(1e-3, 65568, filename='1em3') Z, t, v_error = single_orbit(1e-2, 7205, filename='1em2') Z, t, v_error = single_orbit(1e-1, 2282, filename='1em1') # ============================================================================= # Stk = 1e-2 # tau_p = Stk * tau_f # inv_tau_p = 1.0 / tau_p # # params = np.array([umag, a, inv_tau_p]) # ============================================================================= # ============================================================================= # rs = np.random.RandomState(np.random.MT19937(np.random.SeedSequence(1994))) # # state vector Z = (x, y, u, v) = (position, velocity) # NT = 10000 # NP = 1000 # dt = tau_p / 100 # Z_array = np.zeros((NP, NT, 4)) # Z_array[:, 0, :2] = L * (np.random.random((NP, 2)) - 0.5) # ============================================================================= # error study # ============================================================================= # dt_over_tau_p = np.array([1, 0.5, 0.1, 0.01, 1e-3]) # # X_final = np.array([[3.982866730475551e-05, 3.3439593457412727e-06], # [3.982833215030401e-05, 3.3613026108508203e-06], # [3.9828499940903846e-05, 3.352631146450209e-06], # [3.982866730475904e-05, 3.343959345900221e-06], # [3.982866730476029e-05, 3.343959345906732e-06]]) # # error = np.linalg.norm(X_final[:-1] - X_final[-1], axis=1) / l # # # now for stk = 1e-2 # Stk = 1e-2 # tau_p = Stk * tau_f # inv_tau_p = 1.0 / tau_p # # params = np.array([umag, a, inv_tau_p]) # f = 0.001 # dt = tau_p * f # # careful of floating point precision stuff # NT = int(100 / f) # # ============================================================================= # # Z_array = np.zeros((NT, 4)) # # Z_array[0, :2] = l * np.array([0., 0.45]) # # # # Z_array[0, 2:] = fluid_vel(Z_array[0, 0:2], params) # # # # t = np.arange(0, dt * NT, dt) # # # # Z_array = single_loop(Z_array, dt, params) # # # # plt.figure() # # plt.plot(Z_array[:, 0], Z_array[:, 1]) # # plt.axis('square') # # ============================================================================= # # X_final2 = np.array([[4.028229726663308e-05, 5.876446454043164e-06], # [4.0296529621091614e-05, 5.450448043900599e-06], # [4.03072219530019e-05, 5.108570550497141e-06], # [4.0309542987303906e-05, 5.031521393499652e-06], # [4.030977338079171e-05, 5.023813973204555e-06]]) # # error2 = np.linalg.norm(X_final2[:-1] - X_final2[-1], axis=1) / l # # # ============================================================================= # # plt.figure() # # #plt.scatter(dt_over_tau_p[:-1], error) # # plt.scatter(dt_over_tau_p[:-1], error2) # # plt.xscale('log') # # plt.yscale('log') # # plt.xlabel(r"$\frac{\Delta t_p}{\tau_p}$", fontsize=14) # # plt.ylabel("Normalised Error") # # plt.tight_layout() # # plt.savefig("rk5_cash_timestep.eps") # # ============================================================================= # # # # # # # # # # error study Stk = 0.1 # # dt_over_tau_p = np.array([1, 0.5, 0.1, 0.01, 1e-3]) # # # now for stk = 1e-2 # Stk = 0.1 # tau_p = Stk * tau_f # inv_tau_p = 1.0 / tau_p # # params = np.array([umag, a, inv_tau_p]) # f = 0.001 # dt = tau_p * f # # careful of floating point precision stuff # NT = int(10 / f) # # ============================================================================= # # Z_array = np.zeros((NT, 4)) # # Z_array[0, :2] = l * np.array([0., 0.45]) # # # # Z_array[0, 2:] = fluid_vel(Z_array[0, 0:2], params) # # # # t = np.arange(0, dt * NT, dt) # # # # Z_array = single_loop(Z_array, dt, params) # # # # plt.figure() # # plt.plot(Z_array[:, 0], Z_array[:, 1]) # # plt.axis('square') # # ============================================================================= # # X_final_stk01 = np.array([[4.317595654983959e-05, 2.3031536200832674e-05], # [4.343028807436551e-05, 2.0399795025990722e-05], # [4.358408401176122e-05, 1.8099394996742306e-05], # [4.361360423872948e-05, 1.7557505028696605e-05], # [4.361646265817635e-05, 1.7502823046153542e-05]]) # # error_stk01 = np.linalg.norm(X_final_stk01[:-1] - X_final_stk01[-1], axis=1) / l # ============================================================================= # ============================================================================= # plt.figure() # plt.scatter(dt_over_tau_p[:-1], error_stk01) # plt.xscale('log') # plt.yscale('log') # plt.xlabel(r"$\frac{\Delta t_p}{\tau_p}$", fontsize=14) # plt.ylabel("Normalised Error") # plt.tight_layout() # plt.savefig("rk5_cash_timestep_stk01.eps") # ============================================================================= # ============================================================================= # plt.figure() # plt.scatter(dt_over_tau_p[:-1], error2, label='$Stk = 0.01$') # plt.scatter(dt_over_tau_p[:-1], error_stk01, label='$Stk = 0.1$') # plt.xscale('log') # plt.yscale('log') # plt.xlabel(r"$\frac{\Delta t_p}{\tau_p}$", fontsize=14) # plt.ylabel("Normalised Error") # plt.legend() # plt.tight_layout() # #plt.savefig("rk5_cash_timestep_both.eps") # ============================================================================= # ============================================================================= # # Particle Number + Time Study (Stk = 1e-2) # Stk = 1e-2 # tau_p = Stk * tau_f # inv_tau_p = 1.0 / tau_p # # params = np.array([umag, a, inv_tau_p]) # # rs = np.random.RandomState(np.random.MT19937(np.random.SeedSequence(1994))) # # state vector Z = (x, y, u, v) = (position, velocity) # NT = 300000 # NP = 10 # dt = tau_p / 10 # Z_array = np.zeros((NP, NT, 4)) # Z_array[:, 0, :2] = L * (np.random.random((NP, 2)) - 0.5) # # Z_array[:, 0, 2:] = 2 * umag * (np.random.random((NP, 2)) - 0.5) # # # ============================================================================= # # Z_array = many_loop(Z_array, dt, params) # # # # r = np.sqrt(Z_array[:, :, 0] ** 2 + Z_array[:, :, 1] ** 2) # # # # mean = r.mean(axis=0) # # std = r.std(axis=0) # # skewness = skew(r, axis=0) # # kurt = kurtosis(r, axis=0) # # ============================================================================= # # t_norm = np.arange(0, dt * NT, dt) / tau_p # # mean10 = np.load("D:\\tg_particles\\mean10_vrand.npy") # mean100 = np.load("D:\\tg_particles\\mean100_vrand.npy") # mean1000_1 = np.load("D:\\tg_particles\\mean1000_vrand_part1.npy") # mean1000_2 = np.load("D:\\tg_particles\\mean1000_vrand_part2.npy") # mean1000_3 = np.load("D:\\tg_particles\\mean1000_vrand_part3.npy") # mean1000_4 = np.load("D:\\tg_particles\\mean1000_vrand_part4.npy") # mean1000_5 = np.load("D:\\tg_particles\\mean1000_vrand_part5.npy") # mean1000_6 = np.load("D:\\tg_particles\\mean1000_vrand_part6.npy") # mean1000 = np.concatenate((mean1000_1, mean1000_2, mean1000_3, mean1000_4, # mean1000_5, mean1000_6)) # # std10 = np.load("D:\\tg_particles\\std10_vrand.npy") # std100 = np.load("D:\\tg_particles\\std100_vrand.npy") # std1000_1 = np.load("D:\\tg_particles\\std1000_vrand_part1.npy") # std1000_2 = np.load("D:\\tg_particles\\std1000_vrand_part2.npy") # std1000_3 = np.load("D:\\tg_particles\\std1000_vrand_part3.npy") # std1000_4 = np.load("D:\\tg_particles\\std1000_vrand_part4.npy") # std1000_5 = np.load("D:\\tg_particles\\std1000_vrand_part5.npy") # std1000_6 = np.load("D:\\tg_particles\\std1000_vrand_part6.npy") # std1000 = np.concatenate((std1000_1, std1000_2, std1000_3, std1000_4, # std1000_5, std1000_6)) # # skewness10 = np.load("D:\\tg_particles\\skewness10_vrand.npy") # skewness100 = np.load("D:\\tg_particles\\skewness100_vrand.npy") # skewness1000_1 = np.load("D:\\tg_particles\\skewness1000_vrand_part1.npy") # skewness1000_2 = np.load("D:\\tg_particles\\skewness1000_vrand_part2.npy") # skewness1000_3 = np.load("D:\\tg_particles\\skewness1000_vrand_part3.npy") # skewness1000_4 = np.load("D:\\tg_particles\\skewness1000_vrand_part4.npy") # skewness1000_5 = np.load("D:\\tg_particles\\skewness1000_vrand_part5.npy") # skewness1000_6 = np.load("D:\\tg_particles\\skewness1000_vrand_part6.npy") # skewness1000 = np.concatenate((skewness1000_1, skewness1000_2, skewness1000_3, skewness1000_4, # skewness1000_5, skewness1000_6)) # # kurt10 = np.load("D:\\tg_particles\\kurt10_vrand.npy") # kurt100 = np.load("D:\\tg_particles\\kurt100_vrand.npy") # kurt1000_1 = np.load("D:\\tg_particles\\kurt1000_vrand_part1.npy") # kurt1000_2 = np.load("D:\\tg_particles\\kurt1000_vrand_part2.npy") # kurt1000_3 = np.load("D:\\tg_particles\\kurt1000_vrand_part3.npy") # kurt1000_4 = np.load("D:\\tg_particles\\kurt1000_vrand_part4.npy") # kurt1000_5 = np.load("D:\\tg_particles\\kurt1000_vrand_part5.npy") # kurt1000_6 = np.load("D:\\tg_particles\\kurt1000_vrand_part6.npy") # kurt1000 = np.concatenate((kurt1000_1, kurt1000_2, kurt1000_3, kurt1000_4, # kurt1000_5, kurt1000_6)) # # # ============================================================================= # # np.save("D:\\tg_particles\\mean1000_vrand_part2", mean) # # np.save("D:\\tg_particles\\std1000_vrand_part2", std) # # np.save("D:\\tg_particles\\skewness1000_vrand_part2", skewness) # # np.save("D:\\tg_particles\\kurt1000_vrand_part2", kurt) # # ============================================================================= # # # ============================================================================= # # fig, ax = plt.subplots(sharex=True) # # plt.subplot(221) # # plt.plot(t_norm, mean10 / l, label='10') # # plt.plot(t_norm, mean100 / l, label='100') # # plt.plot(t_norm, mean1000 / l, label='1000') # # #plt.legend(title=r"$N_p$") # # plt.title("Mean") # # plt.xlabel(r"$\frac{t}{\tau_p}$", fontsize=16) # # plt.ylabel(r"$\frac{}{l}$", fontsize=16) # # plt.subplot(222) # # plt.plot(t_norm, std10 / l, label='10') # # plt.plot(t_norm, std100 / l, label='100') # # plt.plot(t_norm, std1000 / l, label='1000') # # #plt.legend(title=r"$N_p$") # # plt.title("Standard Deviation") # # plt.xlabel(r"$\frac{t}{\tau_p}$", fontsize=16) # # plt.ylabel(r"$\frac{M_2^{\frac{1}{2}}}{l}$", fontsize=16) # # plt.subplot(223) # # plt.plot(t_norm, skewness10, label='10') # # plt.plot(t_norm, skewness100, label='100') # # plt.plot(t_norm, skewness1000, label='1000') # # #plt.legend(title=r"$N_p$") # # plt.title("Skewness") # # plt.xlabel(r"$\frac{t}{\tau_p}$", fontsize=16) # # plt.ylabel(r"$\frac{M_3}{M_2^{\frac{3}{2}}}$", fontsize=16) # # plt.subplot(224) # # plt.plot(t_norm, kurt10, label='10') # # plt.plot(t_norm, kurt100, label='100') # # plt.plot(t_norm, kurt1000, label='1000') # # plt.legend(title=r"$N_p$") # # plt.title("Kurtosis") # # plt.xlabel(r"$\frac{t}{\tau_p}$", fontsize=16) # # plt.ylabel(r"$\frac{M_4}{M_2^4}$", fontsize=16) # # fig.tight_layout() # # #plt.savefig("tg_particle_moments_random.eps") # # ============================================================================= # # # normalised by l # # ============================================================================= # # mean_random = np.array([0.7071067811865314, 0.7071067811865319, # # 0.7071067811865323]) # # # normalised by l # # std_random = np.array([1.1869440184542384e-16, ]) # # ============================================================================= # ============================================================================= # Particle Number + Time Study (Stk = 1e-1) # ============================================================================= # Stk = 1e-1 # tau_p = Stk * tau_f # inv_tau_p = 1.0 / tau_p # # params = np.array([umag, a, inv_tau_p]) # # rs = np.random.RandomState(np.random.MT19937(np.random.SeedSequence(1994))) # # state vector Z = (x, y, u, v) = (position, velocity) # NT = 5000 # less time steps needed for same total time as Stk=1e-2 # NP = 10 # dt = tau_p / 10 # # ============================================================================= # # Z_array = np.zeros((NP, NT, 4)) # # Z_array[:, 0, :2] = L * (np.random.random((NP, 2)) - 0.5) # # # # Z_array[:, 0, 2:] = 2 * umag * (np.random.random((NP, 2)) - 0.5) # # # # Z_array = many_loop(Z_array, dt, params) # # # # r = np.sqrt(Z_array[:, :, 0] ** 2 + Z_array[:, :, 1] ** 2) # # # # mean = r.mean(axis=0) # # std = r.std(axis=0) # # skewness = skew(r, axis=0) # # kurt = kurtosis(r, axis=0) # # ============================================================================= # # t_norm = np.arange(0, dt * NT, dt) / tau_p # # mean10 = np.load("D:\\tg_particles\\mean10_stk01_vrand.npy") # mean100 = np.load("D:\\tg_particles\\mean100_stk01_vrand.npy") # mean1000 = np.load("D:\\tg_particles\\mean1000_stk01_vrand.npy") # mean10000 = np.load("D:\\tg_particles\\mean10000_stk01_vrand.npy") # # std10 = np.load("D:\\tg_particles\\std10_stk01_vrand.npy") # std100 = np.load("D:\\tg_particles\\std100_stk01_vrand.npy") # std1000 = np.load("D:\\tg_particles\\std1000_stk01_vrand.npy") # std10000 = np.load("D:\\tg_particles\\std10000_stk01_vrand.npy") # # skewness10 = np.load("D:\\tg_particles\\skewness10_stk01_vrand.npy") # skewness100 = np.load("D:\\tg_particles\\skewness100_stk01_vrand.npy") # skewness1000 = np.load("D:\\tg_particles\\skewness1000_stk01_vrand.npy") # skewness10000 = np.load("D:\\tg_particles\\skewness10000_stk01_vrand.npy") # # kurt10 = np.load("D:\\tg_particles\\kurt10_stk01_vrand.npy") # kurt100 = np.load("D:\\tg_particles\\kurt100_stk01_vrand.npy") # kurt1000 = np.load("D:\\tg_particles\\kurt1000_stk01_vrand.npy") # kurt10000 = np.load("D:\\tg_particles\\kurt10000_stk01_vrand.npy") # # # # ============================================================================= # # np.save("D:\\tg_particles\\mean10000_stk01_vrand", mean) # # np.save("D:\\tg_particles\\std10000_stk01_vrand", std) # # np.save("D:\\tg_particles\\skewness10000_stk01_vrand", skewness) # # np.save("D:\\tg_particles\\kurt10000_stk01_vrand", kurt) # # ============================================================================= # # fig, ax = plt.subplots(sharex=True) # plt.subplot(221) # plt.plot(t_norm, mean10 / l, label='10') # plt.plot(t_norm, mean100 / l, label='100') # plt.plot(t_norm, mean1000 / l, label='1000') # plt.plot(t_norm, mean10000 / l, label='10000') # #plt.legend(title=r"$N_p$") # plt.title("Mean") # plt.xlabel(r"$\frac{t}{\tau_p}$", fontsize=16) # plt.ylabel(r"$\frac{}{l}$", fontsize=16) # plt.subplot(222) # plt.plot(t_norm, std10 / l, label='10') # plt.plot(t_norm, std100 / l, label='100') # plt.plot(t_norm, std1000 / l, label='1000') # plt.plot(t_norm, std10000 / l, label='10000') # #plt.legend(title=r"$N_p$") # plt.title("Standard Deviation") # plt.xlabel(r"$\frac{t}{\tau_p}$", fontsize=16) # plt.ylabel(r"$\frac{M_2^{\frac{1}{2}}}{l}$", fontsize=16) # plt.subplot(223) # plt.plot(t_norm, skewness10, label='10') # plt.plot(t_norm, skewness100, label='100') # plt.plot(t_norm, skewness1000, label='1000') # plt.plot(t_norm, skewness10000, label='10000') # #plt.legend(title=r"$N_p$") # plt.title("Skewness") # plt.xlabel(r"$\frac{t}{\tau_p}$", fontsize=16) # plt.ylabel(r"$\frac{M_3}{M_2^{\frac{3}{2}}}$", fontsize=16) # plt.subplot(224) # plt.plot(t_norm, kurt10, label='10') # plt.plot(t_norm, kurt100, label='100') # plt.plot(t_norm, kurt1000, label='1000') # plt.plot(t_norm, kurt10000, label='10000') # plt.legend(title=r"$N_p$") # plt.title("Kurtosis") # plt.xlabel(r"$\frac{t}{\tau_p}$", fontsize=16) # plt.ylabel(r"$\frac{M_4}{M_2^4}$", fontsize=16) # fig.tight_layout() # #plt.savefig("tg_particle_moments_random_stk01.eps") # =============================================================================