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.
Tuesday, July 23, 2013
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.
Friday, June 21, 2013
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.
Thursday, June 20, 2013
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:
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.
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 (
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.
Thursday, September 20, 2012
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.
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.
Friday, September 7, 2012
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)
Thursday, January 5, 2012
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.
Friday, December 30, 2011
Telemac Benchmark Test
I recently took an old desktop machine of mine and installed GeekoCFD which is an openSuse based distribution with a number of open-source CFD tools including OpenFOAM. I then downloaded and compiled TELEMAC v6p1. Since I last compiled TELEMAC, there is a new Python-based method which does not require Perl. The installation instructions can be found here and are fairly complete. My compilation only used the scalar version of Telemac and not the parallel version. I will be attempting to compile a parallel version in the near future.
I then ran a Telemac3d case to compare the Linux version against my Windows version. My windows machine is a laptop with an Intel Core i3-350M processor (dual core running at 2.26 GHz) with 4 GB of DDR3 memory running 64 bit Windows 7. It was also compiled using the Intel Fortran compiler. My old desktop is running openSuse 11.4 with the KDE desktop (very nice!). It has an Intel Core2 Duo E4500 running at 2.2 GHz and 1 GB of DDR memory. I ran a Telemac3D case and it took 2hrs and 32 minutes on my old desktop and 1 hour on my laptop.
Interesting the difference in the run time difference, as both machines are running dual core processors at nearly the same clock speed - however the newer i3 must be more optimized. The other large difference is that the Linux only has 1GB of memory and also is using the gFortran compiler.
I then ran a Telemac3d case to compare the Linux version against my Windows version. My windows machine is a laptop with an Intel Core i3-350M processor (dual core running at 2.26 GHz) with 4 GB of DDR3 memory running 64 bit Windows 7. It was also compiled using the Intel Fortran compiler. My old desktop is running openSuse 11.4 with the KDE desktop (very nice!). It has an Intel Core2 Duo E4500 running at 2.2 GHz and 1 GB of DDR memory. I ran a Telemac3D case and it took 2hrs and 32 minutes on my old desktop and 1 hour on my laptop.
Interesting the difference in the run time difference, as both machines are running dual core processors at nearly the same clock speed - however the newer i3 must be more optimized. The other large difference is that the Linux only has 1GB of memory and also is using the gFortran compiler.
My next experiment will be to try using Amazon EC2 and sign out a high CPU server, with the eventual goal of being able to cluster a number of instances in the cloud.
Friday, September 30, 2011
Experimenting with OpenFOAM
OpenFOAM is an open source CFD package written in c++ that I have been interested in trying out for some time. I finally sat down and tried my first model using it - a sort of "hello world" example. My goal was to set up my simple dune geometry (dune of 400mm length and 20mm height)and run a simple flow over it. My objective was just to better understand how to define a model geometry in OpenFOAM and the basics for setting up a model run.
OpenFOAM has excellent documentation and also includes a complete set of tutorial models that are great places to start modelling from. I decided to generate this model using the icoFoam solver which solves the incompressible laminar Navier-Stokes equation (more here). Again, in reality this case will be turbulent, but icoFoam is much simpler to set up than a turbulent model.
First step was to copy over and set up the required folder structure. I copied the files from the Cavity example which is one of the tutorials under the icoFoam tutorial directory.
Defining the Mesh
OpenFOAM does not have a GUI which can make things difficult and adds to the learning curve. This is particularly true for building the mesh. For this example I decided to use the BlockMesh mesh generation utility which creates parametric meshes (see here for a further discussion on setting up the mesh).
Before I entered in my mesh I drew it out on paper. This was a key first step since there is no GUI and it can be difficult to visualize what you are doing otherwise. To define my mesh, I first determined all of the vertexes in my problem.
vertices
(
(0 20 0 ) //Vertex 0
(10 20 0 ) //Vertex 1
(10 150 0 ) //Vertex 2
(0 150 0 ) //Vertex 3
(0 20 10 ) //Vertex 4
(10 20 10 ) //Vertex 5
(10 150 10 ) //Vertex 6
(0 150 10 ) //Vertex 7
...
Next, I broke down my model geometry into "blocks" with each block containing 8 vertices. For each block you also define the mesh spacing:
hex (0 1 2 3 4 5 6 7) (5 20 1) simpleGrading (1 1 1)
In total I had 6 blocks that defined my model geometry. Once the blocks are defined, I need to identify patches (boundaries). Patches are defined by a collection of block-faces (e.g. 4 vertices) and have a defined type (e.g. wall or patch). This can be a little confusing - i suggest following through the tutorials in the OpenFoam documentation which gives a better understanding of this - once you get the hang of it, it is quite simple. The patch for the top of my model is given by (note, the name "topWall" is defined by the user):
wall topWall
(
(3 7 6 2)// watch the order of these RHR
(2 6 11 9)
(9 11 15 13)
(13 15 19 17)
(17 19 23 21)
(21 23 27 25)
(25 27 31 29)
)
The definition for the inlet to my model is as follows - note it is a type patch as opposed to a wall :
patch inlet
(
(0 4 7 3)
)
The most important thing to remember is that the vertices that define the boundary need to be in a specific order. The order is based on the right-hand-rule - using your right hand curl your fingers in the direction that you ordered the vertex, your thumb should be pointing out of the model. My full blockMeshDict file looks like this:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 0.001;
vertices
(
(0 20 0 ) //0
(10 20 0 )
(10 150 0 )
(0 150 0 ) //3
(0 20 10 )
(10 20 10 )
(10 150 10 )
(0 150 10 ) //7
(50 0 0 ) //8
(50 150 0 )
(50 0 10 )
(50 150 10 ) //11
(75 0 0 ) // 12
(75 150 0 )
(75 0 10 )
(75 150 10 ) //15
(137.5 2 0 )//16
(137.5 150 0 )
(137.5 2 10 )
(137.5 150 10 )//19
(325 18.40 0 )// 20
(325 150 0 )
(325 18.40 10 )
(325 150 10 )// 23
(390 20 0 ) //24
(390 150 0 )
(390 20 10 )
(390 150 10 ) //27
(400 20 0 )// 28
(400 150 0 )
(400 20 10 )
(400 150 10 )// 31
);
blocks
(
hex (0 1 2 3 4 5 6 7) (5 20 1) simpleGrading (1 1 1) // dx dy dz
hex (1 8 9 2 5 10 11 6) (30 20 1) simpleGrading (1 1 1)
hex (8 12 13 9 10 14 15 11) (20 20 1) simpleGrading (1 1 1)
hex (12 16 17 13 14 18 19 15) (20 20 1) simpleGrading (1 1 1)
hex (16 20 21 17 18 22 23 19) (40 20 1) simpleGrading (1 1 1)
hex (20 24 25 21 22 26 27 23) (10 20 1) simpleGrading (1 1 1)
hex (24 28 29 25 26 30 31 27) (5 20 1) simpleGrading (1 1 1)
);
edges
(
);
patches
(
wall topWall
(
(3 7 6 2)// watch the order of these RHR
(2 6 11 9)
(9 11 15 13)
(13 15 19 17)
(17 19 23 21)
(21 23 27 25)
(25 27 31 29)
)
wall bottomWall
(
(0 1 5 4)// watch the order of these RHR
(1 8 10 5)
(8 12 14 10)
(12 16 18 14)
(16 20 22 18)
(20 24 26 22)
(24 28 30 26)
)
patch inlet // note these are patches!!
(
(0 4 7 3)
)
patch outlet
(
(28 29 31 30)
)
empty frontAndBack
(
(0 3 2 1)// watch the order of these RHR
(4 5 6 7)
(1 2 9 8)
(5 10 11 6)
(8 9 13 12)
(10 14 15 11)
(12 13 17 16)
(14 18 19 15)
(16 17 21 20)
(18 22 23 19)
(20 21 25 24)
(22 26 27 23)
(24 25 29 28)
(26 30 31 27)
)
);
mergePatchPairs
(
);
// ************************************************************************* //
Once defined, I ran the blockMesh utility which generated the mesh. Just another note. When I was generating the mesh, I defined it using one block at a time and ran the meshing utility which made it much easier to identify problems. I then viewed the mesh in ParaView:
My mesh looks ok, however there is still work to be done to adjust the mesh. Once the mesh was built, I set up the initial/boundary conditions. For icoFaom, the velocity and pressure fields need to be defined. I used a simple inlet velocity and defined the outlet pressure.For all of the walls I defined them as being no-slip including the top wall (I will need to look up how to simulate an open channel type flow). Since this is a 2D case, the front and back walls are defined as being empty.
// * * * * * * * Velocity Field * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.5 0 0);
}
outlet
{
type zeroGradient;
}
topWall
{
type fixedValue;
value uniform (0 0 0);
}
bottomWall
{
type fixedValue;
value uniform (0 0 0);
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * Pressure * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
topWall
{
type zeroGradient;
}
bottomWall
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //
I used the model run settings (defined in the controlDict file) the same as for the cavity case. Other than that there was little else to do. I ran the model and viewed the results in ParaView:
Again, I an not interested so much in the results of the model, but a good first step towards understanding how to set up a geometry and run the model.
OpenFOAM has excellent documentation and also includes a complete set of tutorial models that are great places to start modelling from. I decided to generate this model using the icoFoam solver which solves the incompressible laminar Navier-Stokes equation (more here). Again, in reality this case will be turbulent, but icoFoam is much simpler to set up than a turbulent model.
First step was to copy over and set up the required folder structure. I copied the files from the Cavity example which is one of the tutorials under the icoFoam tutorial directory.
Defining the Mesh
OpenFOAM does not have a GUI which can make things difficult and adds to the learning curve. This is particularly true for building the mesh. For this example I decided to use the BlockMesh mesh generation utility which creates parametric meshes (see here for a further discussion on setting up the mesh).
Before I entered in my mesh I drew it out on paper. This was a key first step since there is no GUI and it can be difficult to visualize what you are doing otherwise. To define my mesh, I first determined all of the vertexes in my problem.
vertices
(
(0 20 0 ) //Vertex 0
(10 20 0 ) //Vertex 1
(10 150 0 ) //Vertex 2
(0 150 0 ) //Vertex 3
(0 20 10 ) //Vertex 4
(10 20 10 ) //Vertex 5
(10 150 10 ) //Vertex 6
(0 150 10 ) //Vertex 7
...
Next, I broke down my model geometry into "blocks" with each block containing 8 vertices. For each block you also define the mesh spacing:
hex (0 1 2 3 4 5 6 7) (5 20 1) simpleGrading (1 1 1)
In total I had 6 blocks that defined my model geometry. Once the blocks are defined, I need to identify patches (boundaries). Patches are defined by a collection of block-faces (e.g. 4 vertices) and have a defined type (e.g. wall or patch). This can be a little confusing - i suggest following through the tutorials in the OpenFoam documentation which gives a better understanding of this - once you get the hang of it, it is quite simple. The patch for the top of my model is given by (note, the name "topWall" is defined by the user):
wall topWall
(
(3 7 6 2)// watch the order of these RHR
(2 6 11 9)
(9 11 15 13)
(13 15 19 17)
(17 19 23 21)
(21 23 27 25)
(25 27 31 29)
)
The definition for the inlet to my model is as follows - note it is a type patch as opposed to a wall :
patch inlet
(
(0 4 7 3)
)
The most important thing to remember is that the vertices that define the boundary need to be in a specific order. The order is based on the right-hand-rule - using your right hand curl your fingers in the direction that you ordered the vertex, your thumb should be pointing out of the model. My full blockMeshDict file looks like this:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 0.001;
vertices
(
(0 20 0 ) //0
(10 20 0 )
(10 150 0 )
(0 150 0 ) //3
(0 20 10 )
(10 20 10 )
(10 150 10 )
(0 150 10 ) //7
(50 0 0 ) //8
(50 150 0 )
(50 0 10 )
(50 150 10 ) //11
(75 0 0 ) // 12
(75 150 0 )
(75 0 10 )
(75 150 10 ) //15
(137.5 2 0 )//16
(137.5 150 0 )
(137.5 2 10 )
(137.5 150 10 )//19
(325 18.40 0 )// 20
(325 150 0 )
(325 18.40 10 )
(325 150 10 )// 23
(390 20 0 ) //24
(390 150 0 )
(390 20 10 )
(390 150 10 ) //27
(400 20 0 )// 28
(400 150 0 )
(400 20 10 )
(400 150 10 )// 31
);
blocks
(
hex (0 1 2 3 4 5 6 7) (5 20 1) simpleGrading (1 1 1) // dx dy dz
hex (1 8 9 2 5 10 11 6) (30 20 1) simpleGrading (1 1 1)
hex (8 12 13 9 10 14 15 11) (20 20 1) simpleGrading (1 1 1)
hex (12 16 17 13 14 18 19 15) (20 20 1) simpleGrading (1 1 1)
hex (16 20 21 17 18 22 23 19) (40 20 1) simpleGrading (1 1 1)
hex (20 24 25 21 22 26 27 23) (10 20 1) simpleGrading (1 1 1)
hex (24 28 29 25 26 30 31 27) (5 20 1) simpleGrading (1 1 1)
);
edges
(
);
patches
(
wall topWall
(
(3 7 6 2)// watch the order of these RHR
(2 6 11 9)
(9 11 15 13)
(13 15 19 17)
(17 19 23 21)
(21 23 27 25)
(25 27 31 29)
)
wall bottomWall
(
(0 1 5 4)// watch the order of these RHR
(1 8 10 5)
(8 12 14 10)
(12 16 18 14)
(16 20 22 18)
(20 24 26 22)
(24 28 30 26)
)
patch inlet // note these are patches!!
(
(0 4 7 3)
)
patch outlet
(
(28 29 31 30)
)
empty frontAndBack
(
(0 3 2 1)// watch the order of these RHR
(4 5 6 7)
(1 2 9 8)
(5 10 11 6)
(8 9 13 12)
(10 14 15 11)
(12 13 17 16)
(14 18 19 15)
(16 17 21 20)
(18 22 23 19)
(20 21 25 24)
(22 26 27 23)
(24 25 29 28)
(26 30 31 27)
)
);
mergePatchPairs
(
);
// ************************************************************************* //
My mesh looks ok, however there is still work to be done to adjust the mesh. Once the mesh was built, I set up the initial/boundary conditions. For icoFaom, the velocity and pressure fields need to be defined. I used a simple inlet velocity and defined the outlet pressure.For all of the walls I defined them as being no-slip including the top wall (I will need to look up how to simulate an open channel type flow). Since this is a 2D case, the front and back walls are defined as being empty.
// * * * * * * * Velocity Field * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.5 0 0);
}
outlet
{
type zeroGradient;
}
topWall
{
type fixedValue;
value uniform (0 0 0);
}
bottomWall
{
type fixedValue;
value uniform (0 0 0);
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.5 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * Pressure * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
topWall
{
type zeroGradient;
}
bottomWall
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //
I used the model run settings (defined in the controlDict file) the same as for the cavity case. Other than that there was little else to do. I ran the model and viewed the results in ParaView:
 |
| Velocity Field |
 |
| Pressure Field |
 |
| Vertical Velocity (v) |
Wednesday, August 31, 2011
Structure of the TELEMAC Selafin 3D File
Blue Kenue is a great tool useful for pre and post processing TELEMAC models. However there are times were you want to do some more in depth analysis. I have been looking into developing velocity profiles from my model runs, something which is not possible from Blue Kenue.(EDIT: Blue Kenue does allow the extraction of a vertical profile, simply select a node of your velocity grid (in a 2D view is best) and to bring up the context menu and choose “Extract Vertical Profile” - Thanks Martin Serrer for alerting me to that!) It does have tools to export model runs to ASCII files however some of these tools are limited. There was this discussed on the Open Telemac Forum related to his issue that got me started.
Fortan
The first thing that you will need to do is extract your model results from the Selafin file. The Selafin file is in Binary format which makes things a little complex. The structure of a Selafin file is well described in Appendix 3 of the TELEMAC-2D User Manual. I have not programmed in Fortran since 1990, but I thought I would give it a try.
PROGRAM ReadSelaphin
! http://www.opentelemac.org/index.php?option=com_kunena&func=view&catid=20&id=1162&Itemid=62&lang=en
INTEGER IERR
CHARACTER(LEN=80) :: strTitle
CHARACTER(LEN=32) :: strVar
INTEGER I, NBV1, NBV2, NELEM,NPOIN,NDP, TEMP, ERROR
INTEGER IPARAM(10)
INTEGER, allocatable :: IKLE(:,:)
INTEGER, allocatable :: IPOBO(:)
REAL*4, allocatable :: X(:)
REAL*4, allocatable :: Y(:)
REAL*4 TIME
REAL*4 , allocatable :: RES(:)
OPEN(1, FILE="C:\Telemac\bin\Yu_3d_sm_240940_5.slf", FORM="UNFORMATTED",STATUS='OLD',IOSTAT=IERR)
IF(IERR.NE.0) PRINT*,'ERROR ',IERR
OPEN(2, FILE="C:\Telemac\bin\Sta2.txt", FORM="FORMATTED")
READ(1, END=5, ERR=7) strTitle(1:80)
WRITE(2,*) strTitle(1:80)
!1 record containing the two integers NBV(1) and NBV(2) (NBV(1) the number of variables, NBV(2) with the value of 0),
READ(1, END=5, ERR=7) NBV1, NBV2
WRITE(2,*) NBV1, NBV2
! NBV(1) records containing the names and units of each variable (over 32 characters),
DO 20 j=1, NBV1
READ(1, END=5, ERR=7) strVar(1:32)
WRITE(2,*) strVar(1:32)
20 continue
! 1 record containing the integers table IPARAM (10 integers, of which only 4 are currently being used).
READ(1, END=5, ERR=7) IPARAM(1:10)
WRITE(2,*) IPARAM(1:10)
READ(1, END=5, ERR=7) NELEM,NPOIN,NDP, TEMP
WRITE(2,*) NELEM,NPOIN,NDP, TEMP
WRITE(2,*) "=================IKLE================"
!1 record containing table IKLE (integer array of dimension (NDP,NELEM) which is the connectivity table
allocate(IKLE(NDP,NELEM),stat=ERROR)
if(ERROR.ne.0) THEN
PRINT*, "Could not allocate memory IKLE"
GO TO 1000
end if
READ(1, END=5, ERR=7) IKLE(1:NDP,1:NELEM)
WRITE(2,*) IKLE(1:NDP,1:NELEM)
WRITE(2,*) "=================IPOBO================"
!1 record containing table IPOBO (integer array of dimension NPOIN);
allocate(IPOBO(NPOIN),stat=ERROR)
if(ERROR.ne.0) THEN
PRINT*, "Could not allocate memory IPOBO"
GO TO 1000
end if
READ(1, END=5, ERR=7) IPOBO(1:NPOIN)
WRITE(2,*) IPOBO(1:NPOIN)
WRITE(2,*) "=================X================"
!1 record containing table X (real array of dimension NPOIN containing the abscissas of the points),
allocate(X(NPOIN),stat=ERROR)
if(ERROR.ne.0) THEN
PRINT*, "Could not allocate memory IPOBO"
GO TO 1000
end if
READ(1, END=5, ERR=7) X(1:NPOIN)
WRITE(2,*) X(1:NPOIN)
WRITE(2,*) "=================Y================"
!1 record containing table Y (real array of dimension NPOIN containing the ordinates of the points),
allocate(Y(NPOIN),stat=ERROR)
if(ERROR.ne.0) THEN
PRINT*, "Could not allocate memory IPOBO"
GO TO 1000
end if
READ(1, END=5, ERR=7) Y(1:NPOIN)
WRITE(2,*) Y(1:NPOIN)
! FOR EACH TIMESTEP
allocate(RES(NPOIN),stat=ERROR)
if(ERROR.ne.0) THEN
PRINT*, "Could not allocate memory RES"
GO TO 1000
end if
do
READ(1, END=5, ERR=7) TIME
WRITE(2,*) TIME
! For each record
DO 30 j=1, NBV1
READ(1, END=5, ERR=7) RES(1:NPOIN)
WRITE(2,*) RES(1:NPOIN)
30 continue
end do
5 PRINT*,'END OF fILE 1'
GO TO 1000
7 PRINT*,'ERROR IN FILE 1'
GO TO 1000
999 WRITE(2,*) strTitle(1:80)
1000 CONTINUE
CLOSE(1)
CLOSE(2)
STOP
END
Please excuse any errors with this file - I have run it successfully but have not had the chance to verify the results. This is because I ended up using MATLAB. I found a nice set of scripts for reading and writing Selafin files on the MATLAB website. It has a set of m-files used to read and write the Selafin file. The main ones that I was using were:
MATLAB
I needed to generate the mean flow velocities over a specific timespan so I modified telmean.m to take in a start and end timestep. Next I developed my own script to parse the m.RESULT variable and plot the profiles.
The variables for each time step are stored in an (NPOIN * NBV x 1) array- where NPOIN is the number of points and NBV is the number of variables (e.g. elevation, U, V, W for my model runs). This is a very long array and parsing it into the variables is a somewhat complex task.
In the MATAB scripts, the array is accessible through the variable m.RESULT. The RESULT file for a 3D Selafin file is written out by the vertical layers. Each layer has NPOIN/NPLAN points, where NPLAN is the number of vertical layers. The parser that I wrote ends up creating and populating X,Y,Z,U,V,W arrays which are much more easy to access. It then generates a velocity profile along one of my "channels". There is the code:
function m = telplotprofile(INFILE,CHANNEL,NCCHAN,vectorScale,INTERVAL)
%****M* Telemac/telplotprofile.m
% NAME
% telplotprofile.m
%
% PURPOSE
% Reads a processed Telemac file (e.g. mean of set of time steps) and plots
% the velocity profile. Assumes that the channel mesher has been used.
%
% USAGE
% m = telplotprofile(INFILE,CHANNEL,vectorScale,INTERVAL)
%
% INPUTS
% INFILE - the input file
% CHANNEL - the channel number
% vectorScale - scaling factor for the
% INTERVAL - spacing of the velocity profiles (1 = all, 2 = every
% second)
% NCCHAN - number of cross-channel elements
%
% Example Usage:
% inputFile = 'C:\Working\Civil850\Analysis\Yu_3d_sm_240940_5';
% meanFile = 'C:\Working\Civil850\Analysis\Yu_3d_sm_240940_5_mean.slf';
% telmean2(inputFile,41, 100, meanFile);
% telplotprofile(meanFile,2,9,0.02,3);
%
% author: Patrick Grover
% email: patrick.grover@queensu.ca
% release date: 31-Aug-2011
% For testing run telplotprofile('') from the command line
if isempty(INFILE)
INFILE ='C:\Working\Civil850\Analysis\Yu_3d_sm_240940_5_mean.slf';
CHANNEL = 9; % The channel to plot INPUT
vectorScale = 0.02;
INTERVAL = 1;
NCCHAN = 9; % number of cross-channel elements
end
% Read the output file which now has the averaged values
m = telheadr(INFILE);
m.NPOINCHAN = m.NPOIN / (m.NPLAN*NCCHAN); % the number of points along a profile (downstream
m.NPOINLAY = m.NPOIN / (m.NPLAN); % the number of points on a vertical layer
m.X = zeros(m.NPOINCHAN,1);
m.Y = zeros(m.NPOINCHAN,1);
m.Z = zeros(m.NPOINCHAN,m.NPLAN);
m.U = zeros(m.NPOINCHAN,m.NPLAN);
m.V = zeros(m.NPOINCHAN,m.NPLAN);
m.W = zeros(m.NPOINCHAN,m.NPLAN);
m = telstepr(m,1); %since is an average file - read only the first timestep
% Read in the X and Y values
for i=1 : m.NPOINCHAN
POINT = ((CHANNEL-1)*m.NPOINCHAN) + i;
m.X(i) =m.XYZ(POINT,1);
m.Y(i) =m.XYZ(POINT,2);
%fprintf('\nX = %d\n',m.X(i));
%fprintf('\nY = %d\n',m.Y(i));
end
% now read in the z values
for layer=1:m.NPLAN
for i=1 : m.NPOINCHAN
POINT = ((CHANNEL-1)*m.NPOINCHAN) + i;
%fprintf('\nPOINT = %d\n',POINT);
POINT = POINT + ((layer-1)*m.NPOINLAY);
%fprintf('\nPOINT = %d\n',POINT);
m.Z(i,layer) =m.RESULT(POINT);
%fprintf('\nZ = %d\n',m.Z(i));
end
end
% now read in the U,V,W values
for layer=1:m.NPLAN
for i=1 : m.NPOINCHAN
POINT = ((CHANNEL-1)*m.NPOINCHAN) + i;
%fprintf('\nPOINT = %d\n',POINT);
POINT = POINT + ((layer-1)*m.NPOINLAY);
POINTU = POINT + m.NPOIN;
POINTV = POINT + m.NPOIN*2;
POINTW = POINT + m.NPOIN*3;
m.U(i,layer) = m.RESULT(POINTU);
m.V(i,layer) = m.RESULT(POINTV);
m.W(i,layer) = m.RESULT(POINTW);
end
end
umax = max(max(m.U()));
m.U = m.U ./ umax;
figure; hold on;
for i=1 : INTERVAL: m.NPOINCHAN
bFlag = 1; % a simple flag that indicates if it is the first time through
for layer=1:m.NPLAN
X0 = m.X(i);
Z0 = m.Z(i,layer);
X1 = X0 + (m.U(i,layer))*vectorScale;
Z1 = Z0;
VX = [X0 X1];
VY = [Z0 Z1];
if m.U(i,layer) > 0
plot(VX,VY,'k');
else
plot(VX,VY,'r');
end
if bFlag~=1
BackX = [X0OLD0 X0];
BackZ = [Z0OLD0 Z0];
plot(BackX,BackZ,'k');
FrontX = [X0OLD1 X1];
FrontZ = [Z0OLD1 Z1];
plot(FrontX,FrontZ,'k');
end
bFlag = 0;
X0OLD0 = X0;
X0OLD1 = X1;
Z0OLD0 = Z0;
Z0OLD1 = Z1;
end
end
% Finally plot the surface.
bFlag = 1; % a simple flag that indicates if it is the first time through
for i=1: m.NPOINCHAN
X1 = m.X(i);
Z1 = m.Z(i,1);
if bFlag~=1;
BottomX = [X0 X1];
BottomZ = [Z0 Z1];
fprintf('\nZ = %d\n',BottomZ);
plot(BottomX,BottomZ,'k');
end
X0 =X1;
Z0 = Z1;
bFlag = 0;
end
hold off;
This generates the following velocity profile graph:
Fortan
The first thing that you will need to do is extract your model results from the Selafin file. The Selafin file is in Binary format which makes things a little complex. The structure of a Selafin file is well described in Appendix 3 of the TELEMAC-2D User Manual. I have not programmed in Fortran since 1990, but I thought I would give it a try.
PROGRAM ReadSelaphin
! http://www.opentelemac.org/index.php?option=com_kunena&func=view&catid=20&id=1162&Itemid=62&lang=en
INTEGER IERR
CHARACTER(LEN=80) :: strTitle
CHARACTER(LEN=32) :: strVar
INTEGER I, NBV1, NBV2, NELEM,NPOIN,NDP, TEMP, ERROR
INTEGER IPARAM(10)
INTEGER, allocatable :: IKLE(:,:)
INTEGER, allocatable :: IPOBO(:)
REAL*4, allocatable :: X(:)
REAL*4, allocatable :: Y(:)
REAL*4 TIME
REAL*4 , allocatable :: RES(:)
OPEN(1, FILE="C:\Telemac\bin\Yu_3d_sm_240940_5.slf", FORM="UNFORMATTED",STATUS='OLD',IOSTAT=IERR)
IF(IERR.NE.0) PRINT*,'ERROR ',IERR
OPEN(2, FILE="C:\Telemac\bin\Sta2.txt", FORM="FORMATTED")
READ(1, END=5, ERR=7) strTitle(1:80)
WRITE(2,*) strTitle(1:80)
!1 record containing the two integers NBV(1) and NBV(2) (NBV(1) the number of variables, NBV(2) with the value of 0),
READ(1, END=5, ERR=7) NBV1, NBV2
WRITE(2,*) NBV1, NBV2
! NBV(1) records containing the names and units of each variable (over 32 characters),
DO 20 j=1, NBV1
READ(1, END=5, ERR=7) strVar(1:32)
WRITE(2,*) strVar(1:32)
20 continue
! 1 record containing the integers table IPARAM (10 integers, of which only 4 are currently being used).
READ(1, END=5, ERR=7) IPARAM(1:10)
WRITE(2,*) IPARAM(1:10)
READ(1, END=5, ERR=7) NELEM,NPOIN,NDP, TEMP
WRITE(2,*) NELEM,NPOIN,NDP, TEMP
WRITE(2,*) "=================IKLE================"
!1 record containing table IKLE (integer array of dimension (NDP,NELEM) which is the connectivity table
allocate(IKLE(NDP,NELEM),stat=ERROR)
if(ERROR.ne.0) THEN
PRINT*, "Could not allocate memory IKLE"
GO TO 1000
end if
READ(1, END=5, ERR=7) IKLE(1:NDP,1:NELEM)
WRITE(2,*) IKLE(1:NDP,1:NELEM)
WRITE(2,*) "=================IPOBO================"
!1 record containing table IPOBO (integer array of dimension NPOIN);
allocate(IPOBO(NPOIN),stat=ERROR)
if(ERROR.ne.0) THEN
PRINT*, "Could not allocate memory IPOBO"
GO TO 1000
end if
READ(1, END=5, ERR=7) IPOBO(1:NPOIN)
WRITE(2,*) IPOBO(1:NPOIN)
WRITE(2,*) "=================X================"
!1 record containing table X (real array of dimension NPOIN containing the abscissas of the points),
allocate(X(NPOIN),stat=ERROR)
if(ERROR.ne.0) THEN
PRINT*, "Could not allocate memory IPOBO"
GO TO 1000
end if
READ(1, END=5, ERR=7) X(1:NPOIN)
WRITE(2,*) X(1:NPOIN)
WRITE(2,*) "=================Y================"
!1 record containing table Y (real array of dimension NPOIN containing the ordinates of the points),
allocate(Y(NPOIN),stat=ERROR)
if(ERROR.ne.0) THEN
PRINT*, "Could not allocate memory IPOBO"
GO TO 1000
end if
READ(1, END=5, ERR=7) Y(1:NPOIN)
WRITE(2,*) Y(1:NPOIN)
! FOR EACH TIMESTEP
allocate(RES(NPOIN),stat=ERROR)
if(ERROR.ne.0) THEN
PRINT*, "Could not allocate memory RES"
GO TO 1000
end if
do
READ(1, END=5, ERR=7) TIME
WRITE(2,*) TIME
! For each record
DO 30 j=1, NBV1
READ(1, END=5, ERR=7) RES(1:NPOIN)
WRITE(2,*) RES(1:NPOIN)
30 continue
end do
5 PRINT*,'END OF fILE 1'
GO TO 1000
7 PRINT*,'ERROR IN FILE 1'
GO TO 1000
999 WRITE(2,*) strTitle(1:80)
1000 CONTINUE
CLOSE(1)
CLOSE(2)
STOP
END
Please excuse any errors with this file - I have run it successfully but have not had the chance to verify the results. This is because I ended up using MATLAB. I found a nice set of scripts for reading and writing Selafin files on the MATLAB website. It has a set of m-files used to read and write the Selafin file. The main ones that I was using were:
- telheadr.m - reads the header info from a Selafin file.
- telstepr.m - reads a specified time step
- telmean.m - calculates the mean from all time-steps.
MATLAB
I needed to generate the mean flow velocities over a specific timespan so I modified telmean.m to take in a start and end timestep. Next I developed my own script to parse the m.RESULT variable and plot the profiles.
The variables for each time step are stored in an (NPOIN * NBV x 1) array- where NPOIN is the number of points and NBV is the number of variables (e.g. elevation, U, V, W for my model runs). This is a very long array and parsing it into the variables is a somewhat complex task.
In the MATAB scripts, the array is accessible through the variable m.RESULT. The RESULT file for a 3D Selafin file is written out by the vertical layers. Each layer has NPOIN/NPLAN points, where NPLAN is the number of vertical layers. The parser that I wrote ends up creating and populating X,Y,Z,U,V,W arrays which are much more easy to access. It then generates a velocity profile along one of my "channels". There is the code:
function m = telplotprofile(INFILE,CHANNEL,NCCHAN,vectorScale,INTERVAL)
%****M* Telemac/telplotprofile.m
% NAME
% telplotprofile.m
%
% PURPOSE
% Reads a processed Telemac file (e.g. mean of set of time steps) and plots
% the velocity profile. Assumes that the channel mesher has been used.
%
% USAGE
% m = telplotprofile(INFILE,CHANNEL,vectorScale,INTERVAL)
%
% INPUTS
% INFILE - the input file
% CHANNEL - the channel number
% vectorScale - scaling factor for the
% INTERVAL - spacing of the velocity profiles (1 = all, 2 = every
% second)
% NCCHAN - number of cross-channel elements
%
% Example Usage:
% inputFile = 'C:\Working\Civil850\Analysis\Yu_3d_sm_240940_5';
% meanFile = 'C:\Working\Civil850\Analysis\Yu_3d_sm_240940_5_mean.slf';
% telmean2(inputFile,41, 100, meanFile);
% telplotprofile(meanFile,2,9,0.02,3);
%
% author: Patrick Grover
% email: patrick.grover@queensu.ca
% release date: 31-Aug-2011
% For testing run telplotprofile('') from the command line
if isempty(INFILE)
INFILE ='C:\Working\Civil850\Analysis\Yu_3d_sm_240940_5_mean.slf';
CHANNEL = 9; % The channel to plot INPUT
vectorScale = 0.02;
INTERVAL = 1;
NCCHAN = 9; % number of cross-channel elements
end
% Read the output file which now has the averaged values
m = telheadr(INFILE);
m.NPOINCHAN = m.NPOIN / (m.NPLAN*NCCHAN); % the number of points along a profile (downstream
m.NPOINLAY = m.NPOIN / (m.NPLAN); % the number of points on a vertical layer
m.X = zeros(m.NPOINCHAN,1);
m.Y = zeros(m.NPOINCHAN,1);
m.Z = zeros(m.NPOINCHAN,m.NPLAN);
m.U = zeros(m.NPOINCHAN,m.NPLAN);
m.V = zeros(m.NPOINCHAN,m.NPLAN);
m.W = zeros(m.NPOINCHAN,m.NPLAN);
m = telstepr(m,1); %since is an average file - read only the first timestep
% Read in the X and Y values
for i=1 : m.NPOINCHAN
POINT = ((CHANNEL-1)*m.NPOINCHAN) + i;
m.X(i) =m.XYZ(POINT,1);
m.Y(i) =m.XYZ(POINT,2);
%fprintf('\nX = %d\n',m.X(i));
%fprintf('\nY = %d\n',m.Y(i));
end
% now read in the z values
for layer=1:m.NPLAN
for i=1 : m.NPOINCHAN
POINT = ((CHANNEL-1)*m.NPOINCHAN) + i;
%fprintf('\nPOINT = %d\n',POINT);
POINT = POINT + ((layer-1)*m.NPOINLAY);
%fprintf('\nPOINT = %d\n',POINT);
m.Z(i,layer) =m.RESULT(POINT);
%fprintf('\nZ = %d\n',m.Z(i));
end
end
% now read in the U,V,W values
for layer=1:m.NPLAN
for i=1 : m.NPOINCHAN
POINT = ((CHANNEL-1)*m.NPOINCHAN) + i;
%fprintf('\nPOINT = %d\n',POINT);
POINT = POINT + ((layer-1)*m.NPOINLAY);
POINTU = POINT + m.NPOIN;
POINTV = POINT + m.NPOIN*2;
POINTW = POINT + m.NPOIN*3;
m.U(i,layer) = m.RESULT(POINTU);
m.V(i,layer) = m.RESULT(POINTV);
m.W(i,layer) = m.RESULT(POINTW);
end
end
umax = max(max(m.U()));
m.U = m.U ./ umax;
figure; hold on;
for i=1 : INTERVAL: m.NPOINCHAN
bFlag = 1; % a simple flag that indicates if it is the first time through
for layer=1:m.NPLAN
X0 = m.X(i);
Z0 = m.Z(i,layer);
X1 = X0 + (m.U(i,layer))*vectorScale;
Z1 = Z0;
VX = [X0 X1];
VY = [Z0 Z1];
if m.U(i,layer) > 0
plot(VX,VY,'k');
else
plot(VX,VY,'r');
end
if bFlag~=1
BackX = [X0OLD0 X0];
BackZ = [Z0OLD0 Z0];
plot(BackX,BackZ,'k');
FrontX = [X0OLD1 X1];
FrontZ = [Z0OLD1 Z1];
plot(FrontX,FrontZ,'k');
end
bFlag = 0;
X0OLD0 = X0;
X0OLD1 = X1;
Z0OLD0 = Z0;
Z0OLD1 = Z1;
end
end
% Finally plot the surface.
bFlag = 1; % a simple flag that indicates if it is the first time through
for i=1: m.NPOINCHAN
X1 = m.X(i);
Z1 = m.Z(i,1);
if bFlag~=1;
BottomX = [X0 X1];
BottomZ = [Z0 Z1];
fprintf('\nZ = %d\n',BottomZ);
plot(BottomX,BottomZ,'k');
end
X0 =X1;
Z0 = Z1;
bFlag = 0;
end
hold off;
This generates the following velocity profile graph:
Subscribe to:
Posts (Atom)