! Copyright (c) 2014-2015, 2019-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2015-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

??!------------------------------------------------------------------------------
??! LB derive Macros
??!------------------------------------------------------------------------------
??! Formulation of the velocity calculation from PDF
??! You must provide
??! * an output array vel( 3 )
??! * an input quantity, pdf(stencil%QQ) - state array of an element
??! * an input quantity, densInv - inverse of density to convert mom to vel
??! * an input quantity, stencil in the mus_scheme_layout_type
?? text :: velocity(vel, pdf, densInv, stencil)
    ?vel?(1) = sum( ?pdf? * ?stencil?%cxDirRK(1, :) ) * ?densInv?
    ?vel?(2) = sum( ?pdf? * ?stencil?%cxDirRK(2, :) ) * ?densInv?
    ?vel?(3) = sum( ?pdf? * ?stencil?%cxDirRK(3, :) ) * ?densInv?
?? end text velocity
??!------------------------------------------------------------------------------


??!------------------------------------------------------------------------------
??! Incompressible formulation of the equilibrium for specific direction.
??! Returns a scalar.
??! You must provide
??! * an input quantity rho (density = sum(f_i))
??! * an input quantity vel(1:3) (velocity = sum(c_i*f_i))
??! * an input layout - mus_scheme_layout_type
??! * an input quantity iDir, specific direction
??! * a reference density rho0 (usually =1 )
??! * cs2inv from tem_param_module
??! * The selected real kind rk
??!------------------------------------------------------------------------------
?? text :: EqDir_incomp(eq, rho, vel, layout, iDir, rho0, cs2inv, rk)
  ?eq? = ?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.5_?rk?                        &
    &        - sum(?vel?(:) * ?vel?(:)) * 0.5_?rk? * ?cs2inv? ) )
?? end text EqDir_incomp
??!------------------------------------------------------------------------------


??!------------------------------------------------------------------------------
??! Compressible formulation of the equilibrium for specific direction
??! You must provide
??! * an output scalar fEq
??! * an input quantity rho (density = sum(f_i))
??! * an input quantity vel(1:3) (velocity = sum(c_i*f_i))
??! * an input layout - mus_scheme_layout_type
??! * an input quantity iDir, specific direction
??! It uses cs2inv = 3._rk from tem_param_module
??! and the _rk real kind for double precision
?? text :: EqDir_comp(fEq, rho, vel, layout, iDir, cs2inv, rk)
    ! calculate equilibrium density
    ?fEq? = ?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.5_?rk?                          &
       &        - sum(?vel?(:)*?vel?(:)) * 0.5_?rk? * ?cs2inv? ) )
?? end text EqDir_comp
??!------------------------------------------------------------------------------


??!-----------------------------------------------------------------------------
??! Formulation to derive deviatoric second-order moments from non-equilibrium
??! PDF. Same formulation can be used for both compressible and incompressible 
??! models
?? text :: devSecondMom(nEqTens, nEqPDF, layout)
   ! Now calculate the symmetric deviatoric second-order tensor of 
   ! nonEquilibrium part
   ! the static part cs2 I is usually neglected. 
   ! however, in current implementation it is considered
   ?nEqTens?(1) = sum( (?layout?%fStencil%cxcx(1,:) - cs2) * ?nEqPDF?)  !XX 
   ?nEqTens?(2) = sum( (?layout?%fStencil%cxcx(2,:) - cs2) * ?nEqPDF?)  !YY
   ?nEqTens?(3) = sum( (?layout?%fStencil%cxcx(3,:) - cs2) * ?nEqPDF?)  !ZZ
   ?nEqTens?(4) = sum( (?layout?%fStencil%cxcx(4,:) ) * ?nEqPDF?)  !XY
   ?nEqTens?(5) = sum( (?layout?%fStencil%cxcx(5,:) ) * ?nEqPDF?)  !YZ
   ?nEqTens?(6) = sum( (?layout?%fStencil%cxcx(6,:) ) * ?nEqPDF?)  !XZ
?? end text devSecondMom
??!-----------------------------------------------------------------------------
