tvd_rk2.m

ZEUS is a family of Eulerian (grid based) Magneto-Hydrodynamic codes (MHD) for use in astrophysics,

function tvd_rk2(alldata)

% This is the Navier-Stokes solver called by the GUI "Navier2d.m"
%
% The unsteady 2D incompressible Navier-Stokes equations are integrated
% using a 2nd order finite-volume (FV) method for the spatial discretisation 
% and a 1st/2nd order projection method for the time-stepping.
%
% The FV method is based on a non-staggered, vertex centred arrangement
% using the median dual mesh as the control volumes. 2nd order
% discretisations are developed for the divergence/gradient, Laplacian and
% non-linear operators:
%
%   - Divergence/gradient:  Green's theorem
%   - Laplacian          :  Green's theorem applied to face averaged
%                           gradients (unweighted least squares)
%   - Non-linear         :  Upwind piece-wise linear extrapolation with
%                           the HQUICK slope limiter
%
% The non-staggered arrangement was shown to support the common mesh scale
% pressure oscillation and hence the standard 2nd order incremental
% pressure correction method was modified to incorporate Rhie-Chow
% stabilisation. This reduced the order of the time integration to 1st/2nd.
%
% Unlike many other Navier-Stokes solvers this function makes use of a
% complete LU factorisation for the Poisson pressure equation. Using the
% fill reducing UMFPACK algorithm this approach was shown to be far more
% efficient than standard iterative approaches, although this could well be
% limited to MATLAB implementations and non-moving meshes.
%
% The momentum equations are solved in a semi-implicit manner. The
% non-linear terms are evaluated explicitly via a 2nd order TVD
% Runge-Kutta method. The BiCGSTAB method is used to solve the linear
% systems that arise due to the implicit treatment of the viscous terms.
%
%
%   AUTHOR : Darren Engwirda (2005-2006)
%   VERSION: 2.3 (08/05/2006)
%
% If you would like more information please contact me at:
%
%   [e-mail]
%
% With many useful contributions from Steve Armfield at The University of
% Sydney, Australia.
%
% Naver2d is Copyright (C) 2005-2006 Darren Engwirda. See "copyright.m" for
% full details.
%
%
%   UPDATE 19/04/2006: A HEAP OF CHANGES!! 
%
% (1): The order of the stages in the TVD RK update was switched to put the
% implicit step last (with the pressure correction in the middle). This
% should improve the quality of Neumann BC's for the velocities (outflow,
% symmetry etc).
%
% (2): Inverse distance weighted least squares are now used in the Laplacian.
%
% (3): An advection/diffusion equation for a scalar tracer has been added.
% There is also the option to solve coupled thermal problems where the
% tracer is taken as the temperature. Thermal coupling is done via the
% Oberbeck-Boussinesq quasi-incompressible method.
%
% (4): A module has been added to calculate the [x,y] forces on boundaries
% within the flow, based on the pressure and skin friction effects.
%
% (5): A new pressure outflow boundary condition has been added. The values
% at boundary pressures nodes are extrapolated from the neighbours (setting 
% dP/dn=0). This hopefully should lead to a much less reflective boundary
% condition, although I have noticed some problems when inconsistent IC's
% are used (see readme).
%
% (6): General code structure changes. It's a bit faster and hopefully
% easier to look at??
%
%
%   UPDATE 22/04/2006: BETTER MEMORY OVERHEAD
%
% Based on some failed attempts to use really big meshes (200,000 tri), the
% memory usage has been improved. The most significant change is the use of
% only one Laplacian operator (as opposed to 4 previously). This operator
% is built with Neumann BC's. Dirichlet BC's are imposed approximately
% by adding large terms in the diagonal mass matricies. I have also tried
% to remove a heap of intermediate variables. Anyway, the memory usage is
% now about 1/2 that of previous versions. The maximum mesh size possible is
% computer dependent, although current versions of windows/MATLAB impose a 
% 1.2GB PF limit, independent of the amount of RAM you have...
%
% Big problems often can't fit in physical RAM and make use of virtual
% memory. This is bad. Generally I've seen that the CPU loading decreases
% dramatically (to like 5%) as soon as the memory use spills over into
% virtual. I'm currently trying to think about ways to fix this, but a 
% couple of obvious ones:
%   - get more RAM
%   - don't run any other programs in the background
%   - clear the MATLAB workspace before running navier2d
%
% One possible fix is to switch over to a fully iterative solver for the
% Poisson problem (like an SSOR preconditioned BiCGSTAB). For problems 
% that can use the RAM, I have seen that this can be about 3 times slower 
% than the UMFPACK LU method, BUT, because the LU information would no 
% longer need storage it has the potential to allow bigger meshes to fit 
% in the RAM, giving 100% CPU loads. I could switch over to the fully
% iterative solver once a certain size has been reached...
%
% At the moment, I've decided not to do this, because most computers these
% days have >=1GB of RAM, so MATLAB should run out of PF first. This is not
% always going to be the case (espec for 64bit systems) so if you think you
% are having memory problems, PLEASE CONTACT ME!
%
%
%   UPDATE 02/05/2006: NEW LAPLACIAN (FASTER!!)
%
% (1): A new discretisation of the Laplacian has been implemented. The new method 
% is very similar to a Galerkin FEM, although I (derived) and implemented
% it based on the FVM. The new method uses a much more compact stencil,
% only referencing the nearest-neighbours (as opposed to the next-to-nearest
% like the old method) so the resulting operator is much more sparse. This
% is good because:
%   - Less memory is used to store both the matrix and the LU information.
%   - SMVP's are completed faster.
%   - The LU forward/back solve is completed faster.
%   - The LU decomposition via UMFPACK seems to be stable for lower pivot
%     tolerances (gives sparser LU).
%
% The numerical results seem to be very similar compared to results obtained 
% using the old Laplacian, and based on a mesh refinement study of the steady 
% driven cavity flow the method shows very similar approx 2nd order convergence.
% Overall, the moral of this story is that (basically) identical results
% should now be obtained in about 75% of the CPU time taken previously and
% bigger meshes can be used before running out of RAM/memory.
%
% (2): Better boundary conditions. Extrapolated outflow BC's now use a piecewise
% linear extrapolation for all flow variables. This makes things much smoother,
% i.e the vorticity is now smoothly continuous through the boundary...


% Call the MATLAB memory garbage collection
pack


% Extract data
solver    = alldata.solver;         % Solver type: UVP, tracer, thermal
animation = alldata.animation;      % Animation settings
flag      = alldata.flag;           % Calculate lift/drag forces for edges with flag=true

% Mesh based data
p    = alldata.mesh.p;              % Nodes
t    = alldata.mesh.t;              % Triangulation
bnd  = alldata.mesh.bnd;            % True for bnd nodes
A    = alldata.mesh.A;              % Cell area
e    = alldata.mesh.e;              % Edges
pc   = alldata.mesh.pc;             % Centroids 
e2t  = alldata.mesh.e2t;            % Edge to triangle connectivity
n2n  = alldata.mesh.n2n;            % Node to node connectivity
hnx  = alldata.mesh.hnx;
hny  = alldata.mesh.hny;            
numn = size(p,1);                   % Number of nodes
nume = size(e,1);                   % Number of edges


% Sparse finite volume operators
[L,Gx,Gy,Lp,Up,pp,cp,bc] = build_operators(alldata.mesh,alldata.bc);


% ===================================================================
%                      INTEGRATION SETTINGS
% ===================================================================

nu    = alldata.settings(4);        % Viscosity
mu    = alldata.settings(6);        % Viscosity (tracer)
mmax  = alldata.settings(1);        % Max steps
tmax  = alldata.settings(2);        % Max time
CFL   = alldata.settings(3);        % CFL number
freq  = alldata.settings(5);        % Output frequency
cgtol = 1e-5;                       % Tolerance for BiCGSTAB solver
maxit = 25;                         % Max iters for BiCGSTAB solver  

if solver==3
    alpha = alldata.settings(7);    % Temperature/velocity coupling
else
    alpha = 0;                      % Tracer only - no coupling
end

% Flag BC's 
dbcu     = find(bc(:,1)==1);
dbcv     = find(bc(:,3)==1);
dbcp     = find(bc(:,5) >0);  
dbcs     = find(bc(:,7)==1);  
extrap_p = find(bc(:,5)==2);
extrap_u = extrap_p;                % At this stage only have "extrap"
extrap_v = extrap_p;                % BC's at outflows... I may extend
extrap_s = extrap_p;                % this in the future.

% Initial conditions
U = alldata.init.U;  U(dbcu) = bc(dbcu,2);
V = alldata.init.V;  V(dbcv) = bc(dbcv,4);
S = alldata.init.S;  S(dbcs) = bc(dbcs,8);
P = alldata.init.P;

% Viscous coefficients
diagL = full(diag(L));      % Diagonal terms

% Mass matrcies (use diagonal scaling to get Dirichlet BC's)
M_u = repmat(1,numn,1);  M_u = 2*M_u.*A/nu;  M_u(dbcu) = 1e3 * diagL(dbcu);  M_u(bc(:,1)==0) = 0; 
M_v = repmat(1,numn,1);  M_v = 2*M_v.*A/nu;  M_v(dbcv) = 1e3 * diagL(dbcv);  M_v(bc(:,3)==0) = 0;
M_s = repmat(1,numn,1);  M_s = 2*M_s.*A/mu;  M_s(dbcs) = 1e3 * diagL(dbcs);  M_s(bc(:,7)==0) = 0;


% ===================================================================
%                      NON-LINEAR COEFFICIENTS
% =================================================================== 

% Internal triangles
in = e2t(:,2)>0;

% Edge midpoints
xm = 0.5*(p(e(:,1),1)+p(e(:,2),1));
ym = 0.5*(p(e(:,1),2)+p(e(:,2),2));

% Median edge midpoints
xm1 = 0*xm;  xm1(in) = 0.5*(pc(e2t(in,1),1)+xm(in));  
xm2 = 0*xm;  xm2(in) = 0.5*(pc(e2t(in,2),1)+xm(in));  
ym1 = 0*ym;  ym1(in) = 0.5*(pc(e2t(in,1),2)+ym(in));
ym2 = 0*ym;  ym2(in) = 0.5*(pc(e2t(in,2),2)+ym(in));

% xy from nodes to median edge midpoints
dx = [xm1-p(e(:,1),1), xm1-p(e(:,2),1), xm2-p(e(:,1),1), xm2-p(e(:,2),1)];  
dy = [ym1-p(e(:,1),2), ym1-p(e(:,2),2), ym2-p(e(:,1),2), ym2-p(e(:,2),2)];

% 3rd node in associated triangles
% for each edge
N3 = repmat(0,nume,2);
for k = 1:nume
    if e2t(k,2)>0
        for q = 1:3
            if (t(e2t(k,1),q)~=e(k,1)) && (t(e2t(k,1),q)~=e(k,2))
                N3(k,1) = t(e2t(k,1),q);
            end
            if (t(e2t(k,2),q)~=e(k,1)) && (t(e2t(k,2),q)~=e(k,2))
                N3(k,2) = t(e2t(k,2),q);
            end
        end
    end
end

% COEFFCIENTS FOR LIFT/DRAG CALCS
if any(flag)
    he  = sqrt(hnx(flag,2).^2+hny(flag,2).^2);
    nxe = hnx(flag,2)./he;
    nye = hny(flag,2)./he;
end


% ===================================================================
%                   SETUP ANIMATION FIGURE WINDOWS
% ===================================================================

set(figure(2),'Name','Residuals','Units','Normalized','Position',[0.05 ,0.05 ,0.45,0.35]);
set(figure(3),'Name','UV Quiver','Units','Normalized','Position',[0.525,0.525,0.45,0.35]);
if animation(1)>0
    switch animation(1)
        case 1, set(figure(4),'Name','U Velocity'    ,'Units','Normalized','Position',[0.05,0.525,0.45,0.35]);
        case 2, set(figure(4),'Name','V Velocity'    ,'Units','Normalized','Position',[0.05,0.525,0.45,0.35]);
        case 3, set(figure(4),'Name','Total Velocity','Units','Normalized','Position',[0.05,0.525,0.45,0.35]);
        case 4, set(figure(4),'Name','Pressure'      ,'Units','Normalized','Position',[0.05,0.525,0.45,0.35]);
        case 5, set(figure(4),'Name','abs(Vorticity)','Units','Normalized','Position',[0.05,0.525,0.45,0.35]);
        case 6, set(figure(4),'Name','Tracer'        ,'Units','Normalized','Position',[0.05,0.525,0.45,0.35]);
    end
end
if animation(2)>0
    switch animation(2)
        case 1, set(figure(5),'Name','U Velocity'    ,'Units','Normalized','Position',[0.525,0.05,0.45,0.35]);
        case 2, set(figure(5),'Name','V Velocity'    ,'Units','Normalized','Position',[0.525,0.05,0.45,0.35]);
        case 3, set(figure(5),'Name','Total Velocity','Units','Normalized','Position',[0.525,0.05,0.45,0.35]);
        case 4, set(figure(5),'Name','Pressure'      ,'Units','Normalized','Position',[0.525,0.05,0.45,0.35]);
        case 5, set(figure(5),'Name','abs(Vorticity)','Units','Normalized','Position',[0.525,0.05,0.45,0.35]);
        case 6, set(figure(5),'Name','Tracer'        ,'Units','Normalized','Position',[0.525,0.05,0.45,0.35]);
    end
end
drawnow


% ===================================================================
%                          MAIN LOOP
% ===================================================================

% Memory garbage collection
clear i j pc xm ym xm1 ym1 xm2 ym2 he, pack

% Loop counters
tN     = 0; 
tV     = 0; 
tQ     = 0; 
tstart = cputime;
m      = 0;

% Initial grad(P)
Px = smvp(Gx,P);
Py = smvp(Gy,P);

% Residuals
resx     = repmat(0,mmax,1);
resy     = repmat(0,mmax,1);
time     = repmat(0,mmax,1);
residual = ~any(flag);

% Viscous time-step limit
if solver==1
    dtvisc = 2*min(A)/nu;
else
    dtvisc = 2*min(A)/max(nu,mu);    
end

% Estimate initial dt from BC's
sqrtA  = sqrt(A);
idt    = max(sqrt(U.^2+V.^2)./sqrtA)/CFL;
dt     = min(1/max(idt,eps), dtvisc);


Snew = S;
while (time(m+1)<tmax)&&(m<mmax)
    
    
    %     A NOTE ON THE HQUICK SCHEME: The non-linear eval is based 
    %     on a piecewise-linear upwind reconstruction at each 
    %     median edge midpoint. A "limiter" (HQUICK) is applied to the 
    %     reconstruction to try to make the scheme "positive" 
    %     (non-oscillatory) in regions of sharp gradient while 
    %     non-diffusive in smooth regions. The current implementation 
    %     is a short cut, and does not apply the limiter in a proper 
    %     multidimensional fashion. This means that the scheme is not 
    %     truly TVD and can admit oscillations. The scheme should be 
    %     bounded (TVB??).
    
    
    % TVD RK3 METHOD - 1ST RK STAGE:
 
    
    % ===================================================================
    %                       TRANSPORT EQUATIONS
    % ===================================================================
    tic
    
    if solver>1     %   Solve advection-diffusion equations 
                    %   for the [U,V] velocity components as 
                    %   well as the scalar tracer S.
        
                
        % NON-LINEAR EVAL (Limited piecewise-linear upwind)
        
        % smvp's
        ux = smvp(Gx,U);  uy = smvp(Gy,U);
        vx = smvp(Gx,V);  vy = smvp(Gy,V);
        sx = smvp(Gx,S);  sy = smvp(Gy,S);
        
        % Assemble non-linear vectors
        [Nu,Nv,Ns] = nonlinear_uvs(U,ux,uy,V,vx,vy,S,sx,sy,e,N3,in,hnx,hny,dx,dy);
        
        % Explicit Euler stage 1
        Ss      = S + dt*( (mu*smvp(L,S) - Ns)./A                );  
        Us      = U + dt*( (nu*smvp(L,U) - Nu)./A - Px           );  
        Vs      = V + dt*( (nu*smvp(L,V) - Nv)./A - Py + alpha*S );  
        Ss(bnd) = S(bnd);
        Us(bnd) = U(bnd);
        Vs(bnd) = V(bnd);
        
        % Extrapolate boundary values where necessary
        if any(extrap_u), Us = bnd_extrap(Us,ux,uy,extrap_u,p,n2n); end
        if any(extrap_v), Vs = bnd_extrap(Vs,vx,vy,extrap_v,p,n2n); end
        if any(extrap_s), Ss = bnd_extrap(Ss,sx,sy,extrap_s,p,n2n); end
        
        
    else    %   Only solve advection-diffusion equations 
            %   for the [U,V] velocity components.

        
        % NON-LINEAR EVAL (Limited piecewise-linear upwind)
        
        % smvp's
        ux = smvp(Gx,U);  uy = smvp(Gy,U);
        vx = smvp(Gx,V);  vy = smvp(Gy,V);
        
        % Assemble non-linear vectors
        [Nu,Nv] = nonlinear_uv(U,ux,uy,V,vx,vy,e,N3,in,hnx,hny,dx,dy);
        
        % Explicit Euler stage 1
        Us      = U + dt*( (nu*smvp(L,U) - Nu)./A - Px );  
        Vs      = V + dt*( (nu*smvp(L,V) - Nv)./A - Py ); 
        Us(bnd) = U(bnd);
        Vs(bnd) = V(bnd);