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.
This blog is my personal notes while conducting research for my PHD. My research interests are sediment transport and morphological evolution that occurs during dam break flows.
Thursday, September 20, 2012
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)






