// Fields for data gathering

List< scalarField > allSpeciesN_RU
(
    molecules.pairPotentials().nIds(),
    scalarField (mesh.nCells(), 0.0)
);

List< scalarField > allSpeciesM_RU
(
    molecules.pairPotentials().nIds(),
    scalarField (mesh.nCells(), 0.0)
);

List< vectorField > allSpeciesVelocitySum_RU
(
    molecules.pairPotentials().nIds(),
    vectorField (mesh.nCells(), vector::zero)
);

List< scalarField > allSpeciesVelocityMagSquaredSum_RU
(
    molecules.pairPotentials().nIds(),
    scalarField (mesh.nCells(), 0.0)
);

// Geometric Fields for IO

Info << nl << "Creating fields." << endl;

/*---------------------------------------------------------------------------*\
    Number density
\*---------------------------------------------------------------------------*/

PtrList<volScalarField> allSpeciesRhoN
(
    molecules.pairPotentials().nIds()
);

forAll (allSpeciesRhoN, rN)
{
    Info << "    Creating number density field for "
        << molecules.pairPotentials().idList()[rN] << endl;

    allSpeciesRhoN.set
    (
        rN,
        new volScalarField
        (
            IOobject
            (
                "rhoN_" + molecules.pairPotentials().idList()[rN],
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            mesh,
            dimless/dimVolume,
            "zeroGradient"
        )
    );
    allSpeciesRhoN[rN].internalField() = scalarField (mesh.nCells(), 0.0);
    allSpeciesRhoN[rN].correctBoundaryConditions();
}

Info << "    Creating total number density field" << endl;

volScalarField totalRhoN
(
    IOobject
    (
        "rhoN_total",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimless/dimVolume,
    "zeroGradient"
);
totalRhoN.internalField() = scalarField (mesh.nCells(), 0.0);
totalRhoN.correctBoundaryConditions();

/*---------------------------------------------------------------------------*\
    Mass density
\*---------------------------------------------------------------------------*/

PtrList<volScalarField> allSpeciesRhoM
(
    molecules.pairPotentials().nIds()
);

forAll (allSpeciesRhoM, rM)
{
    Info << "    Creating mass density field for "
        << molecules.pairPotentials().idList()[rM] << endl;

    allSpeciesRhoM.set
    (
        rM,
        new volScalarField
        (
            IOobject
            (
                "rhoM_" + molecules.pairPotentials().idList()[rM],
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            mesh,
            dimMass/dimVolume,
            "zeroGradient"
        )
    );
    allSpeciesRhoM[rM].internalField() = scalarField (mesh.nCells(), 0.0);
    allSpeciesRhoM[rM].correctBoundaryConditions();
}

Info << "    Creating total mass density field" << endl;

volScalarField totalRhoM
(
    IOobject
    (
        "rhoM_total",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimMass/dimVolume,
    "zeroGradient"
);
totalRhoM.internalField() = scalarField (mesh.nCells(), 0.0);
totalRhoM.correctBoundaryConditions();

/*---------------------------------------------------------------------------*\
    Bulk velocity
\*---------------------------------------------------------------------------*/

PtrList<volVectorField> allSpeciesVelocity
(
    molecules.pairPotentials().nIds()
);

forAll (allSpeciesVelocity, v)
{
    Info << "    Creating velocity field for "
        << molecules.pairPotentials().idList()[v] << endl;

    allSpeciesVelocity.set
    (
        v,
        new volVectorField
        (
            IOobject
            (
                "velocity_" + molecules.pairPotentials().idList()[v],
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            mesh,
            dimVelocity,
            "zeroGradient"
        )
    );
    allSpeciesVelocity[v].internalField() =
        vectorField (mesh.nCells(), vector::zero);
    allSpeciesVelocity[v].correctBoundaryConditions();
}

Info << "    Creating total velocity field" << endl;

// volVectorField totalVelocity
// (
//     IOobject
//     (
//         "velocity_total",
//         runTime.timeName(),
//         mesh,
//         IOobject::NO_READ,
//         IOobject::AUTO_WRITE
//     ),
//     mesh,
//     dimVelocity,
//     "zeroGradient"
// );
// totalVelocity.internalField() = vectorField (mesh.nCells(), vector::zero);
// totalVelocity.correctBoundaryConditions();


volVectorField totalVelocity

(
    IOobject
    (

    "velocity_total",
    runTime.timeName(),
    mesh,
    IOobject::NO_READ,
    IOobject::AUTO_WRITE

    ),
    mesh,
    dimensionedVector("zero", dimVelocity, vector::zero)
);

/*---------------------------------------------------------------------------*\
    Kinetic temperature
\*---------------------------------------------------------------------------*/

PtrList<volScalarField> allSpeciesTemperature
(
    molecules.pairPotentials().nIds()
);

forAll (allSpeciesTemperature, t)
{
    Info << "    Creating temperature field for "
        << molecules.pairPotentials().idList()[t] << endl;

    allSpeciesTemperature.set
    (
        t,
        new volScalarField
        (
            IOobject
            (
                "temperature_" + molecules.pairPotentials().idList()[t],
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            mesh,
            dimTemperature,
            "zeroGradient"
        )
    );
    allSpeciesTemperature[t].internalField() = scalarField (mesh.nCells(), 0.0);
    allSpeciesTemperature[t].correctBoundaryConditions();
}

Info << "    Creating total temperature field" << endl;

volScalarField totalTemperature
(
    IOobject
    (
        "temperature_total",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimTemperature,
    "zeroGradient"
);
totalTemperature.internalField() = scalarField (mesh.nCells(), 0.0);
totalTemperature.correctBoundaryConditions();

/*---------------------------------------------------------------------------*\
    Mean kinetic energy
\*---------------------------------------------------------------------------*/


PtrList<volScalarField> allSpeciesMeanKE
(
    molecules.pairPotentials().nIds()
);

forAll (allSpeciesMeanKE, mKE)
{
    Info << "    Creating mean kinetic energy field for "
        << molecules.pairPotentials().idList()[mKE] << endl;

    allSpeciesMeanKE.set
    (
        mKE,
        new volScalarField
        (
            IOobject
            (
                "meanKE_" + molecules.pairPotentials().idList()[mKE],
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            mesh,
            dimensionSet(1, 2, -2, 0, 0, 0, 0),
            "zeroGradient"
        )
    );
    allSpeciesMeanKE[mKE].internalField() = scalarField (mesh.nCells(), 0.0);
    allSpeciesMeanKE[mKE].correctBoundaryConditions();
}

Info << "    Creating total mean kinetic energy field" << endl;

volScalarField totalMeanKE
(
    IOobject
    (
        "meanKE_total",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionSet(1, 2, -2, 0, 0, 0, 0),
    "zeroGradient"
);
totalMeanKE.internalField() = scalarField (mesh.nCells(), 0.0);
totalMeanKE.correctBoundaryConditions();