I am trying to compute the numerical solution of the 2D non linear convection equation
$$\frac{\partial u}{\partial t}+u\nabla u=0 $$
I am using the technique used in this work at page 38-40 .
The case that under consideration is a square domain of size $2\times2$ with Dirichlet boundary condition (BC) at north and south boundaries and periodic BC for the east and west boundaries .
The initial condition is a gaussian function $u_i=exp((x-0.8)^2+(y-1)^2)$ for the two component of $u$, as shown in figure (the quantity plotted is the velocity module $|u|$)
The problem is that during the computation there is an accumulation of the velocity magnitude around the maximum value of the velocity until there is a explosion of the simulation.
The two following picture shows the value of u at two times during the simulation. The first at $t=0.015$ ad the second $t=0.02$
This behaviour doesn't seem very physical to me. On the internet I am reading of the possibility of this kind of simulation of showing shock-like solution. The question is, is this the case ? Or it may depends on some bugs in the code ? Or a first order scheme is not suited for this kind of nonlinear equations ?
The only thing that I can say is, that is not linked to the CFL condition I have tried different to reduce the time step and the results are always the same.
Thank you for your help