Adding Dirichlet condition leads to NAN #2
Description
When I modify the Laplacian matrix code to impose zero Dirichlet condition on the top boundary, the program leads to NAN. I believe the modification is correct assuming that the pressure array values P(x_i,y_j)=Pij are arranged into a Nx*Ny vector in "natural order":
Here is the modified code for the Laplacian matrix creation:
% Creat Laplacian operator for solving pressure Poisson equation
for j=1:ny %number of block matrices along y
for i=1:nx %number of block matrices along x
q = i+(j-1)*nx;
%if(mod(q,1*nx)==0)
if(q > (nx-1)*ny) % if top boundary
L(q,q)=1; % zero Dirichlet condition
else
L(q,q)=-(2*dxi^2+2*dyi^2);
for ii=i-1:2:i+1
if (ii>0 && ii<=nx) % Interior points
L(q,ii+(j-1)*nx)=dxi^2;
else % Neuman conditions on boundary
L(q,q)= ...
L(q,q)+dxi^2;
end
end
for jj=j-1:2:j+1
if (jj>0 && jj<=ny) % Interior points
L(q,i+(jj-1)*nx)=dyi^2;
else % Neuman conditions on boundary
L(q,q)= ...
L(q,q)+dyi^2;
end
end
end
end
end
When Nx = Ny = 5, the Laplacian matrix looks like this, which I believe is correct:
Here is the complete matlab program I wrote following your guide:
clear; clc;
% grid parameters
Lx=1; % x length
Ly=1; % y length
nx=5; % number of grid points along x
ny=5; % number of grid Points along y
%simulation parameters
dt=0.001; % time step
%t_final=0.1; % simulation run time
t_final=0.2; % simulation run time
% physics parameters
nu=0.001; % kinematic Viscosity
rho=1.0; % density
% Index extents
imin=2;
imax=imin+nx-1; % ok
jmin=2; % ok
jmax=jmin+ny-1;
%create mesh sizes
dx=Lx/(nx);
dy=Ly/(ny);
%disp(["dx=",dx]);
%disp(["dy=",dy]);
% Create mesh sizes
dxi=1/dx;
dyi=1/dy;
% Number of timesteps
Nt=t_final/dt;
% preallocate important arrays
p=zeros(imax,jmax); % ok
%p=zeros(nx,ny); % ok
us=zeros(imax+1,jmax+1);
vs=zeros(imax+1,jmax+1);
% R=zeros(imax,1);
R=zeros(nx,1);
u=zeros(imax+1,jmax+1);
v=zeros(imax+1,jmax+1);
L=zeros(nx*ny,nx*ny);
% initial values
t=0; % initial time
u_bot=0; % Initial Velocity for bottom wall
u_top=1; % Initial Velocity for top wall
v_lef=0; % Initial Velocity for left wall
v_rig=0; % Initial Velocity for right wall
% Creat Laplacian operator for solving pressure Poisson equation
for j=1:ny %number of block matrices along y
for i=1:nx %number of block matrices along x
q = i+(j-1)*nx;
%if(mod(q,1*nx)==0)
if(q > (nx-1)*ny) % if top boundary
L(q,q)=1; % zero Dirichlet condition
else
L(q,q)=-(2*dxi^2+2*dyi^2);
for ii=i-1:2:i+1
if (ii>0 && ii<=nx) % Interior points
L(q,ii+(j-1)*nx)=dxi^2;
else % Neuman conditions on boundary
L(q,q)= L(q,q)+dxi^2;
end
end
for jj=j-1:2:j+1
if (jj>0 && jj<=ny) % Interior points
L(q,i+(jj-1)*nx)=dyi^2;
else % Neuman conditions on boundary
L(q,q)= L(q,q)+dyi^2;
end
end
end
end
end
%L(1,:)=0; L(1,1)=1; % since we have Dirichlet condition, we don't need this anymore
disp("L="); disp(num2str(L));
% solver
disp("Starting simulation");
while t <= t_final
% update time
t = t + dt;
disp("t:");
disp(t);
u_top=1;
% u Momentum Predictor
for j = jmin:jmax
for i = imin+1:imax
v_here = 0.25*(v(i-1,j)+v(i-1,j+1)+v(i,j)+v(i,j+1));
us(i,j)=u(i,j)+dt* ...
(nu*(u(i-1,j)-2*u(i,j)+u(i+1,j))*dxi^2 ...
+nu*(u(i,j-1)-2*u(i,j)+u(i,j+1))*dyi^2 ...
-u(i,j)*(u(i+1,j)-u(i-1,j))*0.5*dxi ...
-v_here*(u(i,j+1)-u(i,j-1))*0.5*dyi);
end
end
% v Momentum predictor
for j = jmin+1:jmax
for i = imin:imax
u_here = 0.25*(u(i,j-1)+u(i+1,j-1)+u(i,j)+u(i+1,j));
vs(i,j)=v(i,j)+dt* ...
(nu*(v(i-1,j)-2*v(i,j)+v(i+1,j))*dxi^2 ...
+nu*(v(i,j-1)-2*v(i,j)+v(i,j+1))*dyi^2 ...
-u_here*(v(i+1,j)-v(i-1,j))*0.5*dxi...
-v(i,j)*(v(i,j+1)-v(i,j-1))*0.5*dyi);
end
end
% Compute right-hand side (R) of the Pressure Poisson Equation using
% the predictor velocities (vs) & (us).
n=0;
for j=jmin:jmax
for i=imin:imax
n=n+1;
R(n)=-rho/dt*((us(i+1,j)-us(i,j))*dxi+(vs(i,j+1)-vs(i,j))*dyi);
%if(mod(n,nx)==0)
if(n > (nx-1)*ny) %if top boundary
R(n)=0; % force 0 value
end
%disp(["Time: ",t]);
disp(["R[",n,"]:",R(n)]);
end
end
% Find pressure solution for matrix inversion, Using Laplacian
% Operator (L) and Right hand side of Pressure Poisson Equations (R)
pv=L\R;
disp("pv=");
disp(pv);
n=0;
p=zeros(imax,jmax);
% Convert pressure into mesh
for j=jmin:jmax
for i=imin:imax
n=n+1;
%if(mod(n,nx)==0)
if(n > (nx-1)*ny) % if top boundary
pv(n)=0; % force zero Dirichlet
end
p(i,j)=pv(n);
end
end
disp("p="); disp(num2str(p));
% Corrector of x^(n+1) & y^(n+1) velocity
for j=jmin:jmax
for i=imin+1:imax
u(i,j)=us(i,j)-dt/rho*(p(i-1,j-1)-p(i-2,j-1))*dxi;
end
end
for j=jmin+1:jmax
for i=imin:imax
v(i,j)=vs(i,j)-dt/rho*(p(i-1,j-1)-p(i-1,j-2))*dyi;
end
end
% Set Boundary Conditions
u(:,jmin-1)=u(:,jmin)-2*(u(:,jmin)-u_bot);
u(:,jmax+1)=u(:,jmax)-2*(u(:,jmax)-u_top);
v(imin-1,:)=v(imin,:)-2*(v(imin,:)-v_lef);
v(imax+1,:)=v(imax,:)-2*(v(imax,:)-v_rig);
end;
disp("End of simulation");
disp("Last p="); disp(num2str(p));
disp("Last u="); disp(num2str(u));
disp("Last v="); disp(num2str(v));
disp("Finita!");
It leads to NAN result at around 0.006 s time.