%**************************************************************************
%*******                Simulation for Hookean dumbbell        ************
%*                         shear flow  Oldroyd B                          *
%*                                                                        *
%**************************** Tony Lelievre 30/05/2007 ********************
%**************************************************************************

% run Couette_Oldroyd_B

% The nondimensional equations are:
%
% Re . d_t u = (1- Epsilon) d_x,x u + (Epsilon/We) d_x tau
% u(0,x)=0
% u(t,0)=v
% u(t,1)=0
% d_t tau + (1/We) tau = d_x u
%
% The velocity on the boundary is actually progressively set to v
%
% Warning: Here, tau=\E(PQ) (and not (Epsilon/We)\E(PQ))

% Oldroyd-B

clear all;

% Physical parameters
Re=0.1;
Eps=0.9;
We=0.5;
v=1.;
T=1.; % Maximal time

% Discretization
% Space
I=100;
dx=1/I;
mesh=[0:dx:1];
% Time
N=100;
dt=T/N;

% Matrices

D1=diag(ones(1,I-1),-1);D1=D1(2:I,:);D1=[D1,zeros(I-1,1)];
D2=diag(ones(1,I-1));D2=[zeros(I-1,1),D2,zeros(I-1,1)];
D3=diag(ones(1,I-1),+1);D3=D3(1:(I-1),:);D3=[zeros(I-1,1),D3];

% Mass matrix

M=(1/6)*D1+(2/3)*D2+(1/6)*D3;
M=M.*dx;
M=sparse(M);
MM=M(:,2:I);

% Stiffness matrix
  
A=(-1)*D1+2*D2+(-1)*D3;
A=A./dx;
A=sparse(A);
AA=A(:,2:I);

% Vectors

u=zeros(I+1,1); % Initial velocity
tau=zeros(I,1); % Initial stress: \E(PQ)=0 at t=0
gradtau=zeros(I-1,1);

% Time iterations

BB=Re*MM./dt+(1-Eps)*AA;
CLL=zeros(I+1,1);
for t=dt:dt:T,
  uold=u;
  gradtau=tau(2:I)-tau(1:(I-1));
  if ((t/T)<0.1) 
	CLL(1)=v*10*(t/T);
  else
	 CLL(1)=v ;
  end;
  CL=(Re*M./dt+(1-Eps)*A)*CLL;
  F=(Re*M./dt)*u-CL+(Eps/We)*gradtau;
  u(2:I)=BB\F;
  
  if ((t/T)<0.1) 
	u(1)=v*10*(t/T);
  else
	 u(1)=v;
  end;

  for l=1:I	
	tau(l)=(1-dt/We).*tau(l)+(dt/dx)*(u(l+1)-u(l));
%     tau(l)=(1-dt/We).*tau(l)+dt/dx*(uold(l+1)-uold(l));
  end;
  
  % Drawings
  plot(mesh',u,mesh',[(Eps/We)*tau;(Eps/We)*tau(I)]);
  axis([0 1 -1 1.2]);  
  drawnow;

end;

legend('velocity','stress');