Snip calculating mass flow

1 Description

This describes how to set up the solver in order to calculate the mass fluxes on the boundary. It works in serial and parallel.

It was taken from the thread Calculation of integral boundary values.

2 Preparing the solver

Firstly, you need to #include "SortableList.H" as one of the top level includes in a solver application.

Then #include "buildGlobalBoundaryList.H" once "createMesh.H" and other initialisations have been run in the main function. (This only needs to be run once.)

The result of this is a wordList variable called "globalPatchNames" that is available on each processor, and contains a sorted list of all real patches in the model.

To compute mass fluxes after each timestep, #include "computeMassFlux.H"

3 The files

3.1 buildGlobalBoundaryList.H

 
#include "SortableList.H"
 
wordList buildGlobalBoundaryList(polyMesh &mesh)
{
    int i;
    label numPatches;
    SortableList<word> sortablePatchNames(0);
    List<wordList> gatheredPatchNames(Pstream::nProcs());
 
    Info<< "Building global boundary list" << endl;
 
    // collect all real patches into a wordList
    wordList realPatches(0);
    forAll (mesh.boundaryMesh(), patchIndex)
    {
        const polyPatch& patch = mesh.boundaryMesh()[patchIndex];
        if (typeid(patch) != typeid(processorPolyPatch))
        {
            realPatches.setSize(realPatches.size()+1);
            realPatches[realPatches.size()-1] = patch.name();
        }
    }
 
    // gather real patches from each processor
    gatheredPatchNames[Pstream::myProcNo()] = realPatches;
    Pstream::gatherList(gatheredPatchNames);
 
    if (Pstream::master())
    {
        // store unique patches (using a HashSet)
        HashSet<word> hashSetPatchNames(0);
        forAll (gatheredPatchNames, patchListProc)
        {
            forAll (gatheredPatchNames[patchListProc], patchNameItr)
            {
                hashSetPatchNames.insert(
                    gatheredPatchNames[patchListProc][patchNameItr]
                    );
            }
        }
 
        // store patch names into a SortableList
        sortablePatchNames.setSize(hashSetPatchNames.size());
        i=0;
        for
        (
            HashSet<word>::const_iterator iter = hashSetPatchNames.begin();
            iter != hashSetPatchNames.end();
            ++iter
        )
        {
            sortablePatchNames[i++] = iter.key();
        }
 
        // *** SORT CAUSES SEG-FAULT **********
	// sortablePatchNames.sort();
        // ************************************
 
        // save the number of patches
        numPatches = sortablePatchNames.size();
    }
 
    // scatter the number of patches to each processor
    Pstream::scatter(numPatches);
 
    // reset the list size on all other processors
    if (!Pstream::master())
    {
        sortablePatchNames.setSize(numPatches);
    }
 
    // scatter individual patch names to each processor
    // and rebuild the patch list
    word tmpPatchName;
    for (i=0; i<numPatches; i++)
    {
        if(Pstream::master())
        {
            tmpPatchName = sortablePatchNames[i];
        }
        Pstream::scatter(tmpPatchName);
        if(!Pstream::master())
        {
            sortablePatchNames[i] = tmpPatchName;
        }
    }
 
    // return the global patch names
    return sortablePatchNames;
}

3.2 calculateMassFlux.H

 
{
    scalar netFlux = 0;
    forAll (globalBoundaryList, patchI)
    {
        scalar flux = 0;
        label patchIndex = mesh.boundaryMesh().findPatchID(globalBoundaryList[patchI]);
        if (patchIndex >= 0)
        {
            flux = sum(phi.boundaryField()[patchIndex]);
        }
        reduce(flux, sumOp<scalar>());
        Info<< "Mass flux at " << globalBoundaryList[patchI]
            << " = " << flux << endl;
        netFlux += flux;
    }
    Info<< "Net mass flux = " << netFlux << endl << endl;
}