Расчет расстояний между точками в сетке CFD

Question

У меня есть структурированная сетка, которая имеет около 3e + 6 очков и

Пожалуйста, просмотрите изображение: https://i.stack.imgur.com/eUUkQ.jpg

Каждая точка в физической области (евклидова) имеет индекс ( i, j и k ) в вычислительной области.

Мне нужно перебрать индексы вычислительной области и выполнить вычисления с соответствующими точками.

Например, длина в направлении i при данном индексе будет (псевдокод): длина = vec_len (точка ( i + 1, j, k ) - точка ( i, j, k ))

Также важно выяснить соотношение длин в данной точке. Например, я бы вычислил две соседние длины в направлении i и разделил бы их между собой.

То, что я придумал, занимает слишком много времени и, вероятно, не использует весь потенциал, который может предложить NumPy.

Я сделал ndarray, заполненный нулями, который будет содержать все координаты XYZ сетки.

block_data =numpy.zeros((i_dim,  j_dim, k_dim, 3), dtype='float')

Число 3 соответствует 3 элементам, x, y и z.

Так что, если бы я хотел, чтобы значения z составляли i = 3, j = 7, k = 10 , это было бы:

Z = block_data[3][7][10][2]

Точкой в ​​евклидовом пространстве будет (1,3) ndarray:

point = block_data[i][j][k]

Способ вычисления длины между двумя точками:

numpy.linalg.norm(point2 - point1)

Только расчет длины занимает около 1,5 мс, и я хотел бы рассчитать расстояния во всех точках и во всех направлениях: 3e + 6 * 3.

Я думаю, что есть проблема с подходом, как я строю основной блок ndarray (block_data), потому что это ограничивает меня в расчете только по двум точкам за раз, то есть только по двум маленьким (1,3 ) ndarrays. Если я правильно помню, производить вычисления на небольших массивах не очень эффективно.

Как я могу подойти к проблеме и ускорить время выполнения? Есть ли книга, рекомендуемая для такого рода проблем? Спасибо: -)

Answer 1

Чтобы вычислить евклидово расстояние в массиве такого размера, я рекомендую использовать векторизованный подход:

def euclid_dist(array, direction):
    if direction == 'i':  # make shifted views depending on the direction
        p1 = array[:-1, :, :]
        p2 = array[1:, :, :]
    elif direction == 'j':
        p1 = array[:, :-1, :]
        p2 = array[:, 1:, :]
    elif direction == 'k':
        p1 = array[:, :, :-1]
        p2 = array[:, :, 1:]
    else:
        raise ValueError('direction ' + direction + ' not known.')

    # get euclidean distance for all points in direction:
    euc_dist = (((p1 - p2)*(p1 - p2)).sum(axis=3))**0.5

    return euc_dist

Использование небольшого тестового массива с:

arr = np.random.randint(-20, 20, 5*5*5*3).reshape(5, 5, 5, 3)
eu_i = euclid_dist(arr, 'i')
eu_j = euclid_dist(arr, 'j')
# test some values:
print(eu_i[2, 1, 2] == np.linalg.norm(arr[2, 1, 2] - arr[3, 1, 2]))
# Out 64: True
print(eu_j[1, 1, 1] == np.linalg.norm(arr[1, 1, 1] - arr[1, 2, 1]))
# Out 65: True

Некоторые тайминги для большого массива с 8e6 точками и 24e6 значениями:

big_arr = np.random.rand(200, 200, 200, 3)
%timeit euclid_dist(big_arr, 'i')
# 644 ms ± 57.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Имхо, это довольно быстро для массива такого размера. :) Если я правильно читаю ваши настройки, это примерно в 19000 раз быстрее, чем ваш код.