# Unexpected results when solvind advection-diffusion equation

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 thanks.

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.

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.setPhysicalName(2, fluid_tag, "fluid")

gmsh.model.setPhysicalName(2, source_tag, "source")

gmsh.model.occ.synchronize()
edges = gmsh.model.getBoundary([(2, source[0])], oriented=False)
gmsh.model.mesh.field.setNumbers(1, "EdgesList", [e[1] for e in edges])
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.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]])
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()
``````