tvd_rk2.m

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

        % 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
        
        
    end
    tN = tN + toc;
    
    
    % 2ND RK STAGE
    
    
    % ===================================================================
    %                       TRANSPORT EQUATIONS
    % ===================================================================
    
    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)
        tic
        
        % smvp's
        ux = smvp(Gx,Us);  uy = smvp(Gy,Us);
        vx = smvp(Gx,Vs);  vy = smvp(Gy,Vs);
        sx = smvp(Gx,Ss);  sy = smvp(Gy,Ss);
        
        % Assemble non-linear vectors
        [Nu,Nv,Ns] = nonlinear_uvs(Us,ux,uy,Vs,vx,vy,Ss,sx,sy,e,N3,in,hnx,hny,dx,dy);
        
        tN = tN + toc;
        
        
        % Implicit viscous step (TVD RK2 for non-linear terms)
        tic 
        
        % Mass matrices
        M_u_dt = M_u/dt;
        M_v_dt = M_v/dt;
        M_s_dt = M_s/dt;
        
        % RHS vectors (put BC's in rhs)
        rhss      = ( A.*( S+Ss)/dt                  - Ns )/mu;  
        rhsu      = ( A.*((U+Us)/dt - Px)            - Nu )/nu;   
        rhsv      = ( A.*((V+Vs)/dt - Py + alpha*Ss) - Nv )/nu;  
        rhss(bnd) = M_s_dt(bnd).*bc(bnd,8);
        rhsu(bnd) = M_u_dt(bnd).*bc(bnd,2);
        rhsv(bnd) = M_v_dt(bnd).*bc(bnd,4);
        
        % Solve linear systems
        Us         = myBiCGSTAB(L,M_u_dt,rhsu,Us,cgtol,maxit,diagL);  
        Vs         = myBiCGSTAB(L,M_v_dt,rhsv,Vs,cgtol,maxit,diagL);  
        Snew       = myBiCGSTAB(L,M_s_dt,rhss,Ss,cgtol,maxit,diagL);  
        Us(dbcu)   = U(dbcu);
        Vs(dbcv)   = V(dbcv);
        Snew(dbcs) = S(dbcs);
        
        % 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), Snew = bnd_extrap(Snew,sx,sy,extrap_s,p,n2n); end
        
        tV = tV + toc;
        
        
    else    %   Only solve advection-diffusion equations 
            %   for the [U,V] velocity components.

        
        % NON-LINEAR EVAL (Limited piecewise-linear upwind)
        tic
        
        % smvp's
        ux = smvp(Gx,Us);  uy = smvp(Gy,Us);
        vx = smvp(Gx,Vs);  vy = smvp(Gy,Vs);
        
        % Assemble non-linear vectors
        [Nu,Nv] = nonlinear_uv(Us,ux,uy,Vs,vx,vy,e,N3,in,hnx,hny,dx,dy);
        
        tN = tN + toc;
        
        
        % Implicit viscous step (TVD RK2 for non-linear terms)
        tic 
        
        % Mass matrices
        M_u_dt = M_u/dt;
        M_v_dt = M_v/dt;
        
        % RHS vectors (put BC's in rhs)
        rhsu      = ( A.*((U+Us)/dt - Px) - Nu )/nu;   
        rhsv      = ( A.*((V+Vs)/dt - Py) - Nv )/nu;  
        rhsu(bnd) = M_u_dt(bnd).*bc(bnd,2);
        rhsv(bnd) = M_v_dt(bnd).*bc(bnd,4); 
        
        % Solve linear systems
        Us       = myBiCGSTAB(L,M_u_dt,rhsu,Us,cgtol,maxit,diagL);  
        Vs       = myBiCGSTAB(L,M_v_dt,rhsv,Vs,cgtol,maxit,diagL);  
        Us(dbcu) = U(dbcu);
        Vs(dbcv) = V(dbcv);
        
        % 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
        
        tV = tV + toc;
        
    end  


    % ===================================================================
    %                     PRESSURE POISSON EQUATION
    % ===================================================================
    tic
    
    % Pressure "free" velocity field
    Uh       = Us + dt*Px;  
    Vh       = Vs + dt*Py;  
    Uh(dbcu) = Us(dbcu);
    Vh(dbcv) = Vs(dbcv);
    
    % Divergence
    D       = A.*(smvp(Gx,Uh) + smvp(Gy,Vh)) - dt*smvp(L,P);  
    D(dbcp) = 0;
    
    % Solve Poisson for pressure correction
    Q  = LUsubs(Lp,Up,pp,cp,D/dt); 
    Q  = Q-sum(Q)/numn;     % Set mean=0
    Qx = smvp(Gx,Q);
    Qy = smvp(Gy,Q);
 
    % ===================================================================
    %                        MOMENTUM CORRECTOR
    % ===================================================================
    
    Unew       = Us - dt*Qx;  
    Vnew       = Vs - dt*Qy;  
    Unew(dbcu) = Us(dbcu);
    Vnew(dbcv) = Vs(dbcv);
    
    % Pressure update
    P = P+Q;  Px = Px+Qx;  Py = Py+Qy;
    
    % Extrapolate boundary values where necessary
    if any(extrap_p)
        P  = bnd_extrap(P,Px,Py,extrap_p,p,n2n);
        Px = smvp(Gx,P);  
        Py = smvp(Gy,P);
    end
    
    
    % Extrapolate boundary values where necessary
    if any(extrap_u), Unew = bnd_extrap(Unew,ux,uy,extrap_u,p,n2n); end
    if any(extrap_v), Vnew = bnd_extrap(Vnew,vx,vy,extrap_v,p,n2n); end
    
    
    tQ = tQ + toc;
    
       
    % ===================================================================
    %                          FULL UPDATE
    % ===================================================================
    
    % Residuals
    if residual     % Velocity residuals
        
        resx(m+1) = sum(abs(Unew-U))/(dt*numn*max(norm(Unew,inf),eps));
        resy(m+1) = sum(abs(Vnew-V))/(dt*numn*max(norm(Vnew,inf),eps));
        
    else            % Lift/drag forces
        
        % Average velocity gradients on bnd edges
        uxe = 0.5*(ux(e(flag,1))+ux(e(flag,2))); uye = 0.5*(uy(e(flag,1))+uy(e(flag,2)));
        vxe = 0.5*(vx(e(flag,1))+vx(e(flag,2))); vye = 0.5*(vy(e(flag,1))+vy(e(flag,2)));
        
        % Form dvt/dn (normal gradient of tangential velocity) on bnd edges
        dvt_dn = nxe.*nye.*(uxe-vye) - (nxe.^2).*vxe + (nye.^2).*uye;
        
        % Average pressure on bnd edges
        Pe = 0.5*(P(e(flag,1))+P(e(flag,2)));
        
        % Drag/lift forces
        resx(m+1) = 2*sum( nu*hny(flag,2).*dvt_dn - hnx(flag,2).*Pe );          % Multiply by 2, remember that [hnx,hny]
        resy(m+1) = 2*sum( nu*hnx(flag,2).*dvt_dn + hny(flag,2).*Pe );          % on the boundary is *0.5...   
        
    end
    
    % Updates
    U = Unew;
    V = Vnew;
    S = Snew;
    
    % CFL based step-size
    idt = max(sqrt(U(~bnd).^2+V(~bnd).^2)./sqrtA(~bnd))/CFL;
    dt  = min(1/max(idt,eps), dtvisc);
    
    
    % ===================================================================
    %                           ANIMATION
    % ===================================================================

    if ~mod(m+1,freq)
        
        figure(2)       % Always plot residuals

        semilogy(time,abs(resx),'b',time,abs(resy),'r'), grid on 
        
        legend('x-residual','y-residual'); xlabel(['Time (sec) at step: ',num2str(m+1)])
        
        
        figure(3)       % Always plot quiver
        
        triquiver(U,V,alldata.mesh);
        
        
        if animation(1)>0, figure(4)    % Variable 1
            switch animation(1)
                case 1, H = U;
                case 2, H = V;
                case 3, H = sqrt(U.^2+V.^2);
                case 4, H = P;
                case 5, H = abs(smvp(Gx,V)-smvp(Gy,U));
                case 6, H = S;
            end
            switch animation(3)
                case 1
                    trimesh(t,p(:,1),p(:,2),H);
                case 2
                    trisurf(t,p(:,1),p(:,2),H); shading interp
                case 3
                    tricontour(H,40,alldata.mesh);
            end
            if animation(4)==1
                view(2), axis equal
            else
                view(3)
            end
        end
        if animation(2)>0, figure(5)    % Variable 2
            switch animation(2)
                case 1, H = U;
                case 2, H = V;
                case 3, H = sqrt(U.^2+V.^2);
                case 4, H = P;
                case 5, H = abs(smvp(Gx,V)-smvp(Gy,U));
                case 6, H = S;
            end
            switch animation(3)
                case 1
                    trimesh(t,p(:,1),p(:,2),H);
                case 2
                    trisurf(t,p(:,1),p(:,2),H); shading interp
                case 3
                    tricontour(H,40,alldata.mesh);
            end
            if animation(4)==1
                view(2), axis equal
            else
                view(3)
            end
        end

        drawnow
        
    end
       
    % Advance counters
    m         = m+1;
    time(m+1) = time(m)+dt;
    
end

% Save back to GUI data
alldata.init.U    = U;
alldata.init.V    = V;
alldata.init.S    = S;
alldata.init.P    = P;
alldata.init.W    = smvp(Gx,V)-smvp(Gy,U);
alldata.init.time = time(1:mmax);
alldata.init.resx = resx(1:mmax);
alldata.init.resy = resy(1:mmax);

set(figure(1),'UserData',alldata);


% Print cost stats 
cost = {
       ['Total time: ',num2str(cputime-tstart),' seconds']
       ''
       ['Total steps: ',num2str(m)]
       ''
       ['Non-linear evals: ',num2str(tN),' seconds']
       ''
       ['Viscous evals: ',num2str(tV),' seconds']
       ''
       ['Pressure evals: ',num2str(tQ),' seconds']
       };
       
msgbox(cost,'Done','none')  

return



% SUB-FUNCTIONS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Nu,Nv] = nonlinear_uv(u,ux,uy,v,vx,vy,e,N3,in,hnx,hny,dx,dy);

% Assemble the non-linear vectors for the velocities

% Loop over edges
Nu = 0*u; Nv = Nu;
for k = 1:size(e,1)
    if in(k)
        
        % Nodes
        n1 = e(k,1);  n2 = e(k,2); 
        
        % Nodal velocities
        u1 = u(n1);  u2 = u(n2);  um = 5*(u1+u2)/12;
        v1 = v(n1);  v2 = v(n2);  vm = 5*(v1+v2)/12;
        
        % Velocity coeffs (median midpoint)
        um1 = um + u(N3(k,1))/6;  um2 = um + u(N3(k,2))/6;  
        vm1 = vm + v(N3(k,1))/6;  vm2 = vm + v(N3(k,2))/6;  
        
        % Integral face normal velocities
        un1 = hnx(k,1)*um1 + hny(k,1)*vm1;
        un2 = hnx(k,2)*um2 + hny(k,2)*vm2;
        
        % EDGE 1
        if un1<0
            % Extrapolation
            rGu = dx(k,2)*ux(n2)+dy(k,2)*uy(n2);
            rGv = dx(k,2)*vx(n2)+dy(k,2)*vy(n2);
            % Gradient ratios
            if rGu~=0,  Ru = (um1-u2)/rGu;  else  Ru = 0;  end
            if rGv~=0,  Rv = (vm1-v2)/rGv;  else  Rv = 0;  end
            % HQUICK
            if Ru<=0,  psiu = 0;  else  psiu = 4*Ru*rGu/(3+Ru);  end
            if Rv<=0,  psiv = 0;  else  psiv = 4*Rv*rGv/(3+Rv);  end
            % Upwind flux
            cu = un1*(u2+psiu);
            cv = un1*(v2+psiv);
        else
            % Extrapolation
            rGu = dx(k,1)*ux(n1)+dy(k,1)*uy(n1);
            rGv = dx(k,1)*vx(n1)+dy(k,1)*vy(n1);
            % Gradient ratios
            if rGu~=0,  Ru = (um1-u1)/rGu;  else  Ru = 0;  end
            if rGv~=0,  Rv = (vm1-v1)/rGv;  else  Rv = 0;  end
            % HQUICK
            if Ru<=0,  psiu = 0;  else  psiu = 4*Ru*rGu/(3+Ru);  end
            if Rv<=0,  psiv = 0;  else  psiv = 4*Rv*rGv/(3+Rv);  end
            % Upwind flux
            cu = un1*(u1+psiu);
            cv = un1*(v1+psiv);
        end
        
        % EDGE 2
        if un2<0
            % Extrapolation
            rGu = dx(k,4)*ux(n2)+dy(k,4)*uy(n2);
            rGv = dx(k,4)*vx(n2)+dy(k,4)*vy(n2);
            % Gradient ratios
            if rGu~=0,  Ru = (um2-u2)/rGu;  else  Ru = 0;  end
            if rGv~=0,  Rv = (vm2-v2)/rGv;  else  Rv = 0;  end
            % HQUICK
            if Ru<=0,  psiu = 0;  else  psiu = 4*Ru*rGu/(3+Ru);  end
            if Rv<=0,  psiv = 0;  else  psiv = 4*Rv*rGv/(3+Rv);  end
            % Upwind flux
            cu = cu + un2*(u2+psiu);
            cv = cv + un2*(v2+psiv);
        else
            % Extrapolation
            rGu = dx(k,3)*ux(n1)+dy(k,3)*uy(n1);
            rGv = dx(k,3)*vx(n1)+dy(k,3)*vy(n1);
            % Gradient ratios
            if rGu~=0,  Ru = (um2-u1)/rGu;  else  Ru = 0;  end
            if rGv~=0,  Rv = (vm2-v1)/rGv;  else  Rv = 0;  end
            % HQUICK
            if Ru<=0,  psiu = 0;  else  psiu = 4*Ru*rGu/(3+Ru);  end
            if Rv<=0,  psiv = 0;  else  psiv = 4*Rv*rGv/(3+Rv);  end
            % Upwind flux
            cu = cu + un2*(u1+psiu);
            cv = cv + un2*(v1+psiv);
        end
        
        % NODAL VECTORS