Hello
I am trying to solve Darcy equation for a domain with two-layers using DG formulation.
My governing equations reads as:
\begin{align}
\alpha \mathbf{u} + \mathrm{grad}[p] = \gamma \mathbf{b} \quad \mathrm{in} \; \Omega \
, \mathrm{div}[\mathbf{u}] = 0 \quad \mathrm{on} \; \Omega \
, \mathbf{u} \cdot \widehat{\mathbf{n}} = u_n \quad \mathrm{on} \; \Gamma^v \
, p = p_0 \quad \mathrm{on} \; \Gamma^{p}
\end{align}
and Darcy DG formulation based on B.~Santiago, 2010 CMAME could be written as:
\begin{align}
\int_{\Omega} \alpha \mathbf{u} \cdot \mathbf{w} \; d \Omega
%
+\int_{\Omega} \mathrm{grad} p \cdot \mathbf{w} \; d \Omega
%
- \int_{\Gamma^{int}} [[ p ]] \cdot < \mathbf{w} > \; d \Gamma
%
-\int_{\Omega} \mathbf{u} \cdot \mathrm{grad}[q] \; d \Omega
\end{align}
\begin{align}
+\int_{\Gamma^{int}} [[ q ]] \cdot < \mathbf{u} > \; d \Gamma
%
= \int_{\Omega} \gamma \mathbf{b} \cdot \mathbf{w} \; d \Omega
- \int_{\Gamma^{p}} p_0 \widehat{\mathbf{n}} \cdot \mathbf{w} \; d \Gamma
\end{align}
or equivalently, in the divergence form:
\begin{align}
\int_{\Omega} \alpha \mathbf{u} \cdot \mathbf{w} \; d \Omega
%
-\int_{\Omega} p \mathrm{div}[\mathbf{w}] \; d \Omega
%
+ \int_{\Gamma^{int}} [[ \mathbf{w} ]] < p > \; d \Gamma
%
+\int_{\Omega} \mathrm{div}[\mathbf{u}] q \; d \Omega
\end{align}
\begin{align}
-\int_{\Gamma^{int}} [[ \mathbf{u} ]] < q > \; d \Gamma
%
= \int_{\Omega} \gamma \mathbf{b} \cdot \mathbf{w} \; d \Omega
- \int_{\Gamma^{p}} p_0 \widehat{\mathbf{n}} \cdot \mathbf{w} \; d \Gamma
\end{align}
Herein, $\Gamma^{int}$ denotes the integration on internal edges. $[[ \cdot ]]$ denotes jump operator and $< \cdot >$ denotes average operator. The pressure BC is prescribed strongly in the weak form and velocity BC will be imposed strongly.
Domain:[C.f. Hughes, Masud, Wan, 2006 CAME] We have a square domain divided into upper and lower layers with different $\alpha_t$ and $\alpha_b$, respectively. On top left and top right BCs, $u_n = v_{\mathrm{top}} = \frac{1}{\alpha_t}$ are prescribed. On bottom
left and right BCs, $u_n =v_{\mathrm{bottom}} = \frac{1}{\alpha_b}$ are prescribed. While top and bottom have zero normal velocities, for uniqueness, a reference $p = p_0 = 1.0$ is prescribed on bottom left corner. Here is the code:
from dolfin import *
from mshr import*
#=====================================;
# Create mesh and identify boundary ;
#=====================================;
mesh = RectangleMesh(Point(0.0, 0.0), Point(5.0, 4.0), 50, 50)
print("Plotting a RectangleMesh")
#plot(mesh, title="Rectangle (right/left)")
boundaries = MeshFunction("size_t",mesh, mesh.topology().dim()-1)
domains = CellFunction("size_t", mesh)
facets = FacetFunction("uint",mesh)
#====================================================;
# Define function spaces and mixed (product) space ;
#====================================================;
vSpace = VectorFunctionSpace(mesh,"DG",1)
pSpace = FunctionSpace(mesh,"DG",1)
wSpace = vSpace * pSpace
#===================================;
# Define trial and test functions ;
#===================================;
(v,p) = TrialFunctions(wSpace)
(w,q) = TestFunctions(wSpace)
#======================;
# Define body forces ;
#======================;
rhob= Constant((0.0,0.0))
#============================;
# Define medium properties ;
#============================;
alpha_b = Constant(1.)
alpha_t = Constant(0.1)
p_corner=Constant(1.)
v_b = Constant(1.0)
v_t = Constant(10.0)
v_bottom = Constant(0.0)
v_top = Constant(0.0)
#===========================;
# Identify the boundaries ;
#===========================;
class DownRegion(SubDomain):
def inside(self,x,on_boundary):
return (between(x[0], (0.0, 5.0)) and between(x[1], (0.0, 2.0)))
class UpRegion(SubDomain):
def inside(self,x,on_boundary):
return (between(x[0], (0.0, 5.0)) and between(x[1], (2.0, 4.0)))
class LD(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and (near(x[0],0.0) and between(x[1], (0.0, 2.0)))
class LU(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and (near(x[0],0.0) and between(x[1], (2.0, 4.0)))
class RD(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and (near(x[0],5.0) and between(x[1], (0.0, 2.0)))
class RU(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and (near(x[0],5.0) and between(x[1], (2.0, 4.0)))
class Bottom(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and (near(x[1],0.0) and between(x[0], (0.0, 5.0)))
class Top(SubDomain):
def inside(self,x,on_boundary):
return on_boundary and (near(x[1],4.0) and between(x[0], (0.0, 5.0)))
class Corner(SubDomain):
def inside(self,x,on_boundary):
return x[0] < DOLFIN_EPS and x[1] < DOLFIN_EPS
DownRegion = DownRegion()
UpRegion = UpRegion()
LD = LD()
LU = LU()
RD = RD()
RU = RU()
Bottom = Bottom()
Top = Top()
Corner = Corner()
domains.set_all(0)
DownRegion.mark(domains, 1)
UpRegion.mark(domains, 2)
boundaries.set_all(0)
LD.mark(boundaries,3)
LU.mark(boundaries,9)
RD.mark(boundaries,5)
RU.mark(boundaries,7)
Bottom.mark(boundaries,4)
Top.mark(boundaries,8)
Corner.mark(boundaries,10)
#=================================;
# Dirichlet boundary conditions ;
#=================================;
BC_LD = DirichletBC(wSpace.sub(0).sub(0),v_b,LD)
BC_RD = DirichletBC(wSpace.sub(0).sub(0),v_b,RD)
BC_LU = DirichletBC(wSpace.sub(0).sub(0),v_t,LU)
BC_RU = DirichletBC(wSpace.sub(0).sub(0),v_t,RU)
BC_Top = DirichletBC(wSpace.sub(0).sub(1),Constant(0.0),Top)
BC_Bottom = DirichletBC(wSpace.sub(0).sub(1),Constant(0.0),Bottom)
bcs = [BC_LD,BC_RD,BC_LU,BC_RU,BC_Top,BC_Bottom]
#=====================================;
# Define normal vector and mesh size ;
#=====================================;
dx = Measure("dx")[domains]
ds = Measure("ds")[boundaries]
dS = Measure("dS")[facets]
n = FacetNormal(mesh)
h = CellSize(mesh)
#===========================;
# Define variational form ;
#===========================
a = dot(w, alpha_b * v) * dx(1) +\
dot(w, alpha_t * v) * dx(2) +\
dot(grad(p) , w) * dx -\
dot(v , grad(q)) * dx -\
dot(jump(p,n) , avg(w))* dS +\
dot(jump(q,n) , avg(v)) * dS
L = dot(w,rhob)*dx -\
dot(p_corner * n , w) * ds(10)
#====================;
# Compute solution ;
#====================;
solution = Function(wSpace)
solve(a == L,solution,bcs)
(vsol,psol) = solution.split(True)
#=======================================;
# Dump solution to file in VTK format ;
#=======================================;
file = File('V_DG.pvd')
file << vsol
file = File('P_DG.pvd')
file << psol
The exact solution of the problem is zero vertical velocity everywhere, constant horizontal velocity equal to $\frac{1}{\alpha}$ in each layer, and the same horizontal linear pressure variation in all layers.
However, this code returns zero velocity and pressure fields everywhere.
[I experimented with prescribing velocity BCs weakly or pressure BC strongly, but it returns
unreasonable answers]
I would be grateful if someone can help me.
