require "common"
-- Scaling and initialization settings
verbose = false
if studyError then
  interpolation_method = { method = intp
   , test = true, --testEach = true --, testFluids = true
   } -- , noIntpFromFiner = true }
  a= {1., 0., 0., 0., 0., 0. }
  a= {1.3248373, 0.65673836745, 0.83865622008, 0.217345933, 1.393957623594, 0.9784636222, 1.018284747635}
  init_allElems = true
  startFromZero = true
  iterateIC = false
  tEnd = 0.0
else
  startFromZero = true
  interpolation_method = { method = intp }
  iterateIC = false
  tEnd = 10.  -- run up to which simulation time?
end

label = 'TGV'
if iterateIC then
  relaxation = relaxation..'_fixedVel'
  tEnd = 1.
end
-- choose stencil accroding to dimension
if kernelDim == 3 then
  stencil = 'd3q19'
elseif kernelDim == 2 then
  stencil = 'd2q9'
else
  stencil = 'd1q3'
end

viscosity = 1 / Re
u0     = 1.0 -- initial velocity magnitude
rho0LB = 1.
rho0   = 1.

if scaling == 'acoustic' then
  uLB = 0.05
  dt  = uLB / u0 * dx
  viscLB = viscosity*dt/dx/dx
  omega = 1./(3.*viscLB + 0.5)
else
  omega = 1.9
  viscLB  = 1./6.*(2./omega - 1.)
  dt = viscLB/viscosity*dx*dx
  uLB=u0*dt/dx
end


-- origin is the center of first element
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] + math.pi/3, x0[2] + math.pi/3, x0[3] }
--xProbe = { x0[1], x0[2]+2*dx, x0[3] }

t = 0.
cs2LB = 1./3.
p0 = rho0LB/rho0*dx^2/dt^2*cs2LB
p0 = 0

interval = tEnd/10.
minInt = dt

gamma = 0.0  -- geometry factor
kx = twoPi / length -- wave number in x direction
ky = twoPi / length -- wave number in y direction
k  = math.sqrt(kx^2+ky^2) -- rate of decay coefficient

-- when gamma == 1, amplitude = 1
amplitudeX  = 2.0/sqrt3 * math.sin(gamma + twoPi/3.0)
amplitudeY  = 2.0/sqrt3 * math.sin(gamma - twoPi/3.0)
amplitudeZ  = 2.0/sqrt3 * math.sin(gamma)

-- Analytical solution with respect to x0
if testing == false then
  function pressure(x,y,z,t)
    p = (math.cos(2.0*kx*(x-x0[1]))+math.cos(2.0*ky*(y-x0[2])))*math.exp(-2.*k^2*viscosity*t)
    return (randFluc*math.random()+0.25) * (u0)^2 * p + p0
  end

  -- velocity X
  function velocityX(x,y,z,t)
    -- when gamma == 1, amplitude = 1
    velX = math.sin(kx*(x-x0[1])) * math.cos(ky*(y-x0[2])) * math.exp(-k^2*viscosity*t)
    return (u0+randFluc*math.random()) * velX
  end

  -- velocity Y
  function velocityY(x,y,z,t)
    -- when gamma == 1, amplitude = 1
    velY = math.cos(kx*(x-x0[1])) * math.sin(ky*(y-x0[2])) * math.exp(-k^2*viscosity*t)
    return -(u0+randFluc*math.random()) * velY
  end

  -- velocity Z
  function velocityZ(x,y,z,t)
    return 0
  end

  -- Sxx
  function Sxx(x,y,z,t)
    sxx = (u0+randFluc*math.random()) * math.cos((x-x0[1]))*math.cos((y-x0[2]))
    return sxx * math.exp(-k^2*viscosity*t)
  end
  -- Syy
  function Syy(x,y,z,t)
    syy = (u0+randFluc*math.random()) * math.cos((x-x0[1]))*math.cos((y-x0[2]))
    return -syy * math.exp(-k^2*viscosity*t)
  end
  -- Sxy
  function Sxy(x,y,z,t)
    return  0.0
  end
else
  --originX = 9.*math.pi/7.
  --originY = 6.*math.pi/7.
  --amplitude = 10.
  --halfwidth = length/20.
  function polynomialQuadratic(x,y,z)
    return a[1] + a[2]*x + a[3]*y + a[4]*x^2 + a[5]*y^2 + a[6]*x*y
    end
  function pressure(x,y,z,t)
    return polynomialQuadratic(x,y,z)
    end
  function velocityX(x,y,z,t)
    return polynomialQuadratic(x,y,z)
    end
  function velocityY(x,y,z,t)
    return polynomialQuadratic(x,y,z)
    end
  function Sxx(x,y,z,t)
    return polynomialQuadratic(x,y,z)
    end
  function Syy(x,y,z,t)
    return polynomialQuadratic(x,y,z)
    end
  function Sxy(x,y,z,t)
    return polynomialQuadratic(x,y,z)
    end
end

if dimension == 2 then -- 2D

  -- initial velocity X
  function ic_velocityX(x,y,z)
    return velocityX(x,y,z,0.)
  end

  -- initial velocity Y
  function ic_velocityY(x,y,z)
    return velocityY(x,y,z,0.)
  end

  -- initial velocity Z
  function ic_velocityZ(x,y,z)
    return 0.
  end

  -- initial pressure
  function ic_pressFluc(x,y,z)
--    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

  -- initial Strain rate xx
  function ic_Sxx(x,y,z)
    return Sxx(x,y,z,0.)
  end
  -- initial Strain rate yy
  function ic_Syy(x,y,z)
    return Syy(x,y,z,0.)
  end
  -- initial Strain rate xy
  function ic_Sxy(x,y,z)
    return Sxy(x,y,z,0.)
  end
  -- initial Strain rate zz
  function ic_Szz(x,y,z)
    return 0.
  end
  -- initial Strain rate yz
  function ic_Syz(x,y,z)
    return 0.
  end
  -- initial Strain rate xz
  function ic_Sxz(x,y,z)
    return 0.
  end
else  -- 3D
  function ic_velocityX(x,y,z) -- 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) -- 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) -- 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) -- initial conditions given in physical units, so convert to LB
-- From ATeles return rho0 * (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 -- initial condition functions

function ic_pressure(x,y,z)
  if model == 'lbm_incomp' then
    return ic_pressFluc(x,y,z) --/rho0/dx^2*dt^2/cs2LB + rho0
  else
    return ic_pressFluc(x,y,z)
  end
end

function refPressure(x,y,z,t)
  return pressure(xProbe[1],xProbe[2],xProbe[3],t) -- +rho0LB*cs2LB*dx*dx/dt/dt
end

-- reference velocity magnitude at probe point
function refVelocity(x,y,z,t)
  velX = velocityX(x, y, z, t)^2
  velY = velocityY(x, y, z, t)^2
  velZ = velocityZ(x, y, z, t)^2
  return math.sqrt( velX + velY + velZ )
end
-- reference velocity X at probe point
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_pressure(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/' --mesh_l'..level..'_'..refinement..'/'
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.
sim_control = {
  time_control = {
    max = tEnd,           -- Maximal iteration to reach
    interval = interval }    -- Interval for checking density
}

physics = { dt = dt, rho0 = rho0 }
-- 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 = {
  identify = {
    label = 'global',  -- scheme name 
    kind = model,      -- simulation type of this scheme
    relaxation = relaxation, -- relaxation type (bgk, mrt, ...)
    layout = stencil
  }, -- identify table

-- interpolation method to use for the scheme
interpolation_method =  interpolation_method,

  fluid = { omega = omega, rho0 = rho0LB },
-- Initial condition for each field
initial_condition = {
-- initial with lua function
  pressure  = ic_pressure,--/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
}, -- intitial condition end

variable = {
  -- analytical solution
  { name = 'refPressure',  ncomponents = 1, vartype = 'st_fun', st_fun = pressure },
  { name = 'refVelMag  ',  ncomponents = 1, vartype = 'st_fun', st_fun = refVelocity },
  { name = 'refVelocityX', ncomponents = 1, vartype = 'st_fun', st_fun = velocityX },
  { name = 'refVelocityY', ncomponents = 1, vartype = 'st_fun', st_fun = velocityY },
  { name = 'refSxx',       ncomponents = 1, vartype = 'st_fun', st_fun = Sxx },
  { name = 'refSyy',       ncomponents = 1, vartype = 'st_fun', st_fun = Syy },
  -- error
  { name = 'error_vel', ncomponents = 1, vartype = 'operation',
    operation = {
      kind = 'difference',
      input_varname = {'refVelMag','vel_mag_phy'},
    },
  },
  { name = 'error_p', ncomponents = 1, vartype = 'operation',
    operation = {
      kind = 'difference', 
      input_varname = {'refPressure','pressure_phy'},
    },
  },
},

tracking = {
  { -- probe point
    label = dimension..'D_l'..level..'_'..relaxation..'_probe',
    variable = {
      -- for 2D, diff is inserted here
      -- { name = 'difference', ncomponents = 1, dep = {'refPressure','pressure_phy'} }
      'pressure_phy',
      'refPressure',
      'velocity_phy',
      'refVelocityX',
      'refVelocityY',
      'strain_rate_phy',
      'refSxx',
      'refSyy',
    },
    folder = 'tracking/',
    shape = {kind = 'canoND', object = {origin = xProbe } },
    output={format = 'ascii'},
    time_control = {min = 0, max = tEnd, interval = minInt},
  },
  { -- average kinetic energy over whole domain
   label = dimension..'D_l'..level..'_'..relaxation..'_kE',
   variable = {
     'kinetic_energy_phy',
     'kinetic_energy',
     'density',
   },
   reduction = { 'average', 'average', 'sum'},
   folder = 'tracking/',
   shape = { kind = 'all', },
   output={format = 'ascii'},
   time_control = {min = 0, max = tEnd, interval = minInt},
  },
  { -- whole domain as harvester format
   label = dimension..'D_l'..level..'_'..relaxation..'_Z',
   variable = {
     'pressure_phy',
     'velocity_phy',
     --'shear_stress_phy',
     'strain_rate_phy',
     'refPressure',
     'refVelocityX',
     'refVelocityY',
     'refSxx',
     'refSyy',
   }, -- variable table
   folder = 'tracking/',
   shape = { kind = 'global', },
   output={format = 'vtk'},
   time_control = { min = 0, max = tEnd, interval = interval},
  },
  { -- l2norm error of velocity
   label = dimension..'D_l'..level..'_'..relaxation..'_errVel',
   variable = {
     'error_vel',
   }, -- variable table
   reduction = {'l2norm'},
   folder = 'tracking/',
   shape  = { kind = 'all', },
   output={format = 'ascii'},
   time_control = {min = 0, max = tEnd, interval = interval},
  },
  { -- l2norm error of pressure
   label = dimension..'D_l'..level..'_'..relaxation..'_errPress',
   variable = {
     'error_p',
   },
   reduction = {'l2norm'},
   folder = 'tracking/',
   shape = {kind = 'all',
     object = { 
       origin = x0,
       vec = { { length, 0., 0. }, 
               { 0., length, 0. }
       }, segments = {2*2^maxLevel, 2*2^maxLevel}
     }
   },
   output={format = 'ascii'},
   time_control = {min = 0, max = tEnd, interval = interval},
  },
} -- track table

} -- scheme table

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