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

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

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

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

logging = { level = 1 }

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

-- flow domain [0, 2*pi]
length = twoPi
dx     = length / (2^level)

-- folder = 'mesh/'
-- mesh = folder
mesh = {
  predefined = 'slice',
  origin = { 0, 0, 0 },
  length = length,
  refinementLevel = level,
}
simulation_name = 'TGV_2D'
printRuntimeInfo = false

origin = {0.0, 0.0, 0.0}      -- bounding box origin

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

---------------------------------------------------------------------
-- Simulation setup
io_buffer_size = 32

layout     = 'd2q9'
layout     = 'd3q19'
relaxation = 'bgk'
scheme_label = layout..'_'..relaxation

-- 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.6
  nuLB  = ( 2.0/omega - 1.0 ) / 6.0
  dt  = nuLB / nuPhy * dx * dx
  uLB = u0 * dt / dx
end

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

cs2 = 1.0 / 3.0
p0 = rho0LB / rhoPhy * dx^2 / dt^2 * cs2
if kind == 'lbm_incomp' then
  p0 = 0
end

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/dt) * dt
interval = tEnd / 10.0

-- Analytical solution with respect to x0
function pressure(x,y,z,t)
  -- when kx == ky
  p = (math.cos(2.0*kx*(x-x0[1]))+math.cos(2.0*ky*(y-x0[2]))) * math.exp(-2.0*t/tD)
  return p0 - u0^2 * p * rhoPhy / 4.0
end
function ic_pressure(x,y,z)
  return pressure(x,y,z,0)
end

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

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

-- velocity Z
function ic_velocityZ(x,y,z)
  return 0.0
end

-- velocity
function velocity(x, y, z, t)
  return { velocityX( x,y,z,t ), velocityY( x,y,z,t ), 0.0 }
end

-- Sxx
function Sxx(x,y,z,t)
  sxx = math.sin(kx*(x-x0[1])) * math.sin(ky*(y-x0[2])) * math.exp(-t/tD)
  return u0 * sxx
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
-- Strain
function strain( x,y,z,t )
  return { Sxx(x,y,z,t), Syy(x,y,z,t), Sxy(x,y,z,t), 0, 0, 0 }
end
