Implementation of 2D advection equation with second order upwind scheme?

43 views (last 30 days)
Hi all,
I am trying to numerically discretize a 2D advection equation to model the transport of rocks with thickness (h_debris) on top of glacier ice with velocity components (velx_mod and vely_mod). The equation to be discretized is as follows:
dh/dt = -d(u*h)/dx - d(v*h)/dy = -u(dh/dx) - v(dh/dy) - h(du/dx) - h(dv/dy)
% in the model:
% u = velx_mod (x component of ice surface velocity)
% v = vely_mod (y component of ice surface velocity)
% h = h_debris (thickness of the rock layer)
I therefore did the following for a second order upwind scheme in 2D:
if velx_mod(i,j) < 0
term1_debris_x(i,j) = -velx_mod(i,j).*((-h_debris(i+2,j)+4*h_debris(i+1,j)-3*h_debris(i,j))./(2*deltax_d)) - h_debris(i,j).*((-velx_mod(i+2,j)+4*velx_mod(i+1,j)-3*velx_mod(i,j))./(2*deltax_d));
elseif velx_mod(i,j) >= 0
term1_debris_x(i,j) = -velx_mod(i,j).*((3*h_debris(i,j)-4*h_debris(i-1,j)+h_debris(i-2,j))./(2*deltax_d)) - h_debris(i,j).*((3*velx_mod(i,j)-4*velx_mod(i-1,j)+velx_mod(i-2,j))./(2*deltax_d));
if vely_mod(i,j) < 0
term1_debris_y(i,j) = -vely_mod(i,j).*((-h_debris(i,j+2)+4*h_debris(i,j+1)-3*h_debris(i,j))./(2*deltax_d)) - h_debris(i,j).*((-vely_mod(i,j+2)+4*vely_mod(i,j+1)-3*vely_mod(i,j))./(2*deltax_d));
elseif vely_mod(i,j) >= 0
term1_debris_y(i,j) = -vely_mod(i,j).*((3*h_debris(i,j)-4*h_debris(i,j-1)+h_debris(i,j-2))./(2*deltax_d)) - h_debris(i,j).*((3*vely_mod(i,j)-4*vely_mod(i,j-1)+vely_mod(i,j-2))./(2*deltax_d));
h_debris(i,j) = h_debrisini(i,j) + deltat_d*(term1_debris_x(i,j) + term1_debris_y(i,j));
I am not an expert in this field and this is the first time I implement this advection equation in 2D. Could someone check my code and assure I did the right thing here? Thanks already.
Yoni Verhaegen -WE-1718-
Yoni Verhaegen -WE-1718- on 6 May 2022
Thank you for your input, it really helps me a lot and I appreciate it. So if the velocity field were to change over time, mny original implementation (as stated in the code of my question here), would be correct?

Sign in to comment.

Accepted Answer

Torsten on 6 May 2022
Edited: Torsten on 6 May 2022
Formally, your implementation is correct as it is (I did not check the one-sided difference formulae in detail).
The problems I see are:
  1. the stability of your scheme (usually, 2nd order schemes which are not stabilized (as your scheme) tend to become instable/oscillatory). That's why you should start with first-order upwind or better: with my suggestion from below.
  2. the velocity field. Usually, a velocity field is divergence free to ensure mass-conservation (i.e. the condition du/dx + dv/dy = 0 automatically holds). In this case, your terms h*du/dx + h*dv/dy automatically cancel out and you are left with the equation dh/dt + u*dh/dx + v*dh/dy = 0. You should test this before starting. If you don't get du/dx + dv/dy approximately 0, you must be very careful with the results you get from your simulation.
  3. The boundary condition and where you have to set it.
  4. The implementation of the scheme. Note that one cell apart from an inflow boundary, you must use a value for h beyond the boundary points (h_debris(i+2,j),h_debris(i-2,j),(-h_debris(i,j+2),h_debris(i,j-2)). So, an extrapolation is necessary.
I strongly recommend using clawpack for your task where a second-order upwind implementation for the 2-dimensional advection equation is already implemented:
Yoni Verhaegen -WE-1718-
Yoni Verhaegen -WE-1718- on 6 May 2022
Unfortunately, I am not an expert in this field, so I don't really know how I neede to check and handle these issues. I put a lot of effort in it already to get it where it is now, and for that I am really happy that you have been helping me. Would it be an idea to share the model code somehow so that we can figure it out together?

Sign in to comment.

More Answers (0)


Find more on MATLAB Compiler 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!

Translated by