Advice on solving advection equation

variational formulation
6

Thanks very much. I think I’m still very far away from even understanding what is going wrong… @Smith_Jack the step-9 link is exactly what I tried to do with the penalty term I mentioned earlier. That does not work (or rather does not work as I’d expect, and the penalty weight affects the solution).

I tried to adapt the upwind flux formulation from Lax-Friedrichs Flux for advection equation :

from dolfinx.mesh import create_unit_square
from mpi4py import MPI
import numpy as np
from dolfinx import mesh as dolf_mesh
from dolfinx import fem, io
import ufl
from ufl import dx, grad, inner, jump
from dolfinx.fem.petsc import LinearProblem


mesh = create_unit_square(MPI.COMM_WORLD, 100, 100)
V_e = fem.functionspace(mesh, ("DG", 2))

# An example field `T`
T = fem.Function(V_e)
T.interpolate(lambda x: 0.5 - x[1])


def bot_boundary(x):
    return np.isclose(x[1], 0.0)


def rest_boundary(x):
    return np.logical_not(bot_boundary(x))


bottom_facets = dolf_mesh.locate_entities_boundary(
    mesh, mesh.topology.dim - 1, bot_boundary
)
bottom_vals = np.full(bottom_facets.shape, 1, np.intc)
rest_facets = dolf_mesh.locate_entities_boundary(
    mesh, mesh.topology.dim - 1, rest_boundary
)
rest_vals = np.full(rest_facets.shape, 2, np.intc)

indices = np.hstack((bottom_facets, rest_facets))
values = np.hstack((bottom_vals, rest_vals))
indices, pos = np.unique(indices, return_index=True)
marker = dolf_mesh.meshtags(mesh, mesh.topology.dim - 1, indices, values[pos])
ds = ufl.Measure("ds", subdomain_data=marker, domain=mesh)
dS = ufl.Measure("dS", domain=mesh)
n = ufl.FacetNormal(mesh)

v = ufl.TestFunction(V_e)
u = ufl.TrialFunction(V_e)

beta = grad(T)

a = -inner(beta, grad(v)) * u * dx
a += v * u * inner(beta, n) * ds(2)  # Do we need this?

un = (inner(beta, n) + abs(inner(beta, n)))/2.0
flux = jump(un*u)

a += inner(jump(v), flux)*dS
L = v * T * ds(1)
problem = LinearProblem(a, L, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

uh = problem.solve()
print(uh.x.array)
with io.VTXWriter(mesh.comm, "mvp.bp", [uh]) as file:
    file.write(0.0)

but the resulting vector uh is just [inf inf inf ... inf inf inf]… I’d appreciate any hint of what might be wrong here, or some reference implementation as this seems to be such a standard problem. Thanks very much.