From 0731d93d2a74f08b7ef82c04f6f11b68faa1f66a Mon Sep 17 00:00:00 2001 From: Leo Kotipalo Date: Fri, 5 Feb 2021 00:26:40 +0200 Subject: [PATCH 01/14] Initial commit for implementing divergence-free reconstruction of B --- scripts/divergenceless_reconstruction.py | 176 +++++++++++++++++++++++ 1 file changed, 176 insertions(+) create mode 100644 scripts/divergenceless_reconstruction.py diff --git a/scripts/divergenceless_reconstruction.py b/scripts/divergenceless_reconstruction.py new file mode 100644 index 00000000..475a3952 --- /dev/null +++ b/scripts/divergenceless_reconstruction.py @@ -0,0 +1,176 @@ +import pytools as pt +import numpy as np +from scipy.special import legendre + +P = tuple(map(lambda n: legendre(n), (0, 1, 2, 3, 4))) + +# y and z are simply cyclic rotations of this +def interpolate_x(a, x, y, z): + B_x = 0 + for i in range(5): + for j in range(4): + for k in range(4): + B_x += a[i][j][k] * P[i](x) * P[j](y) * P[k](z) + return B_x + +# Returns moments as [ijk, order, x, y, z] +def solve_moments_from_B(fg_b): + # Cyclically rotate each component and rotate coordinates back + x_moments = np.transpose(solve_x_moments(np.transpose(fg_b[:, :, :, 0], (0, 1, 2))), (0, 1, 2, 3)) + y_moments = np.transpose(solve_x_moments(np.transpose(fg_b[:, :, :, 1], (1, 2, 0))), (0, 3, 1, 2)) + z_moments = np.transpose(solve_x_moments(np.transpose(fg_b[:, :, :, 2], (2, 0, 1))), (0, 2, 3, 1)) + return (x_moments, y_moments, z_moments) + +# Returns moments for x as [order, x, y, z] +# With order (0, y, z, yy, yz, zz, yyy, yyz, yzz, zzz) +# y and z as cyclic rotations +def solve_x_moments(B_x): + x_first = np.gradient(B_x)[1:3] + x_second = [] + for i in range(2): + j = i + 1 + x_second += np.gradient(x_first[i])[j:3] + x_third = [] + for i in range(3): + if i < 2: + j = i+1 + else: + j = 2 + x_third += np.gradient(x_second[i])[j:3] + return np.array([B_x] + x_first + x_second + x_third) + +# Solves a, b, c components with a[x, y, z], b[y, z, x] and c[z, x, y] up to given order +# Input B should be the output of solve_moments_from_B +def solve_coefficients(fg_b, x, y, z, order = 4): + B_moments = solve_moments_from_B(fg_b) + abc = np.zeros((3, 5, 4, 4)) + + # 4th order + if (order > 3): + print('4th order!') + for i in range(3): + j = (i+1) % 3 + k = (i+2) % 3 + coords = (x if i else x + 1, y if j else y + 1, z if k else z + 1) + print(coords) + a = abc[i] + Bx = B_moments[i] + + a[0][3][0] = 1/2 * (Bx[6, x, y, z] + Bx[6][coords]) + a[0][2][1] = 1/2 * (Bx[7, x, y, z] + Bx[7][coords]) + a[0][1][2] = 1/2 * (Bx[8, x, y, z] + Bx[8][coords]) + a[0][0][3] = 1/2 * (Bx[9, x, y, z] + Bx[9][coords]) + + a[1][3][0] = Bx[6, x, y, z] - Bx[6][coords] + a[1][2][1] = Bx[7, x, y, z] - Bx[7][coords] + a[1][1][2] = Bx[8, x, y, z] - Bx[8][coords] + a[1][0][3] = Bx[9, x, y, z] - Bx[9][coords] + + for i in range(3): + j = (i+1) % 3 + k = (i+2) % 3 + a = abc[i] + b = abc[j] + c = abc[k] + + # Should be correct, check later + a[4][0][0] = -1/4 * (b[1][0][3] + c[1][3][0]) + a[3][1][0] = -7/30 * c[1][2][1] + a[3][0][1] = -7/30 * b[1][1][2] + a[2][2][0] = -3/20 * c[1][1][2] + a[2][0][2] = -3/20 * b[1][2][1] + + # 3rd order + if (order > 2): + print('3rd order!') + for i in range(3): + j = (i+1) % 3 + k = (i+2) % 3 + coords = (x if i else x + 1, y if j else y + 1, z if k else z + 1) + a = abc[i] + Bx = B_moments[i] + + a[0][2][0] = 1/2 * (Bx[3, x, y, z] + Bx[3][coords]) - 1/6 * a[2][2][0] + a[0][1][1] = 1/2 * (Bx[4, x, y, z] + Bx[4][coords]) + a[0][0][2] = 1/2 * (Bx[5, x, y, z] + Bx[5][coords]) - 1/6 * a[2][0][2] + + a[1][2][0] = Bx[3, x, y, z] - Bx[3][coords] + a[1][1][1] = Bx[4, x, y, z] - Bx[4][coords] + a[1][0][2] = Bx[5, x, y, z] - Bx[5][coords] + + for i in range(3): + j = (i+1) % 3 + k = (i+2) % 3 + a = abc[i] + b = abc[j] + c = abc[k] + + # Should be correct, check later + a[3][0][0] = -1/3 * (b[1][0][2] + c[1][2][0]) + a[2][1][0] = -1/4 * c[1][1][1] + a[2][0][1] = -1/4 * b[1][1][1] + + # 2nd order + if (order > 1): + print('2nd order!') + for i in range(3): + j = (i+1) % 3 + k = (i+2) % 3 + coords = (x if i else x + 1, y if j else y + 1, z if k else z + 1) + a = abc[i] + Bx = B_moments[i] + + a[0][1][0] = 1/2 * (Bx[1, x, y, z] + Bx[1][coords]) - 1/6 * a[2][1][0] + a[0][0][1] = 1/2 * (Bx[2, x, y, z] + Bx[2][coords]) - 1/6 * a[2][0][1] + + a[1][1][0] = (Bx[1, x, y, z] - Bx[1][coords]) - 1/10 * a[3][1][0] + a[1][0][1] = (Bx[2, x, y, z] - Bx[2][coords]) - 1/10 * a[3][0][1] + + for i in range(3): + j = (i+1) % 3 + k = (i+2) % 3 + a = abc[i] + b = abc[j] + c = abc[k] + + # Should be correct, check later + a[2][0][0] = -1/2 * (b[1][0][1] + c[1][1][0]) - 3/35 * a[4][0][0] - 1/20 * (b[3][0][1] + c[3][1][0]) + + # 1st order + if (order > 0): + print('1st order!') + for i in range(3): + j = (i+1) % 3 + k = (i+2) % 3 + coords = (x if i else x + 1, y if j else y + 1, z if k else z + 1) + a = abc[i] + Bx = B_moments[i] + + a[1][0][0] = (Bx[0, x, y, z] - Bx[0][coords]) - 1/10 * a[3][0][0] + + # 0th order + for i in range(3): + print('0th order!') + j = (i+1) % 3 + k = (i+2) % 3 + coords = (x if i else x + 1, y if j else y + 1, z if k else z + 1) + a = abc[i] + Bx = B_moments[i] + + a[0][0][0] = 1/2 * (Bx[0, x, y, z] + Bx[0][coords]) - 1/6 * a[2][0][0] - 1/70 * a[4][0][0] + + # Check constraint: + test = abc[0][1][0][0] + abc[1][1][0][0] + abc[2][1][0][0] + 1/10 * (abc[0][3][0][0] + abc[1][3][0][0] + abc[2][3][0][0]) + if abs(test) > 0: + print("Something went wrong, sum (17) is " + str(test)) + + return abc + +# Code for testing +#f = pt.vlsvfile.VlsvReader('TEST/bulk.0000057.vlsv') +#fg_b = f.read_fsgrid_variable('fg_b') +#fg_b_vol = f.read_fsgrid_variable('fg_b_vol') +#for i in range(5): +# abc = solve_coefficients(fg_b, 5, 5, 5, i) +# print([interpolate_x(abc[0], 0.5, 0.5, 0.5), interpolate_x(abc[1], 0.5, 0.5, 0.5), interpolate_x(abc[2], 0.5, 0.5, 0.5)]) +#print(fg_b_vol[5, 5, 5]) \ No newline at end of file From 5de4545f646183a825ed5e91534636a44b27903e Mon Sep 17 00:00:00 2001 From: Leo Kotipalo Date: Thu, 11 Feb 2021 15:07:40 +0200 Subject: [PATCH 02/14] Fix coords used for solving coefficients, improve moment code --- scripts/divergenceless_reconstruction.py | 138 ++++++++++++++--------- 1 file changed, 84 insertions(+), 54 deletions(-) diff --git a/scripts/divergenceless_reconstruction.py b/scripts/divergenceless_reconstruction.py index 475a3952..7e96b875 100644 --- a/scripts/divergenceless_reconstruction.py +++ b/scripts/divergenceless_reconstruction.py @@ -1,6 +1,7 @@ import pytools as pt import numpy as np from scipy.special import legendre +from time import perf_counter P = tuple(map(lambda n: legendre(n), (0, 1, 2, 3, 4))) @@ -13,58 +14,70 @@ def interpolate_x(a, x, y, z): B_x += a[i][j][k] * P[i](x) * P[j](y) * P[k](z) return B_x +def interpolate_y(b, x, y, z): + return interpolate_x(b, y, z, x) + +def interpolate_z(c, x, y, z): + return interpolate_x(c, z, x, y) + # Returns moments as [ijk, order, x, y, z] def solve_moments_from_B(fg_b): - # Cyclically rotate each component and rotate coordinates back - x_moments = np.transpose(solve_x_moments(np.transpose(fg_b[:, :, :, 0], (0, 1, 2))), (0, 1, 2, 3)) - y_moments = np.transpose(solve_x_moments(np.transpose(fg_b[:, :, :, 1], (1, 2, 0))), (0, 3, 1, 2)) - z_moments = np.transpose(solve_x_moments(np.transpose(fg_b[:, :, :, 2], (2, 0, 1))), (0, 2, 3, 1)) + x_moments = solve_x_moments(fg_b[:, :, :, 0]) + y_moments = solve_y_moments(fg_b[:, :, :, 1]) + z_moments = solve_z_moments(fg_b[:, :, :, 2]) return (x_moments, y_moments, z_moments) # Returns moments for x as [order, x, y, z] # With order (0, y, z, yy, yz, zz, yyy, yyz, yzz, zzz) # y and z as cyclic rotations def solve_x_moments(B_x): - x_first = np.gradient(B_x)[1:3] - x_second = [] + x_moments = np.zeros((10,) + np.shape(B_x)) + x_moments[0] = B_x + x_moments[1:3] = np.gradient(B_x)[1:3] + start = 3 for i in range(2): j = i + 1 - x_second += np.gradient(x_first[i])[j:3] - x_third = [] + x_moments[start:start+3-j] = np.gradient(x_moments[1+i])[j:3] + start += (3-j) for i in range(3): - if i < 2: - j = i+1 - else: - j = 2 - x_third += np.gradient(x_second[i])[j:3] - return np.array([B_x] + x_first + x_second + x_third) + j = i + 1 if i < 2 else 2 + x_moments[start:start+3-j] = np.gradient(x_moments[3+i])[j:3] + start += (3-j) + return x_moments + +def solve_y_moments(B_y): + return np.transpose(solve_x_moments(np.transpose(B_y, (1, 2, 0))), (0, 3, 1, 2)) + +def solve_z_moments(B_z): + return np.transpose(solve_x_moments(np.transpose(B_z, (2, 0, 1))), (0, 2, 3, 1)) # Solves a, b, c components with a[x, y, z], b[y, z, x] and c[z, x, y] up to given order # Input B should be the output of solve_moments_from_B -def solve_coefficients(fg_b, x, y, z, order = 4): - B_moments = solve_moments_from_B(fg_b) +def solve_coefficients(B_moments, xyz, order = 4): abc = np.zeros((3, 5, 4, 4)) + x = xyz[0] + y = xyz[1] + z = xyz[2] # 4th order if (order > 3): - print('4th order!') for i in range(3): j = (i+1) % 3 k = (i+2) % 3 - coords = (x if i else x + 1, y if j else y + 1, z if k else z + 1) + coords = (x if i else x + 1, y if k else y + 1, z if j else z + 1) print(coords) a = abc[i] Bx = B_moments[i] - a[0][3][0] = 1/2 * (Bx[6, x, y, z] + Bx[6][coords]) - a[0][2][1] = 1/2 * (Bx[7, x, y, z] + Bx[7][coords]) - a[0][1][2] = 1/2 * (Bx[8, x, y, z] + Bx[8][coords]) - a[0][0][3] = 1/2 * (Bx[9, x, y, z] + Bx[9][coords]) + a[0][3][0] = 1/2 * (Bx[6][xyz] + Bx[6][coords]) + a[0][2][1] = 1/2 * (Bx[7][xyz] + Bx[7][coords]) + a[0][1][2] = 1/2 * (Bx[8][xyz] + Bx[8][coords]) + a[0][0][3] = 1/2 * (Bx[9][xyz] + Bx[9][coords]) - a[1][3][0] = Bx[6, x, y, z] - Bx[6][coords] - a[1][2][1] = Bx[7, x, y, z] - Bx[7][coords] - a[1][1][2] = Bx[8, x, y, z] - Bx[8][coords] - a[1][0][3] = Bx[9, x, y, z] - Bx[9][coords] + a[1][3][0] = Bx[6][xyz] - Bx[6][coords] + a[1][2][1] = Bx[7][xyz] - Bx[7][coords] + a[1][1][2] = Bx[8][xyz] - Bx[8][coords] + a[1][0][3] = Bx[9][xyz] - Bx[9][coords] for i in range(3): j = (i+1) % 3 @@ -82,21 +95,20 @@ def solve_coefficients(fg_b, x, y, z, order = 4): # 3rd order if (order > 2): - print('3rd order!') for i in range(3): j = (i+1) % 3 k = (i+2) % 3 - coords = (x if i else x + 1, y if j else y + 1, z if k else z + 1) + coords = (x if i else x + 1, y if k else y + 1, z if j else z + 1) a = abc[i] Bx = B_moments[i] - a[0][2][0] = 1/2 * (Bx[3, x, y, z] + Bx[3][coords]) - 1/6 * a[2][2][0] - a[0][1][1] = 1/2 * (Bx[4, x, y, z] + Bx[4][coords]) - a[0][0][2] = 1/2 * (Bx[5, x, y, z] + Bx[5][coords]) - 1/6 * a[2][0][2] + a[0][2][0] = 1/2 * (Bx[3][xyz] + Bx[3][coords]) - 1/6 * a[2][2][0] + a[0][1][1] = 1/2 * (Bx[4][xyz] + Bx[4][coords]) + a[0][0][2] = 1/2 * (Bx[5][xyz] + Bx[5][coords]) - 1/6 * a[2][0][2] - a[1][2][0] = Bx[3, x, y, z] - Bx[3][coords] - a[1][1][1] = Bx[4, x, y, z] - Bx[4][coords] - a[1][0][2] = Bx[5, x, y, z] - Bx[5][coords] + a[1][2][0] = Bx[3][xyz] - Bx[3][coords] + a[1][1][1] = Bx[4][xyz] - Bx[4][coords] + a[1][0][2] = Bx[5][xyz] - Bx[5][coords] for i in range(3): j = (i+1) % 3 @@ -112,19 +124,18 @@ def solve_coefficients(fg_b, x, y, z, order = 4): # 2nd order if (order > 1): - print('2nd order!') for i in range(3): j = (i+1) % 3 k = (i+2) % 3 - coords = (x if i else x + 1, y if j else y + 1, z if k else z + 1) + coords = (x if i else x + 1, y if k else y + 1, z if j else z + 1) a = abc[i] Bx = B_moments[i] - a[0][1][0] = 1/2 * (Bx[1, x, y, z] + Bx[1][coords]) - 1/6 * a[2][1][0] - a[0][0][1] = 1/2 * (Bx[2, x, y, z] + Bx[2][coords]) - 1/6 * a[2][0][1] + a[0][1][0] = 1/2 * (Bx[1][xyz] + Bx[1][coords]) - 1/6 * a[2][1][0] + a[0][0][1] = 1/2 * (Bx[2][xyz] + Bx[2][coords]) - 1/6 * a[2][0][1] - a[1][1][0] = (Bx[1, x, y, z] - Bx[1][coords]) - 1/10 * a[3][1][0] - a[1][0][1] = (Bx[2, x, y, z] - Bx[2][coords]) - 1/10 * a[3][0][1] + a[1][1][0] = (Bx[1][xyz] - Bx[1][coords]) - 1/10 * a[3][1][0] + a[1][0][1] = (Bx[2][xyz] - Bx[2][coords]) - 1/10 * a[3][0][1] for i in range(3): j = (i+1) % 3 @@ -138,22 +149,20 @@ def solve_coefficients(fg_b, x, y, z, order = 4): # 1st order if (order > 0): - print('1st order!') for i in range(3): j = (i+1) % 3 k = (i+2) % 3 - coords = (x if i else x + 1, y if j else y + 1, z if k else z + 1) + coords = (x if i else x + 1, y if k else y + 1, z if j else z + 1) a = abc[i] Bx = B_moments[i] - a[1][0][0] = (Bx[0, x, y, z] - Bx[0][coords]) - 1/10 * a[3][0][0] + a[1][0][0] = (Bx[0][xyz] - Bx[0][coords]) - 1/10 * a[3][0][0] # 0th order for i in range(3): - print('0th order!') j = (i+1) % 3 k = (i+2) % 3 - coords = (x if i else x + 1, y if j else y + 1, z if k else z + 1) + coords = (x if i else x + 1, y if k else y + 1, z if j else z + 1) a = abc[i] Bx = B_moments[i] @@ -161,16 +170,37 @@ def solve_coefficients(fg_b, x, y, z, order = 4): # Check constraint: test = abc[0][1][0][0] + abc[1][1][0][0] + abc[2][1][0][0] + 1/10 * (abc[0][3][0][0] + abc[1][3][0][0] + abc[2][3][0][0]) - if abs(test) > 0: - print("Something went wrong, sum (17) is " + str(test)) + #if abs(test) > 0: + # print("Something went wrong, sum (17) is " + str(test)) return abc +def center_value(B_moments, xyz, order=4): + abc = solve_coefficients(B_moments, xyz, order) + x = xyz[0] + y = xyz[1] + z = xyz[2] + return [interpolate_x(abc[0], 0.5, 0.5, 0.5), interpolate_y(abc[1], 0.5, 0.5, 0.5), interpolate_z(abc[2], 0.5, 0.5, 0.5)] + +def center_values(B_moments, coords, order=4): + return np.array(map(lambda xyz: center_value(B_moments, xyz, order), coords)) + +def all_center_values(B_moments, order=4): + return + # Code for testing -#f = pt.vlsvfile.VlsvReader('TEST/bulk.0000057.vlsv') -#fg_b = f.read_fsgrid_variable('fg_b') -#fg_b_vol = f.read_fsgrid_variable('fg_b_vol') -#for i in range(5): -# abc = solve_coefficients(fg_b, 5, 5, 5, i) -# print([interpolate_x(abc[0], 0.5, 0.5, 0.5), interpolate_x(abc[1], 0.5, 0.5, 0.5), interpolate_x(abc[2], 0.5, 0.5, 0.5)]) -#print(fg_b_vol[5, 5, 5]) \ No newline at end of file +f = pt.vlsvfile.VlsvReader('/wrk/users/lkotipal/TEST/bulk.0000057.vlsv') +fg_b = f.read_fsgrid_variable('fg_b') +print(np.shape(fg_b)) +fg_b_vol = f.read_fsgrid_variable('fg_b_vol') +print('File read!') +t = perf_counter() +B_moments = solve_moments_from_B(fg_b) +print(f'B_moments solved in {perf_counter() - t} seconds!') +for i in range(5): + t = perf_counter() + print(center_value(B_moments, (50, 50, 50), i)) + #abc = solve_coefficients(B_moments, 5, 5, 5, i) + #print([interpolate_x(abc[0], 0.5, 0.5, 0.5), interpolate_y(abc[1], 0.5, 0.5, 0.5), interpolate_z(abc[2], 0.5, 0.5, 0.5)]) + print(f'Coefficients up to order {i} solved in {perf_counter() - t} seconds!') +print(fg_b_vol[50, 50, 50]) \ No newline at end of file From a0f929640c5f49b79e376ddf9e884a8692e77ce6 Mon Sep 17 00:00:00 2001 From: Leo Kotipalo Date: Sun, 21 Feb 2021 03:11:42 +0200 Subject: [PATCH 03/14] Implement vectorization of solving coefficients. Still needs testing, might have discrepancies with single cell calculation. --- scripts/divergenceless_reconstruction.py | 139 ++++++++++++++++++++++- 1 file changed, 136 insertions(+), 3 deletions(-) diff --git a/scripts/divergenceless_reconstruction.py b/scripts/divergenceless_reconstruction.py index 7e96b875..f07ebb02 100644 --- a/scripts/divergenceless_reconstruction.py +++ b/scripts/divergenceless_reconstruction.py @@ -65,7 +65,6 @@ def solve_coefficients(B_moments, xyz, order = 4): j = (i+1) % 3 k = (i+2) % 3 coords = (x if i else x + 1, y if k else y + 1, z if j else z + 1) - print(coords) a = abc[i] Bx = B_moments[i] @@ -175,6 +174,133 @@ def solve_coefficients(B_moments, xyz, order = 4): return abc +def neighboursum(a, idx): + second = a[1:, 1:, 1:] + if idx == 0: + first = a[:-1, 1:, 1:] + elif idx == 1: + first = a[1:, :-1, 1:] + elif idx == 2: + first = a[1:, 1:, :-1] + return first + second + +def neighbourdiff(a, idx): + second = a[1:, 1:, 1:] + if idx == 0: + first = a[:-1, 1:, 1:] + elif idx == 1: + first = a[1:, :-1, 1:] + elif idx == 2: + first = a[1:, 1:, :-1] + return first - second + +# Solves a, b, c components with a[x, y, z], b[y, z, x] and c[z, x, y] up to given order +# Input B should be the output of solve_moments_from_B +def solve_all_coefficients(B_moments, order = 4): + shp = np.shape(B_moments[0][0]) + abc = np.zeros((3, 5, 4, 4, shp[0] - 1, shp[1] - 1, shp[2] - 1)) + + # 4th order + if (order > 3): + for i in range(3): + a = abc[i] + Bx = B_moments[i] + + a[0][3][0] = 1/2 * neighboursum(Bx[6], i) + a[0][2][1] = 1/2 * neighboursum(Bx[7], i) + a[0][1][2] = 1/2 * neighboursum(Bx[8], i) + a[0][0][3] = 1/2 * neighboursum(Bx[9], i) + + a[1][3][0] = neighbourdiff(Bx[6], i) + a[1][2][1] = neighbourdiff(Bx[7], i) + a[1][1][2] = neighbourdiff(Bx[8], i) + a[1][0][3] = neighbourdiff(Bx[9], i) + + for i in range(3): + j = (i+1) % 3 + k = (i+2) % 3 + a = abc[i] + b = abc[j] + c = abc[k] + + # Should be correct, check later + a[4][0][0] = -1/4 * (b[1][0][3] + c[1][3][0]) + a[3][1][0] = -7/30 * c[1][2][1] + a[3][0][1] = -7/30 * b[1][1][2] + a[2][2][0] = -3/20 * c[1][1][2] + a[2][0][2] = -3/20 * b[1][2][1] + + # 3rd order + if (order > 2): + for i in range(3): + a = abc[i] + Bx = B_moments[i] + + a[0][2][0] = 1/2 * neighboursum(Bx[3], i) - 1/6 * a[2][2][0] + a[0][1][1] = 1/2 * neighboursum(Bx[4], i) + a[0][0][2] = 1/2 * neighboursum(Bx[5], i) - 1/6 * a[2][0][2] + + a[1][2][0] = neighbourdiff(Bx[3], i) + a[1][1][1] = neighbourdiff(Bx[4], i) + a[1][0][2] = neighbourdiff(Bx[5], i) + + for i in range(3): + j = (i+1) % 3 + k = (i+2) % 3 + a = abc[i] + b = abc[j] + c = abc[k] + + # Should be correct, check later + a[3][0][0] = -1/3 * (b[1][0][2] + c[1][2][0]) + a[2][1][0] = -1/4 * c[1][1][1] + a[2][0][1] = -1/4 * b[1][1][1] + + # 2nd order + if (order > 1): + for i in range(3): + a = abc[i] + Bx = B_moments[i] + + a[0][1][0] = 1/2 * neighboursum(Bx[1], i) - 1/6 * a[2][1][0] + a[0][0][1] = 1/2 * neighboursum(Bx[2], i) - 1/6 * a[2][0][1] + + a[1][1][0] = neighbourdiff(Bx[1], i) - 1/10 * a[3][1][0] + a[1][0][1] = neighbourdiff(Bx[2], i) - 1/10 * a[3][0][1] + + for i in range(3): + j = (i+1) % 3 + k = (i+2) % 3 + a = abc[i] + b = abc[j] + c = abc[k] + + # Should be correct, check later + a[2][0][0] = -1/2 * (b[1][0][1] + c[1][1][0]) - 3/35 * a[4][0][0] - 1/20 * (b[3][0][1] + c[3][1][0]) + + # 1st order + if (order > 0): + for i in range(3): + a = abc[i] + Bx = B_moments[i] + + a[1][0][0] = neighbourdiff(Bx[0], i) - 1/10 * a[3][0][0] + + # 0th order + for i in range(3): + a = abc[i] + Bx = B_moments[i] + + a[0][0][0] = 1/2 * neighboursum(Bx[0], i) - 1/6 * a[2][0][0] - 1/70 * a[4][0][0] + + # Check constraint: + test = abc[0][1][0][0] + abc[1][1][0][0] + abc[2][1][0][0] + 1/10 * (abc[0][3][0][0] + abc[1][3][0][0] + abc[2][3][0][0]) + print(np.amax(test)) + #if abs(test) > 0: + # print("Something went wrong, sum (17) is " + str(test)) + + return np.pad(abc, [(0, 0), (0, 0), (0, 0), (0, 0), (0, 1), (0, 1), (0, 1)]) + def center_value(B_moments, xyz, order=4): abc = solve_coefficients(B_moments, xyz, order) x = xyz[0] @@ -195,12 +321,19 @@ def all_center_values(B_moments, order=4): fg_b_vol = f.read_fsgrid_variable('fg_b_vol') print('File read!') t = perf_counter() -B_moments = solve_moments_from_B(fg_b) +B_moments = solve_moments_from_B(fg_b[0:100, 0:100, 0:100]) print(f'B_moments solved in {perf_counter() - t} seconds!') for i in range(5): t = perf_counter() + coeffs = solve_coefficients(B_moments, (50, 50, 50), i) + print(f'Single coefficients up to order {i} solved in {perf_counter() - t} seconds!') print(center_value(B_moments, (50, 50, 50), i)) + + t = perf_counter() + all_coeffs = solve_all_coefficients(B_moments, i) + print(f'All coefficients up to order {i} solved in {perf_counter() - t} seconds!') + + print(all_coeffs[:, :, :, :, 50, 50, 50] - coeffs) #abc = solve_coefficients(B_moments, 5, 5, 5, i) #print([interpolate_x(abc[0], 0.5, 0.5, 0.5), interpolate_y(abc[1], 0.5, 0.5, 0.5), interpolate_z(abc[2], 0.5, 0.5, 0.5)]) - print(f'Coefficients up to order {i} solved in {perf_counter() - t} seconds!') print(fg_b_vol[50, 50, 50]) \ No newline at end of file From 25c16894fe661cac3e29ed145f3ec813c98a3bc9 Mon Sep 17 00:00:00 2001 From: Leo Kotipalo Date: Fri, 5 Mar 2021 14:59:29 +0200 Subject: [PATCH 04/14] Fix neighboursum and neighbourdiff --- scripts/divergenceless_reconstruction.py | 40 +++++++++++++++--------- 1 file changed, 25 insertions(+), 15 deletions(-) diff --git a/scripts/divergenceless_reconstruction.py b/scripts/divergenceless_reconstruction.py index f07ebb02..355caa46 100644 --- a/scripts/divergenceless_reconstruction.py +++ b/scripts/divergenceless_reconstruction.py @@ -175,23 +175,23 @@ def solve_coefficients(B_moments, xyz, order = 4): return abc def neighboursum(a, idx): - second = a[1:, 1:, 1:] + second = a[:-1, :-1, :-1] if idx == 0: - first = a[:-1, 1:, 1:] + first = a[1:, :-1, :-1] elif idx == 1: - first = a[1:, :-1, 1:] + first = a[:-1, 1:, :-1] elif idx == 2: - first = a[1:, 1:, :-1] + first = a[:-1, :-1, 1:] return first + second def neighbourdiff(a, idx): - second = a[1:, 1:, 1:] + second = a[:-1, :-1, :-1] if idx == 0: - first = a[:-1, 1:, 1:] + first = a[1:, :-1, :-1] elif idx == 1: - first = a[1:, :-1, 1:] + first = a[:-1, 1:, :-1] elif idx == 2: - first = a[1:, 1:, :-1] + first = a[:-1, :-1, 1:] return first - second # Solves a, b, c components with a[x, y, z], b[y, z, x] and c[z, x, y] up to given order @@ -303,10 +303,18 @@ def solve_all_coefficients(B_moments, order = 4): def center_value(B_moments, xyz, order=4): abc = solve_coefficients(B_moments, xyz, order) - x = xyz[0] - y = xyz[1] - z = xyz[2] - return [interpolate_x(abc[0], 0.5, 0.5, 0.5), interpolate_y(abc[1], 0.5, 0.5, 0.5), interpolate_z(abc[2], 0.5, 0.5, 0.5)] + return [interpolate_x(abc[0], 0, 0, 0), interpolate_y(abc[1], 0, 0, 0), interpolate_z(abc[2], 0, 0, 0)] + +def ave_B(B_moments, xyz, order=4): + abc = solve_coefficients(B_moments, xyz, order) + a = abc[0] + b = abc[1] + c = abc[2] + return ( + a[0][0][0] - 3/8 * (a[2][0][0] + a[0][2][0] + a[0][0][2]) + 9/64 * (a[2][2][0] + a[2][0][2]) + 15/128 * a[4][0][0], + b[0][0][0] - 3/8 * (b[2][0][0] + b[0][2][0] + b[0][0][2]) + 9/64 * (b[2][2][0] + b[2][0][2]) + 15/128 * b[4][0][0], + c[0][0][0] - 3/8 * (c[2][0][0] + c[0][2][0] + c[0][0][2]) + 9/64 * (c[2][2][0] + c[2][0][2]) + 15/128 * c[4][0][0] + ) def center_values(B_moments, coords, order=4): return np.array(map(lambda xyz: center_value(B_moments, xyz, order), coords)) @@ -326,14 +334,16 @@ def all_center_values(B_moments, order=4): for i in range(5): t = perf_counter() coeffs = solve_coefficients(B_moments, (50, 50, 50), i) + #print(ave_B(B_moments, (50, 50, 50), i)) print(f'Single coefficients up to order {i} solved in {perf_counter() - t} seconds!') - print(center_value(B_moments, (50, 50, 50), i)) + print([interpolate_x(coeffs[0], 0, 0, 0), interpolate_y(coeffs[1], 0, 0, 0), interpolate_z(coeffs[2], 0, 0, 0)]) t = perf_counter() all_coeffs = solve_all_coefficients(B_moments, i) print(f'All coefficients up to order {i} solved in {perf_counter() - t} seconds!') + print([interpolate_x(all_coeffs[0,:,:,:, 50, 50, 50], 0, 0, 0), interpolate_y(all_coeffs[1,:,:,:, 50, 50, 50], 0, 0, 0), interpolate_z(all_coeffs[2,:,:,:, 50, 50, 50], 0, 0, 0)]) - print(all_coeffs[:, :, :, :, 50, 50, 50] - coeffs) + #print(all_coeffs[:, :, :, :, 50, 50, 50] - coeffs) #abc = solve_coefficients(B_moments, 5, 5, 5, i) #print([interpolate_x(abc[0], 0.5, 0.5, 0.5), interpolate_y(abc[1], 0.5, 0.5, 0.5), interpolate_z(abc[2], 0.5, 0.5, 0.5)]) -print(fg_b_vol[50, 50, 50]) \ No newline at end of file +print(fg_b_vol[50, 50, 50]) From 1576e7b81e602c5cc16491e095d20ba7ca74c46b Mon Sep 17 00:00:00 2001 From: Leo Kotipalo Date: Fri, 5 Mar 2021 17:20:45 +0200 Subject: [PATCH 05/14] Split tests away, implement more vectorization --- scripts/divergenceless_reconstruction.py | 55 ++++++++----------- scripts/divergenceless_reconstruction_test.py | 35 ++++++++++++ 2 files changed, 59 insertions(+), 31 deletions(-) create mode 100644 scripts/divergenceless_reconstruction_test.py diff --git a/scripts/divergenceless_reconstruction.py b/scripts/divergenceless_reconstruction.py index 355caa46..a12475fa 100644 --- a/scripts/divergenceless_reconstruction.py +++ b/scripts/divergenceless_reconstruction.py @@ -1,7 +1,6 @@ import pytools as pt import numpy as np from scipy.special import legendre -from time import perf_counter P = tuple(map(lambda n: legendre(n), (0, 1, 2, 3, 4))) @@ -53,6 +52,7 @@ def solve_z_moments(B_z): # Solves a, b, c components with a[x, y, z], b[y, z, x] and c[z, x, y] up to given order # Input B should be the output of solve_moments_from_B +# Might be deprecated. A lot faster than calculating all coefficients but that's only a few seconds anyway for a 100^3 array def solve_coefficients(B_moments, xyz, order = 4): abc = np.zeros((3, 5, 4, 4)) x = xyz[0] @@ -257,6 +257,7 @@ def solve_all_coefficients(B_moments, order = 4): a[2][0][1] = -1/4 * b[1][1][1] # 2nd order + if (order > 1): for i in range(3): a = abc[i] @@ -294,11 +295,12 @@ def solve_all_coefficients(B_moments, order = 4): a[0][0][0] = 1/2 * neighboursum(Bx[0], i) - 1/6 * a[2][0][0] - 1/70 * a[4][0][0] # Check constraint: - test = abc[0][1][0][0] + abc[1][1][0][0] + abc[2][1][0][0] + 1/10 * (abc[0][3][0][0] + abc[1][3][0][0] + abc[2][3][0][0]) - print(np.amax(test)) + #test = abc[0][1][0][0] + abc[1][1][0][0] + abc[2][1][0][0] + 1/10 * (abc[0][3][0][0] + abc[1][3][0][0] + abc[2][3][0][0]) + #print(np.amax(test)) #if abs(test) > 0: # print("Something went wrong, sum (17) is " + str(test)) + # return np.pad(np.transpose(abc, (4, 5, 6, 0, 1, 2, 3)), [(0, 1), (0, 1), (0, 1), (0, 0), (0, 0), (0, 0), (0, 0)]) return np.pad(abc, [(0, 0), (0, 0), (0, 0), (0, 0), (0, 1), (0, 1), (0, 1)]) def center_value(B_moments, xyz, order=4): @@ -316,34 +318,25 @@ def ave_B(B_moments, xyz, order=4): c[0][0][0] - 3/8 * (c[2][0][0] + c[0][2][0] + c[0][0][2]) + 9/64 * (c[2][2][0] + c[2][0][2]) + 15/128 * c[4][0][0] ) +def all_ave_Bs(B_moments, order=4): + abc = solve_all_coefficients(B_moments, order) + a = abc[0] + b = abc[1] + c = abc[2] + return np.transpose(np.array( + ( + a[0][0][0] - 3/8 * (a[2][0][0] + a[0][2][0] + a[0][0][2]) + 9/64 * (a[2][2][0] + a[2][0][2]) + 15/128 * a[4][0][0], + b[0][0][0] - 3/8 * (b[2][0][0] + b[0][2][0] + b[0][0][2]) + 9/64 * (b[2][2][0] + b[2][0][2]) + 15/128 * b[4][0][0], + c[0][0][0] - 3/8 * (c[2][0][0] + c[0][2][0] + c[0][0][2]) + 9/64 * (c[2][2][0] + c[2][0][2]) + 15/128 * c[4][0][0] + ) + ), [1, 2, 3, 0]) + def center_values(B_moments, coords, order=4): - return np.array(map(lambda xyz: center_value(B_moments, xyz, order), coords)) + # Looks scuffed but is faster + abc = np.transpose(np.transpose(solve_all_coefficients(B_moments, order), (4, 5, 6, 0, 1, 2, 3))[coords[:, 0], coords[:, 1], coords[:, 2]], (1, 2, 3, 4, 0)) + return np.transpose(np.array([interpolate_x(abc[0], 0, 0, 0), interpolate_y(abc[1], 0, 0, 0), interpolate_z(abc[2], 0, 0, 0)])) + #return all_center_values(B_moments, order)[coords[:, 0], coords[:, 1], coords[:, 2], :] def all_center_values(B_moments, order=4): - return - -# Code for testing -f = pt.vlsvfile.VlsvReader('/wrk/users/lkotipal/TEST/bulk.0000057.vlsv') -fg_b = f.read_fsgrid_variable('fg_b') -print(np.shape(fg_b)) -fg_b_vol = f.read_fsgrid_variable('fg_b_vol') -print('File read!') -t = perf_counter() -B_moments = solve_moments_from_B(fg_b[0:100, 0:100, 0:100]) -print(f'B_moments solved in {perf_counter() - t} seconds!') -for i in range(5): - t = perf_counter() - coeffs = solve_coefficients(B_moments, (50, 50, 50), i) - #print(ave_B(B_moments, (50, 50, 50), i)) - print(f'Single coefficients up to order {i} solved in {perf_counter() - t} seconds!') - print([interpolate_x(coeffs[0], 0, 0, 0), interpolate_y(coeffs[1], 0, 0, 0), interpolate_z(coeffs[2], 0, 0, 0)]) - - t = perf_counter() - all_coeffs = solve_all_coefficients(B_moments, i) - print(f'All coefficients up to order {i} solved in {perf_counter() - t} seconds!') - print([interpolate_x(all_coeffs[0,:,:,:, 50, 50, 50], 0, 0, 0), interpolate_y(all_coeffs[1,:,:,:, 50, 50, 50], 0, 0, 0), interpolate_z(all_coeffs[2,:,:,:, 50, 50, 50], 0, 0, 0)]) - - #print(all_coeffs[:, :, :, :, 50, 50, 50] - coeffs) - #abc = solve_coefficients(B_moments, 5, 5, 5, i) - #print([interpolate_x(abc[0], 0.5, 0.5, 0.5), interpolate_y(abc[1], 0.5, 0.5, 0.5), interpolate_z(abc[2], 0.5, 0.5, 0.5)]) -print(fg_b_vol[50, 50, 50]) + abc = solve_all_coefficients(B_moments, order) + return np.transpose(np.array([interpolate_x(abc[0], 0, 0, 0), interpolate_y(abc[1], 0, 0, 0), interpolate_z(abc[2], 0, 0, 0)]), (1, 2, 3, 0)) \ No newline at end of file diff --git a/scripts/divergenceless_reconstruction_test.py b/scripts/divergenceless_reconstruction_test.py new file mode 100644 index 00000000..c3ad71d6 --- /dev/null +++ b/scripts/divergenceless_reconstruction_test.py @@ -0,0 +1,35 @@ +from divergenceless_reconstruction import solve_moments_from_B, center_value, center_values, all_center_values, all_ave_Bs +import pytools as pt +import numpy as np +from time import perf_counter +# Code for testing divergenceless_resconstruction.py + +f = pt.vlsvfile.VlsvReader('/wrk/users/lkotipal/TEST/bulk.0000057.vlsv') +fg_b = f.read_fsgrid_variable('fg_b') +print(np.shape(fg_b)) +fg_b_vol = f.read_fsgrid_variable('fg_b_vol') +print('File read!') +t = perf_counter() +B_moments = solve_moments_from_B(fg_b[0:100, 0:100, 0:100]) +print(f'B_moments solved in {perf_counter() - t} seconds!') +print() +print(f'Volumetric b in (50, 50, 50): {fg_b_vol[50, 50, 50]}') +print() +for i in range(5): + t = perf_counter() + print(center_value(B_moments, (50, 50, 50), i)) + print(f'Single coefficients up to order {i} solved in {perf_counter() - t} seconds!') + + t = perf_counter() + print(np.shape(center_values(B_moments, np.array([(50, 50, 50), (51, 51, 51), (54, 54, 64), (43, 65, 42), (43, 65, 32)]), i))) + print(f'Multiple coefficients up to order {i} solved in {perf_counter() - t} seconds!') + + t = perf_counter() + print(all_center_values(B_moments, i)[50, 50, 50]) + print(f'All coefficients up to order {i} solved in {perf_counter() - t} seconds!') + + print(f'Average scaled difference between Analysator and Vlasiator in order {i}:') + print(np.average(((all_ave_Bs(B_moments, i))[:-1, :-1, :-1] - fg_b_vol[:99, :99, :99])/fg_b_vol[:99,:99,:99])) + print(f'Max scaled difference between Analysator and Vlasiator in order {i}:') + print(np.max(abs(((all_ave_Bs(B_moments, i))[:-1, :-1, :-1] - fg_b_vol[:99, :99, :99])/fg_b_vol[:99,:99,:99]))) + print() \ No newline at end of file From 4fa4c89562797b6af4fca482f9cdc54c1db95a2d Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Mon, 19 Jun 2023 16:44:40 +0300 Subject: [PATCH 06/14] Diving for Jacobians. --- scripts/divergenceless_reconstruction.py | 58 +++++++++++++++++++++++- 1 file changed, 57 insertions(+), 1 deletion(-) diff --git a/scripts/divergenceless_reconstruction.py b/scripts/divergenceless_reconstruction.py index a12475fa..7034b195 100644 --- a/scripts/divergenceless_reconstruction.py +++ b/scripts/divergenceless_reconstruction.py @@ -19,6 +19,50 @@ def interpolate_y(b, x, y, z): def interpolate_z(c, x, y, z): return interpolate_x(c, z, x, y) +# y and z are simply cyclic rotations of this +def interpolate_dbxdx(a, x, y, z): + B_x = 0 + for i in range(5): + for j in range(4): + for k in range(4): + B_x += a[i][j][k] * np.polyder(P[i])(x) * P[j](y) * P[k](z) + return B_x + +def interpolate_dbydx(b, x, y, z): + return interpolate_dbxdx(b, y, z, x) + +def interpolate_dbzdx(c, x, y, z): + return interpolate_dbxdx(c, z, x, y) + +def interpolate_dbxdy(a, x, y, z): + B_x = 0 + for i in range(5): + for j in range(4): + for k in range(4): + B_x += a[i][j][k] * P[i](x) * np.polyder(P[j])(y) * P[k](z) + return B_x + +def interpolate_dbydy(b, x, y, z): + return interpolate_dbxdy(b, y, z, x) + +def interpolate_dbzdy(c, x, y, z): + return interpolate_dbxdy(c, z, x, y) + +def interpolate_dbxdz(a, x, y, z): + B_x = 0 + for i in range(5): + for j in range(4): + for k in range(4): + B_x += a[i][j][k] * P[i](x) * P[j](y) * np.polyder(P[k])(z) + return B_x + +def interpolate_dbydz(b, x, y, z): + return interpolate_dbxdz(b, y, z, x) + +def interpolate_dbzdz(c, x, y, z): + return interpolate_dbxdz(c, z, x, y) + + # Returns moments as [ijk, order, x, y, z] def solve_moments_from_B(fg_b): x_moments = solve_x_moments(fg_b[:, :, :, 0]) @@ -339,4 +383,16 @@ def center_values(B_moments, coords, order=4): def all_center_values(B_moments, order=4): abc = solve_all_coefficients(B_moments, order) - return np.transpose(np.array([interpolate_x(abc[0], 0, 0, 0), interpolate_y(abc[1], 0, 0, 0), interpolate_z(abc[2], 0, 0, 0)]), (1, 2, 3, 0)) \ No newline at end of file + return np.transpose(np.array([interpolate_x(abc[0], 0, 0, 0), interpolate_y(abc[1], 0, 0, 0), interpolate_z(abc[2], 0, 0, 0)]), (1, 2, 3, 0)) + +def all_center_values_dbxdi(B_moments, order=4): + abc = solve_all_coefficients(B_moments, order) + return np.transpose(np.array([interpolate_dbxdx(abc[0], 0, 0, 0), interpolate_dbxdy(abc[1], 0, 0, 0), interpolate_dbxdz(abc[2], 0, 0, 0)]), (1, 2, 3, 0)) + +def all_center_values_dbydi(B_moments, order=4): + abc = solve_all_coefficients(B_moments, order) + return np.transpose(np.array([interpolate_dbydx(abc[0], 0, 0, 0), interpolate_dbydy(abc[1], 0, 0, 0), interpolate_dbydz(abc[2], 0, 0, 0)]), (1, 2, 3, 0)) + +def all_center_values_dbzdi(B_moments, order=4): + abc = solve_all_coefficients(B_moments, order) + return np.transpose(np.array([interpolate_dbzdx(abc[0], 0, 0, 0), interpolate_dbzdy(abc[1], 0, 0, 0), interpolate_dbzdz(abc[2], 0, 0, 0)]), (1, 2, 3, 0)) \ No newline at end of file From ba6b371bc96d1fcb96ec2138d4b11860eaabfd7c Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Tue, 20 Jun 2023 11:10:22 +0300 Subject: [PATCH 07/14] Divergence-free derivatives evaluation included. --- scripts/divergenceless_reconstruction.py | 57 +++++++++++++++++++----- 1 file changed, 47 insertions(+), 10 deletions(-) diff --git a/scripts/divergenceless_reconstruction.py b/scripts/divergenceless_reconstruction.py index 7034b195..fc57b6a5 100644 --- a/scripts/divergenceless_reconstruction.py +++ b/scripts/divergenceless_reconstruction.py @@ -29,10 +29,22 @@ def interpolate_dbxdx(a, x, y, z): return B_x def interpolate_dbydx(b, x, y, z): - return interpolate_dbxdx(b, y, z, x) + B_y = 0 + for i in range(5): + for j in range(4): + for k in range(4): + B_y += b[i][j][k] * P[i](y) * P[j](z) * np.polyder(P[k])(x) + return B_y + #return interpolate_dbxdx(b, y, z, x) def interpolate_dbzdx(c, x, y, z): - return interpolate_dbxdx(c, z, x, y) + B_z = 0 + for i in range(5): + for j in range(4): + for k in range(4): + B_z += c[i][j][k] * P[i](z) * np.polyder(P[j])(x) * P[k](y) + return B_z + #return interpolate_dbxdx(c, z, x, y) def interpolate_dbxdy(a, x, y, z): B_x = 0 @@ -43,10 +55,22 @@ def interpolate_dbxdy(a, x, y, z): return B_x def interpolate_dbydy(b, x, y, z): - return interpolate_dbxdy(b, y, z, x) + B_y = 0 + for i in range(5): + for j in range(4): + for k in range(4): + B_y += b[i][j][k] * np.polyder(P[i])(y) * P[j](z) * P[k](x) + return B_y + #return interpolate_dbxdy(b, y, z, x) def interpolate_dbzdy(c, x, y, z): - return interpolate_dbxdy(c, z, x, y) + B_z = 0 + for i in range(5): + for j in range(4): + for k in range(4): + B_z += c[i][j][k] * P[i](z) * P[j](x) * np.polyder(P[k])(y) + return B_z + # return interpolate_dbxdy(c, z, x, y) def interpolate_dbxdz(a, x, y, z): B_x = 0 @@ -57,10 +81,22 @@ def interpolate_dbxdz(a, x, y, z): return B_x def interpolate_dbydz(b, x, y, z): - return interpolate_dbxdz(b, y, z, x) + B_y = 0 + for i in range(5): + for j in range(4): + for k in range(4): + B_y += b[i][j][k] * P[i](y) * np.polyder(P[j])(z) * P[k](x) + return B_y + # return interpolate_dbxdz(b, y, z, x) def interpolate_dbzdz(c, x, y, z): - return interpolate_dbxdz(c, z, x, y) + B_z = 0 + for i in range(5): + for j in range(4): + for k in range(4): + B_z += c[i][j][k] * np.polyder(P[i])(z) * P[j](x) * P[k](y) + return B_z + # return interpolate_dbxdz(c, z, x, y) # Returns moments as [ijk, order, x, y, z] @@ -383,16 +419,17 @@ def center_values(B_moments, coords, order=4): def all_center_values(B_moments, order=4): abc = solve_all_coefficients(B_moments, order) - return np.transpose(np.array([interpolate_x(abc[0], 0, 0, 0), interpolate_y(abc[1], 0, 0, 0), interpolate_z(abc[2], 0, 0, 0)]), (1, 2, 3, 0)) + return np.transpose(np.array([interpolate_x(abc[0], 0, 0, 0), interpolate_y(abc[1], 0, 0, 0), interpolate_z(abc[2], 0, 0, 0)]), (1, 2, 3, 0)) def all_center_values_dbxdi(B_moments, order=4): abc = solve_all_coefficients(B_moments, order) - return np.transpose(np.array([interpolate_dbxdx(abc[0], 0, 0, 0), interpolate_dbxdy(abc[1], 0, 0, 0), interpolate_dbxdz(abc[2], 0, 0, 0)]), (1, 2, 3, 0)) + return np.transpose(np.array([interpolate_dbxdx(abc[0], 0, 0, 0), interpolate_dbxdy(abc[0], 0, 0, 0), interpolate_dbxdz(abc[0], 0, 0, 0)]), (1, 2, 3, 0)) def all_center_values_dbydi(B_moments, order=4): abc = solve_all_coefficients(B_moments, order) - return np.transpose(np.array([interpolate_dbydx(abc[0], 0, 0, 0), interpolate_dbydy(abc[1], 0, 0, 0), interpolate_dbydz(abc[2], 0, 0, 0)]), (1, 2, 3, 0)) + return np.transpose(np.array([interpolate_dbydx(abc[1], 0, 0, 0), interpolate_dbydy(abc[1], 0, 0, 0), interpolate_dbydz(abc[1], 0, 0, 0)]), (1, 2, 3, 0)) def all_center_values_dbzdi(B_moments, order=4): abc = solve_all_coefficients(B_moments, order) - return np.transpose(np.array([interpolate_dbzdx(abc[0], 0, 0, 0), interpolate_dbzdy(abc[1], 0, 0, 0), interpolate_dbzdz(abc[2], 0, 0, 0)]), (1, 2, 3, 0)) \ No newline at end of file + return np.transpose(np.array([interpolate_dbzdx(abc[2], 0, 0, 0), interpolate_dbzdy(abc[2], 0, 0, 0), interpolate_dbzdz(abc[2], 0, 0, 0)]), (1, 2, 3, 0)) + From 76d092fee7a7c459efaff8cea44742945c8e06b4 Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Wed, 21 Jun 2023 08:57:43 +0300 Subject: [PATCH 08/14] Variable order reconstruction with reduced memory footprint. Performance for EGI fsgrid ops ~quadrupled. Still to encapsulate this into a reconstructor class. --- scripts/divergenceless_reconstruction.py | 205 +++++++++++++++-------- 1 file changed, 132 insertions(+), 73 deletions(-) diff --git a/scripts/divergenceless_reconstruction.py b/scripts/divergenceless_reconstruction.py index fc57b6a5..f72388ac 100644 --- a/scripts/divergenceless_reconstruction.py +++ b/scripts/divergenceless_reconstruction.py @@ -3,13 +3,13 @@ from scipy.special import legendre P = tuple(map(lambda n: legendre(n), (0, 1, 2, 3, 4))) - +order = 2 # y and z are simply cyclic rotations of this def interpolate_x(a, x, y, z): B_x = 0 - for i in range(5): - for j in range(4): - for k in range(4): + for i in range(order+1): + for j in range(order): + for k in range(order): B_x += a[i][j][k] * P[i](x) * P[j](y) * P[k](z) return B_x @@ -22,78 +22,78 @@ def interpolate_z(c, x, y, z): # y and z are simply cyclic rotations of this def interpolate_dbxdx(a, x, y, z): B_x = 0 - for i in range(5): - for j in range(4): - for k in range(4): + for i in range(order+1): + for j in range(order): + for k in range(order): B_x += a[i][j][k] * np.polyder(P[i])(x) * P[j](y) * P[k](z) return B_x def interpolate_dbydx(b, x, y, z): B_y = 0 - for i in range(5): - for j in range(4): - for k in range(4): + for i in range(order+1): + for j in range(order): + for k in range(order): B_y += b[i][j][k] * P[i](y) * P[j](z) * np.polyder(P[k])(x) return B_y #return interpolate_dbxdx(b, y, z, x) def interpolate_dbzdx(c, x, y, z): B_z = 0 - for i in range(5): - for j in range(4): - for k in range(4): + for i in range(order+1): + for j in range(order): + for k in range(order): B_z += c[i][j][k] * P[i](z) * np.polyder(P[j])(x) * P[k](y) return B_z #return interpolate_dbxdx(c, z, x, y) def interpolate_dbxdy(a, x, y, z): B_x = 0 - for i in range(5): - for j in range(4): - for k in range(4): + for i in range(order+1): + for j in range(order): + for k in range(order): B_x += a[i][j][k] * P[i](x) * np.polyder(P[j])(y) * P[k](z) return B_x def interpolate_dbydy(b, x, y, z): B_y = 0 - for i in range(5): - for j in range(4): - for k in range(4): + for i in range(order+1): + for j in range(order): + for k in range(order): B_y += b[i][j][k] * np.polyder(P[i])(y) * P[j](z) * P[k](x) return B_y #return interpolate_dbxdy(b, y, z, x) def interpolate_dbzdy(c, x, y, z): B_z = 0 - for i in range(5): - for j in range(4): - for k in range(4): + for i in range(order+1): + for j in range(order): + for k in range(order): B_z += c[i][j][k] * P[i](z) * P[j](x) * np.polyder(P[k])(y) return B_z # return interpolate_dbxdy(c, z, x, y) def interpolate_dbxdz(a, x, y, z): B_x = 0 - for i in range(5): - for j in range(4): - for k in range(4): + for i in range(order+1): + for j in range(order): + for k in range(order): B_x += a[i][j][k] * P[i](x) * P[j](y) * np.polyder(P[k])(z) return B_x def interpolate_dbydz(b, x, y, z): B_y = 0 - for i in range(5): - for j in range(4): - for k in range(4): + for i in range(order+1): + for j in range(order): + for k in range(order): B_y += b[i][j][k] * P[i](y) * np.polyder(P[j])(z) * P[k](x) return B_y # return interpolate_dbxdz(b, y, z, x) def interpolate_dbzdz(c, x, y, z): B_z = 0 - for i in range(5): - for j in range(4): - for k in range(4): + for i in range(order+1): + for j in range(order): + for k in range(order): B_z += c[i][j][k] * np.polyder(P[i])(z) * P[j](x) * P[k](y) return B_z # return interpolate_dbxdz(c, z, x, y) @@ -133,8 +133,8 @@ def solve_z_moments(B_z): # Solves a, b, c components with a[x, y, z], b[y, z, x] and c[z, x, y] up to given order # Input B should be the output of solve_moments_from_B # Might be deprecated. A lot faster than calculating all coefficients but that's only a few seconds anyway for a 100^3 array -def solve_coefficients(B_moments, xyz, order = 4): - abc = np.zeros((3, 5, 4, 4)) +def solve_coefficients(B_moments, xyz): + abc = np.zeros((3, order+1, order, order)) x = xyz[0] y = xyz[1] z = xyz[2] @@ -213,8 +213,12 @@ def solve_coefficients(B_moments, xyz, order = 4): a[0][1][0] = 1/2 * (Bx[1][xyz] + Bx[1][coords]) - 1/6 * a[2][1][0] a[0][0][1] = 1/2 * (Bx[2][xyz] + Bx[2][coords]) - 1/6 * a[2][0][1] - a[1][1][0] = (Bx[1][xyz] - Bx[1][coords]) - 1/10 * a[3][1][0] - a[1][0][1] = (Bx[2][xyz] - Bx[2][coords]) - 1/10 * a[3][0][1] + a[1][1][0] = (Bx[1][xyz] - Bx[1][coords]) + if order > 2: + a[1][1][0] -= 1/10 * a[3][1][0] + a[1][0][1] = (Bx[2][xyz] - Bx[2][coords]) + if order > 2: + a[1][1][0] -= 1/10 * a[3][0][1] for i in range(3): j = (i+1) % 3 @@ -224,7 +228,11 @@ def solve_coefficients(B_moments, xyz, order = 4): c = abc[k] # Should be correct, check later - a[2][0][0] = -1/2 * (b[1][0][1] + c[1][1][0]) - 3/35 * a[4][0][0] - 1/20 * (b[3][0][1] + c[3][1][0]) + a[2][0][0] = -1/2 * (b[1][0][1] + c[1][1][0]) + if order > 3: + a[2][0][0] -= 3/35 * a[4][0][0] + if order > 2: + a[2][0][0] -= 1/20 * (b[3][0][1] + c[3][1][0]) # 1st order if (order > 0): @@ -235,7 +243,9 @@ def solve_coefficients(B_moments, xyz, order = 4): a = abc[i] Bx = B_moments[i] - a[1][0][0] = (Bx[0][xyz] - Bx[0][coords]) - 1/10 * a[3][0][0] + a[1][0][0] = (Bx[0][xyz] - Bx[0][coords]) + if order > 2: + a[1][0][0] -= 1/10 * a[3][0][0] # 0th order for i in range(3): @@ -245,10 +255,14 @@ def solve_coefficients(B_moments, xyz, order = 4): a = abc[i] Bx = B_moments[i] - a[0][0][0] = 1/2 * (Bx[0, x, y, z] + Bx[0][coords]) - 1/6 * a[2][0][0] - 1/70 * a[4][0][0] + a[0][0][0] = 1/2 * (Bx[0, x, y, z] + Bx[0][coords]) + if order > 1: + a[0][0][0] -= 1/6 * a[2][0][0] + if order > 3: + a[0][0][0] -= 1/70 * a[4][0][0] # Check constraint: - test = abc[0][1][0][0] + abc[1][1][0][0] + abc[2][1][0][0] + 1/10 * (abc[0][3][0][0] + abc[1][3][0][0] + abc[2][3][0][0]) + #test = abc[0][1][0][0] + abc[1][1][0][0] + abc[2][1][0][0] + 1/10 * (abc[0][3][0][0] + abc[1][3][0][0] + abc[2][3][0][0]) #if abs(test) > 0: # print("Something went wrong, sum (17) is " + str(test)) @@ -276,9 +290,9 @@ def neighbourdiff(a, idx): # Solves a, b, c components with a[x, y, z], b[y, z, x] and c[z, x, y] up to given order # Input B should be the output of solve_moments_from_B -def solve_all_coefficients(B_moments, order = 4): +def solve_all_coefficients(B_moments): shp = np.shape(B_moments[0][0]) - abc = np.zeros((3, 5, 4, 4, shp[0] - 1, shp[1] - 1, shp[2] - 1)) + abc = np.zeros((3, order+1, order, order, shp[0] - 1, shp[1] - 1, shp[2] - 1)) # 4th order if (order > 3): @@ -346,8 +360,12 @@ def solve_all_coefficients(B_moments, order = 4): a[0][1][0] = 1/2 * neighboursum(Bx[1], i) - 1/6 * a[2][1][0] a[0][0][1] = 1/2 * neighboursum(Bx[2], i) - 1/6 * a[2][0][1] - a[1][1][0] = neighbourdiff(Bx[1], i) - 1/10 * a[3][1][0] - a[1][0][1] = neighbourdiff(Bx[2], i) - 1/10 * a[3][0][1] + a[1][1][0] = neighbourdiff(Bx[1], i) + if order > 2: + a[1][1][0] -= 1/10 * a[3][1][0] + a[1][0][1] = neighbourdiff(Bx[2], i) + if order > 2: + a[1][0][1] -= 1/10 * a[3][0][1] for i in range(3): j = (i+1) % 3 @@ -357,7 +375,11 @@ def solve_all_coefficients(B_moments, order = 4): c = abc[k] # Should be correct, check later - a[2][0][0] = -1/2 * (b[1][0][1] + c[1][1][0]) - 3/35 * a[4][0][0] - 1/20 * (b[3][0][1] + c[3][1][0]) + a[2][0][0] = -1/2 * (b[1][0][1] + c[1][1][0]) + if order > 3: + a[2][0][0] -= 3/35 * a[4][0][0] + if order > 2: + a[2][0][0] -= 1/20 * (b[3][0][1] + c[3][1][0]) # 1st order if (order > 0): @@ -365,14 +387,20 @@ def solve_all_coefficients(B_moments, order = 4): a = abc[i] Bx = B_moments[i] - a[1][0][0] = neighbourdiff(Bx[0], i) - 1/10 * a[3][0][0] + a[1][0][0] = neighbourdiff(Bx[0], i) + if order > 2: + a[1][0][0] -= 1/10 * a[3][0][0] # 0th order for i in range(3): a = abc[i] Bx = B_moments[i] - a[0][0][0] = 1/2 * neighboursum(Bx[0], i) - 1/6 * a[2][0][0] - 1/70 * a[4][0][0] + a[0][0][0] = 1/2 * neighboursum(Bx[0], i) + if order > 1: + a[0][0][0] -= 1/6 * a[2][0][0] + if order > 3: + a[0][0][0] -= 1/70 * a[4][0][0] # Check constraint: #test = abc[0][1][0][0] + abc[1][1][0][0] + abc[2][1][0][0] + 1/10 * (abc[0][3][0][0] + abc[1][3][0][0] + abc[2][3][0][0]) @@ -381,55 +409,86 @@ def solve_all_coefficients(B_moments, order = 4): # print("Something went wrong, sum (17) is " + str(test)) # return np.pad(np.transpose(abc, (4, 5, 6, 0, 1, 2, 3)), [(0, 1), (0, 1), (0, 1), (0, 0), (0, 0), (0, 0), (0, 0)]) + return np.pad(abc, [(0, 0), (0, 0), (0, 0), (0, 0), (0, 1), (0, 1), (0, 1)]) -def center_value(B_moments, xyz, order=4): +def center_value(B_moments, xyz): abc = solve_coefficients(B_moments, xyz, order) return [interpolate_x(abc[0], 0, 0, 0), interpolate_y(abc[1], 0, 0, 0), interpolate_z(abc[2], 0, 0, 0)] -def ave_B(B_moments, xyz, order=4): - abc = solve_coefficients(B_moments, xyz, order) - a = abc[0] - b = abc[1] - c = abc[2] - return ( - a[0][0][0] - 3/8 * (a[2][0][0] + a[0][2][0] + a[0][0][2]) + 9/64 * (a[2][2][0] + a[2][0][2]) + 15/128 * a[4][0][0], - b[0][0][0] - 3/8 * (b[2][0][0] + b[0][2][0] + b[0][0][2]) + 9/64 * (b[2][2][0] + b[2][0][2]) + 15/128 * b[4][0][0], - c[0][0][0] - 3/8 * (c[2][0][0] + c[0][2][0] + c[0][0][2]) + 9/64 * (c[2][2][0] + c[2][0][2]) + 15/128 * c[4][0][0] - ) - -def all_ave_Bs(B_moments, order=4): - abc = solve_all_coefficients(B_moments, order) +def ave_B(B_moments, xyz): + abc = solve_coefficients(B_moments, xyz) a = abc[0] b = abc[1] c = abc[2] - return np.transpose(np.array( - ( + if order > 3: + return ( a[0][0][0] - 3/8 * (a[2][0][0] + a[0][2][0] + a[0][0][2]) + 9/64 * (a[2][2][0] + a[2][0][2]) + 15/128 * a[4][0][0], b[0][0][0] - 3/8 * (b[2][0][0] + b[0][2][0] + b[0][0][2]) + 9/64 * (b[2][2][0] + b[2][0][2]) + 15/128 * b[4][0][0], c[0][0][0] - 3/8 * (c[2][0][0] + c[0][2][0] + c[0][0][2]) + 9/64 * (c[2][2][0] + c[2][0][2]) + 15/128 * c[4][0][0] ) - ), [1, 2, 3, 0]) + elif order > 1: + return ( + a[0][0][0] - 3/8 * (a[2][0][0] + a[0][2][0] + a[0][0][2]) + 9/64 * (a[2][2][0] + a[2][0][2]), + b[0][0][0] - 3/8 * (b[2][0][0] + b[0][2][0] + b[0][0][2]) + 9/64 * (b[2][2][0] + b[2][0][2]), + c[0][0][0] - 3/8 * (c[2][0][0] + c[0][2][0] + c[0][0][2]) + 9/64 * (c[2][2][0] + c[2][0][2]) + ) + else: + return ( + a[0][0][0], + b[0][0][0], + c[0][0][0] + ) -def center_values(B_moments, coords, order=4): +def all_ave_Bs(B_moments): + abc = solve_all_coefficients(B_moments) + a = abc[0] + b = abc[1] + c = abc[2] + if order > 3: + return np.transpose(np.array( + ( + a[0][0][0] - 3/8 * (a[2][0][0] + a[0][2][0] + a[0][0][2]) + 9/64 * (a[2][2][0] + a[2][0][2]) + 15/128 * a[4][0][0], + b[0][0][0] - 3/8 * (b[2][0][0] + b[0][2][0] + b[0][0][2]) + 9/64 * (b[2][2][0] + b[2][0][2]) + 15/128 * b[4][0][0], + c[0][0][0] - 3/8 * (c[2][0][0] + c[0][2][0] + c[0][0][2]) + 9/64 * (c[2][2][0] + c[2][0][2]) + 15/128 * c[4][0][0] + ) + ), [1, 2, 3, 0]) + elif order > 1: + return np.transpose(np.array( + ( + a[0][0][0] - 3/8 * (a[2][0][0] + a[0][2][0] + a[0][0][2]) + 9/64 * (a[2][2][0] + a[2][0][2]), + b[0][0][0] - 3/8 * (b[2][0][0] + b[0][2][0] + b[0][0][2]) + 9/64 * (b[2][2][0] + b[2][0][2]), + c[0][0][0] - 3/8 * (c[2][0][0] + c[0][2][0] + c[0][0][2]) + 9/64 * (c[2][2][0] + c[2][0][2]) + ) + ), [1, 2, 3, 0]) + else: + return np.transpose(np.array( + ( + a[0][0][0], + b[0][0][0], + c[0][0][0] + ) + ), [1, 2, 3, 0]) + +def center_values(B_moments, coords): # Looks scuffed but is faster - abc = np.transpose(np.transpose(solve_all_coefficients(B_moments, order), (4, 5, 6, 0, 1, 2, 3))[coords[:, 0], coords[:, 1], coords[:, 2]], (1, 2, 3, 4, 0)) + abc = np.transpose(np.transpose(solve_all_coefficients(B_moments), (4, 5, 6, 0, 1, 2, 3))[coords[:, 0], coords[:, 1], coords[:, 2]], (1, 2, 3, 4, 0)) return np.transpose(np.array([interpolate_x(abc[0], 0, 0, 0), interpolate_y(abc[1], 0, 0, 0), interpolate_z(abc[2], 0, 0, 0)])) #return all_center_values(B_moments, order)[coords[:, 0], coords[:, 1], coords[:, 2], :] -def all_center_values(B_moments, order=4): - abc = solve_all_coefficients(B_moments, order) +def all_center_values(B_moments): + abc = solve_all_coefficients(B_moments) return np.transpose(np.array([interpolate_x(abc[0], 0, 0, 0), interpolate_y(abc[1], 0, 0, 0), interpolate_z(abc[2], 0, 0, 0)]), (1, 2, 3, 0)) -def all_center_values_dbxdi(B_moments, order=4): - abc = solve_all_coefficients(B_moments, order) +def all_center_values_dbxdi(B_moments): + abc = solve_all_coefficients(B_moments) return np.transpose(np.array([interpolate_dbxdx(abc[0], 0, 0, 0), interpolate_dbxdy(abc[0], 0, 0, 0), interpolate_dbxdz(abc[0], 0, 0, 0)]), (1, 2, 3, 0)) -def all_center_values_dbydi(B_moments, order=4): - abc = solve_all_coefficients(B_moments, order) +def all_center_values_dbydi(B_moments): + abc = solve_all_coefficients(B_moments) return np.transpose(np.array([interpolate_dbydx(abc[1], 0, 0, 0), interpolate_dbydy(abc[1], 0, 0, 0), interpolate_dbydz(abc[1], 0, 0, 0)]), (1, 2, 3, 0)) -def all_center_values_dbzdi(B_moments, order=4): - abc = solve_all_coefficients(B_moments, order) +def all_center_values_dbzdi(B_moments): + abc = solve_all_coefficients(B_moments) return np.transpose(np.array([interpolate_dbzdx(abc[2], 0, 0, 0), interpolate_dbzdy(abc[2], 0, 0, 0), interpolate_dbzdz(abc[2], 0, 0, 0)]), (1, 2, 3, 0)) From 6577dbf082d9d34c1fdd18d2d98dc5dc6b4403ca Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Thu, 10 Aug 2023 13:04:32 +0300 Subject: [PATCH 09/14] Some optimization to solve_moments_from_B --- scripts/divergenceless_reconstruction.py | 33 +++++++++++++++++------- 1 file changed, 23 insertions(+), 10 deletions(-) diff --git a/scripts/divergenceless_reconstruction.py b/scripts/divergenceless_reconstruction.py index f72388ac..7b6c1e91 100644 --- a/scripts/divergenceless_reconstruction.py +++ b/scripts/divergenceless_reconstruction.py @@ -109,19 +109,32 @@ def solve_moments_from_B(fg_b): # Returns moments for x as [order, x, y, z] # With order (0, y, z, yy, yz, zz, yyy, yyz, yzz, zzz) # y and z as cyclic rotations +# 6..9 required for order > 3 +# 3..5 required for order > 2 +# 1..2 required for order > 1 +# 0 required for order = 0 def solve_x_moments(B_x): + if order > 2: + n_moments = 10 + # order = 4 and order = 3 have also different n_moments, but order=2 is what we actually need, so look into this later + elif order == 2: + n_moments = 3 + else: + n_moments = 1 x_moments = np.zeros((10,) + np.shape(B_x)) x_moments[0] = B_x - x_moments[1:3] = np.gradient(B_x)[1:3] - start = 3 - for i in range(2): - j = i + 1 - x_moments[start:start+3-j] = np.gradient(x_moments[1+i])[j:3] - start += (3-j) - for i in range(3): - j = i + 1 if i < 2 else 2 - x_moments[start:start+3-j] = np.gradient(x_moments[3+i])[j:3] - start += (3-j) + if order > 1: + x_moments[1:3] = np.gradient(B_x)[1:3] + if order > 2: + start = 3 + for i in range(2): + j = i + 1 + x_moments[start:start+3-j] = np.gradient(x_moments[1+i])[j:3] + start += (3-j) + for i in range(3): + j = i + 1 if i < 2 else 2 + x_moments[start:start+3-j] = np.gradient(x_moments[3+i])[j:3] + start += (3-j) return x_moments def solve_y_moments(B_y): From 4d3efc8a262952e04c27ee9f35be30eb1d39698f Mon Sep 17 00:00:00 2001 From: Markku Alho Date: Tue, 14 May 2024 15:35:23 +0300 Subject: [PATCH 10/14] Fwd changes from xo_tools --- scripts/divergenceless_reconstruction.py | 41 +++++++++++------------- 1 file changed, 18 insertions(+), 23 deletions(-) diff --git a/scripts/divergenceless_reconstruction.py b/scripts/divergenceless_reconstruction.py index 7b6c1e91..2a5a139c 100644 --- a/scripts/divergenceless_reconstruction.py +++ b/scripts/divergenceless_reconstruction.py @@ -109,32 +109,19 @@ def solve_moments_from_B(fg_b): # Returns moments for x as [order, x, y, z] # With order (0, y, z, yy, yz, zz, yyy, yyz, yzz, zzz) # y and z as cyclic rotations -# 6..9 required for order > 3 -# 3..5 required for order > 2 -# 1..2 required for order > 1 -# 0 required for order = 0 def solve_x_moments(B_x): - if order > 2: - n_moments = 10 - # order = 4 and order = 3 have also different n_moments, but order=2 is what we actually need, so look into this later - elif order == 2: - n_moments = 3 - else: - n_moments = 1 x_moments = np.zeros((10,) + np.shape(B_x)) x_moments[0] = B_x - if order > 1: - x_moments[1:3] = np.gradient(B_x)[1:3] - if order > 2: - start = 3 - for i in range(2): - j = i + 1 - x_moments[start:start+3-j] = np.gradient(x_moments[1+i])[j:3] - start += (3-j) - for i in range(3): - j = i + 1 if i < 2 else 2 - x_moments[start:start+3-j] = np.gradient(x_moments[3+i])[j:3] - start += (3-j) + x_moments[1:3] = np.gradient(B_x)[1:3] + start = 3 + for i in range(2): + j = i + 1 + x_moments[start:start+3-j] = np.gradient(x_moments[1+i])[j:3] + start += (3-j) + for i in range(3): + j = i + 1 if i < 2 else 2 + x_moments[start:start+3-j] = np.gradient(x_moments[3+i])[j:3] + start += (3-j) return x_moments def solve_y_moments(B_y): @@ -505,3 +492,11 @@ def all_center_values_dbzdi(B_moments): abc = solve_all_coefficients(B_moments) return np.transpose(np.array([interpolate_dbzdx(abc[2], 0, 0, 0), interpolate_dbzdy(abc[2], 0, 0, 0), interpolate_dbzdz(abc[2], 0, 0, 0)]), (1, 2, 3, 0)) +def all_center_values_dbidi(B_moments): + abc = solve_all_coefficients(B_moments) + return np.concatenate( + (np.transpose(np.array([interpolate_dbxdx(abc[0], 0, 0, 0), interpolate_dbxdy(abc[0], 0, 0, 0), interpolate_dbxdz(abc[0], 0, 0, 0)]), (1, 2, 3, 0)), + np.transpose(np.array([interpolate_dbydx(abc[1], 0, 0, 0), interpolate_dbydy(abc[1], 0, 0, 0), interpolate_dbydz(abc[1], 0, 0, 0)]), (1, 2, 3, 0)), + np.transpose(np.array([interpolate_dbzdx(abc[2], 0, 0, 0), interpolate_dbzdy(abc[2], 0, 0, 0), interpolate_dbzdz(abc[2], 0, 0, 0)]), (1, 2, 3, 0))), + axis=3) + From 818686b001b70ea14652b267215b65be804d83b8 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Wed, 7 Aug 2024 14:51:03 +0300 Subject: [PATCH 11/14] Localization --- scripts/divergenceless_reconstruction.py | 44 ++++++++++++------------ 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/scripts/divergenceless_reconstruction.py b/scripts/divergenceless_reconstruction.py index 2a5a139c..4cdd73b4 100644 --- a/scripts/divergenceless_reconstruction.py +++ b/scripts/divergenceless_reconstruction.py @@ -268,7 +268,7 @@ def solve_coefficients(B_moments, xyz): return abc -def neighboursum(a, idx): +def neighborsum(a, idx): second = a[:-1, :-1, :-1] if idx == 0: first = a[1:, :-1, :-1] @@ -278,7 +278,7 @@ def neighboursum(a, idx): first = a[:-1, :-1, 1:] return first + second -def neighbourdiff(a, idx): +def neighbordiff(a, idx): second = a[:-1, :-1, :-1] if idx == 0: first = a[1:, :-1, :-1] @@ -300,15 +300,15 @@ def solve_all_coefficients(B_moments): a = abc[i] Bx = B_moments[i] - a[0][3][0] = 1/2 * neighboursum(Bx[6], i) - a[0][2][1] = 1/2 * neighboursum(Bx[7], i) - a[0][1][2] = 1/2 * neighboursum(Bx[8], i) - a[0][0][3] = 1/2 * neighboursum(Bx[9], i) + a[0][3][0] = 1/2 * neighborsum(Bx[6], i) + a[0][2][1] = 1/2 * neighborsum(Bx[7], i) + a[0][1][2] = 1/2 * neighborsum(Bx[8], i) + a[0][0][3] = 1/2 * neighborsum(Bx[9], i) - a[1][3][0] = neighbourdiff(Bx[6], i) - a[1][2][1] = neighbourdiff(Bx[7], i) - a[1][1][2] = neighbourdiff(Bx[8], i) - a[1][0][3] = neighbourdiff(Bx[9], i) + a[1][3][0] = neighbordiff(Bx[6], i) + a[1][2][1] = neighbordiff(Bx[7], i) + a[1][1][2] = neighbordiff(Bx[8], i) + a[1][0][3] = neighbordiff(Bx[9], i) for i in range(3): j = (i+1) % 3 @@ -330,13 +330,13 @@ def solve_all_coefficients(B_moments): a = abc[i] Bx = B_moments[i] - a[0][2][0] = 1/2 * neighboursum(Bx[3], i) - 1/6 * a[2][2][0] - a[0][1][1] = 1/2 * neighboursum(Bx[4], i) - a[0][0][2] = 1/2 * neighboursum(Bx[5], i) - 1/6 * a[2][0][2] + a[0][2][0] = 1/2 * neighborsum(Bx[3], i) - 1/6 * a[2][2][0] + a[0][1][1] = 1/2 * neighborsum(Bx[4], i) + a[0][0][2] = 1/2 * neighborsum(Bx[5], i) - 1/6 * a[2][0][2] - a[1][2][0] = neighbourdiff(Bx[3], i) - a[1][1][1] = neighbourdiff(Bx[4], i) - a[1][0][2] = neighbourdiff(Bx[5], i) + a[1][2][0] = neighbordiff(Bx[3], i) + a[1][1][1] = neighbordiff(Bx[4], i) + a[1][0][2] = neighbordiff(Bx[5], i) for i in range(3): j = (i+1) % 3 @@ -357,13 +357,13 @@ def solve_all_coefficients(B_moments): a = abc[i] Bx = B_moments[i] - a[0][1][0] = 1/2 * neighboursum(Bx[1], i) - 1/6 * a[2][1][0] - a[0][0][1] = 1/2 * neighboursum(Bx[2], i) - 1/6 * a[2][0][1] + a[0][1][0] = 1/2 * neighborsum(Bx[1], i) - 1/6 * a[2][1][0] + a[0][0][1] = 1/2 * neighborsum(Bx[2], i) - 1/6 * a[2][0][1] - a[1][1][0] = neighbourdiff(Bx[1], i) + a[1][1][0] = neighbordiff(Bx[1], i) if order > 2: a[1][1][0] -= 1/10 * a[3][1][0] - a[1][0][1] = neighbourdiff(Bx[2], i) + a[1][0][1] = neighbordiff(Bx[2], i) if order > 2: a[1][0][1] -= 1/10 * a[3][0][1] @@ -387,7 +387,7 @@ def solve_all_coefficients(B_moments): a = abc[i] Bx = B_moments[i] - a[1][0][0] = neighbourdiff(Bx[0], i) + a[1][0][0] = neighbordiff(Bx[0], i) if order > 2: a[1][0][0] -= 1/10 * a[3][0][0] @@ -396,7 +396,7 @@ def solve_all_coefficients(B_moments): a = abc[i] Bx = B_moments[i] - a[0][0][0] = 1/2 * neighboursum(Bx[0], i) + a[0][0][0] = 1/2 * neighborsum(Bx[0], i) if order > 1: a[0][0][0] -= 1/6 * a[2][0][0] if order > 3: From ed7cfab286bad9d5aaeb0eedb933553773f36218 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Wed, 7 Aug 2024 14:52:48 +0300 Subject: [PATCH 12/14] Reference Balsara --- scripts/divergenceless_reconstruction.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/scripts/divergenceless_reconstruction.py b/scripts/divergenceless_reconstruction.py index 4cdd73b4..43dce119 100644 --- a/scripts/divergenceless_reconstruction.py +++ b/scripts/divergenceless_reconstruction.py @@ -1,3 +1,6 @@ +# Reconstruction and interpolation of magntic field from face-average values +# Based on Balsara 2011 (https://doi.org/10.48550/arXiv.0811.2192) + import pytools as pt import numpy as np from scipy.special import legendre From 64ac396024593aaa0f649e8e3fb0bc9635fdd0a0 Mon Sep 17 00:00:00 2001 From: lkotipal Date: Wed, 7 Aug 2024 15:01:46 +0300 Subject: [PATCH 13/14] Better test script --- scripts/divergenceless_reconstruction_test.py | 61 +++++++++++-------- 1 file changed, 36 insertions(+), 25 deletions(-) mode change 100644 => 100755 scripts/divergenceless_reconstruction_test.py diff --git a/scripts/divergenceless_reconstruction_test.py b/scripts/divergenceless_reconstruction_test.py old mode 100644 new mode 100755 index c3ad71d6..7b9aee3a --- a/scripts/divergenceless_reconstruction_test.py +++ b/scripts/divergenceless_reconstruction_test.py @@ -1,35 +1,46 @@ +#!/usr/bin/env python from divergenceless_reconstruction import solve_moments_from_B, center_value, center_values, all_center_values, all_ave_Bs import pytools as pt import numpy as np from time import perf_counter +from sys import argv # Code for testing divergenceless_resconstruction.py -f = pt.vlsvfile.VlsvReader('/wrk/users/lkotipal/TEST/bulk.0000057.vlsv') -fg_b = f.read_fsgrid_variable('fg_b') -print(np.shape(fg_b)) -fg_b_vol = f.read_fsgrid_variable('fg_b_vol') -print('File read!') -t = perf_counter() -B_moments = solve_moments_from_B(fg_b[0:100, 0:100, 0:100]) -print(f'B_moments solved in {perf_counter() - t} seconds!') -print() -print(f'Volumetric b in (50, 50, 50): {fg_b_vol[50, 50, 50]}') -print() -for i in range(5): - t = perf_counter() - print(center_value(B_moments, (50, 50, 50), i)) - print(f'Single coefficients up to order {i} solved in {perf_counter() - t} seconds!') +def main(): + if len(argv) < 2: + print("Usage: ./divergenceless_reconstruction_test.py FILENAME") + exit(1) + filename = argv[1] + f = pt.vlsvfile.VlsvReader(filename) + fg_b = f.read_fsgrid_variable('fg_b') + print(np.shape(fg_b)) + fg_b_vol = f.read_fsgrid_variable('fg_b_vol') + print('File read!') t = perf_counter() - print(np.shape(center_values(B_moments, np.array([(50, 50, 50), (51, 51, 51), (54, 54, 64), (43, 65, 42), (43, 65, 32)]), i))) - print(f'Multiple coefficients up to order {i} solved in {perf_counter() - t} seconds!') + B_moments = solve_moments_from_B(fg_b[0:100, 0:100, 0:100]) + print(f'B_moments solved in {perf_counter() - t} seconds!') + print() + print(f'Volumetric b in (50, 50, 50): {fg_b_vol[50, 50, 50]}') + print() + for i in range(5): + t = perf_counter() + print(center_value(B_moments, (50, 50, 50), i)) + print(f'Single coefficients up to order {i} solved in {perf_counter() - t} seconds!') - t = perf_counter() - print(all_center_values(B_moments, i)[50, 50, 50]) - print(f'All coefficients up to order {i} solved in {perf_counter() - t} seconds!') + t = perf_counter() + print(np.shape(center_values(B_moments, np.array([(50, 50, 50), (51, 51, 51), (54, 54, 64), (43, 65, 42), (43, 65, 32)]), i))) + print(f'Multiple coefficients up to order {i} solved in {perf_counter() - t} seconds!') + + t = perf_counter() + print(all_center_values(B_moments, i)[50, 50, 50]) + print(f'All coefficients up to order {i} solved in {perf_counter() - t} seconds!') + + print(f'Average scaled difference between Analysator and Vlasiator in order {i}:') + print(np.average(((all_ave_Bs(B_moments, i))[:-1, :-1, :-1] - fg_b_vol[:99, :99, :99])/fg_b_vol[:99,:99,:99])) + print(f'Max scaled difference between Analysator and Vlasiator in order {i}:') + print(np.max(abs(((all_ave_Bs(B_moments, i))[:-1, :-1, :-1] - fg_b_vol[:99, :99, :99])/fg_b_vol[:99,:99,:99]))) + print() - print(f'Average scaled difference between Analysator and Vlasiator in order {i}:') - print(np.average(((all_ave_Bs(B_moments, i))[:-1, :-1, :-1] - fg_b_vol[:99, :99, :99])/fg_b_vol[:99,:99,:99])) - print(f'Max scaled difference between Analysator and Vlasiator in order {i}:') - print(np.max(abs(((all_ave_Bs(B_moments, i))[:-1, :-1, :-1] - fg_b_vol[:99, :99, :99])/fg_b_vol[:99,:99,:99]))) - print() \ No newline at end of file +if __name__ == "__main__": + main() \ No newline at end of file From 4c0a4eaaa4b818e50d34512b6f13f76aedfe6dac Mon Sep 17 00:00:00 2001 From: lkotipal Date: Wed, 7 Aug 2024 15:13:33 +0300 Subject: [PATCH 14/14] Exit codes --- scripts/divergenceless_reconstruction_test.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/scripts/divergenceless_reconstruction_test.py b/scripts/divergenceless_reconstruction_test.py index 7b9aee3a..18005871 100755 --- a/scripts/divergenceless_reconstruction_test.py +++ b/scripts/divergenceless_reconstruction_test.py @@ -9,7 +9,7 @@ def main(): if len(argv) < 2: print("Usage: ./divergenceless_reconstruction_test.py FILENAME") - exit(1) + return(1) filename = argv[1] f = pt.vlsvfile.VlsvReader(filename) @@ -42,5 +42,9 @@ def main(): print(np.max(abs(((all_ave_Bs(B_moments, i))[:-1, :-1, :-1] - fg_b_vol[:99, :99, :99])/fg_b_vol[:99,:99,:99]))) print() + return 0 + if __name__ == "__main__": - main() \ No newline at end of file + ret = main() + if(ret): + exit(ret)