in.general_surface

# SETUP --------------------------------------------------------------------------------------

shell "bash -c 'rm -f dumps/*.dat'"

shell "mkdir -p dumps"

units           si
seed            11111
dimension       3

# outflow in +x, periodic y,z. will inject at xlo (x=0)
boundary        o p p

# domain size (m)
variable        xmin equal -1.5
variable        xmax equal 1.1
variable        ymin equal -1.1
variable 	    ymax equal 1.1
variable	    zmin equal -1.1
variable	    zmax equal 1.1

variable	    Lx equal ${xmax}-${xmin}
variable        Ly equal ${ymax}-${ymin}
variable        Lz equal ${zmax}-${zmin}

create_box      ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax}

# Load atmospheric data from NRLMSIS (run: python3 tools/load_atm_data.py <altitude_km>)
include         data/atm.sparta

print "Loaded NRLMSIS data: rho=${rho} nrho=${nrho} T=${T} vx=${vx}"

variable 	    kB equal 1.380649e-23  # J/K
variable 	    d equal 3.7e-10  # m
variable 	    R equal 287.05 # (J / kg*K)
variable 	    lambda equal ${kB}/(sqrt(2.0)*PI*${d}*${d}*${rho}*${R}) # m
variable	    vbar equal sqrt(8*${kB}*${T}/(PI*${rho}/${nrho})) # m/s

print "MEAN FREE PATH = ${lambda} m"
print "MEAN MOLEC VEL = ${vbar} m/s"

# grid resolution

variable	    dx equal ${lambda}/3 # m
variable	    mct equal ${lambda}/${vbar} # s
variable	    mtt equal ${dx}/${vbar} # s
variable	    dt equal ((${mct}<${mtt})*${mct}+(${mct}>=${mtt})*${mtt})/3.0 # min(${mct},${mtt})/3.0

print "CELL SIZE MUST BE LESS THAN ${dx} m"
print "MEAN COLL TIME = ${mct} s"
print "MEAN TRANSIT TIME = ${mtt} s"
print "TIMESTEP MUST BE < ${dt} s" # timestep set at end of script

create_grid	    150 100 50
balance_grid    rcb part

variable	    diag_freq   equal 100 # dump diagnostics every _ timesteps
variable 	    tstep       equal 1.0e-7 # choose small timestep << cell flight time; start 1e-7 s

# SURFACE GEOMETRY ----------------------------------------------------------------------------

read_surf	    surf/AMPT_sat_inlet.surf  group ampt        # create surface group "ampt"

# SPECIES AND MIXTURE -------------------------------------------------------------------------

# set mixture properties from NRLMSIS data
mixture     	atm nrho ${nrho} vstream ${vx} 0.0 0.0 temp ${T}

# particle weighting: set target of 2e5 sim particles, calculate weighting factor
variable        Ns_target equal 2000000.0
variable        Vol       equal ${Lx}*${Ly}*${Lz}
variable        fnum      equal ${nrho}*${Vol}/${Ns_target}
global          fnum ${fnum} # global bc will be queried by each new particle

variable	    parts_per equal ${nrho}/${Ns_target}
print "particles per particle = ${parts_per}"

# global useful for stat query later
global          nrho ${nrho}

# create an initial fill (n 0 -> auto compute # of particles from fnum and nrho)
create_particles	atm n 0

# continuous inflow; inject gas from xlo every step (fix -> run each time step, ID: "in")
fix		        in emit/face atm xlo

# SURFACE COLLISIONS --------------------------------------------------------------------------

                # compute ID, type, surf group, mixture, property (energy flux on surface)
compute         compute_qwall surf ampt atm ke # ke = kinetic + internal energy. W/m^2 *multi-column array for all surface groups

                # running‑average every step so flux is never exactly zero (exponential decay)
fix             flux ave/surf ampt 1 1 1 c_compute_qwall[1] ave running  # running mean, updates each step
fix             fix_Tsurf surf/temp ampt 1 f_flux 300.0 0.9 Tsurf  # Stefan-Boltzmann

# Define collision models
surf_collide 	wall diffuse s_Tsurf 0.9           # define collision model "wall", random dir, use local facet temp, acc
surf_modify 	ampt collide wall                  # attach model to facets

# DIAGNOSTICS ---------------------------------------------------------------------------------

                # per-grid properties: nrho, massrho, pxrho (for flow calcs)
compute         prop_grid grid all atm nrho massrho pxrho
                # per-grid temperature
compute         prop_grid_temp grid all atm temp
                # per-grid pressure (thermal/grid supports press property)
compute         prop_grid_press thermal/grid all atm press
                # per-grid velocity components (for streamlines)
compute         prop_grid_flow grid all atm u v w
                # single grid dump with all properties
dump            dump_grid grid all ${diag_freq} dumps/grid.*.dat id xlo ylo zlo xc yc zc &
                    c_prop_grid[*] c_prop_grid_temp[*] c_prop_grid_press[*] c_prop_grid_flow[*]

                # ID, data type, mixture, every _ steps, filename, columns (particle data)
dump            dump_part particle atm ${diag_freq} dumps/part.*.dat id type x y z vx vy vz

		        # ID, data type, region, every _ steps, filename, columns (facet id, triangle vertices, surface temperature)
dump    	    dump_surf surf ampt ${diag_freq} dumps/surf.*.dat id v1x v1y v1z v2x v2y v2z v3x v3y v3z f_flux[*] s_Tsurf

# Measure simulated flow properties for domain ---------------------

compute         vol_grid  property/grid all vol

# Use grid-style variables to multiply density by volume for each cell
variable        nrho_vol grid c_prop_grid[1]*c_vol_grid
variable        rho_vol  grid c_prop_grid[2]*c_vol_grid
variable        px_vol   grid c_prop_grid[3]*c_vol_grid

compute         sum_vol      reduce sum c_vol_grid
compute         sum_nrho_vol reduce sum v_nrho_vol
compute         sum_rho_vol  reduce sum v_rho_vol
compute         sum_px_vol   reduce sum v_px_vol

variable        nrho_avg equal c_sum_nrho_vol/c_sum_vol
variable        rho_avg  equal c_sum_rho_vol/c_sum_vol
variable        ux_avg   equal c_sum_px_vol/c_sum_rho_vol

fix             boxprops ave/time 1 1 100 v_ux_avg v_rho_avg v_nrho_avg ave running file dumps/box_props.dat

# DRAG CALC (short for caclulator) ------------------------
# sum of surface forces (direct):
                # compute per-surface-element force components (fx,fy,fz) for all ampt surfaces
compute         surfF surf ampt atm fx fy fz
                # time-average per-surf forces (running mean)
fix             surfavg ave/surf ampt 1 1 1 c_surfF[*] ave running
                # reduce (sum) per-surf averaged fx -> global total drag
compute         drag reduce sum f_surfavg[1]

# compute forces for x-normal surfaces (ram drag)
# NOTE: For complex geometry, we skip directional grouping.
# Fill component drags with zeros to keep downstream analysis compatible.
variable        drag_xnorm equal 0.0

# compute forces for y/z-normal surfaces (skin friction)
variable        drag_walls equal 0.0

# write drag components (timestep, total_drag, ram_drag, skin_friction) to file
fix             dragout ave/time 1 1 100 c_drag v_drag_xnorm v_drag_walls file dumps/direct_drag.dat mode scalar

# ----------------------------------------------------------

# cell-averaged (streaming+thermal) temperature
compute         Tbox temp # define a compute Tbox that calculates domain everaged temp
stats           ${diag_freq} # print diagnostics every _ timesteps
stats_style     step cpu np nattempt ncoll c_Tbox c_drag v_drag_xnorm v_drag_walls # print timestep, runtime, particles, collision stats, avg temp, drag components

timestep        ${tstep}
collide		    vss atm vss/air.vss # variable soft sphere model
run             5000

shell "rm -f data/rho.dat data/nrho.dat data/T.dat data/vx.dat"