Advice on solving advection equation

variational formulation
1

Hi, I am trying to solve the following linear equation in a simple 2D domain (like a unit square):

\nabla ( u \nabla T) = 0 \ in \ \Omega, \\ T + u (n, \nabla T) = 0 \ on \ \Gamma_1 \\ \mathrm{(no \ bcs \ on \ } \Gamma_2 == \partial \Omega - \Gamma_1 )

Here u is the unknown scalar field, and T is a known scalar field (not sure that is relevant, but I obtained T as a solution to a different problem).

So far I could only write something like this for the weak formulation:

\int_\Omega - u (\nabla v, \nabla T) dx + \int_{\Gamma_2} v u (n, \nabla T) ds = \int_{\Gamma_2} v T ds

which translates to the code:

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
import ufl
from ufl import dx, grad, inner
from dolfinx.fem.petsc import LinearProblem


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

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


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


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


# Got this from https://github.com/FEniCS/dolfinx/blob/1ccffb0aa545634ea6e160236a7b043dc2a7fdd2/python/test/unit/fem/test_assemble_domains.py#L121
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)

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

a = -u * inner(grad(v), grad(T)) * dx
a += v * u * inner(grad(T), ufl.FacetNormal(mesh)) * ds(2)
L = v * T * ds(1)
problem = LinearProblem(a, L, bcs=[])

uh = problem.solve()

However the results just look rubbish, and in fact the solution sometimes consists of infs depending on my mesh (with petsc_options={"ksp_type": "preonly", "pc_type": "lu"} in the problem definition).

I am really struggling to understand what I am doing wrong and would appreciate any tips of how to solve this. Thank you!

2

Please write this in latex format, ie with $ encapsulation.

It is Also hard to Give you guidance without a minimal example at hand.

3

Thanks @dokken really appreciate that. I updated the latex and added an mvp code (I’m running on dolfinx/dolfinx:stable, ab9e82b923553955b6bc59078c2bf0d0bb9c0cc69e811453aecbca7710921ae6)

4

Apparently adding a term like

a -= div(u * grad(T)) * div(v * grad(T)) * 100 * dx

helps a bit, but I’m not sure how robust this approach is…

Edit: helps == at least I could get it to work for case of a known analytical solution where T = 1 - y. But for non-trivial T, the penalty coefficient affects the solution.

5

Solving the advection equation usually needs some stabilization technique，such as SUPG method，for the details, see step-9 of the deal.II tutorial. Or you can read some books about the application of FEM in fluid dynamics, such as Finite elemenet methods for Flow problems written by Jean Donea and Antonio Huerta

1 Like
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.

7

I am still struggling to learn the fluid dynamics, so I am not sure where the problem is. :joy:
But the variational form I get is the following, there is no \int_{\Gamma_2} vu(n,\nabla T)ds term.
And from my understanding, Discontinuous Galerkin Finite element is used when the solution suddenly change in the domain (step-12), is this your case in your problem ?

image
image770×439 9.16 KB