----------------------- 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
                                                                          
--------------------------------------------------------------------!!!
-- Musubi configuration file. 
require "common"
tracking_fol = 'tracking/'
restart_fol = 'restart/'

-- Simulation name
simulation_name = 'poisson'
mesh = 'mesh/' -- Mesh information
printRuntimeInfo = false

scaling = 'diffusive'
control_routine = 'fast'
io_buffer_size = 10 -- default is 80 MB
logging = {level=10}

tmax = 1000000
-- Time step settigs
sim_control = {
  time_control = { 
    min = 0,
    max = {iter=tmax},
    interval = {iter=tmax/100000}
  } -- time control
 ,abort_criteria = {
    steady_state = true,
    convergence = {
      variable = {'potential_phy'}, 
      shape = {
        kind = 'all',
      },
      time_control = {min = 0, max = {iter=tmax}, interval = {iter=10}},
      reduction = 'average',
      norm='average', nvals = 100, absolute = true,
      condition = { threshold = 1.e-10, operator = '<=' }
    }
  }
} -- simulation control

-- restart 
NOrestart = {
      NOread = restart_fol..'Poisson_lastHeader.lua',
      write = restart_fol,
 }

logging = {level=10}

-- needed to dump variable in physical unit
physics = { dt = dt, rho0=1000.0, temp0=273.0 }

poisson = { 
  potential_diffusivity = 0.167, 
  permittivity=permit, 
  poisson_boltzmann = {
    moleDens0 = moleDens0,
    temp = temp, 
    valence={1.0,-1.0}, 
  }  
}

-- Initial condition 
initial_condition = { potential = ref_pot }

identify = {
  layout='d2q9',
  kind='poisson_boltzmann_nonlinear', 
  relaxation = 'bgk'
}

-- Parameter for analytical solution
k = math.sqrt(2*moleDens0*N_A*valence_sqr*charge^2/(permit*k_b*temp))*length
--print(k)
function analy_pot(x,y,z,t)  
  -- reciprocal of debye length
  ek = math.exp(k)
  eminusk = math.exp(-k)
  term_1 = (ek - 1)/(ek-eminusk)*math.exp(-k*x/length)  
  term_2 = (1-eminusk)/(ek-eminusk)*math.exp(k*x/length)
  return (term_1+term_2)*ref_pot 
end

-- Boundary conditions
boundary_condition = {  
{ label = 'left', 
   kind = 'potential_dirichlet',
   --qVal = 0.5,
   --curved = true,
   potential = ref_pot},
{ label = 'right', 
   kind = 'potential_dirichlet',
   --qVal = 0.5,
   --curved = true,
   potential = ref_pot},
}


variable = {
  {
    name = 'analy_pot',
    ncomponents = 1,
    vartype = 'st_fun',
    st_fun = analy_pot
  },
  {
    name = 'error',
    ncomponents = 1,
    vartype = 'operation',
    operation = {
      kind = 'difference',
      input_varname = {'potential_phy','analy_pot'}
    }
  }
}

-- Tracking
tracking = {
  {
    label = 'probe', 
    folder = tracking_fol,
    variable = {'potential_phy'}, 
    shape = {
            kind = 'canoND', object = {origin ={length*0.5,dx/2.0,zpos}}
    },
    time_control = {min = 0, max = {iter=tmax}, interval = {iter=10}},
    output={format = 'ascii'}
  },
  {
    label = 'line', 
    folder = tracking_fol,
    variable = {'potential_phy','analy_pot','error','electric_field_phy', 'charge_density_phy'}, 
    shape = {
      kind = 'canoND', 
      object = {
        origin ={-dx,dx/2.0,zpos},
        vec = {length+2*dx,0.0,0.0},
        segments = nLength+2
      }
    },
    time_control = {min = {iter=tmax}, max = {iter=tmax}, interval = {iter=tmax}},
    output={format = 'asciiSpatial'}
  },
  {
   label = 'l2error', 
   folder = tracking_fol,
   variable = {'error', 'analy_pot'}, 
   shape = {kind = 'canoND', 
     object = {origin = {-dx,dx/2.0,zpos},
     vec = {{ length+2*dx, 0.0,0.0}},
     segments = nLength+2
     }
   },
   reduction = {'l2norm','sum'},
   time_control = {min = tmax, max = tmax, interval = tmax},
   --time_control = {min = {iter=tmax}, max = {iter=tmax}, interval = {iter=tmax}},
   output={format = 'ascii'}      
  },
}

