Gestion des tableaux d'entrée / sortie de longueur variable, arrêt de la computation lorsqu'une condition est remplie (calculs de la couche limite)

Question

Je suis en train d'essayer d'écrire du code pour calculer le développement de la couche limite sur un profil aérodynamique.

Je commence avec un ensemble de points, S, qui décrivent la forme du profil aérodynamique.

Je lance un solveur non visqueux sur ces points et j'extrais le point de stagnation, qui sera un élément de S.

L'ensemble de points S est maintenant divisé par le point de stagnation en deux ensembles s_u et s_l, qui décrivent la forme du profil aérodynamique sous les couches limites supérieure et inférieure.

Voici mon premier problème. Supposons que j'écrive un composant, BoundaryLayerSolve, qui prend un vecteur s de points dans la couche limite et un vecteur u_e des vitesses de bord de cette couche limite. Les deux s et u_e seront de la même longueur. Si je voulais utiliser ce composant deux fois, une fois pour chaque côté du profil aérodynamique, je devrais connaître a priori la position du point de stagnation, qui ne peut pas être trouvé tant que le modèle n'est pas déjà configuré et exécuté pour le solveur non visqueux. Comment pourrais-je le configurer pour que je puisse gérer ces tailles de tableau d'entrée inconnues?

Maintenant que les tableaux d'entrée s et u_e sont connus pour les couches limites supérieure et inférieure, les couches limites peuvent être calculées. J'utiliserai deux modèles, un pour la région laminaire de la couche limite et un pour la région turbulente. Supposons que ces deux modèles utilisent des calculs totalement différents, et que ces calculs sont suffisamment complexes pour que, pour trouver leurs dérivées partielles analytiques, ils doivent être décomposés en sous-composants.

Les calculs laminaires commencent au point de stagnation et avancent le long de s. À chaque étape, la valeur de l'épaisseur de momentum est calculée. Si le seuil de transition est atteint, les calculs laminaires doivent s'arrêter et des calculs turbulents doivent être utilisés pour le reste de s. Je tiens à souligner qu'il n'est pas possible de faire du laminair pour toute la longueur de s, du turbulent pour toute la longueur de s, puis simplement trancher la partie arrière de la sortie laminar et la partie avant de la sortie turbulente. Le calcul turbulent doit commencer avec les valeurs du calcul laminar au point de transition.

En pseudocode, cela ressemblerait à quelque chose comme:

point_de_transition = 0
pour index, point in enumerate(s):
    point_de_transition += 1
    epaisseur_momentum = laminar(point)
    si epaisseur_momentum > seuil_transition:
        break

pour point in s[point_de_transition+1:]:    
    turbulent(s)

Voici mon deuxième problème. Le point de transition n'est pas connu avant pendant les calculs laminaires, et donc la longueur du tableau de sortie des épaisseurs de momentum pour les calculs laminaires n'est pas connue à priori. Par conséquent, la longueur du tableau de points en entrée pour les calculs turbulents n'est pas non plus connue à priori. Comment puis-je contourner cela?

Voici mon troisième problème. Comment puis-je faire en sorte qu'un composant arrête son calcul une fois qu'une condition est remplie, et transmette le calcul à un autre composant pour terminer le calcul restant du tableau?

---

Pour résumer mes questions:

1. Comment puis-je faire en sorte qu'un composant prenne une tranche de taille initialement inconnue à partir d'un tableau de taille initialement connue?
2. Comment puis-je faire en sorte qu'un composant produise un tableau de taille initialement inconnue, et qu'un deuxième composant prenne en entrée un deuxième tableau de taille initialement inconnue lorsque les deux tableaux inconnus ajoutent jusqu'à la taille d'un tableau initialement connu?
3. Comment puis-je faire en sorte qu'un composant interrompe un calcul une fois qu'une condition est remplie, et transmette le calcul à un deuxième composant pour terminer le calcul restant?

Je sais que ceci est une longue question, et pourrait peut-être être subdivisée. Merci de l'avoir lue, d'avance.

Answer 1

Je vais répondre à vos questions un peu dans le désordre :

Comment puis-je faire en sorte qu'un composant produise un tableau de taille initialement inconnue, et qu'un second composant introduise un second tableau de taille initialement inconnue lorsque les deux tableaux inconnus s'additionnent à la taille d'un tableau initialement connu ?

Depuis OpenMDAO V3.6.0, vous ne pouvez plus le faire. OpenMDAO exige que vous lui indiquiez la taille de toutes les entrées et sorties. Il existe une légère solution de contournement via l'option caractéristique shape_by_conn ce qui permettra aux composants en aval d'être dimensionnés par les comps en amont, mais en fin de compte le début de cette chaîne de dimensionnement doit avoir une valeur donnée au moment de la configuration. Vous ne pouvez donc pas dimensionner dynamiquement les valeurs de sortie pendant l'exécution. Cependant, je vais vous fournir une solution de contournement que vous pouvez utiliser à la place.

Vos questions 1 et 3 peuvent toutes deux être traitées par une astuce appelée surallocation. En général, cela signifie allouer un tableau plus grand que ce dont vous avez réellement besoin, et ajouter une entrée supplémentaire qui vous aide à garder la trace de la quantité à utiliser. La manière exacte d'utiliser cette astuce est un peu spécifique au problème, mais voici un exemple générique qui montre les bases :

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)

Dans cet exemple de jouet, il y a une chaîne de composants identiques et chacun d'eux continue à compter à partir de l'endroit où le dernier s'est arrêté. C'est évidemment un exemple trivial, mais les bases sont là. Dans votre cas, vous auriez des composants différents pour les calculs laminaires et de transition, mais l'idée est la même. Commencez par la partie laminaire, bouclez jusqu'à ce que vous atteigniez la condition trop. Arrêtez-vous là, et notez l'indice auquel vous vous êtes arrêté. Puis transmettez cette information en aval au composant suivant.

Quelques détails méritent d'être notés :

1. Si vous le vouliez, vous pourriez faire les surfaces supérieure et inférieure ensemble dans une chaîne de composants. Il n'est pas nécessaire de le faire, mais moins de composants est généralement plus rapide.
2. Si vous aviez 1000 nœuds sur la surface supérieure, vous pourriez surallouer un tableau de longueur 1000 et le faire passer dans tout le calcul (c'est en gros ce que fait mon exemple). Par ailleurs, vous pouvez savoir que le point de déclenchement de la turbulence se situera toujours entre les 20e et 50e points. Dans ce cas, vous pourriez sur-allouer un tableau de longueur 50 pour la première composante pour les calculs laminaires, et un tableau de longueur 980 pour la seconde composante turbulente. Si vous aviez de grandes grilles, où le chevauchement potentiel était beaucoup plus étroit, cela permettrait d'économiser de la mémoire.

Une remarque sur les produits dérivés :

Tout ce concept implique un mouvement discret d'indices qui est techniquement non différentiable. La façon dont vous faites fonctionner cela dans le contexte des dérivées continues est de supposer que pour toute linéarisation (à chaque fois que vous calculez les dérivées), les valeurs des indices sont fixes (et ne participent donc pas aux dérivées). Dans votre cas, la valeur de l'indice spécifique de la fin du calcul laminaire changera d'une exécution à l'autre, mais à chaque fois que vous voulez connaître les dérivées, vous supposez qu'elle est fixe.

Au sens mathématique strict, si vous utilisiez une optimisation basée sur le gradient, cette hypothèse n'est pas valable. La fonction n'est vraiment pas différentiable. Cependant, dans la pratique, nous pouvons souvent compter sur l'hypothèse que le modèle converge et que l'emplacement de l'indice changeant devient effectivement fixe. Cependant, si cette stabilisation ne se produit pas (peut-être avez-vous une discrétisation vraiment très fine et qu'il y a suffisamment de bruit pour que le point de déclenchement continue à sauter de quelques nœuds quoi qu'il arrive), vous risquez d'avoir des problèmes de convergence étroite.

Je ne pense pas qu'il soit probable que cette question de l'indice discret soit un problème pour votre cas, et je vous encourage donc à essayer ! Cependant, si vous trouvez que le problème est trop bruyant en utilisant cette approche, alors vous devriez envisager de reformuler le problème d'une manière plus continue.

Une formulation plus continue : en utilisant les données discrètes le long de l'aile, définissez une fonction d'interpolation continue (vous pouvez utiliser la fonction d'interpolation continue d'OpenMDAO). MétaModèleStructuré ou utiliser un solveur de Newton pour faire converger un ajustement des moindres carrés vers une fonction de base prédéfinie. ). L'entrée de votre interpolant ne serait pas un emplacement x,y mais plutôt un emplacement paramétrique s le long de la surface de l'aile (c'est-à-dire une longueur de trajectoire à partir d'un point de départ arbitraire). Il y a plusieurs façons de gérer cela, mais pour rester simple, disons que vous choisissez un point discret sur le nez de l'aile comme 0, puis chaque point vers l'arrière sur le dessus de l'aile est positif s=1, s=2, s=3, s=4, etc. Chaque point sur la surface inférieure est négatif s=-1, s=-2, s=-3, s=-4, etc. Bien que j'aie choisi des valeurs entières, il est important de réaliser que l'interpolant est valable pour toute valeur réelle de s entre les limites positives et négatives de votre grille.

Une fois que vous avez cette interpolation, vous pouvez résoudre l'emplacement du point de stagnation dans l'espace continu de l'interpolant (c'est-à-dire trouver la valeur de s dont vous avez besoin). Pour ce faire, vous devrez probablement différencier l'interpolant puis résoudre l'endroit où cette dérivée devient nulle. Il s'agirait d'une deuxième composante implicite, après la première qui fournit les coefficients de votre ajustement lisse.

Ensuite, vous pouvez transmettre les coefficients interpolants et l'emplacement de la stagnation à une composante laminaire, et effectuer vos calculs à rebours à partir de ce point sur les surfaces supérieure et inférieure jusqu'à ce que vous atteigniez le point de transition. Bien que vous puissiez choisir de procéder par étapes discrètes pour produire une grille de données de sortie uniformément espacée, vous souhaitez que les choses restent continues en ce qui concerne le point de transition. Supposons que vous marchiez avec un pas fixe en s jusqu'au premier point qui est au-delà de la condition de transition. À ce moment-là, vous lancerez un processus de recherche par bissection entre les deux derniers emplacements de recherche pour trouver l'emplacement "exact" (dans une certaine tolérance raisonnable) du point de transition dans la zone de transition. s . Pour être clair, les calculs laminaires et turbulents seraient effectués sur leurs propres sous-grilles, qui ont été développées en utilisant la fonction d'interpolation.

Ensuite, vous pourriez transmettre les conditions à la position exacte de la transition au composant suivant en aval (ainsi que son s ) et commencer un autre calcul de marche avec les calculs de turbulence.

Je noterai que ma suggestion est similaire à la réponse donnée par Rob Falck, bien que mon approche ne nécessite pas de relation implicite. Vous pouvez résoudre le problème de ma façon explicite ou de sa façon implicite. Les deux approches présentent des avantages et des inconvénients que je n'aborderai pas ici. Le point essentiel est qu'une formulation véritablement continue est possible, même avec une approche de calcul explicite. L'un des avantages de cette approche est qu'il est beaucoup plus facile de conserver des tableaux de longueur fixe. Vous pouvez simplement allouer n points discrets pour chacune des sous-grilles laminaire et turbulente et laisser l'espacement physique varier un peu lorsque le point de stagnation et le point laminaire/turbulent se déplacent.

En pratique, je ne pense pas que vous deviez essayer d'utiliser l'approche continue. Je pense qu'il est plus correct de la différencier et qu'elle serait capable de converger plus étroitement. Si l'exactitude mathématique est agréable, elle n'est pas toujours pratique. En pratique, je pense que vous pourriez trouver une grille appropriée pour travailler dans un format plus discret et être opérationnel plus rapidement. Si vous constatez que vous avez des problèmes de bruit ou des difficultés à converger, vous pourrez changer d'approche par la suite.

Une chose à laquelle je pense qui pourrait vous rendre sensible au bruit de l'approche discrète serait si vous prévoyez de créer une boucle de rétroaction entre le solveur de couche limite et le solveur inviscide. Certains codes prennent l'épaisseur de la couche limite comme une nouvelle forme d'aile et la renvoient aux calculs inviscides, puis recalculent les distributions de pression de surface et recalculent une nouvelle couche limite sur la surface de l'aile originale. Si vous voulez essayer quelque chose comme ça, le bruit de l'approche discrète vous causera probablement plus de problèmes. Dans ce cas, l'approche plus continue est justifiée.

Answer 2

Pour les questions 1 et 2 : Actuellement, OpenMDAO ne gère pas le dimensionnement dynamique des entrées et des sorties pendant l'exécution, et il n'y a actuellement aucun projet de changer cela. Nous offrons shape_by_conn pour permettre aux variables d'être dimensionnées en fonction de leurs sources/cibles, mais je ne pense pas que cela soit ce que vous voulez ici, puisque les deux côtés sont indéterminés dans votre formulation.

Question 3 : Si nous traitons le problème implicitement, alors nous pouvons forcer la transition entre les deux calculs à se produire au niveau de la jonction entre les régions laminaires et turbulentes. Par exemple, dans Dymos lorsque nous propageons une trajectoire, nous n'utilisons pas de déclencheurs d'événements comme le ferait une simulation de pas de temps traditionnelle. Au lieu de cela, nous utilisons des contraintes au début ou à la fin d'une "phase" de la trajectoire pour forcer la condition de transition à se produire à la jonction.

Mon inclination est d'essayer de formuler le problème de cette manière :

1. Utiliser l'interpolation pour fournir le point sur le profilé d'aile en fonction d'une variable indépendante. S = f(x). Imaginez que le point entouré de rouge glisse continuellement autour du profilé d'aile lorsque x est modifié.

2. Supposer que le point de transition est connu a priori, donc nous avons deux principaux groupes de calcul : laminaires et turbulents. Chaque groupe évalue un certain nombre de points sur le profilé d'aile (N_l et N_u). La valeur du paramètre x qui détermine le point de transition peut "glisser" d'avant en arrière de sorte que jusqu'à ce que la valeur de x du point de transition supposé corresponde à la valeur désirée réelle (en imposant des résidus ou des contraintes au point de transition, avec x en tant que variable implicite ou variable de conception).

3. Au lieu de transmettre les sorties de la portion laminaires à la portion turbulente au point de transition, considérez ces valeurs comme des variables indépendantes dans la section turbulente, puis forcez-les à correspondre aux valeurs à la fin de la section laminaires (soit en mettant en place des résidus, soit en tant que contraintes pour l'optimiseur).

En fonction des détails particuliers de la mise en œuvre, vous pourriez devoir contraindre les valeurs supposées pour obtenir des calculs valides de chaque section. Je ne suis pas non plus sûr de ce que vous utiliseriez comme fonction objective ici, si vous formuliez cela en utilisant une approche d'optimisation au lieu d'un solveur.