Usage examples

This section documents the example scripts shipped with the source code.

Lid-driven cavity flow

The openfoam_cavity_tutorial.py script provides an example of low-level manipulation of OpenFOAM cases. In this example we shall re-create the initial example from the OpenFOAM users’ guide. It’s worth reading over that section first before trying to follow the transliteration below.

Firstly, we need to import some things from the firefish.case module:

from firefish.case import Case, FileName, FileClass

The Case class encapsulates an OpenFOAM case directory. We don’t want to overwrite an existing case and so we write a little convenience wrapper function:

def create_new_case(case_dir):
    # Check that the specified case directory does not already exist
    if os.path.exists(case_dir):
        raise RuntimeError(
            'Refusing to write to existing path: {}'.format(case_dir)
        )

    # Create the case
    return Case(case_dir)

The mutable_data_file() method will return a context manager which can be used to manipulate an OpenFOAM data file. The data file is created if it does not yet exist, its contents are parsed into a dictionary and the dictionary is returned from the context manager. One the context is left the dictionary is re-written to disk.

The upshot of this is that the programmer is insulated from manipulating OpenFOAM data files directly. Let’s write the controlDict file from the tutorial:

def write_initial_control_dict(case):
    # Control dict from tutorial
    control_dict = {
        'application': 'icoFoam',
        'startFrom': 'startTime',
        'startTime': 0,
        'stopAt': 'endTime',
        'endTime': 0.5,
        'deltaT': 0.005,
        'writeControl': 'timeStep',
        'writeInterval': 20,
        'purgeWrite': 0,
        'writeFormat': 'ascii',
        'writePrecision': 6,
        'writeCompression': 'off',
        'timeFormat': 'general',
        'timePrecision': 6,
        'runTimeModifiable': True,
    }

    with case.mutable_data_file(FileName.CONTROL) as d:
        d.update(control_dict)

Well-known file names are available through the FileName class.

The blockMeshDict file is the next one to be created. This is an example of a relatively complex file. The complexity is somewhat hidden by the mapping to and from the Python domain but there is still some subtlety. Notice particularly the rather odd way in which the boundary dictionary is specified:

def write_block_mesh_dict(case):
    block_mesh_dict = {
        'convertToMeters': 0.1,

        'vertices': [
            [0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0],
            [0, 0, 0.1], [1, 0, 0.1], [1, 1, 0.1], [0, 1, 0.1],
        ],

        'blocks': [
            (
                'hex', [0, 1, 2, 3, 4, 5, 6, 7], [20, 20, 1],
                'simpleGrading', [1, 1, 1],
            )
        ],

        'edges': [],

        # Note the odd way in which boundary is defined here as a
        # list of tuples.
        'boundary': [
            ('movingWall', {
                'type': 'wall',
                'faces': [ [3, 7, 6, 2] ],
            }),
            ('fixedWalls', {
                'type': 'wall',
                'faces': [
                    [0, 4, 7, 3],
                    [2, 6, 5, 1],
                    [1, 5, 4, 0],
                ],
            }),
            ('frontAndBack', {
                'type': 'empty',
                'faces': [
                    [0, 3, 2, 1],
                    [4, 5, 6, 7],
                ],
            }),
        ],

        'mergePatchPairs': [],
    }

    with case.mutable_data_file(FileName.BLOCK_MESH) as d:
        d.update(block_mesh_dict)

At this point in the tutorial we’re ready to run the blockMesh command which is one function call away:

case.run_tool('blockMesh')

We’re close to being able to run the icoFoam utility. The transport properties need to be defined:

from firefish.case import Dimension

with case.mutable_data_file(FileName.TRANSPORT_PROPERTIES) as tp:
   tp['nu'] = (Dimension(0, 2, -1, 0, 0, 0, 0), 0.01)

We also need to create the initial conditions. Notice that we have to specify a different class when creating them:

def write_initial_conditions(case):
    # Create the p initial conditions
    p_file = case.mutable_data_file(
        '0/p', create_class=FileClass.SCALAR_FIELD_3D
    )
    with p_file as p:
        p.update({
            'dimensions': Dimension(0, 2, -2, 0, 0, 0, 0),
            'internalField': ('uniform', 0),
            'boundaryField': {
                'movingWall': { 'type': 'zeroGradient' },
                'fixedWalls': { 'type': 'zeroGradient' },
                'frontAndBack': { 'type': 'empty' },
            },
        })

    # Create the U initial conditions
    U_file = case.mutable_data_file(
        '0/U', create_class=FileClass.VECTOR_FIELD_3D
    )
    with U_file as U:
        U.update({
            'dimensions': Dimension(0, 1, -1, 0, 0, 0, 0),
            'internalField': ('uniform', [0, 0, 0]),
            'boundaryField': {
                'movingWall': {
                    'type': 'fixedValue', 'value': ('uniform', [1, 0, 0])
                },
                'fixedWalls': {
                    'type': 'fixedValue', 'value': ('uniform', [0, 0, 0])
                },
                'frontAndBack': { 'type': 'empty' },
            },
        })

Before we can run icoFoam, we must create the mysterious fvSolution file:

def write_fv_solution(case):
    fv_solution = {
        'solvers': {
            'p': {
                'solver': 'PCG',
                'preconditioner': 'DIC',
                'tolerance': 1e-6,
                'relTol': 0,
            },
            'U': {
                'solver': 'smoothSolver',
                'smoother': 'symGaussSeidel',
                'tolerance': 1e-5,
                'relTol': 0,
            },
        },
        'PISO': {
            'nCorrectors': 2,
            'nNonOrthogonalCorrectors': 0,
            'pRefCell': 0,
            'pRefValue': 0,
        }
    }

    with case.mutable_data_file(FileName.FV_SOLUTION) as d:
        d.update(fv_solution)

And also the equally mysterious fvSchemes file:

def write_fv_schemes(case):
    fv_schemes = {
        'ddtSchemes': { 'default': 'Euler' },
        'gradSchemes': { 'default': 'Gauss linear', 'grad(p)': 'Gauss linear' },
        'divSchemes': { 'div(phi,U)': 'Gauss linear', 'default': 'none' },
        'laplacianSchemes': { 'default': 'Gauss linear orthogonal' },
        'interpolationSchemes': { 'default': 'linear' },
        'snGradSchemes': { 'default': 'orthogonal' },
    }

    with case.mutable_data_file(FileName.FV_SCHEMES) as d:
        d.update(fv_schemes)

The example script includes a main() function which performs all of these steps:

def main(case_dir='cavity'):
    # Create a new case file, raising RuntimeError if the directory already
    # exists.
    case = create_new_case(case_dir)

    # Add the information needed by blockMesh.
    write_initial_control_dict(case)
    write_block_mesh_dict(case)

    # At this point there is enough to run blockMesh.
    case.run_tool('blockMesh')

    # Update the physical properties.
    with case.mutable_data_file(FileName.TRANSPORT_PROPERTIES) as tp:
        tp['nu'] = (Dimension(0, 2, -1, 0, 0, 0, 0), 0.01)

    # Write the fvSolution and fvSchemes files.
    write_fv_solution(case)
    write_fv_schemes(case)

    # Write the initial conditions for the p and U fields.
    write_initial_conditions(case)

    # Run the icoFoam application.
    case.run_tool('icoFoam')

After the example script is run, paraFoam may be run to inspect the solution.

Supersonic flow over wedge

The openfoam_supersonic_wedge.py script provides an example of setting up a compressible flow solver in OpenFoam.

As in openfoam_cavity_tutorial.py we set up the OpenFOAM case directory using the firefish.case framework.

For flows with a Mach number above 0.3 compressible effects become non negligible. A compressible solver must therefore be used. In this case we use rhoCentralFoam. The control file must therefore be set accordingly:

def write_control_dict(case, n_iter):
	"""Sets up the control dictionary.
	In this example we use the rhoCentralFoam compressible solver"""

	# Control dict from tutorial
	control_dict = {
		'application': 'rhoCentralFoam',
		'startFrom': 'startTime',
		'startTime': 0,
		'stopAt': 'endTime',
		'endTime': n_iter,
		'deltaT': 0.001,
		'writeControl': 'runTime',
		'writeInterval': 1,
		'purgeWrite': 0,
		'writeFormat': 'ascii',
		'writePrecision': 6,
		'writeCompression': 'off',
		'timeFormat': 'general',
		'timePrecision': 6,
		'runTimeModifiable': True,
		'adjustTimeStep' : 'no',
		'maxCo' : 1,
		'maxDeltaT' : 1e-6,
	}

	with case.mutable_data_file(FileName.CONTROL) as d:
		d.update(control_dict)

The mesh needs to be set up using the blockMeshDict. The mesh consits of three blocks in order to model the upstream and downstream portions as well as the wedge itself. The numbering order in which the vertices are set is very important! We declare a block via:

hex', [0, 7, 2, 1, 8, 15, 10, 9], [40, 40, 1],
            'simpleGrading', [1, 1, 1]

The first line declares which vertices make up the corners of the block. This explanation best describes the order in which the vertices should be listed. The second part of the statement describes the cell density within the block in each of the three directions. The last part is used when we want the mesh density to vary within the block.

We must also set the thermodynamic properties of the gas. In this case the properties have been chosen so that the gas has a ratio of specific heats of 1.4 and that, if the temperature is 1K, then the speed of sound is 1m/s. As this is a commonly used fluid it can be done using the write_standard_thermophysical_properties function in the firefish.case module. We do this via:

write_standard_thermophysical_properties(case, StandardFluid.DIMENSIONLESS_AIR)

As this is a compressible flow we must also set the initial value of the temperatue field. Notice also that for the velocity we have set a slip boundary condition on the solid walls. This is because we are using an inviscid solver. When we move to a viscous solver we must set a no slip boundary condition on the solid walls.

def write_initial_conditions(case):
	"""Sets the initial conditions"""
	# Create the p initial conditions
	p_file = case.mutable_data_file(
		'0/p', create_class=FileClass.SCALAR_FIELD_3D
	)
	with p_file as p:
		p.update({
			'dimensions': Dimension(1, -1, -2, 0, 0, 0, 0),
			'internalField': ('uniform', 1),
			'boundaryField': {
				'inlet' : {'type' : 'fixedValue', 'value' : 'uniform 1'},
				'outlet': {'type': 'zeroGradient'},
				'fixedWalls': {'type': 'zeroGradient'},
				'frontAndBack': {'type': 'empty'},
			},
		})

	# Create the U initial conditions
	U_file = case.mutable_data_file(
		'0/U', create_class=FileClass.VECTOR_FIELD_3D
	)
	with U_file as U:
		U.update({
			'dimensions': Dimension(0, 1, -1, 0, 0, 0, 0),
			'internalField': ('uniform', [2, 0, 0]),
			'boundaryField': {
				'inlet' : {'type' : 'fixedValue',
						   'value' : ('uniform', [2, 0, 0])},
				'outlet': {
					'type': 'zeroGradient'
				},
				'fixedWalls': {
					'type': 'slip'
				},
				'frontAndBack': {'type': 'empty'},
			},
		})
		# Create the T initial conditions
	T_file = case.mutable_data_file(
		'0/T', create_class=FileClass.SCALAR_FIELD_3D
	)
	with T_file as T:
		T.update({
			'dimensions': Dimension(0, 0, 0, 1, 0, 0, 0),
			'internalField': ('uniform', 1),
			'boundaryField': {
				'inlet' : {'type' : 'fixedValue', 'value' : ('uniform', 1)},
				'outlet': {
					'type': 'zeroGradient'
				},
				'fixedWalls': {
					'type': 'zeroGradient'
				},
				'frontAndBack': {'type': 'empty'},
			},
		})

The example script includes a main() function which performs all of these steps. A boolean value can be passed to main() in order to reduce the number of iterations and so speed up automatic testing.

def main(case_dir='wedge', n_iter=10):
    #Try to create new case directory
    case = create_new_case(case_dir)
    # Add the information needed by blockMesh.
    write_control_dict(case, n_iter)
    write_block_mesh_dict(case)
    #we generate the mes1h
    case.run_tool('blockMesh')
    #we prepare the thermophysical and turbulence properties
    write_standard_thermophysical_properties(case, StandardFluid.DIMENSIONLESS_AIR)
    write_turbulence_properties(case)
    #we write fvScheme and fvSolution
    write_fv_schemes(case)
    write_fv_solution(case)
    write_initial_conditions(case)
    case.run_tool('rhoCentralFoam')

After the example script is run, paraFoam may be run to inspect the solution.

Snappy Hex Mesh Example

The snappy_hex_example.py script provides an example of running snappyHexMesh.

In order to be able to run snappyHexMesh we need to set up a control dictionary even though it plays no part in the actual mesh generation process. Likewise, in order to use paraFoam, we need to set up fvSchemes and fvSolution.

For snappyHexMesh to work we must have an underlying block mesh. This is generated in make_block_mesh and follows the same procedure as in previous examples.

For the actual mesh generation we first of all load a geometry using the new Geometry class:

rocket = Geometry(GeometryFormat.STL,'example.stl','example',case)

The idea behind this class is that, when we support more geometry file types in the future, it will abstract away the need to worry about wether something is an STL or OBJ etc.

We next scale and translate the rocket:

rocket.scale(0.5);
rocket.translate([0.5,2,2])

The Geometry class also contains mesh quality settings for this particular geometry.

Now that the geometry has been loaded we use it to initialise the SnappyHexMesh class:

snap = SnappyHexMesh(rocket,4,case)

This creates a new SnappyHexMesh class based on the example geometry and with a surface refinement level of 4. The SnappyHexMesh class automatically sets the mesh generation settings to a set of default values. These can be altered:

snap.snap=True
snap.snapTolerance = 8;
snap.locationToKeep = [0.0012,0.124,0.19]
snap.addLayers=False

Once the mesh generator is set up we can make the mesh via the call:

snap.generate_mesh()

Several things happen all at once here: * Surface features are extracted from the geometry (saving the STL in the trisurfaces directory if it has not already done so) * The mesh quality settings are written to a dictionary * The snappyHexMesh dictionary is written using the attributes of the SnappyHexMesh class * snappyHexMesh is run as a tool within the case directory

Once this has run the mesh can be viewed via paraFoam

Join STL Example

The join_stl_example.py script provides an example of combining multiple STL files into a single geometry and then generating a mesh through snappyHexMesh. It is worth having a look at snappy_hex_example.py first in order to get a more detailed overview on how SnappyHexMesh works.

It is now very straightforwards to generate a mesh made up from multiple .STL files.

Firstly one needs to make a list containing the paths of each STL file:

path_list = ['STLS/nosecone.stl', 'STLS/upperTube.stl', 'STLS/lowerTube.stl', 'STLS/finCan.stl', 'STLS/tail.stl']

One then needs to make a list of the human-readable names that correspond to each file:

part_list = ['nosecone', 'upperTube', 'lowerTube', 'finCan', 'tail']

When examining the mesh in paraFoam or when getting force outputs it will be these names that appear.

Once this has been done the Geometry classes can be loaded by a single call of:

parts = load_multiple_geometries(GeometryFormat.STL,path_list,part_list,case)

This produces a list of firefish.geometry.Geometry objects which can scaled, translated and rotated independently as required using the normal geometry functions.

We now use this list of Geometry objects to initialise SnappyHexMesh:

snap = SnappyHexMesh(parts,4,case)

As in snappy_hex_example.py, we can now alter the settings of Snappy Hex Mesh by altering the attributes of the SnappyHexMesh class. Once we’ve updated these we generate the mesh via:

snap.generate_mesh()

This generates four different meshes: the blank block mesh, the castellated mesh, the snapped mesh and a mesh with layer addition. In order to use the final mesh as the starting point of our simulation we perform some trickery to delete the meshes we don’t want and move the mesh we do want into the constant folder

def getTrueMesh(case):
    #the proper mesh is in the final time directory, delete the one in constant
    os.chdir(case.root_dir_path)
    call (["rm", "-r", "constant/polyMesh"])
    call (["mv", "0.002/polyMesh", "constant/"])
    call (["rm", "-r", "0.001/"])
    call (["rm", "-r", "0.002/"])
    os.chdir("../")

Kinematics Example

The kinematics_example.py script uses the new firefish.kinematics framework to implement the kinematics example in kinematics_basis.py

We firstly define a class that inherits from firefish.kinematics.KinematicBody. We want to model a cylindrical rocket whose principal moments of inertia vary as the rocket burns its motor. In order to vary the prinicpal moments of inertia automatically during the timestepping routine we must overload the update_moi() function.

class CylinderRocket(KinematicBody):
    """A rocket that is a cylinder with  evenly distributed mass"""
    def update_moi(self):
        radius = 0.3
        height = 2
       
        Ixx = (self.mass/12.0)*(3*radius**2 + height**2)
        Iyy = (self.mass/12.0)*(3*radius**2 + height**2)
        Izz = self.mass*radius**2 / 2.0
        
        self.MoI = [Ixx,Iyy,Izz]

In the main function we now undergo the time stepping routine. For each time step we must pass the forces acting on the rocket along with the torques to the routine. These forces and torques must be given in the body coordinate system. The example here burns the motor for fifty seconds and then lets it coast

def main():
    """Run through the simulation with a 50s motor burn"""
    initialMass = 100
    initialInertias = [0, 0, 0]
    rocket = CylinderRocket(initialMass,initialInertias)
    dt =1; duration = 100; gravity = 9.81;
    simulation = KinematicSimulation(rocket,gravity,duration,dt)
    while simulation.tIndex*dt <= duration:

        thrust = 0

        if simulation.tIndex*dt <= 50:
            thrust = 2000.0

        forces = [2.0, 0.0, thrust]
        torques = [0.0, 0.0, 0.0]
        mdot = 0.1

        i = simulation.tIndex
       
        simulation.time_step(forces,torques,mdot)
    
    return simulation.posits[100,2] #z position