-- This file contains important parameters
-- get mesh refine level
if os.getenv('LEVEL') == nil or os.getenv('LEVEL') == "" then
  level = 7
else
  level = tonumber( os.getenv('LEVEL'))
end

-- get Reynolds number
if os.getenv('Re') == nil or os.getenv('Re') == "" then
  Re = 25
else
  Re = tonumber( os.getenv('Re'))
end

-- get Reynolds number
kind = os.getenv('KIND')
if kind == nil or kind == "" then
  kind = 'lbm_incomp'
end

-- whether to visualize global domain
if os.getenv('GLOBAL') == "TRUE" then
  track_global = true
else
  track_global = false
end

-- get scaling method
if os.getenv('SCALING') == "DIFFUSIVE" or os.getenv('SCALING') == "diffusive" then
  scaling = 'diffusive'
else
  scaling = 'acoustic'
end

-- get interpolation method
if os.getenv('INTP') == "LINEAR" or os.getenv('INTP') == "linear" then
  intp  = 'linear'
else
  intp  = 'quadratic'
end

refinementLevel = 1
minLevel = level
minlevel = minLevel
maxLevel = level+refinementLevel

-- some constant parameters
sqrt3 = math.sqrt(3.0)
twoPi = 2.0 * math.pi

-- flow domain [0, 2*pi]
length = twoPi
dx     = length / (2^minLevel)
dxMin  = length / (2^maxLevel)
dxDash = 0.5 * dxMin

folder = 'mesh/'
mesh   = folder
simulation_name = 'TGV_3D'
scheme_label = 'fluid'

origin = {0., 0., 0.}      -- bounding box origin
refine_origin = { length/4 + dxDash, length/4 + dxDash, length/4 + dxDash }
refine_length = length / 2 - dxMin
refine_vec1 = { refine_length,           0.0,            0.0 }
refine_vec2 = {           0.0, refine_length,            0.0 }
refine_vec3 = {           0.0,           0.0,  refine_length }

uc = 0.00
theta = math.pi / 3.0
ucx = uc * math.cos(theta)
ucy = uc * math.sin(theta)

---------------------------------------------------------------------
-- Simulation setup
relaxation = 'bgk'
--relaxation = 'mrt'
layout     = 'd3q19'
interpolation_method = { method = intp }

-- Scaling and initialization settings
u0      = 1.0
rho0LB  = 1.0
rhoPhy  = 1.0
nuPhy   = 1.0 / Re
viscPhy = nuPhy * rhoPhy

-- if scaling == 'acoustic' then
--   -- keep velocity constant, compute omega from dt
--   uLB = 0.05
--   dt  = uLB / u0 * dx
--   nuLB = nuPhy * dt/dx/dx
--   omega = 1.0 / ( 3.0 * nuLB + 0.5)
-- else
-- keep omage constant, compute dt from omage
omega = 1.8
nuLB  = ( 2.0/omega - 1.0 ) / 6.0
dt  = nuLB / nuPhy * dx * dx
uLB = u0 * dt / dx

-- origin is the center of first element
x0 = { 0.5*dx, 0.5*dx, 0.5*dx }
x0 = { 0.0, 0.0, 0.0 }

xProbe = { x0[1] + math.pi/3, x0[2] + math.pi/3, x0[3] }

cs2 = 1./3.
p0 = rho0LB / rhoPhy * dx^2 / dt^2 * cs2

kx = twoPi / length -- wave number in x direction
ky = twoPi / length -- wave number in y direction
tD = 1.0 / nuPhy / (kx^2+ky^2) -- rate of decay coefficient

tEnd = math.ceil(tD/dx) * dx
tEnd = 2.0
interval = tEnd / 10.0

-- Analytical solution with respect to x0
function pressure(x,y,z,t)
  xc = x - x0[1]
  yc = y - x0[2]
  zc = z - x0[3]
  p1 = math.cos(2*xc) * math.cos(2*zc) + 2.0 * math.cos(2*yc)
  p2 = math.cos(2*yc) * math.cos(2*zc) + 2.0 * math.cos(2*xc)
  return (p0 + (p1+p2)/16.0) * math.exp(-2.0*t/tD)
end
function ic_pressure(x,y,z)
  return pressure(x,y,z,0)
end

-- velocity X
function velocityX(x,y,z,t)
  velX = math.sin(x-x0[1]) * math.cos(y-x0[2]) * math.cos(z-x0[3]) * math.exp(-t/tD)
  return u0 * velX
end
function ic_velocityX(x,y,z)
  return velocityX(x,y,z,0)
end

-- velocity Y
function velocityY(x,y,z,t)
  velY = math.cos(x-x0[1]) * math.sin(y-x0[2]) * math.cos(z-x0[3]) * math.exp(-t/tD)
  return u0 * -velY
end
function ic_velocityY(x,y,z)
  return velocityY(x,y,z,0)
end

-- velocity Z
function velocityZ(x,y,z,t)
  return 0
end
function ic_velocityZ(x,y,z)
  return velocityZ(x,y,z,0)
end

--velocity mag
function Velocity(x, y, z, t)
  velX = velocityX( x,y,z,t )
  velY = velocityY( x,y,z,t )
  velZ = velocityY( x,y,z,t )
  return math.sqrt( velX^2 + velY^2 + velZ^2 )
end

-- Sxx
function Sxx(x,y,z,t)
  sxx = math.cos(x-x0[1]) * math.cos(y-x0[2]) * math.cos(z-x0[3])
  return sxx * math.exp(-t/tD)
end
function ic_Sxx(x,y,z)
  return Sxx(x,y,z,0)
end
-- Syy
function Syy(x,y,z,t)
  return -Sxx(x,y,z,t)
end
function ic_Syy(x,y,z)
  return Syy(x,y,z,0)
end
-- Sxy
function Sxy(x,y,z,t)
  return 0.0  -- when kx == ky
end
function ic_Sxy(x,y,z)
  return Sxy(x,y,z,0)
end
-- Sxz
function Sxz(x,y,z,t)
  sxz = -0.5 * math.sin(x-x0[1]) * math.cos(y-x0[2]) * math.sin(z-x0[3])
  return sxz * math.exp(-t/tD)
end
function ic_Sxz(x,y,z)
  return Sxz(x,y,z,0)
end
-- Syz
function Syz(x,y,z,t)
  syz =  0.5 * math.cos(x-x0[1]) * math.sin(y-x0[2]) * math.sin(z-x0[3])
  return syz * math.exp(-t/tD)
end
function ic_Syz(x,y,z)
  return Syz(x,y,z,0)
end

-- Kinetic energy
function kinetic_energy( x, y, z, t )
  velX = velocityX( x,y,z,t )
  velY = velocityY( x,y,z,t )
  return 0.5 * ( velX^2 + velY^2 ) * rhoPhy * dx^3
end
