require "seeder"

nRotations = 1
interpolation_method = 'quadratic'
preIterate=false
model='lbm'
rel='mrt'
cutoff=false
-- set these physical parameters
Re = 150
radius_rot = 0.2*length   -- [m]
circulation = 32 -- [m^2/s]
-- see equation on lower page 37 
-- Skriver Dahl 1996 (PhD)
tPeriod = 8*math.pi^2*radius_rot^2/circulation
tEnd  = nRotations*tPeriod

if preIterate then
  relaxation = rel..'_fixedVel'
else
  if cutoff then 
    relaxation = rel..'_omega_spatial'
  else
    relaxation = rel
  end
end

-- Musubi configuration file. 
rho0LB = 1
cs2LB  = 1./3.
originX = 00.0
originY = 00.0
originZ = 00.0
origin = {originX, originY, originZ}
--[]
-- sponge increase ration
-- value in the sponge (viscosity is set higher than inside the fluid)
w_sponge = 10. 
-- inner (fluid) increase ratio
-- has to match the inner fluid region's value -> 100% = 1.
w_fluid = 1.    
-- sponge should start at a radius of 85% from the origin
cutoff_rmin = 0.85
-- sponge should end  at a radius of 97% from the origin
cutoff_rmax = 0.97
-- start of sponge region (in spherical coordinates)
cutoff_min = cutoff_rmin*length*0.5
-- end of sponge region (in spherical coordinates)
cutoff_max = cutoff_rmax*length*0.5

-- define the spatial function for the sponge
-- this factor will be multiplied with the viscosity in Musubi
-- define this as a multiplicator to viscosity
function spatialFunction(x,y,z,t)
   -- this is a spherical function, so compute the radius from the origin
   r = math.sqrt( (x - originX)^2 + (y-originY)^2)
   --print("x: "..x.."y: "..y.."r:"..r)
   if(  r < cutoff_min) then
     -- if we're inside the fluid region, the viscosity should not be changed -> w_fluid
     res =w_fluid
   elseif( r  >= cutoff_max ) then
     -- if we're inside the fluid region, the viscosity should not be changed -> w_fluid
     res =  w_sponge
   else
     -- if we're between fluid and sponge, apply a linear progression from w_fluid to w_sponge
     slope = (w_fluid-w_sponge)/(cutoff_min-cutoff_max)
     res =  r*slope + w_fluid - cutoff_min*slope
   end
   --print ('result'..res)
   -- return the current value, which is a multiplicator to viscosity
   return res
end

-- Definition in physical units
csPhys = 300    -- m/s
rho0 = 1
Ma = circulation/(4*math.pi*radius_rot*csPhys)
-- Re = circulation*rho0/viscosity see K. Liow 2007 last page
viscosity = circulation*rho0/Re
print( 'rotation Mach number '.. Ma )
--b = radius_rot( cos( omega * t ) + math.i*sin( omega*t) ) 
--MaLB = circulationLB/(4*math.pi()*radius_rotLB*c_sLB)
-- Determine physical quantities
csLB = 1./math.sqrt(3.)
dt = csLB/csPhys*dx
viscLB = viscosity*dt/dx/dx
omega = 1./(3.*viscLB + 0.5)
if model == 'lbm' then
  p0 = rho0LB*cs2LB*dx*dx/dt/dt*rho0 --csPhys
else -- in the incompressible model, use reference pressure 0
  p0 = 0
end
print ('ref press'..p0)

halfwidth = length*0.05
amplitude = 0.0001
background = 0.0
function ic_1Dgauss_pulse(x, y, z, t)
  return p0+amplitude*math.exp(-0.5/(halfwidth^2)*( x - originX )^2)
end
function ic_2Dgauss_pulse(x, y, z, t)
  return p0+amplitude*math.exp(-0.5/(halfwidth^2)*(( x - originX )^2+( y - originY )^2))
end
-- Debug initial condition gaussian pulse 2d
nitial_condition = { 
  pressure  = ic_2Dgauss_pulse, 
  velocityX = 0., velocityY=0., velocityZ=0.} 
-- Initial condition 
ini_pressure = { predefined = "crvpPressure", 
                radius_rot = radius_rot,
                t          = 0,
                cs         = csPhys, 
                p0         = p0,
                cutoff_length = length*0.5,
                cutoff_rmin = cutoff_rmin,
                cutoff_rmax = cutoff_rmax,
                circulation= circulation,
                center = origin}
initial_condition = { 
  pressure  = ini_pressure,
  velocityX = { predefined = "crvpX", 
               cs         = csPhys, 
               radius_rot = radius_rot,
               circulation= circulation,
               cutoff_length = length*0.5,
               cutoff_rmin = cutoff_rmin,
               cutoff_rmax = cutoff_rmax,
               center = origin },
   velocityY = { predefined = "crvpY", 
               cs         = csPhys, 
               radius_rot = radius_rot,
               circulation= circulation,
               cutoff_length = length*0.5,
               cutoff_rmin = cutoff_rmin,
               cutoff_rmax = cutoff_rmax,
               center = origin },
  velocityZ = 0.0 }
-- Simulation name
simulation_name = 'crvp'
mesh = 'mesh/' -- Mesh information
ebug = { 
  debugMode = true, 
  debugFiles = true, 
  dumpDependencies = false
}

fluid = { 
  visc_spatial = { 
    fun = spatialFunction, 
    store = true},
  omega = omega, 
  rho0 = rho0 }

identify = {
  kind = model, 
  layout = stencil, 
  label = 'crvp',  
  relaxation=relaxation 
}

physics = { 
  dt = dt, 
  rho0 = rho0, 
}

bc_type2 = 'wall'
bc_type = 'wall'
bc_type = 'outlet_nrbc'
bc_type2 = bc_type

-- Boundary conditions
boundary_condition = {  
{ label = 'northwest',
   kind = 'wall', pressure = p0, sigma = 0.001, kappa = 1.0, L = length*0.5 }
,{ label = 'east',
   kind = bc_type,  pressure = p0, sigma = 0.001, kappa = 1.0, L = length*0.5 }
,{ label = 'west',
   kind = bc_type2, pressure = p0, sigma = 0.001, kappa = 1.0, L = length*0.5 }
,{ label = 'north',
   kind = bc_type2, pressure = p0, sigma = 0.001, kappa = 1.0, L = length*0.5 }
,{ label = 'south',
   kind = bc_type2, pressure = p0, sigma = 0.001, kappa = 1.0, L = length*0.5 }
}
                     
interval = tEnd/10
-- Time step settigs
sim_control = {
  time_control = {
    max = tEnd,
     interval = interval,
    }
}



-- Tracking              
tracking = {
 { label = 'stat_l'..level, 
 variable = {'vel_mag'}, 
 reduction = {'max'},
 shape = {kind = 'all' },
 time_control =  {
   interval = {iter=20}, 
   min = 0, 
   max =tEnd, 
 },
 format = 'ascii',
 folder = './tracking/'
 },
 { label = 'pressCutX_l'..level, 
 variable = {'pressure','velocity'}, 
 shape = {
   kind = 'canoND', 
   object = {
     origin = {-length*0.5,0.,0.}, 
     vec = {length, 0.,0.}, 
     segments = 2*nElemsMax } 
   },
   time_control = {interval = interval, min =0, max =tEnd},
   format = 'asciiSpatial',
   folder = './tracking/'
 },
 { label = 'hvs_l'..level, 
 variable = {'pressure_phy','velocity_phy', 'density','velocity'}, 
 --shape = {
 --  kind = 'canoND', 
 --  object = {
 --    origin ={-length*0.2,-length*0.15,0.}, 
 --    vec = {
 --      {length*0.7, 0.,0.},
 --      {0., length*0.5,0.}
 --    }, 
 --    segments = 
 --      {2*nElemsMax,nElemsMax} 
 --    } 
 --  },
 shape = {kind = 'all', object = {origin ={1.,0.,0.} } },
 time_control =  {interval = interval, min = 0, max = tEnd }, 
 format = 'harvester',
 folder = './tracking/'
 },
 { label = 'probePressure_l'..level, 
 variable = {'density','pressure_phy'}, 
 shape = {
   { kind = 'canoND', 
     object = {origin ={1.,0.,0.} } }, 
   { kind = 'canoND', 
     object = {origin ={length*0.5*0.5,0.,0.} } }
   },
 time_control =  {interval ={iter= 1}, min = 0, max= tEnd }, 
 format = 'ascii',
 folder = './tracking/'
 }
} 

if preIterate == true then
  restart = { write = 'restart/',
              time_control = { 
                min = tEnd, 
                max = tEnd, 
                interval = 1
              }
            }
end
if restart_read == true then
  restart = { read = 'restart/crvp_lastHeader.lua',
              write = 'restart_read/',
              time_control = { min = 1, max = tEnd, interval =  1},
}
end

-- Output settings
--utput = { active = true, -- VTK output activated?
--  folder = 'output/',     -- Output location
--  vtk = true,             -- VTK output activated?
  --dumpGhosts = true,
  --dumpHalos = true,
--  time_control = {interval = interval, min = 0, max = -1}
--}


if verbose then
  print ("Current level     " .. level)
  print ("Current omega     " .. omega)
  print ("----------------  ")
  print ("Physical units    ")
  print ("Reference length  m " .. lPhys )
  print ("Reference vel     m/s".. csPhys )
  print ("Spatial step size m " .. dx)
  print ("   Time step size s " .. dt)
end

