sqrt3 = math.sqrt(3.0)
twoPi = 2.0 * math.pi
level = 6
Re    = 25
kind  = 'lbm_incomp'
relaxation = 'mrt'
length = twoPi
dx     = length / (2^level)

-- Simulation setup
logging = { level = 3 }
mesh = {
  predefined = 'slice',
  origin = { 0.0, 0.0, 0.0 },
  length = length,
  refinementLevel = level,
}
simulation_name = 'L'..level
layout = 'd2q9'
printRuntimeInfo = false
io_buffer_size = 30

-- Test case parameters
origin = {0.0, 0.0, 0.0}      -- bounding box origin

uc    = 0.0
theta = math.pi / 3.0
ucx = uc * math.cos(theta)
ucy = uc * math.sin(theta)
kx = twoPi / length -- wave number in x direction
ky = twoPi / length -- wave number in y direction
---------------------------------------------------------------------
-- Scaling and initialization settings
u0      = 0.1
rho0LB  = 1.0
cs2     = 1.0 / 3.0
rhoPhy  = 1.0
nuPhy   = 1.0 / Re
viscPhy = nuPhy * rhoPhy

omega = 1.9
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*dx, 1.5*dx }

xProbe = { x0[1] + math.pi/3, x0[2] + math.pi/3, x0[3] }
--p0 = rho0LB / rhoPhy * dx^2 / dt^2 * cs2
p0 = 0
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-ucx*t))+math.cos(2.0*ky*(y-ucy*t))) * 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
function pressureLB(x,y,z,t)
  return pressure(x,y,z,0) * rhoPhy * dt^2 / dx^2
end

-- velocity X
function velocityX(x,y,z,t)
  velX = math.cos(kx*(x-ucx*t)) * math.sin(ky*(y-ucy*t))
  velX = velX * 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-ucx*t)) * math.cos(ky*(y-ucy*t))
  velY = velY * 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-ucx*t)) * math.sin(ky*(y-ucy*t)) * 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

-- Musubi config ---------------------------------------------------
sim_control = {
  time_control = { max = tEnd, interval = tEnd / 10, }
}

physics = { dt = dt, rho0 = rhoPhy }
identify = {
  kind = kind,
  relaxation = relaxation,
  layout = layout,
} -- identify table

fluid = { omega = omega }

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

variable = {
  -- new functions
  { name = 'pressure',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = {fun = pressure}
  },
  { name = 'ic_pressure',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = {fun = ic_pressure}
  },
  { name = 'pressureLB',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = {fun = pressureLB}
  },
  { name = 'velocityX',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = {fun = velocityX}
  },
  { name = 'ic_velocityX',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = {fun = ic_velocityX}
  },
  { name = 'velocityY',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = {fun = velocityY}
  },
  { name = 'ic_velocityY',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = {fun = ic_velocityY}
  },
  { name = 'velocityZ',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = 0.0
  },
  { name = 'ic_velocityZ',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = {fun = ic_velocityZ}
  },
  { name = 'velocity',
    ncomponents = 3,
    vartype = 'operation',
    operation = {
      kind = 'combine',
      input_varname = {
        'velocityX',
        'velocityY',
        'velocityZ'
      },
    },
  },
  { name = 'Sxx',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = {fun = Sxx}
  },
  { name = 'ic_Sxx',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = {fun = ic_Sxx}
  },
  { name = 'Syy',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = {fun = Syy}
  },
  { name = 'ic_Syy',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = {fun = ic_Syy}
  },
  { name = 'Sxy',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = 0.0
  },
  { name = 'ic_Sxy',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = {fun = ic_Sxy}
  },
  { name = 'strain_time',
    ncomponents = 3,
    vartype = 'st_fun',
    st_fun = {const = {0,0,0}}
  },
  { name = 'strain',
    ncomponents = 6,
    vartype = 'operation',
    operation = {
      kind = 'combine',
      input_varname = {'Sxx','Syy','Sxy','strain_time'}
    }
  },
  -- new functions
  -- analytical solution
  { name = 'refP',  ncomponents = 1, st_fun = pressure,  vartype = 'st_fun' },
  { name = 'refV',  ncomponents = 3, st_fun = velocity,  vartype = 'st_fun' },
  { name = 'refS',  ncomponents = 6, st_fun = strain,    vartype = 'st_fun' },
  -- errors
  { name = 'errorP',     ncomponents = 1, vartype = 'operation',
    operation = {
      kind = 'difference',
      input_varname = { 'pressure_phy', 'refP' },
    },
  },
  { name = 'errorV',     ncomponents = 3, vartype = 'operation',
    operation = {
      kind = 'difference',
      input_varname = { 'velocity_phy', 'refV' },
    },
  },
  { name = 'errorS',     ncomponents = 6, vartype = 'operation',
    operation = {
      kind = 'difference',
      input_varname = { 'strain_rate_phy', 'refS' },
    },
  },
}

tracking = {

  { -- l2norm error of velocity mag
    label = 'global_l2norm',
    variable = {
      'errorP', -- 2
      'refP',   -- 3
      'errorV', -- 4 5 6
      'refV',   -- 7 8 9
      'errorS', -- 10 11 12 13 14 15
      'refS',   -- 16 17 18 19 20 21
    },
    reduction = {'l2norm', 'l2norm', 'l2norm', 'l2norm', 'l2norm', 'l2norm'},
    folder = './',
    shape = { kind = 'all' },
    output = { format = 'ascii' },
    time_control = { min = tEnd, max = tEnd, interval = tEnd, },
  },
} -- track table
