height =  0.40
nHeight = 32.0
l_h = 3.0
nLength = nHeight*l_h 
dx = height/nHeight
length = nLength*dx
--length = 2.2
nLength_bnd = nLength+2
--nLength = 200
--dx = length_bnd/(2^level)
--dx_ini = length/nLength 
level = math.ceil(math.log(nLength_bnd)/math.log(2))
length_bnd = (2^level)*dx
--level = 9
--dx = length_bnd/2^level
--nLength = math.ceil((length)/dx)
--nHeight = math.ceil((height)/dx)
dx_half = dx*0.5
zpos = dx_half
sph_pos = {0.2,0.2,zpos}
radius = 0.05
Dia = radius*2.0
offset = height/2.0

bc_type = 'BounceBack'
bc_type = 'moments'
if bc_type == 'moments' then
  dx_offset = dx
end
if bc_type == 'BounceBack' then
  dx_offset = 0.0
  nHeight = nHeight + 1  
end

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

Re = 1
u_in_phy = Re*nu_phy/height
print(u_in_phy)
u_mean_phy = 3.0*u_in_phy/2.0

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

function u_inflow(x,y,z,t) 
  --return 4.0*u_in_L*y*(height-y)/height^2.0
  --return (0.5/(rho0_p*nu_phy))*press_grad*y*(height-y)
  --return (0.5/(rho0_p*nu_phy))*press_grad*(height*0.5-y)*(height*0.5+y)
  --return 16.0 * u_in_phy * y * z * (height-y)*(height-z) / (height^2.0*height^2)
  --return 16.0*u_in_phy*(height*0.5-y)*(height*0.5+y)*(height*0.5-z)*(height*0.5+z)/(height^2.0*height^2)
  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
u_in_phy = u_inflow(0.0, 0.0,0.0)
u_mean_phy = 2.0*u_in_phy/3.0

Re = u_mean_phy*height/nu_phy

nL = math.ceil(Dia/dx)
Re_check = u_mean_phy*Dia/nu_phy
-- 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.5
  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

function omega_spatial(x,y,z,t)
  if x<1.0 then
    return x
  else
    return 1.0
  end
end

-- Analytical solutions
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
