在 python 中使用 modred 模块进行动态模式分解答案

【问题标题】：Dynamic Mode Decomposition using the modred module in python在 python 中使用 modred 模块进行动态模式分解
【发布时间】：2018-03-25 10:06:27
【问题描述】：

我正在尝试使用 python 中的 modred 模块编写用于动态模式分解的代码。我按照 modred 文档编写 CustomVector 类和 CustomVecHandle 类，如下所示：

class CustomVector(mr.Vector):   
    def __init__(self, grids, data_array):
        self.grids = grids
        self.data_array = data_array
        self.weighted_ip = mr.InnerProductTrapz(*self.grids)

    def __add__(self, other):
        """Return a new object that is the sum of self and other"""
        sum_vec = deepcopy(self)   
        sum_vec.data_array = self.data_array + other.data_array
        return sum_vec
    def __mul__(self, scalar):
        """Return a new object that is ``self * scalar`` """
        mult_vec = deepcopy(self)
        mult_vec.data_array = mult_vec.data_array * scalar
        return mult_vec
    def inner_product(self, other): 
        return  self.weighted_ip(self.data_array, other.data_array)

class CustomVecHandle(mr.VecHandle):
    def __init__(self, vec_path, base_handle=None, scale=None):
        mr.VecHandle.__init__(self, base_handle, scale)
        self.vec_path = vec_path

    def _get(self):
        # read in the values
        print ("reading data from {}".format(self.vec_path))
        reader = vtk.vtkPolyDataReader()
        reader.SetFileName(self.vec_path)
        reader.Update()
        grids = dsa.WrapDataObject(reader.GetOutput()).Points
        data_array = dsa.WrapDataObject(reader.GetOutput()).PointData['U']
        return CustomVector(grids, data_array)

    def _put(self, vec):
        print ("writing data to {}".format(self.vec_path))
        common_writer.SetFileName(self.vec_path)   
        U_common_reader[:]=vec.data_array
        ier1 = common_writer.Write() # return 1 means success
        grids_common_reader[:]=vec.grids
        ier2 = common_writer.Write()
        if ier1!=1 or ier2!=1:
            raise Error()
        return

另外一个函数在所有类之外定义为：

def inner_product(v1, v2):
    return v1.inner_product(v2)

我创建了一个列表对象：

vec_handles=[CustomVecHandle(os.path.join(data_root,d,"U_zNormal.vtk")) for d in dirs]

然后我使用下面的代码来计算模式：

myDMD=mr.DMDHandles(inner_product,max_vecs_per_node=50)
myDMD.compute_decomp(vec_handles)

但是当我运行代码时，我收到以下错误：

Traceback (most recent call last):
  File "vtkDMD.py", line 192, in <module>
    eigvals,Rlo_eigvecs,Llo_eigvecs,cm_eigvals,cl_eigvecs = myDMD.compute_decomp(vec_handles)
  File "C:\Python27\lib\site-packages\modred\dmd.py", line 679, in compute_decomp
    self.vec_handles)
  File "C:\Python27\lib\site-packages\modred\vectorspace.py", line 495, in compute_symmetric_inner_product_mat
    IP_burn = self.inner_product(test_vec, test_vec)
  File "vtkDMD.py", line 142, in inner_product
    return v1.inner_product(v2)
  File "vtkDMD.py", line 138, in inner_product
    return  self.weighted_ip(self.data_array, other.data_array)
  File "C:\Python27\lib\site-packages\modred\vectors.py", line 159, in __call__
    return self.inner_product(vec1, vec2)
  File "C:\Python27\lib\site-packages\modred\vectors.py", line 168, in inner_product
    IP = np.trapz(IP, x=grid)
  File "C:\pv54\bin\lib\site-packages\numpy\lib\function_base.py", line 2941, in trapz
    ret = add.reduce(d * (y[slice1]+y[slice2])/2.0, axis)
ValueError: operands could not be broadcast together with shapes (2,)(28330,)

首先我认为这是因为给函数提供了错误的参数：

return CustomVector(grids, data_array)

然后我将grids替换为grids.T，表示网格的转置，但同样的错误再次出现。当我尝试用data_array.T 替换data_array 时，我得到了同样的错误。那么，是否有人有解决此错误的好主意？

【问题讨论】：

标签： python vtk dmd

【解决方案1】：

我在使用相同模块计算 POD 时遇到了同样的问题。您正在处理的两个数组的形状为 (2,) & (28330,)。尝试将它们转换为指定两个维度的 numpy 数组。这样就可以解决问题了。

【讨论】：