Tutoriel discontinu de la méthode Galerkin : dérivation, flux et implémentation de Python

La méthode de Galerkin discontinue, généralement raccourcie en DG, est une technique numérique pour résoudre les équations aux dérivées partielles. Il combine plusieurs propriétés utiles des méthodes de volume fini et d’éléments finis. Comme la méthode du volume fini, la DG peut préserver les quantités localement grâce à des flux soigneusement définis. Comme la méthode des éléments finis, il prend en charge les mailles flexibles et les approximations polynomiales à l’intérieur de chaque élément.

La caractéristique déterminante de DG est que la solution numérique n’a pas besoin de rester continue au-delà des limites des éléments. Chaque élément a sa propre représentation polynomiale. Les éléments voisins communiquent à travers des flux numériques évalués à leurs interfaces partagées.

Cette structure rend la DG utile pour la propagation des ondes, les équations de transport, le flux compressible, les systèmes d’eau peu profonde, l’électromagnétique et les problèmes qui contiennent des chocs ou des interfaces de matériaux. Il est également attrayant pour le calcul parallèle, car de nombreux calculs peuvent être effectués indépendamment à l’intérieur de chaque élément.

Comment la DG diffère des autres méthodes

Une méthode de volume fini stocke généralement une valeur moyenne dans chaque cellule et calcule le flux d’une quantité conservée à travers les faces des cellules. Il fournit une forte conservation locale, mais de nombreux schémas de base de volumes finis utilisent des approximations spatiales d’ordre relativement faible.

Une méthode d’éléments finis continus représente la solution avec des fonctions de base polynomiales qui se connectent en continu au-delà des limites des éléments. Cela offre une flexibilité géométrique et une approximation d’ordre élevé, mais les discontinuités ne peuvent être représentées directement sans traitement supplémentaire.

DG utilise des fonctions de base polynomiales à l’intérieur de chaque élément tout en permettant des valeurs séparées de chaque côté d’une interface. Les flux numériques déterminent la façon dont ces valeurs interagissent. La méthode fournit donc une conservation locale, une approximation d’ordre élevé et un soutien direct pour des solutions discontinues.

L’équation d’advection continue

L’équation d’advection linéaire unidimensionnelle fournit une introduction claire à DG :

∂u/∂t + ∂f(u)/∂x = 0

Pour une vitesse de transport constante a, le flux physique est :

f(u) = a u

L’équation devient alors :

∂u/∂t + a ∂u/∂x = 0

Ce modèle décrit un profil qui se déplace dans le domaine sans changer de forme. Lorsque a est positif, l’information se déplace de gauche à droite. Lorsqu’elle est négative, l’information se déplace de droite à gauche.

Diviser le domaine en éléments

Supposons que le domaine s’étende de zéro à L. La DG la divise en éléments non chevauchants :

Ω = K₁ ∪ K₂ ∪ ... ∪ Kₙ

Au sein de chaque élément, la solution numérique est représentée par un polynôme :

uₕ(x, t) = Σ Uᵢ(t) φᵢ(x)

Les fonctions φᵢ sont des fonctions de base locale et les coefficients Uᵢ sont les degrés de liberté qui changent avec le temps.

Un polynôme de degré zéro stocke une valeur constante dans chaque élément. Une approximation linéaire utilise deux degrés de liberté locaux dans une dimension. Une approximation quadratique en utilise trois.

Contrairement aux éléments finis continus, DG ne force pas le polynôme d’un élément à correspondre au polynôme dans le suivant. À une interface partagée, la solution peut donc avoir une valeur gauche et une valeur correcte.

Dériver la forme faiblement élémentaire

La formulation DG commence par multiplier l’équation régissant par une fonction de test v et intégrant sur un élément :

∫K (∂uₕ/∂t) v dx + ∫K (∂f(uₕ)/∂x) v dx = 0

L’intégration par parties déplace la dérivée spatiale du flux physique vers la fonction de test :

∫K (∂uₕ/∂t) v dx
- ∫K f(uₕ) ∂v/∂x dx
+ f(uₕ) v |∂K = 0

Le terme limite est essentiel. Il décrit le flux à travers les côtés gauche et droit de l’élément.

Étant donné que la solution DG peut être discontinue à une interface, le flux physique n’y est pas défini de manière unique. Un élément fournit une valeur à partir de la gauche, tandis que son voisin fournit une autre valeur à partir de la droite. DG remplace le flux physique ambigu par un flux numérique :

f̂(u⁻, u⁺)

La forme finale faible d’un élément est :

∫K (∂uₕ/∂t) v dx
- ∫K f(uₕ) ∂v/∂x dx
+ f̂R vR
- f̂L vL = 0

Cette équation contient une contribution élément-local de volume et deux contributions d’interface. Les flux numériques sont les seuls termes qui relient directement les éléments voisins.

Pourquoi les flux numériques sont importants

Un flux numérique doit être cohérent. Lorsque les deux valeurs d’interface sont égales, il doit reproduire le flux physique :

f̂(u, u) = f(u)

Il doit également assurer une stabilité appropriée. Pour les problèmes de transport, le flux doit respecter la direction dans laquelle les informations se déplacent.

Le choix de flux contrôle la force avec laquelle les éléments voisins interagissent et la quantité de dissipation numérique est introduite. Un flux inadapté peut créer des oscillations, un lissage excessif ou des vitesses d’onde incorrectes.

Flux au vent

Pour l’advection linéaire avec une vitesse positive, l’information provient du côté gauche d’une interface. Le flux au vent est donc :

f̂(u⁻, u⁺) = a u⁻    when a > 0

Pour la vitesse négative, les informations proviennent de la droite :

f̂(u⁻, u⁺) = a u⁺    when a < 0

Une expression compacte est :

f̂ = a⁺u⁻ + a⁻u⁺

a⁺ = max(a, 0)
a⁻ = min(a, 0)

Le flux au vent est simple, stable et largement utilisé. Il introduit une certaine dissipation numérique, mais la quantité diminue généralement à mesure que le maillage est raffiné ou que l’ordre polynomial augmente.

Flux local de Lax-Friedrichs

Le flux local de Lax-Friedrichs, également appelé flux de Rusanov, est fréquemment utilisé pour les lois de conservation non linéaires :

f̂(u⁻, u⁺)
= 0.5[f(u⁻) + f(u⁺)]
- 0.5 α(u⁺ - u⁻)

Le paramètre α est une estimation de la plus grande vitesse caractéristique à l’interface. Pour une advection linéaire constante, c’est généralement |a|.

La première partie fait la moyenne des flux physiques. La deuxième partie ajoute une dissipation qui aide à contrôler les modes d’interface instables. Ce flux est robuste et facile à mettre en œuvre, bien qu’il puisse être plus diffusant qu’un solveur Riemann approximatif spécialisé.

Solveurs de Riemann

Les systèmes hyperboliques non linéaires peuvent contenir plusieurs ondes se déplaçant à des vitesses différentes. Les exemples incluent les équations d’Euler pour l’écoulement compressible et les équations d’eau peu profonde.

Un solveur de Riemann examine les états gauche et droit à une interface et estime les ondes produites par leur interaction. Les solveurs exacts de Riemann peuvent être coûteux, de sorte que les codes DG pratiques utilisent généralement des méthodes approximatives telles que ROE, HLL, HLLC ou Rusanov.

Le choix correct dépend de l’équation, de la précision souhaitée, des exigences de robustesse et de la capacité de préserver des propriétés physiques importantes.

Le système à matrice semi-discrète

Après avoir sélectionné les fonctions de base et de test, la forme faible peut être exprimée sous la forme d’un système d’équations différentielles ordinaires :

M dU/dt = R(U)

La matrice de masse contient les intégrales des produits de base-fonction :

Mᵢⱼ = ∫K φᵢ φⱼ dx

Le résidu R(U) comprend les contributions de dérivés de volume et de flux numérique.

Étant donné que les fonctions de base DG appartiennent à des éléments individuels, la matrice de masse globale a une structure en bloc-diagonale. Chaque élément contribue à un petit bloc indépendant. Ces blocs peuvent être inversés séparément.

La matrice de masse n’est pas automatiquement diagonale pour chaque base et règle d’intégration. Il devient diagonal ou approximativement diagonal dans les formulations de DG nodales courantes qui utilisent des points d’interpolation et de quadrature correspondants, tels que la collocation de Gauss-Lobatto. D’autres formulations utilisent de petites matrices d’éléments denses.

Choisir les fonctions de base et la quadrature

Les méthodes DG utilisent généralement des bases polynomiales modales ou nodales. Une base modale représente la solution à travers des modes polynomiaux, souvent basés sur des polynômes de Legendre. Une base nodale stocke les valeurs de la solution aux nœuds d’interpolation à l’intérieur de l’élément.

Les formulations nodales sont pratiques car les valeurs d’interface peuvent être obtenues directement lorsque des nœuds sont placés sur les limites des éléments. Les points de Gauss-Lobatto incluent les deux extrémités de l’élément de référence, tandis que les points de Gauss restent à l’intérieur.

La quadrature numérique évalue les intégrales sous la forme faible. La règle de la quadrature doit être suffisamment précise pour le degré polynomial et tous les termes non linéaires. Une quadrature insuffisante peut introduire une erreur et une instabilité d’alias.

Ordre et précision polynomiales

Pour des solutions suffisamment lisses, une méthode DG bien conçue peut atteindre une erreur proportionnelle à environ h^(p+1). L’augmentation du degré polynomial peut donc améliorer la précision sans ajouter plus d’éléments.

L’ordre supérieur n’est pas toujours mieux. Les chocs et les discontinuités nettes peuvent créer des oscillations près du saut. Des limiteurs, une viscosité artificielle, un filtrage ou des techniques de capture de chocs peuvent être nécessaires.

Un solveur DG linéaire dans NumPy

L’exemple éducatif suivant met en œuvre une méthode DG linéaire p = 1 pour une advection unidimensionnelle périodique. Chaque élément contient deux degrés de liberté situés à ses extrémités.

L’implémentation utilise la forme faible de l’élément, un flux au vent et une méthode Runge-Kutta préservant la stabilité de force du troisième ordre.

import numpy as np

# Domain and model parameters
length = 1.0
number_of_elements = 80
velocity = 1.0
final_time = 0.5

element_width = length / number_of_elements
jacobian = element_width / 2.0

# Reference-element mass matrix for linear basis functions
mass_reference = np.array([
    [2.0 / 3.0, 1.0 / 3.0],
    [1.0 / 3.0, 2.0 / 3.0]
])

mass_matrix = jacobian * mass_reference
inverse_mass = np.linalg.inv(mass_matrix)

# S[i, j] = integral(phi_j * derivative(phi_i)) on [-1, 1]
volume_matrix = np.array([
    [-0.5, -0.5],
    [ 0.5,  0.5]
])

left_vector = np.array([1.0, 0.0])
right_vector = np.array([0.0, 1.0])

# Physical coordinates of local DG nodes
left_edges = np.arange(number_of_elements) * element_width
right_edges = left_edges + element_width
coordinates = np.column_stack((left_edges, right_edges))

# Smooth periodic initial condition
solution = (
    0.5
    + 0.5 * np.sin(2.0 * np.pi * coordinates / length)
)

def upwind_flux(left_state, right_state, speed):
    positive_speed = max(speed, 0.0)
    negative_speed = min(speed, 0.0)

    return (
        positive_speed * left_state
        + negative_speed * right_state
    )

def spatial_residual(values):
    residual = np.zeros_like(values)

    for element in range(number_of_elements):
        left_neighbor = (element - 1) % number_of_elements
        right_neighbor = (element + 1) % number_of_elements

        # States on the left interface
        left_inside = values[element, 0]
        left_outside = values[left_neighbor, 1]

        # States on the right interface
        right_inside = values[element, 1]
        right_outside = values[right_neighbor, 0]

        flux_left = upwind_flux(
            left_outside,
            left_inside,
            velocity
        )

        flux_right = upwind_flux(
            right_inside,
            right_outside,
            velocity
        )

        local_rhs = (
            velocity * volume_matrix @ values[element]
            + flux_left * left_vector
            - flux_right * right_vector
        )

        residual[element] = inverse_mass @ local_rhs

    return residual

# Conservative time-step estimate
time_step = 0.1 * element_width / abs(velocity)
current_time = 0.0

while current_time < final_time:
    dt = min(time_step, final_time - current_time)

    # SSP-RK3 stage 1
    stage_one = solution + dt * spatial_residual(solution)

    # SSP-RK3 stage 2
    stage_two = (
        0.75 * solution
        + 0.25 * (
            stage_one
            + dt * spatial_residual(stage_one)
        )
    )

    # SSP-RK3 stage 3
    solution = (
        (1.0 / 3.0) * solution
        + (2.0 / 3.0) * (
            stage_two
            + dt * spatial_residual(stage_two)
        )
    )

    current_time += dt

print("Simulation completed")
print("Final time:", current_time)
print("Minimum value:", solution.min())
print("Maximum value:", solution.max())

Cet exemple est intentionnellement limité à un problème linéaire lisse. Les solveurs de DG de production nécessitent également des traitements aux limites robustes, des mappages multidimensionnels, une quadrature précise, une évaluation de flux non linéaire, des limiteurs et une intégration plus avancée dans le temps.

Intégration du temps et condition LFC

La discrétisation spatiale DG crée un système d’équations différentielles ordinaires. Une méthode explicite telle que Runge-Kutta peut alors faire avancer ce système dans le temps.

Le pas de temps stable maximum dépend de la vitesse d’onde, de la taille des éléments, du degré polynomial, du flux et de l’intégrateur de temps. Une mise à l’échelle commune est :

Δt ∝ h / [(2p + 1)|a|]

La constante de stabilité exacte dépend de la méthode. Il doit être établi par l’analyse, la documentation ou les tests numériques plutôt que traités comme universels.

L’augmentation de l’ordre polynomial réduit normalement le plus grand pas de temps explicitement stable. DG d’ordre élevé peut donc nécessiter plus d’étapes, même lorsqu’il a besoin de moins d’éléments.

DG pour les problèmes dominés par l’advection

Les méthodes de Galerkin continue standard peuvent développer des oscillations lorsque l’advection domine la diffusion. Des méthodes continues stabilisées telles que SUPG modifient les fonctions de test pour ajouter un contrôle le long des lignes de courant.

DG gère le transport à travers les flux d’interface. Au vent ou approximativement, les flux de Riemann introduisent une stabilisation sensible à la direction et préservent la conservation locale.

Le SUPG reste efficace pour de nombreux problèmes dominés par la convection. La DG est attrayante lorsque la conservation locale, les solutions discontinues, les mailles complexes ou l’adaptabilité des éléments sont des exigences centrales.

DG vs FVM et FEM continue

Avantages de la méthode DG

• La conservation locale est intégrée à l’équilibre des flux d’interface.
• Le degré polynomial peut varier d’un élément à l’autre.
• Les discontinuités ne violent pas l’espace d’approximation.
• Les calculs d’éléments sont bien adaptés au matériel parallèle.
• Des maillages complexes et non structurés peuvent être pris en charge.
• Le raffinement du maillage et l’enrichissement polynomial peuvent être combinés.
• Différents flux numériques peuvent être sélectionnés pour différentes équations.

Limites de la DG

DG utilise normalement plus de degrés de liberté qu’une méthode par éléments finis continus du même degré polynomial, car les éléments voisins ne partagent pas les valeurs d’interface.

La méthode nécessite également une conception de flux minutieuse. Les équations hyperboliques, elliptiques et mixtes nécessitent des traitements d’interface différents. Les opérateurs de diffusion exigent des formulations telles que la pénalité intérieure, la DG locale ou des techniques connexes.

Les méthodes explicites de DG d’ordre élevé peuvent avoir des limites de pas de temps restrictives. Les problèmes non linéaires peuvent nécessiter des limiteurs, des flux stables entropie, une préservation de la positivité ou une viscosité artificielle.

Ces fonctionnalités rendent la DG puissante mais plus difficile à implémenter correctement que les schémas de base des différences finies ou des volumes finis.

erreurs de mise en œuvre courantes

Une erreur courante consiste à utiliser le flux physique directement à une interface discontinue sans définir comment les deux valeurs voisines doivent être combinées.

D’autres erreurs fréquentes incluent :

• Utiliser le mauvais signe pour la contribution des limites gauche ou droite
• Appliquer des orientations vectorielles incohérentes
• En supposant que la matrice de masse est toujours diagonale
• Utilisation de la quadrature trop faible pour les termes non linéaires
• Ignorer la dépendance polynomiale-degré de la limite de LFC
• appliquer des conditions aux limites périodiques ou physiques de manière incorrecte
• Utilisation de polynômes d’ordre élevé à proximité de chocs sans limiteur
• Tester uniquement la sortie visuelle au lieu de la convergence et de la conservation

Comment valider un solveur DG

Commencez par un problème qui a une solution analytique connue. L’advection linéaire périodique est utile car le profil exact est simplement décalé de a t.

Exécutez le solveur avec plusieurs résolutions de maillage et mesurez une norme d’erreur. Pour une solution lisse, le taux de convergence observé doit approcher de l’ordre théorique de la méthode.

Vérifiez la conservation en intégrant la solution à travers le domaine. Pour une advection périodique, la masse totale devrait rester presque constante.

Testez les vitesses positives et négatives pour vérifier la direction du vent. Les données initiales constantes doivent rester constantes. Les états limites doivent également être testés séparément avant de passer à des équations non linéaires.

Quand choisir DG

La DG est un choix fort lorsque le problème contient des ondes, des chocs, des interfaces matérielles ou un comportement dominé par les transports. Il est également utile lorsque la conservation locale est essentielle ou lorsque la simulation bénéficie d’une précision d’ordre élevé sur un maillage non structuré.

Une méthode de base du volume fini peut rester plus simple pour les problèmes de conservation de faible ordre. Les éléments finis continus peuvent être plus efficaces pour des problèmes elliptiques ou structurels lisses où des discontinuités ne sont pas attendus.

La méthode numérique doit suivre la structure mathématique de la PDE plutôt que la popularité actuelle ou la disponibilité des logiciels.

Conclusion

La méthode Galerkin discontinue représente la solution avec des polynômes indépendants à l’intérieur de chaque élément. Sa formulation faible produit des termes limites naturels qui sont remplacés par des flux numériques. Ces flux contrôlent la façon dont les informations se déplacent entre les éléments et permettent à la méthode complète de préserver les quantités localement.

DG combine une approximation polynomiale d’ordre élevé, une flexibilité géométrique, une manipulation de discontinuité et une localité d’éléments solide. Ces avantages le rendent utile pour la propagation des ondes, les lois de conservation hyperboliques, les systèmes dominés par la convection et les simulations parallèles de grande taille.

La méthode introduit également une complexité supplémentaire. Les flux, les quadratures, les restrictions de pas de temps, les limiteurs et les conditions aux limites doivent être sélectionnés avec soin. Une mise en œuvre fiable doit être développée progressivement, en commençant par une simple équation linéaire et vérifiée par des tests de convergence et de conservation.

Pour les chercheurs qui comprennent déjà le volume fini ou les méthodes d’éléments finis, la DG fournit une prochaine étape naturelle vers des solveurs PDE localement conservateurs dans le Python scientifique.