Я получаю MemError, когда я интерполирую данные с помощью выражений генератора и итераторов из регулярной сетки на сетке> 2000 раз - PullRequest
Новая
       20

Я получаю MemError, когда я интерполирую данные с помощью выражений генератора и итераторов из регулярной сетки на сетке> 2000 раз

0 голосов
/ 31 октября 2018

это мой первый вопрос здесь, в stackoverflow, потому что я начал писать скрипты с Python3.

Применение

Я создал скрипт Python3, который записывает определение нагрузки подвижного источника тепла для моделирования методом конечных элементов в LS-Dyna. В качестве источника у меня есть дискретизированное трехмерное поле плотности скорости тепловыделения (Вт / см ^ 3), координаты, определяющие сетку конечных элементов и положение центра теплового поля во времени. В качестве вывода я получаю зависящую от времени мощность нагрева, отсортированную после номера элемента для каждого конечного элемента. Это работает уже для разумных размеров (200000 конечных элементов, 3000 местоположений теплового поля, 400000 точек данных в тепловом поле).

Задача

Для больших сеток с конечными элементами (4 000 000 элементов) у меня не хватает памяти (60 ГБ ОЗУ, python3 64Bit). Чтобы проиллюстрировать проблему дальше, я подготовил минимальный пример, который работает сам по себе. Он генерирует некоторые искусственные тестовые данные, конечно-элементную сетку, как я ее использую (на самом деле это не обычная сетка) и итератор для новых мест для применения тепла. 

import numpy as np
import math
from scipy.interpolate import RegularGridInterpolator
def main():
    dataCoordinateAxes,dataArray = makeTestData()
    meshInformationArray = makeSampleMesh()
    coordinates = makeSampleCoordinates()
    interpolateOnMesh(dataCoordinateAxes,dataArray,meshInformationArray,coordinates)

def makeTestData():
    x = np.linspace(-0.02,0.02,300)
    y = np.linspace(-0.02,0.02,300)
    z = np.linspace(-0.005,0.005,4)
    data = f(*np.meshgrid(x,y,z,indexing='ij',sparse=True))
    return (x,y,z),data

def f(x,y,z):
    scaling = 1E18
    sigmaXY = 0.01
    muXY = 0
    sigmaZ = 0.5
    muZ = 0.005
    return weight(x,1E-4,muXY,sigmaXY)*weight(y,1E-4,muXY,sigmaXY)*weight(z,0.1,muZ,sigmaZ)*scaling

def weight(x,dx,mu,sigma):
    result = np.multiply(np.divide(np.exp(np.divide(np.square(np.subtract(x,mu)),(-2*sigma**2))),math.sqrt(2*math.pi*sigma**2.)),dx)
    return result

def makeSampleMesh():
    meshInformation = []
    for x in np.linspace(-0.3,0.3,450):
        for y in np.linspace(-0.3,0.3,450):
            for z in np.linspace(-0.005,0.005,5):
                meshInformation.append([x,y,z])
    return np.array(meshInformation)

def makeSampleCoordinates():
    x = np.linspace(-0.2,0.2,500)
    y = np.sqrt(np.subtract(0.2**2,np.square(x)))
    return (np.array([element[0],element[1],0])for element in zip(x,y))

Затем в этой функции выполняется интерполяция. Я удалил все в цикле for, чтобы изолировать проблему. На самом деле я экспортирую кривую загрузки в файл определенного формата. 

def interpolateOnMesh(dataCoordinateAxes,dataArray,meshInformationArray,coordinates):
    interpolationFunction = RegularGridInterpolator(dataCoordinateAxes, dataArray, bounds_error=False, fill_value=None)
    for finiteElementNumber, heatGenerationCurve in enumerate(iterateOverFiniteElements(meshInformationArray, coordinates, interpolationFunction)):
        pass
    return

def iterateOverFiniteElements(meshInformationArray, coordinates, interpolationFunction):
    meshDataIterator = (np.nditer(interpolationFunction(np.subtract(meshInformationArray,coordinateSystem))) for coordinateSystem in coordinates)
    for heatGenerationCurve in zip(*meshDataIterator):
        yield heatGenerationCurve

if __name__ == '__main__':
    main()

Чтобы выявить проблему, я отслеживал потребление памяти с течением времени. потребление памяти с течением времени Похоже, что итерация по массивам результатов занимает значительное количество памяти.

Вопрос

Существует ли менее ресурсоемкий способ перебора точек данных без потери производительности? Если нет, то я предполагаю, что нарежу массив сетки на куски и интерполирую их по одному.

1 Ответ

0 голосов
/ 02 ноября 2018

Пока единственное решение, которое я нашел, было сократить meshInformationArray. Вот модифицированная функция main(): 

def main():
    dataCoordinateAxes,dataArray = makeTestData()
    meshInformationArray = makeSampleMesh()
    coordinates = makeSampleCoordinates()
    sections = int(meshInformationArray.shape[0] / 100000)
    if sections == 0: sections = 1
    for array in iter(np.array_split(meshInformationArray, sections, axis=0)):
        interpolateOnMesh(dataCoordinateAxes,dataArray,array,coordinates)
...