Detached Eddy Simulation of Flow over a Periodic Dune

Detached Eddy Simulation (DES) of periodic flow (outlet is fed back to the inlet) over a 2-D asymmetric dune. Starting condition (t=0s) is the dune modelled using a Reynolds Averaged Navier-Stokes (RANS) model used to initialize the model. Modelled in OpenFOAM using the Spalart-Allmaras delayed detached eddy simulation (DDES) model. This 13 second simulation took 12 hours to compute running in parallel on 4 CPUs.

Running OpenFoam in Parallel

I cannot believe that I waited so long to try running OpenFOAM in parallel. Perhaps most of my models were small enough not to need parallelization and would likely not have benefited. I am now running a much larger case and run times are going to be significant so the benefits of running in parallel will be much more significant.

Running Environment

I am running OpenFOAM using VirturalBox on a Dell T5500 workstation which has a six CPU Xenon processor  (X5650) with 24 GB of RAM. Lots of power there. I have the virtual machine set up to run 4 CPUs. I am running the GeekoCFD virtual machine created by Alberto Passalacqua which is based on the OpenSUSE 12.3-64bit Linux distribution.

Setting up the Parallel Run

I first copied the decomposeParDict dictionary from the interFoam/ras/dambreak tutorial. I did not modify any of the settings for this first run. I then ran the decomposePar command and promptly received an error "Open RTE was unable to open the hostfile:". I did a quick search and found the following page and implemented the fix. I created a symbolic link:

ln -s /etc /usr/lib64/mpi/gcc/openmpi/etc

I re-ran decomposePar and then ran the command to begin the run:

mpirun -np 4 pisoFoam -parallel >logParallel &

The the model began running.

Generating an OpenFOAM mesh from a point cloud (x,y,z)

I was recently provided a file containing x,y,z points defining the bed elevation extracted from a flume study. I needed a way to efficiently convert this to an OpenFOAM mesh file while having control over the mesh resolution. I am not that familiar with the snappyHexMesh program (yet) but one option would be to convert the point cloud to an STL file. The STL file can then be used to generate a mesh using snappyHexMesh. Note that there is a nice program in Matlab that will perform the conversion from a point cloud to STL for you. Incidentally there is a great tutorial on using snappyHexMesh that can be found here.

Since I am familar with the blockMesh utility and GIS, I decided to try another approach. Using ArcGIS I converted the point cloud first to a TIN and then a DEM. A DEM is essentially a regular matrix of elevations and which makes it very easy to work with as opposed to a randomly spaced point cloud. I then exported the DEM to an ASCII DEM. Finally I wrote a PYTHON script to read the ASCII DEM and write out the blockMeshDict file. The script can be found below. Currently the script requires to user to specify the number of rows and columns as well as the cell size however these could also be read from the ASCII DEM header. The script is nice in that it is easy to make modifications to extract subsets. For example I have a version used to extract a 2D slice of the bed.

There is one issue with the code to be aware of. The code treats each centroid of the DEM as a vertex in the blockMesh which is not quite correct. This may or may not be an issue depending on your problem, but can be corrected. If I get around to correcting this I will post the code.

Here is an image of the mesh I was able to generate for my point cloud:

'''
Created on May 10, 2013

@author: Patrick Grover
'''
from numpy import *

# note that matrix starts at 1,1 and goes up to numrows,numcols
def getvertex(row, col, numrows):
    index = numrows*(col-1) + row
    return index

def getBlock(row, col, numrows,numcols, blockSettings):
    retval =  "\thex ("
    retval = retval + str(getvertex(row,col,numrows)-1) + " "
    retval = retval + str(getvertex(row+1,col,numrows)-1) + " "
    retval = retval + str(getvertex(row+1,col+1,numrows)-1) + " "
    retval = retval + str(getvertex(row,col+1,numrows)-1) + " "    
    
    numvertices = numrows*numcols
    
    retval = retval + str(getvertex(row,col,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row+1,col,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row+1,col+1,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row,col+1,numrows) + numvertices-1) + ") "
    retval = retval + blockSettings     
    return retval

def getTopWall(row, col, numrows,numcols):
    numvertices = numrows*numcols
    retval = "(" + str(getvertex(row,col,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row,col+1,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row+1,col+1,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row+1,col,numrows) + numvertices-1) + ") "
    return retval

def getBottomWall(row, col, numrows,numcols):
    retval = "(" + str(getvertex(row,col,numrows) -1) + " "
    retval = retval + str(getvertex(row+1,col,numrows) -1) + " "
    retval = retval + str(getvertex(row+1,col+1,numrows)-1 ) + " "    
    retval = retval + str(getvertex(row,col+1,numrows)-1 )   + ") "
    return retval

def getInlet(row, col, numrows,numcols):
    numvertices = numrows*numcols
    retval = "(" + str(getvertex(row,col,numrows)-1) + " "
    retval = retval + str(getvertex(row,col,numrows)+ numvertices-1) + " "
    retval = retval + str(getvertex(row+1,col,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row+1,col,numrows) -1)    
    retval = retval + ") "
    return retval

def getOutlet(row, col, numrows,numcols):
    numvertices = numrows*numcols
    retval = "(" + str(getvertex(row,col,numrows)-1) + " "
    retval = retval + str(getvertex(row+1,col,numrows) -1) + " "
    retval = retval + str(getvertex(row+1,col,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row,col,numrows) + numvertices-1)      
    retval = retval + ") "
    return retval

def getLeftWall(row, col, numrows,numcols):
    numvertices = numrows*numcols
    retval = "(" + str(getvertex(row,col,numrows)-1) + " "
    retval = retval + str(getvertex(row,col+1,numrows)-1) + " "
    retval = retval + str(getvertex(row,col+1,numrows) + numvertices-1) + " "  
    retval = retval + str(getvertex(row,col,numrows) + numvertices-1)      
    retval = retval + ") "
    return retval

def getRightWall(row, col, numrows,numcols):
    numvertices = numrows*numcols
    retval = "(" + str(getvertex(row,col,numrows)-1) + " "
    retval = retval + str(getvertex(row,col,numrows)+ numvertices-1) + " "
    retval = retval + str(getvertex(row,col+1,numrows) + numvertices-1) + " "
    retval = retval + str(getvertex(row,col+1,numrows)-1)
    retval = retval + ") "
    return retval


#--------------------------------------------------------------------------------
# BlockMesh Parameters

# Manually set these. Could read these from the ASCII header.
depth = 2.0
numrows = 67
numcols = 1000
cellsize = 0.03   

# Meshing details - describes the mesh
blockSettings = "(2 1 20) simpleGrading (1 1 10)\n"

# Input ASCII DEM

asciiFilePath = '/home/geeko/working/GIS/dunes_3cm_67x1000.asc'
# Path to the blockMeshDict

blockMeshPath = "/home/geeko/working/BlockMesh/TestMesh2/constant/polyMesh/blockMeshDict"

# Read in ASCII file
data = genfromtxt(asciiFilePath,delimiter=' ',skiprows=6,unpack=True, filling_values=-9999)

print data

fo= open(blockMeshPath,'w')
fo.write("FoamFile\n{version\t2.0;\n\tformat\tascii;\n\tclass\tdictionary;\n\tobject\tblockMeshDict;\n}\n\n")

fo.write("convertToMeters 1.0;\n\n")


# Write out vertices
fo.write("vertices\n")
fo.write("(\n")

for y in range(0,numcols):
    for x in range(0,numrows):
        #print "x= " + str(x) + "y= " + str(y)
        if data[y,x]==-9999:
            print "!!!!!!!!!!!!!!!! Error!!!!!!!!!!!!!!!!!!!!!"
        #fo.write("( " + str(x*cellsize) + " " + str(y*cellsize) + " " + str(data[numrows-y-1,x]) + ")\n")
        fo.write("( " + str(x*cellsize) + " " + str(y*cellsize) + " " + str(data[y-1,x]) + ")\n")
for y in range(0,numcols):
    for x in range(0,numrows):
        fo.write("( " + str(x*cellsize) + " " + str(y*cellsize) + " " + str(depth) + ")\n")
    
fo.write(");\n")



# Write out blocks
fo.write("blocks\n")
fo.write("(\n")
for row in range (1,numrows):
    for col in range(1,numcols):
        fo.write(getBlock(row, col, numrows,numcols, blockSettings))
fo.write(");\n")
fo.close

#edges
fo.write("edges\n")
fo.write("(\n")
fo.write(");\n");

# Boundaries
fo.write("boundary\n")
fo.write("(\n")


# Top Wall
fo.write("\ttopWall\n")
fo.write("\t{\n")
fo.write("\ttype symmetryPlane;\n")
fo.write("\tfaces\n")
fo.write("\t(\n")
for row in range (1,numrows):
    for col in range(1,numcols):
        fo.write("\t\t" + getTopWall(row, col, numrows,numcols) + "\n")
fo.write("\t);\n")
fo.write("\t}\n")

# getBottomWall
fo.write("\tbottomWall\n")
fo.write("\t{\n")
fo.write("\ttype wall;\n")
fo.write("\tfaces\n")
fo.write("\t(\n")
for row in range (1,numrows):
    for col in range(1,numcols):
        fo.write("\t\t" + getBottomWall(row, col, numrows,numcols) + "\n")
fo.write("\t);\n")
fo.write("\t}\n")

# Inlet
fo.write("\tinlet\n")
fo.write("\t{\n")
fo.write("\ttype patch;\n")
#fo.write("\ttype cyclic;\n")
#fo.write("\tneighbourPatch outlet;\n")
fo.write("\tfaces\n")
fo.write("\t(\n")
for row in range(1,numrows):
    fo.write("\t\t" + getInlet(row, 1, numrows,numcols) + "\n")
fo.write("\t);\n")
fo.write("\t}\n")

# Outlet
fo.write("\toutlet\n")
fo.write("\t{\n")
fo.write("\ttype patch;\n")
#fo.write("\ttype cyclic;\n")
#fo.write("\tneighbourPatch inlet;\n")
fo.write("\tfaces\n")
fo.write("\t(\n")
for row in range(1,numrows):
    fo.write("\t\t" + getOutlet(row, numcols, numrows,numcols) + "\n")
fo.write("\t);\n")
fo.write("\t}\n")



# Front and Back
fo.write("\tfrontAndBack\n")
fo.write("\t{\n")
fo.write("\ttype empty;\n")
fo.write("\tfaces\n")
fo.write("\t(\n")
fo.write("\t// Left side\n")
for col in range (1,numcols):
    fo.write("\t\t" + getLeftWall(1, col, numrows,numcols) + "\n")
fo.write("\t// Right side\n")
for col in range (1,numcols):
    fo.write("\t\t" + getRightWall(numrows, col, numrows,numcols) + "\n")
fo.write("\t);\n")
fo.write("\t}\n")

fo.write(");\n")


#mergePatchPairs
fo.write("mergePatchPairs\n")
fo.write("(\n")
fo.write(");\n");


print "Done!"
    

Problems with OpenFOAM mapFIelds Utility

While refining a mesh, I was trying to map the internal fields from a courser mesh to a finer mesh using the mapFields utility within OpenFOAM.
FOAM FATAL IO ERROR : size ### is not equal to the given value of ###
This stumped me because all I did was increase the resolution of my mesh, nothing else. Two things corrected the problem. Hirst, I deleted all of the time series folders (e.g. /0, /10, /20...) from the target directory. This fixed part of the problem however I was still getting errors when running mapFields. Looking at the error I noticed that it was having trouble processing the files generated during my post-processing. For example I typically generate the Reynolds Shear Stresses (R) and run the normalized wall spacing, y⁺ (yPlusRAS). I deleted these from the source time-step folder and mapFields ran correctly.

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)

Time-Lapse Removal of Dam

The video below is the controlled removal of the Condit Dam on the White Salmon River. While this is not a proper dam failure, the dynamics of the flow resulting from the breach and the morphological changes in the reservoir are very interesting to watch.