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

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

logging = { level = 3 }

minLevel = level
minlevel = level
maxLevel = level

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

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

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

---------------------------------------------------------------------
-- Simulation setup
printRuntimeInfo = false
kind       = 'lbm'
relaxation = 'bgk'
layout     = 'd3q19'
simulation_name = 'TGV_3D'
scheme_label = layout..'_'..relaxation

-- Scaling and initialization settings
u0      = 1.0
rho0LB  = 1.0
rhoPhy  = 1.0
nuPhy   = 1.0 / Re
omega = 1.98
nuLB  = ( 2.0/omega - 1.0 ) / 6.0
dt  = nuLB / nuPhy * dx * dx

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

cs2 = 1./3.
p0  = 0
tD  = 1.0 / nuPhy / 2.0

tEnd = 20
interval = tEnd / 10.0
tmax  = math.ceil(tEnd/dt)

-- Analytical solution with respect to x0
function pressure(x,y,z,t)
  p1 = math.cos(2*x) * (math.cos(2*z) + 2.0)
  p2 = math.cos(2*y) * (math.cos(2*z) + 2.0)
  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) * math.cos(y) * math.cos(z) * 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) * math.sin(y) * math.cos(z) * 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 = velocityZ( 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)
  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) * math.cos(y) * math.sin(z)
  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) * math.sin(y) * math.sin(z)
  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

sim_control = {
  time_control = {
    max = tEnd,
    interval = tEnd / 10,
  },
}

physics = { dt = dt, rho0 = rhoPhy }

identify = {
  label = scheme_label,
  kind = kind,
  relaxation = relaxation,
  layout = layout,
}

fluid = { omega = omega, rho0 = rhoPhy }

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

tracking = {
  { -- Average Kinetic energy (KE) spatial
   label = 'L'..level..'_ke',
   variable = { 'kinetic_energy_phy' },
   reduction = {'sum'},
   folder = 'tracking/',
   shape = { kind = 'all' },
   output = { format = 'ascii' },
   time_control = { min = 0, max = tEnd, interval = 1/100, },
  },

  -- { -- whole domain as harvester format
  --   label = 'L'..level..'_global',
  --   variable = {
  --     'pressure_phy',
  --     'velocity_phy',
  --     -- 'strain_rate_phy',
  --   }, -- variable table
  --   folder = 'output/',
  --   shape = { kind = 'all' },
  --   output = { format = 'vtk', },
  --   time_control = { min = tEnd, max = tEnd, interval = 1/20, },
  --   -- time_control = { min = 0, max = tEnd, interval = dt, },
  -- },
} -- track table
