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 .

[Sajani]

@Torsten So that’s why you suggested working with CLAWPACK. Do you think there’s any possibility to handle this properly in MATLAB, or is it too difficult compared to CLAWPACK? I need to make the code in MATLAB/Python.

[Torsten]

The shallow water equations are a system of hyperbolic differential equations that are hard to solve numerically - especially if the initial conditions are discontinuous. Since you seem to have little experience solving this kind of equations, I suggested CLAWPACK from which I know it works properly.

I didn't test them, but there are also MATLAB codes on the File Exchange that could help you:

https://uk.mathworks.com/matlabcentral/fileexchange/46475-1d-shallow-water-equations-dam-break?s_tid=prof_contriblnk

https://uk.mathworks.com/matlabcentral/fileexchange/68419-1-d-muscl-solver-for-the-shallow-water-equations

https://uk.mathworks.com/matlabcentral/fileexchange/30126-2d-shallow-water-equation

https://uk.mathworks.com/matlabcentral/fileexchange/47126-solving-shallow-water-equations-using-finite-volume-methods

[William Rose]

@Sajani,

I am not familiar with these equatons. @Torsten clearly is. Before I comment on your proposed code for solving the equations, I want to be sure I understand the equations. I assume Q=Q(x,t)=flow, and = free surface elevation, and A=A(x,t) = area of the flow. Am I correct?

The equations in your code do not match the equations in your original question. As @Torsten said, the update equation for Q in the code is missing the term , which should be present on the right hand side. It is also missing the term , which should be present on the right hand side. Are you assuming that the term equals zero?

The equations in your original question include A, B, and Cd. None of these appear in your code. Please explain what they are, and what the units are. My guess is that B is a constant and equals the channel width, and therefore . But the update equaitons in your code only work dimensionally if Q is flow per unit width, i.e. the update equaitons in the code are dimensionally corect only if Q has units of m^2/s. But if that is the case then what is B in your equaitons? I also infer that Cd is a dimensionless constant that inlcudes effects of fluid density. Perhaps you are assuming this term can be ignored. Are my guesses above correct?

Your code includes h=1. This constant is labelled labelled "depth". Please explain.

Your code includes

% 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;

I assume the variables and constants have the following dimensions:

dt=[s], dx=[m], g=[m/s^2], h=[m], eta=[m], and deta_dx is dimensionless.

The update equations include terms on the right hand side with different dimensions (which is an error), if the units for Q are [m^3/s]. But if the units for Q are [m^2/s] (i.e. Q=flow pre unit width), then the update equations are dimensionally correct. Then h must not be channel width, and then I do not understand the relationship between A, B, and η in your original equaitons. Please explain.

[Sajani]

@William Rose,

Actually, I ignored the nonlinear terms because I just want to check if my steps are correct. Here is the modified code. And these are the bouandary and initial conditions;

Boundary: at & at (head)

Initial: at

L = 12000; %Total length
dx = 500;
Nx = L/dx +1; 
B = 1600;  % lagoon width
H = 1.1;  % mean depth
g = 9.81;
Cd = 0.0022; % drag coefficient
% dt = 0.5*dx/(sqrt(g*H)); % according to the CFL condition
dt = 60;
T = 10*60; 
Nt = round(T/dt)+1;
 
x_eta = linspace(0,L,Nx+1);        
x_Q   = x_eta(1:end-1) + dx/2; 

eta = zeros(Nx+1,Nt); % eta free surface elevation
Q   = zeros(Nx,Nt); % cross sectionally averaged volume flux 

for n = 1:Nt-1
    t = n*dt;
    
    % Boundary conditions
    Q(1,:) = 0.5; % lagoon head, equations give zeros when I put zero here
    eta(end,:) = [0.135083512 0.109083512 0.118083512 0.126083512 0.108083512 0.129083512 0.124083512 0.126083512 0.108083512 0.111083512 0.090083512]; % sea level at mouth
    
    % Update Q (momentum)
    for i = 2:Nx-1
        A(i,n) = B* (H+eta(i,n));
        Q(i,n+1) = Q(i,n) - (dt/dx)*g*A(i,n)*(eta(i+1,n) - eta(i,n))- (Q(i,n)/A(i,n))*(Q(i,n) - Q(i-1,n))/dx - B*Cd*(Q(i,n)/A(i,n))*abs(Q(i,n)/A(i,n));
    end
    
    % Update eta (continuity)
    for i = 2:Nx
        eta(i,n+1) = eta(i,n) - dt/B*(Q(i,n+1) - Q(i-1,n+1))/dx;
    end
end

[Sajani]

@Torsten,

Thank you for your guidance. I’ll go through those codes.

[Torsten]

I didn't check your code carefully, but what immediately strikes me is that you never reference eta(Nx+1,:) . This means that the boundary condition for eta at x=L never "progresses" into the domain.

Further, you compute eta(i,n+1) using Q(:,n+1). You should use an explicit time step here, thus use Q(:,n) instead. Otherwise, your time integration scheme is some mixture of explicit and implicit scheme - hard to say whether it will produce something useful.

And in all sources I searched the non-conservative form of the equations reads

dh/dt + u*dh/dx + h*du/dx = 0
du/dx + u*du/dx = -g*dh/dx + some sink term

Maybe you could explain how this transforms to the system you wrote in your original question.