# -*- coding: utf-8 -*- """ Created on Mon Dec 7 19:14:55 2020 @author: Daniel Powell """ import numpy as np import fvm_2d def test_taylor_green_velocity(): # zero umag returns zero tg = fvm_2d.taylor_green_velocity(1, 1, 0, 1) assert np.allclose(tg, np.array([0.0, 0.0])) # zero a returns zero tg = fvm_2d.taylor_green_velocity(1, 1, 1, 0) assert np.allclose(tg, np.array([0.0, 0.0])) # velocity at (0, 0) is zero tg = fvm_2d.taylor_green_velocity(0, 0, 1, 1) assert np.allclose(tg, np.array([0.0, 0.0])) # velocity at (0, pi/2) is (umag, 0) umag = 1.0 tg = fvm_2d.taylor_green_velocity(0, np.pi / 2, umag, 1) assert np.allclose(tg, np.array([umag, 0.0])) umag = 9.7 tg = fvm_2d.taylor_green_velocity(0, np.pi / 2, umag, 1) assert np.allclose(tg, np.array([umag, 0.0])) def test_taylor_green_ee_velocity(): # trivial umag = 0 a = 1.0 tau_p = 1.0 x = 0.0 y = 0.0 ans = np.array([0.0, 0.0]) u = fvm_2d.taylor_green_velocity(x, y, umag, a) v = fvm_2d.taylor_green_mee_velocity(x, y, umag, a, tau_p, u) assert np.allclose(ans, v) # trivial again umag = 1.0 tau_p = 1.0 ans = np.array([0.0, 0.0]) u = fvm_2d.taylor_green_velocity(x, y, umag, a) v = fvm_2d.taylor_green_mee_velocity(x, y, umag, a, tau_p, u) assert np.allclose(ans, v) # another trivial l = np.pi / a x = l / 3 # arbitrary y = -l / 5 # arbitrary tau_p = 0.0 # should give same result as fluid velocity u = fvm_2d.taylor_green_velocity(x, y, umag, a) v = fvm_2d.taylor_green_mee_velocity(x, y, umag, a, tau_p, u) assert np.allclose(u, v) # still trivial (edge of vortex) umag = 2.0 a = 1.0 tau_p = 0.0 x = 0.0 y = np.pi / 2 ans = np.array([2.0, 0.0]) u = fvm_2d.taylor_green_velocity(x, y, umag, a) v = fvm_2d.taylor_green_ee_velocity(x, y, umag, a, tau_p, u) assert np.allclose(ans, v) # now a non-trivial one (edge of vortex) umag = 2.0 a = 1.0 tau_p = 3.0 x = 0.0 y = np.pi / 2 ans = np.array([2.0, 0.0]) u = fvm_2d.taylor_green_velocity(x, y, umag, a) v = fvm_2d.taylor_green_ee_velocity(x, y, umag, a, tau_p, u) assert np.allclose(ans, v) def test_taylor_green_mee_velocity(): # trivial umag = 0 a = 1.0 tau_p = 1.0 x = 0.0 y = 0.0 ans = np.array([0.0, 0.0]) u = fvm_2d.taylor_green_velocity(x, y, umag, a) v = fvm_2d.taylor_green_mee_velocity(x, y, umag, a, tau_p, u) assert np.allclose(ans, v) # if tau_p = 0, should return the fluid velocity umag = 1.0 tau_p = 0.0 u = fvm_2d.taylor_green_velocity(x, y, umag, a) v = fvm_2d.taylor_green_mee_velocity(x, y, umag, a, tau_p, u) assert np.allclose(u, v) # at the origin with everything = 1 tau_p = 1.0 u = fvm_2d.taylor_green_velocity(x, y, umag, a) v = fvm_2d.taylor_green_mee_velocity(x, y, umag, a, tau_p, u) assert np.allclose(u, v) def taylor_green_mee_velocity_inverse(x, y, umag, a, tau_p, u): """ Computes and returns the 2D Taylor-Green fluid velocity at a point (x, y), for a given umag and a by numerically computing the matrix inverse. u_x = umag cos(ax) sin(ay) u_y = -umag sin(ax) cos(ay) Parameters ---------- x : float x-coordinate (m). y : float y-coordinate (m). umag : float maximum fluid velocity (m/s). a : float flow parameter (/m). tau_p: float particle relaxation time (s). u: array of float local fluid velocity vector (m/s). Returns ------- ee_vel : array of float x and y EE velocity components. """ A = np.array([[1 - (tau_p * umag * a * np.sin(a * x) * np.sin(a * y)), -tau_p * umag * a * np.cos(a * x) * np.cos(a * y)], [tau_p * umag * a * np.cos(a * x) * np.cos(a * y), 1 + (tau_p * umag * a * np.sin(a * x) * np.sin(a * y))]]) A_inv = np.linalg.inv(A) ee_vel = u rel = ((0.5 * tau_p * umag * umag * a) * np.array([np.sin(2 * a * x), np.sin(2 * a * y)])) w = np.matmul(A_inv, rel) ee_vel[0] += w[0] ee_vel[1] += w[1] return ee_vel umag = 0.29283 a = 35521 tau_p = 1e-8 l = np.pi / a x = l / 2 y = l / 2 u = fvm_2d.taylor_green_velocity(x, y, umag, a) v_ans = taylor_green_mee_velocity_inverse(x, y, umag, a, tau_p, u) v = fvm_2d.taylor_green_mee_velocity(x, y, umag, a, tau_p, u) assert np.allclose(v_ans, v) x = 0.2 * l y = -0.3 * l v_ans = taylor_green_mee_velocity_inverse(x, y, umag, a, tau_p, u) v = fvm_2d.taylor_green_mee_velocity(x, y, umag, a, tau_p, u) assert np.allclose(v_ans, v) x = -1.3 * l y = 0.75 * l v_ans = taylor_green_mee_velocity_inverse(x, y, umag, a, tau_p, u) v = fvm_2d.taylor_green_mee_velocity(x, y, umag, a, tau_p, u) assert np.allclose(v_ans, v) def test_get_all_velocities_analytical(): # all velocities zero at the origin u, v, w = fvm_2d.get_all_velocities_analytical(0.0, 0.0, 1.0, 1.0, 1.0) assert np.allclose(u, np.array([0.0, 0.0])) assert np.allclose(v, np.array([0.0, 0.0])) assert np.allclose(w, np.array([0.0, 0.0])) def test_create_uniform_grid(): # take a 2x2 grid on a domain (-1, 1) l = 1.0 N = 2 (X, Y, XF, YF, dx, Area, V) = fvm_2d.create_uniform_grid(N, l) X_ans = np.array([[-0.5, -0.5], [0.5, 0.5]]) Y_ans = X_ans.transpose() XF_ans = np.array([[-1.0, -1.0, -1.0], [0.0, 0.0, 0.0], [1.0, 1.0, 1.0]]) YF_ans = XF_ans.transpose() assert np.allclose(X, X_ans) assert np.allclose(Y, Y_ans) assert np.allclose(XF, XF_ans) assert np.allclose(YF, YF_ans) assert dx == 1.0 assert Area == 1.0 assert V == 1.0 def test_get_compass(): # say we have a 3x3 N = 3 # west, east, north, south lrud = np.array([0, 0, 0, 0], dtype='int') # cell in the middle lrud = fvm_2d.get_compass(1, 1, N) # neighbours must be assert np.all(lrud == np.array([0, 2, 2, 0], dtype='int')) # cell on the east lrud = fvm_2d.get_compass(2, 1, N) # neighbours must be assert np.all(lrud == np.array([1, 0, 2, 0], dtype='int')) # cell on the west lrud = fvm_2d.get_compass(0, 1, N) # neighbours must be assert np.all(lrud == np.array([2, 1, 2, 0], dtype='int')) # cell above lrud = fvm_2d.get_compass(1, 2, N) # neighbours must be assert np.all(lrud == np.array([0, 2, 0, 1], dtype='int')) # cell below lrud = fvm_2d.get_compass(1, 0, N) # neighbours must be assert np.all(lrud == np.array([0, 2, 1, 2], dtype='int')) # cell top west lrud = fvm_2d.get_compass(0, 2, N) # neighbours must be assert np.all(lrud == np.array([2, 1, 0, 1], dtype='int')) # cell top east lrud = fvm_2d.get_compass(2, 2, N) # neighbours must be assert np.all(lrud == np.array([1, 0, 0, 1], dtype='int')) # cell bot west lrud = fvm_2d.get_compass(0, 0, N) # neighbours must be assert np.all(lrud == np.array([2, 1, 1, 2], dtype='int')) # cell bot east lrud = fvm_2d.get_compass(2, 0, N) # neighbours must be assert np.all(lrud == np.array([1, 0, 1, 2], dtype='int')) def test_upwind(): # say we have a flow from west to east # the east coefficient of P should be zero # the west coefficient of the eastern neighbour should be mdot_e N = 2 a_E = np.zeros((N, N)) a_W = np.zeros((N, N)) a_N = np.zeros((N, N)) a_S = np.zeros((N, N)) source = np.zeros((N, N)) l = 1 (X, Y, XF, YF, dx, Area, V) = fvm_2d.create_uniform_grid(N, l) # dummy variables umag = 0 a = 0 tau_p = 0 # just need neighbours u, v, w, neighbours = fvm_2d.velocities_and_neighbours(X, Y, N, umag, a, tau_p) u[0, :, 0] = 1 u[1, :, 0] = 3 Area = 1 alpha = 0 i = 0 j = 0 west, east, north, south = fvm_2d.get_compass(i, j, N) fvm_2d.upwind(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W, a_E, a_S, a_N, alpha, neighbours, source) assert a_E[0, 0] == 0.0 assert a_W[1, 0] > 0.0 # the volumetric flow rate coefficient on that west face should be 2 a_W[1, 0] = 2.0 def test_linear_upwind(): # in a uniform flow, linear should provide the same coefficients as upwind l = 1 N = 20 (X, Y, XF, YF, dx, Area, V) = fvm_2d.create_uniform_grid(N, l) alpha = np.zeros_like(X) umag = 0.0 a = 0.0 tau_p = 0.0 # just need neighbours u, v, w, neighbours = fvm_2d.velocities_and_neighbours(X, Y, N, umag, a, tau_p) u = np.ones((N, N, 2), dtype='float') a_E = np.zeros((N, N)) a_W = np.zeros((N, N)) a_N = np.zeros((N, N)) a_S = np.zeros((N, N)) source = np.zeros((N, N)) a_E_linear = np.zeros((N, N)) a_W_linear = np.zeros((N, N)) a_N_linear = np.zeros((N, N)) a_S_linear = np.zeros((N, N)) source_linear = np.zeros((N, N)) Area = 1 for i in range(N): for j in range(N): west, east, north, south = neighbours[i, j] fvm_2d.upwind(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W, a_E, a_S, a_N, alpha, neighbours, source) fvm_2d.linear_upwind(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W_linear, a_E_linear, a_S_linear, a_N_linear, alpha, neighbours, source_linear) # check coefficients are the same assert np.all(a_E == a_E_linear) assert np.all(a_W == a_W_linear) assert np.all(a_N == a_N_linear) assert np.all(a_S == a_S_linear) # should be no source term for either upwind or linear assert np.all(source == 0.0) assert np.all(source_linear == 0.0) def test_linear_upwind_min_mod(): # in a uniform flow, linear min_mod should provide the same coefficients as # upwind l = 1 N = 20 (X, Y, XF, YF, dx, Area, V) = fvm_2d.create_uniform_grid(N, l) alpha = np.zeros_like(X) umag = 0.0 a = 0.0 tau_p = 0.0 # just need neighbours u, v, w, neighbours = fvm_2d.velocities_and_neighbours(X, Y, N, umag, a, tau_p) u = np.ones((N, N, 2), dtype='float') a_E = np.zeros((N, N)) a_W = np.zeros((N, N)) a_N = np.zeros((N, N)) a_S = np.zeros((N, N)) source = np.zeros((N, N)) a_E_linearm = np.zeros((N, N)) a_W_linearm = np.zeros((N, N)) a_N_linearm = np.zeros((N, N)) a_S_linearm = np.zeros((N, N)) source_linearm = np.zeros((N, N)) Area = 1 for i in range(N): for j in range(N): west, east, north, south = neighbours[i, j] fvm_2d.upwind(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W, a_E, a_S, a_N, alpha, neighbours, source) fvm_2d.linear_upwind_min_mod(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W_linearm, a_E_linearm, a_S_linearm, a_N_linearm, alpha, neighbours, source_linearm) # check coefficients are the same assert np.all(a_E == a_E_linearm) assert np.all(a_W == a_W_linearm) assert np.all(a_N == a_N_linearm) assert np.all(a_S == a_S_linearm) # should be no source term for either upwind or linear assert np.all(source == 0.0) assert np.all(source_linearm == 0.0) # in a linear flow with low flux, should return the same as linear upwind # only care about one cell i = 1 j = 1 west = 0 east = 2 north = 2 south = 0 N = 5 alpha = np.zeros((N, N)) x = np.linspace(0, 4, num=N) X, Y = np.meshgrid(x, x, indexing='ij') u = np.zeros((N, N, 2)) u[:, :, 0] = 10 ** X a_E_linear = np.zeros((N, N)) a_W_linear = np.zeros((N, N)) a_N_linear = np.zeros((N, N)) a_S_linear = np.zeros((N, N)) source_linear = np.zeros((N, N)) a_E_linearm = np.zeros((N, N)) a_W_linearm = np.zeros((N, N)) a_N_linearm = np.zeros((N, N)) a_S_linearm = np.zeros((N, N)) source_linearm = np.zeros((N, N)) fvm_2d.linear_upwind(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W_linear, a_E_linear, a_S_linear, a_N_linear, alpha, neighbours, source_linear) fvm_2d.linear_upwind_min_mod(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W_linearm, a_E_linearm, a_S_linearm, a_N_linearm, alpha, neighbours, source_linearm) assert np.abs(a_E_linear[i, j] - a_E_linearm[i, j]) < 1e-15 assert np.abs(a_W_linear[i, j] - a_W_linearm[i, j]) < 1e-15 assert np.abs(a_N_linear[i, j] - a_N_linearm[i, j]) < 1e-15 assert np.abs(a_S_linear[i, j] - a_S_linearm[i, j]) < 1e-15 assert np.abs(source_linear[i, j] - source_linearm[i, j]) < 1e-15 # a big change in flux should cause the limiter to kick in (only in x-direction) alpha[:, :] = np.logspace(4, 0, num=5) alpha = alpha.transpose() dx = 1e-6 a_E_linear = np.zeros((N, N)) a_W_linear = np.zeros((N, N)) a_N_linear = np.zeros((N, N)) a_S_linear = np.zeros((N, N)) source_linear = np.zeros((N, N)) a_E_linearm = np.zeros((N, N)) a_W_linearm = np.zeros((N, N)) a_N_linearm = np.zeros((N, N)) a_S_linearm = np.zeros((N, N)) source_linearm = np.zeros((N, N)) fvm_2d.linear_upwind(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W_linear, a_E_linear, a_S_linear, a_N_linear, alpha, neighbours, source_linear) fvm_2d.linear_upwind_min_mod(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W_linearm, a_E_linearm, a_S_linearm, a_N_linearm, alpha, neighbours, source_linearm) # as coeffs are always upwind, only the source term (deferred correction) # should be different assert np.abs(a_E_linear[i, j] - a_E_linearm[i, j]) < 1e-15 assert np.abs(a_W_linear[i, j] - a_W_linearm[i, j]) < 1e-15 assert np.abs(a_N_linear[i, j] - a_N_linearm[i, j]) < 1e-15 assert np.abs(a_S_linear[i, j] - a_S_linearm[i, j]) < 1e-15 assert np.abs(source_linear[i, j] - source_linearm[i, j]) > 1e-15 def test_linear_upwind_SUPERBEE(): # in a uniform flow, linear SUPERBEE should provide the same coefficients as # upwind l = 1 N = 20 (X, Y, XF, YF, dx, Area, V) = fvm_2d.create_uniform_grid(N, l) alpha = np.zeros_like(X) umag = 0.0 a = 0.0 tau_p = 0.0 # just need neighbours u, v, w, neighbours = fvm_2d.velocities_and_neighbours(X, Y, N, umag, a, tau_p) u = np.ones((N, N, 2), dtype='float') a_E = np.zeros((N, N)) a_W = np.zeros((N, N)) a_N = np.zeros((N, N)) a_S = np.zeros((N, N)) source = np.zeros((N, N)) a_E_linearS = np.zeros((N, N)) a_W_linearS = np.zeros((N, N)) a_N_linearS = np.zeros((N, N)) a_S_linearS = np.zeros((N, N)) source_linearS = np.zeros((N, N)) Area = 1 for i in range(N): for j in range(N): west, east, north, south = neighbours[i, j] fvm_2d.upwind(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W, a_E, a_S, a_N, alpha, neighbours, source) fvm_2d.linear_upwind_SUPERBEE(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W_linearS, a_E_linearS, a_S_linearS, a_N_linearS, alpha, neighbours, source_linearS) # check coefficients are the same assert np.all(a_E == a_E_linearS) assert np.all(a_W == a_W_linearS) assert np.all(a_N == a_N_linearS) assert np.all(a_S == a_S_linearS) # should be no source term for either upwind or linear assert np.all(source == 0.0) assert np.all(source_linearS == 0.0) # in a linear flow with low flux, should return the same as linear upwind # only care about one cell i = 1 j = 1 west = 0 east = 2 north = 2 south = 0 N = 5 alpha = np.zeros((N, N)) x = np.linspace(0, 4, num=N) X, Y = np.meshgrid(x, x, indexing='ij') u = np.zeros((N, N, 2)) u[:, :, 0] = 10 ** X a_E_linear = np.zeros((N, N)) a_W_linear = np.zeros((N, N)) a_N_linear = np.zeros((N, N)) a_S_linear = np.zeros((N, N)) source_linear = np.zeros((N, N)) a_E_linearS = np.zeros((N, N)) a_W_linearS = np.zeros((N, N)) a_N_linearS = np.zeros((N, N)) a_S_linearS = np.zeros((N, N)) source_linearS = np.zeros((N, N)) fvm_2d.linear_upwind(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W_linear, a_E_linear, a_S_linear, a_N_linear, alpha, neighbours, source_linear) fvm_2d.linear_upwind_SUPERBEE(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W_linearS, a_E_linearS, a_S_linearS, a_N_linearS, alpha, neighbours, source_linearS) assert np.abs(a_E_linear[i, j] - a_E_linearS[i, j]) < 1e-15 assert np.abs(a_W_linear[i, j] - a_W_linearS[i, j]) < 1e-15 assert np.abs(a_N_linear[i, j] - a_N_linearS[i, j]) < 1e-15 assert np.abs(a_S_linear[i, j] - a_S_linearS[i, j]) < 1e-15 assert np.abs(source_linear[i, j] - source_linearS[i, j]) < 1e-15 # a big change in flux should cause the limiter to kick in (only in x-direction) alpha[:, :] = np.logspace(4, 0, num=5) alpha = alpha.transpose() dx = 1e-6 a_E_linear = np.zeros((N, N)) a_W_linear = np.zeros((N, N)) a_N_linear = np.zeros((N, N)) a_S_linear = np.zeros((N, N)) source_linear = np.zeros((N, N)) a_E_linearS = np.zeros((N, N)) a_W_linearS = np.zeros((N, N)) a_N_linearS = np.zeros((N, N)) a_S_linearS = np.zeros((N, N)) source_linearS = np.zeros((N, N)) fvm_2d.linear_upwind(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W_linear, a_E_linear, a_S_linear, a_N_linear, alpha, neighbours, source_linear) fvm_2d.linear_upwind_SUPERBEE(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W_linearS, a_E_linearS, a_S_linearS, a_N_linearS, alpha, neighbours, source_linearS) # as coeffs are always upwind, only the source term (deferred correction) # should be different assert np.abs(a_E_linear[i, j] - a_E_linearS[i, j]) < 1e-15 assert np.abs(a_W_linear[i, j] - a_W_linearS[i, j]) < 1e-15 assert np.abs(a_N_linear[i, j] - a_N_linearS[i, j]) < 1e-15 assert np.abs(a_S_linear[i, j] - a_S_linearS[i, j]) < 1e-15 assert np.abs(source_linear[i, j] - source_linearS[i, j]) > 1e-15 def test_QUICK(): # in a uniform flow, QUICK should provide the same coefficients as upwind l = 1 N = 20 (X, Y, XF, YF, dx, Area, V) = fvm_2d.create_uniform_grid(N, l) alpha = np.zeros_like(X) umag = 0.0 a = 0.0 tau_p = 0.0 # just need neighbours u, v, w, neighbours = fvm_2d.velocities_and_neighbours(X, Y, N, umag, a, tau_p) u = np.ones((N, N, 2), dtype='float') a_E = np.zeros((N, N)) a_W = np.zeros((N, N)) a_N = np.zeros((N, N)) a_S = np.zeros((N, N)) source = np.zeros((N, N)) a_E_QUICK = np.zeros((N, N)) a_W_QUICK = np.zeros((N, N)) a_N_QUICK = np.zeros((N, N)) a_S_QUICK = np.zeros((N, N)) source_QUICK = np.zeros((N, N)) Area = 1 for i in range(N): for j in range(N): west, east, north, south = neighbours[i, j] fvm_2d.upwind(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W, a_E, a_S, a_N, alpha, neighbours, source) fvm_2d.QUICK(i, j, fvm_2d.volumetric_flow_cd_uniform, u, Area, X, Y, dx, umag, a, west, east, south, north, a_W_QUICK, a_E_QUICK, a_S_QUICK, a_N_QUICK, alpha, neighbours, source_QUICK) # check coefficients are the same assert np.all(a_E == a_E_QUICK) assert np.all(a_W == a_W_QUICK) assert np.all(a_N == a_N_QUICK) assert np.all(a_S == a_S_QUICK) # should be no source term for either upwind or QUICK assert np.all(source == 0.0) assert np.all(source_QUICK == 0.0) def test_central_difference_uniform(): # saves pasting these into every function x = np.array([1.0, 3.0]) X, Y = np.meshgrid(x, x, indexing='ij') xf = np.array([0.0, 2.0, 4.0]) XF, YF = np.meshgrid(xf, xf) # fluid velocity u = np.zeros((2, 2, 2), dtype='float') for i in range(2): for j in range(2): u[i, j, 0] = X[i, j] # u_x = x u[i, j, 1] = 0.0 # u_y = 0 # we will test the XF = 2.0 face P_i = 0 n_i = 1 P_j = 0 n_j = 0 # first with constant phi ans = 1.0 phi = np.ones_like(X) * 1.0 cd = fvm_2d.central_difference_uniform(P_i, n_i, P_j, n_j, phi) assert np.abs(ans - cd) < 1e-14 # then with linear variation of phi from 1.0 to 2.0 phi[1, :] = 2.0 ans = 1.5 cd = fvm_2d.central_difference_uniform(P_i, n_i, P_j, n_j, phi) assert np.abs(ans - cd) < 1e-14 def test_volumetric_flow_cd_uniform(): # saves pasting these into every function x = np.array([1.0, 3.0]) X, Y = np.meshgrid(x, x, indexing='ij') xf = np.array([0.0, 2.0, 4.0]) XF, YF = np.meshgrid(xf, xf) # fluid velocity u = np.zeros((2, 2, 2), dtype='float') for i in range(2): for j in range(2): u[i, j, 0] = X[i, j] # u_x = x u[i, j, 1] = 0.0 # u_y = 0 dx = x[1] - x[0] # m # we will test the XF = 2.0 face P_i = 0 n_i = 1 P_j = 0 n_j = 0 # linear velocity gradient: 1m/s at P and 3m/s at E u = np.zeros((2, 2, 2)) umag = 1 a = 1 u[0, :, 0] = 1 u[1, :, 0] = 3 Area = 1 ans = 2.0 * Area alpha = 0.0 vdot = fvm_2d.volumetric_flow_cd_uniform(P_i, P_j, n_i, n_j, u, alpha, 'east', Area, X, Y, dx, umag, a) assert np.abs(ans - vdot) < 1e-14 # should give the same result as top east neighbour vdot = fvm_2d.volumetric_flow_cd_uniform(P_i, P_j, n_i, 1, u, alpha, 'east', Area, X, Y, dx, umag, a) assert np.abs(ans - vdot) < 1e-14 def test_volumetric_flow_exact(): dx = 2.0 umag = 1.0 a = 1.0 Area = dx * dx # zero at the origin on the north face x = np.array([0.0, 0.0]) y = np.array([-5.0, 5.0]) X, Y = np.meshgrid(x, y, indexing='ij') u = np.zeros((2, 2, 2)) alpha = 0.0 vdot = fvm_2d.volumetric_flow_exact(0, 0, 0, 1, u, alpha, 'north', Area, X, Y, dx, umag, a) assert np.abs(vdot) < 1e-14 def test_charge_flux(): # test that zero volume fraction returns zero dx = 2.0 umag = 1.0 a = 1.0 N = 2 l = np.pi / a tau_p = 1.0 (X, Y, XF, YF, dx, Area, V) = fvm_2d.create_uniform_grid(N, l) u, v, w, neighbours = fvm_2d.velocities_and_neighbours(X, Y, N, umag, a, tau_p) alpha = np.zeros((N, N)) flux = fvm_2d.charge_flux(0, 0, 0, 1, u, alpha, 'north', Area, X, Y, dx, umag, a) assert np.abs(flux) < 1e-14 alpha = np.ones((N, N)) flux = fvm_2d.charge_flux(0, 0, 0, 1, u, alpha, 'north', Area, X, Y, dx, umag, a) assert np.abs(flux) > 1e-14 def test_source_uniform_grid_explicit(): # lets have zero volume fraction everywhere alpha = np.zeros((3, 3), dtype='float') source = np.zeros_like(alpha) S_p = np.zeros_like(alpha) w = np.ones((3, 3, 2)) i = 1 j = 1 rho_f = 1.0 Area = 1.0 west, east, north, south = fvm_2d.get_compass(i, j, 3) fvm_2d.source_uniform_grid_explicit(i, j, east, west, north, south, w, alpha, Area, source, S_p, limit=0.0) assert np.all(np.abs(source) < 1e-14) assert np.all(np.abs(S_p) < 1e-14) # uniform volume fraction distribution for the TG flow. The source term then # has an analytical solution given by # S = -rho_f tau_p umag a^2 alpha dV (cos(2ax) + cos(2ay)) def source_analytical(x, y, rho_f, tau_p, umag, a, alpha0, V): return -rho_f * tau_p * umag * (a ** 2) * alpha0 * V * (np.cos(2 * a * x) + np.cos(2 * a * y)) # check that the discretised solution approaches the analytical solution # as the grid is refined N_array = np.array([32, 64, 128, 256, 512, 1024]) norms = np.zeros_like(N_array, dtype='float') for norm, N in enumerate(N_array): x = np.linspace(-np.pi, np.pi, num=N) X, Y = np.meshgrid(x, x, indexing='ij') alpha0 = 1.0 alpha = np.ones_like(X) * alpha0 tau_p = 1.0 rho_f = 1.0 a = 1.0 umag = 1.0 dx = x[1] - x[0] Area = dx ** 2 V = dx ** 3 u, v, w, neighbours = fvm_2d.velocities_and_neighbours(X, Y, N, umag, a, tau_p) ans = np.zeros_like(X) S_p = np.zeros_like(X) disc = np.zeros_like(X) for i in range(N): for j in range(N): west, east, north, south = neighbours[i, j] ans[i, j] = source_analytical(X[i, j], Y[i, j], rho_f, tau_p, umag, a, alpha0, V) fvm_2d.source_uniform_grid_explicit(i, j, east, west, north, south, w, alpha, Area, disc, S_p, limit=0.0) assert np.all(S_p == 0.0) norms[norm] = np.linalg.norm(ans - disc) # ensure that the norm decreases with increasing grid points assert np.all(norms[1:] < norms[:-1]) # ensure that the finest mesh is 'close' to the analytical solution assert norms[-1] < 1e-5 def test_source_uniform_grid_implicit(): # lets have zero volume fraction everywhere alpha = np.zeros((3, 3), dtype='float') S_c = np.zeros_like(alpha) S_p = np.zeros_like(alpha) w = np.ones((3, 3, 2)) i = 1 j = 1 rho_f = 1.0 Area = 1.0 west, east, north, south = fvm_2d.get_compass(i, j, 3) fvm_2d.source_uniform_grid_implicit(i, j, east, west, north, south, w, alpha, Area, S_c, S_p, limit=0.0) source = S_c + (S_p * 0) assert np.all(np.abs(source) < 1e-14) # uniform volume fraction distribution for the TG flow. The source term then # has an analytical solution given by # S = -rho_f tau_p umag a^2 alpha dV (cos(2ax) + cos(2ay)) def source_analytical(x, y, rho_f, tau_p, umag, a, alpha0, V): return -rho_f * tau_p * umag * (a ** 2) * alpha0 * V * (np.cos(2 * a * x) + np.cos(2 * a * y)) # check that the discretised solution approaches the analytical solution # as the grid is refined N_array = np.array([32, 64, 128, 256, 512, 1024]) norms = np.zeros_like(N_array, dtype='float') for norm, N in enumerate(N_array): x = np.linspace(-np.pi, np.pi, num=N) X, Y = np.meshgrid(x, x, indexing='ij') alpha0 = 1.0 alpha = np.ones_like(X) * alpha0 tau_p = 1.0 rho_f = 1.0 a = 1.0 umag = 1.0 dx = x[1] - x[0] Area = dx ** 2 V = dx ** 3 u, v, w, neighbours = fvm_2d.velocities_and_neighbours(X, Y, N, umag, a, tau_p) ans = np.zeros_like(X) S_c = np.zeros_like(X) S_p = np.zeros_like(X) disc = np.zeros_like(X) for i in range(N): for j in range(N): west, east, north, south = neighbours[i, j] ans[i, j] = source_analytical(X[i, j], Y[i, j], rho_f, tau_p, umag, a, alpha0, V) fvm_2d.source_uniform_grid_implicit(i, j, east, west, north, south, w, alpha, Area, S_c, S_p, limit=0.0) disc[i, j] = S_c[i, j] + (S_p[i, j] * alpha0) norms[norm] = np.linalg.norm(ans - disc) # ensure that the norm decreases with increasing grid points assert np.all(norms[1:] < norms[:-1]) # ensure that the finest mesh is 'close' to the analytical solution assert norms[-1] < 1e-5 return N_array, norms def test_calculate_a_P(): # say a 2x2 domain N = 2 V = None dt = None a_E = np.array([[1.0, 2.0], [3.0, 4.0]]) a_W = np.array([[2.0, 3.0], [4.0, 5.0]]) a_N = np.array([[3.0, 4.0], [5.0, 6.0]]) a_S = np.array([[4.0, 5.0], [6.0, 7.0]]) S_p = np.array([[-0.1, -0.2], [-0.3, -0.4]]) a_P = np.zeros_like(a_E) # a_P should equal the sum of neighbours minus S_p ans = np.array([[10.1, 14.2], [18.3, 22.4]]) a_P = fvm_2d.calculate_a_P(a_P, a_E, a_W, a_N, a_S, S_p, N, V, dt) assert np.allclose(ans, a_P) def test_Gauss_Seidel_point_by_point(): N = 2 # essentially simulating convection from west to east # expecting a uniform distribution as solution x = np.array([0.0, 0.0]) X, Y = np.meshgrid(x, x, indexing='ij') u, v, w, neighbours = fvm_2d.velocities_and_neighbours(X, Y, N, 0, 1, 1) a_P = np.array([[1.0, 1.0], [1.0, 1.0]]) a_E = np.zeros_like(a_P) a_N = np.zeros_like(a_P) a_S = np.zeros_like(a_P) a_W = np.array([[0, 0], [1, 1]]) b = np.array([[1.0, 1.0], [0.0, 0.0]]) # put them all together for Ac=b A = np.array([[1.0, 0, 0, 0], [0, 1, 0, 0], [0, -1, 1, 0], [0, 0, -1, 1]]) # answer = (1, 1, 1, 1) ans = np.linalg.solve(A, b.ravel()) # initial guess alpha = np.zeros_like(a_P) disc = fvm_2d.Gauss_Seidel_point_by_point(alpha, a_P, a_E, a_W, a_N, a_S, b, N, neighbours) assert np.allclose(ans, disc.ravel()) def test_Gauss_Seidel_point_by_point_relax(): N = 2 # essentially simulating convection from west to east # expecting a uniform distribution as solution x = np.array([0.0, 0.0]) X, Y = np.meshgrid(x, x, indexing='ij') u, v, w, neighbours = fvm_2d.velocities_and_neighbours(X, Y, N, 0, 1, 1) a_P = np.array([[1.0, 1.0], [1.0, 1.0]]) a_E = np.zeros_like(a_P) a_N = np.zeros_like(a_P) a_S = np.zeros_like(a_P) a_W = np.array([[0, 0], [1, 1]]) b = np.array([[1.0, 1.0], [0.0, 0.0]]) # put them all together for Ac=b A = np.array([[1.0, 0, 0, 0], [0, 1, 0, 0], [0, -1, 1, 0], [0, 0, -1, 1]]) # answer = (1, 1, 1, 1) ans = np.linalg.solve(A, b.ravel()) # initial guess alpha = np.zeros_like(a_P) relax = 1.0 disc = fvm_2d.Gauss_Seidel_point_by_point_relax(alpha, a_P, a_E, a_W, a_N, a_S, b, N, neighbours, relax) assert np.allclose(ans, disc.ravel()) def test_sources_converge_to_same_solution(): # test using a low Stokes number (N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p) = fvm_2d.initialise(N=128, Stk=1e-3, alpha_0=5e-7, rho_f=1.0, rho_p=2650.0, umag=1.543, a=84765.0, flux_func=fvm_2d.volumetric_flow_exact, source_func=fvm_2d.source_uniform_grid_implicit, diff_func=fvm_2d.diffusion_coefficient_none) (alpha_imp, residual_imp, alpha_max_imp, alpha_min_imp, alpha_tot_imp, mean_alpha, moment_2, moment_3, moment_4, a_P, a_E, a_W, a_N, a_S, source, S_p) = fvm_2d.it_resi(N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p, 10000, relax=1.0, limit=0.0, flux_func=fvm_2d.volumetric_flow_exact, source_func=fvm_2d.source_uniform_grid_implicit, diff_func=fvm_2d.diffusion_coefficient_none, res_criterion=1e-14) # check it did converge assert residual_imp[-1] < 1e-14 # now calculate the same solution with an explicit source (N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p) = fvm_2d.initialise(N=128, Stk=1e-3, alpha_0=5e-7, rho_f=1.0, rho_p=2650.0, umag=1.543, a=84765, source_func=fvm_2d.source_uniform_grid_explicit, diff_func=fvm_2d.diffusion_coefficient_none) (alpha_exp, residual_exp, alpha_max_exp, alpha_min_exp, alpha_tot_exp, mean_alpha, moment_2, moment_3, moment_4, a_P, a_E, a_W, a_N, a_S, source, S_p) = fvm_2d.it_resi(N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p, 10000, relax=1.0, limit=0.0, source_func=fvm_2d.source_uniform_grid_explicit, diff_func=fvm_2d.diffusion_coefficient_none, res_criterion=1e-14) # check it did converge assert residual_exp[-1] < 1e-14 # check solutions are the same (super close) assert np.allclose(alpha_imp, alpha_exp, atol=1e-6) assert np.abs(alpha_max_imp[-1] - alpha_max_exp[-1]) < 2e-7 assert np.abs(alpha_min_imp[-1] - alpha_min_exp[-1]) < 1e-7 # total mass fraction differs a bit more assert np.abs(alpha_tot_imp[-1] - alpha_tot_exp[-1]) < 1e-3 def test_volume_fraction_to_mass_fraction(): alpha = 1e-7 rho_f = 1.0 rho_p = 1000.0 c_ans = 9.999e-5 c = fvm_2d.volume_fraction_to_mass_fraction(alpha, rho_p, rho_f) assert np.abs(c - c_ans) < 1e-10 def test_mass_fraction_to_volume_fraction(): c = 9.999001099790129e-05 rho_f = 1.0 rho_p = 1000.0 alpha_ans = 1e-7 alpha = fvm_2d.mass_fraction_to_volume_fraction(c, rho_p, rho_f) assert np.abs(alpha - alpha_ans) < 1e-10 def test_volume_fraction_to_number_density(): alpha = 0.0 d_p = 1e-7 assert np.abs(fvm_2d.volume_fraction_to_number_density(alpha, d_p)) < 1e-15 def test_number_density_to_volume_fraction(): n = 0.0 d_p = 1e-7 assert np.abs(fvm_2d.number_density_to_volume_fraction(n, d_p)) < 1e-15 def test_hertzian_max_contact_area(): # zero particle mass should return zero v_rel = 1.0 d_p = 1e-6 rho_p = 0.0 A_max = fvm_2d.hertzian_max_contact_area(v_rel, d_p, rho_p) assert A_max == 0.0 # zero particle velocity should return zero v_rel = 0.0 d_p = 1e-6 rho_p = 2650.0 A_max = fvm_2d.hertzian_max_contact_area(v_rel, d_p, rho_p) assert A_max == 0.0 def test_hertzian_contact_time(): # zero particle diameter should return zero v_rel = 1.0 d_p = 0.0 rho_p = 2650.0 t_c = fvm_2d.hertzian_contact_time(v_rel, d_p, rho_p) assert t_c == 0.0 # zero particle density should return zero v_rel = 1.0 d_p = 1e-6 rho_p = 0.0 t_c = fvm_2d.hertzian_contact_time(v_rel, d_p, rho_p) assert t_c == 0.0 # if these quantities equal unity v_rel = 1.0 d_p = 1.0 rho_p = 1.0 t_c = fvm_2d.hertzian_contact_time(v_rel, d_p, rho_p) I = 1.4716375921623521 # integral (hertzian_integral.py) E = 7.13e10 # Young's modulus (Pa) nu = 0.17 # Poisson's ratio ans = I * (((5 * np.pi * np.sqrt(2) / 4) * ((1 - (nu ** 2)) / E)) ** 0.4) assert np.abs(t_c - ans) < 1e-14 def test_ireland_charge_transfer(): # if A_c, sigma_0 or t_c is zero, it should return zero A_c = 0.0 sigma_0 = 1.0 tau = 1.0 t_c = 1.0 assert fvm_2d.ireland_charge_transfer(A_c, sigma_0, tau, t_c) == 0.0 A_c = 1.0 sigma_0 = 0.0 tau = 1.0 t_c = 1.0 assert fvm_2d.ireland_charge_transfer(A_c, sigma_0, tau, t_c) == 0.0 A_c = 1.0 sigma_0 = 1.0 tau = 1.0 t_c = 0.0 assert fvm_2d.ireland_charge_transfer(A_c, sigma_0, tau, t_c) == 0.0 # simple check (solution = 1 - (1 / e)) A_c = 1.0 sigma_0 = 1.0 tau = 1.0 t_c = 1.0 Delta_Q = fvm_2d.ireland_charge_transfer(A_c, sigma_0, tau, t_c) ans = 1 - np.exp(-1) assert np.abs(Delta_Q - ans) < 1e-14 def test_taylor_green_ee_strain_rate(): # if tau_p=0, we should just have the fluid strain rate magnitude a = 1.0 umag = 1.0 tau_p = 0.0 # this is zero at the origin x = 0.0 y = 0.0 S_mag_v = fvm_2d.taylor_green_ee_strain_rate(x, y, umag, a, tau_p) assert np.abs(S_mag_v) < 1e-10 # maxima occur at (x, y) = (n * pi, m * pi) # with the strain rate magnitude equalling 2 * umag * a x = np.pi / 2 y = np.pi / 2 S_mag_v = fvm_2d.taylor_green_ee_strain_rate(x, y, umag, a, tau_p) assert np.abs(S_mag_v - 2.0) < 1e-10 # for non-zero tau_p, the origin will provide a non-zero strain rate x = 0.0 y = 0.0 tau_p = 1.0 S_mag_v = fvm_2d.taylor_green_ee_strain_rate(x, y, umag, a, tau_p) assert np.abs(S_mag_v) > 1e-10 # test at (x, y) = (0, 0) # analytical solution is 2 tau_p umag^2 a^2 def zero_analytical(umag, a, tau_p): return 2 * tau_p * (umag ** 2) * (a ** 2) x = 0.0 y = 0.0 umag = 1.0 a = 1.0 tau_p = 1.0 S_mag_v = fvm_2d.taylor_green_ee_strain_rate(x, y, umag, a, tau_p) assert np.abs(S_mag_v - zero_analytical(umag, a, tau_p)) < 1e-15 umag = 15.7 a = 1e-7 tau_p = 1e-6 S_mag_v = fvm_2d.taylor_green_ee_strain_rate(x, y, umag, a, tau_p) assert np.abs(S_mag_v - zero_analytical(umag, a, tau_p)) < 1e-15 # test at (x, y) = (pi / 2a, pi / 2a) # analytical solution is 2 sqrt( umag^2 a^2 (1 + tau_p^2 umag^2 a^2)) def pi_over_two_analytical(umag, a, tau_p): return 2 * np.sqrt(((umag ** 2) * (a ** 2)) + ((tau_p ** 2) * (umag ** 4) * (a ** 4))) umag = 1.0 a = 1.0 tau_p = 1.0 x = np.pi / (2 * a) y = np.pi / (2 * a) S_mag_v = fvm_2d.taylor_green_ee_strain_rate(x, y, umag, a, tau_p) assert np.abs(S_mag_v - pi_over_two_analytical(umag, a, tau_p)) < 1e-15 umag = 15.7 a = 1e-7 tau_p = 1e-6 x = np.pi / (2 * a) y = np.pi / (2 * a) S_mag_v = fvm_2d.taylor_green_ee_strain_rate(x, y, umag, a, tau_p) assert np.abs(S_mag_v - pi_over_two_analytical(umag, a, tau_p)) < 1e-15 def test_particle_mean_free_path(): # expected to be proportional to the particle diameter d_p = 0.0 alpha = 1e-7 assert np.abs(fvm_2d.particle_mean_free_path(d_p, alpha)) < 1e-15 def test_packing_function(): # should return unity if alpha = 0 assert fvm_2d.packing_function(0.0) == 1.0 # should increase with increasing volume fraction alpha = np.logspace(-10, -0.5, num=100) g0 = np.zeros_like(alpha) for i in range(len(alpha)): g0[i] = fvm_2d.packing_function(alpha[i]) assert np.all(g0[1:] >= g0[:-1]) def independent_of_density(): """ This test checks that the same solution is reached, regardless of the choice of particle density. This is expected as the particles have constant density which will mean that rho_p should cancel in each term of the governing equation. Returns ------- None. """ # checking with explicit source term rho_p = 1.0 (N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p) = fvm_2d.initialise(N=128, Stk=1e-3, alpha_0=5e-7, rho_f=1.0, rho_p=rho_p, umag=1.543, a=84765.0, source_func=fvm_2d.source_uniform_grid_explicit, diff_func=fvm_2d.diffusion_coefficient_none) (alpha1, residual1, alpha_max1, alpha_min1, alpha_tot1, a_P, a_E, a_W, a_N, a_S, source, S_p) = fvm_2d.it_resi(N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p, 10000, relax=1.0, limit=0.0, source_func=fvm_2d.source_uniform_grid_explicit, diff_func=fvm_2d.diffusion_coefficient_none) rho_p = 2650.0 (N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p) = fvm_2d.initialise(N=128, Stk=1e-3, alpha_0=5e-7, rho_f=1.0, rho_p=rho_p, umag=1.543, a=84765.0, source_func=fvm_2d.source_uniform_grid_explicit, diff_func=fvm_2d.diffusion_coefficient_none) (alpha2, residual2, alpha_max2, alpha_min2, alpha_tot1, a_P, a_E, a_W, a_N, a_S, source, S_p) = fvm_2d.it_resi(N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p, 10000, relax=1.0, limit=0.0, source_func=fvm_2d.source_uniform_grid_explicit, diff_func=fvm_2d.diffusion_coefficient_none) # check both solutions are the same assert np.allclose(alpha, alpha2) # now check with the implicit source term rho_p = 1.0 (N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p) = fvm_2d.initialise(N=128, Stk=1e-3, alpha_0=5e-7, rho_f=1.0, rho_p=rho_p, umag=1.543, a=84765.0, source_func=fvm_2d.source_uniform_grid_implicit, diff_func=fvm_2d.diffusion_coefficient_none) (alpha1, residual1, alpha_max1, alpha_min1, alpha_tot1, a_P, a_E, a_W, a_N, a_S, source, S_p) = fvm_2d.it_resi(N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p, 10000, relax=1.0, limit=0.0, source_func=fvm_2d.source_uniform_grid_implicit, diff_func=fvm_2d.diffusion_coefficient_none) rho_p = 2650.0 (N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p) = fvm_2d.initialise(N=128, Stk=1e-3, alpha_0=5e-7, rho_f=1.0, rho_p=rho_p, umag=1.543, a=84765.0, source_func=fvm_2d.source_uniform_grid_implicit, diff_func=fvm_2d.diffusion_coefficient_none) (alpha2, residual2, alpha_max2, alpha_min2, alpha_tot1, a_P, a_E, a_W, a_N, a_S, source, S_p) = fvm_2d.it_resi(N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p, 10000, relax=1.0, limit=0.0, source_func=fvm_2d.source_uniform_grid_implicit, diff_func=fvm_2d.diffusion_coefficient_none) # check both solutions are the same assert np.allclose(alpha, alpha2) def test_first_order_backward_euler(): V = 1.0 dt = 1.0 ans = 1.0 coeff = fvm_2d.first_order_backward_euler(V, dt) assert np.abs(ans - coeff) < 1e-15 def test_first_order_unsteady_progresses_single_time_step(): """ Checks that the unsteady solver is only solving for a single time step and not continuing through time Returns ------- None. """ (N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p, dt) = fvm_2d.initialise_unsteady(N=128, Stk=1e-3, alpha_0=5e-7, rho_f=1.0, rho_p=2650.0, umag=1.543, a=84765, source_func=fvm_2d.source_uniform_grid_implicit, diff_func=fvm_2d.diffusion_coefficient_collisions, time_func=fvm_2d.first_order_backward_euler, dt = 1e-6) # 100 iterations on a single time step (implicit source) Nit = 100 NT = 1 alpha, residual, alpha_max, alpha_min, alpha_tot, a_P, a_E, a_W, a_N, a_S, source, S_p = fvm_2d.time_iterate(N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p, Nit, NT, dt, relax=1.0, limit=0.0, source_func=fvm_2d.source_uniform_grid_implicit, diff_func=fvm_2d.diffusion_coefficient_none, time_func=fvm_2d.first_order_backward_euler) original = np.ones_like(alpha) * 5e-7 assert np.allclose(alpha, original) # 100 iterations on a single time step (explicit source) (N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p, dt) = fvm_2d.initialise_unsteady(N=128, Stk=1e-3, alpha_0=5e-7, rho_f=1.0, rho_p=2650.0, umag=1.543, a=84765, source_func=fvm_2d.source_uniform_grid_explicit, diff_func=fvm_2d.diffusion_coefficient_collisions, time_func=fvm_2d.first_order_backward_euler, dt = 1e-6) alpha, residual, alpha_max, alpha_min, alpha_tot, a_P, a_E, a_W, a_N, a_S, source, S_p = fvm_2d.time_iterate(N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p, Nit, NT, dt, relax=1.0, limit=0.0, source_func=fvm_2d.source_uniform_grid_explicit, diff_func=fvm_2d.diffusion_coefficient_none, time_func=fvm_2d.first_order_backward_euler) original = np.ones_like(alpha) * 5e-7 assert np.allclose(alpha, original) def test_unsteady_converges_to_steady_state(): # steady state solution (N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p) = fvm_2d.initialise(N=128, Stk=1e-3, alpha_0=5e-7, rho_f=1.0, rho_p=2650.0, umag=1.543, a=84765, source_func=fvm_2d.source_uniform_grid_implicit, diff_func=fvm_2d.diffusion_coefficient_none) (alpha_steady, residual_steady, alpha_max_steady, alpha_min_steady, alpha_tot_steady, mean_alpha, moment_2, moment_3, moment_4, a_P, a_E, a_W, a_N, a_S, source, S_p) = fvm_2d.it_resi(N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p, 7000, relax=1.0, limit=0.0, source_func=fvm_2d.source_uniform_grid_implicit, diff_func=fvm_2d.diffusion_coefficient_none) # now iterate in time dt = 1e-4 (N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p, dt) = fvm_2d.initialise_unsteady(N=128, Stk=1e-3, alpha_0=5e-7, rho_f=1.0, rho_p=2650.0, umag=1.543, a=84765, source_func=fvm_2d.source_uniform_grid_implicit, diff_func=fvm_2d.diffusion_coefficient_collisions, time_func=fvm_2d.first_order_backward_euler, dt=dt) alpha, residual, alpha_max, alpha_min, alpha_tot, a_P, a_E, a_W, a_N, a_S, source, S_p = fvm_2d.time_iterate(N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p, 10, 1000, dt, relax=1.0, limit=0.0, source_func=fvm_2d.source_uniform_grid_implicit, diff_func=fvm_2d.diffusion_coefficient_none, time_func=fvm_2d.first_order_backward_euler) # check key numbers are very similar assert np.abs(alpha_max[-1] - alpha_max_steady[-1]) < 1e-7 assert np.abs(alpha_min[-1] - alpha_min_steady[-1]) < 1e-7 # harder one to satisfy assert np.linalg.norm(alpha - alpha_steady) < 2e-7 def test_source_linear_analytical(): """ Checks the source term discretisation result against analytical results. These analytical results are due to linear variations, therefore the discretisation result is expected to be exact. Returns ------- None. """ # cell centred at P with the 8 surrounding cells # uniform phi, linear relative velocity gradient in x-direction phi_0 = 1.0 phi = np.ones((3, 3), dtype='float') * phi_0 S_c = np.zeros_like(phi) S_p = np.zeros_like(phi) w = np.zeros((3, 3, 2), dtype='float') w[0, :, 0] = -1 w[2, :, 0] = 1 i = 1 j = 1 west = 0 south = 0 east = 2 north = 2 Area = 1.0 V = 1.0 fvm_2d.source_uniform_grid_explicit(i, j, east, west, north, south, w, phi, Area, S_c, S_p, limit=0.0) assert S_p[1, 1] == 0.0 S_true = -phi_0 * V assert np.abs(S_c[1, 1] - S_true) < 1e-14 # test the y-direction # uniform phi, linear relative velocity gradient in y-direction S_c = np.zeros_like(phi) S_p = np.zeros_like(phi) w = np.zeros((3, 3, 2), dtype='float') w[:, 0, 1] = -1 w[:, 2, 1] = 1 fvm_2d.source_uniform_grid_explicit(i, j, east, west, north, south, w, phi, Area, S_c, S_p, limit=0.0) assert S_p[1, 1] == 0.0 S_true = -phi_0 * V assert np.abs(S_c[1, 1] - S_true) < 1e-14 # test both x and y linear velocity gradient (each unity again) S_c = np.zeros_like(phi) S_p = np.zeros_like(phi) w[0, :, 0] = -1 w[2, :, 0] = 1 fvm_2d.source_uniform_grid_explicit(i, j, east, west, north, south, w, phi, Area, S_c, S_p, limit=0.0) assert S_p[1, 1] == 0.0 S_true = -2 * phi_0 * V assert np.abs(S_c[1, 1] - S_true) < 1e-14 # check with different velocity gradients (du/dx=2, dv/dy=3) S_c = np.zeros_like(phi) S_p = np.zeros_like(phi) w[0, :, 0] = -2 w[2, :, 0] = 2 w[:, 0, 1] = -3 w[:, 2, 1] = 3 fvm_2d.source_uniform_grid_explicit(i, j, east, west, north, south, w, phi, Area, S_c, S_p, limit=0.0) assert S_p[1, 1] == 0.0 S_true = -5 * phi_0 * V assert np.abs(S_c[1, 1] - S_true) < 1e-14 # check with velocity gradients which should cause the source to cancel S_c = np.zeros_like(phi) S_p = np.zeros_like(phi) w[0, :, 0] = 3 w[2, :, 0] = -3 w[:, 0, 1] = -3 w[:, 2, 1] = 3 fvm_2d.source_uniform_grid_explicit(i, j, east, west, north, south, w, phi, Area, S_c, S_p, limit=0.0) assert S_p[1, 1] == 0.0 S_true = 0.0 assert S_c[1, 1] == S_true # check with constant relative velocity, but linear volume fraction S_c = np.zeros_like(phi) S_p = np.zeros_like(phi) w = np.ones((3, 3, 2), dtype='float') phi[0, :] = 0.0 phi[2, :] = 2.0 fvm_2d.source_uniform_grid_explicit(i, j, east, west, north, south, w, phi, Area, S_c, S_p, limit=0.0) assert S_p[1, 1] == 0.0 S_true = -1.0 assert S_c[1, 1] == S_true def test_no_source_leaves_flow_field_unchanged(): # check upwind (N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p) = fvm_2d.initialise(N=128, Stk=1e-3, alpha_0=5e-7, rho_f=1.0, rho_p=2650.0, umag=1.543, a=84765, flux_func=fvm_2d.volumetric_flow_exact, source_func=fvm_2d.source_zero, diff_func=fvm_2d.diffusion_coefficient_none, conv_func=fvm_2d.upwind) (alpha, residual, alpha_max, alpha_min, alpha_tot, mean_alpha, moment_2, moment_3, moment_4, a_P, a_E, a_W, a_N, a_S, source, S_p) = fvm_2d.it_resi(N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p, 10000, relax=1.0, limit=0.0, flux_func=fvm_2d.volumetric_flow_exact, source_func=fvm_2d.source_zero, diff_func=fvm_2d.diffusion_coefficient_none, conv_func=fvm_2d.upwind, a_nb_func=fvm_2d.calculate_a_neighbours, a_P_func=fvm_2d.calculate_a_P, res_criterion=1e-14) assert np.linalg.norm(alpha - 5e-7) < 1e-14 # check QUICK (N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p) = fvm_2d.initialise(N=128, Stk=1e-3, alpha_0=5e-7, rho_f=1.0, rho_p=2650.0, umag=1.543, a=84765, flux_func=fvm_2d.volumetric_flow_exact, source_func=fvm_2d.source_zero, diff_func=fvm_2d.diffusion_coefficient_none, conv_func=fvm_2d.QUICK) (alpha, residual, alpha_max, alpha_min, alpha_tot, mean_alpha, moment_2, moment_3, moment_4, a_P, a_E, a_W, a_N, a_S, source, S_p) = fvm_2d.it_resi(N, umag, a, X, Y, dx, Area, V, alpha, u, v, w, neighbours, a_P, a_E, a_W, a_N, a_S, source, S_p, rho_f, rho_p, tau_p, d_p, 10000, relax=1.0, limit=0.0, flux_func=fvm_2d.volumetric_flow_exact, source_func=fvm_2d.source_zero, diff_func=fvm_2d.diffusion_coefficient_none, conv_func=fvm_2d.QUICK, a_nb_func=fvm_2d.calculate_a_neighbours, a_P_func=fvm_2d.calculate_a_P, res_criterion=1e-14) assert np.linalg.norm(alpha - 5e-7) < 1e-14 def test_turbulent_diffusion_coefficient(): # origin # should be zero at for zero fluid velocity x_f = 0.0 y_f = 0.0 phi_f = 1.0 umag = 0.0 a = 1.0 tau_p = 1.0 d_p = 1.0 dx = 1.0 D = fvm_2d.diffusion_coefficient_turb(x_f, y_f, phi_f, umag, a, tau_p, d_p, dx) assert D == 0.0 # should be zero for zero tau_p umag = 1.0 tau_p = 0.0 D = fvm_2d.diffusion_coefficient_turb(x_f, y_f, phi_f, umag, a, tau_p, d_p, dx) assert D == 0.0 # should be zero for zero a tau_p = 1.0 a = 0.0 D = fvm_2d.diffusion_coefficient_turb(x_f, y_f, phi_f, umag, a, tau_p, d_p, dx) assert D == 0.0 # at the origin, S_mag_v = 2 tau_p umag^2 a^2 # if everything equals unity, S_mag_v = 2, D = 0.02 a = 1.0 D = fvm_2d.diffusion_coefficient_turb(x_f, y_f, phi_f, umag, a, tau_p, d_p, dx) assert np.abs(D - 0.02) < 1e-15 # if umag equals 2, S_mag_v = 8, D = 0.08 umag = 2.0 D = fvm_2d.diffusion_coefficient_turb(x_f, y_f, phi_f, umag, a, tau_p, d_p, dx) assert np.abs(D - 0.08) < 1e-15 # if tau_p equals 3, S_mag_v = 6, D = 0.06 umag = 1.0 tau_p = 3.0 D = fvm_2d.diffusion_coefficient_turb(x_f, y_f, phi_f, umag, a, tau_p, d_p, dx) assert np.abs(D - 0.06) < 1e-15 # if a equals 4, S_mag_v = 32, D = 0.32 tau_p = 1.0 a = 4.0 D = fvm_2d.diffusion_coefficient_turb(x_f, y_f, phi_f, umag, a, tau_p, d_p, dx) assert np.abs(D - 0.32) < 1e-15 # if params = 1 but dx = 5, S_mag_v = 2, D = 0.5 (prop to dx^2) a = 1.0 dx = 5.0 D = fvm_2d.diffusion_coefficient_turb(x_f, y_f, phi_f, umag, a, tau_p, d_p, dx) assert np.abs(D - 0.5) < 1e-15 # should have the same answer independent of phi_f phi_f = 0.0 assert np.abs(fvm_2d.diffusion_coefficient_turb(x_f, y_f, phi_f, umag, a, tau_p, d_p, dx) - 0.5) < 1e-15 phi_f = 1e4 assert np.abs(fvm_2d.diffusion_coefficient_turb(x_f, y_f, phi_f, umag, a, tau_p, d_p, dx) - 0.5) < 1e-15 phi_f = -17 assert np.abs(fvm_2d.diffusion_coefficient_turb(x_f, y_f, phi_f, umag, a, tau_p, d_p, dx) - 0.5) < 1e-15 # should have the same answer independent of d_p (assuming constant tau_p) d_p = 0.0 assert np.abs(fvm_2d.diffusion_coefficient_turb(x_f, y_f, phi_f, umag, a, tau_p, d_p, dx) - 0.5) < 1e-15 d_p = 1e4 assert np.abs(fvm_2d.diffusion_coefficient_turb(x_f, y_f, phi_f, umag, a, tau_p, d_p, dx) - 0.5) < 1e-15 d_p = -17 assert np.abs(fvm_2d.diffusion_coefficient_turb(x_f, y_f, phi_f, umag, a, tau_p, d_p, dx) - 0.5) < 1e-15