OFtutorial5.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"

    // For a case being run in parallel, the domain is decomposed into several
    // processor meshes. Each of them is run in a separate process and holds
    // instances of objects like mesh, U or p just as in a single-threaded (serial)
    // computation. These will have different sizes, of course, as they hold
    // fewer elements than the whole, undecomposed, mesh.
    // Pout is a stream to which each processor can write, unlike Info which only
    // gets used by the head process (processor0)
    Pout << "Hello from processor " << Pstream::myProcNo() << "! I am working on "
         << mesh.C().size() << " cells" << endl;

    // To exchange information between processes, special OpenMPI routines need
    // to be called.

    // This goes over each cell in the subdomain and integrates their volume.
    scalar meshVolume(0.);
    forAll(mesh.V(),cellI)
        meshVolume += mesh.V()[cellI];

    // Add the values from all processes together
    Pout << "Mesh volume on this processor: " << meshVolume << endl;
    reduce(meshVolume, sumOp<scalar>());
    Info << "Total mesh volume on all processors: " << meshVolume
        // Note how the reudction operation may be done in place without defning
        // a temporary variable, where appropriate.
         << " over " << returnReduce(mesh.C().size(), sumOp<label>()) << " cells" << endl;
    // During the reduction stage, different operations may be carried out, summation,
    // described by the sumOp template, being one of them.
    // Other very useful operations are minOp and maxOp.
    // Note how the type
    // of the variable must be added to make an instance of the template, here
    // this is done by adding <scalar> in front of the brackets.
    // Custom reduction operations are easy to implement but need fluency in
    // object-oriented programming in OpenFOAM, so we'll skip this for now.

    // Spreading a value across all processors is done using a scatter operation.
    Pstream::scatter(meshVolume);
    Pout << "Mesh volume on this processor is now " << meshVolume << endl;

    // It is often useful to check the distribution of something across all
    // processors. This may be done using a list, with each element of it
    // being written to by only one processor.
    List<label> nInternalFaces (Pstream::nProcs()), nBoundaries (Pstream::nProcs());
    nInternalFaces[Pstream::myProcNo()] = mesh.Cf().size();
    nBoundaries[Pstream::myProcNo()] = mesh.boundary().size();

    // The list may then be gathered on the head node as
    Pstream::gatherList(nInternalFaces);
    Pstream::gatherList(nBoundaries);
    // Scattering a list is also possbile
    Pstream::scatterList(nInternalFaces);
    Pstream::scatterList(nBoundaries);

    // It can also be useful to do things on the head node only
    // (in this case this is meaningless since we are using Info, which already
    // checks this and executes on the head node).
    // Note how the gathered lists hold information for all processors now.
    if (Pstream::master())
    {
        forAll(nInternalFaces,i)
            Info << "Processor " << i << " has " << nInternalFaces[i]
                 << " internal faces and " << nBoundaries[i] << " boundary patches" << endl;
    }

    // As the mesh is decomposed, interfaces between processors are turned
    // into patches, meaning each subdomain sees a processor boundary as a
    // boundary condition.
    forAll(mesh.boundary(),patchI)
        Pout << "Patch " << patchI << " named " << mesh.boundary()[patchI].name() << endl;

    // When looking for processor patches, it is useful to check their type,
    // similarly to how one can check if a patch is of empty type
    forAll(mesh.boundary(),patchI)
    {
        const polyPatch& pp = mesh.boundaryMesh()[patchI];
        if (isA<processorPolyPatch>(pp))
            Pout << "Patch " << patchI << " named " << mesh.boundary()[patchI].name()
                 << " is definitely a processor boundary!" << endl;
    }

    // ---
    // this is an example implementation of the code from tutoral 2 which
    // has been adjusted to run in parallel. Each difference is highlighted
    // as a NOTE.

    // It is conventional in OpenFOAM to move large parts of code to separate
    // .H files to make the code of the solver itself more readable. This is not
    // a standard C++ practice, as header files are normally associated with
    // declarations rather than definitions.
    // A very common include, apart from the setRootCase, createTime, and createMesh,
    // which are generic, is createFields, which is often unique for each solver.
    // Here we've moved all of the parts of the code dealing with setting up the fields
    // and transport constants into this include file.
    #include "createFields.H"

    // pre-calculate geometric information using field expressions rather than
    // cell-by-cell assignment.
	const dimensionedVector originVector("x0", dimLength, vector(0.05,0.05,0.005));
    volScalarField r (mag(mesh.C()-originVector));
    // NOTE: we need to get a global value; convert from dimensionedScalar to scalar
	const scalar rFarCell = returnReduce(max(r).value(), maxOp<scalar>());
    scalar f (1.);

	Info<< "\nStarting time loop\n" << endl;

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        // assign values to the field;
        // sin function expects a dimensionless argument, hence need to convert
        // current time using .value().
        // r has dimensions of length, hence the small value being added to it
        // needs to match that.
        // Finally, the result has to match dimensions of pressure, which are
        // m^2 / s^-2/
		p = Foam::sin(2.*constant::mathematical::pi*f*runTime.time().value())
            / (r/rFarCell + dimensionedScalar("small", dimLength, 1e-12))
            * dimensionedScalar("tmp", dimensionSet(0, 3, -2, 0, 0), 1.);

        // NOTE: this is needed to update the values on the processor boundaries.
        // If this is not done, the gradient operator will get confused around the
        // processor patches.
        p.correctBoundaryConditions();

        // calculate velocity from gradient of pressure
		U = fvc::grad(p)*dimensionedScalar("tmp", dimTime, 1.);
		runTime.write();
	}

    Info<< "End\n" << endl;

    return 0;
}

// ************************************************************************* //