处理可变长度input/output数组，满足条件时终止计算（边界层计算）
Handling variable length input/output arrays, terminating computation when a condition is met (Boundary layer calculations)
我正在尝试编写代码来计算机翼上的边界层发展。
我从描述机翼形状的一组点 S 开始。
我运行 对这些点进行非粘性求解并提取驻点，这将是S.
的一个元素
点集S现在被驻点分成两组su和sl，描述了上下边界层下翼型的形状。
这是我的第一个问题。假设我写了一个组件 BoundaryLayerSolve，它采用边界层中点的矢量 s 和矢量 ue 边界层的边缘速度。 s 和 ue 的长度相同。如果我想使用这个组件两次，一次用于机翼的每一侧，我需要先验地知道停滞点的位置，直到模型已经建立并且 运行 for无粘求解器。我该如何设置才能处理这些未知的输入数组大小？
现在输入的s和ue数组是已知的上下边界层，可以计算边界层。我将使用两种模型，一种用于边界层的层流区域，一种用于湍流区域。假设这两个模型使用完全不同的计算，并且这些计算非常复杂，以至于为了找到它们的解析偏导数，它们必须被分解成子组件。
层流计算从停滞点开始，沿着 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)
这是我的第二个问题。直到在层流计算期间才知道过渡点，因此层流计算的输出 momentum_thickness 数组的长度是先验未知的。因此，用于湍流计算的输入点阵列的长度也不是先验已知的。我该如何解决这个问题？
这是我的第三个问题。如何让一个组件在满足条件后终止其计算，并传递给另一个组件以完成对数组其余部分的计算？
收集我的问题：
- 如何让一个组件从一个初始大小已知的数组中取出一个初始大小未知的切片？
- 当两个未知数组加起来等于初始已知数组的大小时，如何让一个组件输出一个初始未知大小的数组，并让第二个组件输入第二个初始未知大小的数组？
- 如何让组件在满足条件时终止计算，并将计算传递给第二个组件以完成剩余的计算？
我知道这是一个很长的问题，可能会被分解。提前感谢您的阅读。
对于问题 1 和 2：目前 OpenMDAO 不处理执行期间输入和输出的动态大小，并且目前没有任何更改此设置的计划。我们提供 shape_by_conn 以允许变量根据其 sources/targets 进行整形，但我认为这不是您想要的，因为双方在您的公式中都未确定。
问题 3：如果我们隐含地处理这个问题，那么我们可以强制两个计算之间的过渡发生在层流和湍流区域之间的交界处。例如，在 Dymos 中，当我们传播轨迹时，我们不会像典型的时间步进模拟那样使用事件触发器。相反，我们在轨迹“阶段”的开始或结束处使用约束来强制转换条件发生在交界处。
我倾向于尝试以这种方式表述问题：
使用插值法提供机翼上的点作为某个独立变量的函数。 S = f(x)。想象一下，当 x 发生变化时，用红色圈出的点在机翼周围不断滑动。
假设过渡点是先验已知的，那么我们有两个主要的计算组：层流和湍流。每组评估机翼上的一些点（N_l 和 N_u）。确定过渡点的参数 x 的值可以来回“滑动”，直到 x 的假定过渡点值与实际期望值匹配（通过在过渡点强制执行残差或约束，使用 x作为隐式变量或设计变量）。
与其在过渡点将层流部分的输出馈入湍流部分，不如将这些值视为湍流部分中的自变量，然后强制它们与流段末尾的值相匹配层流部分（通过设置残差，或作为优化器的约束）。
根据实施的具体细节，您可能必须绑定假定值才能从每个部分获得有效计算。如果您使用优化方法而不是求解器来制定此函数，我也不确定您将在此处用作 objective 函数。
我要回答你的问题有点乱：
How can I have a component output an initially unknown size array, and have a second component input a second initially unknown size array when both unknown arrays add up to the size of an initially known array?
从 OpenMDAO V3.6.0 开始，您不能这样做。 OpenMDAO 要求您告诉它所有输入和输出的大小。通过 shape_by_conn feature 有一个轻微的解决方法，它可以让下游组件由上游组件调整大小，但最终 sizing-chain 的开始需要在设置时具有给定值。因此，您不能在 运行 时间内动态调整输出值的大小。但是，我将为您提供一个您可以改用的解决方法。
你的问题 1 和 3 都可以用一个叫做 over-allocation 的技巧来解决。
通常，这意味着分配一个比您实际需要的更大的数组，并添加一个额外的输入来帮助您跟踪使用了多少。具体如何使用这个技巧有点问题，但这是一个显示基础知识的通用示例：
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', VectorChain(vec_size=VEC_SIZE, step_size=STEP_SIZE))
c2 = p.model.add_subsystem('c2', VectorChain(vec_size=VEC_SIZE, 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 个节点，则可以 over-allocate 一个长度为 1000 的数组并将其传递给整个计算（这基本上就是我的示例所做的）。或者，您可能知道湍流跳变点总是在第 20 点和第 50 点之间。在这种情况下，您可以 over-allocate 一个长度为 50 的数组用于层流计算的第一个分量，一个长度为 980 的数组用于第二个湍流分量。如果您有大网格，其中潜在的重叠要窄得多，这将节省一些内存。
关于衍生品的注意事项：
整个概念涉及指数的离散运动，从技术上讲是 non-differentiable。在连续导数上下文中进行这项工作的方法是假设对于任何线性化（任何时候计算导数）索引值都是固定的（因此不参与导数）。
在您的情况下，层流计算结束的特定索引值将从一个 运行 变为下一个，但任何时候您想知道导数，您都认为它是固定的。
在严格的数学意义上，如果您使用基于梯度的优化，则此假设无效。函数确实不可微分。然而，在实践中，我们通常可以指望假设设计收敛并且更改索引的位置有效固定。然而，如果这种稳定没有发生（也许你有一个 really-really 精细的离散化，并且有足够的噪声无论如何都可以让跳闸点在几个节点周围跳跃）那么你可能会遇到紧密收敛的问题.
我认为这个 discrete-index 问题不太可能成为您的问题，所以我鼓励您试一试！但是，如果您确实发现使用这种方法时问题过于嘈杂，那么您应该 re-formulating 以更连续的方式考虑问题。
更连续的公式：
使用沿机翼的离散数据，定义连续插值函数（您可以使用 OpenMDAO 的 MetaModelStructured component, or 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 的任何实际值都有效，介于网格的正负限制之间。
一旦你有了那个插值，你就可以在插值的连续 space 中求解停滞点位置（即找到你需要的 s 的值）。为此，您可能需要对插值进行微分，然后求解该导数变为零的位置。这将是第二个隐式组件，紧随第一个输出平滑拟合系数的组件之后。
然后你可以将插值系数和停滞点的位置传递给一个层流分量，并从该点在上下表面上向后执行你的行进计算ace 直到你到达过渡点。
尽管您可能会选择分步进行以生成漂亮均匀的 spaced 数据输出网格，但您仍然希望在过渡点方面保持连续性。假设您在 s 中以固定步长前进，直到超出转换条件的第一个点。那时，您将在最后两个搜索位置之间开始二分搜索过程，以找到 s 中过渡点的“精确”（在一些合理的公差范围内）位置。需要明确的是，层流和湍流计算将自行执行 sub-grids，这些计算是使用插值函数开发的。
然后您可以将确切过渡位置的条件传递到下游的下一个组件（连同其 s 位置）并开始另一个带有湍流计算的行进计算。
我会注意到我的建议与 Rob Falck 给出的答案类似，尽管我的方法不需要隐式关系。你可以用我的显式方式或他的隐式方式来设置问题。两者各有优缺点，这里就不赘述了。关键是即使使用显式计算方法，也可以实现真正连续的公式。这种方法的一个好处是它更容易保持固定长度的数组。您可以为每个层流和湍流分配 n 个离散点 sub-grids 并让物理间距随着驻点和 laminar/turbulent 点移动而稍微变化。
在实践中，我认为您不必尝试使用连续方法。我认为区分它更正确，并且能够更紧密地收敛。虽然数学正确性很好，但并不总是实用。在实践中，我认为您可以找到一个合适的网格以更离散的格式工作，并且可以更快地启动和 运行ning。如果你发现你有噪音问题或难以收敛，那么你可以稍后改变方法。
我能想到的一件事可能会让您对离散方法的噪音敏感，如果您计划在边界层求解器和无粘性求解器之间创建一个 feed-back 循环。一些代码会将边界层厚度作为新的翼型形状并将其传递回无粘性计算，然后 re-compute 表面压力分布并重新计算原始翼型表面上的新边界层。如果你打算尝试这样的事情，离散方法的噪音可能会给你带来更多麻烦。在这种情况下，更连续的方法是必要的。
我正在尝试编写代码来计算机翼上的边界层发展。
我从描述机翼形状的一组点 S 开始。
我运行 对这些点进行非粘性求解并提取驻点，这将是S.
的一个元素点集S现在被驻点分成两组su和sl，描述了上下边界层下翼型的形状。
这是我的第一个问题。假设我写了一个组件 BoundaryLayerSolve，它采用边界层中点的矢量 s 和矢量 ue 边界层的边缘速度。 s 和 ue 的长度相同。如果我想使用这个组件两次，一次用于机翼的每一侧，我需要先验地知道停滞点的位置，直到模型已经建立并且 运行 for无粘求解器。我该如何设置才能处理这些未知的输入数组大小？
现在输入的s和ue数组是已知的上下边界层，可以计算边界层。我将使用两种模型，一种用于边界层的层流区域，一种用于湍流区域。假设这两个模型使用完全不同的计算，并且这些计算非常复杂，以至于为了找到它们的解析偏导数，它们必须被分解成子组件。
层流计算从停滞点开始，沿着 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)
这是我的第二个问题。直到在层流计算期间才知道过渡点，因此层流计算的输出 momentum_thickness 数组的长度是先验未知的。因此，用于湍流计算的输入点阵列的长度也不是先验已知的。我该如何解决这个问题？
这是我的第三个问题。如何让一个组件在满足条件后终止其计算，并传递给另一个组件以完成对数组其余部分的计算？
收集我的问题：
- 如何让一个组件从一个初始大小已知的数组中取出一个初始大小未知的切片？
- 当两个未知数组加起来等于初始已知数组的大小时，如何让一个组件输出一个初始未知大小的数组，并让第二个组件输入第二个初始未知大小的数组？
- 如何让组件在满足条件时终止计算，并将计算传递给第二个组件以完成剩余的计算？
我知道这是一个很长的问题，可能会被分解。提前感谢您的阅读。
对于问题 1 和 2：目前 OpenMDAO 不处理执行期间输入和输出的动态大小，并且目前没有任何更改此设置的计划。我们提供 shape_by_conn 以允许变量根据其 sources/targets 进行整形，但我认为这不是您想要的，因为双方在您的公式中都未确定。
问题 3：如果我们隐含地处理这个问题，那么我们可以强制两个计算之间的过渡发生在层流和湍流区域之间的交界处。例如，在 Dymos 中，当我们传播轨迹时，我们不会像典型的时间步进模拟那样使用事件触发器。相反，我们在轨迹“阶段”的开始或结束处使用约束来强制转换条件发生在交界处。
我倾向于尝试以这种方式表述问题：
使用插值法提供机翼上的点作为某个独立变量的函数。 S = f(x)。想象一下，当
x发生变化时，用红色圈出的点在机翼周围不断滑动。假设过渡点是先验已知的，那么我们有两个主要的计算组：层流和湍流。每组评估机翼上的一些点（
N_l和N_u）。确定过渡点的参数x的值可以来回“滑动”，直到 x 的假定过渡点值与实际期望值匹配（通过在过渡点强制执行残差或约束，使用x作为隐式变量或设计变量）。与其在过渡点将层流部分的输出馈入湍流部分，不如将这些值视为湍流部分中的自变量，然后强制它们与流段末尾的值相匹配层流部分（通过设置残差，或作为优化器的约束）。
根据实施的具体细节，您可能必须绑定假定值才能从每个部分获得有效计算。如果您使用优化方法而不是求解器来制定此函数，我也不确定您将在此处用作 objective 函数。
我要回答你的问题有点乱：
How can I have a component output an initially unknown size array, and have a second component input a second initially unknown size array when both unknown arrays add up to the size of an initially known array?
从 OpenMDAO V3.6.0 开始，您不能这样做。 OpenMDAO 要求您告诉它所有输入和输出的大小。通过 shape_by_conn feature 有一个轻微的解决方法，它可以让下游组件由上游组件调整大小，但最终 sizing-chain 的开始需要在设置时具有给定值。因此，您不能在 运行 时间内动态调整输出值的大小。但是，我将为您提供一个您可以改用的解决方法。
你的问题 1 和 3 都可以用一个叫做 over-allocation 的技巧来解决。 通常，这意味着分配一个比您实际需要的更大的数组，并添加一个额外的输入来帮助您跟踪使用了多少。具体如何使用这个技巧有点问题，但这是一个显示基础知识的通用示例：
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', VectorChain(vec_size=VEC_SIZE, step_size=STEP_SIZE))
c2 = p.model.add_subsystem('c2', VectorChain(vec_size=VEC_SIZE, 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 个节点，则可以 over-allocate 一个长度为 1000 的数组并将其传递给整个计算（这基本上就是我的示例所做的）。或者，您可能知道湍流跳变点总是在第 20 点和第 50 点之间。在这种情况下，您可以 over-allocate 一个长度为 50 的数组用于层流计算的第一个分量，一个长度为 980 的数组用于第二个湍流分量。如果您有大网格，其中潜在的重叠要窄得多，这将节省一些内存。
关于衍生品的注意事项：
整个概念涉及指数的离散运动，从技术上讲是 non-differentiable。在连续导数上下文中进行这项工作的方法是假设对于任何线性化（任何时候计算导数）索引值都是固定的（因此不参与导数）。 在您的情况下，层流计算结束的特定索引值将从一个 运行 变为下一个，但任何时候您想知道导数，您都认为它是固定的。
在严格的数学意义上，如果您使用基于梯度的优化，则此假设无效。函数确实不可微分。然而，在实践中，我们通常可以指望假设设计收敛并且更改索引的位置有效固定。然而，如果这种稳定没有发生（也许你有一个 really-really 精细的离散化，并且有足够的噪声无论如何都可以让跳闸点在几个节点周围跳跃）那么你可能会遇到紧密收敛的问题.
我认为这个 discrete-index 问题不太可能成为您的问题，所以我鼓励您试一试！但是，如果您确实发现使用这种方法时问题过于嘈杂，那么您应该 re-formulating 以更连续的方式考虑问题。
更连续的公式： 使用沿机翼的离散数据，定义连续插值函数（您可以使用 OpenMDAO 的 MetaModelStructured component, or 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 的任何实际值都有效，介于网格的正负限制之间。
一旦你有了那个插值，你就可以在插值的连续 space 中求解停滞点位置（即找到你需要的 s 的值）。为此，您可能需要对插值进行微分，然后求解该导数变为零的位置。这将是第二个隐式组件，紧随第一个输出平滑拟合系数的组件之后。
然后你可以将插值系数和停滞点的位置传递给一个层流分量，并从该点在上下表面上向后执行你的行进计算ace 直到你到达过渡点。 尽管您可能会选择分步进行以生成漂亮均匀的 spaced 数据输出网格，但您仍然希望在过渡点方面保持连续性。假设您在 s 中以固定步长前进，直到超出转换条件的第一个点。那时，您将在最后两个搜索位置之间开始二分搜索过程，以找到 s 中过渡点的“精确”（在一些合理的公差范围内）位置。需要明确的是，层流和湍流计算将自行执行 sub-grids，这些计算是使用插值函数开发的。
然后您可以将确切过渡位置的条件传递到下游的下一个组件（连同其 s 位置）并开始另一个带有湍流计算的行进计算。
我会注意到我的建议与 Rob Falck 给出的答案类似，尽管我的方法不需要隐式关系。你可以用我的显式方式或他的隐式方式来设置问题。两者各有优缺点，这里就不赘述了。关键是即使使用显式计算方法，也可以实现真正连续的公式。这种方法的一个好处是它更容易保持固定长度的数组。您可以为每个层流和湍流分配 n 个离散点 sub-grids 并让物理间距随着驻点和 laminar/turbulent 点移动而稍微变化。
在实践中，我认为您不必尝试使用连续方法。我认为区分它更正确，并且能够更紧密地收敛。虽然数学正确性很好，但并不总是实用。在实践中，我认为您可以找到一个合适的网格以更离散的格式工作，并且可以更快地启动和 运行ning。如果你发现你有噪音问题或难以收敛，那么你可以稍后改变方法。
我能想到的一件事可能会让您对离散方法的噪音敏感，如果您计划在边界层求解器和无粘性求解器之间创建一个 feed-back 循环。一些代码会将边界层厚度作为新的翼型形状并将其传递回无粘性计算，然后 re-compute 表面压力分布并重新计算原始翼型表面上的新边界层。如果你打算尝试这样的事情，离散方法的噪音可能会给你带来更多麻烦。在这种情况下，更连续的方法是必要的。