##  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.
##

# Main program: execute via the command "ruby unsteady.rb"

require 'lb_lattice.rb'
require 'lb_collision.rb'


# Simulation of a flow in a channel, past a cylindrical obstacle
class Channel
    
    # Constructor setup of the lattice, initialization of the constants
    def initialize(omega, u_max, l_x, l_y)
        @omega = omega
        @u_max = u_max
        @l_x   = l_x
        @l_y   = l_y
        
        @lattice =
            LB::Lattice.new( @l_x, @l_y, LB::Collision::LBGK.new(omega) )
    end #initialize

    # Initialize the boundaries, the obstacle and the initial condition
    def iniGeometry(obst_r, obst_x, obst_y)
        # Bounce-back on upper and lower wall
        @l_x.times { |x|
            @lattice.attribute_collision(
                x, 0, LB::Collision::BounceBack.instance )
            @lattice.attribute_collision(
                x, @l_y-1, LB::Collision::BounceBack.instance )
        }

        # Definition of the obstacle, and setup of initial values
        @l_x.times { |x|
            1.upto(@l_y-2){ |y|
                if (x-obst_x)**2 + (y-obst_y)**2 < obst_r**2
                    @lattice.attribute_collision(
                        x, y, LB::Collision::BounceBack.instance )
                else
                    u = [ computePoiseuille(y), 0.0 ]
                    @lattice.ini_equilibrium(x, y, 1, u)
                end
            }
        }

        # Poiseuille profile on the inlet and outlet. Currently, the 
        # only implemented boundary condition is the "Equilibrium BC"
        # at which the f's are set to their equilibrium distribution value.
        # A more accurate BC needs still  to be implemented.
        1.upto(@l_y-2) { |y|
            u = [ computePoiseuille(y), 0.0 ]
            cell = LB::Collision::EquilibriumBoundary.new(u, 1)
            @lattice.attribute_collision( 0, y, cell )
            @lattice.attribute_collision( @l_x-1, y, cell )
        }
    end #iniGeometry

    # Execute n_steps iteration steps (collision+streaming)
    def iterate(n_steps)
        n_steps.times { |t|
            puts t
            @lattice.step
        }
    end #iterate

    # Write velocity norm to a file
    def print_out(fname)
        File.open(fname, "w") { |file|
            @lattice.norm_u.each { |column| file.puts column.join(" ") }
        }
    end #print_out

    private

    # Computation of Poiseuille profile
    def computePoiseuille(y)
        y_   = y-0.5
        l_y_ = @l_y-2
        4 * @u_max / (l_y_*l_y_) * (l_y_*y_ - y_*y_)
    end

end #class Channel


# Main Program
##############

if __FILE__ == $0
  # Constants
  omega    = 1.5
  u_max    = 0.02
  n_steps  = 30
  l_x, l_y = 100, 40
  obst_r   = l_y/10+1
  obst_x   = l_x/5
  obst_y   = l_y/2

  # Simulation
  channel = Channel.new(omega, u_max, l_x, l_y)
  channel.iniGeometry(obst_r, obst_x, obst_y)
  channel.iterate(n_steps)

  # Output
  channel.print_out("u.dat")
end
