If the Stokes-Darcy problem is dirichlet boundary condition, that is the speed of Stokes domain and the pressure of Darcy region are known. According to the code that previously dealt with Neumann boundary conditions we made the following modifications. However, that’s not what we want. Anybody know why?
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from fealpy.tools.show import showmultirate, show_error_table
from ngsolve import *
from netgen.geom2d import SplineGeometry
from numpy import pi
import netgen.gui
ngsglobals.msg_level=1
errorType = [‘||us - us_h||_{\Omega,0}’,
‘||ud - ud_h||_{\Omega, 0}’
]
maxit = 1
errorMatrix = np.zeros((2, maxit), dtype=np.float)
NDof = np.zeros(maxit, dtype=np.float)
def topBottom():
geometry = SplineGeometry()
point coordinates …
pnts = [ (0,0), (1,0), (1,0.5), (1,1), (0,1), (0, 0.5)]
pnums = [geometry.AppendPoint(*p) for p in pnts]
start-point, end-point, boundary-condition, domain on left side, domain on right side:
lines = [ (pnums[0],pnums[1],“gammaD”,2,0), (pnums[1],pnums[2],“gammaD”,2,0),
(pnums[2],pnums[3],“gammaS”,1,0), (pnums[3],pnums[4],“gammaS”,1,0),
(pnums[4],pnums[5],“gammaS”,1,0),
(pnums[5],pnums[0],“gammaD”,2,0), (pnums[5],pnums[2],“gammaSD”,1,2)]
for p1,p2,bc,left,right in lines:
geometry.Append( [“line”, p1, p2], bc=bc, leftdomain=left, rightdomain=right)
return geometry
def GetInterfaceIndicator(mesh):
dom_indicator = GridFunction(L2(mesh,order=0))
dom_indicator.Set(CoefficientFunction([1,2]))
ffes = FacetFESpace(mesh,order=0)
dom_ind_average = GridFunction(ffes)
mf = BilinearForm(ffes,symmetric=True)
mf += SymbolicBFI(ffes.TrialFunction()ffes.TestFunction(),element_boundary=True)
mf.Assemble()
ff = LinearForm(ffes)
ff += SymbolicLFI(dom_indicatorffes.TestFunction(),element_boundary=True)
ff.Assemble()
dom_ind_average.vec.data = mf.mat.Inverse() * ff.vec
interface_indicator = BitArray(ffes.ndof)
for i in range(len(dom_ind_average.vec)):
if (dom_ind_average.vec[i] > 1.75) or (dom_ind_average.vec[i] < 1.25):
interface_indicator[i] = False
else:
interface_indicator[i] = True
print(sum(interface_indicator))
return interface_indicator
def GetFaceIndicator(mesh):
dom_indicator = GridFunction(L2(mesh,order=0))
dom_indicator.Set(CoefficientFunction([1,2]))
ffes = FacetFESpace(mesh,order=0)
dom_ind_average = GridFunction(ffes)
mf = BilinearForm(ffes,symmetric=True)
mf += SymbolicBFI(ffes.TrialFunction()ffes.TestFunction(),element_boundary=True)
mf.Assemble()
ff = LinearForm(ffes)
ff += SymbolicLFI(dom_indicatorffes.TestFunction(),element_boundary=True)
ff.Assemble()
dom_ind_average.vec.data = mf.mat.Inverse() * ff.vec
stokes_indicator = BitArray(ffes.ndof)
darcy_indicator = BitArray(ffes.ndof)
for i in range(len(dom_ind_average.vec)):
if (dom_ind_average.vec[i] > 1.75): # exclude darcy part
stokes_indicator[i] = False
else:
stokes_indicator[i] = True
if (dom_ind_average.vec[i] < 1.25): # exclude darcy part
darcy_indicator[i] = False
else:
darcy_indicator[i] = True
#print(“facets marked as stokes facets:”, stokes_indicator)
#print(“facets marked as darcy facets:”, darcy_indicator)
return stokes_indicator, darcy_indicator
def tang(vec):
return vec - (vec*n)*n
stokesId = [1,0]
darcyId = [0,1]
nu = 1.
kinv = 1.
gamma = (1.+4.pipi)/2.
hmax = 1.0/8
u1 = sin(pix)exp(y/2.)
u2 = cos(pix)exp(y/2.)
uexS = CoefficientFunction((-1./2./pi/piu1,1./piu2))
uexD = CoefficientFunction((-2.u1,1./piu2))
pexS = CoefficientFunction(-1./piu2)
pexD = CoefficientFunction(-2./piu2)
cst1 = 1.+nu*(-.5+1./8./pi/pi)
cst2 = -.5/pi+nu*(pi-.25/pi)
cst3 = 0.5/pi-2.pi
stokesSource = CoefficientFunction((cst1u1,cst2u2))
darcySource = CoefficientFunction(cst3u2)
dcFlag = False
n = specialcf.normal(2)
h = specialcf.mesh_size
order = 1
orderp = order-1
penalty = 10order(order+1)/h
InitialMeshSize = 1
for i in range(maxit):
hmax = InitialMeshSize/2**i
mesh = Mesh(topBottom().GenerateMesh (maxh=hmax))
dom1 = GridFunction(L2(mesh,order=0))
dom2 = GridFunction(L2(mesh,order=0))
dom1.Set(CoefficientFunction(stokesId))
dom2.Set(CoefficientFunction(darcyId))
V1 = HDiv(mesh, order=order, dirichlet="gammaS")
FIXME define trace space only on stokes domain
V2 = FESpace ("vectorfacet", mesh, order=order, dirichlet="gammaS",definedon = [1], flags
={"highest_order_dc": dcFlag})
V3 = L2(mesh, order=orderp)
V = FESpace([V1,V2, V3], flags={"dgjumps": True})
u,uhat, p = V.TrialFunction()
v,vhat, q = V.TestFunction()
crossTerm = -div(u)*q-div(v)*p
du = grad(u)
dv = grad(v)
gradu = CoefficientFunction ( (du,), dims=(2,2) )
gradv = CoefficientFunction ( (dv,), dims=(2,2) )
epsu = 0.5*(gradu+gradu.trans)
epsv = 0.5*(gradv+gradv.trans)
interface_indicator = GetInterfaceIndicator(mesh)
a = BilinearForm(V)
dg stokes
stokesCell = 2*nu*InnerProduct(epsu,epsv)
stokesFacet1 = 2*nu*InnerProduct ( epsu*n, tang(vhat-v) )
stokesFacet2 = 2*nu*InnerProduct ( epsv*n, tang(uhat-u) )
stokesFacet3 = penalty*nu*InnerProduct (tang(vhat-v),tang(uhat-u) )
Hdiv darcy
darcyCell = kinv*u*v
interface (only stokes part is needed)
interfaceTerm = gamma*sqrt(kinv)*(tang(uhat)*tang(vhat))
DEBUGGING
Stokes is correct
a += SymbolicBFI ( stokesCell+dom2*darcyCell + crossTerm, definedon=[1])
a += SymbolicBFI ( darcyCell + crossTerm, definedon = [2])
a += SymbolicBFI ( (stokesFacet1+stokesFacet2+stokesFacet3), element_boundary=True, definedon=[1])
interfacebfi = SymbolicBFI (interfaceTerm, VOL, skeleton=True)
interfacebfi.SetDefinedOnElements(interface_indicator)
a += interfacebfi
a.Assemble()
c = Preconditioner(type="direct", bf=a, flags = { "inverse" : "pardiso" } )
f = LinearForm(V)
f += SymbolicLFI(dom1*stokesSource*v-dom2*darcySource*q)
f += -pexD*(v.Trace()*n)*ds(definedon=mesh.Boundaries("gammaD"))
f += -pexS*(v.Trace()*n)*ds(definedon=mesh.Boundaries("gammaS"))
f.Assemble()
u = GridFunction(V)
u.components[0].Set(dom1*uexS,definedon=mesh.Boundaries("gammaS"))
u.components[1].Set(uexS,definedon=mesh.Boundaries("gammaS"))
bvp = BVP(bf=a,lf=f,gf=u,pre=c, maxsteps=1000, prec=1e-11)
errS = dom1*(u.components[0]-uexS)*(u.components[0]-uexS)
errD = dom2*(u.components[0]-uexD)*(u.components[0]-uexD)
erru1 = sqrt( Integrate(errS,mesh,order = 2*order) )
erru2 = sqrt( Integrate(errD,mesh,order = 2*order) )
errorMatrix[0, i] = erru1
errorMatrix[1, i] = erru2
show_error_table(NDof, errorType, errorMatrix)
plt.show()