The mesh should correspond to whatever geometry you are simulating, with the disk B cut out. To generate the square mesh you mentioned, you could use the following in gmsh:
SetFactory("OpenCASCADE");
L = 20; // size of the box
R = 2; // radius of the disk
Rectangle(1) = {-L/2, -L/2, 0, L, L, 0};
Disk(2) = {0, 0, 0, R, R};
BooleanDifference{ Surface{1}; Delete; }{ Surface{2}; Delete; }
MeshSize {6,7,8,9} = L/10;
MeshSize {5} = R/10;
Physical Surface(1) = {1};
Physical Curve(2) = {5};
Physical Curve(3) = {6,7,8,9};
with the command
gmsh mesh.geo -2 -format msh4 -o mesh.msh
to generate the mesh, and the following Python script to convert it to XDMF format:
import meshio
msh = meshio.read("mesh.msh")
msh.prune_z_0()
triangle_data = msh.cell_data_dict["gmsh:physical"]["triangle"]
line_data = msh.cell_data_dict["gmsh:physical"]["line"]
triangle_cells = msh.cells_dict["triangle"]
line_cells = msh.cells_dict["line"]
triangle_mesh =meshio.Mesh(points=msh.points,
cells=[("triangle", triangle_cells)],
cell_data={"cell_marker":[triangle_data]})
line_mesh = meshio.Mesh(points=msh.points, cells={"line": line_cells},
cell_data={"line_marker":[line_data]})
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("bnd_markers.xdmf", line_mesh)
Your PDE can be solved, starting from some initial condition P(x,0), using the following code. In the code below, you will see that I have considered a homogeneous condition on the external boundary of G, as well as three different inhomogeneous initial conditions. The solution for all three of these initial conditions eventually settles to the same steady-state, namely the solution to
(u \cdot \nabla)P - 2 \Delta P = 0
Based on this solution, it does not appear that [-10, 10] \times [-10, 10] is a large enough domain to accurately represent the entire x_1 x_2-plane, so you may wish to consider experimenting with the mesh. In addition, this solution does not satisfy the final condition P(x,T)=0 that you have specified, so you should double-check the formulation of the problem.
from fenics import *
import numpy as np
T = 1000.0 # final time
N = 100
dt = T/N
Dt = Constant(dt) # timestep
# Create mesh and define function space
mesh = Mesh()
with XDMFFile("mesh.xdmf") as xdmf:
xdmf.read(mesh)
V = FunctionSpace(mesh, 'P', 2)
# Load boundary markers
mvc = MeshValueCollection('size_t', mesh, mesh.topology().dim() - 1)
with XDMFFile("bnd_markers.xdmf") as xdmf:
xdmf.read(mvc)
mf = MeshFunction('size_t', mesh, mvc)
# Define boundary condition
P_disc = Constant(1.0)
P_outer = Constant(0.0)
bc_disc = DirichletBC(V, P_disc, mf, 2)
bc_outer = DirichletBC(V, P_outer, mf, 3)
bcs = [bc_disc, bc_outer]
# Define initial value
# P_old = interpolate(Expression('0.', degree=2),V)
# P_old = interpolate(Expression('1.', degree=2),V)
P_old = interpolate(Expression('x[0]*x[0] + x[1]*x[1] - 4', degree=2),V)
#Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
gamma = Constant(1.5)
f1 = Expression(('1.','1.'), degree=0)
f2 = Expression(('4*pow(x[0],3)/(pow(x[0],4)+pow(x[1],4)+DOLFIN_EPS)',
'4*pow(x[1],3)/(pow(x[0],4)+pow(x[1],4)+DOLFIN_EPS)'),
degree=4)
f = f1 + gamma*f2
F = (u*v + Dt*inner(f, grad(u))*v + 2*Dt*inner(grad(u), grad(v)) - P_old*v)*dx
a, L = lhs(F), rhs(F)
#Time-stepping
P = Function(V)
P.rename("P", "P")
P.assign(P_old)
t = 0.0
xdmf = XDMFFile("result_Pparab_initial.xdmf")
# xdmf.write(P, t)
while t < T:
#Update current time
t += dt
#Compute solution
solve(a == L, P, bcs)
#Plot solution
xdmf.write(P, t)
#Update previous solution
P_old.assign(P)
xdmf.close()
