%% 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 );