Before we begin

About this presentation

• This presentation is about boring topics
  • Programming
  • Numerical evaluation
• It tries to compensate by being flashy
  • Needs a HTML5-capable browser
    • Press F for fullscreen. ESC for overview
• It demonstrate the methods using only a "toy project"
  • But that way all the files can later be made available

Where I work

• ICE Stroemungsforschung
  • Located in Leoben, Austria
  • CFD Consulting
    • Development with OpenFOAM
    • Also using ClosedSource

The rocket

DLR Water rocket

The question:

How much water is the optimum?

• Too much water: air "lost" pressure before all water is driven out
• Too little water: air is still pressurized when all the water is "gone"

We'd like to find the ideal spot with OpenFOAM

• Ideal means: rocket flies highest

What we use

We need a snake

Not that one

We need a cool snake

BTW: that was Snake Plissken from John Carpenters Escape from New York

Python - the friendly snake

Python is ..

• An object-oriented scripting language with a simple syntax

• Used as scripting language in a number of software packages
  • Paraview
  • Salome (Preprocessor)
• Comes with a number of libraries included
  • Text processing, Internet stuff, etc
• A number of libraries for scientific/numeric purposes exist
  • numpy and matplotlib give functionality equivalent to MATLAB
  • scipy adds numerical methods like ODE-integration and optimization
  • pandas adds a layer for data analysis

PyFoam

• .. is a library to help working with OpenFOAM
  • Manipulate cases
  • Control the execution of a run
    • "on the fly" analysis of the solver output (includes plotting)
  • Work with the resulting data files

• There is a number of utilities bundled with it
  • This is the part that is used most
  • … but it has already been discussed in other presentations
• http://openfoamwiki.net/index.php/Contrib_PyFoam
  • Download and further information
  • Announcements on the Message-board and Twitter @swakPyFoam

Preparation

How the case is set up

• Mesh generation
  • blockMesh for the basic mesh (half bottle)
  • makeAxialMesh to use the axial symmetry of the bottle
    • We assume that the bottle is always oriented vertically
• Dimensions of the bottle
  • Radius: 5 centimeter, height: 20 centimeter
  • Bottleneck width: 1 centimeter
• As a solver compressibleInterFoam is used
  • Handles two-phase flow and compressed air
• All cases are set up the same way except for two parameters:
  • Initialized with water up to a specified height
  • Initial pressure

The ParsedParameterFile

This class is the "workhorse" of PyFoam

If setFieldsDict looks like this

pInit 5e5;
fillHeight 0.15;

defaultFieldValues
(
    volScalarFieldValue alphawater 1
    volScalarFieldValue p_rgh $pInit
    volScalarFieldValue p $pInit
);

Then the values can be easily changed like this

from PyFoam.RunDictionary.ParsedParameterFile import ParsedParameterFile
setF=ParsedParameterFile(path.join(self.casedir(),"system","setFieldsDict"))
setF["pInit"]=7e5
setF["fillHeight"]=0.2*0.5
setF.writeFile()

Calling utility from a script

• Every PyFoam-utility is implemented as a class
  • Results of the execution are stored in the object and can be accessed
• For instance if we call a utility on the command-line

pyFoamRunner.py --clear --progress simpleFoam

we can run it from a Python-script like this

from PyFoam.Applications.Runner import Runner
data=Runner(args=["--clear","--progress",
           "simpleFoam","-case",caseName]).getData()

and the variable data then holds the analyzed data from that run

A script to send to the cluster

def setup(self,parameters):
    self.execute("rm -rf 0")
    self.execute("cp -r 0.orig 0")
    self.templateFile("constant/polyMesh/blockMeshDict")
    self.foamRun("blockMesh")
    self.foamRun("makeAxialMesh",foamArgs=["-overwrite"])
    self.foamRun("collapseEdges",foamArgs=["-overwrite"])
    setF=ParsedParameterFile(path.join(self.casedir(),"system","setFieldsDict"))
    setF["pInit"]=parameters["p"]
    setF["fillHeight"]=parameters["frac"]*0.2
    setF.writeFile()
    self.foamRun("setFields")

A typical run

The initial conditions

Lift-off

Note: scaling of pressure is different

Full thrust

Break-through

The end

First analysis

Output analysis of PyFoam

• PyFoam analysis the screen output of the solver
  • Automatically extracts standard information like residuals (and allows plotting them)
  • Allows specification of additional patterns to look for
• Data is stored
  • Full timelines for quick plotting
  • A database of the last values

Looking at the acceleration

We've added an expression in swak4Foam to output the current acceleration of the bottle

$a_{thrust}=\frac{\dot m_{out}*v_{out}}{m(t)}$

like this

Expression acceleration :  min=20.443693

In a file customRegexp we add

bottleAcceleration {
               theTitle "Acceleration of the bottle";
               expr "Expression acceleration :  min=(%f%)";
               titles ( a );
}

and PyFoam extracts these values and writes them to disc

Acceleration of bottle

Current position of the bottle

swak4Foam also integrates

$\dot v=a_{thrust}-g$

$\dot x=v$

to

• get the position of the bottle
• the maximum height

PyFoam extracts that data for us as well ….

Position of bottle

(blue line is position without gravity)

Handling other data

PyFoam is not picky:

• it also works with data generated by OpenFOAM:
  • sets generated by sample
  • timelines generated by
    • probes
    • other functionObjects (did I mention swak4Foam?)
• Processes that data
  • through utilities
    • generate plots
    • writes them into files readable by Excel and friends
  • passes the data on to Python-libraries like numpy and pandas
  • makes it available to ParaView …

PyFoam can read the analyzed data into Paraview

Doing the variation

• Using two loops we vary the conditions
  • 10 different initial pressures
  • 19 filling fractions
• Next step would be going through 190 cases by hand and get
  • the maximum height
  • other data that might be of interest
• This is a tedious and error-prone process
  • Do we want to do this?

Answer:

No

Extracting the data

First we build a database from the cases

pyFoamAddCaseDataToDatabase.py runDatabase.db BottleRocket*/PyFoamRunner.compressibleInterFoam.analyzed/pickledData --create

Then we dump the data to an Excel-file (not that we're going to use that):

pyFoamDumpRunDatabaseToCSV.py runDatabase.db rawRuns.xls --excel --disable=".+Text"  --panda --interactive-after

And get dropped to an ipython-shell:

%pylab
du=self.dump
p=du.pivot("parameters//frac","parameters//p")
p["analyzed//Custom//maxHeight//maximum height"].plot(marker="o")
savefig("heightWithoutDrag.png")

Later we could put that into a script ….

Overview of the variation (maximum height of flight)

Function of initial pressure and filling ratio

Hold on

Over 100 meter?

This can't be right

Fixing the height

What did we forget?

• Check the velocities:
  • They are 50 m/s or more
• The acceleration due to the drag force plays a role

$a_{drag}=-\frac{1}{2} \frac{\rho_{air}}{m(t)} \vec v \left|\vec v\right|A$

Modification of the differential equations

$\dot v=a_{drag}+a_{thrust}-g$

Would make a recalculation of all 190 cases necessary.

But: we already have all the necessary data $a_{thrust}(t)$ and $m(t)$

Enter scipy

Write a script that reads the data

from PyFoam.Applications.RedoPlot import RedoPlot
d=RedoPlot(args=["--pickle-file",
         path.join(d,"PyFoamRunner.compressibleInterFoam.analyzed","pickledPlots"),
         "--numpy"]).plotNumpy

Builds interpolated function from it

from scipy.interpolate import interp1d
from numpy import insert

t=insert(d["bottleAcceleration"]["time"],0,[0])
a=insert(d["bottleAcceleration"]["a"],0,[0])
aFunc=interp1d(t,a,bounds_error=False,fill_value=0)

And integrates the new ODE

from scipy.integrate import odeint
def func(x,t):
    dragForce=0.5*drag*area*x[1]*abs(x[1])/mFunc(t)
    return [x[1],
        aFunc(t)-grav-dragForce]

sol=odeint(func,[0,0],t)

… and put it into a loop over all cases

Plotting the data

• The maximum height is added to the case data

data=Unpickler(open(pickPath)).load()
my={"dragHeight":sol[:,0].max(),
     peakTime":t[sol[:,0].argmax()]}
data["analyzed"]["Custom"]["drag"]=my
p=Pickler(open(pickPath+".addedVel","w"))
p.dump(data)

• And plotted the way the old data was
  • For some cases the ODE did not work: we're filtering this data out

d=self.dump
p=d[d["analyzed//Custom//drag//dragHeight"]<100].pivot("parameters//frac","parameters//p")
p["analyzed//Custom//drag//dragHeight"].plot(marker="o")
plot(p["analyzed//Custom//drag//dragHeight"].idxmax(),p["analyzed//Custom//drag//dragHeight"].max(),'r',linewidth=4)

Heights with drag force

What we learn from this variation

• Higher pressure means higher peak height
  • Big surprise
  • But we could try higher pressures than the recommended 6 bar
• The optimal filling ratio is approximately $0.3$
  • Rises with pressure
  • This is in tune with the DLR recommendation of "fill one third of the bottle with water"

The advantage of scripts

• One case was set up
  • 190 cases were run
• No manual post-processing was required
  • Automatic data collection removes the possibility of human error
• The whole setup can be rerun with a more realistic CAD-bottle (meshed with snappyHexMesh)
  • Minimal setup time, maximal CPU-usage
• PyFoam supports the setup, running and data collection

FoamCommandCenter

Another application of PyFoam

• PyFoam can serve as the basis for more complex applications

• One example is the FoamCommandCenter
  • A in-house application that
    • Lists OpenFOAM-runs in the network (the have to be controlled by PyFoam)
    • Allows checking their state (output and plots)
    • Control the runs (stop them, edit files)
    • … all from the browser
  • It is written with django
    • A Python-framework for Web-applications
  • PyFoam does the OpenFOAM-specific stuff

Overview of runs

One run

Latest solver output

The plots

Goodbye

Material in this presentation

• Some of the demonstrated features are not yet in the released version of PyFoam
  • New release in the next couple of weeks
• The demonstration case is available by request
  • Although you can easily reproduce it

Goodbye

Questions?

This presentation

Was written with

reveal.js

a HTML5-framework for presentation

org-mode

Text outliner mode (and more) in

emacs

Text editor of kings

Bonus: XKCD-style plots always look serious