#!/usr/bin/env bash
set -eo pipefail
source /opt/openfoam13/etc/bashrc
mkdir -p results constant/triSurface
cp input/Al_Rihla_2022.stl constant/triSurface/Al_Rihla_2022.stl
for app in blockMesh snappyHexMesh checkMesh foamRun foamToVTK foamListTimes; do
  command -v "$app" >/dev/null || { echo "Missing application: $app" >&2; exit 2; }
done
foamVersion | tee results/environment.txt
blockMesh > results/blockMesh.log 2>&1
snappyHexMesh > results/snappyHexMesh.log 2>&1
checkMesh 2>&1 | tee results/checkMesh.log
grep -q 'Mesh OK.' results/checkMesh.log || exit 3

foamRun 2>&1 | tee results/steady_startup.log
rm -rf postProcessing
latest=$(foamListTimes -latestTime)
rm -f "$latest/uniform/time"
end_time=$(python3 -c "print(f'{float("$latest") + 0.05:.8f}')")
sed "s/END_TIME_VALUE/$end_time/" system/controlDict.transient > system/controlDict
cp system/fvSchemes.transient system/fvSchemes
cp system/fvSolution.transient system/fvSolution

foamRun 2>&1 | tee results/urans.log
cp -r postProcessing results/
foamToVTK -fields '(U p)' -time '300:' > results/foamToVTK.log 2>&1
mv VTK results/VTK
cp -r constant/polyMesh results/polyMesh

python3 - <<PY
import json, pathlib, re
force_files=list(pathlib.Path('results/postProcessing').rglob('forceCoeffs.dat'))
log=pathlib.Path('results/urans.log').read_text(errors='replace')
times=[float(x) for x in re.findall(r'^Time = ([0-9.eE+-]+)',log,re.M)]
summary={
 "case":"Al Rihla refined-mesh URANS non-rotating orthogonal-correction baseline revision",
 "speed_m_s":25.0,
 "diameter_m":0.22,
 "air_kinematic_viscosity_m2_s":1.5e-5,
 "reynolds_number":25.0*0.22/1.5e-5,
 "model":"transient incompressible URANS kOmegaSST",
 "mesh_cells":263891,
 "transient_requested_duration_s":0.05,
 "transient_start_time":float("$latest"),
 "transient_final_time":times[-1] if times else None,
 "force_file":str(force_files[0]) if force_files else None,
 "transient_velocity_convection":"bounded Gauss upwind",
 "max_courant_number":0.35,
 "limitations":["No prism layers","Short validation duration","First-order velocity convection used for robustness","Orthogonal pressure/gradient treatment used to avoid correction-path crashes","No mesh/time-step independence yet"]
}
pathlib.Path('results/summary.json').write_text(json.dumps(summary,indent=2)+'\n')
print(json.dumps(summary,indent=2))
PY
