#!/bin/bash
cd ${0%/*} || exit 1    # Run from this directory

# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions

plotFile="postProcessing/bubblesPlot/0/LagrangianFieldValue.dat"

# Controls
maxTimeStepFraction="1e300 1e300 1e300 0.5 0.25"
nSubSteps="1 1 1 2 4"
nCorrectors="0 1 2 1 1"
pColour="web-blue"
Tcolours="dark-red red orange-red orange gold"
TanalyticColour="web-green"

# Construct substitutions and titles from the control lists
makeSubstitutions() { for arg in ${!1}; do echo $1=$arg; done }
substitutionss=$(paste -d ',' \
    <(makeSubstitutions maxTimeStepFraction) \
    <(makeSubstitutions nCorrectors))
makeTitle() { for arg in ${!1}; do echo "$arg $2"; done }
titles=$(paste -d ', ' \
    <(makeTitle nSubSteps "Sub Steps" | sed 's/1/No/g') \
    /dev/null \
    <(makeTitle nCorrectors "Corrector" | sed 's/^0/No/g;s/^[^1].*/&s/g'))

# Mesh and run for various parameters
runApplication blockMesh
for substitutions in $substitutionss
do
    runApplication -a foamDictionary \
        system/Lagrangian/bubbles/LagrangianSolution \
        -set "$substitutions"

    runApplication -s $substitutions foamRun

    cp --update=none $plotFile $plotFile.$substitutions
done

# Get the thermodynamic properties
R="6.02214e+26*1.38065e-23"
W=$(foamDictionary constant/Lagrangian/bubbles/physicalProperties \
    -entry mixture/specie/molWeight -value)
Cv=$(foamDictionary constant/Lagrangian/bubbles/physicalProperties \
    -entry mixture/thermodynamics/Cv -value)
gamma="(1+$R/$W/$Cv)"

# Get the initial thermodynamic state
p0=$(grep -m 1 -v "^#" "$plotFile" | awk '{print $4}')
T0=$(grep -m 1 -v "^#" "$plotFile" | awk '{print $5}')

# Analytic temperature for an adiabatic expansion
Tanalytic="($T0*(\$4/$p0)**(($gamma - 1)/$gamma))"

gnuplot << EOF

set terminal postscript eps color enhanced size 4,4
set output "expansion.eps"

substitutionss = split("${substitutionss//"
"/\\\n}", "\n")
titles = split("${titles//"
"/\\\n}", "\n")
Tcolours = split("$Tcolours")

set multiplot layout 2,1 margins 0.17, 0.83, 0.1, 0.95 spacing 0.1

set xlabel 'Time [s]'
set ytics nomirror
set ylabel 'Pressure [kPa]'
set y2tics
set y2label 'Temperature [K]'
set key bottom left Left reverse

plot "$plotFile" us 1:(\$4/1e3) w l lc '$pColour' t 'Pressure', \
    for [ i = 1:|substitutionss| ] \
    "$plotFile.".substitutionss[i] us 1:5 axis x1y2 w l \
    lc rgb Tcolours[i] t 'Temperature, '.titles[i], \
    "$plotFile" us 1:($Tanalytic) axis x1y2 w l \
    lc '$TanalyticColour' dt '-' t 'Analytic Temperature'

set ytics mirror
set ylabel 'Temperature Error [%]'
unset y2tics
unset y2label
set format y '%g'
set log y
set key bottom right Right noreverse

plot [][1e-5:] for [ i = 1:|substitutionss| ] \
    "$plotFile.".substitutionss[i] us 1:(100*abs(\$5/$Tanalytic - 1)) w l \
    lc rgb Tcolours[i] t titles[i]

unset multiplot

EOF

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