--
-- G. I. Taylor and a. E. Green. 
-- Mechanism of the Production of Small Eddies from Large Ones. 
-- Proceedings of the Royal Society A: Mathematical, 
-- Physical and Engineering Sciences, 158(895):499–521, February 1937.
dimension = 3
Re = 10.
omega = 1.8

-- Some parameters for the functions given below
-- Domain Size: L = 2pi
length = 2.*math.pi
label = 'TGV'
level  = 4
relaxation = 'mrt'
if dimension == 3 then
  stencil = 'd3q19'
else
  stencil = 'd2q9'
end

viscosity = 1./Re
rho0LB  = 1.
rho0    = 1.
viscLB  = 1./6.*(2./omega - 1.)
dx = getdxFromLevel( {len_bnd=length, level=level})
dt = getdtFromVisc( {dx = dx, nu_p = viscosity, nu_l = viscLB} )
print ('time step dt'..dt)
xProbe = { 1.518436,3.089233,0.}

t = 0.
cs2LB = 1./3.
-- Time step settings
tmax =  4000 -- some helping variable to use below
interval = tmax/10

-- size factor: chosen to be 0. in Brachet et al.
-- M. Brachet, D. Meiron, S. Orszag, B. Nickel, R. Morf, and U. Frisch, 
-- Small-scale structure of the taylor-green vortex, J. Fluid. Mech., 130 (1983), pp. 411–452.
-- in taylor&Green's definition there is this spatial factor gamma, which we always set to 0.
gamma = 0.
u0=1.
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)
print('dt/dx'..dt/dx)
print('dx   '..dx   )
print('dt   '..dt   )

function pressure(x,y,z,t)
-- for domain -pi < x pi
--  return -0.25*u0^2*math.cos(kL*(x+y))*math.cos(kL*(x-y))*math.exp(-2.*k^2*viscosity*t)
-- for domain 0 < x 2pi
  return 0.25*u0^2*(math.cos(2.*x)+math.cos(2.*y))*math.exp(-2.*k^2*viscosity*t)
  end
function velocityX(x,y,z,t)
  return  u0*math.sin(x)*math.cos(y)*math.exp(-k^2*viscosity*t)
  end
function velocityY(x,y,z,t)
  return -u0*math.cos(x)*math.sin(y)*math.exp(-k^2*viscosity*t)
  end

-- Initial conditions for the Taylor Green Vortex
--
if true then --dimension == 2 then -- 2D
  function ic_velocityX(x,y,z) -- initial conditions given in physical units, so convert to LB
    return  velocityX(x,y,z,0.)*dt/dx
    end
  function ic_velocityY(x,y,z) -- initial conditions given in physical units, so convert to LB
    return  velocityY(x,y,z,0.)*dt/dx
    end
  function ic_velocityZ(x,y,z) -- initial conditions given in physical units, so convert to LB
    return 0.
    end
  function ic_pressFluc(x,y,z) -- initial conditions given in physical units, so convert to LB
--    res = (rho0/4.*(math.cos(2.*2.*math.pi/length*x) + math.cos(2.*2.*math.pi/length*y)))*math.exp(-4.*4.*math.pi^2/length^2*viscosity*t)
--    return res
    return pressure(x,y,z,0.)
    end
else  -- 3D
  -- 3d initialization taken from p. 192
  -- [1] D. Fauconnier, “Development of a Dynamic Finite Difference Method for Large-Eddy Simulation,” pp. 1–296, Nov. 2008.
  function ic_velocityX(x,y,z) -- initial conditions given in physical units, so convert to LB
    res = amplitudeX*math.sin(x)*math.cos(y)*math.cos(z)*dt/dx
    return res
  end
  function ic_velocityY(x,y,z) -- initial conditions given in physical units, so convert to LB
    res =  amplitudeY*math.cos(x)*math.sin(y)*math.cos(z)*dt/dx
    return res
  end
  function ic_velocityZ(x,y,z) -- initial conditions given in physical units, so convert to LB
    return 0. --amplitudeZ*math.cos(x)*math.cos(y)*math.sin(z)*dt/dx
  end
  function ic_pressFluc(x,y,z) -- initial conditions given in physical units, so convert to LB
    return ((1-math.cos(2.*gamma ))/24.*(
    math.cos(2.*x)*math.cos(2.*y) + 2.*math.cos(2.*z)
    ) 
    + (2.+math.cos(2.*gamma)+math.sqrt(3.)*math.sin(2.*gamma))/48.*(
    math.cos(2.*x)*math.cos(2.*z) + 2.*math.cos(2.*y)
    )
    + (2.+math.cos(2.*gamma)-math.sqrt(3.)*math.sin(2.*gamma))/48.*(
    math.cos(2.*y)*math.cos(2.*z) + 2.*math.cos(2.*x)
    )
    )
  end
end
function ic_density(x,y,z)
  return ic_pressFluc(x,y,z)/rho0/dx^2*dt^2/cs2LB + rho0
  end

function refPressure(x,y,z,t)
  return pressure(x,y,z,t)+rho0LB*cs2LB*dx*dx/dt/dt
  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.
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

-- Define the interpolation method to use between different refinement levels
-- available are: average, eqneq, quadratic, quadraticCompact
interpolation_method = 'average' -- default: average

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

control_routine = 'multischeme'  -- another option is 'benchmark'

-- Debug options
debug = { debugMode = true,  -- activate global debug modus
          debugFiles = true, -- open the debugFiles for each process dbgOut0000**.out
          dumpDependencies = true -- write the dependencies between levels to disk
          }
function test(x,y,z,t)
  return 1.
  end
-- scheme model for single fluid simulation
scheme = {
  physics = { dt = dt, rho0 = rho0 },
  identify = {
      label = 'global',  -- scheme name 
      kind = 'lbm_incomp',      -- simulation type of this scheme
      relaxation = relaxation, -- relaxation type (bgk, mrt, ...)
-- Scheme layout
-- This describes the stencil to use in the simulation.
      layout = stencil},
  -- field which defines fluid or species
  -- Single fluid simulation
--  field = {
--     label = 'fluid',
-- Define the Lattice Boltzmann fluid properties
     fluid = { omega = omega, rho0 = rho0LB },
-- Initial condition for each field
-- Define initial conditions per variable.
-- They might be constants, predefined functions (described with a table), or
-- lua functions.
initial_condition = { 
-- initial with lua function
     density   = ic_density, --/rho0/dx^2*dt^2/cs2LB,
     velocityX = ic_velocityX, --*dt/dx,
     velocityY = ic_velocityY, --*dt/dx,
     velocityZ = ic_velocityZ 
     }

--,add_variable = {{'spacetime',1,ic_1Dgauss_pulse2}}
, add_variable = { 
    { 'refPressure', 1, test },
--    { name = 'velocityX', ncomponents = 1, spacetime = velocityX },
  }
, tracking = {{
   label = label..dimension..'D_'..relaxation..'_point',
   variable = { --{'pressure'}, {'velocity'},-- {'kinetic_energy_phy'},
--    { 'spacetime' },
    { name = 'difference', ncomponents = 1, dep = {'refPressure','pressure'} }
--    { name = 'difference', ncomponents = 1, dep = {'velocity_phy','velocityX'} },
-- possible names are: pdfs, density, velocity, velMag,
--                     equilFromState, equilFromMacro, shearstress, wss,
--                     spacetime, difference
--    {'velocity',sub_index={1}}
              }, -- variable table
   folder = 'tracking/',
   shape = {kind = 'canoND', object = {origin = xProbe } },
   format = 'ascii',
   time = {min = 0, max = tmax, interval = 2},
  }
  ,  
{
   label = label..dimension..'D_'..relaxation..'_average',
   variable = { {'density'},{'pressure_phy'}, {'kinetic_energy_phy'},
                {'velMag_phy'},{'velocity_phy'},{'velMag'},{'velocity'}
--    { name = 'spacetime', ncomponents = 1, spacetime = ic_1Dgauss_pulse2 },
--    { 'spacetime' },
--    { name = 'difference', ncomponents = 1, dep = {'density','spacetime'} },
--    {'velocity',sub_index={1}}
              }, -- variable table
   reduction = {'sum', 'sum', 'sum','max','max','max','max'},
   folder = 'tracking/',
   shape = {kind = 'all' },
   format = 'ascii',
   time = {min = 0, max = tmax, interval = 2},
  }
}

,
output = { active = true, -- VTK output activated?
  folder = 'output/',     -- Output location
  vtk = true,             -- VTK output activated?
  time = { -- defining when and how often to do outputs
          min = 0,      -- starting iteration
          max = tmax,     -- closing iteration (negative = forever)
          interval = interval} -- output after each iteration
  } 
}
