% Testing the ice equations in Lagrangian form.
% The equations are u_t = -p_x; rho_t = u_x; with the constraints rho>=0, p>=0.
% The domain is periodic on [0,1].
% Initial conditions are rho constant at 0.5 and u being a sin.

% We are here minimizing the one-norm of p at each time step, subject to the 
% above constraints, using a simplex method.

% The inputs are spatial discretization dx, ratio of dt/dx=mu, and max time.

% The outputs are matrices of velocites u, pressures p, and densities rho at 
% each time and spatial step.

% Complementary program needed is simplex.m, which in turn requires pivot.m


% Written by helga schaffrin, oct. 2003


function [u,p,rho,x,x1,t]=ice_simplex(dx,mu,T);

% First we need to set up the problem

dt=mu*dx;               % Temporal step size

J=floor(1/dx);          % Number of spatial steps
N=floor(T/dt);          % Number of temporal steps

x=[0:dx:1-dx];          % Spatial grid for velocity
x1=x+dx/2;              % Spatial grid for rho and p
t=[0:dt:N*dt];          % Temporal grid


% Initialization

u=zeros(N+1,J);
rho=zeros(N+1,J);
p=zeros(N+1,J);


% Initialization for items used in the simplex minimization

v=zeros(1,2*J);         % Placeholder for (p,q), whereby q has J elements and contains
                        % the slack variables needed to put the problem
                        % into standard simplex form; q corresponds to the updated
                        % values of rho.
p_norm=0;               % Placeholder which temporarily records the 1-norm of p

A=zeros(J,2*J);         % Matrix and vector expressing the rho>=0 constraints in
b=zeros(J,1);           % standard simplex form Av=b
c=zeros(1,2*J);         % To express 1-norm(p) as c*v

% Defining A
% The left half is a tridiagonal submatrix with two entries in the other
% corners; the right half (corresponding to the slack variables) is an 
% identity submatrix.

A(:,J+1:2*J) = eye(J);
A(:,1:J) = -2*mu^2*eye(J);
A(1:J-1,2:J) = A(1:J-1,2:J) + mu^2*eye(J-1);
A(2:J,1:J-1) = A(2:J,1:J-1) + mu^2*eye(J-1);
A(J,1) = mu^2;
A(1,J) = mu^2;

% Defining c

c(1:J)=1;


% Initial conditions

rho(1,:) = 0.5;
u(1,:) = sin(2*pi*x);


% Time stepping
% This happens in two steps: (1) Using the expression for the advanced
% rho without constraints, we find p to enforce the constraint, keeping the
% norm of p at a minimum (simplex method).  As a by-product, we also get the
% advanced values of rho.  (2) Substituting the values found for p, we find u 
% at the new time level.
% And the game starts over...

for n=1:N

    % Step (1)
    % In order to minimize norm(p), we will use the simplex program, which
    % requires as input the c, d, A, and b from the problem in standard
    % simplex form (i.e. z=c*v+d, constraints: Av=b, v>=0).  d=0; A, c are
    % defined above; so it remains to find b, which changes with time.
    
    b(1:J-1) = u(n,2:J)' - u(n,1:J-1)';
    b(J) = u(n,1) - u(n,J);
    b = mu*b + rho(n,:)';
    
    [v,p_norm] = simplex(c, 0, A, b);
    
    % The first J entries of x are p; the others are rho(n+1,:).
    
    p(n+1,:) = v(1:J);
    rho(n+1,:) = v(J+1:2*J);
    
    % Step (2)
    
    u(n+1,2:J) = u(n,2:J) - mu*(p(n+1,2:J) - p(n+1,1:J-1));
    u(n+1,1) = u(n,1) - mu*(p(n+1,1) - p(n+1,J));
    n
end