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.