Unexpected results when solvind advection-diffusion equation

General
1

Hi all,

I am trying to solve a simple advection-diffusion equation using dolfinx, but my results are not converging or it is not solving correctly…
My domain is a simple square with an inner square inside that acts as a source, which I define as explained in this tutorial. The initial condition is just 0 everywhere.

f = fem.function(Q)
f.x.array[:] = 0
cells = ct.find(source_tag)
J.x.array[cells] = 1.0

The functions and solver are set as in the tutorial of the heat equation, with the addition of a vector function space for the velocity of the fluid:

W = fem.VectorFunctionSpace(domain, ("CG", 2))

class FluidVelocity():
    def __init__(self, t, max_speed=1.0):
        self.t = t
        self.max_speed = max_speed

    def __call__(self, x):
        values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType)
        values[0] = self.max_speed * ((x[1] + 2) / 4)**2
        return values

v_fluid = fem.Function(W)
v_fluid.name = "v_fluid"
fluid_velocity = FluidVelocity(t)
v_fluid.interpolate(fluid_velocity)

Lastly, the variational formulation looks like

a = u * v * ufl.dx + dt * D * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx \
+ dt * v * ufl.dot(v_fluid, ufl.grad(u)) * ufl.dx
L = (u_n + dt * f) * v * ufl.dx
bilinear_form = fem.form(a)
linear_form = fem.form(L)

A = fem.petsc.assemble_matrix(bilinear_form, bcs=[bc])
A.assemble()
b = fem.petsc.create_vector(linear_form)

The simulation runs fine but then the result is just numerical zero everywhere in the mesh, which suggest that the source is not being integrated correctly or else I have no idea.

Any help would be appreciated :slight_smile: thanks.

2

Please add a fully reproducible code, as it would be quite a lot of guesswork involved to reproduce your results with the current snippets. Also try to pin down your problem to a single time step, no need to have multiple time steps.

3

I hope it is more clear now, thank you.

import gmsh
import numpy as np
import tqdm

from mpi4py import MPI
from petsc4py import PETSc

from dolfinx import fem, mesh, io, plot
import ufl

# Define temporal parameters
t = 0 # Start time
T = 0.01 # Final time
num_steps = 1
dt = T / num_steps # time step size

# Define mesh
x0, y0 = -2, -2
L = 4
l = 0.2

gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
gmsh.initialize()

fluid_tag, source_tag = 11, 22
source, fluid = [], []
res_min = 0.05

if mesh_comm.rank == model_rank:
    box = gmsh.model.occ.addRectangle(x0, y0, 0, L, L)
    obstacle = gmsh.model.occ.addRectangle(0.1, 0.1, 0, l, l)
    
    whole = gmsh.model.occ.fragment([(gdim, box)], [(gdim, obstacle)])
    gmsh.model.occ.synchronize()
    
    volumes = gmsh.model.occ.getEntities(2)
    
    for dim, tag in volumes:
        mass = gmsh.model.occ.getMass(dim, tag)
        
        if np.isclose(mass, l * l):
            source.append(tag)
        
        else:
            fluid.append(tag)
    
    gmsh.model.addPhysicalGroup(2, fluid, fluid_tag)
    gmsh.model.setPhysicalName(2, fluid_tag, "fluid")
    
    gmsh.model.addPhysicalGroup(2, source, source_tag)
    gmsh.model.setPhysicalName(2, source_tag, "source")
    
    gmsh.model.occ.synchronize()
    edges = gmsh.model.getBoundary([(2, source[0])], oriented=False)
    gmsh.model.mesh.field.add("Distance", 1)
    gmsh.model.mesh.field.setNumbers(1, "EdgesList", [e[1] for e in edges])
    gmsh.model.mesh.field.add("Threshold", 2)
    gmsh.model.mesh.field.setNumber(2, "IField", 1)
    gmsh.model.mesh.field.setNumber(2, "LcMin", l / 3)
    gmsh.model.mesh.field.setNumber(2, "LcMax", 2 * l)
    gmsh.model.mesh.field.setNumber(2, "DistMin", 2 * l)
    gmsh.model.mesh.field.setNumber(2, "DistMax", 6 * l)
    gmsh.model.mesh.field.add("Min", 5)
    gmsh.model.mesh.field.setNumbers(5, "FieldsList", [2])
    gmsh.model.mesh.field.setAsBackgroundMesh(5)

gmsh.option.setNumber("Mesh.Algorithm", 7)
gmsh.model.mesh.generate(gdim) 
gmsh.model.mesh.optimize("Netgen")

domain, ct, ft = io.gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
ct.name = 'Cell markers'
ft.name = 'Facet markers'
gmsh.finalize()

# # Function spaces
V = fem.FunctionSpace(domain, ("CG", 1)) 
W = fem.VectorFunctionSpace(domain, ("CG", 2)) 
Q = fem.FunctionSpace(domain, ("DG", 0))


# initial condition
def initial_condition(x, a=5):
    return np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)


# Fluid velocity
class FluidVelocity():
    def __init__(self, t, max_speed=1.0):
        self.t = t
        self.max_speed = max_speed

    def __call__(self, x):
        values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType)
        values[0] = self.max_speed * ((x[1] + 2) / 4)**2
        return values


# Function for the "previous" timestep
u_n = fem.Function(V)
u_n.name = "u_n"
u_n.interpolate(initial_condition)

# Function for the present timestep
uh = fem.Function(V)
uh.name = "uh"
uh.interpolate(initial_condition)

# Function for the fluid velocity
v_fluid = fem.Function(W)
v_fluid.name = "v_fluid"
fluid_velocity = FluidVelocity(t)
v_fluid.interpolate(fluid_velocity)

# Trial and Test for the variational problem
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)

# BCs
fdim = domain.topology.dim - 1
facets = mesh.locate_entities_boundary(domain, fdim, lambda x: np.full(x.shape[1], True))
dofs = fem.locate_dofs_topological(V, fdim, facets)
bc = fem.dirichletbc(PETSc.ScalarType(0), dofs, V)

## Source and diffusion coefficient
f = fem.Function(Q)
f.x.array[:] = 0
cells = ct.find(source_tag)
f.x.array[cells] = 0.5
D = fem.Constant(domain, PETSc.ScalarType(0.5))

# Variational problem
a = u * v * ufl.dx + dt *D* ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx + dt * v * ufl.dot(v_fluid, ufl.grad(u)) * ufl.dx
L = (u_n + dt * f) * v * ufl.dx
bilinear_form = fem.form(a)
linear_form = fem.form(L)

A = fem.petsc.assemble_matrix(bilinear_form, bcs=[bc])
A.assemble()
b = fem.petsc.create_vector(linear_form)

# Solver
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

xdmf = io.XDMFFile(domain.comm, "result.xdmf", "w")
xdmf.write_mesh(domain)
xdmf.write_function(uh, t)

progress = tqdm.tqdm(total=num_steps)
for i in range(num_steps):
    progress.update(1)
    t += dt

    # Update the right hand side reusing the initial vector
    with b.localForm() as loc_b:
        loc_b.set(0)
    fem.petsc.assemble_vector(b, linear_form)
    
    # Apply Dirichlet boundary condition to the vector
    fem.petsc.apply_lifting(b, [bilinear_form], [[bc]])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    fem.petsc.set_bc(b, [bc])

    # Solve linear problem
    solver.solve(b, uh.vector)
    uh.x.scatter_forward()

    # Update solution at previous time step (u_n)
    u_n.x.array[:] = uh.x.array

    # Write solution to file
    xdmf.write_function(uh, t)

xdmf.close()