--This file contains common parameters

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

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

simulation_name = 'nonNewtonian'
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
R = height / 2

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

-- For Grid convergence,
-- keep Re = U_aver * Ly / viscosity to be constant
omega = 1.9
velMaxLB = ( 2 * Re / omega - Re ) / 4 / ly
tmaxLB  = 5000 * ( 4.^(level - 6 ))

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

if Re == 10 then
  velMaxPL = 0.00446
elseif Re == 20 then
  velMaxPL = 0.00892
end

-- 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
p0Phy = 0
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

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