这是一个暂时保存的代码
%This program complished the process of 2D Eulerian advection with method
%of MIC(allocating memory temporarily).
clear all;clf;
% Set parameters
Nt = 1;
output_interval = 1;
% Length of model
Lx = 500000;
Ly = 500000;
x_mid = Lx/2;
y_mid = Ly/2;
% Sun of point number of model
Nx = 51;
Ny = 51;
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 = zeros(Ny1,Nx1);
vy = zeros(Ny1,Nx1);
p = zeros(Ny1,Nx1);
sxx = zeros(Ny1,Nx1);
syy = zeros(Ny1,Nx1);
exx = zeros(Ny1,Nx1);
eyy = zeros(Ny1,Nx1);
% Vx at Pnodes
pvx = zeros(Ny1,Nx1);
pvy = zeros(Ny1,Nx1);
% Gravity acceleration
g_y = 10;
g_x = 0;
% Basic node
ETA = zeros(Ny,Nx);%density
ETAVP = zeros(Ny,Nx);
MIU = zeros(Ny,Nx);
RHO = zeros(Ny,Nx);%viscosity
sxy = zeros(Ny,Nx);
syx = zeros(Ny,Nx);
% Stokes-equ and continuity-equ
% Temperature-equ
% Amount of markers.
Nx_marker = 200;
Ny_marker = 300;
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); % Viscosity
mium = zeros(1,marknum); % shear mudulus
syym = zeros(1,marknum);
sxxm = zeros(1,marknum);
sxym = 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 = sparse(Nx1*Ny1, Nx1*Ny1);% Setup of model
% Mantle
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;
miu_mantle = 5e10;
% Plume
den_plume = 3200;
vis_plume = 1e+20;
k_plume = 2;%W/(m*K)
Cp_plume = 1100;%J/(Kg*K)
t_plume = 1800;%K
aef_plume = 2e-5;%1/K
Hr_plume = 3e-8;
miu_plume = 3e10;% Striky air
den_striky = 1;
vis_striky = 1e+17;
k_striky = 300;%W/(m*K)
Cp_striky = 1100;%J/(Kg*K)
t_striky = 273;%K
aef_striky = 0;%1/K
Hr_striky = 0;
radius = 100000;
miu_striky = 2e10;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^3ETAP = zeros(Ny1,Nx1);
ETAPVP = zeros(Ny1,Nx1);
RHOP = zeros(Ny1,Nx1);
MIUP = zeros(Ny1,Nx1);
% Parameters and matrixs of exercise 9.3
srxy = zeros(Ny,Nx);
dsxy = zeros(Ny,Nx);
srxx = zeros(Ny1,Nx1);
dsxx = zeros(Ny1,Nx1);
pxyxy = zeros(Ny1,Nx1);
pxxxx = zeros(Ny1,Nx1);
Hsij = zeros(Ny1,Nx1);
T = 1300;%K
aef93 = 3e-5;
Haij = zeros(Ny1,Nx1);Hrij = zeros(Ny1,Nx1);
Aefij = zeros(Ny1,Nx1);
% % DEFINE dt_tem
% kappa=max(max(Km))/min(min(DenCpm));
% dtexp = min(dx,dy)^2/(3*kappa);
% if(vx_matrix(1,1)~=0)
% dtexp=min(dtexp,abs(dx/max(max(vx_matrix)))); % Limitation for horizontal advection timestep
% end
% if(vy_matrix(1,1)~=0)
% dtexp=min(dtexp,abs(dy/max(max(vy_matrix)))); % Limitation for vertical advection timestep
% end
% dt_tem = 0.5*dtexp;
% subgrdifcoe = 1;m = 1;
for jm=1:1:Nx_markerfor im=1:1:Ny_marker% Define marker coordinatesxm(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 propertiesrmark=((xm(m)-x_mid)^2+(ym(m)-y_mid)^2)^0.5;if(rmark<=radius)rhom(m) = den_plume; % Plume densityetam(m) = vis_plume; % Plume viscosityKm(m) = k_plume; % Plume K DenCpm(m) = 3.2e+6; % Plume den*CpTm0(m) = t_plume; %Plume T Aefm(m) = aef_plume;Hrm(m) = Hr_plume;mium(m) = miu_plume;elseif(ym(m)<=0.2*Ly)rhom(m) = den_striky; % Sticky densityetam(m) = vis_striky; % Sticky viscosityKm(m) = k_striky; % Sticky K DenCpm(m) = 3.3e+6; % Sticky den*CpTm0(m) = t_striky; % Sticky TAefm(m) = aef_striky;Hrm(m) = Hr_striky;mium(m) = miu_striky;elserhom(m) = den_mantle; % Mantle densityetam(m) = vis_mantle; % Mantle viscosityKm(m) = k_mantle; % Plume K DenCpm(m) = 3.3e+6; % Plume den*CpTm0(m) = t_mantle; %Plume T Aefm(m) = aef_mantle;Hrm(m) = Hr_mantle;mium(m) = miu_mantle;endsxxm(m) = 10^6;sxym(m) = 0;% Update marker counterm=m+1;end
end
Kcont = 1e+21/dx;% Unknown_number is (Ny+1)*(Nx+1)*3
unknowns = (Nx+1)*(Ny+1)*3;
coe = sparse(unknowns,unknowns);
R = zeros(unknowns, 1);% Boundary condition : free slip = -1 ; no slip = 1
bc_left = -1;
bc_right = -1;
bc_top = -1;
bc_bottom = -1;dxymax = 0.5; % Define dt for stokes equation.
dt = 1e+10; % initial timestepdtkoef=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 NODEw_m_x_node = zeros(Ny,Nx);eta_w_m = zeros(Ny,Nx);miu_w_m = zeros(Ny,Nx);sxy_w_m = zeros(Ny,Nx);% rho in vy and vx nodesRHOYSUM=zeros(Ny1,Nx1);WTYSUM=zeros(Ny1,Nx1);RHOXSUM=zeros(Ny1,Nx1);WTXSUM=zeros(Ny1,Nx1);% K in vy and vx nodesKYSUM=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 markesj=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;endif(i<1)i=1;elseif(i>Ny-1)i=Ny-1;end% Compute distancesdis_x = xm(m)-x(j);dis_y = ym(m)-y(i);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);% | |% * * * * * Discard the m-point on the right.w_m_x_node(i,j) = w_m_x_node(i,j) + wtmij;eta_w_m(i,j) = eta_w_m(i,j) + etam(m).*wtmij;miu_w_m(i,j) = miu_w_m(i,j) + wtmij/mium(m);sxy_w_m(i,j) = sxy_w_m(i,j) + sxym(m)*wtmij;% i+1,jw_m_x_node(i+1,j) = w_m_x_node(i+1,j) + wtmi1j;eta_w_m(i+1,j) = eta_w_m(i+1,j) + etam(m).*wtmi1j;miu_w_m(i+1,j) = miu_w_m(i+1,j) + wtmi1j/mium(m);sxy_w_m(i+1,j) = sxy_w_m(i+1,j) + sxym(m)*wtmi1j;% i,j+1w_m_x_node(i,j+1) = w_m_x_node(i,j+1) + wtmij1;eta_w_m(i,j+1) = eta_w_m(i,j+1) + etam(m).*wtmij1;miu_w_m(i,j+1) = miu_w_m(i,j+1) + wtmij1/mium(m);sxy_w_m(i,j+1) = sxy_w_m(i,j+1) + sxym(m)*wtmij1;% i+1,j+1w_m_x_node(i+1,j+1) = w_m_x_node(i+1,j+1) + wtmi1j1;eta_w_m(i+1,j+1) = eta_w_m(i+1,j+1) + etam(m).*wtmi1j1;miu_w_m(i+1,j+1) = miu_w_m(i+1,j+1) + wtmi1j1/mium(m);sxy_w_m(i+1,j+1) = sxy_w_m(i+1,j+1) + sxym(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;endif(i<1)i=1;elseif(i>Ny-1)i=Ny-1;end% Compute distancesdis_x=xm(m)-xvy(j);dis_y=ym(m)-yvy(i);% Compute weightswtmij=(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);% Update properties% i,j Node% rho and k in vy nodesRHOYSUM(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 NodeRHOYSUM(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 NodeRHOYSUM(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 NodeRHOYSUM(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;endif(i<1)i=1;elseif(i>Ny)i=Ny;end% Compute distancesdis_x=xm(m)-xvx(j);dis_y=ym(m)-yvx(i);% Compute weightswtmij=(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);% Update properties% i,j NodeRHOXSUM(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 NodeRHOXSUM(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 NodeRHOXSUM(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 NodeRHOXSUM(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:Nxfor i = 1:Nyif(w_m_x_node(i,j)>0)% Not need RHO on the basic nodes.Just Y nodes.ETA(i,j) = eta_w_m(i,j)/w_m_x_node(i,j);ETAVP(i,j) = ETA(i,j);%etaVp = eta in viscosity-elastic matreials.MIU(i,j) = w_m_x_node(i,j)/miu_w_m(i,j);sxy(i,j) = sxy_w_m(i,j)/w_m_x_node(i,j);syx(i,j) = sxy(i,j);endendend% Computing Density and K value at the VY and VX-nodes.for j=1:1:Nx1for i=1:1:Ny1% VY-nodes:(Nx1 * Ny) absent Ny1=0 !!!if(WTYSUM(i,j)>0)% rou_vRHOVY(i,j)=RHOYSUM(i,j)/WTYSUM(i,j);% K_vKVY(i,j)=KYSUM(i,j)/WTYSUM(i,j);end% VX-nodes:(Nx * Ny1)absent Nx1=0 !!!if(WTXSUM(i,j)>0)% rou_hRHOVX(i,j)=RHOXSUM(i,j)/WTXSUM(i,j);% K_hKVX(i,j) = KXSUM(i,j)/WTXSUM(i,j);endendend% Interpolate ETA from L-markers to P-nodes.w_m_p_node = zeros(Ny1,Nx1);% weights except TETAp_w_m = zeros(Ny1,Nx1);% Viscosity at P nodesxxp_w_m = zeros(Ny1,Nx1);syyp_w_m = zeros(Ny1,Nx1);MIUp_w_m = zeros(Ny1,Nx1);% miu at P nodeRHOp_w_m = zeros(Ny1,Nx1);% Density at P nodeHr_w_m = zeros(Ny1,Nx1); % Hr at P nodeAef_w_m = zeros(Ny1,Nx1); % Aef at P nodeK_w_m_t = zeros(Ny1,Nx1); % K at P nodeDenCp_w_m_t = zeros(Ny1,Nx1);% Den*Cp at P nodeT_w_m_t = zeros(Ny1,Nx1); % T at P nodew_m_t_T = zeros(Ny1,Nx1); % weight of T% BACKUP of K,HC and T.kht_backup = kht;hct_backup = hct;T_backup = T0;for m=1:1:marknum% Reserve original coordinate.xa = xm(m);ya = ym(m);% 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 distancesdis_x=xm(m)-xp(j);dis_y=ym(m)-yp(i);% Compute weightswtmij=(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 velocityw_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;sxxp_w_m(i,j)=sxxp_w_m(i,j)+sxxm(m)*wtmij;MIUp_w_m(i,j)=MIUp_w_m(i,j)+wtmij/mium(m);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;% KDenCp_w_m_t(i,j) = DenCp_w_m_t(i,j) + DenCpm(m).*wtmij;% Den*CpT_w_m_t(i,j) = T_w_m_t(i,j) + Tm0(m).*DenCpm(m)*wtmij;% Tw_m_t_T(i,j) = w_m_t_T(i,j) + DenCpm(m)*wtmij;% i+1,j Nodew_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;sxxp_w_m(i+1,j)=sxxp_w_m(i+1,j)+sxxm(m)*wtmi1j;MIUp_w_m(i+1,j)=MIUp_w_m(i+1,j)+wtmi1j/mium(m);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);% Tw_m_t_T(i+1,j) = w_m_t_T(i+1,j) + DenCpm(m)*wtmi1j;% i,j+1 Nodew_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;sxxp_w_m(i,j+1)=sxxp_w_m(i,j+1)+sxxm(m)*wtmij1;MIUp_w_m(i,j+1)=MIUp_w_m(i,j+1)+wtmij1/mium(m);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);% Tw_m_t_T(i,j+1) = w_m_t_T(i,j+1) + DenCpm(m)*wtmij1;% i+1,j+1 Nodew_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;sxxp_w_m(i+1,j+1)=sxxp_w_m(i+1,j+1)+sxxm(m)*wtmi1j1;MIUp_w_m(i+1,j+1)=MIUp_w_m(i+1,j+1)+wtmi1j1/mium(m);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);% Tw_m_t_T(i+1,j+1) = w_m_t_T(i+1,j+1) + DenCpm(m)*wtmi1j1;endfor j = 1:Nx1for i = 1:Ny1if(w_m_p_node(i,j)>0)ETAP(i,j) = ETAp_w_m(i,j)/w_m_p_node(i,j);ETAPVP(i,j) = ETAP(i,j);%etaVP = eta in viscosity-elastic materials.MIUP(i,j) = w_m_p_node(i,j)/MIUp_w_m(i,j);sxx(i,j) = sxxp_w_m(i,j)/w_m_p_node(i,j);syy(i,j) = -sxx(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);%Khct(i,j) = DenCp_w_m_t(i,j)/w_m_p_node(i,j);%Den*Cpelsekht(i,j) = kht_backup(i,j);hct(i,j) = hct_backup(i,j);endtd(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);%TelseT0(i,j) = T_backup(i,j);endendend% Why??% Apply boundary condition after interpolate T0 from markers(any t).% top and left CtbcT0(1,2:Nx1-1) = 2*273-T0(2,2:Nx1-1);T0(Ny1,2:Nx1-1) = 2*1500-T0(Ny,2:Nx1-1);% left and right IbcT0(:,1) = T0(:,2);T0(:,Nx1) = T0(:,Nx);% Solve Stokes-equ and Continuity-equ.for j=1:1:Nx1for i=1:1:Ny1% Define global indexes in algebraic spacekvx=((j-1)*Ny1+i-1)*3+1; % Vxkvy=kvx+1; % Vykpm=kvx+2; % P% Vx equation External pointsif(i==1 || i==Ny1 || j==1 || j==Nx || j==Nx1)% Boundary Condition% 1*Vx=0coe(kvx,kvx)=1; % Left partR(kvx)=0; % Right part% Top boundaryif(i==1 && j>1 && j<Nx)coe(kvx,kvx+3)=bc_top; % Left partend% Bottom boundaryif(i==Ny1 && j>1 && j<Nx)coe(kvx,kvx-3)=bc_bottom; % Left partendelse% 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 pointsETA1=ETA(i-1,j);ETA2=ETA(i,j);ETAP1=ETAP(i,j);ETAP2=ETAP(i,j+1);% Density gradientsdRHOdx=(RHOVX(i,j+1)-RHOVX(i,j-1))/2/dx;dRHOdy=(RHOVX(i+1,j)-RHOVX(i-1,j))/2/dy;% Modify parametersETAPMIUij1 = 2*dt*MIUP(i,j+1)*ETAPVP(i,j+1)/(MIUP(i,j+1)*dt+ETAPVP(i,j+1));ETAPMIUij = 2*dt*MIUP(i,j)*ETAPVP(i,j)/(MIUP(i,j)*dt+ETAPVP(i,j)); ETAMIUi1j = 2*dt*MIU(i-1,j)*ETAVP(i-1,j)/(MIU(i-1,j)*dt+ETAVP(i-1,j));ETAMIUij = 2*dt*MIU(i,j)*ETAVP(i,j)/(MIU(i,j)*dt+ETAVP(i,j));ETAPij1 = ETAPVP(i,j+1)/(MIUP(i,j+1)*dt+ETAPVP(i,j+1));ETAPij = ETAPVP(i,j)/(MIUP(i,j)*dt+ETAPVP(i,j)); ETAij = ETAVP(i,j)/(MIU(i,j)*dt+ETAVP(i,j));ETAi1j = ETAVP(i-1,j)/(MIU(i-1,j)*dt+ETAVP(i-1,j));% Left partcoe(kvx,kvx-Ny1*3)=ETAPMIUij /2/dx^2 ; % Vx1(i,j-1)coe(kvx,kvx-3)=ETAMIUi1j/2/dy^2 ; % Vx2(i-1,j)coe(kvx,kvx)=-ETAPMIUij1/2/dx^2 - ETAPMIUij/2/dx^2 ...-ETAMIUij/2/dy^2-ETAMIUi1j/2/dy^2-g_x*dt*dRHOdx; % Vx3(i,j)coe(kvx,kvx+3)=ETAMIUij/2/dy^2 ; % Vx4(i+1,j)coe(kvx,kvx+Ny1*3)=ETAPMIUij1/2/dx^2 ; % Vx5(i,j+1)coe(kvx,kvy)=ETAPMIUij/2/dx/dy-ETAMIUij/2/dy/dx -dRHOdy*g_x*dt/4 ; % Vy2(i,j)coe(kvx,kvy+Ny1*3)=-ETAPMIUij1/2/dx/dy + ETAMIUij/2/dy/dx -dRHOdy*g_x*dt/4 ; % Vy4(i,j+1)coe(kvx,kvy-3)=-ETAPMIUij/2/dx/dy +ETAMIUi1j/2/dx/dy -dRHOdy*g_x*dt/4 ; % Vy1(i-1,j)coe(kvx,kvy+Ny1*3-3)=ETAPMIUij1/2/dx/dy -ETAMIUi1j/2/dx/dy -dRHOdy*g_x*dt/4 ; % Vy3(i-1,j+1)coe(kvx,kpm)=Kcont/dx; % P1(i,j)coe(kvx,kpm+Ny1*3)=-Kcont/dx; % P2(i,j+1)% Right part !!!R(kvx)=-RHOVX(i,j)*g_x - sxx(i,j+1)*ETAPij1/dx + sxx(i,j)*ETAPij/dx...-sxy(i,j)*ETAij/dy + sxy(i-1,j)*ETAi1j/dy ;end% Vy equation External pointsif(j==1 || j==Nx1 || i==1 || i==Ny || i==Ny1)% Boundary Condition% 1*Vy=0coe(kvy,kvy)=1; % Left partR(kvy)=0; % Right part% Left boundaryif(j==1 && i>1 && i<Ny)coe(kvy,kvy+3*Ny1)=bc_left; % Left partend% Right boundaryif(j==Nx1 && i>1 && i<Ny)coe(kvy,kvy-3*Ny1)=bc_right; % Left partendelse% 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 pointsETA1=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;ETAPMIUi1j = 2*dt*MIUP(i+1,j)*ETAPVP(i+1,j)/(MIUP(i+1,j)*dt+ETAPVP(i+1,j));ETAPMIUij = 2*dt*MIUP(i,j)*ETAPVP(i,j)/(MIUP(i,j)*dt+ETAPVP(i,j));ETAMIUij1 = 2*dt*MIU(i,j-1)*ETAVP(i,j-1)/(MIU(i,j-1)*dt+ETAVP(i,j-1));ETAMIUij = 2*dt*MIU(i,j)*ETAVP(i,j)/(MIU(i,j)*dt+ETAVP(i,j)); ETAPi1j = ETAPVP(i+1,j)/(MIUP(i+1,j)*dt+ETAPVP(i+1,j));ETAPij = ETAPVP(i,j)/(MIUP(i,j)*dt+ETAPVP(i,j)); ETAij = ETAVP(i,j)/(MIU(i,j)*dt+ETAVP(i,j));ETAij1 = ETAVP(i,j-1)/(MIU(i,j-1)*dt+ETAVP(i,j-1));% Left partcoe(kvy,kvy-Ny1*3)=ETAMIUij1/2/dx^2 ; % Vy1(i,j-1)coe(kvy,kvy-3)=ETAPMIUij/2/dy^2 ; % Vy2(i-1,j)coe(kvy,kvy)=-ETAPMIUij/2/dy^2-ETAPMIUi1j/2/dy^2 ...-ETAMIUij/2/dx^2-ETAMIUij1/2/dx^2 - dRHOdy*g_y*dt; % Vy3(i,j)coe(kvy,kvy+3)=ETAPMIUi1j/2/dy^2 ; % Vy4(i+1,j)coe(kvy,kvy+Ny1*3)=ETAMIUij/2/dx^2 ; % Vy5(i,j+1)coe(kvy,kvx)=ETAPMIUij/2/dy/dx - ETAMIUij/2/dx/dy -dRHOdx*g_y*dt/4; %Vx3(i,j)coe(kvy,kvx+3)=-ETAPMIUi1j/2/dy/dx + ETAMIUij/2/dx/dy -dRHOdx*g_y*dt/4; %Vx4(i+1,j)coe(kvy,kvx-Ny1*3)=-ETAPMIUij/2/dy/dx + ETAMIUij1/2/dx/dy -dRHOdx*g_y*dt/4; %Vx1(i,j-1)coe(kvy,kvx+3-Ny1*3)=ETAPMIUi1j/2/dy/dx - ETAMIUij1/2/dx/dy -dRHOdx*g_y*dt/4; %Vx2(i+1,j-1)coe(kvy,kpm)=Kcont/dy; % P1coe(kvy,kpm+3)=-Kcont/dy; % P2% Right partR(kvy)=-RHOVY(i,j)*g_y - syy(i+1,j)*ETAPi1j/dy + syy(i,j)*ETAPij/dy...- syx(i,j)*ETAij/dx + syx(i,j-1)*ETAij1/dx;end% P equation External pointsif(i==1 || j==1 || i==Ny1 || j==Nx1 ||...(i==2 && j==2))% Boundary Condition% 1*P=0coe(kpm,kpm)=1; % Left partR(kpm)=0; % Right part% Real BCif(i==2 && j==2)coe(kpm,kpm)=1*Kcont; %Left partR(kpm)=1e+9; % Right partendelse% Internal points: continuity eq.% dVx/dx+dVy/dy=0% Vy1% |% Vx1--P--Vx2% |% Vy2%% Left partcoe(kpm,kvx-Ny1*3)=-1/dx; % Vx1coe(kpm,kvx)=1/dx; % Vx2coe(kpm,kvy-3)=-1/dy; % Vy1coe(kpm,kvy)=1/dy; % Vy2% Right partR(kpm)=0;endendend% Computing the solution U:include vx,vy and pU = coe \ R;% Decompose the solution of global matrix. for j=1:1:Nx1for i=1:1:Ny1% Define global indexes in algebraic spacekvx=((j-1)*Ny1+i-1)*3+1; % Vxkvy=kvx+1; % Vykpm=kvx+2; % P% Reload solutionvx(i,j)=U(kvx);vy(i,j)=U(kvy);p(i,j)=U(kpm)*Kcont;endend% Restrict the time step limiting max_marker_disp by 0.5 grid step% Define dtm for transport of markers.% Why dt*1.2??dt = dt*dtkoef;maxvx=max(max(abs(vx)));maxvy=max(max(abs(vy)));if(dt*maxvx>dxymax*dx)dt=dxymax*dx/maxvx;endif(dt*maxvy>dxymax*dy)dt=dxymax*dy/maxvy;end% Computing strain rate(xy) and deviatoric(xy) stress components.for j = 1:Nxfor i = 1:Ny%ETAVP = ETA,which plastic is absent.Z = MIU(i,j)*dt/(MIU(i,j)*dt+ETAVP(i,j));srxy(i,j) = (1/2)*( (vx(i+1,j)-vx(i,j))/dy ...+ (vy(i,j+1)-vy(i,j))/dx ) ;% srxy is strain rate,sxy is stressxydsxy(i,j) = 2*ETAVP(i,j)*srxy(i,j)*Z + sxy(i,j)*(1-Z);endend% Computing strain rate(xx) and deviatoric(xx) stress components.for j = 2:Nxfor i = 2:Ny Z = MIUP(i,j)*dt/(MIUP(i,j)*dt+ETAPVP(i,j));srxx(i,j) = (1/2)*( (vx(i,j)-vx(i,j-1))/dx ...- (vy(i,j)-vy(i-1,j))/dy );% srxx is strain rate,sxx is stressxxdsxx(i,j) = 2*ETAPVP(i,j)*srxx(i,j)*Z + sxx(i,j)*(1-Z);endenddsxxij = dsxx - sxx;dsxyij = dsxy - sxy;% Computing subgrid diffusion for markersif (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-NODESw_Pnodes = zeros(Ny1,Nx1);w_nodes = zeros(Ny,Nx);DWsxxM = zeros(Ny1,Nx1);% Basic nodeDWsxyM = zeros(Ny,Nx);% 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 nodej=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;endif(i<1)i=1;elseif(i>Ny)i=Ny;end% Compute distancesdis_x=xm(m)-xp(j);dis_y=ym(m)-yp(i);% Compute weightswtmij=(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)sxxm0nodal= sxx(i,j)*wtmij + sxx(i+1,j)*wtmi1j ...+sxx(i,j+1)*wtmij1 + sxx(i+1,j+1)*wtmi1j1;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;endif(i<1)i=1;elseif(i>Ny-1)i=Ny-1;end% Compute distancesdis_x=xm(m)-x(j);dis_y=ym(m)-y(i);% Compute weightswtmij=(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)sxym0nodal= sxy(i,j)*wtmij + sxy(i+1,j)*wtmi1j ...+sxy(i,j+1)*wtmij1 + sxy(i+1,j+1)*wtmi1j1;% Eq.10.14dsxxm = sxxm0nodal - sxxm(m);dsxym = sxym0nodal - sxym(m);% Page209tdiff = etam(m)/mium(m);% Computing subgrid diffusion% Why??dve = 0.5;%209,0<dve<1sdif=-dve*dt/tdiff;if(sdif<-30)sdif=-30;end% Page162 Eq.10.16dsxxmsub=dsxxm*(1-exp(sdif));dsxymsub=dsxym*(1-exp(sdif));% Eq.10.19% Tcorrected = Tm0 + dTsubgrid (+ dTremain)sxxm(m) = sxxm(m) + dsxxmsub;sxym(m) = sxym(m) + dsxymsub;% Interpolate Basic node(dTijsub) from L-markers(dTmsub)% Eq.10.17% Define i,j indexes for the upper left nodej=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;endif(i<1)i=1;elseif(i>Ny)i=Ny;end% Compute distancesdis_x=xm(m)-xp(j);dis_y=ym(m)-yp(i);% Compute weightswtmij=(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);DWsxxM(i,j) = DWsxxM(i,j) + dsxxmsub*wtmij;w_Pnodes(i,j) = w_Pnodes(i,j) + wtmij;DWsxxM(i+1,j) = DWsxxM(i+1,j) + dsxxmsub*wtmi1j;w_Pnodes(i+1,j) = w_Pnodes(i+1,j) + wtmi1j;DWsxxM(i,j+1) = DWsxxM(i,j+1) + dsxxmsub*wtmij1;w_Pnodes(i,j+1) = w_Pnodes(i,j+1) + wtmij1;DWsxxM(i+1,j+1) = DWsxxM(i+1,j+1) + dsxxmsub*wtmi1j1;w_Pnodes(i+1,j+1) = w_Pnodes(i+1,j+1) + wtmi1j1;% Define i,j indexes for the upper left nodej=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;endif(i<1)i=1;elseif(i>Ny-1)i=Ny-1;end% Compute distancesdis_x=xm(m)-x(j);dis_y=ym(m)-y(i);% Compute weightswtmij=(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);DWsxyM(i,j) = DWsxyM(i,j) + dsxymsub*wtmij;w_nodes(i,j) = w_nodes(i,j) + wtmij;DWsxyM(i+1,j) = DWsxyM(i+1,j) + dsxymsub*wtmi1j;w_nodes(i+1,j) = w_Pnodes(i+1,j) + wtmi1j;DWsxyM(i,j+1) = DWsxyM(i,j+1) + dsxymsub*wtmij1;w_nodes(i,j+1) = w_nodes(i,j+1) + wtmij1;DWsxyM(i+1,j+1) = DWsxyM(i+1,j+1) + dsxymsub*wtmi1j1;w_nodes(i+1,j+1) = w_nodes(i+1,j+1) + wtmi1j1;enddsxxijsub = zeros(Ny1,Nx1);dsxyijsub = zeros(Ny,Nx);for j = 1:Nx1for i = 1:Ny1if(w_Pnodes(i,j)>0)% Compting dTijsub% Eq.10.17dsxxijsub(i,j) = DWsxxM(i,j)/w_Pnodes(i,j);if(i~=Ny1 && j ~=Nx1)dsxyijsub(i,j) = DWsxyM(i,j)/w_nodes(i,j);endendendend% Eq.10.18% Update dTijdsxxijremain = dsxxij - dsxxijsub;dsxyijremain = dsxyij - dsxyijsub;end% Interpolate L-markers(dTremain) from P-node(dsxxijremain)for m = 1:marknum% Define i,j indexes for the upper left nodej=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;endif(i<1)i=1;elseif(i>Ny)i=Ny;end% Compute distancesdis_x=xm(m)-xp(j);dis_y=ym(m)-yp(i);% Compute weightswtmij=(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)dsxxmremain=dsxxijremain(i,j)*wtmij + dsxxijremain(i+1,j)*wtmi1j ...+dsxxijremain(i,j+1)*wtmij1 + dsxxijremain(i+1,j+1)*wtmi1j1;if(t>1)sxxm(m) = sxxm(m) + dsxxmremain;else% Interpolate new temperature for the 1st timestepsxxm(m) = dsxx(i,j)*wtmij + dsxx(i+1,j)*wtmi1j ...+dsxx(i,j+1)*wtmij1 + dsxx(i+1,j+1)*wtmi1j1;end% Define i,j indexes for the upper left nodej=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;endif(i<1)i=1;elseif(i>Ny-1)i=Ny-1;end% Compute distancesdis_x=xm(m)-x(j);dis_y=ym(m)-y(i);% Compute weightswtmij=(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)dsxymremain=dsxyijremain(i,j)*wtmij + dsxyijremain(i+1,j)*wtmi1j ...+dsxyijremain(i,j+1)*wtmij1 + dsxyijremain(i+1,j+1)*wtmi1j1;% Tcorrected = Tm0(Tm0+ dTsubgrid) + dTremain% Eq.10.19if(t>1)sxym(m) = sxym(m) + dsxymremain;else% Interpolate new temperature for the 1st timestepsxym(m) = dsxy(i,j)*wtmij + dsxy(i+1,j)*wtmi1j ...+dsxy(i,j+1)*wtmij1 + dsxy(i+1,j+1)*wtmi1j1;endendfor j = 2:Nxfor i = 2:Ny pxyxy(i,j) = (1/4)*(dsxy(i,j)^2/ETAVP(i,j)+dsxy(i-1,j)^2/ETAVP(i-1,j)...+ dsxy(i,j-1)/ETAVP(i,j-1)+dsxy(i-1,j-1)/ETAVP(i-1,j-1));% Shear heating HsHsij(i,j) = dsxx(i,j)^2/ETAPVP(i,j)+ pxyxy(i,j);endend% Adiabatic heating Hafor j = 2:Nx1for i = 2:Ny1kvy = ( vy(i,j)+vy(i-1,j) )/2;kvx = ( vx(i,j)+vx(i,j-1) )/2;kRHOY = ( RHOVY(i,j)+RHOVY(i-1,j) )/2;kRHOX = ( RHOVX(i,j)+RHOVX(i,j-1) )/2;Haij(i,j) = T0(i,j)*Aefij(i,j)*RHOP(i,j)*(g_y*kvy+g_x*kvx);endend% Computing average velocity components at cell centres.for j = 1:Nx1for i = 1:Ny1% Computing internal Vx and Vy at the P nodes.if (j~=1 && i ~=1 && j~=Nx1 && i~=Ny1)% Vx at internal(P) nodes.pvx(i,j) = ( vx(i,j) + vx(i,j-1) )./2; % Vy at internal nodes.pvy(i,j) = ( vy(i,j) + vy(i-1,j) )./2; endendend% Why it is free-slip ??% Applying free-slip boundary condition for velocity in the P nodes .for j = 1:Nx1for i = 1:Ny1% j=1 j=2 j=3 ... j=Nx_node+1% i=1% i=2% i=3% ...% i=Ny_node+1if ( i==1 || i==Ny )% Top and Bottom% Vxif(j>1 && j<Nx1)pvx(i,j) = pvx(i+1,j);end% Vypvy(i,j) = pvy(i+1,j);endif (j==1 || j==Nx)% Right and Left% Vyif(i>1 && i< Ny1)pvy(i,j) = pvy(i,j+1);end% Vxpvx(i,j) = pvx(i,j+1);endendendfor titer = 1:1:4% Ctbc% Why???% Up and downfor j = 2:Nxi = 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) = 273*2;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) = 1500*2;end% left and rightfor i = 1:Ny1j = 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 pointsfor j = 2:Nxfor i = 2:Nyk = (j-1)*Ny1 + i; % k is the indexif vx(i,j)>0dtdx = (T0(i,j)-T0(i,j-1))/dx;elseif vx(i,j)<0dtdx = (T0(i,j+1)-T0(i,j))/dx;endif vy(i,j)>0dtdy = (T0(i,j)-T0(i-1,j))/dy;elseif vy(i,j)<0dtdy = (T0(i+1,j)-T0(i,j))/dy;end% KVX == Kh;KVY == Kvkhij = KVX(i,j);khij1 = KVX(i,j-1);kvij = KVY(i,j);kvi1j = KVY(i-1,j);% Not 2* why?? page177coet(k, k) = (khij+khij1)/dx^2+hct(i,j)/dt...+(kvij+kvi1j)/dy^2;% middle T3 i,jcoet(k, k-1) = -(kvi1j)/(2.*dy.^2); % up T2 i-1,jcoet(k, k+1) = -(kvij)/(2.*dy.^2); % down T4 i+1,j coet(k, k-Ny1) = -(khij1)/(dx.^2); % left T1 i,j-1coet(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)*hct(i,j)/dt + Hrij(i,j) + Haij(i,j) + Hsij(i,j); %- hctdt(i,j)*(vx(i,j)*dtdx+vy(i,j)*dtdy);endend%direct method:right martix is divided by coe martix on the leftu = coet \ bt;% transform vector to grids% U has Ny1 rows and Nx1 columnsU = reshape(u, Ny1, Nx1);T1 = U;% Computing dTi,j = dTi,j(subgrid) + dTi,j(remaining)% Eq10.13dTij = T1-T0;% Adjust dt for must condition:DTmax,and max number of iteration % is 2 % Define dtT for temperature equation.% dtT <= dtm(markers instead of stokes, visible in LINE-727)maxDTcurrent = max(max(abs(dTij)));if(titer<4 && maxDTcurrent>DTmax)dt=dt/maxDTcurrent*DTmax;elsebreak;endend% Computing subgrid diffusion for markersif (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-NODESw_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 nodej=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;endif(i<1)i=1;elseif(i>Ny)i=Ny;end% Compute distancesdis_x=xm(m)-xp(j);dis_y=ym(m)-yp(i);% Compute weightswtmij=(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.14dTm = Tm0nodal - Tm0(m);% Page163tdiff = DenCpm(m)/Km(m)/(2/dx^2+2/dy^2);% Computing subgrid diffusion% Why??sdif=-subgrdifcoe*dt/tdiff;if(sdif<-30)sdif=-30;end% Page162 Eq.10.16dTmsub=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.17DWTM(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;enddTijsub = zeros(Ny1,Nx1);for j = 1:Nx1for i = 1:Ny1if(w_Pnodes(i,j)>0)% Compting dTijsub% Eq.10.17dTijsub(i,j) = DWTM(i,j)/w_Pnodes(i,j);endendend% Eq.10.18% Update dTijdTijremain = dTij - dTijsub;end% Interpolate L-markers(dTremain) from P-node(dTijremain)for m = 1:marknum% Define i,j indexes for the upper left nodej=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;endif(i<1)i=1;elseif(i>Ny)i=Ny;end% Compute distancesdis_x=xm(m)-xp(j);dis_y=ym(m)-yp(i);% Compute weightswtmij=(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.19if(t>1)Tm0(m) = Tm0(m) + dTmremain;else% Interpolate new temperature for the 1st timestepTm0(m) = T1(i,j)*wtmij + T1(i+1,j)*wtmi1j ...+T1(i,j+1)*wtmij1 + T1(i+1,j+1)*wtmi1j1;endend% Figuresubplot(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;colorbartitle(['log10(ETA) after ',num2str(t),' deltat'])% Hasubplot(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);% Rotate stress by Analytical formulas.w = zeros(Ny,Nx);wm = zeros(1,marknum);for j = 1:Nxfor i = 1:Ny%ETAVP = ETA,which plastic is absent.w(i,j) = (1/2)*( (vy(i,j+1)-vy(i,j))/dx ...- (vx(i+1,j)-vx(i,j))/dy ) ;endend% Interpolate wij from basic node to L-markersfor m=1:1:marknumj=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;endif(i<1)i=1;elseif(i>Ny-1)i=Ny-1;end% Compute distancesdis_x=xm(m)-x(j);dis_y=ym(m)-y(i);% Compute weightswtmij=(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 vxm and vym velocitywm(m) = w(i,j)*wtmij+w(i+1,j)*wtmi1j...+ w(i,j+1)*wtmij1+w(i+1,j+1)*wtmi1j1;angleAB = wm(m)*dt;sxymnew = (1/2)*(sxxm(m)-syym(m))*sin(2*angleAB)+sxym(m)*cos(2*angleAB);sxxmnew = sxxm(m)*(cos(angleAB))^2+syym(m)*(sin(angleAB))^2-sxym(m)*sin(2*angleAB) ;sxym(m) = sxymnew;sxxm(m) = sxxmnew;end% 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;endif(i<1)i=1;elseif(i>Ny)i=Ny;end% Compute distancesdis_x=xm(m)-xp(j);dis_y=ym(m)-yp(i);% Compute weightswtmij=(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 vxm and vym velocityvxm(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 nodej=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 distancesdis_x=xm(m)-xvx(j);dis_y=ym(m)-yvx(i);% Compute weightswtmij=(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 velocityvxcij = (1/3)*(2*vx(i,j));vxci1j = (1/3)*(2*vx(i+1,j));vxcij1 = (1/3)*(2*vx(i,j+1));vxci1j1 = (1/3)*(2*vx(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 nodej=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;endif(i<1)i=1;elseif(i>Ny-1)i=Ny-1;end% Compute distancesdis_x=xm(m)-xvy(j);dis_y=ym(m)-yvy(i);% Compute weightswtmij=(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 velocityvycij = (1/3)*(2*vy(i,j));vyci1j = (1/3)*(2*vy(i+1,j));vycij1 = (1/3)*(2*vy(i,j+1));vyci1j1 = (1/3)*(2*vy(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) = vxm(rk)*dt/2 + xa;ym(m) = vym(rk)*dt/2 + ya;elseif(rk == 3)xm(m) = vxm(rk)*dt + xa;ym(m) = vym(rk)*dt + ya;endend% Return original coordinate.xm(m) = xa;ym(m) = ya;% Move markers with effective velocity field.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
