??!------------------------------------------------------------------------------
??! LB derive Macros
??!------------------------------------------------------------------------------
??! Incompressible formulation of the equilibrium 
??! You must provide 
??! * an output array fEq( layout%fStencil%QQ )
??! * an input quantity rho (density = sum(f_i))
??! * an input quantity vel(1:3) (velocity = sum(c_i*f_i))
??! * a stencil in the mus_scheme_layout_type
??! * a reference density rho0 (usually =1 )
??! It uses cs2inv = 3._rk from tem_param_module
??! and the _rk real kind for double precision
?? text :: Eq_incomp(fEq, rho, vel, layout, rho0)
    do iDir = 1, ?layout?%fStencil%QQ
      ! calculate equilibrium density
      ?fEq?( iDir ) = ?layout?%weight( iDir )                    &
        &           *( ?rho? +  ?rho0?*( cs2inv*                 &
        & sum( ?layout?%fStencil%cxDirRK(:,iDir)*?vel?(:))     &
        & + (sum( ?layout?%fStencil%cxDirRK(:,iDir)*?vel?(:))  &
        & *  sum( ?layout?%fStencil%cxDirRK(:,iDir)*?vel?(:))) &
        & *cs2inv*cs2inv*0.5d0                                   &
        & - sum( ?vel?(:)*?vel?(:))*0.5_rk*cs2inv ))
    enddo
?? end text Eq_incomp
??!------------------------------------------------------------------------------
??!
??!------------------------------------------------------------------------------
??! Compressible formulation of the equilibrium 
??! You must provide 
??! * an output array fEq( layout%fStencil%QQ )
??! * an input quantity rho (density = sum(f_i))
??! * an input quantity vel(1:3) (velocity = sum(c_i*f_i))
??! * a stencil in the mus_scheme_layout_type
??! It uses cs2inv = 3._rk from tem_param_module
??! and the _rk real kind for double precision
?? text :: Eq_comp(fEq, rho, vel, layout)
    do iDir = 1, ?layout?%fStencil%QQ
      ! calculate equilibrium density
      ?fEq?( iDir ) = ?layout?%weight( iDir )                    &
        &           * ?rho?*( 1._rk +  ( cs2inv*                 &
        & sum( ?layout?%fStencil%cxDirRK(:,iDir)*?vel?(:))     &
        & + (sum( ?layout?%fStencil%cxDirRK(:,iDir)*?vel?(:))  &
        & *  sum( ?layout?%fStencil%cxDirRK(:,iDir)*?vel?(:))) &
        & *cs2inv*cs2inv*0.5d0                                   &
        & - sum( ?vel?(:)*?vel?(:))*0.5_rk*cs2inv ))
    enddo
?? end text Eq_comp
??!------------------------------------------------------------------------------
??!
??!
??!------------------------------------------------------------------------------
??! Multispecies liquid mixture formulation of the equilibrium 
??! You must provide 
??! * an output array fEq( layout%fStencil%QQ )
??! * an input quantity rho(nFields) of all species
??! * an input quantity moleFrac(nFields) of all species
??! * an input quantity vel(3, nFields)
??! * a stencil in the mus_scheme_layout_type
??! * number of fields: nFields
??! * molecular weight ratio of all species: phi(nFields)
??! * resitivity coefficients of all species: resi_coeff(nFields, nFields)
??! * diffusivity omega: omega_diff 
??! * inverse of mixture%paramB: paramBInv
??! * current field ID
??! * theta_eq
??! * weight0Inv, inverse of lattice weight at rest position
??! It uses cs2inv = 3._rk, t2cs2inv, t2cs4inv from tem_param_module
?? text :: Eq_MSLiquid(fEq, iField, rho, moleFrac, vel, layout )
     
     eqVel(:) = ?vel?(:,?iField?)
     do ifld = 1, nFields
       eqVel(:) = eqVel(:) + resi_coeff(ifld,?iField?) * phi(?iField?) &
         &      * ?molefrac?(ifld) * ( ?vel?(:,ifld) - ?vel?(:,?iField?) ) &
         &      * paramBInv
     end do

     ! total mass density inverse
     totMass_densInv = 1.0_rk/sum(?rho?)

     ! mass averaged mixture velocity
     velAvg(1) = dot_product( ?rho?(:), ?vel?(1, :) )*totMass_densInv
     velAvg(2) = dot_product( ?rho?(:), ?vel?(2, :) )*totMass_densInv
     velAvg(3) = dot_product( ?rho?(:), ?vel?(3, :) )*totMass_densInv

     ! velocity in quadratic term of equilibrium
     velQuadTerm(:) = theta_eq*velAvg(:) + (1.0_rk-theta_eq)*eqVel(:)
     
     ! Calculate the square of velocity
     usq = dot_product( velQuadTerm,velQuadTerm ) * t2cs2inv

     do iDir = 1, ?layout?%fStencil%QQ
       ! Velocity times lattice unit velocity
       ucx = dot_product( ?layout?%fStencil%cxDirRK(:, iDir), eqVel )

       ucxQuadTerm = dot_product( ?layout?%fStencil%cxDirRK(:, iDir), &
         &                        velQuadTerm )

       ! calculate equilibrium
       ?fEq?(iDir) = ?layout?%weight( iDir )                          &
         &         * ?rho?(?iField?) * ( phi(?iField?)  + ucx * cs2inv &
         &         + ucxQuadTerm * ucxQuadTerm * t2cs4inv              &
         &         - usq )
     end do ! iDir
     
     ?fEq?(?layout?%stencil%restPosition) = &
       & ?layout?%weight( ?layout?%fStencil%restPosition ) * ?rho?(?iField?)&
       & * ( weight0Inv + (1.0_rk - weight0Inv)*phi(?iField?) - usq )
?? end text Eq_MSLiquid
??!------------------------------------------------------------------------------
??!
??!
??!------------------------------------------------------------------------------
??! Multispecies gas mixture (Asinari model) formulation of the equilibrium 
??! You must provide 
??! * an output array fEq( layout%fStencil%QQ )
??! * an input quantity rho(nFields) of all species
??! * an input quantity massFrac(nFields) of all species
??! * an input quantity vel(3, nFields)
??! * a stencil in the mus_scheme_layout_type
??! * number of fields: nFields
??! * resitivity coefficients of all species: resi_coeff(nFields, nFields)
??! * current field ID
??! * theta_eq
??! * weight0Inv, inverse of lattice weight at rest position
??! It uses cs2inv = 3._rk, t2cs2inv, t2cs4inv from tem_param_module
?? text :: Eq_MSGas(fEq, iField, rho, massFrac, vel, layout )
     
    eqVel(:) = ?vel?(:,?iField?)
    do ifld = 1, nFields
      eqVel(:) = eqVel(:) + chi(ifld, ?iField?)                       &
        &      * ?massfrac?(ifld) * ( ?vel?(:,ifld) - ?vel?(:,?iField?) ) 
    end do

     ! Calculate the square of velocity
     usq = dot_product( eqVel,eqVel ) * t2cs2inv

     do iDir = 1, ?layout?%fStencil%QQ
       ! Velocity times lattice unit velocity
       ucx = dot_product( ?layout?%fStencil%cxDirRK(:, iDir), eqVel )

       ! calculate equilibrium
       ?fEq?(iDir) = ?layout?%weight( iDir )                          &
         &         * ?rho?(?iField?) * ( phi(?iField?)  + ucx * cs2inv &
         &         + ucx * ucx * t2cs4inv - usq )
     end do ! iDir
     
     ?fEq?(?layout?%stencil%restPosition) = &
       & ?layout?%weight( ?layout?%fStencil%restPosition ) * ?rho?(?iField?)&
       & * ( weight0Inv + (1.0_rk - weight0Inv)*phi(?iField?) - usq )
?? end text Eq_MSGas
??!------------------------------------------------------------------------------
??!
??!
??!------------------------------------------------------------------------------
??! derive velocity of species of untransformed pdf from transformed pdf
??! by solving equation system of nFields.
??! You must provide
??! * an output array: momOut(3, nFields)
??! * an input quantity mole fraction of all species: moleFrac(nFields)
??! * first moments of transformed pdf momIn(3, nFields)
??! It uses:
??! * number of fields: nFields
??! * molecular weight ratio of all species: phi(nFields)
??! * resitivity coefficients of all species: resi_coeff(nFields, nFields)
??! * diffusivity omega: omega_diff 
??! * inverse of mixture%paramB: paramBInv
?? text :: momentumFromMacroLSE( momOut, moleFrac, momIn )
   ! build up the equation system for momentum
    matrixA = 0.0_rk
    do iField = 1, nFields
      ! set diagonal part
      matrixA( iField, iField ) = 1.0_rk
      do iField_2 = 1, nFields
        matrixA( iField, iField ) = matrixA( iField, iField )      &
          & + omega_diff * 0.5_rk * resi_coeff( iField, iField_2 ) &
          & * phi( iField ) * ?moleFrac?( iField_2 ) * paramBInv
      end do
      ! set nonDiagonal
      do iField_3 = 1, nFields
        matrixA( iField, iField_3 ) = matrixA( iField, iField_3 ) &
          & - omega_diff * resi_coeff( iField, iField_3 )         &
          & * phi( iField_3) * ?moleFrac?( iField ) * paramBInv
      end do
    end do
      
    invA = invert_matrix( matrixA )

    ! momentum of all species
    ?momOut?( 1, : ) = matmul( invA, ?momIn?( 1, : ) )
    ?momOut?( 2, : ) = matmul( invA, ?momIn?( 2, : ) )
    ?momOut?( 3, : ) = matmul( invA, ?momIn?( 3, : ) )

?? end text momentumFromMacroLSE
??!------------------------------------------------------------------------------
??!
??!
??!------------------------------------------------------------------------------
??! derive multispecies equilibrium velocity of given iField 
??! from molefraction and velocity of untransformed pdfs for liquid mixture
??! You must provide
??! * an output array: eqVel(3)
??! * Current field ID: iField
??! * an input quantity mole fraction of all species: moleFrac(nFields)
??! * velocity of untransformed pdf vel(3, nFields)
??! It uses:
??! * number of fields: nFields
??! * molecular weight ratio of all species: phi(nFields)
??! * resitivity coefficients of all species: resi_coeff(nFields, nFields)
??! * inverse of mixture%paramB: paramBInv
?? text :: eqVel_fromMacroMSLiquid( eqVel, iField, moleFrac, vel )
    ?eqVel?(:) = ?vel?(:,?iField?)
    do ifld = 1, nFields
      ?eqVel?(:) = ?eqVel?(:) + resi_coeff(ifld,?iField?) * phi(?iField?) &
        &        * ?molefrac?(ifld) * ( ?vel?(:,ifld) - ?vel?(:,?iField?) ) &
        &        * paramBInv
    end do
?? end text eqVel_fromMacroMSLiquid
??!------------------------------------------------------------------------------
??!
??!
??!------------------------------------------------------------------------------
??! derive multispecies equilibrium velocity of given iField 
??! from massfraction and velocity of untransformed pdfs for gas mixture
??! You must provide
??! * an output array: eqVel(3)
??! * Current field ID: iField
??! * an input quantity mass fraction of all species: massFrac(nFields)
??! * velocity of untransformed pdf vel(3, nFields)
??! It uses:
??! * number of fields: nFields
??! * compute chi matrix before. Look deriveEquilMSGas_FromState
?? text :: eqVel_fromMacroMSGas( eqVel, iField, massFrac, vel )
    ?eqVel?(:) = ?vel?(:,?iField?)
    do ifld = 1, nFields
      ?eqVel?(:) = ?eqVel?(:) + chi(ifld, ?iField?)                       &
        &        * ?massfrac?(ifld) * ( ?vel?(:,ifld) - ?vel?(:,?iField?) ) 
    end do
?? end text eqVel_fromMacroMSGas
??!------------------------------------------------------------------------------
