!
!  MIXLBM.f90 
!
!  FUNCTIONS:
!	MIXLBM						- Main
!	UpdateLatticeData			- Update local distribution function f in P(i,j)
!	HydrodynamicMoments			- Compute hydrodynamic moments of f
!	EquilibriumDistribution		- Compute the equilibrium distribution function feq
!	D(M)						- Fick diffusivity
!	B(M1,M2)					- Stefan-Maxwell resistivity
!	writeDAT					- write results to files
!
!****************************************************************************
!
!  PROGRAM: MIXLBM
!
!  PURPOSE: To compare different implementations of species coupling 
!			(Fick model, Maxwell-Stefan) in LBM framework
!
!****************************************************************************

program MIXLBM

	implicit none

! Local variables: vectors
	real(8),dimension(:,:,:,:),allocatable :: fd,fd_new
	real(8),dimension(:,:),allocatable :: f
	real(8),dimension(:,:,:,:),allocatable :: rsigma,psigma,uxsigma,uysigma
	real(8),dimension(:,:),allocatable :: r,p,ux,uy
	real(8),dimension(:),allocatable :: phi
	real(8),dimension(:,:),allocatable :: BCsigma
	real(8),dimension(:,:),allocatable :: Output
	real(8),dimension(:),allocatable :: MM

	real(8),dimension(0:8) :: feq

! Local variables: scalars
	real(8) :: rho_new,ux_new,uy_new

	integer :: nx,ny,nt
	integer :: species,md
	integer :: i,j,k,s,t

! Values
	species = 3;
	nx = 100;
	ny = 10;
	nt = 30;

! Physical model
! 0- Ideal Fick model
! 1- Fick model
! 2- Maxwell-Stefan model
	md = 2;

! Allocate memory
	allocate( fd(1:nx,1:ny,0:8,1:species) )
	allocate( fd_new(1:nx,1:ny,0:8,1:species) )
	allocate( f(0:8,1:species) )
	allocate( rsigma(1:nx,1:ny,1:species,0:nt) )
	allocate( psigma(1:nx,1:ny,1:species,0:nt) )
	allocate( uxsigma(1:nx,1:ny,1:species,0:nt) )
	allocate( uysigma(1:nx,1:ny,1:species,0:nt) )
	allocate( r(1:nx,1:ny) )
	allocate( p(1:nx,1:ny) )
	allocate( ux(1:nx,1:ny) )
	allocate( uy(1:nx,1:ny) )
	allocate( phi(1:species) )
	allocate( BCsigma(1:4,1:species) )
	allocate( MM(1:species) )
	allocate( Output(1:nx,1:ny) )

! Molecular weight
	MM(1) = 1.0d0;
	MM(2) = 2.0d0;
	MM(3) = 3.0d0;		

! R0*T = 2.0 in lattice units, where R0 universal gas constant
	phi(1) = 1.0d0/MM(1);
	phi(2) = 1.0d0/MM(2);
	phi(3) = 1.0d0/MM(3);

! Boundary conditions
	BCsigma(:,1) = (/0.45d0,	-1.0d0,		0.55d0,	-1.0d0/);
	BCsigma(:,2) = (/0.50d0,	-1.0d0,		0.50d0,	-1.0d0/);
	BCsigma(:,3) = (/0.55d0,	-1.0d0,		0.45d0,	-1.0d0/);

!============================================================================
! Initialization
!============================================================================
	do i = 1,nx,1
		do j = 1,ny,1

!	Species 1
			if (i<int(real(nx)/2.0d0)) then
				psigma(i,j,1,0)  = 0.55d0;
			else
				psigma(i,j,1,0)  = 0.45d0;
			endif

!	Species 2
			psigma(i,j,2,0)  = 0.50d0;

!	Species 3
			if (i<int(real(nx)/2.0d0)) then
				psigma(i,j,3,0)  = 0.45d0;
			else
				psigma(i,j,3,0)  = 0.55d0;
			endif

!	Simple initialization of distribution functions (it may be inaccurate)
			do s = 1,species,1

				rsigma(i,j,s,0)  = psigma(i,j,s,0)*3.0d0/phi(s);
				uxsigma(i,j,s,0) = 0.0d0;
				uysigma(i,j,s,0) = 0.0d0;

				call EquilibriumDistribution(rsigma(i,j,s,0),phi(s),uxsigma(i,j,s,0),uysigma(i,j,s,0),feq)
				fd(i,j,:,s)		= feq(:);
				fd_new(i,j,:,s) = feq(:);

			enddo

		enddo
	enddo

!============================================================================
! Solver
!============================================================================
	write(*,'("Iter: "I"  Variable:"ES)') 0,psigma(4,4,3,t)

	do t = 1,nt,1
		do i = 1,nx,1
			do j = 1,ny,1

				call UpdateLatticeData(md,nx,ny,species,fd,i,j,phi,MM,BCsigma,f);

				do s = 1,species,1

					fd_new(i,j,:,s) = f(:,s);
					
					call HydrodynamicMoments(f(:,s),rsigma(i,j,s,t),uxsigma(i,j,s,t),uysigma(i,j,s,t))
					psigma(i,j,s,t) = phi(s)*rsigma(i,j,s,t)/3.0d0;

				enddo
			enddo
		enddo

		fd(:,:,:,:) = fd_new(:,:,:,:);

		write(*,'("Iter: "I"  Variable:"ES)') t, psigma(4,4,3,t)

	enddo

!============================================================================
! Output
!============================================================================
	call writeDAT(nx,ny,psigma(:,:,1,nt),'p1.csv');
	call writeDAT(nx,ny,psigma(:,:,2,nt),'p2.csv');
	call writeDAT(nx,ny,psigma(:,:,3,nt),'p3.csv');

	call writeDAT(nx,ny,rsigma(:,:,1,nt),'r1.csv');
	call writeDAT(nx,ny,rsigma(:,:,2,nt),'r2.csv');
	call writeDAT(nx,ny,rsigma(:,:,3,nt),'r3.csv');

	call writeDAT(nx,ny,uxsigma(:,:,1,nt),'u1.csv');
	call writeDAT(nx,ny,uxsigma(:,:,2,nt),'u2.csv');
	call writeDAT(nx,ny,uxsigma(:,:,3,nt),'u3.csv');

	deallocate( fd )
	deallocate( fd_new )
	deallocate( f )
	deallocate( rsigma )
	deallocate( psigma )
	deallocate( uxsigma )
	deallocate( uysigma )
	deallocate( r )
	deallocate( p )
	deallocate( ux )
	deallocate( uy )
	deallocate( phi )
	deallocate( BCsigma )
	deallocate( MM )
	deallocate( Output )

end program MIXLBM

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!============================================================================
! Subroutines / Functions
!============================================================================
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

subroutine HydrodynamicMoments(f,rho,ux,uy)

	implicit none

! Inlet / Outlet
	real(8),dimension(0:8),intent(in) :: f
	real(8),intent(out) :: rho,ux,uy

! Local variables
	real(8),dimension(0:8) :: zitax,zitay
	real(8) :: qx,qy
	integer :: i

! zitax of D2Q9 lattice
    zitax(:) = (/0.0d0, 1.0d0, 0.0d0,-1.0d0, 0.0d0, 1.0d0,-1.0d0,-1.0d0, 1.0d0/)
! zitay of D2Q9 lattice
    zitay(:) = (/0.0d0, 0.0d0, 1.0d0, 0.0d0,-1.0d0, 1.0d0, 1.0d0,-1.0d0,-1.0d0/)

	rho = 0.0d0;
	qx  = 0.0d0;
	qy  = 0.0d0;

	do i = 0,8,1
		rho = rho + f(i);
		qx  = qx  + zitax(i)*f(i);
		qy  = qy  + zitay(i)*f(i);
	enddo	

	ux = qx/rho;
	uy = qy/rho;

end subroutine HydrodynamicMoments

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

subroutine EquilibriumDistribution(rho,phi,ux,uy,feq)
!
! D2Q9 lattice

	implicit none

! Inlet / Outlet
	real(8),intent(in) :: rho,phi,ux,uy
	real(8),dimension(0:8),intent(out) :: feq

	feq(0) = 4.0d0/9.0d0*rho*(	(9.0d0-5.0d0*phi)/4.0d0							-3.0d0/2.0d0*(ux**2+uy**2) );

	feq(1) = 1.0d0/9.0d0*rho*(	phi	+3.0d0*(ux)		+9.0d0/2.0d0*(ux)**2		-3.0d0/2.0d0*(ux**2+uy**2) );
	feq(2) = 1.0d0/9.0d0*rho*(	phi	+3.0d0*(uy)		+9.0d0/2.0d0*(uy)**2		-3.0d0/2.0d0*(ux**2+uy**2) );
	feq(3) = 1.0d0/9.0d0*rho*(  phi	+3.0d0*(-ux)	+9.0d0/2.0d0*(-ux)**2		-3.0d0/2.0d0*(ux**2+uy**2) );
	feq(4) = 1.0d0/9.0d0*rho*(	phi	+3.0d0*(-uy)	+9.0d0/2.0d0*(-uy)**2		-3.0d0/2.0d0*(ux**2+uy**2) );

	feq(5) = 1.0d0/36.0d0*rho*(	phi+3.0d0*(ux+uy)	+9.0d0/2.0d0*(ux+uy)**2		-3.0d0/2.0d0*(ux**2+uy**2) );
	feq(6) = 1.0d0/36.0d0*rho*(	phi+3.0d0*(-ux+uy)	+9.0d0/2.0d0*(-ux+uy)**2	-3.0d0/2.0d0*(ux**2+uy**2) );
	feq(7) = 1.0d0/36.0d0*rho*(	phi+3.0d0*(-ux-uy)	+9.0d0/2.0d0*(-ux-uy)**2	-3.0d0/2.0d0*(ux**2+uy**2) );
	feq(8) = 1.0d0/36.0d0*rho*(	phi+3.0d0*(ux-uy)	+9.0d0/2.0d0*(ux-uy)**2		-3.0d0/2.0d0*(ux**2+uy**2) );

end subroutine EquilibriumDistribution

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

subroutine UpdateLatticeData(model,nx,ny,species,fd,i,j,phi,MMs,BCs,newfdPs)
!
! update distribution function of all the species in P(i,j)
! * for every direction k it identifies I(i+...,j+...)
! * if I is a point out of the domain, which border ?
! * if I is bulk, read f(k,s) in I and compute rsigma, ux sigma, uysigma in I
! * if I is border, then;
!	- read f(k,s) in P and assume rsigma, ux sigma, uysigma in I (wall / source)
!	- compute f(k,s) in I as feq(k,s)
! * compute barycentric r,p,ux,uy (Fick model)
! * compute barycentric rstar,pstar,uxstar,uystar (Maxwell-Stefan model)
! * compute local equilibrium feq(k) for species s
! * Collision Step: compute new value at the new time step f^+(k,s) in I
! * Streaming Step: take component BB(k) of f^+(k,s) in I and stream it in P
!	BB(k) is an operator which produce the direction opposite to direction k

	implicit none

! Internal Functions
	real(8) :: B,D

! Inlet / Outlet
	integer,intent(in) :: model,nx,ny,species
! distribution function fd(x,y,v)
	real(8),dimension(1:nx,1:ny,0:8,1:species),intent(in) :: fd
! generic point P(i,j)
	integer,intent(in) :: i,j
! rekionship between pressure and density
	real(8),dimension(1:species),intent(in) :: phi
! molecular weights
	real(8),dimension(1:species),intent(in) :: MMs
! boundary conditions:
! BC(i)<0 Neumann,	 wall
! BC(i)>0 Dirichlet, rho = BC(i)
! 1, 0, -1, 0
! 0, 1, 0, -1
	real(8),dimension(1:4,1:species),intent(in) :: BCs

! new distribution function in P(i,j)
	real(8),dimension(0:8,1:species),intent(out) :: newfdPs

! Local variables
	real(8),dimension(:),allocatable :: rsigma,psigma,uxsigma,uysigma
	real(8),dimension(:),allocatable :: uxstar,uystar
	real(8),dimension(:,:),allocatable :: f
	real(8),dimension(:),allocatable :: lambda

	real(8),dimension(0:8) :: feq,fcoll
	integer,dimension(0:8,1:2) :: Incr
	integer,dimension(0:8) :: BB

	real(8) :: r,p,ux,uy
	real(8) :: vx,vy
	real(8) :: temp,temp1,temp2
	real(8) :: MM
	integer :: k,s,vs,ik,iik,iI,jI
	integer :: BCdirection

! Allocate memory
	allocate( f(0:8,1:species) )
	allocate( rsigma(1:species) )
	allocate( psigma(1:species) )
	allocate( uxsigma(1:species) )
	allocate( uysigma(1:species) )
	allocate( lambda(1:species) )
	allocate( uxstar(1:species) )
	allocate( uystar(1:species) )

	Incr(:,1) = (/0, 1, 0, -1, 0, 1, -1, -1, 1/) ! x-axis
	Incr(:,2) = (/0, 0, 1, 0, -1, 1, 1, -1, -1/) ! y-axis

	BB(:) = (/0,3,4,1,2,7,8,5,6/) ! bounce-back

	do k=0,8,1
		iI = i + Incr(k,1);
		jI = j + Incr(k,2);

		do s=1,species,1

			rsigma(s)  = 0.0d0;
			psigma(s)  = 0.0d0;
			uxsigma(s) = 0.0d0;
			uysigma(s) = 0.0d0;

		! Boundary Conditions
			BCdirection = 0;		
			if (iI>nx) then
			! BC-1 
				BCdirection = 1;
 			elseif (jI>ny) then 
			! BC-2 
				BCdirection = 2;
			elseif (iI<1) then 
			! BC-3
				BCdirection = 3;
			elseif (jI<1) then 
			! BC-4
				BCdirection = 4;
			endif

			if (BCdirection == 0) then
			! bulk domain

				do ik=0,8,1
					f(ik,s) = fd(iI,jI,ik,s);
				enddo

				call HydrodynamicMoments(f(:,s),rsigma(s),uxsigma(s),uysigma(s))
				psigma(s)  = rsigma(s)*phi(s)/3.0d0;
			else
			! boundary layer

				if (BCs(BCdirection,s)<0) then
					! wall: specify ux, uy
					do ik=0,8,1
						f(ik,s) = fd(i,j,ik,s);
					enddo
					call HydrodynamicMoments(f(:,s),rsigma(s),uxsigma(s),uysigma(s))				
					uxsigma(s) = 0.0d0;
					uysigma(s) = 0.0d0;
					psigma(s)  = rsigma(s)*phi(s)/3.0d0;
					call EquilibriumDistribution(rsigma(s),phi(s),uxsigma(s),uysigma(s),f(:,s))
				else
					! source: specify rho
					do ik=0,8,1
						f(ik,s) = fd(i,j,ik,s);
					enddo
					call HydrodynamicMoments(f(:,s),rsigma(s),uxsigma(s),uysigma(s))				
					psigma(s) = BCs(BCdirection,s);				
					rsigma(s) = 3.0d0*psigma(s)/phi(s);
					call EquilibriumDistribution(rsigma(s),phi(s),uxsigma(s),uysigma(s),f(:,s))
				endif

			endif

		enddo

		r  = 0.0d0;
		p  = 0.0d0;

		do s=1,species,1
			r = r + rsigma(s);
			p = p + psigma(s);
		enddo

! Barycentric quantities
		ux = 0.0d0;
		uy = 0.0d0;

		do s=1,species,1
			ux = ux + ( rsigma(s)/r )*uxsigma(s);
			uy = uy + ( rsigma(s)/r )*uysigma(s);
		enddo

! Molar quantities
		vx = 0.0d0;
		vy = 0.0d0;

		do s=1,species,1
			vx = vx + ( psigma(s)/p )*uxsigma(s);
			vy = vy + ( psigma(s)/p )*uysigma(s);
		enddo

! Molecular weight for mixture
		temp = 0.0d0;

		do s=1,species,1
			temp = temp + ( rsigma(s)/r )/MMs(s);
		enddo
		MM = 1/temp;

! Modified quantities (Maxwell-Stefan)
		do s=1,species,1
			uxstar(s) = uxsigma(s);
			uystar(s) = uysigma(s);

			do vs=1,species,1
				uxstar(s) = uxstar(s) + (MM**2)/(MMs(s)*MMs(vs))* &
					B(MMs(s),MMs(vs))/B(MMs(s),MMs(s))* &
					( rsigma(vs)/r )*(uxsigma(vs)-uxsigma(s));
				uystar(s) = uystar(s) + (MM**2)/(MMs(s)*MMs(vs))* &
					B(MMs(s),MMs(vs))/B(MMs(s),MMs(s))* &
					( rsigma(vs)/r )*(uysigma(vs)-uysigma(s)); 	
			enddo

		enddo

		do s=1,species,1

! Select your model
			if (model==0) then
			! Ideal Fick model
				call EquilibriumDistribution(rsigma(s),phi(s),0.0d0,0.0d0,feq)
				lambda(s) = phi(s)/( 3.0d0*D(MMs(s)) );
			elseif (model==1) then
			! Fick model
				call EquilibriumDistribution(rsigma(s),phi(s),vx,vy,feq)
				lambda(s) = phi(s)/( 3.0d0*D(MMs(s)) );
			elseif (model==2) then
			! Maxwell-Stefan model
				call EquilibriumDistribution(rsigma(s),phi(s),uxstar(s),uystar(s),feq)
				lambda(s) = p*B(MMs(s),MMs(s))/r;
			endif

! Collision Step
			do ik=0,8,1
				fcoll(ik) = f(ik,s)+lambda(s)*( feq(ik)-f(ik,s) );
			enddo

! Streaming Step
! only one component in BB(k) direction is streamed
! BB(k) is an operator which produce the direction
! opposite to direction k
			newfdPs(BB(k),s) = fcoll(BB(k));

		enddo

	enddo

	deallocate( f )
	deallocate( rsigma )
	deallocate( psigma )
	deallocate( uxsigma )
	deallocate( uysigma )
	deallocate( lambda )
	deallocate( uxstar )
	deallocate( uystar )

end subroutine UpdateLatticeData

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

real(8) function B(m1,m2)
!
! Binary resistivity of Maxwell-Stefan
! for couple (sigma , varsigma) species
! ** it depends on both species !! **

	implicit none

! Inlet/Outlet
	real(8),intent(in) :: m1,m2

	B = (1.0d0/m1+1.0d0/m2)**(-0.5d0);

	return

end function B

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

real(8) function D(m)
!
! Diffusivity of Fick
! for ( sigma ) species
! ** it depends on one species !! **

	implicit none

! Inlet/Outlet
	real(8),intent(in) :: m

	D = 1.d0/m;

	return

end function D

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

subroutine writeDAT(nx,ny,M,name)

	implicit none

! Inlet
	integer,intent(in) :: nx,ny
	real(8),dimension(1:nx,1:ny),intent(in) :: M
	character(8),intent(in) :: name

! Local variables
	integer :: i,j
	integer :: IO_Stat_Number1=-1

	open(UNIT=1,FILE=name,IOSTAT=IO_Stat_Number1)

	if (IO_Stat_Number1==0) then
		print *,'Saving ',name,' ...'

		do j=ny,1,-1
			do i=1,nx,1
				if (i<nx) then
					write(1,'($F",")') M(i,j)
				else
					write(1,'($F)') M(i,j)
				endif
			enddo
			write(1,*)
		enddo

		print *, '...done' 
	endif
	
	close(UNIT=1)

end subroutine writeDAT

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&