Self gravitation:LINE451-486.
%This program complished the process of 2D Eulerian advection with method
%of MIC(allocating memory temporarily).
clear all;clf;
% Set parameters
Nt = 2000;
output_interval = 1;
% Length of model
Lx = 10000000;
Ly = 10000000;
xmid = Lx/2;
ymid = Ly/2;
% Sun of point number of model
Nx = 101;
Ny = 101;
Nx1 = Nx+1;
Ny1 = Ny+1;
% Interval of E-staggered-nodes
dx = Lx / (Nx-1);
dy = Ly / (Ny-1);
% Basic nodes.
x = 0:dx:Lx;
y = 0:dy:Ly;
% Vx staggered nodes.(Nx_node * Ny_node+1)
xvx = 0:dx:Lx;
yvx = -dy/2 : dy : Ly+dy/2;
% Vy staggered nodes.(Nx_node+1 * Ny_node)
xvy = -dx/2 : dx : Lx+dx/2;
yvy = 0:dy:Ly;
% P staggered nodes.(Nx_node+1 * Ny_node+1)
xp = -dx/2 : dx : Lx+dx/2;
yp = -dy/2 : dy : Ly+dy/2;
% Velocity and Pressure matrix
vx_matrix = zeros(Ny1,Nx1);
vy_matrix = zeros(Ny1,Nx1);
p_matrix = zeros(Ny1,Nx1);
% Vx at Pnodes
pvx = zeros(Ny1,Nx1);
pvy = zeros(Ny1,Nx1);
% Gravity acceleration
g_y = zeros(Ny1,Nx1);
g_x = zeros(Ny1,Nx1);
ETA = zeros(Ny,Nx);%density
RHO = zeros(Ny,Nx);%viscosity
% Stokes-equ and continuity-equ
% Temperature-equ
% Amount of markers.
Nx_marker = 4*(Nx-1);
Ny_marker = 4*(Ny-1);
marknum=Nx_marker*Ny_marker;
% Step-length among L-marker points
x_m_step = Lx/Nx_marker;
y_m_step = Ly/Ny_marker;
% Marker points.
% Stokes-equ and continuity-equ
xm=zeros(1,marknum); % Horizontal coordinates, m
ym=zeros(1,marknum);
rhom=zeros(1,marknum); % Density, kg/m^3
etam=zeros(1,marknum);
% Temperature-equ
Km=zeros(1,marknum); % K, kg/m^3 kt
DenCpm=zeros(1,marknum); % den*Cp hct
tdm = zeros(1,marknum); % K/den/Cp td
Tm0 = zeros(1,marknum); % temperature
Aefm = zeros(1,marknum); % K/den/Cp td
Hrm = zeros(1,marknum); % temperature
% Define matrix
hct = zeros(Ny1,Nx1); %Den*Cp
kht = zeros(Ny1,Nx1); %K
T0 = zeros(Ny1,Nx1); %T
td = zeros(Ny1,Nx1); %K/Den*Cp
bt = zeros(Nx1*Ny1, 1);
coet = zeros(Nx1*Ny1, Nx1*Ny1);
bg = zeros(Nx1*Ny1, 1);
coeg = zeros(Nx1*Ny1, Nx1*Ny1);
% Setup of model
% "Squared planet"
den_mantle = 3300;
vis_mantle = 1e+21;
k_mantle = 3;%W/(m*K)
Cp_mantle = 1000;%J/(Kg*K)
t_mantle = 1500;%K
aef_mantle = 3e-5;%1/K
Hr_mantle = 2e-8;
% "Low density core"
den_plume = 3000;
vis_plume = 1e+17;
k_plume = 2;%W/(m*K)
Cp_plume = 1000;%J/(Kg*K)
t_plume = 1800;%K
aef_plume = 2e-5;%1/K
Hr_plume = 2e-7;
% "Striky space"
den_striky = 0;
vis_striky = 1e+17;
k_striky = 3000;%W/(m*K)
Cp_striky = 1100;%J/(Kg*K)
t_striky = 273;%K
aef_striky = 0;%1/K
Hr_striky = 0;
% radius = 100000;
RHOVY=zeros(Ny1,Nx1); %kg/m^3
RHOVX=zeros(Ny1,Nx1); %kg/m^3
KVY=zeros(Ny1,Nx1); %kg/m^3
KVX=zeros(Ny1,Nx1); %kg/m^3
ETAP = zeros(Ny1,Nx1);
RHOP = zeros(Ny1,Nx1);
% Parameters and matrixs of exercise 9.3
srxy = zeros(Ny1,Nx1);
dsxy = zeros(Ny1,Nx1);
srxx = zeros(Ny,Nx1);
dsxx = zeros(Ny,Nx1);
pxyxy = zeros(Ny1,Nx1);
pxxxx = zeros(Ny1,Nx1);
Hsij = zeros(Ny1,Nx1);
Haij = zeros(Ny1,Nx1);
Hrij = zeros(Ny1,Nx1);
Aefij = zeros(Ny1,Nx1);
m = 1;
for jm=1:1:Nx_marker
for im=1:1:Ny_marker
% Define marker coordinates
xm(m)=x_m_step/2+(jm-1)*x_m_step+(rand-0.5)*x_m_step;
ym(m)=y_m_step/2+(im-1)*y_m_step+(rand-0.5)*y_m_step;
% Marker properties
% "Squared planet"
if(xm(m)>=0.2*Lx && xm(m)<=0.8*Lx && ym(m)>=0.2*Ly && ym(m)<=0.8*Ly)
rhom(m) = den_mantle; % Mantle density mantle == planet
etam(m) = vis_mantle; % Mantle viscosity
Km(m) = k_mantle; % Mantle K
DenCpm(m) = den_mantle*Cp_mantle; % Mantle den*Cp
Tm0(m) = t_mantle; % Mantle T
Aefm(m) = aef_mantle;
Hrm(m) = Hr_mantle;
% "Low density core"
if(xm(m)>=0.35*Lx && xm(m)<=0.65*Lx && ym(m)>=0.35*Ly && ym(m)<=0.65*Ly)
rhom(m) = den_plume; % Plume density
etam(m) = vis_plume; % Plume viscosity
Km(m) = k_plume; % Plume K
DenCpm(m) = den_plume*Cp_plume; % Plume den*Cp
Tm0(m) = t_plume; %Plume T
Aefm(m) = aef_plume;
Hrm(m) = Hr_plume;
end
% "Striky space"
else
rhom(m) = den_striky; % Sticky density
etam(m) = vis_striky; % Sticky viscosity
Km(m) = k_striky; % Sticky K
DenCpm(m) = 3.3e+6; % Sticky den*Cp
Tm0(m) = t_striky; % Sticky T
Aefm(m) = aef_striky;
Hrm(m) = Hr_striky;
end
% Update marker counter
m=m+1;
end
end
Kcont = 1e+21/dx;
% Unknown_number is (Ny+1)*(Nx+1)*3
unknown_number = (Nx+1)*(Ny+1)*3;
coe = sparse(unknown_number,unknown_number);
R = zeros(unknown_number, 1);
% Boundary condition : free slip = -1 ; no slip = 1
bc_left = -1;
bc_right = -1;
bc_top = -1;
bc_bottom = -1;
dxymax = 0.5;
dt = 1e+10; % initial timestep
dtkoef=1.2; % timestep increment
DTmax=20; % max temperature change per time step, K
subgrdifcoe = 1;
figure(1);
for t = 1:Nt
% "drunken sailor" instability appear if dt = 0 at the next line.
% dt = 0;
% BASCI NODE
w_m_x_node = zeros(Ny,Nx);
vis_w_m = zeros(Ny,Nx);
% rho in vy and vx nodes
RHOYSUM=zeros(Ny1,Nx1);
WTYSUM=zeros(Ny1,Nx1);
RHOXSUM=zeros(Ny1,Nx1);
WTXSUM=zeros(Ny1,Nx1);
% K in vy and vx nodes
KYSUM=zeros(Ny1,Nx1);
WKYSUM=zeros(Ny1,Nx1);
KXSUM=zeros(Ny1,Nx1);
WKXSUM=zeros(Ny1,Nx1);
% Interpolate viscosity from markers to BASIC NODES.
% Searching the coordinate of upper-left node of every L-marker.
for m=1:1:marknum
% BASIC NODE Interpolate den and vis from markes
j=fix((xm(m)-x(1))/dx)+1;
i=fix((ym(m)-y(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx-1)
j=Nx-1;
end
if(i<1)
i=1;
elseif(i>Ny-1)
i=Ny-1;
end
% Compute distances
dis_x = xm(m)-x(j);
dis_y = ym(m)-y(i);
wtmij = (1-dis_x/dx)*(1-dis_y/dy);
wtmi1j = (dis_x/dx)*(1-dis_y/dy);
wtmij1 = (1-dis_x/dx)*(dis_y/dy);
wtmi1j1 = (dis_x/dx)*(dis_y/dy);
% | |
% * * * * * Discard the m-point on the right.
w_m_x_node(i,j) = w_m_x_node(i,j) + wtmij;
vis_w_m(i,j) = vis_w_m(i,j) + etam(m).*wtmij;
% i+1,j
w_m_x_node(i+1,j) = w_m_x_node(i+1,j) + wtmi1j;
vis_w_m(i+1,j) = vis_w_m(i+1,j) + etam(m).*wtmi1j;
% i,j+1
w_m_x_node(i,j+1) = w_m_x_node(i,j+1) + wtmij1;
vis_w_m(i,j+1) = vis_w_m(i,j+1) + etam(m).*wtmij1;
% i+1,j+1
w_m_x_node(i+1,j+1) = w_m_x_node(i+1,j+1) + wtmi1j1;
vis_w_m(i+1,j+1) = vis_w_m(i+1,j+1) + etam(m).*wtmi1j1;
% Interpolate Density and K to VY-NODE
% Vy staggered nodes.(Nx1 * Ny)
j=fix((xm(m)-xvy(1))/dx)+1;
i=fix((ym(m)-yvy(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx)
j=Nx;
end
if(i<1)
i=1;
elseif(i>Ny-1)
i=Ny-1;
end
% Compute distances
dxmj=xm(m)-xvy(j);
dymi=ym(m)-yvy(i);
% Compute weights
wtmij=(1-dxmj/dx)*(1-dymi/dy);
wtmi1j=(1-dxmj/dx)*(dymi/dy);
wtmij1=(dxmj/dx)*(1-dymi/dy);
wtmi1j1=(dxmj/dx)*(dymi/dy);
% Update properties
% i,j Node
% rho and k in vy nodes
RHOYSUM(i,j)=RHOYSUM(i,j)+rhom(m)*wtmij;
WTYSUM(i,j)=WTYSUM(i,j)+wtmij;
KYSUM(i,j)=KYSUM(i,j)+Km(m)*wtmij;% k in vy nodes
% i+1,j Node
RHOYSUM(i+1,j)=RHOYSUM(i+1,j)+rhom(m)*wtmi1j;
WTYSUM(i+1,j)=WTYSUM(i+1,j)+wtmi1j;
KYSUM(i+1,j)=KYSUM(i+1,j)+Km(m)*wtmi1j;
% i,j+1 Node
RHOYSUM(i,j+1)=RHOYSUM(i,j+1)+rhom(m)*wtmij1;
WTYSUM(i,j+1)=WTYSUM(i,j+1)+wtmij1;
KYSUM(i,j+1)=KYSUM(i,j+1)+Km(m)*wtmij1;
% i+1,j+1 Node
RHOYSUM(i+1,j+1)=RHOYSUM(i+1,j+1)+rhom(m)*wtmi1j1;
WTYSUM(i+1,j+1)=WTYSUM(i+1,j+1)+wtmi1j1;
KYSUM(i+1,j+1)=KYSUM(i+1,j+1)+Km(m)*wtmi1j1;
% Interpolate Density and K to VX-NODE
% Vx staggered nodes.(Nx * Ny1)
j=fix((xm(m)-xvx(1))/dx)+1;
i=fix((ym(m)-yvx(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx-1)
j=Nx-1;
end
if(i<1)
i=1;
elseif(i>Ny)
i=Ny;
end
% Compute distances
dxmj=xm(m)-xvx(j);
dymi=ym(m)-yvx(i);
% Compute weights
wtmij=(1-dxmj/dx)*(1-dymi/dy);
wtmi1j=(1-dxmj/dx)*(dymi/dy);
wtmij1=(dxmj/dx)*(1-dymi/dy);
wtmi1j1=(dxmj/dx)*(dymi/dy);
% Update properties
% i,j Node
RHOXSUM(i,j)=RHOXSUM(i,j)+rhom(m)*wtmij;
WTXSUM(i,j)=WTXSUM(i,j)+wtmij;
KXSUM(i,j)=KXSUM(i,j)+Km(m)*wtmij;% k in vx nodes
% i+1,j Node
RHOXSUM(i+1,j)=RHOXSUM(i+1,j)+rhom(m)*wtmi1j;
WTXSUM(i+1,j)=WTXSUM(i+1,j)+wtmi1j;
KXSUM(i+1,j)=KXSUM(i+1,j)+Km(m)*wtmi1j;
% i,j+1 Node
RHOXSUM(i,j+1)=RHOXSUM(i,j+1)+rhom(m)*wtmij1;
WTXSUM(i,j+1)=WTXSUM(i,j+1)+wtmij1;
KXSUM(i,j+1)=KXSUM(i,j+1)+Km(m)*wtmij1;
% i+1,j+1 Node
RHOXSUM(i+1,j+1)=RHOXSUM(i+1,j+1)+rhom(m)*wtmi1j1;
WTXSUM(i+1,j+1)=WTXSUM(i+1,j+1)+wtmi1j1;
KXSUM(i+1,j+1)=KXSUM(i+1,j+1)+Km(m)*wtmi1j1;
end
% Computing the value of density and viscosity of BASIC-NODE.
for j = 1:Nx
for i = 1:Ny
if(w_m_x_node(i,j)>0)
% Not need RHO on the basic nodes.Just Y nodes.
ETA(i,j) = vis_w_m(i,j)/w_m_x_node(i,j);
end
end
end
% Computing Density and K value at the VY and VX-nodes.
for j=1:1:Nx1
for i=1:1:Ny1
% VY-nodes:(Nx1 * Ny) absent Ny1=0 !!!
if(WTYSUM(i,j)>0)
% rou_v
RHOVY(i,j)=RHOYSUM(i,j)/WTYSUM(i,j);
% K_v
KVY(i,j)=KYSUM(i,j)/WTYSUM(i,j);
end
% VX-nodes:(Nx * Ny1)absent Nx1=0 !!!
if(WTXSUM(i,j)>0)
% rou_h
RHOVX(i,j)=RHOXSUM(i,j)/WTXSUM(i,j);
% K_h
KVX(i,j) = KXSUM(i,j)/WTXSUM(i,j);
end
end
end
% Interpolate ETA from L-markers to P-nodes.
w_m_p_node = zeros(Ny1,Nx1);% weights except T
ETAp_w_m = zeros(Ny1,Nx1);% Viscosity at P node
RHOp_w_m = zeros(Ny1,Nx1);% Density at P node
Hr_w_m = zeros(Ny1,Nx1); % Hr at P node
Aef_w_m = zeros(Ny1,Nx1); % Aef at P node
K_w_m_t = zeros(Ny1,Nx1); % K at P node
DenCp_w_m_t = zeros(Ny1,Nx1);% Den*Cp at P node
T_w_m_t = zeros(Ny1,Nx1); % T at P node
w_m_t_T = zeros(Ny1,Nx1); % weight of T
for m=1:1:marknum
% Reserve original coordinate.
% Interpolate Hr and aef from L-markers to P-staggered nodes.
% P staggered nodes.(Nx1 * Ny1)
j=fix((xm(m)-xp(1))/dx)+1;
i=fix((ym(m)-yp(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx)
j=Nx;
end
% P staggered nodes.(Nx1 * Ny1)
if(i<1)
i=1;
elseif(i>Ny)
i=Ny;
end
% Compute distances
dis_x=xm(m)-xp(j);
dis_y=ym(m)-yp(i);
% Compute weights
wtmij=(1-dis_x/dx)*(1-dis_y/dy);
wtmi1j=(dis_x/dx)*(1-dis_y/dy);
wtmij1=(1-dis_x/dx)*(dis_y/dy);
wtmi1j1=(dis_x/dx)*(dis_y/dy);
% Compute vx velocity
w_m_p_node(i,j)=w_m_p_node(i,j)+wtmij;
ETAp_w_m(i,j)=ETAp_w_m(i,j)+etam(m)*wtmij;
RHOp_w_m(i,j)=RHOp_w_m(i,j)+rhom(m)*wtmij;
Hr_w_m(i,j)=Hr_w_m(i,j)+Hrm(m)*wtmij;
Aef_w_m(i,j) = Aef_w_m(i,j) + Aefm(m)*wtmij;
K_w_m_t(i,j) = K_w_m_t(i,j) + Km(m).*wtmij;% K
DenCp_w_m_t(i,j) = DenCp_w_m_t(i,j) + DenCpm(m).*wtmij;% Den*Cp
T_w_m_t(i,j) = T_w_m_t(i,j) + Tm0(m).*DenCpm(m)*wtmij;% T
w_m_t_T(i,j) = w_m_t_T(i,j) + DenCpm(m)*wtmij;
% i+1,j Node
w_m_p_node(i+1,j)=w_m_p_node(i+1,j)+wtmi1j;
ETAp_w_m(i+1,j)=ETAp_w_m(i+1,j)+etam(m)*wtmi1j;
RHOp_w_m(i+1,j)=RHOp_w_m(i+1,j)+rhom(m)*wtmi1j;
Hr_w_m(i+1,j)=Hr_w_m(i+1,j)+Hrm(m)*wtmi1j;
Aef_w_m(i+1,j) = Aef_w_m(i+1,j) + Aefm(m)*wtmi1j;
K_w_m_t(i+1,j) = K_w_m_t(i+1,j) + Km(m).*wtmi1j;
DenCp_w_m_t(i+1,j) = DenCp_w_m_t(i+1,j) + DenCpm(m).*wtmi1j;
T_w_m_t(i+1,j) = T_w_m_t(i+1,j) + Tm0(m).*wtmi1j*DenCpm(m);% T
w_m_t_T(i+1,j) = w_m_t_T(i+1,j) + DenCpm(m)*wtmi1j;
% i,j+1 Node
w_m_p_node(i,j+1)=w_m_p_node(i,j+1)+wtmij1;
ETAp_w_m(i,j+1)=ETAp_w_m(i,j+1)+etam(m)*wtmij1;
RHOp_w_m(i,j+1)=RHOp_w_m(i,j+1)+rhom(m)*wtmij1;
Hr_w_m(i,j+1)=Hr_w_m(i,j+1)+Hrm(m)*wtmij1;
Aef_w_m(i,j+1) = Aef_w_m(i,j+1) + Aefm(m)*wtmij1;
K_w_m_t(i,j+1) = K_w_m_t(i,j+1) + Km(m).*wtmij1;
DenCp_w_m_t(i,j+1) = DenCp_w_m_t(i,j+1) + DenCpm(m).*wtmij1;
T_w_m_t(i,j+1) = T_w_m_t(i,j+1) + Tm0(m).*wtmij1*DenCpm(m);% T
w_m_t_T(i,j+1) = w_m_t_T(i,j+1) + DenCpm(m)*wtmij1;
% i+1,j+1 Node
w_m_p_node(i+1,j+1)=w_m_p_node(i+1,j+1)+wtmi1j1;
ETAp_w_m(i+1,j+1)=ETAp_w_m(i+1,j+1)+etam(m)*wtmi1j1;
RHOp_w_m(i+1,j+1)=RHOp_w_m(i+1,j+1)+rhom(m)*wtmi1j1;
Hr_w_m(i+1,j+1)=Hr_w_m(i+1,j+1)+Hrm(m)*wtmi1j1;
Aef_w_m(i+1,j+1) = Aef_w_m(i+1,j+1) + Aefm(m)*wtmi1j1;
K_w_m_t(i+1,j+1) = K_w_m_t(i+1,j+1) + Km(m).*wtmi1j1;
DenCp_w_m_t(i+1,j+1) = DenCp_w_m_t(i+1,j+1) + DenCpm(m).*wtmi1j1;
T_w_m_t(i+1,j+1) = T_w_m_t(i+1,j+1) + Tm0(m).*wtmi1j1*DenCpm(m);% T
w_m_t_T(i+1,j+1) = w_m_t_T(i+1,j+1) + DenCpm(m)*wtmi1j1;
end
for j = 1:Nx1
for i = 1:Ny1
if(w_m_p_node(i,j)>0)
ETAP(i,j) = ETAp_w_m(i,j)/w_m_p_node(i,j);
RHOP(i,j) = RHOp_w_m(i,j)/w_m_p_node(i,j);
Hrij(i,j) = Hr_w_m(i,j)/w_m_p_node(i,j);
Aefij(i,j) = Aef_w_m(i,j)/w_m_p_node(i,j);
kht(i,j) = K_w_m_t(i,j)/w_m_p_node(i,j);%K
hct(i,j) = DenCp_w_m_t(i,j)/w_m_p_node(i,j);%Den*Cp
% else
% kht(i,j) = kht_backup(i,j);
% hct(i,j) = hct_backup(i,j);
end
td(i,j) = kht(i,j)/hct(i,j);%K/(Den*Cp)
if(w_m_t_T(i,j)>0)
T0(i,j) = T_w_m_t(i,j)/w_m_t_T(i,j);%T
% else
% T0(i,j) = T_backup(i,j);
end
end
end
% Apply boundary condition after interpolate T0 from markers(any t).
% top and left Ibc
T0(1,2:Nx) = T0(2,2:Nx);
T0(Ny1,2:Nx) = T0(Ny,2:Nx);
% left and right Ibc
T0(:,1) = T0(:,2);
T0(:,Nx1) = T0(:,Nx);
% Solve the possion equation of gravitational potential(Φ)
% coefficients of boundary
% up and down
K = 2/3;% circular
G = 6.672e-11;% (Nm^2)/(kg^2),gravitational constant.
for j = 1:Nx1
for i = 1:Ny1
k = (j-1)*Ny1 + i;
distance = ((xp(j)-xmid)^2 + (yp(i)-ymid)^2)^0.5;
if(distance > Lx/2 || i==1 || i==Ny1 || j==1 || j==Nx1)
% Graviational = 0
coeg(k, k) = 1;
bg(k,1) = 0;
else
% Distribution of nodes.
% j-1 j j+1
% i-1 2Φ(i-1,j)
% i 1Φ(i,j-1) 3Φ(i,j) 5Φ(i,j+1)
% i+1 4Φ(i+1,j)
coeg(k, k) = -2/(dx.^2)-2/(dy.^2); % middleΦ(i,j) 3
coeg(k, k-1) = 1/(dy.^2); % up Φ(i-1,j) 2
coeg(k, k+1) = 1/(dy.^2); % down Φ(i+1,j) 4
coeg(k, k-Ny1) = 1/(dx.^2); % left Φ(i,j-1) 1
coeg(k, k+Ny1) = 1/(dx.^2); % right Φ(i,j+1) 5
bg(k,1) = 4*K*pi*G*RHOP(i,j);
end
end
end
ug = coeg\bg;
Ug = reshape(ug,Ny1,Nx1);
for j = 1:Nx
for i = 1:Ny
g_x(i,j) = -(Ug(i,j+1)-Ug(i,j))/dx;
g_y(i,j) = -(Ug(i+1,j)-Ug(i,j))/dy;
end
end
% Mechanical Solution
% Solve Stokes-equ and Continuity-equ.
for j=1:1:Nx1
for i=1:1:Ny1
% Define global indexes in algebraic space
kvx=((j-1)*Ny1+i-1)*3+1; % Vx
kvy=kvx+1; % Vy
kpm=kvx+2; % P
% Vx equation External points
if(i==1 || i==Ny1 || j==1 || j==Nx || j==Nx1)
% Boundary Condition
% 1*Vx=0
coe(kvx,kvx)=1; % Left part
R(kvx)=0; % Right part
% Top boundary
if(i==1 && j>1 && j<Nx)
coe(kvx,kvx+3)=bc_top; % Left part
end
% Bottom boundary
if(i==Ny1 && j>1 && j<Nx)
coe(kvx,kvx-3)=bc_bottom; % Left part
end
else
% Internal points: x-Stokes eq.
% ETA*(d2Vx/dx^2+d2Vx/dy^2)-dP/dx=0
% Vx2
% |
% Vy1 | Vy3
% |
% Vx1-P1-Vx3-P2-Vx5
% |
% Vy2 | Vy4
% |
% Vx4
%
% Viscosity points
ETA1=ETA(i-1,j);
ETA2=ETA(i,j);
ETAP1=ETAP(i,j);
ETAP2=ETAP(i,j+1);
% Density gradients
dRHOdx=(RHOVX(i,j+1)-RHOVX(i,j-1))/2/dx;
dRHOdy=(RHOVX(i+1,j)-RHOVX(i-1,j))/2/dy;
% Left part
coe(kvx,kvx-Ny1*3)=2*ETAP1/dx^2; % Vx1
coe(kvx,kvx-3)=ETA1/dy^2; % Vx2
coe(kvx,kvx)=-2*(ETAP1+ETAP2)/dx^2-(ETA1+ETA2)/dy^2-dRHOdx*g_x(i,j)*dt; % Vx3
coe(kvx,kvx+3)=ETA2/dy^2; % Vx4
coe(kvx,kvx+Ny1*3)=2*ETAP2/dx^2; % Vx5
coe(kvx,kvy)=-ETA2/dx/dy-dRHOdy*g_x(i,j)*dt/4; % Vy2
coe(kvx,kvy+Ny1*3)=ETA2/dx/dy-dRHOdy*g_x(i,j)*dt/4; % Vy4
coe(kvx,kvy-3)=ETA1/dx/dy-dRHOdy*g_x(i,j)*dt/4; % Vy1
coe(kvx,kvy+Ny1*3-3)=-ETA1/dx/dy-dRHOdy*g_x(i,j)*dt/4; % Vy3
coe(kvx,kpm)=Kcont/dx; % P1
coe(kvx,kpm+Ny1*3)=-Kcont/dx; % P2
% Right part !!!
R(kvx)=-RHOVX(i,j)*g_x(i,j);
end
% Vy equation External points
if(j==1 || j==Nx1 || i==1 || i==Ny || i==Ny1)
% Boundary Condition
% 1*Vy=0
coe(kvy,kvy)=1; % Left part
R(kvy)=0; % Right part
% Left boundary
if(j==1 && i>1 && i<Ny)
coe(kvy,kvy+3*Ny1)=bc_left; % Left part
end
% Right boundary
if(j==Nx1 && i>1 && i<Ny)
coe(kvy,kvy-3*Ny1)=bc_right; % Left part
end
else
% Internal points: y-Stokes eq.
% ETA*(d2Vy/dx^2+d2Vy/dy^2)-dP/dy=-RHO*gy
% Vy2
% |
% Vx1 P1 Vx3
% |
% Vy1----Vy3----Vy5
% |
% Vx2 P2 Vx4
% |
% Vy4
%
% Viscosity points
% Viscosity points
ETA1=ETA(i,j-1);
ETA2=ETA(i,j);
ETAP1=ETAP(i,j);
ETAP2=ETAP(i+1,j);
dRHOdx=(RHOVY(i,j+1)-RHOVY(i,j-1))/2/dx;
dRHOdy=(RHOVY(i+1,j)-RHOVY(i-1,j))/2/dy;
% Left part
coe(kvy,kvy-Ny1*3)=ETA1/dx^2; % Vy1
coe(kvy,kvy-3)=2*ETAP1/dy^2; % Vy2
coe(kvy,kvy)=-2*(ETAP1+ETAP2)/dy^2-...
(ETA1+ETA2)/dx^2 - dRHOdy*g_y(i,j)*dt; % Vy3
coe(kvy,kvy+3)=2*ETAP2/dy^2; % Vy4
coe(kvy,kvy+Ny1*3)=ETA2/dx^2; % Vy5
coe(kvy,kvx)=-ETA2/dx/dy -dRHOdx*g_y(i,j)*dt/4; %Vx3
coe(kvy,kvx+3)=ETA2/dx/dy -dRHOdx*g_y(i,j)*dt/4; %Vx4
coe(kvy,kvx-Ny1*3)=ETA1/dx/dy -dRHOdx*g_y(i,j)*dt/4; %Vx1
coe(kvy,kvx+3-Ny1*3)=-ETA1/dx/dy -dRHOdx*g_y(i,j)*dt/4; %Vx2
coe(kvy,kpm)=Kcont/dy; % P1
coe(kvy,kpm+3)=-Kcont/dy; % P2
% Right part
R(kvy)=-RHOVY(i,j)*g_y(i,j);
end
% P equation External points
if(i==1 || j==1 || i==Ny1 || j==Nx1 ||...
(i==2 && j==2))
% Boundary Condition
% 1*P=0
coe(kpm,kpm)=1; % Left part
R(kpm)=0; % Right part
% Real BC
if(i==2 && j==2)
coe(kpm,kpm)=1*Kcont; %Left part
R(kpm)=1e+9; % Right part
end
else
% Internal points: continuity eq.
% dVx/dx+dVy/dy=0
% Vy1
% |
% Vx1--P--Vx2
% |
% Vy2
%
% Left part
coe(kpm,kvx-Ny1*3)=-1/dx; % Vx1
coe(kpm,kvx)=1/dx; % Vx2
coe(kpm,kvy-3)=-1/dy; % Vy1
coe(kpm,kvy)=1/dy; % Vy2
% Right part
R(kpm)=0;
end
end
end
% Computing the solution U:include vx,vy and p
U = coe \ R;
% Decompose the solution of global matrix.
for j=1:1:Nx1
for i=1:1:Ny1
% Define global indexes in algebraic space
kvx=((j-1)*Ny1+i-1)*3+1; % Vx
kvy=kvx+1; % Vy
kpm=kvx+2; % P
% Reload solution
vx_matrix(i,j)=U(kvx);
vy_matrix(i,j)=U(kvy);
p_matrix(i,j)=U(kpm)*Kcont;
end
end
% Computing strain rate(xy) and deviatoric(xy) stress components.
for j = 1:Nx
for i = 1:Ny
srxy(i,j) = (1/2)*( (vx_matrix(i+1,j)-vx_matrix(i,j))/dy+(vy_matrix(i,j+1)-vy_matrix(i,j))/dx ) ;
dsxy(i,j) = 2*ETA(i,j)*srxy(i,j);
end
end
% Computing strain rate(xx) and deviatoric(xx) stress components.
for j = 2:Nx
for i = 2:Ny
srxx(i,j) = ( vx_matrix(i,j)-vx_matrix(i,j-1) )/dx;
dsxx(i,j) = 2*ETAP(i,j)*srxx(i,j);
end
end
for j = 2:Nx
for i = 2:Ny
if(i~=1 && j~=1)
kxyxyij = srxy(i,j)*dsxy(i,j);
kxyxyi1j = srxy(i-1,j)*dsxy(i-1,j);
kxyxyij1 = srxy(i,j-1)*dsxy(i,j-1);
kxyxyi1j1 = srxy(i-1,j-1)*dsxy(i-1,j-1);
pxyxy(i,j) = (1/4)*(kxyxyij+kxyxyi1j+kxyxyij1+kxyxyi1j1);
kxxxxij = srxx(i,j)*dsxx(i,j);
kxxxxi1j = srxx(i-1,j)*dsxx(i-1,j);
kxxxxij1 = srxx(i,j-1)*dsxx(i,j-1);
kxxxxi1j1 = srxx(i-1,j-1)*dsxx(i-1,j-1);
pxxxx(i,j) = (1/4)*(kxxxxij+kxxxxi1j+kxxxxij1+kxxxxi1j1);
end
% Shear heating Hs
Hsij(i,j) = 2*( pxyxy(i,j) + pxxxx(i,j) ) ;
end
end
% Adiabatic heating Ha
for j = 2:Nx1
for i = 2:Ny1
kvy = ( vy_matrix(i,j)+vy_matrix(i-1,j) )/2;
kvx = ( vx_matrix(i,j)+vx_matrix(i,j-1) )/2;
GXP=(g_x(i,j)+g_x(i,j-1))/2;
GYP=(g_y(i,j)+g_y(i-1,j))/2;
Haij(i,j) = T0(i,j)*Aefij(i,j)*RHOP(i,j)*(GYP*kvy+GXP*kvx);
end
end
% Computing average velocity components at cell centres.
for j = 2:Nx
for i = 2:Ny
% Computing internal Vx and Vy at the P nodes.
% Vx at internal(P) nodes.
pvx(i,j) = ( vx_matrix(i,j) + vx_matrix(i,j-1) )./2;
% Vy at internal nodes.
pvy(i,j) = ( vy_matrix(i,j) + vy_matrix(i-1,j) )./2;
end
end
% Why it is free-slip ??
% Applying free-slip boundary condition for velocity in the P nodes .
% j=1 j=2 j=3 ... j=Nx_node+1
% i=1
% i=2
% i=3
% ...
% i=Ny_node+1
% Apply BC
% Top
pvx(1,2:Nx-1)=-bc_top*pvx(2,2:Nx-1);
% Bottom
pvx(Ny1,2:Nx-1)=-bc_bottom*pvx(Ny,2:Nx-1);
% Left
pvx(:,1)=-pvx(:,2);
% Right
pvx(:,Nx1)=-pvx(:,Nx);
% Apply BC
% Left
pvy(2:Ny-1,1)=-bc_left*pvy(2:Ny-1,2);
% Right
pvy(2:Ny-1,Nx1)=-bc_right*pvy(2:Ny-1,Nx); % Free slip
% Top
pvy(1,:)=-pvy(2,:);
% Bottom
pvy(Ny1,:)=-pvy(Ny,:);
% Restrict the time step limiting max_marker_disp by 0.5 grid step
% Define time steps.
dt = dt*dtkoef;
maxvx=max(max(abs(vx_matrix)));
maxvy=max(max(abs(vy_matrix)));
if(dt*maxvx>dxymax*dx)
dt=dxymax*dx/maxvx;
end
if(dt*maxvy>dxymax*dy)
dt=dxymax*dy/maxvy;
end
khtdt = kht;
hctdt = hct;
for titer = 1:1:2
% Ctbc
% Why???
% Up and down
for j = 2:Nx
i = 1;
k = (j-1)*Ny1 + i;
% T(1,j) - T(2,j) = 0;
coet(k, k) = 1;
coet(k, k+1) = -1;
bt(k, 1) = 0;
i = Ny1;
k = (j-1)*Ny1 + i;
% T(Ny,j) - T(Ny-1,j) = 0;
coet(k, k) = 1;
coet(k, k-1) = -1;
bt(k, 1) = 0;
end
% left and right
for i = 1:Ny1
j = 1;
k = (j-1)*Ny1 + i;
coet(k, k) = 1;
coet(k, k+Ny1) = -1;
bt(k, 1) = 0;
j = Nx1;
k = (j-1)*Ny1 + i;
coet(k, k) = 1;
coet(k, k-Ny1) = -1;
bt(k, 1) = 0;
end
% coefficients of interval points
for j = 2:Nx
for i = 2:Ny
k = (j-1)*Ny1 + i; % k is the index
if vx_matrix(i,j)>0
dtdx = (T0(i,j)-T0(i,j-1))/dx;
elseif vx_matrix(i,j)<0
dtdx = (T0(i,j+1)-T0(i,j))/dx;
end
if vy_matrix(i,j)>0
dtdy = (T0(i,j)-T0(i-1,j))/dy;
elseif vy_matrix(i,j)<0
dtdy = (T0(i+1,j)-T0(i,j))/dy;
end
% KVX == Kh;KVY == Kv
khij = KVX(i,j);
khij1 = KVX(i,j-1);
kvij = KVY(i,j);
kvi1j = KVY(i-1,j);
% Not 2* why?? page177
coet(k, k) = (khij+khij1)/dx^2+hctdt(i,j)/dt...
+(kvij+kvi1j)/dy^2;% middle T3 i,j
coet(k, k-1) = -(kvi1j)/(2.*dy.^2); % up T2 i-1,j
coet(k, k+1) = -(kvij)/(2.*dy.^2); % down T4 i+1,j
coet(k, k-Ny1) = -(khij1)/(dx.^2); % left T1 i,j-1
coet(k, k+Ny1) = -(khij)/(dx.^2); % right T5 i,j+1
% !!!Note that there are not advection term because MIC have
% marker transport.!!!
% Adding H(P-nodes)
bt(k, 1) = T0(i,j)*hctdt(i,j)/dt + Hrij(i,j) + Haij(i,j) + Hsij(i,j); %- hctdt(i,j)*(vx(i,j)*dtdx+vy(i,j)*dtdy);
end
end
%direct method:right martix is divided by coe martix on the left
u = coet \ bt;
%transform vector to grids
% U has Ny1 rows and Nx1 columns
U = reshape(u, Ny1, Nx1);
T1 = U;
% Computing dTi,j = dTi,j(subgrid) + dTi,j(remaining)
% Eq10.13
dTij = T1-T0;
maxDTcurrent = max(max(abs(dTij)));
if(titer<2 && maxDTcurrent>DTmax)
dt=dt/maxDTcurrent*DTmax;
else
break;
end
end
% Computing subgrid diffusion for markers
if (subgrdifcoe>0)
% Existence of subgrid diffusion.
% Page162 How to correct the problem ?
% operating consistent subgrid diffusion .
% Clear subgrid temperature changes(dTijsub) for nodes
% From marker to P-NODES
w_Pnodes = zeros(Ny1,Nx1);
DWTM = zeros(Ny1,Nx1);
% Interpolate Tm0(nodal) from E-nodes(T0) to L-markers.
for m = 1:marknum
% Interpolate T to P-NODES(Nx*Ny)
% Define i,j indexes for the upper left node
j=fix((xm(m)-xp(1))/dx)+1;
i=fix((ym(m)-yp(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx)
j=Nx;
end
if(i<1)
i=1;
elseif(i>Ny)
i=Ny;
end
% Compute distances
dis_x=xm(m)-xp(j);
dis_y=ym(m)-yp(i);
% Compute weights
wtmij=(1-dis_x/dx)*(1-dis_y/dy);
wtmi1j=(1-dis_x/dx)*(dis_y/dy);
wtmij1=(dis_x/dx)*(1-dis_y/dy);
wtmi1j1=(dis_x/dx)*(dis_y/dy);
% Interpolate Basic node(Tm0nodal) from L-markers(T0i,j)
Tm0nodal= T0(i,j)*wtmij + T0(i+1,j)*wtmi1j ...
+T0(i,j+1)*wtmij1 + T0(i+1,j+1)*wtmi1j1;
% Eq.10.14
dTm = Tm0nodal - Tm0(m);
% Page163
tdiff = DenCpm(m)/Km(m)/(2/dx^2+2/dy^2);
% Computing subgrid diffusion
sdif=-subgrdifcoe*dt/tdiff;
if(sdif<-30)
sdif=-30;
end
% Page162 Eq.10.16
dTmsub=dTm*(1-exp(sdif));
% Eq.10.19
% Tcorrected = Tm0 + dTsubgrid (+ dTremain)
Tm0(m) = Tm0(m) + dTmsub;
% Interpolate Basic node(dTijsub) from L-markers(dTmsub)
% Eq.10.17
DWTM(i,j) = DWTM(i,j) + dTmsub*DenCpm(m)*wtmij;
w_Pnodes(i,j) = w_Pnodes(i,j) + DenCpm(m)*wtmij;
DWTM(i+1,j) = DWTM(i+1,j) + dTmsub*DenCpm(m)*wtmi1j;
w_Pnodes(i+1,j) = w_Pnodes(i+1,j) + DenCpm(m)*wtmi1j;
DWTM(i,j+1) = DWTM(i,j+1) + dTmsub*DenCpm(m)*wtmij1;
w_Pnodes(i,j+1) = w_Pnodes(i,j+1) + DenCpm(m)*wtmij1;
DWTM(i+1,j+1) = DWTM(i+1,j+1) + dTmsub*DenCpm(m)*wtmi1j1;
w_Pnodes(i+1,j+1) = w_Pnodes(i+1,j+1) + DenCpm(m)*wtmi1j1;
end
dTijsub = zeros(Ny1,Nx1);
for j = 1:Nx1
for i = 1:Ny1
if(w_Pnodes(i,j)>0)
% Compting dTijsub
% Eq.10.17
dTijsub(i,j) = DWTM(i,j)/w_Pnodes(i,j);
end
end
end
% Eq.10.18
% Update dTij
dTijremain = dTij - dTijsub;
end
% Interpolate L-markers(dTremain) from P-node(dTijremain)
for m = 1:marknum
% Define i,j indexes for the upper left node
j=fix((xm(m)-xp(1))/dx)+1;
i=fix((ym(m)-yp(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx)
j=Nx;
end
if(i<1)
i=1;
elseif(i>Ny)
i=Ny;
end
% Compute distances
dis_x=xm(m)-xp(j);
dis_y=ym(m)-yp(i);
% Compute weights
wtmij=(1-dis_x/dx)*(1-dis_y/dy);
wtmi1j=(1-dis_x/dx)*(dis_y/dy);
wtmij1=(dis_x/dx)*(1-dis_y/dy);
wtmi1j1=(dis_x/dx)*(dis_y/dy);
% Interpolate Basic node(Tm0nodal) from L-markers(T0i,j)
dTmremain=dTijremain(i,j)*wtmij + dTijremain(i+1,j)*wtmi1j ...
+dTijremain(i,j+1)*wtmij1 + dTijremain(i+1,j+1)*wtmi1j1;
% Tcorrected = Tm0(Tm0+ dTsubgrid) + dTremain
% Eq.10.19
if(t>1)
Tm0(m) = Tm0(m) + dTmremain;
else
% Interpolate new temperature for the 1st timestep
Tm0(m) = T1(i,j)*wtmij + T1(i+1,j)*wtmi1j ...
+T1(i,j+1)*wtmij1 + T1(i+1,j+1)*wtmi1j1;
end
end
% Figure
subplot(1,3,1);
colormap('Jet');
pcolor(x,y,log10(ETA));
hold on;
%P nodes:(Nx1 * Ny1).
quiver(xp(1:2:Nx1),yp(1:2:Ny1),pvx(1:2:Ny1,1:2:Nx1),pvy(1:2:Ny1,1:2:Nx1),'w')
xlabel('Horizontal(m)')
ylabel('Vertical(m)')
caxis([17 21]);
set(gca,'xaxislocation','top');
set (gca,'YDir','reverse')
shading interp;
axis ij image;
colorbar
title(['log10(ETA) after ',num2str(t),' deltat'])
% Ha
subplot(1,3,2);
colormap('Jet');
pcolor(xp,yp,Hrij);
xlabel('Horizontal(m)')
ylabel('Vertical(m)')
set(gca,'xaxislocation','top');
set (gca,'YDir','reverse')
shading interp;
colorbar;
axis ij image;
title(['Ha after ',num2str(t),' deltat'])
subplot(1,3,3);
pcolor(xp,yp,T1);colorbar;
xlabel('Horizontal(Km)');
ylabel('Vertical(Km)');
zlabel('Gravitational potential');
title(['(Ibc)(MIC) Temperature after',num2str(t),' deltat'])
set(gca,'xaxislocation','top');
set (gca,'YDir','reverse')
axis ij image;
shading interp;
pause(0.01);
% The classical Runge-Kutta scheme,Fourth-order in sapce First-order in time
% Define Vxa(Vya)1,Vxb(Vyb)2,Vxc(Vyc)3,Vxd(Vyd)4.
vxm=zeros(4,1);
vym=zeros(4,1);
% Interpolate vx and vy components from E-nodes to L-markers.
for m=1:1:marknum
% Reserve original coordinate.
xa = xm(m);
ya = ym(m);
for rk = 1:1:4
% Interpolate vx-p and vy-p from P-staggered nodes to L-markers.
% P staggered nodes.(Nx1 * Ny1)
j=fix((xm(m)-xp(1))/dx)+1;
i=fix((ym(m)-yp(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx)
j=Nx;
end
if(i<1)
i=1;
elseif(i>Ny)
i=Ny;
end
% Compute distances
dis_x=xm(m)-xp(j);
dis_y=ym(m)-yp(i);
% Compute weights
wtmij=(1-dis_x/dx)*(1-dis_y/dy);
wtmi1j=(dis_x/dx)*(1-dis_y/dy);
wtmij1=(1-dis_x/dx)*(dis_y/dy);
wtmi1j1=(dis_x/dx)*(dis_y/dy);
% Compute vxm and vym velocity
vxm(rk) = pvx(i,j)*wtmij+pvx(i+1,j)*wtmi1j...
+pvx(i,j+1)*wtmij1+pvx(i+1,j+1)*wtmi1j1;
vym(rk) = pvy(i,j)*wtmij+pvy(i+1,j)*wtmi1j...
+pvy(i,j+1)*wtmij1+pvy(i+1,j+1)*wtmi1j1;
% Interpolate vx from Vx-staggered nodes to L-markers.
% Define i,j indexes for the upper left node
j=fix((xm(m)-xvx(1))/dx)+1;
i=fix((ym(m)-yvx(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx-1)
j=Nx-1;
end
% Vx staggered nodes.(Nx * Ny+1)
if(i<1)
i=1;
elseif(i>Ny)
i=Ny;
end
% Compute distances
dis_x=xm(m)-xvx(j);
dis_y=ym(m)-yvx(i);
% Compute weights
wtmij=(1-dis_x/dx)*(1-dis_y/dy);
wtmi1j=(dis_x/dx)*(1-dis_y/dy);
wtmij1=(1-dis_x/dx)*(dis_y/dy);
wtmi1j1=(dis_x/dx)*(dis_y/dy);
% Compute vx velocity
vxcij = (1/3)*(2*vx_matrix(i,j));
vxci1j = (1/3)*(2*vx_matrix(i+1,j));
vxcij1 = (1/3)*(2*vx_matrix(i,j+1));
vxci1j1 = (1/3)*(2*vx_matrix(i+1,j+1));
vxm(rk) = (1/3)*vxm(rk) + (vxcij*wtmij+vxci1j*wtmi1j...
+vxcij1*wtmij1+vxci1j1*wtmi1j1);
% Interpolate vy
% Define i,j indexes for the upper left node
j=fix((xm(m)-xvy(1))/dx)+1;
i=fix((ym(m)-yvy(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx)
j=Nx;
end
if(i<1)
i=1;
elseif(i>Ny-1)
i=Ny-1;
end
% Compute distances
dis_x=xm(m)-xvy(j);
dis_y=ym(m)-yvy(i);
% Compute weights
wtmij=(1-dis_x/dx)*(1-dis_y/dy);
wtmi1j=(1-dis_x/dx)*(dis_y/dy);
wtmij1=(dis_x/dx)*(1-dis_y/dy);
wtmi1j1=(dis_x/dx)*(dis_y/dy);
% Compute vx velocity
vycij = (1/3)*(2*vy_matrix(i,j));
vyci1j = (1/3)*(2*vy_matrix(i+1,j));
vycij1 = (1/3)*(2*vy_matrix(i,j+1));
vyci1j1 = (1/3)*(2*vy_matrix(i+1,j+1));
vym(rk) = (1/3)*vym(rk) + (vycij*wtmij+vyci1j*wtmi1j...
+vycij1*wtmij1+vyci1j1*wtmi1j1);
if(rk == 1 || rk == 2)
xm(m) = vx_matrix(rk)*dt/2 + xa;
ym(m) = vy_matrix(rk)*dt/2 + ya;
elseif(rk == 3)
xm(m) = vx_matrix(rk)*dt + xa;
ym(m) = vy_matrix(rk)*dt + ya;
end
end
% Return original coordinate.
xm(m) = xa;
ym(m) = ya;
% Move markers
vxefft = (1/6)*(vxm(1)+2*vxm(2)+2*vxm(3)+vxm(4));
vyefft = (1/6)*(vym(1)+2*vym(2)+2*vym(3)+vym(4));
xm(m)=xm(m)+dt*vxefft;
ym(m)=ym(m)+dt*vyefft;
end
end
Initial model：
Model after 61 dt：
文章来源：https://uudwc.com/A/LmGWJ
文章来源地址https://uudwc.com/A/LmGWJ
