Problem with 3D Meshgrid - Not understanding why some planes are 0 which results in my dx becoming 0
Show older comments
Hello,
For context: Trying to make a Navier Stokes solver in Matlab, starting with 1D just to see if it will work.
To start I am using the meshgrid do discretize the spatial domain into grids. This is done with the following code:
L = 1 % Length of the cube side
n = 100 % Num of gridpoints
x = linspace(0, L, n);
y = linspace(0, L, n);
z = linspace(0, L, n);
[X, Y, Z] = meshgrid(x,y,z);
Now this works and I get three 100x100x100 matrices for X, Y and Z. The problem Im having is understanding how it actually works in generating each individual element. I understand that it reallocates the points in a grid formation from the input values, which is how I want it. However when testing it with this code:
for t=1:T
for i = 2:n-1
for j = 2:n-1
for k = 2:n-1
dx = X(i+1,j,k) - X(i,j,k);
dy = Y(i+1,j,k) - Y(i,j,k);
dz = Z(i+1,j,k) - Z(i,j,k);
if dx <= 1e-10
continue;
else
disp(dx)
end
end
end
u_convection = -0.5/dx * (U(i+1)^2 - U(i-1)^2);
u_diffiusion = nu/dx^2 * (U(i+1) - 2*U(i) + U(i-1));
U_new(i) = U(i) + dt*(u_convection + u_diffiusion);
end
end
U = U_new;
end
The dx terms result in 0 for all the iterations. At this point im just trying for 1D just to make sure it works so dont mind the nested loops, those are for later usage.
From what I can see the X matrix has alot of concecutive values, like X(1,88,56) through to X(99,88,56) all have the value 0.8788. Thus it is logical that dx results in 0.
Could somone please explain why we have the same values in the same "vector" of X, as I dont really understand the documentation of gridmesh(). Also if anybody has any ideas to make the dx term work that would be appriacated as well.
Thanks in advance!
Accepted Answer
For clarity, consider an example in which dx,dy,dz differ, and in which the number of elements differ for each dimension.
x=0:3; y=10:2:14; z=20:3:23; % dx=1, dy=2, dz=3
[X,Y,Z]=meshgrid(x,y,z)
Size of X:
size(X)
The storage order of X is illustrated by
reshape(X,[1,24])
The output of the preceding command extends beyond the window, for me, but you get the idea.
2 Comments
[Edit: Fix mistake in equation for dcdtConv. I used dt where I should have used dx.]
You use c()^2 in the convection term, but it should be just c().
For the 1-D convection-diffusion equaiton, using an array with 3 spatial dimensions, I define a 4-D array for concentration (3 spatial + 1 time). I do not use meshgrid.
dx=0.05; dy=0.5; dz=0.5; % [cm]
D=1; % diffusion constant [cm^2/s]
v=20; % velocity [cm/s]
Stability criterion for diffusion equation, using the explicit forward method:
dt<=(dx)^2/(2*D),
where D=diffusion constant. Therefore we choose dt<=0.00125:
dt=0.001; % [s]
x=0:dx:1; y=0:dy:1; z=0:dz:1; % position [cm]
t=0:dt:0.1; % time [s]
nx=length(x); ny=length(y); nz=length(z); nt=length(t);
c=zeros(nx,ny,nz,nt); % concentration [cm^-3]
% Next: make c=1 at t=0 and x=0.5:
c(11,:,:,1)=ones(ny,nz);
We have initialized the system.
Solve by looping through time and space. I loop from i=2:nx-1 so that there will not be errors at c(i-1) or at c(i+1). This enforces boundary conditions: c=0 at i=1 and at i=nx, at all times. The solution is the same at all points in each y-z plane within the 4D array.
for h=2:nt
for i=2:nx-1
dcdtDiff=D*(c(i-1,:,:,h-1)-2*c(i,:,:,h-1)+c(i+1,:,:,h-1))/(dx^2);
dcdtConv=v*(c(i-1,:,:,h-1)-c(i+1,:,:,h-1))/(2*dx);
dcdt=dcdtDiff+dcdtConv;
c(i,:,:,h)=c(i,:,:,h-1)+dt*dcdt;
end
end
Plot results:
figure
colors=['k','r','g','b','c','m'];
for i=1:6
plot(x,c(:,1,1,i),'.-','Color',colors(i))
hold on
legstr{i}=sprintf('t=%.3f',t(i));
end
legend(legstr); grid on
xlabel('Position (cm)'); ylabel('Concentration (1/cm^3)')
title('Concentration vs Position and Time')
The solution above is with D=1 cm^2/s and v=20 cm/s. The combined effects of diffusion and convection are evident in the plot.
This little demo shows that it works.
William Rose
on 7 May 2024
@Elias Larsson, Just as there is a stability criterion for the diffusion equation, there is also a stability criterion for the convection equation: the Peclet number, Pe=vL/D, where L=dx, should be <=1, if you use the explicit scheme shown above.
The solution in my post above uses D=1, v=20, dx=0.05. Therefore Pe=1, and the solution appears stable. But if you use v=50, then Pe=2.5, and you will see that the numerical solution, using the algorithm above, is unstable.
More Answers (1)
Cris LaPierre
on 6 May 2024
Edited: Cris LaPierre
on 6 May 2024
To simplify, think of a 3D grid with the following coordinates
- X location can either be 1, 2, or 3
- Y location can either be 6, 7, or 8
- Z location can either be 10 or 20
Each dimension of a 3D array corresponds to one of the 3D axes. meshgrid creates a grid of all combinations of these coordinates in (x,y,z) space and returns the corresponding coordinates in X, Y and Z.
- colums correspond to location along X axis
- rows correspond to location along Y axis
- sheets correspond to location along Z axis
You'll see that each column in X contains the same value because it corresponds to the same X location.
Each row in Y contains the same value because it corresponds to the same Y location.
Each sheet in Z contains the same value because they correspond to the same Z location.
x = 1:3;
y = 6:8;
z = [10 20];
[X, Y, Z] = meshgrid(x,y,z)
Another way to think about this is to plot the X, Y, Z values and see how they are arranged in a 3D coordiate system.
scatter3(X(:),Y(:),Z(:),'filled')
xlabel('X'); ylabel('Y'); zlabel('Z');
So an indexing operation of (:,1,1) for each variable will return the coordinates of (1,6,10), (1,7,10) and (1,8,10), or the 3 dots on the Y axis.
2 Comments
Elias
on 7 May 2024
William Rose
on 7 May 2024
@Elias, you're welcome. Good luck with your work.
Categories
Find more on Image Arithmetic in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
