Plotting and extracting data along a wall (patch) using ParaView

On thing I often need to do is to examine different flow properties along a wall from my model runs in OpenFOAM. A very simple way to do this is to use ParaView.

Step 1. Open up the model run in ParaView as normal. Under the Properties, Mesh Regions, select the region that you are interested in (e.g. the name of the Patch from defined in your mesh model). Advance the model to the time you are interested in.

Step 2. From the Filters menu, select Alphabetical, Extract Block. In the properties for the filter, select the name of the patch you are interested in (make sure that you selected the patch to be active in Step 1, otherwise it will not appear). Click Apply, only the patch you selected should be visible.

Step 3. You can export the data to a csv file at this point. Make sure that the Extract Block object is selected in the Pipeline Browser and under the Menu select File and Save Data. A Save File dialog box will open allowing you to select the location where to save the data as a CSV file.

Step 4. To plot the data, again have the Extract Block selected in the Pipeline Browser and select Plot Data from the Filters/Alphabetical menu. Click Apply from the Object Inspector and the graph should appear. Click into the graph to make it active. Select the Display tab on the Object Inspector and you should see a number of selections. Under Attribute Mode, I typically select Cell Data as opposed to Point Data. Under the Line Series section, all of the flow variables are available to toggle on/off and select the legend.

Auto-Calibration using PyFoam

I have started using OpenFOAM for the simulation of flows over dunes. The free-surface  is modelled using a symmetry boundary condition or "rigid-lid". Flow over the dune is therefore created by adding a pressure difference between the inlet and outlet boundaries which are periodic (using the fan boundary). Unfortunately this makes a problem for calibrating the models as the pressure difference must be adjusted for each model. Fortunately, I have been playing with pyFoam and wrote some code to auto-calibrate the pressure difference. PyFoam is a great library for automating runs with OpenFOAM and makes life much easier once you get a hang of it. I had some problems running some utilities from PyFoam (using the UtilityRunner class), but just shelled out the commands instead. The code could definitely be improved but it does the trick. Just a quick comment on what the code is doing...The calibration is performed by matching the free surface velocity at the crest of the dune (U0_Target). After the first run of the model, the pressure difference is adjusted by a percent difference in the appropriate direction. On subsequent runs, linear regression using the two previous results is used to make a guess at the correct pressure difference.

#! /usr/bin/env python

# See http://www.openfoamworkshop.org/2009/4th_Workshop/0_Feature_Presentations/OFW4_2009_Gschaider_PyFoam.pdf slide 48

from PyFoam.Execution.UtilityRunner import UtilityRunner
from PyFoam.Execution.BasicRunner import BasicRunner
from PyFoam.RunDictionary.SolutionDirectory import SolutionDirectory
from PyFoam.RunDictionary.ParsedParameterFile import ParsedParameterFile
from os import path
import os
import time
import csv

case = "kOmega_15_dy2_roughwall"
U0_Target = 0.44

dire = SolutionDirectory(case,archive=None,paraviewLink=False)

dp0 = 0;
dp1 = 0;
U0_0 =0;
U0_1 =0;


for run in range(0,5):
    # Clear out the old results
    print '---------------------- Start --------------------------------------'
    print 'We\'re on time %d' % (run)
    print "Clearing out the old results"
    dire.clearResults()

    # Read the pressure
    pressureFile = ParsedParameterFile( path.join( case ,"0/p"))
    POutlet0 = pressureFile["boundaryField"]["outlet"]["f"][0]
    PInlet0 = pressureFile["boundaryField"]["inlet"]["f"][0]

    print "The inlet p is %f", PInlet0
    print "The outlet p is %f", POutlet0

    #pressureFile.close()

    # Run the model
    print "Starting run of model......."
    BasicRunner( argv =["pimpleFoam","-case ", dire.name ], silent=True).start()
    print "Done running model"


    # Sleep for 30 s
    time.sleep(1)

    # Run sample
    #UtilityRunner(argv=["sample", "-case", dire.name], silent=False).start()

    #pUtil=UtilityRunner(argv=["sample",".",dire.name],silent=False,logname="sample")
    #pUtil.start()

    #from subprocess import call
    # Cannot seem to call sample using the UtilityRunner. Just simply 
    sample_str = "sample -case "
    sample_str += case
    os.system(sample_str)
    time.sleep(5)

    # Get the U0 value from L1
    L1UfilePath = path.join( case ,"sets/200","lineL1_U.xy")  
    with open(L1UfilePath, 'rb') as f:
 mycsv = csv.reader(f,delimiter='\t')
 mycsv = list(mycsv)
 U0 = float(mycsv[15][1])
 print "L1 U is %f", U0
    # We check what run we are on in order to imrpove accuracy of iterations beyond 2
    if run == 0:
      U0_0 = U0
      dp0 = POutlet0
      if U0 < U0_Target:
 print "U0 is too low %f, increasing pressure" % (U0)
 dpNew = POutlet0*1.01
      else:
 print "U0 is too high %f, decreasing pressure" % (U0)
 dpNew = POutlet0*0.99
    elif run == 1:
      U0_1 = U0
      dp1 = POutlet0
      if U0 < U0_Target:
 print "U0 is too low %f, increasing pressure" % (U0)
 dpNew = POutlet0*1.01
      else:
 print "U0 is too high %f, decreasing pressure" % (U0)
 dpNew = POutlet0*0.99
    else:
      
      # Update the variables
      U0_0 = U0_1
      dp0 = dp1
      
      U0_1 = U0
      dp1 = POutlet0
      
      dpNew = dp0 + (U0_Target - U0_0)*((dp1-dp0)/(U0_1 - U0_0))
    
 
    # Update the pressure file with the new pressure
    pressureFile = ParsedParameterFile( path.join( case ,"0/p"))
    pressureFile["boundaryField"]["outlet"]["f"][0] = dpNew
    pressureFile["boundaryField"]["inlet"]["f"][0] = dpNew


    pressureFile.writeFile()

    print "Old P ", POutlet0, " new P ", dpNew; 
    
    #pressureFile.closeFile()
    print '---------------------- End Run --------------------------------------'
print "Done!"
os.system("R -case ", case)
os.system("wallShearStress -case ", case)
os.system("yPlusRAS -case ", case)