Does anyone have a MATLAB code example for a staggered grid (for 1D/2D problems)?

I’m working on solving the shallow water equations using a staggered grid in MATLAB, and I’m looking for example codes or guidance for a staggered grid. Specifically, I want to understand how to set up the grid (interfaces and centers), apply boundary/initial conditions. These are my equations,

[Torsten]

I know the use of staggered grids to compute velocity and pressure for the Navier-Stokes equations.

Why do you think this is necessary for the shallow water equations ?

I suggest working with CLAWPACK because the equations are quite difficult to handle numerically:

https://www.clawpack.org/riemann_book/html/Shallow_water.html

[Sajani]

@Torsten Thank you for your response. I’m using the staggered grid because many references have applied it previously, and I’m trying to implement a simple version in MATLAB before moving to more advanced. If you have a MATLAB code for staggered grids applied to the Navier–Stokes equations, could you please share it? I’ll also check out CLAWPACK as you suggested.

Answer 1

[William Rose]

@Sajani,

Here is a start, just to show how the grids set up in 2D:

dx=0.1; Lx=1;
dy=0.1; Ly=1;
[PX,PY]=meshgrid(0:dx:Lx,0:dy:Ly);
[UX,UY]=meshgrid(dx/2:dx:Lx-dx/2,0:dy:Ly);
[VX,VY]=meshgrid(0:dx:Lx,dy/2:dy:Ly-dy/2);

Display the grids

figure;
plot(PX(:),PY(:),'r.',UX(:),UY(:),'g.',VX(:),VY(:),'b.')
legend('Pressure','U','V'); xlabel('X'); ylabel('Y'); title('Staggered Grids')

For 3D, add a grid for W (z-component of velocity) which is offset by dz/2.

Depending on your boundary conditions, you might want the U and V grid opoints to lie exactly on both boundaries. For example, in a tank, U and V are zero at the edges, so place the U and V grid points on the tank edges, where the respective velocities must be zero. The boundary conditions will include: {U=0 at X=0 and at X=Lx}; {V=0 at Y=0 and at Y=Ly}. One way to make such a grid is shown below.

[UX,UY]=meshgrid(0:dx:Lx,dy/2:dy:Ly-dy/2);
[VX,VY]=meshgrid(dx/2:dx:Lx-dx/2,0:dy:Ly);
[PX,PY]=meshgrid(dx/2:dx:Lx-dx/2,dy/2:dy:Ly-dy/2);

Display the grids

figure;
plot(PX(:),PY(:),'r.',UX(:),UY(:),'g.',VX(:),VY(:),'b.')
legend('Pressure','U','V'); xlabel('X'); ylabel('Y'); title('Staggered Grids')

For the specific equations on the staggered grids, consult the journal articles to which you referred.

[William Rose]

@Sajani,

[Edit: Fix punctuation. No changes to code.]

They key features of the examples above are:

• The U grid points have the same y-values as the pressure points, and are half-offset in the x-direction.
• The V grid points have the same x-values as the pressure points, and are half-offset in the y-direction.

In 3D:

• The U grid points have the same y-values and z-values as the pressure points, and are half-offset in the x-direction.
• The V grid points have the same x-values and z-values as the pressure points, and are half-offset in the y-direction.
• The W grid points have the same x-values and y-values as the pressure points, and are half-offset in the z-direction.

Code for the 3D case, like the second 2D example in my answer above.

dx=0.1; Lx=1;
dy=0.1; Ly=1.1;  % see comment about grid dimensions below * 
dz=0.1; Lz=1.2;
[UX,UY,UZ]=meshgrid(0:dx:Lx,dy/2:dy:Ly-dy/2,dz/2:dz:Lz-dz/2);
[VX,VY,VZ]=meshgrid(dx/2:dx:Lx-dx/2,0:dy:Ly,dz/2:dz:Lz-dz/2);
[WX,WY,WZ]=meshgrid(dx/2:dx:Lx-dx/2,dy/2:dy:Ly-dy/2,0:dz:Lz);
[PX,PY,PZ]=meshgrid(dx/2:dx:Lx-dx/2,dy/2:dy:Ly-dy/2,dz/2:dz:Lz-dz/2);

The 3D grid above will allow easy enforcement of boundary conditions: {velocity=0 normal to the tank edges, at the edges}, because there are velocity nodes of the desired direction on the appropriate edges. This would work for a rectangular tank that is closed on top and full of water. For a rectangular tank that is open on top, you might want a layer of pressure nodes at z=0 (the surface), so you can enforce the boundary condition {P=0 at z=0 (on the surface)}. To do this, modify the code above slightly:

[UX,UY,UZ]=meshgrid(0:dx:Lx,dy/2:dy:Ly-dy/2,0:dz:Lz);
[VX,VY,VZ]=meshgrid(dx/2:dx:Lx-dx/2,0:dy:Ly,0:dz:Lz);
[WX,WY,WZ]=meshgrid(dx/2:dx:Lx-dx/2,dy/2:dy:Ly-dy/2,dz/2:dz:Lz+dz/2);
[PX,PY,PZ]=meshgrid(dx/2:dx:Lx-dx/2,dy/2:dy:Ly-dy/2,0:dz:Lz);

In this case, the W nodes at the bottom of the tank (where W=0) are at level z=Lz+dz/2.

* When I am making 2D or 3D meshes, I like my grids to have slightly different numbers of elements in each direction. If they're the same size in all directions, one can mix up which direction is plotted as which, and there probably won't be an error message. If they're different, therre is more likely to be an error message to alert me to the problem.

[Sajani]

@William Rose Thank you for your explanation. I understand it clearly now. What do you think about this code for the 1-D shallow water equations? I’m mainly looking for guidance on how to update Q and η. (or p,u, and v in your code)

L  = 10;         % domain length
Nx = 100;        % number of cells
dx = L/Nx;       
g  = 9.81;       
h  = 1.0;        % depth
dt = 0.01;       % time step
Nt = 200;        % number of time steps

% Grids
x_eta = dx/2:dx:L-dx/2;   % eta at cell centers
x_Q   = 0:dx:L;           % Q at interfaces

eta = zeros(size(x_eta)); 
Q   = zeros(size(x_Q));

% Time stepping 
for n = 1:Nt
    % Compute spatial derivatives
    deta_dx = diff(eta)/dx;             
    dQ_dx   = diff(Q)/dx;               
    
    % Update Q (momentum eq)
    Q(2:end-1) = Q(2:end-1) - dt*g*h*deta_dx;
    
    % Update eta (continuity eq)
    eta = eta - dt*dQ_dx;
    
    % Simple boundary conditions (closed walls)
    Q(1) = 0; 
    Q(end) = 0;
end

[Torsten]

You will come into trouble if you want to include the missing dQ/dx term in the equation for Q because dQ/dx is given at dx/2:dx:L-dx/2 while Q is given at 0:dx:L .