1D upwind DG implementation does not work propertly.

[afgalindo]

Hello, I am implementing a simple upwind Discontinuous Galerkin method for the one-dimensional transport equation, with periodic boundary conditions in (x).

$$ u_t + u_x = 0, \quad u(x, 0) = \sin(x) $$

For time integration, I am using the TVD-RK 3 method, which is a third-order TVD (Total Variation Diminishing) Runge-Kutta scheme. This method is applied to the transport equation in a manner similar to what is described in the Firedrake tutorial: DG Advection Demo.

Just for the sake of completeness, each element has the form (I_i = [x_{i-1/2}, x_{i+1/2}]). Taking into account that we are using an upwind flux and periodic boundary conditions, the weak form has the form:

$$ \int_{I} u_t, \varphi dx = \int_{I} u \varphi_x dx + \sum_{i=1}^N u(x_{i+1/2}^-) [\varphi]_{i+1/2} $$

I wrote the implementation below for this problem up to final time (T = 0.1), and I always get the wrong result and even blow ups. I have checked the math and the code many times and I cannot seem to find what is going on. Also in many cases the boundary conditions are not preserved. We have checked this with other colleagues and we were wondering if you guys can check this for us.

Code:

from firedrake import *
import matplotlib.pyplot as plt
import numpy as np
from math import pi

if __name__ == "__main__":
    ####################################################################################
    # Mesh parameters 
    ####################################################################################
    left  = 0.0      # Left e``ndpoint of the interval
    right = 2.0*pi  # Right endpoint of the interval
    length = right - left  # Length of the interval
    Nx    = 100        #Number of physical space elements.
    Delta_x = length/Nx  # Mesh size
    mesh   = PeriodicIntervalMesh(Nx, length)  # Periodic mesh for the interval [left, right] 
    x     = SpatialCoordinate(mesh)  # Spatial coordinate
    ####################################################################################
    # Finite element space
    ####################################################################################
    b  = 0  # polynomial degree
    V  = FunctionSpace(mesh, 'DG', 1)
    ################################################################################### 
    # Time discretization parameters
    ###################################################################################
    T   = 0.01        # Final time
    CFL = 0.01  # Courant-Friedrichs-Lewy condition for stability
    dt = CFL*Delta_x/(2.0*(b+1)+1)  # Time step size as a constant for Firedrake  #In principle this should be a function of the mesh size dx.
    dtc = Constant(dt)  # Time step size as a constant for Firedrake
    ####################################################################################
    # Initial condition
    ####################################################################################
    u = Function(V).interpolate(sin(x[0]))  # Function to hold the solution
    u_init = Function(V).assign(u)  # Initial condition
    ######################################################################################
    # Define the variational problem
    ######################################################################################
    du_trial = TrialFunction(V)     # Trial function du/dt.
    phi      = TestFunction(V)      # Test function
    a        = du_trial * phi * dx  # Left-hand side bilinear form 
    L1 = dtc*(u* phi.dx(0)*dx+(u('-')*jump(phi)*dS))  # Right-hand side linear form
    ######################################################################################
    # Time integration: TVD RK3 method
    ######################################################################################  
    # Assemble the Left Hand side matrix 
    du = Function(V)  # Increment in the solution at each time step
    # Variables for RK 
    u1 = Function(V)  # First stage solution
    u2 = Function(V)  # Second stage solution
    # For Rk3 
    L2 = replace(L1, {u: u1})
    L3 = replace(L1, {u: u2})


    # Parameters for solving 
    params = {'ksp_type': 'preonly', 'pc_type': 'bjacobi', 'sub_pc_type': 'ilu'}
    prob1  = LinearVariationalProblem(a, L1, du)
    solv1  = LinearVariationalSolver(prob1, solver_parameters=params)  # Solver for the first stage
    prob2  = LinearVariationalProblem(a, L2, du)
    solv2  = LinearVariationalSolver(prob2, solver_parameters=params)  # Solver for the second stage
    prob3  = LinearVariationalProblem(a, L3, du)
    solv3  = LinearVariationalSolver(prob3, solver_parameters=params)  # Solver for the third stage
    # Create a time loop
    t=0.0
    ######################################################################################   
    while t < T:  # Time loop until final time T
        if t + dt > T:  # Adjust the time step to not exceed T
            dtc.assign(T - t)
        # ---First stage---
        solv1.solve()
        u1.assign(u + du)  # Update the first stage solution
        # ---Second stage---
        solv2.solve()
        u2.assign(0.75*u + 0.25*(u1+du))  # Update the second stage solution
        # ---Third stage---
        solv3.solve()
        u.assign((1.0/3.0)*u + (2.0/3.0)*(u2+du))  # Update the final solution

        t+= float(dtc)  # Increment time
    ######################################################################################
    # Plot the final solution
    ######################################################################################
    u_exact = Function(V).interpolate(sin(x[0]-T))  # Exact solution at final time
    print("Final time:", t)

    # Create fine grid for plotting
    num_plot_points = 100
    x_vals = np.linspace(left, right, num_plot_points)

    u_vals = np.array([u.at(xi,) for xi in x_vals])
    u_exact_vals = np.array([u_exact.at(xi,) for xi in x_vals])

    # Plot
    plt.figure(figsize=(8, 5))
    plt.plot(x_vals, u_vals, label="DG solution", linewidth=2)
    plt.plot(x_vals, u_exact_vals, "--", label="Exact solution", linewidth=2)
    plt.xlabel("x")
    plt.ylabel("u(x)")
    plt.title(f"DG vs Exact at t = {T:.2f}")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()

[stephankramer]

This is incorrect:

L1 = dtc*(u* phi.dx(0)*dx+(u('-')*jump(phi)*dS))

Remember that Firedrake is a general, unstructured FEM framework in which the orientation of facets (the boundaries between cells) are arbitrary. There is therefore no expectation that '-' refers to either the left or right-hand side of the cell boundary in 1D. If you follow the example more closely, and use:

n = FacetNormal(mesh)
un = 0.5*(n[0]+abs(n[0]))
L1 = dtc*(u* phi.dx(0)*dx-(un('+')*u('+')-un('-')*u('-'))*(phi('+')-phi('-'))*dS)  # Right-hand side linear form

your example seems to work correctly. Note that I have simplified the expression for un - you have no "advective velocity" but if you simply use vel=[1,] i.e. the 1D vector field that has just 1 component with value one, and is always pointing to the right, you get dot(vel, n)=n[0] the first component of the normal. The way this works in 1D, is that on the left of the facet (which can be either '+' or '-'), it uses the normal going out of the cell on the left, which is therefore pointing to the right with n[0]=1, and vice verse, and therefore un is 1 on the "upwind" (left) side and 0 on the right. Then the other thing that needs changing is, you used jump(u) but that's not the difference between right and left, it's actually a vector which has an orientation - if you write out the two cases where '+' and '-' are either left/right or right/left, you'll see that the above gives the right expression in both cases.

[afgalindo]

Thank you so much ! This was very helpful!