使用python中的modred模块进行动态模式分解

Question

我正在尝试使用python中的modred模块编写动态模式分解代码。我编写了CustomVector类和CustomVecHandle类，其中包含以下modred文档：

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替换data_array.T时，我得到了同样的错误。那么，有没有人有一个好主意来解决这个错误？

Answer 1

我在使用相同模块计算POD时遇到了同样的问题。您正在处理的两个数组的形状为（2，）和（28330，）。尝试将它们转换为numpy数组，其中指定了两个维度。这将解决问题。