calcProduction
{
    name            calcProduction;
    type            coded;
    libs            ("libutilityFunctionObjects.so");
    enabled         true;
    log             true;
    timeStart       $startTimeFieldAverage;
    timeEnd         $endTime;
    executeControl  adjustableRunTime;
    executeInterval $writeInterval;
    writeControl    adjustableRunTime;
    writeInterval   $writeInterval;
    
    // Basic includes
    codeOptions
    #{ -I$(LIB_SRC)/OpenFOAM/lnInclude -I$(LIB_SRC)/finiteVolume/lnInclude #}; 

    // --- Input(s) ---
    UMeanName       UMean;
    RMeanName       UPrime2Mean;

    // --- Output(s) ---
    OutputName      P_ij_RSTE;

    // Code for calculating Production
    codeWrite
    #{
        #include "fvCFD.H"

        // Input-Output namelist
        const word uMeanNameIn = this->dict().lookupOrDefault<word>("UMeanName", "UMean");
        const word RMeanNameIn = this->dict().lookupOrDefault<word>("RMeanName", "UPrime2Mean");
        const word prodOutNm   = this->dict().lookupOrDefault<word>("OutputName", "P_ij_RSTE");

        Info<< "Calculating Production Term for FO: " << this->name() << " at t=" << mesh().time().timeName() << endl;

        // ----- Reading Input variables ----- //

        // Lookup UMean
        if (!mesh().foundObject<volVectorField>(uMeanNameIn)) { FatalErrorInFunction << "Field " << uMeanNameIn << " not found." << exit(FatalError); }
        const volVectorField& UMean = mesh().lookupObject<volVectorField>(uMeanNameIn);
        
        // Lookup UPrime2Mean
        if (!mesh().foundObject<volSymmTensorField>(RMeanNameIn)) { FatalErrorInFunction << "Field " << RMeanNameIn << " not found." << exit(FatalError); }
        const volSymmTensorField& RMean = mesh().lookupObject<volSymmTensorField>(RMeanNameIn);
        
        Info<< "  Inputs: " << UMean.name() << ", " << RMean.name() << endl;

        // ----- Preparing Output ----- //

        // Calculating velocity gradient tensor
        volTensorField gradUMean = fvc::grad(UMean);

        // Constructing production as a volSymmTensor 
        volSymmTensorField P_ij_RSTE
        (
            IOobject(prodOutNm, mesh().time().timeName(), mesh(), IOobject::NO_READ, IOobject::AUTO_WRITE),
            -2.0 * symm(RMean & gradUMean)
        );


        // Writing Output
        P_ij_RSTE.correctBoundaryConditions();
        P_ij_RSTE.write();

        Info<< "  Written field: " << P_ij_RSTE.name() << endl;
        Info<< endl;
    #};
}