Shallow Water DAE - Index may be greater than 1 ?

Dear all,

I am trying to solve a DAE for the 1D nonlinear shallow water (SW) system of equations

where are the velocity and height field respectively.

The intial conditions are given by and in the spatial domain

. The boundary conditions are and .

The solution of this dam-break problem will produce a left-going rarefaction wave and a right-propagating shock wave.

In order to better capture the shock wave I want to solve simultaneously the SW system with a moving mesh PDE of equation , where ξ is a uniform computational mesh in [0,1] and is a density function.

By rewriting as function of I obtain the following DAE system of index 1

Now, I have discretized the system in space by finite difference method with the moving mesh in the space domain and fixed mesh in the computational domain for

I have in total a system of equations of this form

where and .

I am trying to use the MATLAB ode solver ode23t/ode15s for the system with a singular diagonal Mass matrix ( 0-th in the last N entries). Initially, the physical mesh is supposed to be uniform and v is consistent with the algebraic constraint.

After running the script I obtain the following error :

Error using daeic12 (line 76)
This DAE appears to be of index greater than 1.

Do you have idea why MATLAB does not recognize the system as index 1 ? Is it then more reasonable to write the system without the auxiliary variable v and changing the entries of the mass matrix (no longer diagonal) ?

I have also tried to solve the problem with an implicit solver ode15i by discarding the auxiliary variable v but I get the same error message when I call the decic operator for computing consistent initial conditions.

Thank you for your help.

[Bill Greene]

The error message says "appears to be" but really should say "is." I don't understand your mesh moving function so can't suggest what the problem might be. What is τ? What is Δ?

ode15s can sometimes solve DAE with index>1 if you provide a set of "consistent" initial conditions. That is, your y0 must satisfy your equations at . If it does, then the routine where this error message appears is bypassed.

[Sa_math]

Do you mean that the index is greater than 1 because I have to differentiate the N components of v to obtain a full system of ode ?

I tried to compute consistent initial conditions using fsolve but the solver stopped prematurely

for j = 2:(nx - 1)  

    v(j) = (0.5*(nx-1)^2/(tau*rho(j)))*((x(j+1)-x(j))*(rho(j+1)+rho(j)) ...
             -(x(j)-x(j-1))*(rho(j)+rho(j-1)));
         
end

% Starting guess for physical variables
y0 = [reshape(s',2*nx,1);x;v];


y0 = fsolve(@(v)  fSWE(0,v,nx),y0);

Solver stopped prematurely.

fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 40000 (the default value).

I think that the problem is that I also don't know the value of u'(0) and h'(0). Is there another routine to bypass this errror ?

Thanks.

[Bill Greene]

To be index-1, , must be non-singular. That is specifically what daeic12 is testing for when it issues the error message. The problem is at the end points where your code is

% endpoints are fixed
  yp(3*nx+1) =  0.0;
  yp(nn) = 0.0;  

If you rewrite your code to compute for ALL points in the mesh and remove these so-called "boundary conditions", you can resolve the cause of this error message.

[Sa_math]

Thank you for clarifying this point. Removing the conditions at the boundary I don't bump into any error. Though, the ode23t/15s sub-routine is selecting a really tiny time step in order to meet the integration tolerance, even if this is set of order e-2 in odeset.

Is there a free reference that explains how this method work? I would like to understant if this issue might be due to the Finite Difference Scheme for the SWEs. In that case, I will attempt a different method (FVM, FEM) for solving the system.

[Bill Greene]

Links to downloadable copies of the references for ode15s are at the bottom of the documentation page: ode15s.

If you post your updated m-code, I can run that and perhaps comment further.

If you describe your mesh-moving function in more detail, I might be able to say more about your approach.

I don't think that simply switching to a FVM or FEM will, by itself, improve the behavior of the solution. Looking at your FDM equations, it looks like you are using a central difference approximation to the spatial deriviatives. This is known to create spurious oscillations unless you have a very fine mesh. If you use a backward difference approximation instead (known as "upwind" in the literature) you will eliminate the oscillations but at a cost of excessive smoothing of the solution.

[Sa_math]

Thank you for the reference, Bill. You can find the updated code in attachment.

I am using the moving mesh PDE because I want to reduce the mesh size in proximity of the shock front and consequently the error by linear interpolation. This is realized by the arc-length monitor function in the definition of the PDE. The user-specified parameter τ tweaks the level of adaptability of the physical mesh to the SW system: The smaller τ, the faster the mesh will adapt to the solution.

The upwind scheme does not produce better results than the centered finite difference scheme. I just notice that at the first time step the velocity field starts increasing on the left boundary of the domain and it propagates rightwards before resetting again to 0. I should expect instead that the velocity field starts raising at the center of my domain, where the dam breaks. The same behaviour is obtained for ode15i, so I suppose that there might be instabilities due to the FDM and a FVM might resolve the problem.

Any further help is highly appreciated. Thanks.