Rotation of a Gaussian hill

This example is taken from Section 13.2 in Kuzmin and Turek (2004). The (non-dimensional) equation we are solving is a linear convection-diffusion equation in 2D, given as

\[\partial_t u + \boldsymbol{b} \cdot \nabla u - \kappa \Delta u = 0, \qquad \forall\, \boldsymbol{x} \in \Omega = [-1,+1]^2, \quad t \in [t_0, t_f],\]

where \(t_0 = \pi/2\), \(t_f = 5\pi/2\), with a diffusion constant of \(\kappa = 0.001\) and a spatially-dependent wind velocity profile: \(\boldsymbol{b} = (-y, x)^T\).

The exact solution to this benchmark is a time-dependent Gaussian of the form

\[u(x,y,t) = \frac{1}{4 \pi \kappa t} e^{-r^2/(4 \kappa t)},\]

where \(r^2 = (x - \widehat{x})^2 + (y - \widehat{y})^2\) and the center of the pulse is time-dependent, as such

\[x(t) = x_0 \cos(t) - y_0 \sin(t), \qquad y(t) = -x_0 \sin(t) + y_0 \cos(t),\]

with \((x_0, y_0) = (0, 0.5)\) denoting the initial position of the pulse at \(t = t_0\).

[1] Kuzmin, Dmitri, and Stefan Turek.
“High-resolution FEM-TVD schemes based on a fully multidimensional flux limiter.”
Journal of Computational Physics, 198.1 (2004): 131-158.
[1]:
# Import preliminary modules.
import numpy as np
import ngsolve as ngs
from ngsolve.webgui import Draw
from netgen.occ import OCCGeometry, WorkPlane

from dream.scalar_transport import Initial, transportfields, ScalarTransportSolver, FarField
[2]:
# Define the grid.

def CreateSimpleGrid(ne, lx, ly):

    # Select a common element size.
    h0 = min( lx, ly )/float(ne)

    # Generate a simple rectangular geometry.
    domain = WorkPlane().RectangleC(lx, ly).Face()

    # Assign the name of the internal solution in the domain.
    domain.name = 'internal'

    # For convenience, extract and name each of the edges consistently.
    bottom = domain.edges[0]; bottom.name = 'bottom'
    right  = domain.edges[1]; right.name  = 'right'
    top    = domain.edges[2]; top.name    = 'top'
    left   = domain.edges[3]; left.name   = 'left'

    # Initialize a rectangular 2D geometry.
    geo = OCCGeometry(domain, dim=2)

    # Discretize the domain.
    mesh = ngs.Mesh(geo.GenerateMesh(maxh=h0, quad_dominated=True))

    # Return our fancy grid.
    return mesh

[3]:
# Generate the grid.

# Number of elements per dimension.
nElem1D = 20

# Dimension of the rectangular domain.
xLength = 2.0
yLength = 2.0

# Generate a simple grid.
mesh = CreateSimpleGrid(nElem1D, xLength, yLength)

# Message output detail from netgen.
ngs.ngsglobals.msg_level = 0
ngs.SetNumThreads(4)
Draw(mesh)
[3]:
BaseWebGuiScene
[4]:
# Define analytic solution.
def get_analytic_solution(t, k):

    # Initial pulse location.
    x0 = 0.0
    y0 = 0.5

    # Pulse center trajectory.
    xc =  x0*ngs.cos(t) - y0*ngs.sin(t)
    yc = -x0*ngs.sin(t) + y0*ngs.cos(t)

    # Radial distance (time-dependent).
    r2 = (ngs.x-xc)**2 + (ngs.y-yc)**2
    # Variance of this pulse.
    s2 = get_variance_pulse(t, k)

    # Return the analytic solution.
    return ( 1.0/(s2*ngs.pi) ) * ngs.exp( -r2/(4.0*k*t) )

def get_variance_pulse(t, k):
    return 4.0*k*t

[5]:
# Solver configuration: Scalar transport equation.
cfg = ScalarTransportSolver(mesh)
cfg.time = "transient"
cfg.time.timer.interval = (ngs.pi/2, 5*ngs.pi/2)
cfg.time.timer.step = 0.01

cfg.convection_velocity = (-ngs.y, ngs.x)
cfg.diffusion_coefficient = 1.0e-03
cfg.is_inviscid = False

cfg.riemann_solver = "lax_friedrich"
cfg.fem = "dg" # NOTE, by default, DG is used.
cfg.fem.interior_penalty_coefficient = 10.0
cfg.fem.order = 4
cfg.fem.scheme = "implicit_euler"



U0 = transportfields()
U0.phi = get_analytic_solution(cfg.time.timer.interval[0], cfg.diffusion_coefficient)

Ubc = transportfields()
Ubc.phi = 0.0

cfg.bcs['left|top|bottom|right'] = FarField(fields=Ubc)
cfg.dcs['internal'] = Initial(fields=U0)
[6]:
cfg.initialize()

with ngs.TaskManager():
    # Solve the transport equation.
    cfg.solve()

cfg.io.draw(cfg.fem.get_fields())

dream.time      (INFO) | dg implicit_euler | t: 7.85
[7]:
# Decorator for the actual simulation.

# Extract time configuration.
t0, tf = cfg.time.timer.interval
dt = cfg.time.timer.step.Get()
nt = int(round((tf - t0) / dt))

# Define a decorator for the Gaussian hill routine.
def gaussian_hill_routine(label):
    def decorator(func):
        def wrapper(*args, draw_solution=False, **kwargs):

            # Insert options here.
            func(*args, **kwargs)
            cfg.fem.order = 4
            cfg.fem.interior_penalty_coefficient = 2.0

            # Allocate the necessary data.
            cfg.initialize()

            # Get a reference to the numerical solution and the analytic solution.
            uh = cfg.fem.get_fields("phi").phi
            ue = get_analytic_solution(cfg.time.timer.t, cfg.diffusion_coefficient)

            if draw_solution:
                # cfg.io.draw({"phi": Uh.phi})
                # cfg.io.draw({"Exact[phi]": Uex.phi},  min=0.0, max=10.0)
                cfg.io.draw({"Diff[phi]": (ue - uh)}, min=-0.01, max=0.01)

            # Integration order (for post-processing).
            qorder = 10

            # Data for book-keeping information.
            data = np.zeros((nt, 3), dtype=float)

            with ngs.TaskManager():
                # Time integration loop  (actual simulation).
                for i, t in enumerate(cfg.time.start_solution_routine(True)):

                    # Get the variance, based on the analytic solution.
                    s2e = get_variance_pulse(cfg.time.timer.t.Get(), cfg.diffusion_coefficient.Get())

                    # Compute centroid
                    xh = ngs.Integrate(ngs.x * uh, mesh, order=qorder)
                    yh = ngs.Integrate(ngs.y * uh, mesh, order=qorder)

                    # Compute variance.
                    r2h = (ngs.x - xh)**2 + (ngs.y - yh)**2
                    s2h = ngs.Integrate(r2h * uh, mesh, order=(qorder+2) )

                    # Compute the normalized variance error.
                    var_dif = s2h/s2e - 1.0

                    # Compute the L2-norm of the error.
                    err = np.sqrt( ngs.Integrate( (ue-uh)**2, mesh, order=qorder) )

                    # Store data: time and error metrics.
                    data[i] = [cfg.time.timer.t.Get(), var_dif, err]
            return data
        wrapper.label = label
        return wrapper
    return decorator
[8]:
# Specialized routines.

@gaussian_hill_routine("implicit_euler(hdg)")
def implicit_euler_hdg():
    cfg.fem = "hdg"
    cfg.fem.scheme = "implicit_euler"

@gaussian_hill_routine("implicit_euler(dg)")
def implicit_euler_dg():
    cfg.fem = "dg"
    cfg.fem.scheme = "implicit_euler"

@gaussian_hill_routine("bdf2(hdg)")
def bdf2_hdg():
    cfg.fem = "hdg"
    cfg.fem.scheme = "bdf2"
@gaussian_hill_routine("bdf2(dg)")
def bdf2_dg():
    cfg.fem = "dg"
    cfg.fem.scheme = "bdf2"

@gaussian_hill_routine("sdirk22(hdg)")
def sdirk22_hdg():
    cfg.fem = "hdg"
    cfg.fem.scheme = "sdirk22"
@gaussian_hill_routine("sdirk22(dg)")
def sdirk22_dg():
    cfg.fem = "dg"
    cfg.fem.scheme = "sdirk22"

@gaussian_hill_routine("sdirk33(hdg)")
def sdirk33_hdg():
    cfg.fem = "hdg"
    cfg.fem.scheme = "sdirk33"
@gaussian_hill_routine("sdirk33(dg)")
def sdirk33_dg():
    cfg.fem = "dg"
    cfg.fem.scheme = "sdirk33"

@gaussian_hill_routine("imex_rk_ars443(dg)")
def imex_rk_ars443_dg():
    cfg.fem = "dg"
    cfg.fem.scheme = "imex_rk_ars443"

@gaussian_hill_routine("ssprk3(dg)")
def ssprk3_dg():
    cfg.fem = "dg"
    cfg.fem.scheme = "ssprk3"

@gaussian_hill_routine("crk4(dg)")
def crk4_dg():
    cfg.fem = "dg"
    cfg.fem.scheme = "crk4"

@gaussian_hill_routine("explicit_euler(dg)")
def explicit_euler_dg():
    cfg.fem = "dg"
    cfg.fem.scheme = "explicit_euler"
[9]:
# Figure style and properties.
from matplotlib import pyplot as plt
from matplotlib.ticker import MultipleLocator, FuncFormatter
from fractions import Fraction


# Assign colors to base method names
method_colors = {
    "implicit_euler": "tab:red",
    "bdf2": "tab:blue",
    "sdirk22": "tab:green",
    "sdirk33": "tab:purple"
}

# Helper function that sets up a generic figure.
def setup_fig(xlabel, ylabel, title):
    fig, ax = plt.subplots(figsize=(15, 6))
    ax.set_xlim(t0, tf)
    ax.tick_params(axis='both', labelsize=16)
    ax.set_xlabel(xlabel, fontsize=14)
    ax.set_ylabel(ylabel, fontsize=16)
    ax.set_title(title, fontsize=16)
    ax.grid(True)
    ax.xaxis.set_major_locator(MultipleLocator(base=np.pi/2))
    ax.xaxis.set_major_formatter(FuncFormatter(pi_formatter))
    return fig,ax

# Helper for formatting ticks as multiples of pi.
def pi_formatter(x, pos):
    frac = Fraction(x / np.pi).limit_denominator(8)
    if frac.numerator == 0:
        return "$0$"
    elif frac.denominator == 1:
        return rf"${frac.numerator}\pi$"
    else:
        return rf"${frac.numerator}\pi/{frac.denominator}$"

# Helper function that specifies the plot style properties.
def get_plot_style(label, method_colors):
    if label.endswith("(hdg)"):
        base_name = label.replace("(hdg)", "").strip()
        linestyle = "-"
        marker = 'o'
    elif label.endswith("(dg)"):
        base_name = label.replace("(dg)", "").strip()
        linestyle = "--"
        marker = '^'
    else:
        base_name = label.strip()
        linestyle = ":"
        marker = '.'

    color = method_colors.get(base_name, "black")

    return {
        "label": label,
        "linestyle": linestyle,
        "marker": marker,
        "markersize": 8,
        "markevery": 50,
        "color": color
    }

[10]:
# Run simulation(s).
routines = [implicit_euler_hdg,
                      bdf2_hdg,
                   sdirk22_hdg,
                   sdirk33_hdg]
# Figure 1: Linear plot of the relative error in the variance.
fig1, ax1 = setup_fig("Time", r"$s^2_h / s^2_e - 1$", "Time Evolution of Relative Variance Error")

# Figure 2: logaraithmic plot of the relative error in the L2-norm.
fig2, ax2 = setup_fig("Time", r"$\| u_e - u_h \|$", r"Time Evolution of Relative L2-norm $\| u_e - u_h \|$")

# Run each simulation.
for routine in routines:
    print(f"Running {routine.label}...")
    data = routine()
    style = get_plot_style(routine.label, method_colors)

    ax1.plot(data[:, 0], data[:, 1], **style)
    ax2.semilogy(data[:, 0], data[:, 2], **style)

# Add IMEX-RK for DG only.
data = imex_rk_ars443_dg()
ax1.plot(data[:, 0], data[:, 1], label=imex_rk_ars443_dg.label, linestyle='--', color='k')
ax2.semilogy(data[:, 0], data[:, 2], label=imex_rk_ars443_dg.label, linestyle='--', color='k')

# Final touches.
ax1.legend(frameon=False)
ax2.legend(loc="upper right", frameon=False)
fig1.tight_layout()
fig2.tight_layout()
plt.show()
Running implicit_euler(hdg)...
dream.time      (INFO) | hdg implicit_euler | t: 7.85
Running bdf2(hdg)...
dream.time      (INFO) | hdg sdirk22 | t: 1.65
Running sdirk22(hdg)...
dream.time      (INFO) | hdg sdirk33 | t: 1.63
Running sdirk33(hdg)...
dream.time      (INFO) | dg imex_rk_ars443 | t: 7.85
../_images/examples_gaussian_hill_farfield_10_8.png
../_images/examples_gaussian_hill_farfield_10_9.png
[ ]: