--This file contains common parameters

level = tonumber(os.getenv("level"))
if level == nil then
  level = 8
end
Reynolds = tonumber(os.getenv("Re"))
if Reynolds == nil then
  Reynolds = 10
end


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 = ly / 2    -- channel radius
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  = 10000 * ( 4.^(level - 6 ))

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

-- nonNewtonian parameters
n = tonumber(os.getenv("n"))
if n == nil then
  n = 0.5
end
velMaxPL = 0.0019259
k = viscPhy

-- 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, t)
  return velMaxPL * ( 1 - (math.abs(y)/R)^((n+1)/n) )
end

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

function strainInit(x, y, z, t)
  return strainPL(x,y,z, t) * 2.0
end

-- Wall shear stress for power law
function wssPL(x, y, z, t)
  return k * ( (2.0 * math.abs(strainPL(x,y,z, t)))^n )
end

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