-- This is the double-shear layer test case, commonly used as in 
-- [1]  P. Dellar, "Bulk and shear viscosities in lattice Boltzmann equations," Phys. Rev. E, vol. 64, no. 3, p. 031203, Aug. 2001.
-- [2]  D. Ricot, S. Marie, P. Sagaut, and C. Bailly, "Lattice Boltzmann method with selective viscosity filter," J. Comp. Phys., vol. 228, pp. 4478–4490, 2009.
require 'seeder'
-- General test case geometry and setup settings
kernelDim = 3

nIntervals = 10
intp = 'quadratic'
intp = 'linear'
scaling = 'diffusive'
scaling = 'acoustic'
relaxation = 'mrt'
model='lbm'
model='lbm_incomp'
-- random fluctuation
perturb = 0.05
w = 0.05

-- Scaling and initialization settings
verbose = false
interpolation_method = { method = intp }
tEnd = 0.1   -- run up to which simulation time?
label = 'shearlayer'
if kernelDim == 3 then
  stencil = 'd3q19'
elseif kernelDim == 2 then
  stencil = 'd2q9'
else
  stencil = 'd1q3'
end

length = 1.
viscosity = 0.0001 
dx = length/2.^level 
u0LB=0.05
omega0=1.9
csLB = 1./math.sqrt(3.)
cs0 = 343
Ma = 0.04
u0 = Ma*cs0

if scaling == 'acoustic' then
  dt = csLB/cs0*dx
  viscLB = viscosity*dt/dx/dx
  omega = 1./(3.*viscLB + 0.5)
else
  omega = omega0
  viscLB  = 1./6.*(2./omega - 1.)
  dt = viscLB/viscosity*dx*dx  
  uLB=u0*dt/dx
end 
if verbose then
  print('Scaling  '..scaling..' on level '..level)
  print('  Re       '..Re)
  print('  omega    '..omega)
  print('  uLB      '..uLB)
  print('  dx       '..dx)
  print('  dt       '..dt)
end 

rho0LB  = 1.
rho0    = 1.
if model == 'lbm_incomp' then
  p0 = 0.
else
  p0 = rho0*cs0*cs0
end

cs2LB = 1./3.

interval = tEnd/nIntervals
minInt = 0.0

  function velocityX(x,y,z,t)
    if( y <= 0.5 ) then
      velX=u0*math.tanh(4.*(y-0.25)/w)
    else
      velX=u0*math.tanh(4.*(0.75-y)/w)
    end 
    return  velX
    end
  function velocityY(x,y,z,t)
    return u0*perturb*math.sin(2.*math.pi*(x+0.25))
    end


simulation_name = label
mesh = 'mesh/'

-- Actual time definitions
-- times can be given in terms of iterations or simulation time.
sim_control = {
  time_control = {
    max = tEnd,
    interval = interval},
  }


-- Debug options
ebug = { 
  debugMode = true,  -- activate global debug modus
  debugFiles = true, -- open the debugFiles for each process dbgOut0000**.out
  umpTreeIDs = true,
  umpAuxLists = true,
  heckSteps = true,
  debugRestart = true,
  umpState = true,
--  debugDependencies = true,
--  dumpDependencies = true -- write the dependencies between levels to disk
}

-- scheme model for single fluid simulation
scheme = {
  physics = { dt = dt, rho0 = rho0, setPhysicalQuantities = true },
  identify = {
      label = 'global',  -- scheme name 
      kind = model,      -- simulation type of this scheme
      relaxation = relaxation, -- relaxation type (bgk, mrt, ...)
      layout = stencil},
   fluid = { omega = omega, 
     rho0 = rho0LB },
-- Initial condition for each field
initial_condition = { 
-- initial with lua function
     pressure  = p0,--/rho0/dx^2*dt^2/cs2LB,
     velocityX = velocityX, --*dt/dx,
     velocityY = velocityY, --*dt/dx,
     velocityZ = 0.
--      Sxx = ic_Sxx,
--      Syy = ic_Syy,
--      Szz = ic_Szz,
--      Sxy = ic_Sxy,
--      Syz = ic_Syz,
--      Sxz = ic_Sxz
     }

, 
tracking = {{
   label = label..dimension..'D_l'..level..'_'..relaxation..'_press',
   variable = { 
    'pressure_phy', 
    'density', 
    'velocity', 
    'wss', 
    'pressure', 
              }, 
   folder = 'tracking/',
   shape = {kind = 'canoND', object = {origin = {0.5,0.5,dx/2.} } },
   format = 'ascii',
   time_control = {min = 1, max = tEnd, interval = minInt},
  },
  
  {
   label = label..dimension..'D_l'..level..'_'..relaxation..'_full',
   variable = { 
     'pressure', 'velocity',
     'pressure_phy', 
     'velocity_phy'
   }, -- variable table
   folder = 'tracking/',
   shape = {kind = 'all',
     object = { 
       origin = {0.,0.,0.},
       vec = { { length*0.5, 0., 0. }, 
               { 0., 0., length*0.5 }, 
               { 0., length*0.5, 0. }
       }, segments = {2^level, 2^level, 2^level}
     }
   },
   format = 'harvester',
   time_control = {min=1, max =tEnd, interval = interval},
  }
}
 
}


