# 2D Riemann Problem: Dimensional Splitting with Front Tracking

In this example, we consider the 2D scalar conservation law

with Riemann initial data, that is, with an intial function that has a constant value in each of the four quadrants and discontinuities along the x- and y-axes.

To solve the one-dimensional problems arising as part of the operator splitting, we will use front tracking followed by a projection step onto underlying regular grid.

## Initial setup

N = 100; h = 2/N; x = -1:h:1; y = 0.5*(x(1:end-1)+x(2:end)); T = 1.0; delta = 0.1/sqrt(N); [X,Y]=meshgrid(y,y);

## Exact solutions

We consider two different setups. In both setups u0=1/4 in the second and fourth quadrant. Setting u0=1/2 in the first quadrant and u0=-1/2 in the third quadrant gives a self-similar wave pattern where the four original constant states are separated by two pairs of rarefaction waves. For each pair, the two rarefaction waves meet in a sharp kink along the line y = x. Switching the values in the first and third quadrant gives a wave pattern that consists of the four constant states separated by six shocks forming two triple points at (3t/8,−t/8) and (−t/8, 3t/8).

subplot(2,2,1), pcolor(X,Y,rpsola(X,Y,0)), shading flat; colorbar, title('Problem A') subplot(2,2,2), pcolor(X,Y,rpsolb(X,Y,0)), shading flat; colorbar, title('Problem B') subplot(2,2,3), pcolor(X,Y,rpsola(X,Y,1)), shading flat; colorbar subplot(2,2,4), pcolor(X,Y,rpsolb(X,Y,1)), shading flat; colorbar

## Riemann problem B: shocks

First, we study the error for the second Riemann problem as a function of the number of operator splitting steps

u0 = rpsolb(X,Y,0); utrue = rpsolb(X,Y,1); rel = sum(sum(abs(utrue))); for i=1:4, Nstep = 4^(i-1); u = DimSplit(u0,x,x,'burgerRsol','burgerRsol',delta,Nstep,T,'wbar'); feil = sum(sum(abs(utrue-u)))/rel; subplot(2,2,i); pcolor(y,y,u), axis equal image; colorbar, shading flat title(sprintf('Grid: %dx%d Steps: %d', N, N, Nstep)); xlabel(sprintf('L1 error: %4.2f %%', 100*feil)); p=get(gca,'position'); p([3 4])=p([3 4])+0.03; set(gca,'position',p,'XTickLabel',[]); end;

From the figure, we see that the error is a convex function of the number of splitting steps. With few splitting steps, the error is dominated by the splitting error and the strong self-sharpening effects present near the shocks will counteract the smearing introduced by the projection operator that brings the front tracking solution back from a possibly irregular grid and back to the underlying Cartesian grid. As the number of time steps increases, the number of projections increases and the self-sharpening effects of diminish so that the smearing errors become dominant. Because the two error mechanisms work in different directions, a minimum error is observed for an intermediate number of splitting steps.

## Convergence test

We perform a grid-refinement study for the two Riemann problems, fixing the CFL number to 10 and report the error as a logarithmic plot

nu=10; errA=ones(1,6); errB=errA; for i=1:6, N = 8*2^i; h = 2/N; x = -1:h:1; y = 0.5*(x(1:end-1)+x(2:end)); delta = 0.1/sqrt(N); [X,Y] = meshgrid(y,y); % Riemann problem A u = DimSplit(rpsola(X,Y,0), x,x,'burgerRsol','burgerRsol', ... delta, round(T/(nu*h)), T,'wbar'); utrue = rpsola(X,Y,1); rel = sum(sum(abs(utrue))); errA(i) = sum(sum(abs(utrue-u)))/rel; % Riemann problem B u = DimSplit(rpsolb(X,Y,0), x,x,'burgerRsol','burgerRsol', ... delta, round(T/(nu*h)), T,'wbar'); utrue = rpsolb(X,Y,1); rel = sum(sum(abs(utrue))); errB(i) = sum(sum(abs(utrue-u)))/rel; end clf, semilogy(4:9,errA,'-o',4:9,errB,'x-') xlabel('log2 of grid size'); ylabel('Relative L1 error'); legend('Problem A','Problem B');

As pointed out above, the strong self-sharpening effects present near the shocks in Riemann problem B will counteract the smearing introduced by the projection operator leading to a convergence rate of approximately one. For the smooth solution resulting from Riemann problem A, the self-sharpening effects are much weaker, thereby giving a lower convergence rate.