[CIG-CS] Sinker example in Matlab
Walter Landry
walter at geodynamics.org
Wed Jul 27 14:42:17 PDT 2011
Hello Everyone,
Recently, a paper came out that looked into discretization errors for
the finite difference method
Duretz, T., D. A. May, T. V. Gerya, and P. J. Tackley (2011),
Discretization errors and free surface stabilization in the finite
difference and marker-in-cell method for applied geodynamics: A
numerical study, Geochem. Geophys. Geosyst., 12, Q07004,
doi:10.1029/2011GC003567.
One benchmark that they did was a circular inclusion under pure shear.
This usually has the same problem as we saw with the sinking block
benchmark, where there are consistent overshoots and undershoots of
the pressure near the material interface.
However, Duretz et. al. claimed to get a convergent solution, so I
thought that they might have gotten around this issue. Reading the
paper, the only difference that I could see that really matters
between my finite difference code and theirs is that they used a
direct solver.
Modifying my code to use a direct solver would be difficult.
Fortunately, as part of Taras Gerya's book, he also distributes Matlab
code that solves Stokes using finite differences and a direct solver.
He even solves a version of the sinker benchmark. So I was able to
modify that to reproduce the form of the sinker benchmark that I had
applied to the deal.II and Samrai codes. I have attached a picture
showing the instantaneous solution for the pressure.
I still see the overshoot and undershoot. So I am not sure what is
going on with Duretz et. al.
In any case, one result of all this is that I now have a Matlab code
that solves the instantaneous sinker problem. I am attaching it in
case anyone is interested. I have only used it in Octave, though it
should work in Matlab.
Cheers,
Walter Landry
walter at geodynamics.org
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Matlab_2D_pressure.png
Type: image/png
Size: 12322 bytes
Desc: not available
Url : http://geodynamics.org/pipermail/cig-cs/attachments/20110727/596925bb/attachment-0001.png
-------------- next part --------------
% Mechanical benchmark for a falling square block;
% solution of 2D Stokes continuity and advection equations
% with finite differences and marker-in-cell technique
% using external function
% Stokes_Continuity_solver_ghost.m
% Staggered Grid with external velocity ghost nodes
%
% vx vx vx
%
% vy +---vy---+---vy---+ vy
% | | |
% vx P vx P vx
% | | |
% vy +---vy---+---vy---+ vy
% | | |
% vx P vx P vx
% | | |
% vy +---vy---+---vy---+ vy
%
% vx vx vx
%
% Lines show basic grid
% Basic (density) nodes are shown with +
% Ghost nodes shown outside the basic grid
% are used for boundary conditions
% Clearing all variables and arrays
clear;
% Clearing figures
clf;
% Acceleration of Gravity, m/s^2
g=1;
% Rock density (kg/m3), viscosity (Pa s)
% Medium
MRHO(1)=0.00;
META(1)=1;
% Block
MRHO(2)=0.03;
META(2)=1e+2;
% Model size, m
xsize=1;
ysize=1;
% Velocity Boundary condition specified by bleft,bright,btop,bbot
% (1=free slip -1=no slip) are implemented from ghost nodes
% directly into Stokes and continuity equations
bleft=1;
bright=1;
btop=1;
bbottom=1;
% Pressure in the upermost, leftmost (first) cell
prfirst=0;
% Defining resolution
xnum=63;
ynum=63;
% Defining gridsteps
xstp=xsize./(xnum-1);
ystp=ysize./(ynum-1);
% Defining number of markers and steps between them in the horizontal and vertical direction
xmx=5; %number of markers per cell in horizontal direction
ymy=5; %number of markers per cell in vertical direction
mxnum=(xnum-1)*xmx; %total number of markers in horizontal direction
mynum=(ynum-1)*ymy; %total number of markers in vertical direction
mxstep=xsize/mxnum; %step between markers in horizontal direction
mystep=ysize/mynum; %step between markers in vertical direction
% Creating markers arrays
MX=zeros(mynum,mxnum); % X coordinate
MY=zeros(mynum,mxnum); % Y coordinate
MI=zeros(mynum,mxnum); % Type
% Defining intial position of markers
% Defining lithological structure of the model
for xm = 1:1:mxnum
for ym = 1:1:mynum
MX(ym,xm)=xm*mxstep-mxstep/2;
MY(ym,xm)=ym*mystep-mystep/2;
MI(ym,xm)=1;
% Density, viscosity structure definition for block
% Relative distances for the marker inside the grid
dx=MX(ym,xm)/xsize;
dy=MY(ym,xm)/ysize;
if(dx>=0.4 && dx<=0.6 && dy>=0.1 && dy<=0.3)
MI(ym,xm)=2;
end
end
end
% Rock type, density and viscosity arrays
typ1 = zeros(ynum,xnum);
etas1 = zeros(ynum,xnum);
etan1 = zeros(ynum-1,xnum-1);
rho1 = zeros(ynum,xnum);
% Backup rock type, density and viscosity arrays
typ0 = typ1;
etas0 = etas1;
etan0 = etan1;
rho0 = rho1;
% Clear rock type, density and viscosity arrays
typ1 = zeros(ynum,xnum);
etas1 = zeros(ynum,xnum);
etan1 = zeros(ynum-1,xnum-1);
rho1 = zeros(ynum,xnum);
% Clear wights for basic nodes
wtnodes=zeros(ynum,xnum);
% Clear wights for etas
wtetas=zeros(ynum,xnum);
% Clear wights for etan
wtetan=zeros(ynum-1,xnum-1);
% Interpolating parameters from markers to nodes
for xm = 1:1:mxnum
for ym = 1:1:mynum
% xn rho(xn,yn)--------------------rho(xn+1,yn)
% ? ^ ?
% ? ? ?
% ? dy ?
% ? ? ?
% ? v ?
% ?<----dx--->o Mrho(xm,ym) ?
% ? ?
% ? ?
% xn+1 rho(xn,yn+1)-------------------rho(xn+1,yn+1)
%
% Define indexes for upper left node in the cell where the marker is
% !!! SUBTRACT 0.5 since int16(0.5)=1
xn=double(int16(MX(ym,xm)./xstp-0.5))+1;
yn=double(int16(MY(ym,xm)./ystp-0.5))+1;
if (xn<1)
xn=1;
end
if (xn>xnum-1)
xn=xnum-1;
end
if (yn<1)
yn=1;
end
if (yn>ynum-1)
yn=ynum-1;
end
% Define normalized distances from marker to the upper left node;
dx=MX(ym,xm)./xstp-xn+1;
dy=MY(ym,xm)./ystp-yn+1;
% Add density to 4 surrounding nodes
rho1(yn,xn)=rho1(yn,xn)+(1.0-dx).*(1.0-dy).*MRHO(MI(ym,xm));
wtnodes(yn,xn)=wtnodes(yn,xn)+(1.0-dx).*(1.0-dy);
rho1(yn+1,xn)=rho1(yn+1,xn)+(1.0-dx).*dy.*MRHO(MI(ym,xm));
wtnodes(yn+1,xn)=wtnodes(yn+1,xn)+(1.0-dx).*dy;
rho1(yn,xn+1)=rho1(yn,xn+1)+dx.*(1.0-dy).*MRHO(MI(ym,xm));
wtnodes(yn,xn+1)=wtnodes(yn,xn+1)+dx.*(1.0-dy);
rho1(yn+1,xn+1)=rho1(yn+1,xn+1)+dx.*dy.*MRHO(MI(ym,xm));
wtnodes(yn+1,xn+1)=wtnodes(yn+1,xn+1)+dx.*dy;
% Add shear viscosity etas() and rock
% type typ() to 4 surrounding nodes only
% using markers located at <=0.5
% gridstep distances from nodes
if(dx<=0.5 && dy<=0.5)
etas1(yn,xn)=etas1(yn,xn)+(1.0-dx).*(1.0-dy).*META(MI(ym,xm));
typ1(yn,xn)=typ1(yn,xn)+(1.0-dx).*(1.0-dy).*MI(ym,xm);
wtetas(yn,xn)=wtetas(yn,xn)+(1.0-dx).*(1.0-dy);
end
if(dx<=0.5 && dy>=0.5)
etas1(yn+1,xn)=etas1(yn+1,xn)+(1.0-dx).*dy.*META(MI(ym,xm));
typ1(yn+1,xn)=typ1(yn+1,xn)+(1.0-dx).*dy.*MI(ym,xm);
wtetas(yn+1,xn)=wtetas(yn+1,xn)+(1.0-dx).*dy;
end
if(dx>=0.5 && dy<=0.5)
etas1(yn,xn+1)=etas1(yn,xn+1)+dx.*(1.0-dy).*META(MI(ym,xm));
typ1(yn,xn+1)=typ1(yn,xn+1)+dx.*(1.0-dy).*MI(ym,xm);
wtetas(yn,xn+1)=wtetas(yn,xn+1)+dx.*(1.0-dy);
end
if(dx>=0.5 && dy>=0.5)
etas1(yn+1,xn+1)=etas1(yn+1,xn+1)+dx.*dy.*META(MI(ym,xm));
typ1(yn+1,xn+1)=typ1(yn+1,xn+1)+dx.*dy.*MI(ym,xm);
wtetas(yn+1,xn+1)=wtetas(yn+1,xn+1)+dx.*dy;
end
% Add normal viscosity etan() to the
% center of current cell
etan1(yn,xn)=etan1(yn,xn)+(1.0-abs(0.5-dx)).*(1.0-abs(0.5-dy)).*META(MI(ym,xm));
wtetan(yn,xn)=wtetan(yn,xn)+(1.0-abs(0.5-dx)).*(1.0-abs(0.5-dy));
end
end
% Computing Viscosity, density, rock type for nodal points
for i=1:1:ynum;
for j=1:1:xnum;
% Density
if (wtnodes(i,j)~=0)
% Compute new value interpolated from markers
rho1(i,j)=rho1(i,j)./wtnodes(i,j);
else
% If no new value is interpolated from markers old value is used
rho1(i,j)=rho0(i,j);
end
% Shear viscosity and type
if (wtetas(i,j)~=0)
% Compute new value interpolated from markers
etas1(i,j)=etas1(i,j)./wtetas(i,j);
typ1(i,j)=typ1(i,j)./wtetas(i,j);
else
% If no new value is interpolated from markers old value is used
etas1(i,j)=etas0(i,j);
typ1(i,j)=typ0(i,j);
end
% Normal viscosity
if (i<ynum && j<xnum)
if (wtetan(i,j)~=0)
% Compute new value interpolated from markers
etan1(i,j)=etan1(i,j)./wtetan(i,j);
else
% If no new value is interpolated from markers old value is used
etan1(i,j)=etan0(i,j);
end
end
end
end
% Computing right part of Stokes (RX, RY) and Continuity (RC) equation
% vx, vy, P
vx1=zeros(ynum+1,xnum);
vy1=zeros(ynum,xnum+1);
pr1=zeros(ynum-1,xnum-1);
% Right parts of equations
RX1=zeros(ynum+1,xnum);
RY1=zeros(ynum,xnum+1);
RC1=zeros(ynum-1,xnum-1);
% Grid points cycle
for i=1:1:ynum;
for j=1:1:xnum;
% Right part of x-Stokes Equation
if(j>1 && i>1 && j<xnum)
RX1(i,j)=0;
end
% Right part of y-Stokes Equation
if(j>1 && i>1 && i<ynum)
RY1(i,j)=-g*(rho1(i,j)+rho1(i,j-1))/2;
end
end
end
% Solving of Stokes and Continuity equations on nodes
% and computing residuals
% by calling function Stokes_Continuity_solver_ghost()
[vx1,resx1,vy1,resy1,pr1,resc1]=Stokes_Continuity_solver_ghost(prfirst,etas1,etan1,xnum,ynum,xstp,ystp,RX1,RY1,RC1,bleft,bright,btop,bbottom);
% Defining scale for Stokes residuals from y-Stokes equation
% dSIGMAij/dj-dP/di=-RHO*gi=0 => Stokes scale=abs(RHO*gi)
stokesscale= MRHO(1)*g;
% Defining scale for Continuity residuals from y-Stokes equation
% dvx/dx+dvy/dy=0 can be transformed to 2ETA(dvx/dx+dvy/dy)/dx=0
% which is similar to dSIGMAij/dj and has scale given above
% therefore continuity scale = scale=abs(RHO*gi/ETA*dx)
continscale= MRHO(1)*g/META(1)*xstp;
figure(1), clf;
pcolor(pr1);
# shading interp;
axis tight;
zlabel('Pressure')
colorbar();
# % Plotting Residuals for x-Stokes as surface
# figure(1), clf;
# subplot(2,3,1)
# resx0=resx1/stokesscale;
# surf(resx0);
# shading interp;
# axis tight;
# zlabel('residual x-Stokes')
# % colorbar;
# % Plotting Residuals for y-Stokes as surface
# subplot(2,3,2)
# resy0=resy1/stokesscale;
# surf(resy0);
# shading interp;
# axis tight;
# zlabel('residual Y-stokes')
# % colorbar;
# % Plotting Residuals for Continuity as surface
# subplot(2,3,3)
# resc0=resc1/continscale;
# surf(resc0);
# shading interp;
# axis tight;
# zlabel('residual continuity')
# % colorbar;
# % Plotting vx
# subplot(2,3,4)
# pcolor(vx1);
# shading interp;
# axis ij;
# title('Vx, m/s')
# colorbar;
# % Plotting vy
# subplot(2,3,5)
# pcolor(vy1);
# shading interp;
# axis ij;
# title('Vy, m/s')
# colorbar;
# % Plotting density
# subplot(2,3,6)
# pcolor(pr1);
# shading interp;
# axis ij;
# colorbar;
# title(['Density (kg/m^3) Step=',num2str(0),' Myr=',num2str(0)]);
# % Stop for 0.1 second
drawnow % draw above figure without delay
-------------- next part --------------
% Function Stokes_Continuity_solver_ghost()
% This function formulates and solves
% Stokes and Continuity equations defined on 2D regularly spaced staggered grid
% with specified resolution (xnum, ynum) and gridsteps (xstp, ystp)
% given distribution of right parts for all equations (RX,RY,RC) on the grid
% and given variable shear (etas) and normal (etan) viscosity distributions
% pressure is normalized relative to given value (prnorm) in the first cell
%
% Velocity Boundary condition specified by bleft,bright,btop,bbot
% (1=free slip -1=no slip) are implemented from ghost nodes
% directly into Stokes and continuity equations
%
% Function returns solution for velocity and pressure (vx,vy,pr)
% and distribution of residuals (resx,resy,resc)
function[vx,resx,vy,resy,pr,resc]=Stokes_Continuity_solver_ghost(prnorm,etas,etan,xnum,ynum,xstp,ystp,RX,RY,RC,bleft,bright,btop,bbottom)
%
% Staggered Grid with external velocity ghost nodes
%
% vx vx vx
%
% vy +---vy---+---vy---+ vy
% | | |
% vx P vx P vx
% | | |
% vy +---vy---+---vy---+ vy
% | | |
% vx P vx P vx
% | | |
% vy +---vy---+---vy---+ vy
%
% vx vx vx
%
% Lines show basic grid
% Basic (density) nodes are shown with +
% Ghost nodes shown outside the basic grid
% are used for boundary conditions
% Poisson-like equations koefficients
xkf=1/xstp^2;
xkf2=2/xstp^2;
ykf=1/ystp^2;
ykf2=2/ystp^2;
xykf=1/(xstp*ystp);
% Koefficient for scaling pressure
pscale=2*etan(1)/(xstp+ystp);
% Horizontal shift index
ynum3=(ynum-1)*3;
% Creating matrix
L=sparse((xnum-1)*(ynum-1)*3,(xnum-1)*(ynum-1)*3);
R=zeros((xnum-1)*(ynum-1)*3,1);
% Solving of Stokes and continuity equations on nodes
for i=1:1:ynum-1
for j=1:1:xnum-1
% Indexes for P,vx,vy
ivx=((j-1)*(ynum-1)+(i-1))*3+1;
ivy=ivx+1;
ipr=ivx+2;
% x-Stokes equation dSIGMAxx/dx+dSIGMAxy/dy-dP/dx=RX
if (j<xnum-1)
% x-Stokes equation stensil
% +-------------------- -+----------------------+
% | | |
% | | |
% | vx(i-1,j) |
% | | |
% | | |
% +-----vy(i-1,j)---etas(i,j+1)---vy(i-1,j+1)---+
% | | |
% | | |
% vx(i,j-1) pr(i,j) vx(i,j) P(i,j+1) vx(i,j+1)
% | etan(i,j) | etan(i,j+1) |
% | | |
% +------vy(i,j)---etas(i+1,j+1)---vy(i,j+1)----+
% | | |
% | | |
% | vx(i+1,j) |
% | | |
% | | |
% +-------------------- -+----------------------+
% Right Part
R(ivx,1)=RX(i+1,j+1);
% Computing Current x-Stokes coefficients
% Central Vx node
L(ivx,ivx)=-xkf2*(etan(i,j+1)+etan(i,j))-ykf*(etas(i+1,j+1)+etas(i,j+1));
% Left Vx node
if (j>1)
ivxleft=ivx-ynum3;
L(ivx,ivxleft)=xkf2*etan(i,j);
end
% Right Vx node
if (j<xnum-2)
ivxright=ivx+ynum3;
L(ivx,ivxright)=xkf2*etan(i,j+1);
end
% Top Vx node
if (i>1)
ivxtop=ivx-3;
L(ivx,ivxtop)=ykf*etas(i,j+1);
else
L(ivx,ivx)=L(ivx,ivx)+btop*ykf*etas(i,j+1);
end
% Bottom Vx node
if (i<ynum-1)
ivxbottom=ivx+3;
L(ivx,ivxbottom)=ykf*etas(i+1,j+1);
else
L(ivx,ivx)=L(ivx,ivx)+bbottom*ykf*etas(i+1,j+1);
end
% Top Left Vy node
if (i>1)
ivytopleft=ivx-3+1;
L(ivx,ivytopleft)=xykf*etas(i,j+1);
end
% Top Right Vy node
if (i>1)
ivytopright=ivx-3+1+ynum3;
L(ivx,ivytopright)=-xykf*etas(i,j+1);
end
% Bottom Left Vy node
if (i<ynum-1)
ivybottomleft=ivx+1;
L(ivx,ivybottomleft)=-xykf*etas(i+1,j+1);
end
% Bottom Right Vy node
if (i<ynum-1)
ivybottomright=ivx+1+ynum3;
L(ivx,ivybottomright)=xykf*etas(i+1,j+1);
end
% Left P node
iprleft=ivx+2;
L(ivx,iprleft)=pscale/xstp;
% Right P node
iprright=ivx+2+ynum3;
L(ivx,iprright)=-pscale/xstp;
% Ghost Vx_parameter=0 used for numbering
else
L(ivx,ivx)=1;
R(ivx,1)=0;
end
% y-Stokes equation dSIGMAyy/dy+dSIGMAyx/dx-dP/dy=RY
if (i<ynum-1)
% y-Stokes equation stensil
% +-------------------- -+-------vy(i-1,j)------+----------------------+
% | | | |
% | | | |
% | vx(i,j-1) P(i,j) vx(i,j) |
% | | etan(i,j) | |
% | | | |
% +-----vy(i,j-1)---etas(i+1,j)---vy(i,j)--etas(i+1,j+1)---vy(i,j+1)---+
% | | | |
% | | | |
% | vx(i+1,j-1) P(i+1,j) vx(i+1,j) |
% | | etan(i+1,j) | |
% | | | |
% +----------------------+-------vy(i+1,j)------+----------------------+
%
% Right Part
R(ivy,1)=RY(i+1,j+1);
% Computing Current y-Stokes coefficients
% Central Vy node
L(ivy,ivy)=-ykf2*(etan(i+1,j)+etan(i,j))-xkf*(etas(i+1,j+1)+etas(i+1,j));
% Top Vy node
if(i>1)
ivytop=ivy-3;
L(ivy,ivytop)=ykf2*etan(i,j);
end
% Bottom Vy node
if(i<ynum-2)
ivybottom=ivy+3;
L(ivy,ivybottom)=ykf2*etan(i+1,j);
end
% Left Vy node
if(j>1)
ivyleft=ivy-ynum3;
L(ivy,ivyleft)=xkf*etas(i+1,j);
else
L(ivy,ivy)=L(ivy,ivy)+bleft*xkf*etas(i+1,j);
end
% Right Vy node
if(j<xnum-1)
ivyright=ivy+ynum3;
L(ivy,ivyright)=xkf*etas(i+1,j+1);
else
L(ivy,ivy)=L(ivy,ivy)+bright*xkf*etas(i+1,j+1);
end
% Top left Vx node
if (j>1)
ivxtopleft=ivy-1-ynum3;
L(ivy,ivxtopleft)=xykf*etas(i+1,j);
end
% Bottom left Vx node
if (j>1)
ivxbottomleft=ivy-1+3-ynum3;
L(ivy,ivxbottomleft)=-xykf*etas(i+1,j);
end
% Top right Vx node
if (j<xnum-1)
ivxtopright=ivy-1;
L(ivy,ivxtopright)=-xykf*etas(i+1,j+1);
end
% Bottom right Vx node
if (j<xnum-1)
ivxbottomright=ivy-1+3;
L(ivy,ivxbottomright)=xykf*etas(i+1,j+1);
end
% Top P node
iprtop=ivy+1;
L(ivy,iprtop)=pscale/ystp;
% Bottom P node
iprbottom=ivy+1+3;
L(ivy,iprbottom)=-pscale/ystp;
% Ghost Vy_parameter=0 used for numbering
else
L(ivy,ivy)=1;
R(ivy,1)=0;
end
% Continuity equation dvx/dx+dvy/dy=RC
if (j>1 || i>1)
% Continuity equation stensil
% +-----vy(i-1,j)--------+
% | |
% | |
% vx(i,j-1) pr(i,j) vx(i,j)
% | |
% | |
% +------vy(i,j)---------+
%
% Right Part
R(ipr,1)=RC(i,j);
% Computing Current Continuity coefficients
% Left Vx node
if (j>1)
ivxleft=ipr-2-ynum3;
L(ipr,ivxleft)=-pscale/xstp;
end
% Right Vx node
if (j<xnum-1)
ivxright=ipr-2;
L(ipr,ivxright)=pscale/xstp;
end
% Top Vy node
if (i>1)
ivytop=ipr-1-3;
L(ipr,ivytop)=-pscale/ystp;
end
% Bottom Vy node
if (i<ynum-1)
ivybottom=ipr-1;
L(ipr,ivybottom)=pscale/ystp;
end
% Pressure definition for the upper left node
else
L(ipr,ipr)=2*pscale/(xstp+ystp);
R(ipr,1)=2*prnorm/(xstp+ystp);
end
end
end
% Solve matrix
S=L\R;
% Reload solution
for i=1:1:ynum-1
for j=1:1:xnum-1
% Indexes for P,vx,vy
ivx=((j-1)*(ynum-1)+(i-1))*3+1;
ivy=ivx+1;
ipr=ivx+2;
% Reload Vx
if (j<xnum-1)
vx(i+1,j+1)=S(ivx);
end
% Reload Vy
if (i<ynum-1)
vy(i+1,j+1)=S(ivy);
end
% Reload P
pr(i,j)=S(ipr)*pscale;
end
end
% Apply vx boundary conditions
% Left Boundary
vx(:,1)=0;
% Right Boundary
vx(:,xnum)=0;
% Upper Boundary
vx(1,:)=btop*vx(2,:);
% Lower Boundary
vx(ynum+1,:)=bbottom*vx(ynum,:);
% Apply vy boundary conditions
% Upper Boundary
vy(1,:)=0;
% Lower Boundary
vy(ynum,:)=0;
% Left Boundary
vy(:,1)=bleft*vy(:,2);
% Right Boundary
vy(:,xnum+1)=bright*vy(:,xnum);
% Computing residuals
for i=1:1:ynum+1;
for j=1:1:xnum+1;
% x-Stokes equation dSIGMAxx/dx+dSIGMAxy/dy-dP/dx=RX
if (j<xnum+1)
% vx-Boundrary conditions
if (i==1 || i==ynum+1 || j==1 || j==xnum);
resx(i,j)=0;
% x-Stokes equation
else
% Computing Current x-Stokes residual
% dSIGMAxx/dx-dP/dx
resx(i,j)=RX(i,j)-(xkf2*(etan(i-1,j)*(vx(i,j+1)-vx(i,j))-etan(i-1,j-1)*(vx(i,j)-vx(i,j-1)))-(pr(i-1,j)-pr(i-1,j-1))/xstp);
% dSIGMAxy/dy
resx(i,j)=resx(i,j)-(etas(i,j)*((vx(i+1,j)-vx(i,j))*ykf+(vy(i,j+1)-vy(i,j))*xykf)-etas(i-1,j)*((vx(i,j)-vx(i-1,j))*ykf+(vy(i-1,j+1)-vy(i-1,j))*xykf));
end
end
% y-Stokes equation dSIGMAyy/dy+dSIGMAyx/dx-dP/dy=RY
if (i<ynum+1)
% vy-Boundrary conditions
if (i==1 || i==ynum || j==1 || j==xnum+1);
resy(i,j)=0;
%y-Stokes equation
else
% Computing current residual
% dSIGMAyy/dy-dP/dy
resy(i,j)=RY(i,j)-(ykf2*(etan(i,j-1)*(vy(i+1,j)-vy(i,j))-etan(i-1,j-1)*(vy(i,j)-vy(i-1,j)))-(pr(i,j-1)-pr(i-1,j-1))/ystp);
% dSIGMAxy/dx
resy(i,j)=resy(i,j)-(etas(i,j)*((vy(i,j+1)-vy(i,j))*xkf+(vx(i+1,j)-vx(i,j))*xykf)-etas(i,j-1)*((vy(i,j)-vy(i,j-1))*xkf+(vx(i+1,j-1)-vx(i,j-1))*xykf));
end
end
% Continuity equation dvx/dx+dvy/dy=RC
if (i<ynum && j<xnum);
% Computing current residual
resc(i,j)=RC(i,j)-((vx(i+1,j+1)-vx(i+1,j))/xstp+(vy(i+1,j+1)-vy(i,j+1))/ystp);
end
end
end
More information about the CIG-CS
mailing list