tvd_rk2.m

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

        Nu(n1) = Nu(n1)+cu;  Nu(n2) = Nu(n2)-cu;
        Nv(n1) = Nv(n1)+cv;  Nv(n2) = Nv(n2)-cv;
        
    end
end

return


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

% Assemble the non-linear vectors for the velocities and the scalar field

% Loop over edges
Nu = 0*u; Nv = Nu; Ns = 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;
        s1 = s(n1);  s2 = s(n2);  sm = 5*(s1+s2)/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;
        sm1 = sm + s(N3(k,1))/6;  sm2 = sm + s(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);
            rGs = dx(k,2)*sx(n2)+dy(k,2)*sy(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
            if rGs~=0,  Rs = (sm1-s2)/rGs;  else  Rs = 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
            if Rs<=0,  psis = 0;  else  psis = 4*Rs*rGs/(3+Rs);  end
            % Upwind flux
            cu = un1*(u2+psiu);
            cv = un1*(v2+psiv);
            cs = un1*(s2+psis);
        else
            % Extrapolation
            rGu = dx(k,1)*ux(n1)+dy(k,1)*uy(n1);
            rGv = dx(k,1)*vx(n1)+dy(k,1)*vy(n1);
            rGs = dx(k,1)*sx(n1)+dy(k,1)*sy(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
            if rGs~=0,  Rs = (sm1-s1)/rGs;  else  Rs = 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
            if Rs<=0,  psis = 0;  else  psis = 4*Rs*rGs/(3+Rs);  end
            % Upwind flux
            cu = un1*(u1+psiu);
            cv = un1*(v1+psiv);
            cs = un1*(s1+psis);
        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);
            rGs = dx(k,4)*sx(n2)+dy(k,4)*sy(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
            if rGs~=0,  Rs = (sm2-s2)/rGs;  else  Rs = 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
            if Rs<=0,  psis = 0;  else  psis = 4*Rs*rGs/(3+Rs);  end
            % Upwind flux
            cu = cu + un2*(u2+psiu);
            cv = cv + un2*(v2+psiv);
            cs = cs + un2*(s2+psis);
        else
            % Extrapolation
            rGu = dx(k,3)*ux(n1)+dy(k,3)*uy(n1);
            rGv = dx(k,3)*vx(n1)+dy(k,3)*vy(n1);
            rGs = dx(k,3)*sx(n1)+dy(k,3)*sy(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
            if rGs~=0,  Rs = (sm2-s1)/rGs;  else  Rs = 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
            if Rs<=0,  psis = 0;  else  psis = 4*Rs*rGs/(3+Rs);  end
            % Upwind flux
            cu = cu + un2*(u1+psiu);
            cv = cv + un2*(v1+psiv);
            cs = cs + un2*(s1+psis);
        end
        
        % NODAL VECTORS
        Nu(n1) = Nu(n1)+cu;  Nu(n2) = Nu(n2)-cu;
        Nv(n1) = Nv(n1)+cv;  Nv(n2) = Nv(n2)-cv;
        Ns(n1) = Ns(n1)+cs;  Ns(n2) = Ns(n2)-cs;
        
    end
end

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = bnd_extrap(f,fx,fy,vec,p,n2n) 

% Set nodal values in VEC by average linear extrapolation 
% from neighbours

for iter = 1:4          % Coupled, so do some iters to converge
    for k = 1:length(vec)
        % Take average of neighbours
        z = 1; veck = vec(k); sumf = f(veck);
        while n2n(veck,z)>0
            % Neighbour
            nn = n2n(veck,z);
            % Extrapolation
            sumf = sumf + f(nn) + 0.5*( (p(veck,1)-p(nn,1))*(fx(veck)+fx(nn)) ...
                                      + (p(veck,2)-p(nn,2))*(fy(veck)+fy(nn)) ); 
            % Counter
            z = z+1;
        end
        f(veck) = sumf/z;
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,iter] = myBiCGSTAB(V,M_dt,b,x,tol,maxit,diagV)

% BiCGSTAB iterative method for the viscous terms in the transport equations. 
% Calls mex "smvp.c" for a fast matrix-vector product. Jacobi
% pre-conditioned.

% Scale tolerance wrt RHS
tol = tol*norm(b,inf);

% MAIN LOOP
if tol==0                   % All zero rhs has all zero solution
    
    x = 0*x;  iter = 0;
    
else                        % BiCGSTAB method
    
    % Initial residual
    r     = x;              % Pre-alloc (important for speed...)
    r     = b - M_dt.*x + smvp(V,x);
    rdash = r'- sqrt(eps);  % Fudge by sqrt(eps) - gives better convergence??
    
    % Initialise
    a    = 0; p = r;
    v    = r; w = 1;
    rhol = 1;

    % Jacobi preconditioner
    M = 1./(M_dt-diagV);
    
    for iter = 1:maxit
    
        rho = rdash*r;

        rholw = rhol*w;
        if rholw==0 
            error('BiCGSTAB failure')
        end

        % Search vector
        p = r + (rho*a)*(p - w*v)/rholw;

        % HALF ITERATE STEP        
        % Jacobi Preconditioner
        ph = M.*p;     
        v  = M_dt.*ph-smvp(V,ph);        
        a  = rho/(rdash*v);
    
        % Residual for half iterate
        s = r - a*v;
    
        % FULL ITERATE STEP
        % Jacobi Preconditioner
        sh = M.*s;   
        t  = M_dt.*sh-smvp(V,sh);
        tt = t';
        w  = (tt*s)/(tt*t);
    
        % Solution & residual at full iterate
        x = x + a*ph + w*sh;
        r = s - w*t;
    
        % Break iteration if converged
        if norm(r,inf)<=tol
            return
        end
    
        % For the next step
        rhol = rho;
    
    end
    % iter==maxit to get here
    disp('WARNING: BiCGSTAB did not converge. Try decreasing the CFL number.')
end

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function tricontour(Hn,num,data)

% This is the guts of my TRICONTOUR function from the FEX.

p    = data.p;      % Nodes
e    = data.e;      % Edges
t    = data.t;      % Triangulation
pt   = data.pc;     % Centroids
pe   = data.pe;     % Edge midpoints
e2t  = data.e2t;    % Edge to triangle connectivity    
eINt = data.eINt;   % Triangles as edges
numt = size(t,1);   % Number of triangles

% Endpoint nodes in edges
e1 = e(:,1); e2 = e(:,2);

% Interpolate to centroids and edge midpoints
Ht = (Hn(t(:,1))+Hn(t(:,2))+Hn(t(:,3)))/3;   
He = (Hn(e1)+Hn(e2))/2;  

% Data increment
dH = (max(Ht)-min(Ht))/(num+1);

% MAIN LOOP
x   = []; y = []; z = [];
in  = false(numt,1);
lev = max(Ht);
vec = 1:numt;
old = in;
for v = 1:num       % Loop over contouring levels

    % Find centroid values >= current level
    lev   = lev-dH;
    i     = vec(Ht>=lev);
    i     = i(~old(i));     % Don't need to check triangles from higher levels
    in(i) = true;

    % Locate boundary edges in group
    bnd  = [i; i; i];
    next = 1;
    for k = 1:length(i)
        ct    = i(k);
        count = 0;
        for q = 1:3     % Loop through edges in ct
            ce = eINt(ct,q);
            if ~in(e2t(ce,1)) || ((e2t(ce,2)>0)&&~in(e2t(ce,2)))    
                bnd(next) = ce;     % Found bnd edge
                next      = next+1;
            else
                count = count+1;    % Count number of non-bnd edges in ct
            end
        end
        if count==3                 % If 3 non-bnd edges ct must be in middle of group
            old(ct) = true;         % & doesn't need to be checked for the next level
        end
    end
    bnd = bnd(1:next-1);
    
    % Place nodes approximately on contours by interpolating across bnd edges
    numb  = next-1;
    penew = pe;
    cc    = repmat(0,2*numb,2);
    for k = 1:numb
        
        ce = bnd(k);        % Current edge
        t1 = e2t(ce,1);     % Associated
        t2 = e2t(ce,2);     % triangles
        
        if t2>0     % Move node based on neighbouring centroids
            dx = pt(t2,1)-pt(t1,1);
            dy = pt(t2,2)-pt(t1,2);
            r  = (lev-Ht(t1))/(Ht(t2)-Ht(t1));
        else        % Move node based on internal centroid & bnd midpoint
            dx = pe(ce,1)-pt(t1,1);
            dy = pe(ce,2)-pt(t1,2);
            r  = (lev-Ht(t1))/(He(ce)-Ht(t1));
        end
        
        % New node position
        penew(ce,1) = pt(t1,1)+r*dx;
        penew(ce,2) = pt(t1,2)+r*dy;
        
        % Do a temp connection between adjusted node & endpoint nodes in
        % ce so that the connectivity between neighbouring adjusted nodes
        % can be determined
        m         = 2*k-1;
        cc(m,1)   = e1(ce);
        cc(m,2)   = ce;
        cc(m+1,1) = e2(ce);
        cc(m+1,2) = ce;
        
    end
    
    % Sort connectivity to place connected edges in sucessive rows
    [j,i] = sort(cc(:,1));
    cc    = cc(i,1:2);
    
    % Connect adjacent adjusted nodes
    k    = 1;
    next = 1;
    while k<(2*numb)
        if cc(k,1)==cc(k+1,1)
            cc(next,1) = cc(k,2);
            cc(next,2) = cc(k+1,2);
            next       = next+1;
            k          = k+2;       % Skip over connected edge
        else
            k = k+1;                % Node has only 1 connection - will be picked up above
        end
    end
    cc = cc(1:next-1,1:2);
    
    % xzy data
    xc = [penew(cc(:,1),1),penew(cc(:,2),1)];
    yc = [penew(cc(:,1),2),penew(cc(:,2),2)];
    zc = lev + 0*xc;
    
    % Add contours to list
    x = [x; xc];
    y = [y; yc];
    z = [z; zc];

end

% Draw contours
newplot, patch('Xdata',x','Ydata',y','Zdata',z','Cdata',z','facecolor','none','edgecolor','flat');

axis([min(p(:,1)),max(p(:,1)),min(p(:,2)),max(p(:,2))])

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function triquiver(u,v,data)

% Draw a quiver plot of the velocity field, but do NOT scale the length of
% the arrows with magnitude. Instead colour the arrows with magnitude. Also
% plot a few contours of total velocity using tricontour.

x = data.p(:,1);
y = data.p(:,2);
A = data.A;

uv    = sqrt(u.^2+v.^2);                % Total velocity  
scale = 0.75*sqrt(A)./max(uv,eps);      % Scale arrow length with cell area

xend   = x+scale.*u;
yend   = y+scale.*v;
xarrow = 0.3*x+0.7*xend;
yarrow = 0.3*y+0.7*yend;

newplot, patch('xdata'    ,[x, xend; xarrow+0.2*scale.*v, xend; xarrow-0.2*scale.*v, xend]', ...
               'ydata'    ,[y, yend; yarrow-0.2*scale.*u, yend; yarrow+0.2*scale.*u, yend]', ...
               'cdata'    ,[uv,uv; uv,uv; uv,uv]'                                          , ... 
               'facecolor','none'                                                          , ...
               'edgecolor','flat');

axis equal, axis off, hold on

tricontour(uv,10,data); hold off

return