height =  0.41
width = height
nHeight = 16.
l_h = 3.0
nLength = nHeight*l_h 
dx = height/nHeight
length = nLength*dx
nLength_bnd = nLength+2
level = math.ceil(math.log(nLength_bnd)/math.log(2))
length_bnd = (2^level)*dx

dx_half = dx*0.5
zpos = dx_half
offset = 0.0 --height/2.0
dx_eps = length_bnd/2^10

qVal = false
all_outlet = false
side_allslip = false

--flow parameters
nu_phy = 1e-3 --m^2/s
rho0_p = 1.0
rho0_L = 1.0
cs2 = 1.0/3.0

Re = 20
u_mean_phy = Re*nu_phy/height
u_in_phy = 9.0*u_mean_phy/4.0

-- set true for acoustic scaling, false for diffusive scaling
acoustic_scaling = false

if acoustic_scaling == true then
--acoustic scaling
  u_in_L = 0.08
  dt = u_in_L*dx/u_in_phy
  u_mean_L = 2.0*u_in_L/3.0
  nu_L = nu_phy*dt/dx^2.
  omega = 1.0/(3.0*nu_L+0.5)
else
--diffusive scaling
  omega = 1.7
  nu_L = (1.0/omega-0.5)/3.0
  --nu_L = 0.0500
  --omega = 1.0/(3.0*nu_L+0.5)
  dt = nu_L*dx^2/nu_phy
  u_mean_L = u_mean_phy*dt/dx
  u_in_L = 3.0*u_mean_L/2.0
end

--press_drop = 8*u_in_phy*rho0_p*nu_phy*length/height^2
ref_press = 1.0 --( rho0_L*cs2 ) * ( rho0_p * (dx*dx) / (dt*dt) )
--press_grad = 1e-3
press_drop = 32*u_in_phy*rho0_p*nu_phy*length/height^2/height^2
press_grad = press_drop/length
press_inlet = ref_press+press_drop

cutoff_min = length/5.0
cutoff_max = length
-- max omega cut off ratio
w_max = 1.0
-- min omega cut off ratio
w_min = 10.0
function omega_spatial(x,y,z)
  if ( x < cutoff_min ) then
    res = w_max
  elseif (x >= cutoff_max)  then
    res = w_min
  else
    slope = (w_max-w_min)/(cutoff_min-cutoff_max)
    res =  x*slope + w_max - cutoff_min*slope
  end
  return res
end

-- Analytical solutions
--=====================

-- velocity
-------------------------------------------------------------------------------------------------------
-- function taken from "Benchmark Computations of Laminar Flow Around a Cylinder" by Schaefer and Turek
-------------------------------------------------------------------------------------------------------

function u_inflow_simple(x,y,z,t) 
  return (0.5/(rho0_p*nu_phy))*press_grad*y*z*(height-y)*(height-z)
 
 ----------------------
 -- for 'offset ~= 0.':
 ----------------------
 
 -- return (0.5/(rho0_p*nu_phy))*press_grad*(height*0.5-y)*(height*0.5+y)*(height*0.5-z)*(height*0.5+z)

end
--------------------------------------------------------------------------------------
-- function taken from "Shear stress in lattice Boltzmann simulations" by Timm Krueger
--------------------------------------------------------------------------------------

sum_1 = 0
for n = 1,14,2 do 
  a = n*math.pi*width/(2.0*height)
  sum_1 = sum_1 + ( 1.0/(n^3.) * math.sin(n*math.pi/2.0) / (( 1. + math.exp(-2.*a) ) / (2.*math.exp(-a))) )
end

E = (height^2.0)/4.0 - 8.0*(height^2.0)/(math.pi^3.0) * sum_1

function inner(y,z,o)
   return (1./(o*1.) * (math.sin(o*math.pi*y/width)*math.sin(1.*math.pi*z/height)) / 
       ( ((o/width)^2.) + ((1./height)^2.) )) + 
       (1./(o*3.) * (math.sin(o*math.pi*y/width)*math.sin(3.*math.pi*z/height)) / 
       ( ((o/width)^2.) + ((3./height)^2.) )) + 
       (1./(o*5.) * (math.sin(o*math.pi*y/width)*math.sin(5.*math.pi*z/height)) / 
       ( ((o/width)^2.) + ((5./height)^2.) )) + 
       (1./(o*7.) * (math.sin(o*math.pi*y/width)*math.sin(7.*math.pi*z/height)) / 
       ( ((o/width)^2.) + ((7./height)^2.) )) + 
       (1./(o*9.) * (math.sin(o*math.pi*y/width)*math.sin(9.*math.pi*z/height)) / 
       ( ((o/width)^2.) + ((9./height)^2.) )) + 
       (1./(o*11.) * (math.sin(o*math.pi*y/width)*math.sin(11.*math.pi*z/height)) / 
       ( ((o/width)^2.) + ((11./height)^2.) )) +
       (1./(o*13.) * (math.sin(o*math.pi*y/width)*math.sin(13.*math.pi*z/height)) / 
       ( ((o/width)^2.) + ((13./height)^2.) )) 
end    

function foo1(y,z,m)
    s = m*math.pi*(y-width/2.0)/height
    t = m*math.pi*width/(2.0*height)
    return 1.0/(m^3) * ((1. + math.exp(-2.*s)) / (2.*math.exp(-s))) / 
            ((1. + math.exp(-2.*t)) / (2.*math.exp(-t))) * math.sin(m*math.pi*z/height)
end    

function sum_2(y,z)
    return foo1(y,z,1) +foo1(y,z,3) +foo1(y,z,5) +foo1(y,z,7) +foo1(y,z,9) +foo1(y,z,11) +foo1(y,z,13) +foo1(y,z,15) 
end

function u_inflow_complex_wall(x,y,z,t)
    return u_in_phy/E * 32./(math.pi^4.) * ( inner(y,z,1) 
                                                  + inner(y,z,3) 
                                                  + inner(y,z,5)
                                                  + inner(y,z,7)
                                                  + inner(y,z,9)
                                                  + inner(y,z,11)
                                                  + inner(y,z,13) )
end

function u_inflow_complex_inner(x,y,z,t)
    return u_in_phy/E * ( z*(height-z) - 8.0*(height^2.0)/(math.pi^3)*sum_2(y,z) )
end
------------------------------------------------------------------------------------------------------------------
-- Analytical solutions
--function u_inflow(x,y,z,t) 
--  return (0.5/(rho0_p*nu_phy))*press_grad*(height*0.5-y)*(height*0.5+y)*(height*0.5-z)*(height*0.5+z)
--end

function wssA(x, y, z,t)
  return nu_phy * 2 * math.abs(y) * u_in_phy / (height*0.5)^2.0
end
function strainRateA(x, y, z, t)
  return -1*(y) * u_in_phy / R / R
end
function pressA(x,y,z, t)
  return ref_press+press_drop*(length-x)/length
end
-- X velocity
velocityX = u_inflow_simple
--            { predefined = 'combined',
--            temporal= {predefined='smooth', min_factor = 0.5, max_factor=1.0, 
--            from_time=0, to_time=1}, 
--            spatial = u_mean_phy
--            }
