Before and after the case

New PyFoam Features

or: Comparing solvers without any claim for correctness

Bernhard F.W. Gschaider

Wien, 2014-06-16

"The lesser the people know about how sausages and laws are made, the better they sleep in the night"

Otto v. Bismarck

Bismarck knew nothing about presentations.

This presentation is about how it was made

Who is this Bismarck anyway?

Overview

The aim of this presentation is to demonstrate two new features of PyFoam

  • Semi-automatic case setup
  • Working with IPython-notebooks

Structure of the presentation

  • What is PyFoam
  • Our "project"
  • Setting up the case
  • Getting the results

And who is PyFoam?

  • PyFoam is a Python-library
  • Aimed at working with OpenFOAM
    • Manipulates case files
    • Starts runs and analyzes their output
  • There are utilities implemented with this library
    • This is the part that people usually call "PyFoam"
  • Most parts of PyFoam written by the speaker
    • Which doesn't mean that he knows how to use it properly

Further information about PyFoam

Updates

Are usually announced:

  • On the Wiki-page
  • On the MessageBoard
  • On Twitter @swakPyFoam

The case

p-U-coupled solver

  • The next version of foam-extend has a p-U-coupled solver pUCoupledFoam
    • Sets up one big matrix with pressure and all three velocity components
    • Solves them all at the same time
  • I never worked with it before
    • But people who did are quite excited
    • Supposedly it converges faster and is more stable

The case

  • There are tutorial-cases that come with pUCoupledFoam
    • Never trust tutorials provided with the solver: Their purpose is to make the solver look good
  • Instead a well-known case is used that also comes with (Open)FOAM
    • pitzDaily: a backward facing step for which measurment data by Pitz and Daily exists
In [4]: 
from IPython.display import Image
In [6]: 
Image(filename="pitzDailyOverview.png")
Out[6]:

What we're comparing with

We like both types of music: we like Country and Western

Jake & Elwood Blues

  • There are currently two major forks of the code originally called Foam
    • OpenFOAM by ESI-OpenCFD
    • Foam-extend by Hrvoje Jasak (Wikki Ltd) and community contributors
  • Foam-extend has block-coupled matrixes for some time now
    • Until now the examples were conjugate heat transfer and coupled transport equations
  • OpenFOAM has a similar facility for approximately a year
    • There is now an option to solve the U-Eqn in regular solvers coupled
    • According to release notes current versions have improved numerics

We're going to compare:

  • coupled31 : pUCoupledFoam from foam-extend-3.1
  • segregated31: simpleFoam from foam-extend-3.1
  • segregated23: simpleFoam from OpenFOAM-2.3
  • coupled23: simpleFoam from OpenFOAM-2.3 with the coupled option for the U-Eqn

What we're comparing

  • The base case will be the pitzDaily from foam-extend-3.1
  • For OpenFOAM-2.3 the fvSchemes from OF-2.3 will be used
    • Because this is what is probably the best there
  • For pUCoupledFoam the fvSolution from the tutorial case was used

Numbers that are compared

  • Residuals
    • Are calculated the same way. So they should be comparable
  • Pressure drop between inlet and outlet
  • Detachment point of the circulation behind the step
    • The place on the lower wall where the flow changes direction
  • Profiles defined on specific locations
    • Experimental data is measured from the original paper

Detachment point

In [3]: 
Image(filename='pitzDailyDetachmentPoint.png')
Out[3]:

The measured profiles

In [7]: 
Image(filename="pitzDailyMeasurments.png")
Out[7]:

The problem of slighly different cases

  • The pitzDaily-cases are similar on OpenFOAM-2.3and foam-extend-3.1
    • But not the same
  • Some of the differences are even incompatible
  • One solution:
    • Set up case for one distro
    • Run
    • Comment out stuff that doesn't fit the other distro
    • Run
    • Afterwards change files back
  • Or:
    • Copy case
    • Adapt for other distro
    • If things (BC etc) change: remember to change in both cases consistently
  • Both ways of working are rather tedious
    • Would be nice to have a way to have both in one file
      • Something like #ifdef but more flexible

New utility: pyFoamPrepareCase.py

  • Case setup usually involves these steps:
    1. Getting original files
    2. Editing files (boundary conditions, material parameters)
    3. Creating the mesh
    4. Setting initial conditions
  • The current "standard" Allrun is only good at 3 and 4
    • Unless you use changeDictionary
  • The pyFoamPrepareCase.py utility tries to automatize this
    • Based on templates

Workflow of pyFoamPrepareCase.py

  1. Clear old data from the case
    • If a folder 0.orgis present remove that too
  2. For every file with the extension .template do template expansion
  3. Create a mesh
    • If a file meshCreate.sh is found execute this
    • Otherwise execute blockMesh
  4. Copy every foo.org (file or directory) to foo
  5. Do template expansion for every file with the extension .postTemplate
  6. Execute caseSetup.sh if present

The scripts

That is the easy part: the case is so simple we don't need no scripts

And you know how to write scripts, right?

Template expansion

  • Template expansion has been in PyFoam for a long time
  • Big leap forward was the inclusion of pyratemp(http://www.simple-is-better.org/template/pyratemp.html)
    • A light-weight templating library
      • That's why it is part of PyFoam- no need to install
    • Allows loops and branches
    • The version in PyFoam is modified and allows things that on the original application field (webservers) would be insanely dangerous
  • PyFoam-specific extensions
    • Text between |- and -| is evaluated as a Python-expression and the result is inserted
    • Lines starting with $$ are discarded
      • Rest of the line is executed as a Python-command (usually variable declarations)

Inline expansion

  • The pitzDaily-case that comes with OpenFOAM specifies $10 \frac{m}{s}$ as the inlet velocity
    • The Pitz-Daily paper says it's $13.3 \frac{m}{s}$
    • We want to be flexible

In 0.org/U.template:

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (|- UIn -| 0 0);
    }

is expanded to (assuming UIn is $13.3$)

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (13.3 0 0);
    }

Branches and loops

  • pyratemp-specific commands start with <!--( and end with )-->
    • Similar to HTML-comments
  • Between these tags is a "Python-like" command
  • <!--(if expr)--> starts a block that is kept if expr is True
    • elif and else are also possible
  • Block ends with <!--(end)-->

For details see documentation

Branches

In fvSolution.template incompatible relaxation factors are set:

relaxationFactors
{
<!--(if pUCoupled and not of23)-->
    U               0.9;
    k               0.95;
    epsilon         0.95;
<!--(elif not of23)-->
    p               0.3;
    U               0.7;
    k               0.7;
    epsilon         0.7;
    R               0.7;
    nuTilda         0.7;
<!--(else)-->
    fields
    {
        p               0.3;
    }
    equations
    {
        U               0.7;
        k               0.7;
        epsilon         0.7;
        R               0.7;
        nuTilda         0.7;
    }
<!--(end)-->
}

pUCoupledFoam expansion

If pUCoupled is set to True and of23 is False this becomes:

relaxationFactors
{
    U               0.9;
    k               0.95;
    epsilon         0.95;
}

Where come the parameters from?

  • Parameters are specified when calling pyFoamPrepareCase.py
    • As a string with values with the option --values-string
    • In a parameter file specified with --parameter-file
      • Parameter-files are regular OpenFOAM-dictionaries
  • These options can be specified multiple times
    • Last value wins

For example parameters.coupled31:

#include "parameters.basic"

pUCoupled on;

of23 off;

This parameter-set is used with

pyFoamPrepareCase.py . --parameter-file=parameteres.segregated31

Other preparations

  • Add an expression for the pressure drop
    • With swak4Foam
  • Add an expression for the detachment point
    • Again: swak4Foam
  • Add sample lines on the places specified

Running the cases

  • Now we run the cases
  • Preferably with pyFoamPlotRunner.py

Post-processing in IPython

  • OpenFOAM produces two kinds of data sets
    • Field data
      • Paraview is quite OK for this
    • Boring numbers from evaluations, samples etc
      • Paraview is not so good for this
      • Excel sucks by definition
      • Gnuplot is OK but cryptic for some

Is there an OpenSource-solution for this?

What is IPython

  • python is a programming language
    • Interpreted
    • Simple syntax
    • Wide range of numeric software
      • numpy for matrices
      • sympy for symbolic compuations
      • scipy for numerical algorithms
      • matplotlib for plotting
      • pandas for convenient handling of data-sets
  • ipython is a shell for python
    • Autocompletion, debugging-support, etc
  • ipython-notebooks are a browser-based interface to this
    • Cell-based:
      • Enter one input-cell
      • Directly below it is an output cell
    • Allows annotations
      • Inclusion of images
      • Equation conversions
    • Newer versions allow GUI-elements
  • All this together is a computing interface similar to Mathematica-notebooks or MATLAB
  • You can even write presentations in it

IPython and PyFoam

  • It was always possible to call PyFoam from an ipython-notebook
    • Even import data for pandas
  • But it was tedious:
    • Long library command calls
    • Frequently changing back to the shell ("What is the name of the sample-file again?")
    • Data-sets with different resolutions had to be resampled to be compared

New features

  • Wrapper for the pandas DataFrame that adds convenience features
    • Resampling
    • Weighted averaging
  • A Case-class for easily calling functionality to import data
    • Also adds info about the case to the notebook
  • A Notebook class to write notebooks from a python-program
    • Also allows storing data into the notebook
      • This is experimental and has performance issues

All this is integrated into one utility

pyFoamIPythonNotebook.py

  • Main application (it knows other tricks too)
    • Create an IPython-notebook based on an (Open)Foam-case
      • Only a starter. The user must "fill it with life"
      • But most repetitive work is already done

In the notebook

  • Most important imports to make working easy
    • Commented what they do
  • Case-object that points to the current case
    • Basic reports about the case
  • GUI-elements to select data that was found
    • Data is stored in Python-variables
      • may be stored in the notebook as well

Demo

Now this is the moment where Bernhard says "I'll change to the shell and show you"

He'll show:

  • How the notebook is set up
  • Opening it
  • Some simple commands

If you read this presentation on the internet: wait 4 minutes before continuing

Evaluations

This section shows some example evaluation. First we set up the environment (remember: most of this is done by pyIPythonNotebook.pyfor you)

In [14]: 
%matplotlib inline
import matplotlib
matplotlib.rcParams['figure.figsize'] = (10.0, 8.0)
from PyFoam.Wrappers.Pandas import PyFoamDataFrame as DataFrame
from PyFoam.IPython.Case import Case

Now we get us the four cases

In [15]: 
case=Case('/Volumes/Foam/Cases/Pfau2014/pitzDaily_coupled31')
case23=Case('/Volumes/Foam/Cases/Pfau2014/pitzDaily_coupled23')
caseSeg=Case('/Volumes/Foam/Cases/Pfau2014/pitzDaily_segregated31')
caseSeg23=Case('/Volumes/Foam/Cases/Pfau2014/pitzDaily_segregated23')

About a case

In [16]: 
case.size(level=3)
Out[16]:

Size of the case

Faces 49180
Points 25012
Cells 12225
In [17]: 
case.boundaryConditions(level=3)
Out[17]:

Table of boundary conditions for t = 0

frontAndBack inlet lowerWall outlet upperWall
Patch Type empty patch wall patch wall
Length 24450 30 250 57 223
R empty fixedValue zeroGradient zeroGradient zeroGradient
U empty fixedValue fixedValue zeroGradient fixedValue
epsilon empty fixedValue epsilonWallFunction zeroGradient epsilonWallFunction
epsilon.old empty fixedValue epsilonWallFunction zeroGradient epsilonWallFunction
k empty fixedValue kqRWallFunction zeroGradient kqRWallFunction
k.old empty fixedValue kqRWallFunction zeroGradient kqRWallFunction
nuTilda empty fixedValue zeroGradient zeroGradient zeroGradient
nut empty calculated nutWallFunction calculated nutWallFunction
p empty zeroGradient zeroGradient fixedValue zeroGradient

Reading the data

In [18]: 
dataCouple=case.pickledPlots('PyFoamRunner.pUCoupledFoam.analyzed/pickledPlots')
In [19]: 
dataCouple23=case23.pickledPlots('PyFoamRunner.simpleFoam.analyzed/pickledPlots')
In [20]: 
dataSeg=caseSeg.pickledPlots('PyFoamRunner.simpleFoam.analyzed/pickledPlots')
In [21]: 
dataSeg23=caseSeg23.pickledPlots('PyFoamRunner.simpleFoam.analyzed/pickledPlots')

Which datasets are there?

In [22]: 
dataCouple.keys()
Out[22]:
['continuity', 'deltaP', 'linear', 'detachement', 'iterations', 'execution']

Simple plotting

In [23]: 
dataCouple["iterations"].plot()
Out[23]:
<matplotlib.axes.AxesSubplot at 0x10d19aad0>

Comparing residuals

This looks worse than it is (renaming datasets for merging):

In [24]: 
linCouple=DataFrame(dataCouple["linear"].select(lambda k:(k.find("_")<0 and k!="Up[2]"),1))
In [25]: 
linCouple23=DataFrame(dataCouple23["linear"].select(lambda k:(k.find("_")<0 and k.find("Uz")<0),1).rename(columns=lambda n:n+"_cou23"))
In [26]: 
linSeg=DataFrame(dataSeg["linear"].select(lambda k:(k.find("_")<0),1).rename(columns=lambda n:n+"_seg"))
In [27]: 
linSeg23=DataFrame(dataSeg23["linear"].select(lambda k:(k.find("_")<0),1).rename(columns=lambda n:n+"_seg23"))

Comparing coupled and segregated

In [28]: 
linSeg.addData(linCouple).plot(logy=True)
Out[28]:
<matplotlib.axes.AxesSubplot at 0x10d15ba50>

Lets have a closer look:

In [29]: 
linSeg.addData(linCouple).plot(logy=True,xlim=(0,400))
Out[29]:
<matplotlib.axes.AxesSubplot at 0x10d0e55d0>

Seems like coupled converges really faster

Comparing OpenFOAM-2.3 and foam-extend-3.1 (segregated)

In [30]: 
linSeg.addData(linSeg23).plot(logy=True)
Out[30]:
<matplotlib.axes.AxesSubplot at 0x10d0c6f10>

I don't dare to say who is winning

Comparing coupled and segregated for 2.3

In [31]: 
linSeg23.addData(linCouple23).plot(logy=True)
Out[31]:
<matplotlib.axes.AxesSubplot at 0x10cfee2d0>

There is a difference

It converges fast, but is it fast?

Caution: These numbers were not gathered under serious benchmarking conditions (there were other things going on on that machine. So the numbers may be skewed)

The calculation times per second are accumulated to get the total calculation time:

In [32]: 
DataFrame(dataCouple["execution"].rename(columns=lambda n:n+" coupled")).addData(
   dataSeg["execution"].rename(columns=lambda n:n+" segregated")).cumsum().plot()
Out[32]:
<matplotlib.axes.AxesSubplot at 0x10d00ab50>

This looks bad for coupled but he converged near $t=550$ and segregated is not nearly finished.

Solution properties

Let's see how the properties of the flow converge

Detachment point

Clean unfit data:

In [33]: 
dataSeg["detachement"][dataSeg["detachement"]<-1000]=float("NaN")

Plot:

In [34]: 
DataFrame(dataCouple["detachement"].rename(columns={"x":"coupled"})).addData(
   dataSeg["detachement"].rename(columns={"x":"segregated"})).plot()
Out[34]:
<matplotlib.axes.AxesSubplot at 0x10d054a50>

Seems that converged long before the residuals

Pressure difference

Renaming data for merging

In [36]: 
deltaPData=DataFrame(dataCouple["deltaP"].rename(columns={"avg":"coupled"})).addData(
   dataSeg["deltaP"].rename(columns={"avg":"segregated"}))
In [37]: 
deltaPData.plot()
Out[37]:
<matplotlib.axes.AxesSubplot at 0x10d059590>

Strong oscillations in the beginning. Try without them:

In [38]: 
deltaPData.plot(xlim=(100,1000),ylim=(-20,0))
Out[38]:
<matplotlib.axes.AxesSubplot at 0x10d9e2090>

This too converges faster

But how big are the oscillations?

In [39]: 
deltaPData.describe()
Out[39]:
coupled segregated
count 999.000000 999.000000
mean -9.092461 -5.414749
std 7.353295 309.862965
min -167.096000 -4124.230000
25% -9.077380 -10.544450
50% -9.075930 -9.775140
75% -9.075830 -8.965875
max 117.733000 6731.340000
integral -9137.696661 -8770.111473
valid length 998.000000 998.000000
weighted average -9.156009 -8.787687

11 rows × 2 columns

More than one order of magnitude smaller for coupled!

Comparing profiles

Reading the profile data from the experiment:

In [41]: 
refU1=case.sampleTime('./experimentReference','experiment1','0')
refU3=case.sampleTime('./experimentReference','experiment3','0')
refU6=case.sampleTime('./experimentReference','experiment6','0')

Reading the sample data from the cases

In [51]: 
sample1Couple=case.sampleTime('./experimentSample','experiment1','1000')
sample3Couple=case.sampleTime('./experimentSample','experiment3','1000')
sample6Couple=case.sampleTime('./experimentSample','experiment6','1000')
sample1Couple=sample1Couple.rename(columns=lambda n:n.replace("experiment1_t=1000","coupled"))
sample3Couple=sample3Couple.rename(columns=lambda n:n.replace("experiment3_t=1000","coupled"))
sample6Couple=sample6Couple.rename(columns=lambda n:n.replace("experiment6_t=1000","coupled"))
In [52]: 
sample1Seg=caseSeg.sampleTime('./experimentSample','experiment1','1000')
sample3Seg=caseSeg.sampleTime('./experimentSample','experiment3','1000')
sample6Seg=caseSeg.sampleTime('./experimentSample','experiment6','1000')
sample1Seg=sample1Seg.rename(columns=lambda n:n.replace("experiment1_t=1000","segregated"))
sample3Seg=sample3Seg.rename(columns=lambda n:n.replace("experiment3_t=1000","segregated"))
sample6Seg=sample6Seg.rename(columns=lambda n:n.replace("experiment6_t=1000","segregated"))
In [53]: 
sample1Seg23=caseSeg23.sampleTime('./postProcessing/experimentSample','experiment1','1000')
sample3Seg23=caseSeg23.sampleTime('./postProcessing/experimentSample','experiment3','1000')
sample6Seg23=caseSeg23.sampleTime('./postProcessing/experimentSample','experiment6','1000')
sample1Seg23=sample1Seg23.rename(columns=lambda n:n.replace("experiment1_t=1000","segregated 2.3"))
sample3Seg23=sample3Seg23.rename(columns=lambda n:n.replace("experiment3_t=1000","segregated 2.3"))
sample6Seg23=sample6Seg23.rename(columns=lambda n:n.replace("experiment6_t=1000","segregated 2.3"))

The sample plots

In [48]: 
DataFrame(DataFrame(DataFrame(sample1Couple["coupled magU"]).addData(
   refU1,mergeIndex=True)).addData(sample1Seg["segregated magU"],mergeIndex=True)).addData(
sample1Seg23["segregated 2.3 magU"],mergeIndex=True).plot()
Out[48]:
<matplotlib.axes.AxesSubplot at 0x10db3e150>
In [54]: 
DataFrame(DataFrame(DataFrame(sample3Couple["coupled magU"]).addData(
   refU3,mergeIndex=True)).addData(sample3Seg["segregated magU"],mergeIndex=True)).addData(
sample3Seg23["segregated 2.3 magU"],mergeIndex=True).plot()
Out[54]:
<matplotlib.axes.AxesSubplot at 0x10dec8050>
In [55]: 
DataFrame(DataFrame(DataFrame(sample6Couple["coupled magU"]).addData(
   refU6,mergeIndex=True)).addData(sample6Seg["segregated magU"],mergeIndex=True)).addData(
sample6Seg23["segregated 2.3 magU"],mergeIndex=True).plot()
Out[55]:
<matplotlib.axes.AxesSubplot at 0x10dcbc950>

Comparing the simulation results

Coupled vs segregated

In [56]: 
(sample1Couple["coupled magU"]-sample1Seg["segregated magU"]).plot()
Out[56]:
<matplotlib.axes.AxesSubplot at 0x10deed450>

2.3 vs 3.1

In [58]: 
(sample1Seg23["segregated 2.3 magU"]-sample1Seg["segregated magU"]).plot()
Out[58]:
<matplotlib.axes.AxesSubplot at 0x10df38c50>

Conclusion on the sample data

  • All 4 solvers show similar results
  • Results are quite different from the original data
    • Possible reasons:
      • Data only very roughly measured
      • Viscosity of fluid not know

Conclusion

  • Two new utilities for PyFoam were demonstrated
    • One to automatically set up parameterized cases
    • The other to create notebooks for the evaluation of OpenFOAM-cases
  • Use of IPython-notebooks as a platform for numeric evaluations
    • Well suited to document the work
    • Maybe not the best tool to write presentations
      • But very nice to do commented evaluations
  • The pUCoupledFoam solver was compared with simpleFoam
    • Results are only preliminary
    • Seems to converge much faster
      • How this translates into CPU-time needs to be investigated further

Where to get it

Errata

  • During the presentation it was suggested from the audience to use the inlet width $25.4 mm$ as $D$ and calculate the actual viscosity from the Reynolds-number provided
    • Using that the viscosity is 50% higher than in the original tutorials
  • The presentation has not been updated to reflect that
In []: