%% Program to implement multi-species RK3 algorithm described in
%  "Computational fluctuating fluid dynamics", J.B. Bell, A. Garcia, and
%  S. Williams, to appear in ESAIM: Mathematical Modelling and Numerical Analysis.
%  MATLAB program written by Anton de la Fuente
clear all;
set(0,'defaulttextinterpreter','latex')
disp(' ')
%% Set physical parameters and system constants
k_B = 1.38054e-16;  % Boltzmann's constant
gamma = 5/3;        % Specific heat ratio
gammaM1 = gamma-1;
mass1 = 6.63e-23;
massR = 3;
mass0 = massR*mass1;
CmassEq = .5;
diam0 = 3.66e-8;    % Particle diameter
diam1 = 3.66e-8;

Rgas0 = k_B/mass0;  % Gas constant
Rgas1 = k_B/mass1;
Cv0 = 1.5*Rgas0;    % Specific heat per particle
Cv1 = 1.5*Rgas1;
Cp0 = gamma*Cv0;
Cp1 = gamma*Cv1;
rhoEq = 1.78e-3;    % Equilibrium mass density
tempEq = 273;       % Equilibrium temperature
presEq = rhoEq*((1-CmassEq)*Rgas0 + CmassEq*Rgas1)*tempEq;
cSound = sqrt(gamma*presEq/rhoEq); % Sound speed
velEq = 0;%cSound;                    % Equilibrium fluid velocity (can be zero)
Jeq = rhoEq*velEq;                 % Equilibrium momentum density
Eeq = presEq/gammaM1 + 0.5*rhoEq*velEq^2;    % Energy density


nxCells = 20;  
nyCells = 20;  
nzCells = 20;    

nSteps = 100000;

% Auxiliary variables for calculating transport coefficients
Cn = 1.016;
Ck = 1.025;
Cd = 1.133;
Cs = 1.299;
massx = 2*mass0*mass1/(mass0+mass1);
diamx = (diam0 + diam1)/2;
auxmpm = (mass0 + mass1)^2/(4*mass0*mass1);
auxmmm = (mass0 - mass1)^2/(2*mass0*mass1);
U0 = 4/15 - 17/60*mass0/mass1 + auxmmm;
U1 = 4/15 - 17/60*mass1/mass0 + auxmmm;

alpha1 = 0.25 * (sqrt(7)+1);                  % Interpolation coefficient
alpha2 = 0.25 * (sqrt(7)-1);                  % Interpolation coefficient

%% Initialize the grid values
dx = 2.7e-6;   dy = dx;   dz = dx;            % Grid spacing 
cellVol = dx*dy*dz;                           % Volume of single grid point
rho  = rhoEq *ones(nxCells,nyCells,nzCells);     % Mass density (\rho)
rho1  = rho*CmassEq; 
Jx = Jeq*ones(nxCells,nyCells,nzCells);          % Momentum density (J_x)
Jy = zeros(nxCells,nyCells,nzCells);             % Momentum density (J_y)
Jz = zeros(nxCells,nyCells,nzCells);             % Momentum density (J_z)
E  = Eeq *ones(nxCells,nyCells,nzCells);         % Energy density (E)

if nxCells == 1
    ip = 1; im = 1; ipp = 1;
else
    ip  = [2:nxCells, 1];                      % i+1, with periodic b.c.
    im  = [nxCells, 1:(nxCells-1)];            % i-1, with periodic b.c.
    ipp = [3:nxCells, 1, 2];                   % i+2, with periodic b.c.
end

if nyCells == 1
    jp = 1; jm = 1; jpp = 1;
else
    jp  = [2:nyCells, 1];                      % i+1, with periodic b.c.
    jm  = [nyCells, 1:(nyCells-1)];            % i-1, with periodic b.c.
    jpp = [3:nyCells, 1, 2];                   % i+2, with periodic b.c.
end

if nzCells == 1
    kp = 1; km = 1; kpp = 1;
else
    kp  = [2:nzCells, 1];                      % i+1, with periodic b.c.
    km  = [nzCells, 1:(nzCells-1)];            % i-1, with periodic b.c.
    kpp = [3:nzCells, 1, 2];                   % i+2, with periodic b.c.
end


%% Initialize statistics variables
rhoMean  = zeros(nxCells,nyCells,nzCells);     rhoVar = zeros(nxCells,nyCells,nzCells);
JxMean   = zeros(nxCells,nyCells,nzCells);     JxVar  = zeros(nxCells,nyCells,nzCells);
Emean  = zeros(nxCells,nyCells,nzCells);       Evar = zeros(nxCells,nyCells,nzCells);
rho1Mean  = zeros(nxCells,nyCells,nzCells);    rho1Var = zeros(nxCells,nyCells,nzCells);
rhoJxCov = zeros(nxCells,nyCells,nzCells);     
rhoECov = zeros(nxCells,nyCells,nzCells); 
JxECov = zeros(nxCells,nyCells,nzCells); 
rhoAvgInit = mean(rho(:));  JxAvgInit  = mean(Jx(:));  EavgInit = mean(E(:));

%%%%%% Main loop %%%%%%
dt = 1.0e-12;   % Time step
dtOdx = dt/dx;  dtOdy = dt/dy;  dtOdz = dt/dz;
cCoeff = sqrt(2)*sqrt(1/(dt*cellVol));
sCoeffii = sqrt(2)*sqrt(4*k_B/(3*dt*cellVol));% Momentum noise coefficient for S_{ij} when i = j
sCoeffij = sqrt(2)*sqrt(k_B/(dt*cellVol));    % Momentum noise coefficient for S_{ij} when i \neq j
qCoeff = sqrt(2)*sqrt(k_B/(dt*cellVol));      % Energy noise coefficient

nSkip = nSteps/4; nSamp = 0;    			  % Start taking statistics after nSkip steps
tic;    % Start MATLAB timer
for iStep=1:nSteps

    for iStage=1:3  % RK3 has three stages
        
        %% Set current values for this stage
        switch( iStage ) 
            case {1}
                rho_ = rho;    rho1_ = rho1;    
                Jx_ = Jx;      Jy_ = Jy;     Jz_ = Jz;
                E_ = E;
            case {2}
                rho_ = rho_13; rho1_ = rho1_13;       
                Jx_ = Jx_13;   Jy_ = Jy_13;   Jz_ = Jz_13; 
                E_ = E_13;
            case {3}
                rho_ = rho_23;  rho1_ = rho1_23;     
                Jx_ = Jx_23;   Jy_ = Jy_23;   Jz_ = Jz_23;
                E_ = E_23;
        end

        %% Compute conserved variables on half-grid by special interpolation
        %% Notation:  rho_   = \rho_{i,j,k}     (Full grid) 
        %%            rho_ip = \rho_{i+1/2,j,k} (Half grid)
        %%            rho_jp = \rho_{i,j+1/2,k}
        %%            rho_kp = \rho_{i,j,k+1/2}
        rho_ip = alpha1*( rho_(ip,:,:) + rho_ ) - alpha2*( rho_(ipp,:,:) + rho_(im,:,:) );
        rho_jp = alpha1*( rho_(:,jp,:) + rho_ ) - alpha2*( rho_(:,jpp,:) + rho_(:,jm,:) );
        rho_kp = alpha1*( rho_(:,:,kp) + rho_ ) - alpha2*( rho_(:,:,kpp) + rho_(:,:,km) );
        
        rho1_ip = alpha1*( rho1_(ip,:,:) + rho1_ ) - alpha2*( rho1_(ipp,:,:) + rho1_(im,:,:) );
        rho1_jp = alpha1*( rho1_(:,jp,:) + rho1_ ) - alpha2*( rho1_(:,jpp,:) + rho1_(:,jm,:) );
        rho1_kp = alpha1*( rho1_(:,:,kp) + rho1_ ) - alpha2*( rho1_(:,:,kpp) + rho1_(:,:,km) );
        
        Jx_ip = alpha1*( Jx_(ip,:,:) + Jx_ ) - alpha2*( Jx_(ipp,:,:) + Jx_(im,:,:) );
        Jx_jp = alpha1*( Jx_(:,jp,:) + Jx_ ) - alpha2*( Jx_(:,jpp,:) + Jx_(:,jm,:) );
        Jx_kp = alpha1*( Jx_(:,:,kp) + Jx_ ) - alpha2*( Jx_(:,:,kpp) + Jx_(:,:,km) );
        
        Jy_ip = alpha1*( Jy_(ip,:,:) + Jy_ ) - alpha2*( Jy_(ipp,:,:) + Jy_(im,:,:) );
        Jy_jp = alpha1*( Jy_(:,jp,:) + Jy_ ) - alpha2*( Jy_(:,jpp,:) + Jy_(:,jm,:) );
        Jy_kp = alpha1*( Jy_(:,:,kp) + Jy_ ) - alpha2*( Jy_(:,:,kpp) + Jy_(:,:,km) );
        
        Jz_ip = alpha1*( Jz_(ip,:,:) + Jz_ ) - alpha2*( Jz_(ipp,:,:) + Jz_(im,:,:) );
        Jz_jp = alpha1*( Jz_(:,jp,:) + Jz_ ) - alpha2*( Jz_(:,jpp,:) + Jz_(:,jm,:) );
        Jz_kp = alpha1*( Jz_(:,:,kp) + Jz_ ) - alpha2*( Jz_(:,:,kpp) + Jz_(:,:,km) );
        
        E_ip = alpha1*( E_(ip,:,:) + E_ ) - alpha2*( E_(ipp,:,:) + E_(im,:,:) );
        E_jp = alpha1*( E_(:,jp,:) + E_ ) - alpha2*( E_(:,jpp,:) + E_(:,jm,:) );
        E_kp = alpha1*( E_(:,:,kp) + E_ ) - alpha2*( E_(:,:,kpp) + E_(:,:,km) );

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% Compute hydrodynamic variables %%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Cmass_  = rho1_./rho_;
        Cnum_   = massR*Cmass_./(1. - (1 - massR)*Cmass_);
        
        velX_   = Jx_   ./ rho_  ;                                 
        velX_ip = Jx_ip ./ rho_ip;  
        velX_jp = Jx_jp ./ rho_jp;
        velX_kp = Jx_kp ./ rho_kp;
        
        velY_   = Jy_   ./ rho_  ;
        velY_ip = Jy_ip ./ rho_ip;                          
        velY_jp = Jy_jp ./ rho_jp;
        velY_kp = Jy_kp ./ rho_kp;
        
        velZ_   = Jz_   ./ rho_  ;
        velZ_ip = Jz_ip ./ rho_ip;                          
        velZ_jp = Jz_jp ./ rho_jp;
        velZ_kp = Jz_kp ./ rho_kp;
        
        vel2_ = velX_.^2 + velY_.^2 + velZ_.^2; 
        vel2_ip = velX_ip.^2 + velY_ip.^2 + velZ_ip.^2;      
        vel2_jp = velX_jp.^2 + velY_jp.^2 + velZ_jp.^2;
        vel2_kp = velX_kp.^2 + velY_kp.^2 + velZ_kp.^2;
       
        pres_ = (E_ - 0.5*(rho_.*vel2_))*gammaM1;
        pres_ip = (E_ip - 0.5*(rho_ip.*vel2_ip))*gammaM1;
        pres_jp = (E_jp - 0.5*(rho_jp.*vel2_jp))*gammaM1;
        pres_kp = (E_kp - 0.5*(rho_kp.*vel2_kp))*gammaM1;
        
        temp_ = (E_./rho_ - 0.5*vel2_)./(Cv0*(1-Cmass_) + Cv1.*Cmass_);
        temp_ip = (temp_(ip,:,:) + temp_)/2;
        temp_jp = (temp_(:,jp,:) + temp_)/2;
        temp_kp = (temp_(:,:,kp) + temp_)/2;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% Compute transport coefficients %%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        sqrtT_ = sqrt(temp_);
        
        %% Viscosity
        n0_ = 5/(16*diam0^2) * sqrt(mass0*k_B/pi)*sqrtT_;
        n1_ = 5/(16*diam1^2) * sqrt(mass1*k_B/pi)*sqrtT_;      
        nx_ = 5/(16*diamx^2) * sqrt(massx*k_B/pi)*sqrtT_;
      
        mC2_ = (1 - Cnum_).^2;  C2_ = Cnum_.^2;  CmC2_ = 2*Cnum_ - 2*C2_ ;      
        Xn_ = mC2_./n0_ + CmC2_./nx_ + C2_./n1_;    
        Yn_ = 3/5*( mC2_./n0_*massR + CmC2_.*nx_./(n0_.*n1_)*auxmpm + C2_./(n1_*massR) );  
        Zn_ = 3/5*( mC2_*massR + CmC2_.* (auxmpm*(nx_./n0_ + nx_./n1_) - 1) + C2_/massR );

        visc_ = Cn*(1 + Zn_)./(Xn_ + Yn_);
        visc_ip = (visc_(ip,:,:) + visc_)/2;
        visc_jp = (visc_(:,jp,:) + visc_)/2;
        visc_kp = (visc_(:,:,kp) + visc_)/2;
        
        %% Thermal Conductivity
        k0_ = 15*k_B/(4*mass0)*n0_;        
        k1_ = 15*k_B/(4*mass1)*n1_;       
        kx_ = 15*k_B/(4*massx)*nx_;
        
        Uy_ = 4/15*auxmpm*kx_.^2./(k0_.*k1_) - 17/60 + 13/16*auxmmm;        
        Uz_ = 4/15*(auxmpm*(kx_./k0_ + kx_./k1_) - 1) - 17/60;

        Xk_ = mC2_./k0_ + CmC2_./kx_ + C2_./k1_;
        Yk_ = mC2_.*U0./k0_ + CmC2_.*Uy_./kx_ + C2_.*U1./k1_;        
        Zk_ = mC2_.*U0 + CmC2_.*Uz_ + C2_.*U1;
  
        thermC_ = Ck*(1 + Zk_)./(Xk_ + Yk_);
        
        %% Mass diffusion coefficient
        D_ = 6*Cd*nx_./5;
        
        %% Thermal diffusion coefficient
        S0_ = (mass0+mass1)/(2*mass1)*kx_./k0_ - 15/8*(mass1-mass0)/mass0 - 1;        
        S1_ = (mass0+mass1)/(2*mass0)*kx_./k1_ + 15/8*(mass1-mass0)/mass1 - 1;
        kt_ = Cs * CmC2_./(12.*kx_) .* (S1_.*Cnum_ - S0_.*(1-Cnum_))./(Xk_ + Yk_);
        
        %% Baro diffusion coefficient
        kp_ = (mass0-mass1)*Cmass_.*(1-Cmass_) .* ((1-Cmass_)/mass0 + Cmass_/mass1);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% Compute gradients of hydrodynamic variables %%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        dvelXdx_ip = ( velX_(ip,:,:) - velX_ )/dx;
        dvelYdx_ip = ( velY_(ip,:,:) - velY_ )/dx;
        dvelZdx_ip = ( velZ_(ip,:,:) - velZ_ )/dx;
        
        dvelXdx_jp = ( velX_(ip,jp,:) + velX_(ip,:,:) - velX_(im,jp,:) - velX_(im,:,:) )/(4*dx);
        dvelYdx_jp = ( velY_(ip,jp,:) + velY_(ip,:,:) - velY_(im,jp,:) - velY_(im,:,:) )/(4*dx);
        dvelZdx_jp = ( velZ_(ip,jp,:) + velZ_(ip,:,:) - velZ_(im,jp,:) - velZ_(im,:,:) )/(4*dx);
        
        dvelXdx_kp = ( velX_(ip,:,kp) + velX_(ip,:,:) - velX_(im,:,kp) - velX_(im,:,:) )/(4*dx);
        dvelYdx_kp = ( velY_(ip,:,kp) + velY_(ip,:,:) - velY_(im,:,kp) - velY_(im,:,:) )/(4*dx);
        dvelZdx_kp = ( velZ_(ip,:,kp) + velZ_(ip,:,:) - velZ_(im,:,kp) - velZ_(im,:,:) )/(4*dx);

		dvelXdy_ip = ( velX_(ip,jp,:) + velX_(:,jp,:) - velX_(ip,jm,:) - velX_(:,jm,:) )/(4*dy);
        dvelYdy_ip = ( velY_(ip,jp,:) + velY_(:,jp,:) - velY_(ip,jm,:) - velY_(:,jm,:) )/(4*dy);
        dvelZdy_ip = ( velZ_(ip,jp,:) + velZ_(:,jp,:) - velZ_(ip,jm,:) - velZ_(:,jm,:) )/(4*dy);
        
        dvelXdy_jp = ( velX_(:,jp,:) - velX_ )/dy;
        dvelYdy_jp = ( velY_(:,jp,:) - velY_ )/dy;
        dvelZdy_jp = ( velZ_(:,jp,:) - velZ_ )/dy;
        
		dvelXdy_kp = ( velX_(:,jp,kp) + velX_(:,jp,:) - velX_(:,jm,kp) - velX_(:,jm,:) )/(4*dy);
        dvelYdy_kp = ( velY_(:,jp,kp) + velY_(:,jp,:) - velY_(:,jm,kp) - velY_(:,jm,:) )/(4*dy);
        dvelZdy_kp = ( velZ_(:,jp,kp) + velZ_(:,jp,:) - velZ_(:,jm,kp) - velZ_(:,jm,:) )/(4*dy);
		
		dvelXdz_ip = ( velX_(ip,:,kp) + velX_(:,:,kp) - velX_(ip,:,km) - velX_(:,:,km) )/(4*dz);
        dvelYdz_ip = ( velY_(ip,:,kp) + velY_(:,:,kp) - velY_(ip,:,km) - velY_(:,:,km) )/(4*dz);
        dvelZdz_ip = ( velZ_(ip,:,kp) + velZ_(:,:,kp) - velZ_(ip,:,km) - velZ_(:,:,km) )/(4*dz);
        
		dvelXdz_jp = ( velX_(:,jp,kp) + velX_(:,:,kp) - velX_(:,jp,km) - velX_(:,:,km) )/(4*dz);
        dvelYdz_jp = ( velY_(:,jp,kp) + velY_(:,:,kp) - velY_(:,jp,km) - velY_(:,:,km) )/(4*dz);
        dvelZdz_jp = ( velZ_(:,jp,kp) + velZ_(:,:,kp) - velZ_(:,jp,km) - velZ_(:,:,km) )/(4*dz);
		
        dvelXdz_kp = ( velX_(:,:,kp) - velX_ )/dz;
        dvelYdz_kp = ( velY_(:,:,kp) - velY_ )/dz;
        dvelZdz_kp = ( velZ_(:,:,kp) - velZ_ )/dz;
        
        dtempdx = ( temp_(ip,:,:) - temp_ )/dx;
        dtempdy = ( temp_(:,jp,:) - temp_ )/dy;
        dtempdz = ( temp_(:,:,kp) - temp_ )/dz;
        
        dpresdx = ( pres_(ip,:,:) - pres_ )/dx;
        dpresdy = ( pres_(:,jp,:) - pres_ )/dy;
        dpresdz = ( pres_(:,:,kp) - pres_ )/dz;
        
        dCmassdx = ( Cmass_(ip,:,:) - Cmass_ )/dx;
        dCmassdy = ( Cmass_(:,jp,:) - Cmass_ )/dy;
        dCmassdz = ( Cmass_(:,:,kp) - Cmass_ )/dz;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% Compute stochastic fluxes %%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        if iStage == 1
            A_ = Cmass_.*(1-Cmass_).*(mass1*(1-Cmass_) + mass0.*Cmass_);
            Cx = cCoeff * sqrt(D_(ip,:,:).*A_(ip,:,:) + D_.*A_) .* randn(nxCells,nyCells,nzCells);
            Cy = cCoeff * sqrt(D_(:,jp,:).*A_(:,jp,:) + D_.*A_) .* randn(nxCells,nyCells,nzCells);
            Cz = cCoeff * sqrt(D_(:,:,kp).*A_(:,:,kp) + D_.*A_) .* randn(nxCells,nyCells,nzCells);
        
            SxAux = sqrt(visc_(ip,:,:).*temp_(ip,:,:) + visc_.*temp_);
            Sxx = sCoeffii * SxAux .* randn(nxCells,nyCells,nzCells);
            Syx = sCoeffij * SxAux .* randn(nxCells,nyCells,nzCells);
            Szx = sCoeffij * SxAux .* randn(nxCells,nyCells,nzCells);

            SyAux = sqrt(visc_(:,jp,:).*temp_(:,jp,:) + visc_.*temp_);
            Sxy = sCoeffij * SyAux .* randn(nxCells,nyCells,nzCells);
            Syy = sCoeffii * SyAux .* randn(nxCells,nyCells,nzCells);
            Szy = sCoeffij * SyAux .* randn(nxCells,nyCells,nzCells);

            SzAux = sqrt(visc_(:,:,kp).*temp_(:,:,kp) + visc_.*temp_);
            Sxz = sCoeffij * SzAux .* randn(nxCells,nyCells,nzCells);
            Syz = sCoeffij * SzAux .* randn(nxCells,nyCells,nzCells);
            Szz = sCoeffii * SzAux .* randn(nxCells,nyCells,nzCells);

            Qx = qCoeff * sqrt(thermC_(ip,:,:).*temp_(ip,:,:).^2 + thermC_.*temp_.^2) .* randn(nxCells,nyCells,nzCells);
            Qy = qCoeff * sqrt(thermC_(:,jp,:).*temp_(:,jp,:).^2 + thermC_.*temp_.^2) .* randn(nxCells,nyCells,nzCells);
            Qz = qCoeff * sqrt(thermC_(:,:,kp).*temp_(:,:,kp).^2 + thermC_.*temp_.^2) .* randn(nxCells,nyCells,nzCells);
  %      end
        
        %%%%%%%%%%%%%%%%%%%%%%%%
        %% Dissipative Fluxes %%
        %%%%%%%%%%%%%%%%%%%%%%%%
        
        %% Stress-energy tensor
        Txx = visc_ip .*(4/3* dvelXdx_ip - 2/3* (dvelYdy_ip + dvelZdz_ip) );
        Tyx = visc_ip .*(dvelXdy_ip + dvelYdx_ip);
        Tzx = visc_ip .*(dvelXdz_ip + dvelZdx_ip);
        
        Txy = visc_jp .*(dvelYdx_jp + dvelXdy_jp);
        Tyy = visc_jp .*(4/3* dvelYdy_jp - 2/3* (dvelXdx_jp + dvelZdz_jp) );
        Tzy = visc_jp .*(dvelYdz_jp + dvelZdy_jp);
        
        Txz = visc_kp .*(dvelZdx_kp + dvelXdz_kp);
        Tyz = visc_kp .*(dvelZdy_kp + dvelYdz_kp);
        Tzz = visc_kp .*(4/3* dvelZdz_kp - 2/3* (dvelXdx_kp + dvelYdy_kp) );
        
        %% Mass diffusion flux
        jx = (D_(ip,:,:) + D_)/2.*(dCmassdx + ...
            (kt_(ip,:,:)./temp_(ip,:,:) + kt_./temp_)/2.*dtempdx + ...
            (kp_(ip,:,:)./pres_(ip,:,:) + kp_./pres_)/2.*dpresdx); 
        jy = (D_(:,jp,:) + D_)/2.*(dCmassdy + ...
            (kt_(:,jp,:)./temp_(:,jp,:) + kt_./temp_)/2.*dtempdy + ...
            (kp_(:,jp,:)./pres_(:,jp,:) + kp_./pres_)/2.*dpresdy);
        jz = (D_(:,:,kp) + D_)/2.*(dCmassdz + ...
            (kt_(:,:,kp)./temp_(:,:,kp) + kt_./temp_)/2.*dtempdz + ...
            (kp_(:,:,kp)./pres_(:,:,kp) + kp_./pres_)/2.*dpresdz);
        
        G_ip = (Cp1 - Cp0).*temp_ip + 1/2*( kt_(ip,:,:) .* ...
            k_B.*temp_(ip,:,:)./(Cmass_(ip,:,:).*(1-Cmass_(ip,:,:)).*...
            (mass1*(1-Cmass_(ip,:,:))+mass0*Cmass_(ip,:,:))) + ...
            kt_.*k_B.*temp_./(Cmass_.*(1-Cmass_).*...
            (mass1*(1-Cmass_)+mass0*Cmass_)) );
        
        G_jp = (Cp1 - Cp0).*temp_jp + 1/2*( kt_(:,jp,:) .* ...
            k_B.*temp_(:,jp,:)./(Cmass_(:,jp,:).*(1-Cmass_(:,jp,:)).*...
            (mass1*(1-Cmass_(:,jp,:))+mass0*Cmass_(:,jp,:))) + ...
            kt_.*k_B.*temp_./(Cmass_.*(1-Cmass_).*...
            (mass1*(1-Cmass_)+mass0*Cmass_)) );
        
        G_kp = (Cp1 - Cp0).*temp_kp + 1/2*( kt_(:,:,kp) .* ...
            k_B.*temp_(:,:,kp)./(Cmass_(:,:,kp).*(1-Cmass_(:,:,kp)).*...
            (mass1*(1-Cmass_(:,:,kp))+mass0*Cmass_(:,:,kp))) + ...
            kt_.*k_B.*temp_./(Cmass_.*(1-Cmass_).*...
            (mass1*(1-Cmass_)+mass0*Cmass_)) );

        
        %%%%%%%%%%%%%%%%%%%%%%%%%%
        %% Compute total fluxes %%
        %%%%%%%%%%%%%%%%%%%%%%%%%%
        rhoFlux_ip = -Jx_ip;   
        rhoFlux_jp = -Jy_jp;   
        rhoFlux_kp = -Jz_kp;
        
        rho1Flux_ip = -rho1_ip.*velX_ip + jx + Cx;
        rho1Flux_jp = -rho1_jp.*velY_jp + jy + Cy;
        rho1Flux_kp = -rho1_kp.*velZ_kp + jz + Cz;
        
        JxFlux_ip = -Jx_ip .* velX_ip - pres_ip + Txx + Sxx;
        JyFlux_ip = -Jy_ip .* velX_ip + Tyx + Syx;
        JzFlux_ip = -Jz_ip .* velX_ip + Tzx + Szx;
        
        JxFlux_jp = -Jx_jp .* velY_jp + Txy + Sxy;
        JyFlux_jp = -Jy_jp .* velY_jp - pres_jp + Tyy + Syy;
        JzFlux_jp = -Jz_jp .* velY_jp + Tzy + Szy;
        
        JxFlux_kp = -Jx_kp .* velZ_kp + Txz + Sxz;
        JyFlux_kp = -Jy_kp .* velZ_kp + Tyz + Syz;
        JzFlux_kp = -Jz_kp .* velZ_kp - pres_kp + Tzz + Szz;
        
        Eflux_ip = -(E_ip + pres_ip).*velX_ip + (thermC_(ip,:,:)+thermC_)/2.*dtempdx + ...
            velX_ip.*Txx + velY_ip.*Tyx + velZ_ip.*Tzx + G_ip.*jx + Qx + ...
            velX_ip.*Sxx + velY_ip.*Syx + velZ_ip.*Szx + G_ip.*Cx;
        
        Eflux_jp = -(E_jp + pres_jp).*velY_jp + (thermC_(:,jp,:)+thermC_)/2.*dtempdy + ...
            velX_jp.*Txy + velY_jp.*Tyy + velZ_jp.*Tzy + G_jp.*jy + Qy + ...
            velX_jp.*Sxy + velY_jp.*Syy + velZ_jp.*Szy + G_jp.*Cy;
        
        Eflux_kp = -(E_kp + pres_kp).*velZ_kp + (thermC_(:,:,kp)+thermC_)/2.*dtempdz + ...
            velX_kp.*Txz + velY_kp.*Tyz + velZ_kp.*Tzz + G_kp.*jz + Qz + ...
            velX_kp.*Sxz + velY_kp.*Syz + velZ_kp.*Szz + G_kp.*Cz;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% Compute update for this stage %%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        switch( iStage ) 
            case {1}
                rho_13  = rho +  dtOdx * ( rhoFlux_ip - rhoFlux_ip(im,:,:) )+ ...
                                 dtOdy * ( rhoFlux_jp - rhoFlux_jp(:,jm,:) )+ ...
                                 dtOdz * ( rhoFlux_kp - rhoFlux_kp(:,:,km) );
                rho1_13 = rho1 + dtOdx * ( rho1Flux_ip - rho1Flux_ip(im,:,:) )+ ...
                                 dtOdy * ( rho1Flux_jp - rho1Flux_jp(:,jm,:) )+ ...
                                 dtOdz * ( rho1Flux_kp - rho1Flux_kp(:,:,km) );
                Jx_13   = Jx  +  dtOdx * ( JxFlux_ip  - JxFlux_ip(im,:,:) ) + ...
                                 dtOdy * ( JxFlux_jp  - JxFlux_jp(:,jm,:) ) + ...
                                 dtOdz * ( JxFlux_kp  - JxFlux_kp(:,:,km) );
                Jy_13   = Jy  +  dtOdx * ( JyFlux_ip  - JyFlux_ip(im,:,:) ) + ...
                                 dtOdy * ( JyFlux_jp  - JyFlux_jp(:,jm,:) ) + ...
                                 dtOdz * ( JyFlux_kp  - JyFlux_kp(:,:,km) );
                Jz_13   = Jz  +  dtOdx * ( JzFlux_ip  - JzFlux_ip(im,:,:) ) + ...
                                 dtOdy * ( JzFlux_jp  - JzFlux_jp(:,jm,:) ) + ...
                                 dtOdz * ( JzFlux_kp  - JzFlux_kp(:,:,km) );
                E_13    = E  +   dtOdx * ( Eflux_ip - Eflux_ip(im,:,:) ) + ...
                                 dtOdy * ( Eflux_jp - Eflux_jp(:,jm,:) ) + ...
                                 dtOdz * ( Eflux_kp - Eflux_kp(:,:,km) );          
            case {2}
                rho_23  = 0.75*rho + 0.25*rho_13 +   0.25*dtOdx * ( rhoFlux_ip - rhoFlux_ip(im,:,:) ) + ...
                                                     0.25*dtOdy * ( rhoFlux_jp - rhoFlux_jp(:,jm,:) ) + ...
                                                     0.25*dtOdz * ( rhoFlux_kp - rhoFlux_kp(:,:,km) );                                  
                rho1_23 = 0.75*rho1 + 0.25*rho1_13 + 0.25*dtOdx * ( rho1Flux_ip - rho1Flux_ip(im,:,:) ) + ...
                                                     0.25*dtOdy * ( rho1Flux_jp - rho1Flux_jp(:,jm,:) ) + ...
                                                     0.25*dtOdz * ( rho1Flux_kp - rho1Flux_kp(:,:,km) );                               
                Jx_23   = 0.75*Jx  + 0.25*Jx_13  +   0.25*dtOdx * ( JxFlux_ip  - JxFlux_ip(im,:,:) ) + ...
                                                     0.25*dtOdy * ( JxFlux_jp  - JxFlux_jp(:,jm,:) ) + ...
                                                     0.25*dtOdz * ( JxFlux_kp  - JxFlux_kp(:,:,km) );
                Jy_23   = 0.75*Jy  + 0.25*Jy_13  +   0.25*dtOdx * ( JyFlux_ip  - JyFlux_ip(im,:,:) ) + ...
                                                     0.25*dtOdy * ( JyFlux_jp  - JyFlux_jp(:,jm,:) ) + ...
                                                     0.25*dtOdz * ( JyFlux_kp  - JyFlux_kp(:,:,km) );    
                Jz_23   = 0.75*Jz  + 0.25*Jz_13  +   0.25*dtOdx * ( JzFlux_ip  - JzFlux_ip(im,:,:) ) + ...
                                                     0.25*dtOdy * ( JzFlux_jp  - JzFlux_jp(:,jm,:) ) + ...
                                                     0.25*dtOdz * ( JzFlux_kp  - JzFlux_kp(:,:,km) );
                E_23    = 0.75*E + 0.25*E_13  +      0.25*dtOdx * ( Eflux_ip - Eflux_ip(im,:,:) ) + ...
                                                     0.25*dtOdy * ( Eflux_jp - Eflux_jp(:,jm,:) ) + ...
                                                     0.25*dtOdz * ( Eflux_kp - Eflux_kp(:,:,km) );                             
            case {3}
                rho = (1/3)*rho + (2/3)*rho_23 +    (2/3)*dtOdx * ( rhoFlux_ip - rhoFlux_ip(im,:,:) ) + ...
                                                    (2/3)*dtOdy * ( rhoFlux_jp - rhoFlux_jp(:,jm,:) ) + ...
                                                    (2/3)*dtOdz * ( rhoFlux_kp - rhoFlux_kp(:,:,km) );
                rho1 = (1/3)*rho1 + (2/3)*rho1_23 + (2/3)*dtOdx * ( rho1Flux_ip - rho1Flux_ip(im,:,:) ) + ...
                                                    (2/3)*dtOdy * ( rho1Flux_jp - rho1Flux_jp(:,jm,:) ) + ...
                                                    (2/3)*dtOdz * ( rho1Flux_kp - rho1Flux_kp(:,:,km) );
                Jx   = (1/3)*Jx  + (2/3)*Jx_23  +   (2/3)*dtOdx * ( JxFlux_ip  - JxFlux_ip(im,:,:) ) + ...
                                                    (2/3)*dtOdy * ( JxFlux_jp  - JxFlux_jp(:,jm,:) ) + ...
                                                    (2/3)*dtOdz * ( JxFlux_kp  - JxFlux_kp(:,:,km) ); 
                Jy   = (1/3)*Jy  + (2/3)*Jy_23  +   (2/3)*dtOdx * ( JyFlux_ip  - JyFlux_ip(im,:,:) ) + ...
                                                    (2/3)*dtOdy * ( JyFlux_jp  - JyFlux_jp(:,jm,:) ) + ...
                                                    (2/3)*dtOdz * ( JyFlux_kp  - JyFlux_kp(:,:,km) );
                Jz   = (1/3)*Jz  + (2/3)*Jz_23  +   (2/3)*dtOdx * ( JzFlux_ip  - JzFlux_ip(im,:,:) ) + ...
                                                    (2/3)*dtOdy * ( JzFlux_jp  - JzFlux_jp(:,jm,:) ) + ...
                                                    (2/3)*dtOdz * ( JzFlux_kp  - JzFlux_kp(:,:,km) );                             
                E    = (1/3)*E  + (2/3)*E_23    +   (2/3)*dtOdx * ( Eflux_ip - Eflux_ip(im,:,:) ) + ...
                                                    (2/3)*dtOdy * ( Eflux_jp - Eflux_jp(:,jm,:) ) + ...
                                                    (2/3)*dtOdz * ( Eflux_kp - Eflux_kp(:,:,km) );                                          
        end
    end % End of RK3 stages loop

	%% Accumulate statistics
    if( iStep > nSkip )
        nSamp = nSamp + 1;
        rhoMean = rhoMean + rho;             rhoVar = rhoVar + rho.^2;
        JxMean = JxMean + Jx;             	 JxVar = JxVar + Jx.^2;
        Emean = Emean + E;                   Evar = Evar + E.^2;
        rho1Mean = rho1Mean + rho1;          rho1Var = rho1Var + rho1.^2;
        rhoJxCov = rhoJxCov + rho.*Jx; 
        rhoECov = rhoECov + rho.*E; 
        JxECov = JxECov + Jx.*E; 
    end
    
    %% Periodically display diagnostics
    if( rem(iStep,(nSteps/20)) < 1 )
        timeElapsed = toc;
        fprintf('Finished %g of %g samples; about %g seconds to completion\n',...
            iStep,nSteps,timeElapsed/iStep*(nSteps-iStep)); 
    end
end % End of Main Loop
fprintf('Total time elapsed: %g seconds\n\n',toc)


%% Normalize statistics
rhoMean = rhoMean/nSamp;
rhoVar = rhoVar/nSamp - rhoMean.^2;
JxMean = JxMean/nSamp;
JxVar = JxVar/nSamp - JxMean.^2;
Emean = Emean/nSamp;
Evar = Evar/nSamp - Emean.^2;
rho1Mean = rho1Mean/nSamp;
rho1Var = rho1Var/nSamp - rho1Mean.^2;
rhoJxCov = rhoJxCov/nSamp - rhoMean.*JxMean;
rhoECov = rhoECov/nSamp - rhoMean.*Emean;
JxECov = JxECov/nSamp - JxMean.*Emean;

%% Verify conservation
rhoAvgFinal = mean(rho(:));  
JxAvgFinal  = mean(Jx(:));    
EavgFinal = mean(E(:));
fprintf('Conserved variables:     (Initial)  /  (Final)\n');
fprintf('Density            :  %12.6g  /  %g \n',rhoAvgInit,rhoAvgFinal);
fprintf('Momentum           :  %12.6g  /  %g \n',JxAvgInit,JxAvgFinal);
fprintf('Energy             :  %12.6g  /  %g \n',EavgInit,EavgFinal);

%% Compute theory values for variances and print results
zeta = 1 + (massR-1)^2/massR*CmassEq*(1-CmassEq);
Neq = (CmassEq*rhoEq/mass1 + (1-CmassEq)*rhoEq/mass0)*cellVol;
rhoVarTh = zeta*rhoEq^2/Neq;
rho1VarTh = (CmassEq + (1-CmassEq)/massR)*CmassEq*rhoEq^2/Neq;
JxVarTh = rhoVarTh*velEq^2 + rhoEq*k_B*tempEq/cellVol;
EvarTh = gamma/gammaM1^2*presEq^2/Neq + gamma/gammaM1*rhoEq*k_B*tempEq*velEq^2/cellVol + ...
    zeta/4*rhoEq^2*velEq^4/Neq;
rhoJxCovTh = rhoVarTh*velEq;
rhoECovTh = 1/gammaM1*rhoEq*k_B*tempEq/cellVol + zeta/2*rhoEq^2/Neq*velEq^2;
JxECovTh = gamma/gammaM1*rhoEq*k_B*tempEq/cellVol*velEq + zeta/2*rhoEq^2/Neq*velEq^3;

fprintf('\nVariances and Covariances\n');
fprintf('           (Measured)   : (Theory)     (Percent Error)\n');
fprintf('rho     =  %-13.6g: %-13.6g %-7.4g %%\n', mean(rhoVar(:)),  rhoVarTh,  100*(mean(rhoVar(:))  - rhoVarTh) /rhoVarTh );
fprintf('Jx      =  %-13.6g: %-13.6g %-7.4g %%\n', mean(JxVar(:)),   JxVarTh,   100*(mean(JxVar(:))   -JxVarTh)   /JxVarTh );
fprintf('Energy  =  %-13.6g: %-13.6g %-7.4g %%\n', mean(Evar(:)),    EvarTh,    100*(mean(Evar(:))    -EvarTh)    /EvarTh );
fprintf('rho1    =  %-13.6g: %-13.6g %-7.4g %%\n', mean(rho1Var(:)), rho1VarTh, 100*(mean(rho1Var(:)) -rho1VarTh) /rho1VarTh );
fprintf('Rho-Jx  =  %-13.6g: %-13.6g %-7.4g %%\n', mean(rhoJxCov(:)),rhoJxCovTh,100*(mean(rhoJxCov(:))-rhoJxCovTh)/rhoJxCovTh );
fprintf('Rho-E   =  %-13.6g: %-13.6g %-7.4g %%\n', mean(rhoECov(:)), rhoECovTh, 100*(mean(rhoECov(:)) -rhoECovTh) /rhoECovTh );
fprintf('Jx-E    =  %-13.6g: %-13.6g %-7.4g %%\n', mean(JxECov(:)),  JxECovTh,  100*(mean(JxECov(:))  -JxECovTh)  /JxECovTh );