处理可变长度的输入/输出数组，满足条件时终止计算边界层计算

如何解决处理可变长度的输入/输出数组，满足条件时终止计算边界层计算

我正在尝试编写代码来计算翼型上的边界层发展。

我从一组描述机翼形状的点开始，S。

我在这些点上运行一个无粘性求解器并提取驻点，这将是 S 的一个元素。

点集合 S 现在被驻点分为两组 s_u 和 s_l，描述了上下边界层下方的翼型形状。

这是我的第一个问题。假设我编写了一个组件 BoundaryLayerSolve，它采用边界层中点的向量 s 和边缘速度的向量 u_e 到它边界层。 s 和 u_e 的长度相同。如果我想两次使用这个组件，对机翼的每一侧使用一次，我需要先验地知道停滞点的位置，直到模型已经建立并为无粘性求解器运行后才能找到。我该如何设置才能处理这些未知的输入数组大小？

现在输入 s 和 u_e 数组对于上边界层和下边界层是已知的，可以计算边界层。我将使用两种模型，一种用于边界层的层流区域，一种用于湍流区域。假设这两个模型使用完全不同的计算，并且这些计算足够复杂，以至于为了找到它们的解析偏导数，它们必须分解为子组件。

层流计算从停滞点开始，沿着s前进。在每一步，计算动量厚度的值。如果达到过渡阈值，则必须停止层流计算，并且必须在剩余的s 中使用湍流计算。我要强调的是，不是对所有 s 进行层流，对所有 s 进行湍流，然后简单地切掉层流输出的后部，然后湍流输出的前部。湍流计算必须从过渡点层流计算的值开始。

在伪代码中，这将类似于：

transition_point = 0
for index,point in enumerate(s):
    transition_point += 1
    momentum_thickness = laminar(point)
    if momentum_thickness > transition_threshold:
        break

for point in s[transition_point+1:]:    
    turbulent(s)

这是我的第二个问题。过渡点直到层流计算期间才知道，因此层流计算的输出动量_厚度数组的长度是先验未知的。因此，用于湍流计算的输入点阵列的长度也不是先验已知的。我该如何解决这个问题？

这是我的第三个问题。如何在满足条件后让组件终止其计算，并将其传递给另一个组件以完成对数组其余部分的计算？

收集我的问题：

如何让组件从初始已知大小的数组中获取初始未知大小的切片？
当两个未知数组加起来等于初始已知数组的大小时，如何让组件输出初始未知大小的数组，并让第二个组件输入第二个初始未知大小的数组？
如何在满足条件时让组件终止计算，并将计算传递给第二个组件以完成剩余的计算？

我明白这是一个很长的问题，可能会被分解。感谢您提前阅读。

解决方法

我来回答你的问题有点乱：

当两个未知数组加起来等于初始已知数组的大小时，如何让组件输出初始未知大小的数组，并让第二个组件输入第二个初始未知大小的数组？

从 OpenMDAO V3.6.0 开始，您不能这样做。 OpenMDAO 要求您告诉它所有输入和输出的大小。通过 shape_by_conn feature 有一个轻微的解决方法，它将让下游组件由上游组件调整大小，但最终该调整链的开始需要在设置时具有给定的值。因此，您不能在运行时动态调整输出值的大小。但是，我将为您提供一个您可以使用的变通方法。

您的问题 1 和 3 都可以通过一种称为过度分配的技巧来解决。通常，这意味着分配比您实际需要的更大的数组，并添加一个额外的输入来帮助您跟踪要使用的数量。具体如何使用这个技巧有点问题，但这里是一个展示基础知识的通用示例：

import numpy as np
import openmdao.api as om


class VectorChain(om.ExplicitComponent): 

    def initialize(self):
        self.options.declare('vec_size',types=int)
        self.options.declare('step_size',types=int)

    def setup(self): 
        vec_size = self.options['vec_size']
        step_size = self.options['vec_size']

        # this is the index of the last valid value in the data array
        self.add_input('in_idx',shape=1,val=0) 
        # NOTE: though this could be done as a discrete variable,# that will confuse OpenMDAO's derivative system,so just pass it as a float. 

        # actual data array
        self.add_input('in_vec',shape=vec_size,val=np.zeros(vec_size))

        self.add_output('out_idx',shape=1)
        self.add_output('out_vec',shape=vec_size)

        # NOTES: You do **NOT** define derivatives wrt the idx variable!! 
        self.declare_partials('out_vec','in_vec',method='CS')

    def compute(self,inputs,outputs): 

        in_vec = inputs['in_vec']
        out_vec = outputs['out_vec']
        
        out_vec[:] = in_vec.copy()

        # always use the first two entries to
        # indicate the start/end of the valid data
        i_start_idx = int(inputs['in_idx'])
        i_end_idx = i_start_idx + self.options['step_size']
        i_size = i_end_idx - i_start_idx

        if i_start_idx == 0: 
            out_vec[0] = 1 # get the counting started
            i_start_idx = 1

        # note: careful to account for the open end of the  
        # interval when computing the end_idx value
        # print(self.pathname)
        for i in range(i_start_idx,i_start_idx+self.options['step_size']): 
            out_vec[i] = out_vec[i-1] + 1

        o_end_idx = i_start_idx + i_size

        outputs['out_idx'] = o_end_idx


if __name__ == "__main__": 

    p = om.Problem()

    VEC_SIZE = 12
    STEP_SIZE = 3

    c0 = p.model.add_subsystem('c0',VectorChain(vec_size=VEC_SIZE,step_size=STEP_SIZE))
    c1 = p.model.add_subsystem('c1',step_size=STEP_SIZE))
    c2 = p.model.add_subsystem('c2',step_size=STEP_SIZE))

    p.model.connect('c0.out_idx','c1.in_idx')
    p.model.connect('c0.out_vec','c1.in_vec')

    p.model.connect('c1.out_idx','c2.in_idx')
    p.model.connect('c1.out_vec','c2.in_vec')


    p.setup()

    p.run_model()

    p.model.list_outputs(print_arrays=True)

在这个玩具示例中，它们是一连串相同的组件，每个组件都从最后一个停止的地方开始计数。这显然是一个微不足道的例子，但基础知识就在那里。在您的情况下，层流和过渡计算会有不同的组件，但想法是相同的。从层流部分开始，循环直到达到 trop 条件。停在那里，并注意您停在的索引。然后将该信息向下游传递给下一个组件。

这里有一些细节值得注意：

如果您愿意，可以在一个组件链中同时制作上表面和下表面。不需要这样做，但组件越少通常速度越快。
如果上表面有 1000 个节点，则可以过度分配长度为 1000 的数组并将其传递给整个计算（这基本上就是我的示例所做的）。或者，您可能知道湍流跳闸点总是在第 20 和第 50 点之间。在这种情况下，您可以为层流计算的第一个分量过度分配长度为 50 的数组，为第二个湍流分量分配一个长度为 980 的数组。如果您有大网格，其中潜在的重叠要窄得多，这将节省一些内存。

关于衍生品的说明：

整个概念涉及指数的离散运动，这在技术上是不可微的。在连续导数上下文中进行这项工作的方法是假设对于任何线性化（任何时候计算导数），索引值都是固定的（因此不参与导数）。在您的情况下，层流计算结束的特定索引值会从一次运行更改为下一次运行，但任何时候您想知道导数，您都假设它是固定的。

在严格的数学意义上，如果您使用基于梯度的优化，则此假设无效。函数确实是不可微的。然而，在实践中，我们通常可以指望设计收敛和变化索引的位置变得有效固定的假设。然而，如果这种稳定没有发生（也许你有一个非常非常好的离散化并且有足够的噪声来保持跳变点无论如何都跳几个节点）那么你可能会遇到紧密收敛的问题。

我认为这个离散索引问题不太可能对您的情况造成影响，因此我鼓励您尝试一下！但是，如果您确实发现使用这种方法时问题过于嘈杂，那么您应该考虑以更连续的方式重新制定问题。

更连续的公式： 使用沿翼型的离散数据，定义连续插值函数（您可以使用 OpenMDAO 的 MetaModelStructured 组件，或 use a newton-solver to converge a least-squares fit to a pre-defined basis function）。插值的输入不是 x,y 位置，而是沿翼型表面的参数位置 s（即从任意起点开始的路径长度）。有几种方法可以解决这个问题，但为了简单起见，假设您在机翼前端选择一个离散点为 0，然后机翼顶部向后的每个点都是正数 s=1,s=2， s=3,s=4,etc. 下表面的每个点都是负数 s=-1,s=-2,s=-3,s=-4 etc. 虽然我选择了整数值，但它很重要意识到插值对网格的正负限制之间的任何 s 实际值都有效。

一旦你有了那个插值，你就可以在插值的连续空间中求解驻点位置（即找到你需要的 s 的值）。为此，您可能需要对插值进行微分，然后求解该导数在哪里变为零。这将是第二个隐式组件，在第一个输出平滑拟合系数之后。

然后，您可以将插值系数和停滞位置传递给层流分量，并从上下表面上的该点向后执行行进计算，直到到达过渡点。尽管您可能选择以离散的步骤前进以生成良好的均匀分布的数据输出网格，但您仍然希望在过渡点方面保持连续。假设您以 s 中的固定步骤前进，直到超出转换条件的第一个点。那时，您将在最后两个搜索位置之间开始二分搜索过程，以在 s 中找到过渡点的“精确”（在一些合理的容差范围内）位置。需要说明的是，层流和湍流计算将在它们自己的子网格上进行，这些子网格是使用插值函数开发的。

然后，您可以将准确转换位置的条件传递到下游的下一个组件（连同其s 位置），并使用湍流计算开始另一个行进计算。

我会注意到我的建议类似于 Rob Falck 给出的答案，尽管我的方法不需要隐含关系。你可以用我明确的方式或他的隐含方式来设置问题。两者各有优缺点，这里就不赘述了。关键是，即使使用明确的计算方法，也可以实现真正连续的公式化。这种方法的一个好处是它更容易保持固定长度的数组。您可以为每个层流和湍流子网格分配 n 个离散点，并让物理间距随着驻点和层流/湍流点的移动而略有不同。

在实践中，我认为您不必尝试使用连续方法。我认为区分它更正确，并且能够更紧密地收敛。虽然数学正确性很好，但并不总是实用。在实践中，我认为您可以找到一个合适的网格，以更离散的格式工作并更快地启动和运行。如果您发现有噪音问题或难以收敛，那么您可以稍后切换方法。

我能想到的一件事可能会让您对离散方法的噪声敏感，那就是您是否计划在边界层求解器和无粘性求解器之间创建反馈循环。一些代码将边界层厚度作为新的翼型形状并将其传递回无粘性计算，然后重新计算表面压力分布并重新计算原始翼型表面上的新边界层。如果您打算尝试这样的事情，离散方法的噪音可能会给您带来更多麻烦。在这种情况下，需要采用更持续的方法。

对于问题 1 和 2：目前 OpenMDAO 在执行期间不处理输入和输出的动态大小调整，并且目前没有任何计划来改变这一点。我们提供 shape_by_conn 以允许根据其来源/目标来塑造变量，但我认为这不是您在这里想要的，因为在您的公式中双方都未确定。

问题 3：如果我们隐式地处理问题，那么我们可以强制两个计算之间的转换发生在层流区域和湍流区域之间的连接处。例如，在 Dymos 中，当我们传播轨迹时，我们不会像典型的时间步进模拟那样使用事件触发器。相反，我们在轨迹“阶段”的开始或结束处使用约束来强制转换条件发生在连接处。

我倾向于尝试以这种方式来表述问题：

使用插值提供翼型上的点作为某个自变量的函数。 S = f(x)。想象一下，当 x 发生变化时，红色圆圈中的点围绕机翼连续滑动。
假设过渡点是先验已知的，所以我们有两个主要的计算组：层流和湍流。每组评估机翼上的一些点（N_l 和 N_u）。确定过渡点的参数 x 的值可以来回“滑动”，直到 x 的假定过渡点值与实际期望值匹配（通过在过渡点强制残差或约束，用 {{ 1}} 作为隐式变量或设计变量）。
不是将输出从层流部分馈入过渡点处的湍流部分，而是将这些值视为湍流部分中的自变量，然后强制它们与层流部分末端的值（或通过设置残差，或作为优化器的约束）。

根据实现的特定细节，您可能必须限制假定值才能从每个部分中获得有效的计算。如果您使用优化方法而不是求解器来制定目标函数，我也不确定您会在这里使用什么作为目标函数。

处理可变长度的输入/输出数组，满足条件时终止计算边界层计算

如何解决处理可变长度的输入/输出数组，满足条件时终止计算边界层计算

解决方法

相关推荐