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

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!

Answer 1

[William Rose]

@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.

[William Rose]

@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]

@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.

Answer 2

[Cris LaPierre]

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.

[Elias]

Alright its begining to make sense, and you actuallt gave me an idea! Big thanks for the response!

[William Rose]

@Elias, you're welcome. Good luck with your work.