OpenFOAM postprocessing: a few shortcuts

Run-time postprocessing

For pump performance calculation (and in many other cases as well), one usually needs data about pressure rise (head), flow, and power. If cavitation is an issue, maximum and minimum values for pressure come in handy, also maximum and minimum velocity can tell about convergence or numerical issues. All of these values are best calculated at runtime so it’s best to put the settings in controlDict.functionObjects.

// impeller forces for torque/power/viscous losses/axial and radial forces/mechanical losses
forces_impeller
{
    type    forces;
    libs    ("libforces.so");
    patches (impeller_walls);

    rho     rhoInf;
    rhoInf  985; // depends on temperature, also check viscosity in constant/transportProperties

    CofR (0 0 0);
}

// pressure difference between inlet and outlet for head calculation
delta_p
{
    type    fieldValueDelta;
    libs    ("libfieldFunctionObjects.so");
    operation   subtract;

    region1
    {
        type        surfaceFieldValue;
        operation   average;
        fields      (p);
        writeFields no;
        regionType  patch;
        name        outlet;
    }

    region2
    {
        type        surfaceFieldValue;
        operation   average;
        fields      (p);
        writeFields no;
        regionType  patch;
        name        inlet;
    }
}

// flow (is specified in 0/U but here it's added for an additional check)
Q
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");

    writeFields     false;
    surfaceFormat   foam;

    regionType      patch;
    name            outlet;

    operation       areaIntegrate;

    fields
    (
        U
    );
}

// maximum values: maximum velocity should be around impeller circumferential velocity
// difference between minimum pressure and inlet pressure is NPSH
maxValues
{
    type        fieldMinMax;
    libs        ("libfieldFunctionObjects.so");

    writeToFile no;
    log         yes;
    location    no;
    mode        magnitude;
    fields      (U p);
}

Values are calculated after every time step and results are saved in postProcessing/ directory.

Post-postprocessing

There is also a handful of other quantities of interest like yPlus for mesh quality check and Q and vorticity for turbulence/losses visualisation. The problem with these quantities are that the whole field gets saved each timestep. If your writeInterval in controlDict is not 1 OpenFOAM makes a mess with timestep directories. Now every timestep has a directory, some only have yPlus and Q fields and others have all calculated quantities. This makes working with Paraview difficult so it’s better to postprocess these fields separately.

mpirun -np 4 pimpleFoam -parallel -postProcess -latestTime -func yPlus
mpirun -np 4 pimpleFoam -parallel -postProcess -latestTime -func Q
mpirun -np 4 pimpleFoam -parallel -postProcess -latestTime -func vorticity

The postProcessing/ directory

If for any reason the simulation is interrupted and then restarted from latestTime (check your controlDict), new files are added to postProcessing directory. Every added file gets an additional timestep timestamp added to its name. With many restarts the files multiply and make gathering data difficult.

Gathering data

This is a script for reading postprocessing data from all created files mentioned above. It gathers the values in a single Python dictionary and pickles it to a file. It can then be read in Paraview with programmable filter/source for visualization or read from another Python script for plotting and so on.

from subprocess import Popen, PIPE
import os

import re

import pickle

# case times
times = []

# data to be written
data = {
    # [time]: {
    #   delta_p: [number],
    #   Q: [number]
    #   forces_impeller: [F_pressure, F_viscous, M_pressure, M_viscous]
    # }   
}

for t in times:
    data[t] = {}

def read_pp_data(pp_object, parse_function):
    for froot, _, fnames in os.walk(os.path.join('postProcessing', pp_object)):
        for fname in fnames:
            filename = os.path.join(froot, fname)
            if not os.path.isfile(filename):
                continue

            with open(filename, 'r') as f:
                print filename
                
                # skip lines with comment
                for line in f:
                    if line[0] == '#':
                        continue

                    # timestep
                    time = line.split()[0]

                    # add to an existing entry or create a new one
                    if time not in times:
                        times.append(time)
                        data[time] = {}

                    # add data to collection
                    data[time][pp_object] = parse_function(line)

# parsing files: each value has different format
split_re = re.compile('[ \t\n\r()]+') # split: whitespace, parenthesis

def parse_delta_p(line): # delta_p function object
    values = line.split()

    return float(values[1])


def parse_Q(line): # surfaceFieldValue areaIntegrate (flow)
    return float(split_re.split(line)[2])

def parse_forces_impeller(line): # forces functionObject
    # forces format:
    # 0                   1 2 3    4 5 6   7 8 9            10 11 12 13 14 15 16 17 18
    # Time         forces(pressure viscous porous)	moments(pressure viscous porous)
    forces = split_re.split(line)

    # interested in: pressure and viscous forces and moments
    indexes = (1, 2, 3, 4, 5, 6, 10, 11, 12, 13, 14, 15)
    return [float(forces[i]) for i in indexes]
    

read_pp_data('delta_p', parse_delta_p)
read_pp_data('Q', parse_Q)
read_pp_data('forces_impeller', parse_forces_impeller)

# change text numbers to floats
for k in data.keys():
    data[float(k)] = data[k]
    del data[k]

# save the data to be read in paraView (?)
with open('postProcessing.pickle', 'wb') as f:
    pickle.dump(data, f, pickle.HIGHEST_PROTOCOL)

Plotting collected data

Loading a pickled dictionary is a piece of cake. I usually cut the first few results away since timesteps usually don’t converge and therefore values are unphysical.

import pickle
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats

import collect_pp # import the above script to avoid running it manually

# start only after transients decay
t_start = 0.001 # here interesting data starts

# pickled data format:
# 0.000444: {
#   'Q': -0.000222930392,
#   'delta_p': 200.0, 
#   'forces_impeller': [-70.6219451, -0.307142376, 0.116242584, 0.0130158996, -0.0167682449, 0.00424048528, -0.320480923, -0.00105645189, -0.00185268198, -0.0417415834, -4.51203576e-05, -0.00011419315]
# },

with open('postProcessing.pickle', 'rb') as f:
    data = pickle.load(f)

# times
times = data.keys()
times.sort()

n_start = np.searchsorted(np.array(times), t_start)
times = times[n_start:]

p_H = [] # head
p_Q = [] # flow
p_Mp = [] # pressure torque
p_Mv = [] # viscous torque
p_Fx = [] # axial force

for t in times:
    try:
        p_H.append(data[t]['delta_p']) 
        p_Q.append(-data[t]['Q']) 
        Mp.append(-data[t]['forces_impeller'][6]) 
        Mv.append(-data[t]['forces_impeller'][9]) 
        Fx.append(data[t]['forces_impeller'][0]) 
        
        # time
        p_times.append(t)
    except KeyError:
        continue

Filtering stuff

Ocassionally a timestep doesn’t converge and so calculated forces/moments/pressures are way too high or too low. That ruins the day for graphing and these values should be filtered. Actually, it would be better to find out why calculation doesn’t converge in the first place but I haven’t had the time yet for a proper ‘debugging’ session. So I use a little trim utility from scipy:

# filter out the outliers; weird unexplained results that occur randomly
# first, determine which are they. 
def filter_data(values):
    # use trimmed values for average or it doesn't 'mean' anything
    avg = np.average(scipy.stats.trim_mean(values, 0.1))
    
    return np.array(np.where(np.abs(values - avg) > abs(0.5*avg)))[0]

i_filter = np.concatenate((
        filter_data(p_H),
        filter_data(p_Q),
        filter_data(p_Mp),
        filter_data(p_Mv),
        filter_data(p_Fx)))

i_filter = list(set(i_filter))

# delete all data with given indexes
p_times = np.delete(p_times, i_filter)
p_H = np.delete(p_H, i_filter)
p_Q = np.delete(p_Q, i_filter)
p_Mp = np.delete(p_Mp, i_filter)
p_Mv = np.delete(p_Mv, i_filter)
p_Fx = np.delete(p_Fx, i_filter)

Now you’re ready to do whatever with the data. Like plotting:

# plot results
fig, ax1 = plt.subplots()

ax1.plot(p_times, p_H, 'b-')
ax1.set_xlabel('Time (s)')

ax1.set_ylabel('H [dm]', color='b')
ax1.tick_params('y', colors='b')

ax2 = ax1.twinx()

ax2.plot(p_times, p_Mp, 'r.')
ax2.set_ylabel('M [Nm]', color='r')
ax2.tick_params('y', colors='r')

Enjoy!

4 Comments

  1. Very good material, thanks!!!!!!

  2. Thanks for the material. I am getting issues with running delta p. Please have a look at the error. Also find my control dict code below. Any idea on how to resolve the issue.

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    –> FOAM Warning :
    Unknown function type fieldValueDelta

    Valid function types :

    86
    (
    AMIWeights
    BilgerMixtureFraction
    CourantNo
    Curle
    DESModelRegions
    DMD
    LambVector
    Lambda2
    MachNo
    ObukhovLength
    PecletNo
    Q
    XiReactionRate
    add
    blendingFactor
    columnAverage
    components
    continuityError
    ddt
    ddt2
    derivedFields
    div
    enstrophy
    externalCoupled
    extractEulerianParticles
    fieldAverage
    fieldCoordinateSystemTransform
    fieldExtents
    fieldMinMax
    flowType
    flux
    fluxSummary
    forceCoeffs
    forces
    grad
    heatTransferCoeff
    histogram
    interfaceHeight
    limitFields
    log
    mag
    magSqr
    mapFields
    momentum
    momentumError
    multiFieldValue
    multiply
    nearWallFields
    particleDistribution
    patchProbes
    pow
    pressure
    probes
    processorField
    proudmanAcousticPower
    psiReactionThermoMoleFractions
    psiReactionsSensitivityAnalysis
    psiSpecieReactionRates
    randomise
    readFields
    reference
    regionSizeDistribution
    rhoReactionThermoMoleFractions
    rhoReactionsSensitivityAnalysis
    rhoSpecieReactionRates
    setFlow
    sets
    stabilityBlendingFactor
    streamFunction
    streamLine
    subtract
    surfaceDistance
    surfaceFieldValue
    surfaceInterpolate
    surfaces
    turbulenceFields
    valueAverage
    volFieldValue
    vorticity
    wallBoundedStreamLine
    wallHeatFlux
    wallShearStress
    writeCellCentres
    writeCellVolumes
    yPlus
    zeroGradient
    )

    From static Foam::autoPtr<Foam::functionObject> Foam::functionObject::New(const Foam::word&, const Foam::Time&, const Foam::dictionary&)
    in file db/functionObjects/functionObject/functionObject.C at line 117.

    –> loading function object ‘delta_p’

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    /——————————–– C++ –———————————-\
    | ========= | |
    | \ / F ield | OpenFOAM: The Open Source CFD Toolbox |
    | \ / O peration | Version: v2106 |
    | \ / A nd | Website: http://www.openfoam.com |
    | \/ M anipulation | |
    *—————————————————————————*/
    FoamFile
    {
    version 2.0;
    format ascii;
    class dictionary;
    location “system”;
    object controlDict;
    }
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    application simpleFoam;

    startFrom latestTime;

    startTime 0;

    stopAt endTime;

    endTime 1.8;

    deltaT 0.003;

    writeControl timeStep;

    writeInterval 200;//100;

    purgeWrite 0;

    writeFormat ascii;

    writePrecision 6;

    writeCompression off;

    timeFormat general;

    timePrecision 6;

    runTimeModifiable true;

    libs
    (
    “libforces.so” // this is for forces
    );

    functions
    {
    //#includeFunc residuals

    forces
    {
    type forces;
    //functionObjectLibs ("libforces.so"); //Found [v1612] 'functionObjectLibs' entry
    libs ("libforces.so"); //In new
    patches (impeller); // sum the forces and moments on those patches
    writeControl timeStep;
    writeInterval 1;
    p p;
    U U;
    log true;
    rhoInf 998.2;
    rho rhoInf;
    CofR (0 0 0.094289); //centreOfRotation
    }

    wallShearStress
    {
    type wallShearStress;
    libs ("libfieldFunctionObjects.so");
    // Optional entries
    patches (impeller); // sum the forces and moments on those patch
    executeControl timeStep;
    writeControl writeTime;
    }

    // pressure difference between inlet and outlet for head calculation
    delta_p
    {
    type fieldValueDelta;
    libs (fieldFunctionObjects);
    operation subtract;

    region1
    {
    type surfaceFieldValue;
    operation average;
    fields (p);
    writeFields no;
    regionType patch;
    name outlet;
    }

    region2
    {
    type surfaceFieldValue;
    operation average;
    fields (p);
    writeFields no;
    regionType patch;
    name inlet;
    }
    }
    // flow (is specified in 0/U but here it's added for an additional check)
    Q
    {
    type surfaceFieldValue;
    libs ("libfieldFunctionObjects.so");

    writeFields false;
    surfaceFormat foam;

    regionType patch;
    name outlet;

    operation areaIntegrate;

    fields
    (
    U
    );
    }

    // maximum values: maximum velocity should be around impeller circumferential velocity
    // difference between minimum pressure and inlet pressure is NPSH
    maxValues
    {
    type fieldMinMax;
    libs ("libfieldFunctionObjects.so");

    writeToFile no;
    log yes;
    location no;
    mode magnitude;
    fields (U p);
    }

    }

    // ************************************************************************* //

Leave a Reply to Tom AlexCancel reply