Неожиданное переключение координат при изменении положения объекта в трехмерной сетке после выравнивания и интеграции. ⇐ Python

[Anonymous]

Я работаю над обработкой и анализом данных поля скоростей из файлов VTI с использованием Python. Рабочий процесс включает в себя чтение файлов VTI, извлечение данных о скорости, вычисление траекторий объектов, создание траекторий в виде карт, а затем выполнение различных вычислений. Однако я столкнулся с проблемой, когда на картах создается впечатление, что оси смещены, что приводит к неверным расчетам дельты между соседними точками сетки.
Я использую интеграцию np.meshgrid объектов. или частицы изменяют форму и сглаживаются для интеграции, а затем принимают желаемую форму в зависимости от оси домена. Я не смог найти точку, в которой переключались оси x и z. Я мог сказать, что они переключались из-за того, что разница между одной точкой и соседом была равна нулю вдоль желаемой оси, но вдоль противоположной оси, поэтому z для x и x для z я получал ожидаемые координатные расстояния. По оси Y не возникло никаких проблем.
Кроме того, я заметил, что если я изменю x_grid, скажем, на 1–11 вместо 0–10, то рассчитанная там дельта будет такой, как и ожидалось. , то есть форма, которая, по моему мнению, подходит. Включив «order = F» для изменения формы, чтобы следовать порядку индексов, подобному Фортрану, во время процесса изменения формы сохраняется правильное выравнивание осей.
Почему изменение формы не работает как ожидается, если я не буду использовать индексацию Fortran? Или я неправильно понял или испортил изменение формы на другом этапе интеграции?
Мне удалось воссоздать проблему в следующем фрагменте:
import numpy as np
from scipy.integrate import solve_ivp

x_grid = np.linspace(0, 1, 10)
y_grid = np.linspace(0, 1, 10)
z_grid = np.linspace(0, 1, 10)
resolution = 0.01

domain_x = np.arange(x_grid.min(), x_grid.max() + resolution, resolution)
domain_y = np.arange(y_grid.min(), y_grid.max() + resolution, resolution)
domain_z = np.arange(z_grid.min(), z_grid.max() + resolution, resolution)

domain_shape = (len(domain_x), len(domain_y), len(domain_z))

#initial positions
initial_positions = np.array(np.meshgrid(domain_x, domain_y, domain_z, indexing='ij')).T.reshape(-1, 3)

def path_calc(t, y):
return y
# Solving
t_span = (0, 1)
times = np.linspace(0, 1, 10)
sol = solve_ivp(
path_calc,
t_span,
initial_positions.flatten(),
method='RK45',
t_eval=times,
vectorized=True
)
trajectories = sol.y.T.reshape(len(times), -1, 3)

def compute_maps(trajectories, domain_shape):
positions_final = trajectories[-1]

expected_num_positions = np.prod(domain_shape)
actual_num_positions = positions_final.shape[0]

if expected_num_positions != actual_num_positions:
raise ValueError(f"Expected number of positions {expected_num_positions} does not match actual number {actual_num_positions}")

map_x = positions_final[:, 0].reshape(domain_shape, order='F')
map_y = positions_final[:, 1].reshape(domain_shape, order='F')
map_z = positions_final[:, 2].reshape(domain_shape, order='F')

map_x_non = positions_final[:, 0].reshape(domain_shape)
map_y_non = positions_final[:, 1].reshape(domain_shape)
map_z_non = positions_final[:, 2].reshape(domain_shape)

maps = {
'map_x': map_x,
'map_y': map_y,
'map_z': map_z,
'map_x_non': map_x_non,
'map_y_non': map_y_non,
'map_z_non': map_z_non
}

return maps

maps = compute_maps(trajectories, domain_shape)

delta_x_calc = np.abs(maps['map_x'][1, 0, 0] - maps['map_x'][0, 0, 0])
delta_y_calc = np.abs(maps['map_y'][0, 1, 0] - maps['map_y'][0, 0, 0])
delta_z_calc = np.abs(maps['map_z'][0, 0, 1] - maps['map_z'][0, 0, 0])
delta_x_calc_non = np.abs(maps['map_x_non'][1, 0, 0] - maps['map_x_non'][0, 0, 0])
delta_y_calc_non = np.abs(maps['map_y_non'][0, 1, 0] - maps['map_y_non'][0, 0, 0])
delta_z_calc_non = np.abs(maps['map_z_non'][0, 0, 1] - maps['map_z_non'][0, 0, 0])


Подробнее здесь: https://stackoverflow.com/questions/787 ... d-grid-aft