Dear all,
I have implemented standard finite element code for coupled Stokes-Darcy problem (Here I am attaching the link of the pdf where the model problem has been described [https://drive.google.com/drive/folders/1yjBxjZfF5rSFCa0V0SbETqHlfIeIje4Y?usp=drive_link]). It is running well. But Discontinuous Galerkin(DG) SIPG code for aforementioned problem doesn’t not providing desired results. The problem contains interface condition which occurs coupling between Stokes and Darcy equation. But I think the DG code is not approximating the interface properly and taking zero value at the interface. But in my model problem the specific boundary condition at the interface is not prescribed. Below, I have attached the DG code for the Stokes-Darcy problem.
real nu = 1./1.;
verbosity = 0;
macro div(ux,uy) (dx(ux) + dy(uy))//
int NN = 32;
border D1(t=0, 1){x=t; y=0; label=1;}
border D2(t=0, 1){x=1; y=t; label=2;}
border D3(t=0, 1){x=1-t; y=1; label=3;}//interface of darcy region
border D4(t=0, 1){x=0; y=1-t; label=4;}
border N1(t=0, 1){x=t; y=1; label=5;} //interface of Stokes region
border N2(t=1, 2){x=1; y=t; label=6;}
border N3(t=0, 1){x=1-t; y=2; label=7;}
border N4(t=0, 1){x=0; y=2-t; label=8;}
mesh Th1 = buildmesh( N1(NN) + N2(NN)+ N3(NN) + N4(NN));
mesh Th2 = buildmesh( D1(NN) + D2(NN) + D3(NN) + D4(NN));
mesh Th = Th1+Th2; //mesh
plot(Th);
fespace Vh(Th1,P1dc); // Discontinous P1 finite element
Vh ux, uy, uu1, uu2, dux,duy, vx, vy;
fespace Ph1(Th1, P0); // Discontinous P0 finite element
Ph1 p1, pp1,q1,dp1;
fespace Ph2(Th2, P1dc);// Discontinous P1 finite element
Ph2 p2, pp2,q2,dp2;
real pena=50.; // a paramater to add penalisation
Vh Uex1= x^2*(y-1)^2+y; //exact solution
Vh Uex2= -2x(y-1)^3/(3)+2-pisin(pix); //exact solution
Ph1 p1ex= (2-pisin(pix))sin(piy/(2)); //exact solution
Ph2 p2ex= (2-pisin(pix))(1-y-cos(piy));//exact solution
func U1b= x^2*(y-1)^2+y; //boundary condition
func U2b= -2x(y-1)^3/(3)+2-pisin(pix); //boundary condition
func P22b= (2-pisin(pix))(1-y-cos(piy)); //boundary condition
// Define right hand side for Stokes and Darcy
func f1 = (-2x^2-2(y-1)^2-(pi^2)cos(pix)sin(piy/(2.)));
func f2 = (4x(y-1)-(pi^3)sin(pix)+(2-pisin(pix))(pi/2)cos(piy/(2.)));
func f3 = -(pi^3)sin(pix)(1-y-cos(piy))-(2-pisin(pix))(pi^2)cos(piy);
solve Stokes([dux,duy,dp1], [vx,vy,q1])=
int2d(Th1)(
nu*(dx(dux)dx(vx)
+ dy(dux)dy(vx)
+ dx(duy)dx(vy)
+ dy(duy)dy(vy))
- div(dux,duy) * q1
- dp1div(vx,vy)
- 1e-4dp1q1)
+ intalledges(Th1)(
(nu(mean(dx(dux))N.xjump(vx)+mean(dy(dux))N.yjump(vx)+mean(dx(duy))N.xjump(vy) +mean(dy(duy))N.yjump(vy) +mean(dx(vx))N.xjump(dux)+mean(dy(vx))N.yjump(dux)+mean(dx(vy))N.xjump(duy)+mean(dy(vy))N.yjump(duy)
+ (pena/lenEdge)(jump(dux)jump(vx)+jump(duy)jump(vy)))
-(jump(vx)N.x+jump(vy)N.y)mean(dp1)
-(jump(dux)N.x+jump(duy)N.y)mean(q1)
))
+int1d(Th1,6,7,8)(nu(
(dx(vx))N.x(U1b)+(dy(vx))N.y(U1b)+(dx(vy))N.x(U2b)+(dy(vy))N.y(U2b)
- (pena/lenEdge)((U1b)(vx)+(U2b)(vy)))
-((U1b)N.x+(U2b)N.y)(q1) )
+int1d(Th1,5)((dp2(vx0+vy(-1)))) //interface
+int1d(Th1,5)((dux(1)+duy0)(vx*(1)+vy0)) //interface
-int2d(Th1)([f1,f2]'[vx,vy]);
solve DAR(dp2, q2)=
int2d(Th2)(
dx(dp2)*dx(q2)
+ dy(dp2)*dy(q2)
)
+ intalledges(Th2)((
((mean(dx(dp2))*N.x*jump(q2)+mean(dy(dp2))*N.y*jump(q2)
+mean(dx(q2))*N.x*jump(dp2)+mean(dy(q2))*N.y*jump(dp2)
+ (pena/lenEdge)*(jump(dp2)*jump(q2))) /nTonEdge)
))
+int1d(Th2,1,2,4)((
(dx(q2))*N.x*(P22b)+(dy(q2))*N.y*(P22b)
- (pena/lenEdge)*((P22b)*(q2)))
)
-int1d(Th2,3)((dux*0+duy*(-1))*q2) //interface
-int2d(Th2)(f3*q2)
;
plot([Uex1,Uex2],wait=1, cmm="exact solution");
plot([dux,duy],wait=1, cmm="approximate solution");
plot(p2ex,wait=1, cmm="exact solution");
plot(dp2,wait=1, cmm="approximate solution");
Please, help me out to resolve this problem. Thanks in advance.