Hi,
I want to solve the Poisson equation on a domain like this:
Where I want the normal component of the gradient of the solution to bei zero on red and apply different constant Dirichlet BCs on blue and white.
When I project the negative gradient of the solution I obtain following field:
What I noticed is that the velocity spike at the node between the inflow and free slip boundaries increases with increasing mesh density. My question is whether this is physical, or something is wrong somewhere.
an mwe:
from dolfin import *
from matplotlib import pyplot as plt
class Boundary(SubDomain):
def inside(self,x, on_boundary):
return on_boundary
class Inflow(SubDomain):
def inside(self,x, on_boundary):
return (.3 <= x[0] <= .7) and x[1] < .001 and on_boundary
class Top(SubDomain):
def inside(self,x, on_boundary):
return near(x[1],1) and on_boundary
mesh = UnitSquareMesh(25,25,diagonal='left')
CG1 = FunctionSpace(mesh,'CG',1)
DG0 = VectorFunctionSpace(mesh,'DG',0)
p = TrialFunction(CG1)
q = TestFunction(CG1)
n = FacetNormal(mesh)
sub_domains = MeshFunction('size_t', mesh, 1)
sub_domains.set_all(0)
Boundary().mark(sub_domains, 3)
Top().mark(sub_domains, 2)
Inflow().mark(sub_domains, 1)
ds = Measure('ds', domain = mesh, subdomain_data = sub_domains)
a = dot(grad(p), grad(q))*dx
L = Constant(0)*q*dx
bc1 = DirichletBC(CG1, Constant(1),sub_domains,1)
bc2 = DirichletBC(CG1, Constant(0),sub_domains,2)
bc = [bc1, bc2]
p = Function(CG1)
solve(a == L, p,bc,solver_parameters={'linear_solver': 'mumps'})
u = project(grad(-p), DG0)
plt.figure(dpi=400)
c = plot(p)
plot(mesh, linewidth=.2)
plt.colorbar(c)
plt.savefig("p.png")
plt.figure(dpi=400)
c = plot(u)
plot(mesh, linewidth=.2)
plt.colorbar(c)
plt.savefig("u.png")