-- This file contains important parameters
--
level = tonumber(os.getenv("LEVEL"))
if level == nil then
  level = 6
end
Re = 25
kind = 'lbm_incomp'
scaling = 'diffusive'

printRuntimeInfo = true
logging = { level = 10 }
-- null_device = "/dev/null"

-- some constant parameters
twoPi = 2.0 * math.pi

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

mesh = {
  predefined = 'cube',
  origin = {0, 0, 0},
  length = length,
  refinementLevel = level,
}
simulation_name = 'TGV_3D'
scheme_label = 'fluid'

---------------------------------------------------------------------
-- Simulation setup
relaxation = 'bgk'
layout     = 'd3q19'

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

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

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

-- 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,t)
  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,t)
  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,t)
  return velocityY(x,y,z,0)
end

-- velocity Z
function velocityZ(x,y,z,t)
  return 0
end
function ic_velocityZ(x,y,z,t)
  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,t)
  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,t)
  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,t)
  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,t)
  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,t)
  return Syz(x,y,z,0)
end
