Problem with 3D Meshgrid - Not understanding why some planes are 0 which results in my dx becoming 0

2 次查看（过去 30 天）
Elias
Elias 2024-5-6
评论： William Rose 2024-5-7
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!

请先登录，再回答此问题。

采纳的回答

William Rose
William Rose 2024-5-6
Ran in:
@Elias Larsson,
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)
X =
X(:,:,1) = 0 1 2 3 0 1 2 3 0 1 2 3 X(:,:,2) = 0 1 2 3 0 1 2 3 0 1 2 3
Y =
Y(:,:,1) = 10 10 10 10 12 12 12 12 14 14 14 14 Y(:,:,2) = 10 10 10 10 12 12 12 12 14 14 14 14
Z =
Z(:,:,1) = 20 20 20 20 20 20 20 20 20 20 20 20 Z(:,:,2) = 23 23 23 23 23 23 23 23 23 23 23 23
Size of X:
size(X)
ans = 1x3
3 4 2
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
The storage order of X is illustrated by
reshape(X,[1,24])
ans = 1x24
0 0 0 1 1 1 2 2 2 3 3 3 0 0 0 1 1 1 2 2 2 3 3 3
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
The output of the preceding command extends beyond the window, for me, but you get the idea.
  2 个评论
William Rose
William Rose 2024-5-7
编辑：William Rose 2024-5-7
Ran in:
@Elias Larsson,
[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
William Rose 2024-5-7
@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.

请先登录，再进行评论。

更多回答（1 个）

Cris LaPierre
Cris LaPierre 2024-5-6
编辑：Cris LaPierre 2024-5-6
Ran in:
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)
X =
X(:,:,1) = 1 2 3 1 2 3 1 2 3 X(:,:,2) = 1 2 3 1 2 3 1 2 3
Y =
Y(:,:,1) = 6 6 6 7 7 7 8 8 8 Y(:,:,2) = 6 6 6 7 7 7 8 8 8
Z =
Z(:,:,1) = 10 10 10 10 10 10 10 10 10 Z(:,:,2) = 20 20 20 20 20 20 20 20 20
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.

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by