Discontinuous Galerkin-Methoden-Tutorial: Ableitung, Flüsse und Python-Implementierung

Die diskontinuierliche Galerkin-Methode, die normalerweise auf DG verkürzt wird, ist eine numerische Technik zur Lösung von partiellen Differentialgleichungen. Es kombiniert mehrere nützliche Eigenschaften von Finite-Volumen- und Finite-Elemente-Methoden. Wie bei der Finite-Volumen-Methode kann die DG Mengen lokal durch sorgfältig definierte Flussmittel konservieren. Wie bei der Finite-Elemente-Methode unterstützt es flexible Netze und polynomische Approximationen in jedem Element.

Das definierende Merkmal von DG ist, dass die numerische Lösung nicht über Elementgrenzen hinweg kontinuierlich bleiben muss. Jedes Element hat seine eigene Polynomdarstellung. Benachbarte Elemente kommunizieren über numerische Flüsse, die an ihren gemeinsamen Schnittstellen bewertet werden.

Diese Struktur macht DG nützlich für Wellenausbreitung, Transportgleichungen, komprimierbare Strömung, Flachwassersysteme, Elektromagnete und Probleme, die Schocks oder Materialgrenzflächen enthalten. Es ist auch attraktiv für das parallele Rechnen, da viele Berechnungen in jedem Element unabhängig durchgeführt werden können.

Wie DG sich von anderen Methoden unterscheidet

Eine Finite-Volumen-Methode speichert normalerweise einen Mittelwert in jeder Zelle und berechnet den Fluss einer konservierten Menge durch Zellflächen. Es bietet eine starke lokale Erhaltung, aber viele grundlegende endliche Volumenschemata verwenden räumliche Approximationen mit relativ niedriger Ordnung.

Eine kontinuierliche Finite-Elemente-Methode stellt die Lösung mit polynomialen Basisfunktionen dar, die sich kontinuierlich über Elementgrenzen verbinden. Dies bietet geometrische Flexibilität und Approximation höherer Ordnung, aber Diskontinuitäten können ohne zusätzliche Behandlung nicht direkt dargestellt werden.

DG verwendet in jedem Element Polynom-Basisfunktionen und erlaubt gleichzeitig separate Werte auf beiden Seiten einer Schnittstelle. Numerische Flüsse bestimmen, wie diese Werte interagieren. Das Verfahren bietet daher lokale Konservierung, Approximation hoher Ordnung und direkte Unterstützung für diskontinuierliche Lösungen.

die kontinuierliche Advektionsgleichung

Die eindimensionale lineare Advektionsgleichung bietet eine klare Einführung in die DG:

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

Für konstante Transportgeschwindigkeit a ist der physikalische Fluss:

f(u) = a u

Die Gleichung wird dann:

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

Dieses Modell beschreibt ein Profil, das sich durch die Domäne bewegt, ohne die Form zu ändern. Wenn a positiv ist, wandern die Informationen von links nach rechts. Wenn es negativ ist, werden Informationen von rechts nach links übertragen.

Unterteilen der Domäne in Elemente

Angenommen, die Domäne erstreckt sich von null bis L. DG teilt es in nicht überlappende Elemente:

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

Innerhalb jedes Elements wird die numerische Lösung durch ein Polynom dargestellt:

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

Die Funktionen φᵢ sind lokale Basisfunktionen, und die Koeffizienten Uᵢ sind die Freiheitsgrade, die sich mit der Zeit ändern.

Ein Polynom von Grad Null speichert in jedem Element einen konstanten Wert. Eine lineare Approximation verwendet zwei lokale Freiheitsgrade in einer Dimension. Eine quadratische Approximation verwendet drei.

Im Gegensatz zu kontinuierlichen finiten Elementen zwingt DG das Polynom von einem Element nicht, um das Polynom im nächsten übereinzustimmen. An einer gemeinsamen Schnittstelle kann die Lösung daher einen linken Wert und einen rechten Wert haben.

Ableiten der elementweise schwachen Form

Die GD-Formulierung beginnt mit der Multiplikation der maßgeblichen Gleichung mit einer Testfunktion v und der Integration über ein Element:

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

Die Integration nach Teilen verschiebt die räumliche Ableitung vom physikalischen Fluss zur Testfunktion:

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

Der Grenzbegriff ist wesentlich. Es beschreibt den Fluss durch die linke und rechte Seite des Elements.

Da die DG-Lösung an einer Schnittstelle diskontinuierlich sein kann, ist der physikalische Fluss dort nicht eindeutig definiert. Ein Element liefert einen Wert von links, während sein Nachbar einen anderen Wert von rechts angibt. DG ersetzt den mehrdeutigen physikalischen Fluss durch einen numerischen Fluss:

f̂(u⁻, u⁺)

Die letzte schwache Form für ein Element ist:

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

Diese Gleichung enthält einen element-lokalen Volumenbeitrag und zwei Schnittstellenbeiträge. Die numerischen Flüsse sind die einzigen Begriffe, die benachbarte Elemente direkt verbinden.

Warum numerische Flüsse wichtig sind

Ein numerischer Fluss muss konsistent sein. Wenn die beiden Schnittstellenwerte gleich sind, sollte der physikalische Fluss reproduziert werden:

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

Es muss auch eine angemessene Stabilität bieten. Bei Transportproblemen sollte der Fluss die Richtung berücksichtigen, in die sich Informationen bewegen.

Die Flusswahl steuert, wie stark benachbarte Elemente interagieren und wie viel numerische Dissipation eingeführt wird. Ein ungeeigneter Fluss kann Schwingungen, übermäßige Glättung oder falsche Wellengeschwindigkeiten verursachen.

Aufwind-Fluss

Für die lineare Advektion mit positiver Geschwindigkeit kommt die Information von der linken Seite einer Schnittstelle. Der Aufwindfluss ist daher:

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

Für die negative Geschwindigkeit kommen Informationen von rechts:

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

Ein kompakter Ausdruck ist:

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

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

Der Aufwindfluss ist einfach, stabil und weit verbreitet. Es führt eine gewisse numerische Dissipation ein, aber die Menge nimmt im Allgemeinen ab, wenn das Netz verfeinert wird oder die Polynomreihenfolge zunimmt.

Lokaler LAX-FRIEDRICHS FLUX

Der lokale Lax-Friedrichs-Fluss, auch Rusanov-Fluss genannt, wird häufig für nichtlineare Erhaltungsgesetze verwendet:

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

Der Parameter α ist eine Schätzung der größten Kennzahl an der Schnittstelle. Für konstante lineare Advektion ist es normalerweise |a| .

Der erste Teil mittelt die physikalischen Flüsse. Der zweite Teil fügt eine Dissipation hinzu, die hilft, instabile Schnittstellenmodi zu steuern. Dieser Fluss ist robust und einfach zu implementieren, kann jedoch diffusiver sein als ein spezialisierter ungefährer Riemann-Solver.

Riemann-Solver

Nichtlineare hyperbolische Systeme können mehrere Wellen enthalten, die sich mit unterschiedlichen Geschwindigkeiten bewegen. Beispiele sind die Euler-Gleichungen für die komprimierbare Strömung und die Flachwassergleichungen.

Ein Riemann-Solver untersucht die linken und rechten Zustände an einer Grenzfläche und schätzt die Wellen, die durch ihre Wechselwirkung erzeugt werden. Exakte Riemann-Solver können teuer sein, daher verwenden praktische DG-Codes üblicherweise ungefähre Methoden wie Roe, HLL, HLLC oder Rusanov.

Die richtige Wahl hängt von der Gleichung, der gewünschten Genauigkeit, den Robustheitsanforderungen und der Fähigkeit ab, wichtige physikalische Eigenschaften beizubehalten.

Das halbdiskrete Matrixsystem

Nach der Auswahl von Basis- und Testfunktionen kann die schwache Form als System gewöhnlicher Differentialgleichungen ausgedrückt werden:

M dU/dt = R(U)

Die Massenmatrix enthält Integrale von Basisfunktionsprodukten:

Mᵢⱼ = ∫K φᵢ φⱼ dx

Das Residuum R(U) umfasst die Volumenableitung und den numerischen Flussbeitrag.

Da die DG-Basisfunktionen zu einzelnen Elementen gehören, hat die globale Massenmatrix eine Blockdiagonalstruktur. Jedes Element trägt einen kleinen unabhängigen Block bei. Diese Blöcke können separat invertiert werden.

Die Massenmatrix ist nicht automatisch diagonal für jede Basis- und Integrationsregel. Es wird diagonal oder annähernd diagonal in gängigen Knoten-DG-Formulierungen, die passende Interpolations- und Quadraturpunkte wie Gauss-Lobatto-Kollokation verwenden. Andere Formulierungen verwenden kleine dichte Elementmatrizen.

Auswahl von Basisfunktionen und Quadratur

DG-Methoden verwenden üblicherweise modale oder nodale Polynombasen. Eine modale Basis stellt die Lösung durch Polynommodi dar, die häufig auf Legendre-Polynomen basieren. Eine Knotenbasis speichert Lösungswerte an Interpolationsknoten im Element.

Knotenformulierungen sind praktisch, da Schnittstellenwerte direkt erhalten werden können, wenn Knoten an Elementgrenzen platziert werden. Gauss-Lobatto-Punkte umfassen beide Endpunkte des Referenzelements, während Gauß-Punkte darin bleiben.

Die numerische Quadratur bewertet die Integrale in der schwachen Form. Die Quadraturregel muss für den Polynomgrad und alle nichtlinearen Terme ausreichend genau sein. Unzureichende Quadratur kann zu Aliasing-Fehlern und Instabilität führen.

Polynomreihenfolge und Genauigkeit

Für ausreichend reibungslose Lösungen kann ein gut ausgelegtes DG-Verfahren einen Fehler proportional zu ungefähr h^(p+1) erzielen. Eine Erhöhung des Polynomgrads kann daher die Genauigkeit verbessern, ohne mehr Elemente hinzuzufügen.

Höhere Ordnung ist nicht immer besser. Schocks und scharfe Diskontinuitäten können Schwingungen in der Nähe des Sprunges verursachen. Begrenzer, künstliche Viskosität, Filterung oder Schockerfassungstechniken können erforderlich sein.

Ein linearer DG-Solver in Numpy

Das folgende Lehrbeispiel implementiert eine lineare p = 1-DG-Methode für die periodische eindimensionale Advektion. Jedes Element enthält zwei Freiheitsgrade an seinen Endpunkten.

Die Implementierung verwendet das Element Schwache Form, einen Aufwindfluss und eine dritte Ordnung, die stabilitätserhaltende Runge-Kutta-Methode.

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())

Dieses Beispiel ist absichtlich auf ein glattes lineares Problem beschränkt. Produktions-DG-Solver benötigen außerdem robuste Grenzbehandlungen, mehrdimensionale Zuordnungen, genaue Quadratur, nichtlineare Flussmittelbewertung, Begrenzer und eine erweiterte Zeitintegration.

Zeitintegration und die CFL-Bedingung

Die räumliche DG-Diskretisierung schafft ein System gewöhnlicher Differentialgleichungen. Eine explizite Methode wie Runge-Kutta kann dieses System dann rechtzeitig voranbringen.

Der maximale stabile Zeitschritt hängt von der Wellengeschwindigkeit, der Elementgröße, dem Polynomgrad, dem Fluss und dem Zeitintegrator ab. Eine gängige Skalierung ist:

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

Die genaue Stabilitätskonstante ist methodenabhängig. Es sollte durch Analyse, Dokumentation oder numerische Tests festgelegt und nicht als universell behandelt werden.

Das Erhöhen der Polynomreihenfolge reduziert normalerweise den größten expliziten Zeitschritt. Hochwertige DG erfordert daher möglicherweise mehr Schritte, auch wenn sie weniger Elemente benötigt.

DG für advektionsdominierte Probleme

Standard-kontinuierliche Galerkin-Methoden können Schwingungen entwickeln, wenn die Advektion die Diffusion dominiert. Stabilisierte kontinuierliche Methoden wie supg modifizieren die Testfunktionen, um die Kontrolle entlang von Streamlines hinzuzufügen.

DG übernimmt den Transport durch Schnittstellenflüsse. Aufwärtswind oder ungefähre Riemann-Flüsse führen richtungsbewusste Stabilisierung ein und bewahren die lokale Erhaltung.

SUPG bleibt bei vielen reibungslosen Konvektionsproblemen wirksam. DG ist attraktiv, wenn lokale Konservierung, diskontinuierliche Lösungen, komplexe Netze oder elementbezogene Adaptivität zentrale Anforderungen sind.

DG vs. FVM und kontinuierliche FEM

Vorteile der DG-Methode

• Die lokale Erhaltung ist in die Schnittstellenflussbilanz integriert.
• Der Polynomgrad kann zwischen den Elementen variieren.
• Diskontinuitäten verletzen nicht den Approximationsraum.
• Elementberechnungen eignen sich gut für parallele Hardware.
• Komplexe und unstrukturierte Netze können unterstützt werden.
• Netzverfeinerung und Polynomanreicherung können kombiniert werden.
• Für verschiedene Gleichungen können unterschiedliche numerische Flüsse ausgewählt werden.

Einschränkungen der DG

DG verwendet normalerweise mehr Freiheitsgrade als eine kontinuierliche Finite-Elemente-Methode mit demselben Polynomgrad, da benachbarte Elemente keine Schnittstellenwerte teilen.

Die Methode erfordert auch eine sorgfältige Flussentwürfe. Hyperbolische, elliptische und gemischte Gleichungen erfordern unterschiedliche Grenzflächenbehandlungen. Diffusionsoperatoren benötigen Formulierungen wie Innenstrafe, lokale DG oder verwandte Techniken.

Explizite DG-Methoden höherer Ordnung können restriktive Zeitschrittbeschränkungen haben. Nichtlineare Probleme können Begrenzer, entropiestabile Flüsse, Positivitätserhaltung oder künstliche Viskosität erfordern.

Diese Funktionen machen die DG leistungsfähig, aber schwieriger, sie richtig zu implementieren als grundlegende Finite-Differenz- oder Finite-Volumen-Schemata.

Häufige Implementierungsfehler

Ein häufiger Fehler besteht darin, den physischen Fluss direkt an einer diskontinuierlichen Schnittstelle zu verwenden, ohne zu definieren, wie die beiden benachbarten Werte kombiniert werden sollen.

Andere häufige Fehler sind:

• Verwenden des falschen Vorzeichens für den linken oder rechten Grenzbeitrag
• Anwenden inkonsistenter Normal-Vektor-Orientierungen
• Angenommen, die Massenmatrix ist immer diagonal
• Verwenden von Quadratur, die für nichtlineare Begriffe zu schwach ist
• Ignorieren der Polynom-Grad-Abhängigkeit der CFL-Grenze
• periodische oder physikalische Randbedingungen falsch anwenden
• Verwenden höherwertiger Polynome in der Nähe von Schocks ohne Begrenzer
• Testen nur visuelle Ausgabe anstelle von Konvergenz und Konservierung

So validieren Sie einen DG-Solver

Beginnen Sie mit einem Problem, das eine bekannte analytische Lösung hat. Eine periodische lineare Advektion ist nützlich, da das genaue Profil einfach um a t verschoben wird.

Führen Sie den Solver mit mehreren Netzauflösungen aus und messen Sie eine Fehlernorm. Für eine glatte Lösung sollte sich die beobachtete Konvergenzrate der theoretischen Reihenfolge der Methode nähern.

Überprüfen Sie die Konservierung durch Integration der Lösung über die Domäne. Für die periodische Advektion sollte die Gesamtmasse nahezu konstant bleiben.

Testen Sie sowohl positive als auch negative Geschwindigkeiten, um die Aufwindrichtung zu überprüfen. Konstante Anfangsdaten sollten konstant bleiben. Grenzzustände sollten auch separat getestet werden, bevor sie zu nichtlinearen Gleichungen übergehen.

Wann wählen Sie DG

DG ist eine gute Wahl, wenn das Problem Wellen, Schocks, Materialschnittstellen oder transportdominiertes Verhalten enthält. Es ist auch nützlich, wenn die lokale Konservierung unerlässlich ist oder wenn die Simulation von der Genauigkeit hoher Ordnung auf einem unstrukturierten Netz profitiert.

Eine grundlegende Methode des endlichen Volumens kann für Konservierungsprobleme niedriger Ordnung einfacher bleiben. Kontinuierliche endliche Elemente können für glatte elliptische oder strukturelle Probleme effizienter sein, wenn keine Unterbrechungen erwartet werden.

Die numerische Methode sollte der mathematischen Struktur der PDE folgen und nicht die aktuelle Popularität oder die Verfügbarkeit der Software.

Schlussfolgerung

Die diskontinuierliche Galerkin-Methode stellt die Lösung mit unabhängigen Polynomen in jedem Element dar. Seine schwache Formulierung erzeugt natürliche Grenzbegriffe, die durch numerische Flüsse ersetzt werden. Diese Flüsse steuern, wie sich Informationen zwischen Elementen bewegen und die vollständige Methode, um die Mengen lokal zu erhalten.

DG kombiniert Polynom-Approximation höherer Ordnung, geometrische Flexibilität, Diskontinuitätshandhabung und starke Elementlokalität. Diese Vorteile machen es wertvoll für Wellenausbreitung, hyperbolische Erhaltungsgesetze, konvektionsdominierte Systeme und große parallele Simulationen.

Das Verfahren führt auch zusätzliche Komplexität ein. Flüsse, Quadratur, Zeitschrittbeschränkungen, Begrenzer und Randbedingungen müssen sorgfältig ausgewählt werden. Eine zuverlässige Implementierung sollte schrittweise entwickelt werden, beginnend mit einer einfachen linearen Gleichung und durch Konvergenz- und Konservierungstests verifiziert werden.

Für Forscher, die bereits endliche Volumen- oder Finite-Elemente-Methoden verstehen, bietet DG einen natürlichen nächsten Schritt zu lokal konservativen PDE-Solvern in wissenschaftlicher Python.