----------------------- PLEASE READ THIS ---------------------------!!!

-- This input file is set up to run for regression check
-- Please make sure you DO NOT MODIFY AND PUSH it to the repository
                                                                              
--------------------------------------------------------------------!!!   



require "seeder"
-- General test case geometry and setup settings
kernelDim = 3
dimension = 2
Re     = 25
periodic = true

-- random fluctuation
randFluc = 0.00
-- enforce equal number of elements on all processes
balance = 'none'            

-- Scaling and initialization settings
scaling = 'acoustic'
verbose = false
startFromZero = true
iterateIC = false
tEnd = 0.52-- run up to which simulation time?
interpolation_method = 'quadratic'
relaxation = 'bgk'
spatial = true
label = 'TGV'
maxlevel = level+refinement
if spatial then
  relaxation = relaxation..'_omega_spatial'
end
if kernelDim == 3 then
  stencil = 'd3q19'
elseif kernelDim == 2 then
  stencil = 'd2q9'
else
  stencil = 'd1q3'
end

u0=1.
length = 2.*math.pi
viscosity = 1/Re
dx = length/2.^level 
u0LB=0.05
omega0=1.9

if scaling == 'acoustic' then
  uLB=u0LB
  dt = uLB/u0*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.

x0      = { 0.5*dx, 0.5*dx, 0.5*dx }
if periodic then
  x0[3] = 1.5*dx
end 
xProbe = { x0[1], 2.^(level-3)*dx+x0[2], x0[3] }
--xProbe = { x0[1], x0[2]+2*dx, x0[3] }

t = 0.
cs2LB = 1./3.

interval = tEnd/10.
minInt = 0.0

gamma = 0.
kx=2*math.pi/length
ky=2*math.pi/length
k=math.sqrt(kx^2+ky^2)
kL=kx
amplitudeX  =  2./math.sqrt(3.)*math.sin(gamma+2.*math.pi/3.)
amplitudeY  =  2./math.sqrt(3.)*math.sin(gamma-2.*math.pi/3.)
amplitudeZ  =  2./math.sqrt(3.)*math.sin(gamma)


  originX = origin[1]+ 0.25*length
  originY = origin[2]+ 0.25*length
  amplitude = 10.
  halfwidth = length/50.
  function pressure(x,y,z,t)
    return amplitude*math.exp(-0.5/(halfwidth^2)*(( x - originX )^2+( y - originY )^2))
    end
  function velocityX(x,y,z,t)
    return  0.
    end
  function velocityY(x,y,z,t)
    return 0.
    end

function refVelocity(x,y,z,t)
  return math.sqrt(velocityX(xProbe[1],xProbe[2],xProbe[3],t)^2 + velocityY(xProbe[1],xProbe[2],xProbe[3],t)^2)
  end
function refVelocityX(x,y,z,t)
  return velocityX(xProbe[1],xProbe[2],xProbe[3],t)
  end
if verbose then
  print('  evaluate press at probe '..pressure(xProbe[1], xProbe[2], xProbe[3], 0.))
  print('  evaluate dens at probe '..res)
  print('  evaluate velX  at probe '..velocityX(xProbe[1], xProbe[2], xProbe[3], 0.))
  print('  evaluate velY  at probe '..velocityY(xProbe[1], xProbe[2], xProbe[3], 0.))
end 

simulation_name = label
-- Defining the mesh can take two forms:
-- * simply an string specifying the directory to read the mesh from, if it is
--   indeed an directory there has to be a trailing "/", otherwise the string
--   is simply a prefix to all mesh files.
-- * a table for a predefined mesh. Right now only the full 'cube' without
--   any boundary conditions (periodic in all directions) is available as a
--   predefined mesh.
if refinement > 0 or periodic then
  mesh = 'mesh/'
else
  mesh = { predefined='cube',    -- use the predefined full cube
           origin = {0.,0.,0.},  -- origin of the cube
           length = length,      -- length of the cube
           refinementLevel = level } -- refinement level to resolve the cube
end 

-- Actual time definitions
-- times can be given in terms of iterations or simulation time.
time = {useIterations = false, -- Timings given in iterations? default: true
        max = tEnd,           -- Maximal iteration to reach
        interval = interval }    -- Interval for checking density

-- Debug options
ebug = { 
  verbosity = 100,
  debugMode = true,  -- activate global debug modus
  debugFiles = true, -- open the debugFiles for each process dbgOut0000**.out
  debugDependencies = true,
  dumpTreeIDs = true,
  dumpAuxLists = true,
  checkSteps = true,
  ebugRestart = true,
  dumpState = true,
--  debugDependencies = true,
--  dumpDependencies = true -- write the dependencies between levels to disk
}

-- outer omega cutoff ratio
w_min = 100.
-- inner omega cutoff ratio
w_max = 1.    
cutoff_rmin = 0.25
cutoff_rmax = 0.35
cutoff_min = cutoff_rmin*length*0.5
cutoff_max = cutoff_rmax*length*0.5


function spatialFunction(x,y,z,t)
   if(  y < cutoff_min) then
     res = w_max
   elseif( y  >= cutoff_max ) then
     res =  w_min
   else
     slope = (w_max-w_min)/(cutoff_min-cutoff_max)
     res =  y*slope + w_max - cutoff_min*slope
   end
   --print ('result'..res)
   return res
end


-- scheme model for single fluid simulation
scheme = {
  physics = { dt = dt, rho0 = rho0, setPhysicalQuantities = true },
  identify = {
      label = 'global',  -- scheme name 
      kind = 'lbm_incomp',      -- simulation type of this scheme
      relaxation = relaxation, -- relaxation type (bgk, mrt, ...)
      layout = stencil},
   fluid = { 
     visc_spatial = { 
       predefined = 'lua_fun',
       fun = spatialFunction,
       visc_spatial = spatialFunction, 
       store = true},
     omega = omega, 
     rho0 = rho0LB },
-- Initial condition for each field
initial_condition = { 
-- initial with lua function
     pressure   = pressure, --/rho0/dx^2*dt^2/cs2LB,
     velocityX = velocityX, --*dt/dx,
     velocityY = velocityY, --*dt/dx,
     velocityZ = 0.
     }

, 
tracking = {{
   label = 'press',
   variable = { 
    'pressure_phy', 
    'density', 
    'velocity', 
    'wss', 
    'pressure', 
              }, 
   folder = 'tracking/',
   shape = {kind = 'canoND', object = {origin = {0.,length*0.25, origin[3]+dx+shift+dxDash } }},
   format = 'ascii',
   time = {min = 0, max = tEnd, interval = minInt},
  }
  ,{
   label = 'hvs',
   variable = { 
     'pressure_phy', 
     'velocity_phy'
   }, -- variable table
   folder = 'tracking/',
   shape = {kind = 'all',
     object = { 
       origin = x0,
       vec = { { length, 0., 0. }, 
               { 0., length, 0. }
       }, segments = {2^maxlevel, 2^maxlevel}
     }
   },
   format = 'harvester',
   time = {min = 0, max = tEnd, interval = interval},
  }
--  ,{
--   label = label..dimension..'D_l'..level..'_'..relaxation..'_Y',
--   variable = { 
--     'pressure', 'velocity',
--     'pressure_phy', 
--     'velocity_phy'
--   }, -- variable table
--   folder = 'tracking/',
--   shape = {kind = 'canoND',
--     object = { 
--       origin = {0.,0.,0.},
--       vec = { { length, 0., 0. }, 
--               { 0., 0., length }
--       }, segments = {2^level, 2^level}
--     }
--   },
--   format = 'harvester',
--   time = {min = 0, max = tEnd, interval = interval},
--  }
--  ,{
--   label = label..dimension..'D_l'..level..'_'..relaxation..'_full',
--   variable = { 
--     'pressure', 'velocity',
--     'pressure_phy', 
--     'velocity_phy'
--   }, -- variable table
--   folder = 'tracking/',
--   shape = {kind = 'canoND',
--     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 = {min = 0, max = tEnd, interval = interval},
--  }
}

,
}

if walls then
  print("Walls")
 -- Boundary conditions
scheme["boundary_condition"] = {  
 { label = 'walls',
    kind = 'wall' }
 }
end

