CasADi Integrator setup for transport equation

Hello there,

I am using Matlab R2023a and the CasADi (https://web.casadi.org/) 3.6.4 version for Windows.

What I want to do with it is to solve a transport/traffic equation with the algorithm of CasADi of the following form:

ρ: density

: constant velocity

h: spatial step size

So you can see that this is an ODE with in which I already discretized the spatial differentiation. For describing my problem I set the initial condition of rho to a gaussian and the boundary condition on the "left side" to stay constant. I setup a grid with 40 grid points and want to integrate over one time step. In the following you see the code that I used until now:

addpath(genpath("casadi-3.6.4-windows64-matlab2018b"));

import casadi.*

%% Integrator creation

% symbolic variables for setting up the integrator

rho = MX.sym('rho',[params.sim.gridPoints,1]);

dxi = MX.sym('dxi');

u_dot = MX.sym('u_dot');

v_const = MX.sym('v_const');

% model equations

x =[rho]; % state

p =[u_dot;v_const;dxi] % parameters

diffEquation =[u_dot; % first grid point bc of boundary condition

-v_const*(rho(2:end)-rho(1:end-1))/dxi; % all other points

];

% summing up the dynamics

mySys = struct('x',x,'p',p,'ode',diffEquation);

% create integrator

TranspIntegrator=integrator('TranspIntegrator','collocation',mySys);

disp(TranspIntegrator)

%% Problem Solution

u_dot = 0; % this is the left boundary condition. With u_dot=0 the most left rho value should stay the same.

v_const = 100; % chosen constant velocity

dxi = 1; % spatial step size

% create initial condition

xAxes = linspace(0,4,40); % spatial vector for creating initial condition

rho_0 = (10+(250-40)*gaussmf(xAxes,[0.5 2]))'; % gaussian initial condition

disp(x_0)

r = TranspIntegrator('x0',rho_0,'p',[u_dot;v_const;dxi]);

disp(r.xf)

The displayed result looks like this:

[10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704, 10.0704]

Regarding the fact that I have a constant velocity in a transport equation I would expect my rho values to shift along the x Axes in positiv direction. What the results are giving me are 40 values that are in direct realtion to rho_0(1) and u_dot (when I changed this value there where slight changes in the result).

I guess I missed something in the implementation of the solving algorithm but I could not find it. So I am asking for your help. If I forgot some crutial information I will deliver it afterwards. Anyway I already want to thank you for any help. Best Regards!

Answer 1

[Torsten]

With a velocity of 100 and a length of 4 for the region of integration, it takes at most 4/100 sec until "rho" gets constant and equal to its value at the inlet. What is the time for which you received the result vector ?

[Cedric]

Hey Torsten,

Thanks a lot for your help, I think this was a typical case of tunnel vision :D

As you already guessed right I recived the result vector for an iteration over one second.

I changed the v_const to 4 and plotted my result for the first 3 time steps. These results also fit into my expectation. If anyone is interested the picture of the plots is in the appendix.

So thanks again! :)

[Torsten]

But if v = 4 and t = 1, why does the curve only advance by dx = 0.5 in positive x-direction ? Or does t = 1 mean t = 0.125 ?

[Cedric]

I think that I am still having some issues with the units.

In the upper example I had dxi=1. And chose a totally different space in between my xAxes for my initial conditions what does not make sense.

dxi = 1; % spatial step size

% create initial condition

xAxes = linspace(0,4,40); % spatial vector for creating initial condition

I changed these conditions again to the following ones:

%% solve problem

u_dot = 0; % this is the left boundary condition. With u_dot=0 the most left rho value should stay the same.

v_const = 1;

length = 4;

gridPoints = 40;

dxi = length/gridPoints;

% create initial condition

xAxes = linspace(0,length,gridPoints); % spatial vector

As you can see I also changed v_const to 1.

Plotting this leads to the results in the appendix (I shifted the initial conditions to make the effect of the changes clear). I think now it makes more sense.

[Torsten]

Yes, the speed of transport looks correct now. But the dissipation is far too strong. The curve in its form should remain - it should only be translated to the right. But the peak gets lower and the whole curve smears out. This is due to your first-order discretization with Euler backwards in space. You will either need many more grid points or an upwind discretization of second order to keep your curve in its original form.

[Cedric]

I see your point and will consider these steps in my ongoing work. Thanks for your great help :)

If there will come up some more interesting things up I will post theme here.