kind = os.getenv("KIND")
if kind == nil then
  kind = 'lbm'
end

relaxation = os.getenv("RELAXATION")
if relaxation == nil then
  relaxation = 'bgk'
end

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

Re = 25
simulation_name = 'tgv3D'
scaling = 'diffusive'
printRuntimeInfo = false
logging = { level = 1 }

sim_control = {
  time_control = {
    max = { iter = 2048 },
    interval = { iter = 2048 },
  }
}

-- paramsters
twoPi   = 2.0 * math.pi
length  = twoPi
dx      = length / (2^level)
rho0LB  = 1.0
rhoPhy  = 1.0
nuPhy   = 1.0 / Re
viscPhy = nuPhy * rhoPhy
tD = 1.0 / nuPhy / 2

omega = 1.8
nuLB  = ( 2.0/omega - 1.0 ) / 6.0
dt  = nuLB / nuPhy * dx * dx
cs2 = 1./3.
p0 = rho0LB / rhoPhy * dx^2 / dt^2 * cs2

-- Analytical solution. origin is assumed to be 0 0 0
function pressure(x,y,z,t)
  p1 = math.cos(2*x) * math.cos(2*z) + 2.0 * math.cos(2*y)
  p2 = math.cos(2*y) * math.cos(2*z) + 2.0 * math.cos(2*x)
  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) * math.cos(y) * math.cos(z) * math.exp(-t/tD)
  return 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) * math.sin(y) * math.cos(z) * math.exp(-t/tD)
  return -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) * math.cos(y) * math.cos(z)
  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) * math.cos(y) * math.sin(z)
  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) * math.sin(y) * math.sin(z)
  return syz * math.exp(-t/tD)
end
function ic_Syz(x,y,z,t)
  return Syz(x,y,z,0)
end

mesh = {
  predefined = 'cube',
  origin = {0, 0, 0},
  length = length,
  refinementLevel = level,
}

physics = { dt = dt, rho0 = rhoPhy }

scheme = {
  identify = {
    label = kind..'_'..relaxation,
    kind = kind,
    relaxation = relaxation,
    layout = layout,
  }, -- identify table

  fluid = { omega = omega, rho0 = rhoPhy },

  -- Initial condition for each field
  initial_condition = {
    pressure  = ic_pressure,
    velocityX = ic_velocityX,
    velocityY = ic_velocityY,
    velocityZ = 0.0,
    Sxx = ic_Sxx,
    Syy = ic_Syy,
    Szz = 0.0,
    Sxy = 0.0,
    Syz = ic_Syz,
    Sxz = ic_Sxz,
  }, -- intitial condition end

} -- scheme table
