----------------------- 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"
folderRead = 'restart_read/'
folderWrite = 'restart/'
mode = 'write'
if os.getenv('mode') == 'read' then
  mode = 'read'
end 
-- General test case geometry and setup settings
kernelDim = 3
dimension = 2
Re     = 25
periodic = true
testing = false
-- 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'
omm_reduced = false
ontrol_routine= 'debug'
nit_allTimeLayers = true

label = 'TGV'
maxlevel = level+refinement
if iterateIC then
  relaxation = relaxation..'_fixedVel'
  tEnd = 14.
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 = dt

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)

if testing == false then
  function pressure(x,y,z,t)
    return (randFluc*math.random()+0.25)*(u0)^2*(math.cos(2.*(x-x0[1]))+math.cos(2.*(y-x0[2])))*math.exp(-2.*k^2*viscosity*t)
    end
  function velocityX(x,y,z,t)
    return  (u0+randFluc*math.random())*math.sin((x-x0[1]))*math.cos((y-x0[2]))*math.exp(-k^2*viscosity*t)
    end
  function velocityY(x,y,z,t)
    return -(u0+randFluc*math.random())*math.cos((x-x0[1]))*math.sin((y-x0[2]))*math.exp(-k^2*viscosity*t)
    end
  function Sxx(x,y,z,t)
    return
    2.*(u0+randFluc*math.random())*math.cos((x-x0[1]))*math.cos((y-x0[2])) --*viscosity 
    --*math.exp(-k^2*viscosity*t)
    end
  function Sxy(x,y,z,t)
    return  0.*viscosity
    end
  function Syy(x,y,z,t)
    return
    -2.*(u0+randFluc*math.random())*math.cos((x-x0[1]))*math.cos((y-x0[2])) --*viscosity
    --*math.exp(-k^2*viscosity*t)
    end
else
  originX = 9.*math.pi/7.
  originY = 6.*math.pi/7.
  amplitude = 10.
  halfwidth = length/20.
  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
end

if dimension == 2 then
  function ic_velocityX(x,y,z,t) 
    return  velocityX(x,y,z,0.)
    end
  function ic_velocityY(x,y,z,t) 
    return  velocityY(x,y,z,0.)
    end
  function ic_velocityZ(x,y,z,t) -- initial conditions given in physical units, so convert to LB
    return 0.
    end
  function ic_pressFluc(x,y,z,t) -- 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
  function ic_Sxx(x,y,z,t)
    return Sxx(x,y,z,0.)
    end
  function ic_Syy(x,y,z,t)
    return Syy(x,y,z,0.)
    end
  function ic_Sxy(x,y,z,t)
    return Sxy(x,y,z,0.)
    end
   function ic_Szz(x,y,z,t)
    return 0.
    end
   function ic_Syz(x,y,z,t)
    return 0.
    end
   function ic_Sxz(x,y,z,t)
    return 0.
    end
 
else  -- 3D
  function ic_velocityX(x,y,z,t) -- initial conditions given in physical units, so convert to LB
    -- for gamma = 0 amplitude is +1
    res = amplitudeX*math.sin(x-x0[1])*math.cos(y-x0[2])*math.cos(z-x0[3])--*dt/dx
    return res
  end
  function ic_velocityY(x,y,z,t) -- initial conditions given in physical units, so convert to LB
    -- for gamma = 0 amplitude is -1
    res =  amplitudeY*math.cos(x-x0[1])*math.sin(y-x0[2])*math.cos(z-x0[3])--*dt/dx
    return res
  end
  function ic_velocityZ(x,y,z,t) -- initial conditions given in physical units, so convert to LB
    -- for gamma = 0 this is 0
    return amplitudeZ*math.cos(x-x0[1])*math.cos(y-x0[2])*math.sin(z-x0[3]) --*dt/dx
  end
  function ic_pressFluc(x,y,z,t) -- initial conditions given in physical units, so convert to LB
-- From ATeles return rho0 * math.pow(u0,2)/16.*(math.cos(2.*x)+math.cos(2.*y))*(math.cos(2.*z)+2.)
return ((1-math.cos(2.*gamma ))/24.*(
math.cos(2.*(x-x0[1]))*math.cos(2.*(y-x0[2])) + 2.*math.cos(2.*(z-x0[3])))
+ (2.+math.cos(2.*gamma)+math.sqrt(3.)*math.sin(2.*gamma))/48.*(
math.cos(2.*(x-x0[1]))*math.cos(2.*(z-x0[3])) + 2.*math.cos(2.*(y-x0[2])))
  + (2.+math.cos(2.*gamma)-math.sqrt(3.)*math.sin(2.*gamma))/48.*(
  math.cos(2.*(y-x0[2]))*math.cos(2.*(z-x0[3])) + 2.*math.cos(2.*(x-x0[1]))))
  end
  function ic_Sxx(x,y,z,t)
    return 2.*(u0+randFluc*math.random())*math.cos((x-x0[1]))*math.cos((y-x0[2]))*math.cos((z-x0[3]))
    end
  function ic_Syy(x,y,z,t)
    return -2.*(u0+randFluc*math.random())*math.cos((x-x0[1]))*math.cos((y-x0[2])) *math.cos((z-x0[3])) 
    end
  function ic_Szz(x,y,z,t)
    return  0.
    end
  function ic_Sxy(x,y,z,t)
    return  0.
    end
  function ic_Syz(x,y,z,t)
    return  (u0+randFluc*math.random())*math.cos((x-x0[1]))*math.sin((y-x0[2])) *math.sin((z-x0[3]))
    end
  function ic_Sxz(x,y,z,t)
    return  -(u0+randFluc*math.random())*math.sin((x-x0[1]))*math.cos((y-x0[2])) *math.sin((z-x0[3]))
    end
end
function ic_density(x,y,z,t)
  --if scaling == 'diffusive' then
    return ic_pressFluc(x,y,z) --/rho0/dx^2*dt^2/cs2LB + rho0
  --else
  --  return ic_pressFluc(x,y,z) + rho0
  --end 
  end

function refPressure(x,y,z,t)
  return pressure(xProbe[1],xProbe[2],xProbe[3],t) -- +rho0LB*cs2LB*dx*dx/dt/dt
  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.))
  res = ic_density(xProbe[1], xProbe[2], xProbe[3])
  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.
-- here: timeControl is given in physical simulation time
sim_control = { 
  time_control = { max = tEnd, -- maximum simulation time to reach 
                   interval = interval }   -- Interval for checking density
}

logging = {level=10}
-- physics table to convert unit from physics to lattice and vice versa.
-- if this table is present, all input paramters to musubi like
-- fluid/species property, initial_condition, boundary condition (amplitude
-- and transient times )must defined in physical units.
physics = { dt = dt, rho0 = rho0 }

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

-- scheme model for single fluid simulation
scheme = {
  identify = {
      label = 'global',  -- scheme name 
      kind = 'lbm_incomp',      -- simulation type of this scheme
      relaxation = relaxation, -- relaxation type (bgk, mrt, ...)
      layout = stencil},
   fluid = { omega = omega, 
     rho0 = rho0 },
-- Initial condition for each field
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, 
     Sxx = ic_Sxx,
     Syy = ic_Syy,
     Szz = ic_Szz,
     Sxy = ic_Sxy,
     Syz = ic_Syz,
     Sxz = ic_Sxz
     }

, variable = { 
    { name = 'refPressure', 
      ncomponents = 1, 
      vartype = 'st_fun',
      st_fun = refPressure },
    { name ='refVelMag  ',  
      ncomponents = 1, 
      vartype = 'st_fun',
      st_fun = refVelocity },
    { name = 'refVelocityX', 
      ncomponents = 1, 
      vartype = 'st_fun',
      st_fun = refVelocityX },
  }
, 
tracking = {{
   label = label..dimension..'D_l'..level..'_'..relaxation..'_press',
   variable = { 
    'pressure_phy', 
    'refPressure', 
    'density', 
    'velocity', 
    'wss', 
    'pressure', 
              }, 
   folder = 'tracking/',
   shape = {kind = 'canoND', object = {origin = xProbe } },
   format = 'ascii',
   time_control = {min = 0, max = tEnd, interval = minInt},
  }
  ,  
{
   label = label..dimension..'D_l'..level..'_'..relaxation..'_vel',
   variable = { 
    'velocity_phy', 
              }, 
   folder = 'tracking/',
   shape = {kind = 'canoND', object = {origin = xProbe } },
   format = 'ascii',
   time_control = {min = 0, max = tEnd, interval = minInt},
  }
  ,{
   label = label..dimension..'D_l'..level..'_'..relaxation..'_kE',
   variable = { 'kinetic_energy_phy', 'kinetic_energy', 'density'
                              }, -- variable table
   reduction = {'average', 'average', 'sum'},
   folder = 'tracking/',
   shape = {kind = 'all',
     object = { 
       origin = {0.,0.,length/2.},
       vec = { { length, 0., 0. }, 
               { 0., length, 0. }
       }, segments = { 2^level, 2^level }
     }
   },
   format = 'ascii',
   time_control = {min = 0, max = tEnd, interval = minInt},
  }  
  ,{
   label = label..dimension..'D_l'..level..'_'..relaxation..'_Z',
   variable = { 
     'pressure', 'velocity',
     'pressure_phy', 
     'velocity_phy'
   }, -- variable table
   folder = 'tracking/',
   shape = {kind = 'canoND',
     object = { 
       origin = x0,
       vec = { { length, 0., 0. }, 
               { 0., length, 0. }
       }, segments = {2^maxlevel, 2^maxlevel}
     }
   },
   format = 'harvester',
   time_control = {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_control = {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_control = {min = 0, max = tEnd, interval = interval},
--  }
}

 
}
if dimension == 2 then
  table.insert( scheme.tracking[1].variable, 1, 
    { name = 'difference', ncomponents = 1, dep = {'refPressure','pressure_phy'} })
  table.insert( scheme.tracking[2].variable, 1, 
    { name = 'difference', ncomponents = 1, dep = {'refVelMag','velMag_phy'} }
     )
end 
if iterateIC then
  restart = {
    write = './restart/', 
    time_control = {min = tEnd, max = tEnd }
  }
else 
  if startFromZero == false then
    restart = {
      read = './restart/'..label..'_lastHeader.lua'
    }
  end 
end
-- Read or write?
if mode == 'read' then
  print ('   READ MODE ')
  print ('   File: '..folderWrite..label..'_header_0.000E+00.lua')
restart = { read = folderWrite..label..'_header_0.000E+00.lua',
            write = folderRead,
            time_control = { min = 0, max = tEnd, interval = {iter=1}}}
elseif mode == 'write' then
  print ('   WRITE MODE ')
restart = { 
            write = folderWrite,
            time_control = { min = 0, max = tEnd, interval = {iter=1}}}
end

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

