--This file contains common parameters

--level = os.getenv("level")
--Reynolds = os.getenv("Re")
level = 8
Reynolds = 10
nNwtn = "carreau_yasuda"

diffusive = true  -- diffusive or acoustic scaling?
usePeriodic = true -- true = 2D channel, false = 3D channel
useObstacle = false -- not used yet
multi_level = 0   -- 0 = single-level, 1 = multi-level
printRuntimeInfo = false
scaling = 'diffusive'
identify_kind = 'lbm'

mesh = 'mesh/' -- Mesh files location

length = 8.0 * 0.01    -- physical length of channel
dx = length / 2.^level
height = length / 8

-- Do not change these values.
-- They are taken as a base to calculate the position of the periodic plane
ly = height / dx       -- number of elements for channel height
lx = length / dx - 2   -- number of elements for channel height
R = height / 2

-- For Grid convergence,
-- keep Re = U_aver * Ly / viscosity to be constant
omega = 1.9
velMaxLB = ( 2 * Reynolds / omega - Reynolds ) / 4 / ly
tmaxLB  = 1024 * 4

-- Physical parameters
densityPhy = 1050
viscPhy = 0.0035
nuPhy = viscPhy / densityPhy

-- nonNewtonian parameters
n = os.getenv("n")
n = 0.5
k = viscPhy
--velMaxPL = 0.001779
if Reynolds == 10 then
  velMaxPL = 0.00446
elseif Reynolds == 20 then
  velMaxPL = 0.00892
end
velMaxPL = 0.00892

-- reference LB quantities
rho0 = 1.
cs2  = 1./3.
p0LB = rho0*cs2
nuLB = ( 1.0 / omega - 0.5 ) / 3.0   -- viscosity
dt = nuLB * dx * dx / nuPhy
tmax = tmaxLB

velMaxPhy = velMaxLB * dx / dt
p0Phy = p0LB * densityPhy * dx * dx / dt / dt -- pressure at outlet
tmaxPhy = tmaxLB * dt

dp_dx = 2 * velMaxPhy * viscPhy / R / R
dp = dp_dx * length

-- Analytical solutions for Power Law nonNewtonian fluid
-- velocity for power law (PL)
function velPL(x, y, z)
  return velMaxPL * ( 1 - (math.abs(y)/R)^((n+1)/n) )
end

-- strain rate for power law
-- todo: add a sign function
function strainPL(x, y, z)
  return -0.5 * (velMaxPL/R) * ((n+1)/n) * ( (math.abs(y)/R)^(1/n) ) * math.abs(y)/y
end

-- Wall shear stress for power law
function wssPL(x, y, z)
  -- strainRate = strainPL(x, y, z)
  -- shearRate  = math.sqrt( 2.0 * strainRate * strainRate )
  -- visc = k * (shearRate^(n-1))
  -- return 2.0 * visc * strainRate
  return k * ( (2.0 * math.abs(strainPL(x,y,z)))^n )
end

-- shear Rate for power law
function shearRatePL(x, y, z)
  return 2.0 * math.abs( strainPL(x,y,z) )
end

-- analytical solutions for Newtonian Flow (just for reference)
-- Wall shear stress
function wssA(x, y, z)
  return viscPhy * 2 * math.abs(y) * velMaxPhy / R / R
end
-- strain rate
function strainRateA(x, y, z)
  return -y * velMaxPhy / R / R
end
-- shear rate
function shearRateA(x, y, z)
  return math.abs( 2.0 * strainRateA(x,y,z) )
end
-- velocity
function velA(x, y, z)
  return velMaxPhy * ( 1 - (y)^2 / R^2 )
end
-- pressure
function pressureA(x, y, z)
  return ( x - (-length/2) ) * dp_dx + p0Phy
end
