Convergence (Quantitative)

telemac2d.py [STUDY-NAME].cas -s

Ce chapitre utilise les fichiers de simulation de Telemac steady 2d tutorial, avec une définition modifiée de l’étape temporelle et des périodes d’impression:

/ steady2d-conv.cas
TIME STEP : 1.
NUMBER OF TIME STEPS : 10000
GRAPHIC PRINTOUT PERIOD : 50
LISTING PRINTOUT PERIOD : 50

En outre, la simulation a été ré-exécutée avec le drapeau -s, qui écrit la liste complète dans un fichier nommé [FILE-NAME].cas_YEAR-MM-DD-HHhMMminSSs.sortie dans le répertoire de simulation:

telemac2d.py steady2d-conv.cas -s

La direction .cas et les fichiers .sortie peuvent être téléchargés à partir des dépôts hydro-informatiques.com:

• téléchargez stable2d-conv.cas

• télécharger stable2d-conv.cas 2023-07-26-18h41min26s.sortie

Extraire et vérifier les données Flux¶

Les modèles de carnets TELEMAC Jupyter (HOMETEL/notebooks/ > data manip/extraction/.ipynb* ou workshops/exo fluxes.ipynb) fournissent des indications pour extraire les données des résultats de simulation; toutefois, les modèles ne constituent pas un cadre directement applicable pour évaluer la convergence de masse aux limites en fonction de NUMBER OF TIME STEPS. À cette fin, hydro-informatique.com maintient le paquet Python léger pythomac (version ≥ 3.0.0 est décrit ici). Le paquet ne nécessite que numpy, pandas et matplotlib (see the Python installation guide) et court en dehors de l’environnement Python TELEMAC. Deux options d’installation sont disponibles:

pip install pythomac

git clone https://github.com/hydro-informatics/pythomac.git
pip install -e pythomac

La fonction centrale est pythomac.extract_fluxes(). Il localise la liste la plus récente .sortie à côté du fichier de pilotage, analyse le solde de volume et le flux signé imprimé pour chaque limite de liquide à chaque liste imprimée (les formats classiques THERE IS n LIQUID BOUNDARIES et TELEMAC v9 NUMBER OF LIQUID BOUNDARIES: listing sont reconnus), et écrit dans le répertoire de simulation:

• extracted-fluxes.csv - la série chronologique du volume dans le domaine et du flux à travers chaque limite liquide;

• flux-convergence.png - un tracé des magnitudes du flux au cours du temps de simulation (facultatif, plotting=True).

The function returns the extracted series as a pandas.DataFrame indexed by simulation time; the working directory of the calling process is not modified. The implementation can be inspected in flux_analyst.py on GitHub, and the complete API documentation is available at https://pythomac.readthedocs.io.

Pour appliquer la fonction, copiez le code suivant dans un nouveau script Python appelé, par exemple, example_flux_convergence.py, situé dans le répertoire où la simulation stable2d a été effectuée (ou téléchargez exemple flux convergence.py):

# example_flux_convergence.py

from pathlib import Path
from pythomac import extract_fluxes

simulation_dir = str(Path(__file__).parents[1])
telemac_cas = "steady2d.cas"

fluxes_df = extract_fluxes(
    model_directory=simulation_dir,
    cas_name=telemac_cas,
    plotting=True
)

Exécutez le script Python depuis un terminal (ou Anaconda Prompt) dans le répertoire de simulation :

python example_flux_convergence.py

Le script place dans le dossier de simulation :

• le fichier CSV extracted-fluxes.csv (télécharger); et

• la courbe de convergence des flux (flux-convergence.png) à travers les limites du modèle (voir Fig. 1), ce qui indique qualitativement que les flux ont approché la convergence après environ 6000-7000 étapes.

Identifier la convergence¶

Pour déterminer si et quand les flux limites convergent, le déséquilibre relatif des flux est évalué à chaque heure d’impression $t$ comme suit:

où $Q_{i,t}$ et $Q_{j,t}$ = les flux sortants et entrants franchissent les limites du modèle à l’heure $t$, respectivement. Les magnitudes de flux $|\cdot|$ sont nécessaires parce que TELEMAC signale les flux limites avec une convention de signe (flux positif, flux négatif); le bilan massique correspond donc à $|Q_{i,t}| = |Q_{j,t}|$, de sorte que $\varepsilon_{t} \to 0$ à la convergence, et la normalisation par l’afflux $|Q_{j,t}|$ rendes $\varepsilon_{t}$ sans dimension. Dans une simulation stable et régulière, le rapport des déséquilibres de flux consécutifs approche une constante de convergence $c_{\varepsilon}$ égale à l’unité avec un temps croissant:

La combinaison du taux de convergence (ou ordre) $\iota$ et de la constante de convergence $c_{\varepsilon}$ indique:

• convergence linéaire si $\iota$ = 1 ** et** $c_{\varepsilon} \in ]0, 1[$;

• convergence lente sublinéaire si $\iota$ = 1 ** et** $c_{\varepsilon}$ = 1;

• rapide convergence superlinéaire si $\iota$ > 1 et $c_{\varepsilon} \in ]0, 1]$; et

• divergence si $\iota$ = 1 ** et** $c_{\varepsilon}$ > 1, ou $\iota$ < 1.

Le moment où une simulation stable peut être considérée comme ayant atteint un état stable est identifié par le début de la convergence sublinéaire ($\iota$ = 1 et $c_{\varepsilon}$ = 1); c’est-à-dire le moment $t$ au-delà duquel chaque étape supplémentaire $t+1$ améliore la précision du modèle seulement de façon insignifiante (le terme insignificative est quantifié dans le section below). Dans l’hypothèse où le modèle converge sous une forme ou une autre, régler $c_{\varepsilon}$ = 1 rendement $\iota(t)$ en fonction de $\varepsilon_{t}$ et $\varepsilon_{t+1}$:

\début{align} \label{estimation convergence} \frac{\varepsilon {t+1}}{\varepsilon^{\iota(t)}} {t}} &=c {\varepsilon} & \Leftrightarrow \ \iota(t) &= \frac{1}{c {\varepsilon}} \cdot \log {\varepsilon {t}}\varepsilon {t+1} & \overbrace{\Longleftrightarrow}^{c {\varepsilon} = 1}\ \iota(t) &= \log {\varepsilon {t}}\varepsilon {t+1} & \end{align}

These relations are implemented in the pythomac.calculate_convergence() function, which returns a pandas.DataFrame with the columns "Relative imbalance" ($\varepsilon_{t+1}$, Equation (1)) and "Convergence rate" ($\iota(t)$), indexed by simulation time. Its core reads:

import numpy as np
import pandas as pd


def calculate_convergence(series_1, series_2, conv_constant=1.):
    # relative flux imbalance epsilon_t = ||Q_in| - |Q_out|| / |Q_in|; the magnitudes |.|
    #   are needed because Telemac reports outflow negative, so that balance -> epsilon -> 0
    epsilon = np.abs(np.abs(series_1) - np.abs(series_2)) / np.abs(series_1)
    # derive epsilon at t and t+1
    epsilon_t0 = epsilon[:-1]  # cut off last element
    epsilon_t1 = epsilon[1:]   # cut off element zero
    # return the relative imbalance and the convergence rate iota as a pandas DataFrame
    return pd.DataFrame({
        "Relative imbalance": epsilon_t1,
        "Convergence rate": np.emath.logn(epsilon_t0, epsilon_t1) / conv_constant,
    })

Pour calculer $\iota(t)$ (nom variable Python: iota_t) avec la fonction ci-dessus, modifier le exemple flux convergence.py script Python comme suit:

# example_flux_convergence.py

# ...
# add to header:
from pythomac import calculate_convergence

# calculate fluxes_df (see above code block)
fluxes_df = [...]

# back-calculate the printout spacing (in simulation seconds) from the flux index
timestep_in_cas = int(max(fluxes_df.index.values) / (len(fluxes_df.index.values) - 1))

# calculate iota (t) with the calculate_convergence function
iota_t = calculate_convergence(
    series_1=fluxes_df["Fluxes Boundary 1"][1:],  # remove first zero-entry
    series_2=fluxes_df["Fluxes Boundary 2"][1:],  # remove first zero-entry
    cas_timestep=timestep_in_cas,
    plot_dir=simulation_dir,
)

Le taux de convergence résultant $\iota(t)$ est tracé à Fig. 2 pour le steady 2d tutorial avec les périodes d’impression modifiées de 50 secondes et un temps de simulation total de 10000 secondes.

Temps de simulation optimal dérivé¶

Pour économiser le temps de calcul, l’étape dans laquelle les flux d’entrée et de sortie ont convergé présente un intérêt pratique. Les flux tracés dans Fig. 1 et le taux de convergence dans Fig. 2 suggèrent qualitativement que la simulation s’est stabilisée après environ 6000 secondes (étapes de temps). L’extrémité locale dans les deux chiffres près de 4000 étapes dans le temps marque l’interaction des fronts mouillants se propageant à partir des limites amont et aval (voir le animation in the steady 2d tutorial); la convergence monotonique se fixe seulement par la suite.

Parce qu’un jugement purement visuel de convergence est subjectif, un critère objectif est adopté : la longueur optimale de simulation est le plus petit temps $t$ au-delà duquel le déséquilibre relatif du flux $\varepsilon_{t}$ (Equation (1)) reste en permanence en dessous d’une tolérance cible $\varepsilon_{tar}$. Les tolérances de $\varepsilon_{tar}$ = 10$^{-4}$ sont généralement acceptables pour les essais d’étalonnage préliminaires, alors que la validation et les essais d’initialisation à chaud justifient des valeurs plus petites (10$^{-6}$ ou moins). Comme l’illustre Fig. 2, le déséquilibre peut tomber temporairement sous la tolérance et remonter à nouveau (ici près de 4000 pas de temps, lorsque le front en amont franchit la limite en aval); seul le dernier passage permanent est pertinent. L’implémentation algorithmique détecte donc la dernière fois à laquelle $\varepsilon_{t} \geq \varepsilon_{tar}$ et désigne l’impression subséquente comme temps de convergence. Ce critère est implémenté à pythomac.get_convergence_time(), qui retourne l’index d’impression du passage à niveau permanent, ou numpy.nan (avec un avertissement) si la tolérance n’est jamais maintenue. Modifier le script exemple flux convergence.py comme suit:

# example_flux_convergence.py

# ...
# add to header:
from pythomac import get_convergence_time

# calculate fluxes_df and iota_t (see above code blocks)
fluxes_df = [...]
iota_t = [...]

# identify the printout index from which the relative flux imbalance stays
# permanently below the target tolerance (epsilon_tar)
convergence_time_iteration = get_convergence_time(
    relative_imbalance=iota_t["Relative imbalance"],
    convergence_precision=1.0E-4
)

if not str(convergence_time_iteration).lower() == "nan":
    print("The simulation converged after {0} simulation seconds ({1}th printout).".format(
            str(timestep_in_cas * convergence_time_iteration), str(convergence_time_iteration)))

The simulation converged after 6000 simulation seconds (120th printout).

Avec le temps de convergence établi, le mot clé NUMBER OF TIME STEPS dans le fichier de pilotage .cas peut être réduit en conséquence, par exemple:

/ steady2d-conv.cas
TIME STEP : 1.
NUMBER OF TIME STEPS : 6000
GRAPHIC PRINTOUT PERIOD : 50
LISTING PRINTOUT PERIOD : 50

Dépannage des instabilités et des divergences¶

Si une simulation régulière ne parvient pas à obtenir des flux stables, ou si les flux divergent, vérifiez que toutes les limites sont bien définies selon la section des projecteurs sur boundary conditions, et consultez le workflow dans la section sur mass conservation.