##  Lattice Boltzmann sample, written in Ruby
##
##  Main author: Jean-Luc Falcone
##  Co-author: Jonas Latt
##  Copyright (C) 2006 University of Geneva
##  Address: Jean-Luc Falcone, Rue General Dufour 24,
##           1211 Geneva 4, Switzerland 
##
##  This program is free software; you can redistribute it and/or
##  modify it under the terms of the GNU General Public License
##  as published by the Free Software Foundation; either version 2
##  of the License, or (at your option) any later version.
##
##  This program is distributed in the hope that it will be useful,
##  but WITHOUT ANY WARRANTY; without even the implied warranty of
##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##  GNU General Public License for more details.
##
##  You should have received a copy of the GNU General Public 
##  License along with this program; if not, write to the Free 
##  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
##  Boston, MA  02110-1301, USA.
##

require 'singleton.rb'
require 'lb_util.rb'

module LB

# Various collision steps are implemented, depending on the dynamics of a
# cell: BGK, bounce-back, equilibrium value, ...
module Collision
  
    # BGK bulk dynamics
    class LBGK
      
        # Constructor: initialize the local relaxation parameter, related
        # to the fluid viscosity
        def initialize(omega)
            @omega = omega
            @omega_complement = 1 - omega
        end #initialize(omega)
        
        # Computation of the local density
        def rho(f)
            sum = 0
            f.each { |val| sum += val }
            sum
        end #def rho(f)
        
        # Computation of the local flow velocity
        def u(f,rho=rho(f))
            u0 = u1 = 0
            f.each_with_index { |f_i,i|
                u0 += C[i][0]*f_i
                u1 += C[i][1]*f_i
            }
            [u0/rho, u1/rho]
        end #def u(f)
        
        # Computation of the local equilibrium
        def f_eq(f,i)
            rho = rho(f)
            u = u(f,rho)
            equilibrium(i, rho, u)
        end #f_eq(f,i)

        # Local collision step
        def update(f)
            f_prime = Array.new
            f.each_with_index {|f_i,i|
                f_prime << @omega_complement*f_i + @omega*f_eq(f,i)
            }
            f_prime   
        end #update(f)

        private

        include D2Q9

    end #class LBGK

    # Bounce-back for no-slip boundary conditions
    class BounceBack

	 include D2Q9
        
        # The object can be used as a singleton, as it has no
        # internal state
        include Singleton
          
        # Computation of the local density (dummy value)
        def rho(f)
            1.0
        end 

        # Computation of the local flow velocity (dummy value)
        def u(f)
            [0.0,0.0]
        end

        # Local collision step
        def update(f)
             f_prime = Array.new
             f.each_index {|i|  
                 f_prime << f[Opposite[i]]              
             } 
             f_prime
        end #update(f)

    end #class BounceBack
      
    # Equilibrium distribution for Dirichlet velocity boundaries.
    # This BC is only accurate at zeroth order in the hierarchy
    # of the Chapman-Enskog expansion
    class EquilibriumBoundary

	include D2Q9

        # Constructor: initialize the value for velocity and density
        def initialize(u0, rho0)
            @u0 = u0
            @rho0 = rho0
            @f0 = compute_f0 
        end #initialize(u0, rho0)
          
        # Computation of the local density
        def rho(f)
            @rho0
        end #def rho(f)
        
        # Computation of the local flow velocity
        def u(f)
            @u0
        end #def u(f)

        # Computation of the local equilibrium
        def f_eq(f,i)
            equilibrium(i, @rho0, @u0)
        end #f_eq(f,i)

        # Local collision step
        def update(f)
            @f0.dup
        end
        
        private
        
        def compute_f0
            f0 = []
            0.upto(8) { |i|
              f0 << f_eq(nil,i)
            }  
            f0
        end

    end #class EquilibriumBoundary
      
end # module Collision

end # module LB
