/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2022 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 "polyTopoChangeMap.H"
#include "processorFvPatchField.H"

// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

template<class GeoField>
void Foam::fvMeshDistribute::printFieldInfo(const fvMesh& mesh)
{
    HashTable<const GeoField*> flds
    (
        mesh.objectRegistry::lookupClass<GeoField>()
    );

    forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
    {
        const GeoField& fld = *iter();

        Pout<< "Field:" << iter.key() << " internal size:" << fld.size()
            << endl;

        forAll(fld.boundaryField(), patchi)
        {
            Pout<< "    " << patchi
                << ' ' << fld.boundaryField()[patchi].patch().name()
                << ' ' << fld.boundaryField()[patchi].type()
                << ' ' << fld.boundaryField()[patchi].size()
                << endl;
        }
    }
}


template<class T, class Mesh>
void Foam::fvMeshDistribute::saveBoundaryFields
(
    PtrList<FieldField<fvsPatchField, T>>& bflds
) const
{
    // Save whole boundary field

    typedef GeometricField<T, fvsPatchField, Mesh> fldType;

    HashTable<const fldType*> flds
    (
        static_cast<const fvMesh&>(mesh_).objectRegistry::lookupClass<fldType>()
    );

    bflds.setSize(flds.size());

    label i = 0;

    forAllConstIter(typename HashTable<const fldType*>, flds, iter)
    {
        const fldType& fld = *iter();

        bflds.set(i, fld.boundaryField().clone().ptr());

        i++;
    }
}


template<class T, class Mesh>
void Foam::fvMeshDistribute::mapBoundaryFields
(
    const polyTopoChangeMap& map,
    const PtrList<FieldField<fvsPatchField, T>>& oldBflds
)
{
    // Map boundary field

    const labelList& oldPatchStarts = map.oldPatchStarts();
    const labelList& faceMap = map.faceMap();

    typedef GeometricField<T, fvsPatchField, Mesh> fldType;

    HashTable<fldType*> flds
    (
        mesh_.objectRegistry::lookupClass<fldType>()
    );

    if (flds.size() != oldBflds.size())
    {
        FatalErrorInFunction
            << abort(FatalError);
    }

    label fieldi = 0;

    forAllIter(typename HashTable<fldType*>, flds, iter)
    {
        fldType& fld = *iter();
        typename fldType::Boundary& bfld =
            fld.boundaryFieldRef();

        const FieldField<fvsPatchField, T>& oldBfld = oldBflds[fieldi++];

        // Pull from old boundary field into bfld.

        forAll(bfld, patchi)
        {
            fvsPatchField<T>& patchFld = bfld[patchi];
            label facei = patchFld.patch().start();

            forAll(patchFld, i)
            {
                label oldFacei = faceMap[facei++];

                // Find patch and local patch face oldFacei was in.
                forAll(oldPatchStarts, oldPatchi)
                {
                    label oldLocalI = oldFacei - oldPatchStarts[oldPatchi];

                    if (oldLocalI >= 0 && oldLocalI < oldBfld[oldPatchi].size())
                    {
                        patchFld[i] = oldBfld[oldPatchi][oldLocalI];
                    }
                }
            }
        }
    }
}


template<class T>
void Foam::fvMeshDistribute::initMapExposedFaces
(
    PtrList<Field<T>>& iflds
) const
{
    HashTable<const SurfaceField<T>*> flds
    (
        static_cast<const fvMesh&>(mesh_).lookupClass<SurfaceField<T>>()
    );

    iflds.setSize(flds.size());

    label fieldi = 0;

    forAllConstIter(typename HashTable<const SurfaceField<T>*>, flds, iter)
    {
        iflds.set(fieldi, Field<T>(mesh_.nFaces()));

        const SurfaceField<T>& fld = *iter();

        SubList<T>(iflds[fieldi], fld.primitiveField().size()) =
            fld.primitiveField();

        forAll(fld.boundaryField(), patchi)
        {
            const fvsPatchField<T>& pfld = fld.boundaryField()[patchi];

            SubList<T>(iflds[fieldi], pfld.size(), pfld.patch().start()) =
                pfld;
        }

        fieldi++;
    }
}


template<class T>
void Foam::fvMeshDistribute::mapExposedFaces
(
    const polyTopoChangeMap& map,
    const PtrList<Field<T>>& oldFlds
)
{
    HashTable<SurfaceField<T>*> flds
    (
        mesh_.objectRegistry::lookupClass<SurfaceField<T>>()
    );

    label fieldi = 0;

    forAllIter(typename HashTable<SurfaceField<T>*>, flds, iter)
    {
        SurfaceField<T>& fld = *iter();

        const Field<T>& oldFld = oldFlds[fieldi];

        const bool negateIfFlipped = isFlux(fld);

        forAll(fld.boundaryField(), patchi)
        {
            fvsPatchField<T>& patchFld = fld.boundaryFieldRef()[patchi];

            forAll(patchFld, i)
            {
                const label facei = patchFld.patch().start()+i;
                const label oldFacei = map.faceMap()[facei];

                if (oldFacei < map.nOldInternalFaces())
                {
                    if (negateIfFlipped && map.flipFaceFlux().found(facei))
                    {
                        patchFld[i] = flipOp()(oldFld[oldFacei]);
                    }
                    else
                    {
                        patchFld[i] = oldFld[oldFacei];
                    }
                }
                else
                {
                    patchFld[i] = oldFld[oldFacei];
                }
            }
        }

        fieldi++;
    }
}


template<class GeoField>
void Foam::fvMeshDistribute::correctProcessorPatchFields()
{
    typedef processorFvPatchField<typename GeoField::value_type>
        processorPatchFieldType;

    HashTable<GeoField*> flds
    (
        mesh_.objectRegistry::lookupClass<GeoField>()
    );

    forAllIter(typename HashTable<GeoField*>, flds, iter)
    {
        GeoField& fld = *iter();

        typename GeoField::Boundary& bfld = fld.boundaryFieldRef();

        if
        (
            Pstream::defaultCommsType == Pstream::commsTypes::blocking
         || Pstream::defaultCommsType == Pstream::commsTypes::nonBlocking
        )
        {
            label nReq = Pstream::nRequests();

            forAll(bfld, patchi)
            {
                if (bfld[patchi].coupled() || isA<processorPatchFieldType>(bfld[patchi]))
                {
                    bfld[patchi].initEvaluate(Pstream::defaultCommsType);
                }
            }

            // Block for any outstanding requests
            if
            (
                Pstream::parRun()
             && Pstream::defaultCommsType == Pstream::commsTypes::nonBlocking
            )
            {
                Pstream::waitRequests(nReq);
            }

            forAll(bfld, patchi)
            {
                if (bfld[patchi].coupled() || isA<processorPatchFieldType>(bfld[patchi]))
                {
                    bfld[patchi].evaluate(Pstream::defaultCommsType);
                }
            }
        }
        else if (Pstream::defaultCommsType == Pstream::commsTypes::scheduled)
        {
            const lduSchedule& patchSchedule =
                mesh_.globalData().patchSchedule();

            forAll(patchSchedule, patchEvali)
            {
                if (bfld[patchEvali].coupled() || bfld[patchEvali].coupled())
                {
                    if (patchSchedule[patchEvali].init)
                    {
                        bfld[patchSchedule[patchEvali].patch]
                            .initEvaluate(Pstream::commsTypes::scheduled);
                    }
                    else
                    {
                        bfld[patchSchedule[patchEvali].patch]
                            .evaluate(Pstream::commsTypes::scheduled);
                    }
                }
            }
        }
    }
}


template<class GeoField>
void Foam::fvMeshDistribute::sendFields
(
    const label domain,
    const wordList& fieldNames,
    const fvMeshSubset& subsetter,
    Ostream& toNbr
)
{
    // Send fields. Note order supplied so we can receive in exactly the same
    // order.
    // Note that field gets written as entry in dictionary so we
    // can construct from subdictionary.
    // (since otherwise the reading as-a-dictionary mixes up entries from
    // consecutive fields)
    // The dictionary constructed is:
    //  volScalarField
    //  {
    //      p {internalField ..; boundaryField ..;}
    //      k {internalField ..; boundaryField ..;}
    //  }
    //  volVectorField
    //  {
    //      U {internalField ...  }
    //  }

    // volVectorField {U {internalField ..; boundaryField ..;}}

    toNbr << GeoField::typeName << token::NL << token::BEGIN_BLOCK << token::NL;
    forAll(fieldNames, i)
    {
        if (debug)
        {
            Pout<< "Subsetting field " << fieldNames[i]
                << " for domain:" << domain << endl;
        }

        // Send all fieldNames. This has to be exactly the same set as is
        // being received!
        const GeoField& fld =
            subsetter.baseMesh().lookupObject<GeoField>(fieldNames[i]);

        tmp<GeoField> tsubfld = subsetter.interpolate(fld);

        toNbr
            << fieldNames[i] << token::NL << token::BEGIN_BLOCK
            << tsubfld
            << token::NL << token::END_BLOCK << token::NL;
    }
    toNbr << token::END_BLOCK << token::NL;
}


template<class GeoField>
void Foam::fvMeshDistribute::receiveFields
(
    const label domain,
    const wordList& fieldNames,
    typename GeoField::Mesh& mesh,
    PtrList<GeoField>& fields,
    const dictionary& fieldDicts
)
{
    if (debug)
    {
        Pout<< "Receiving fields " << fieldNames
            << " from domain:" << domain << endl;
    }

    fields.setSize(fieldNames.size());

    forAll(fieldNames, i)
    {
        if (debug)
        {
            Pout<< "Constructing field " << fieldNames[i]
                << " from domain:" << domain << endl;
        }

        fields.set
        (
            i,
            new GeoField
            (
                IOobject
                (
                    fieldNames[i],
                    mesh.thisDb().time().timeName(),
                    mesh.thisDb(),
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
                mesh,
                fieldDicts.subDict(fieldNames[i])
            )
        );
    }
}


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

Create time

Create mesh for time = 0

Selecting fvMeshDistributor loadBalancer
Selecting distributor zoltan
Selecting solver multicomponentFluid
Selecting thermodynamics package 
{
    type            hePsiThermo;
    mixture         multiComponentMixture;
    transport       sutherland;
    thermo          janaf;
    energy          sensibleEnthalpy;
    equationOfState perfectGas;
    specie          specie;
}

Selecting turbulence model type laminar
Selecting laminar stress model Stokes
No MRF models present

Courant Number mean: 2.73409e-05 max: 0.000175681
Selecting combustion model laminar
Selecting chemistry solver 
{
    solver          ode;
    method          chemistryModel;
}

odeChemistryModel: Number of species = 5
chemistryModel: Number of species = 5 and reactions = 1
Selecting ODE solver seulex
    using integrated reaction rate
Selecting thermophysical transport type laminar
Selecting default laminar thermophysical transport model unityLewisFourier

PIMPLE: No convergence criteria found


PIMPLE: Operating solver in transient mode with 1 outer corrector
PIMPLE: Operating solver in PISO mode



Starting time loop

Courant Number mean: 2.73409e-05 max: 0.000175681
deltaT = 1.11111e-06
Time = 1.11111e-06s

diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab:  Solving for O2, Initial residual = 1, Final residual = 1.26231e-10, No Iterations 1
DILUPBiCGStab:  Solving for H2O, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab:  Solving for CH4, Initial residual = 1, Final residual = 1.26283e-10, No Iterations 1
DILUPBiCGStab:  Solving for CO2, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab:  Solving for h, Initial residual = 1, Final residual = 1.2625e-10, No Iterations 1
DICPCG:  Solving for p, Initial residual = 1, Final residual = 0.0291404, No Iterations 1
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 1.242e-06, global = -6.6435e-07, cumulative = -6.6435e-07
DICPCG:  Solving for p, Initial residual = 0.00828587, Final residual = 3.63402e-07, No Iterations 4
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 5.44339e-11, global = -2.80532e-11, cumulative = -6.64378e-07
ExecutionTime = 0.097289 s  ClockTime = 0 s

Courant Number mean: 3.71025e-05 max: 0.000195819
deltaT = 1.26984e-06
Time = 2.38095e-06s

diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab:  Solving for O2, Initial residual = 0.240887, Final residual = 4.03063e-11, No Iterations 1
DILUPBiCGStab:  Solving for H2O, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab:  Solving for CH4, Initial residual = 0.240844, Final residual = 4.0334e-11, No Iterations 1
DILUPBiCGStab:  Solving for CO2, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab:  Solving for h, Initial residual = 0.296454, Final residual = 1.4082e-10, No Iterations 1
DICPCG:  Solving for p, Initial residual = 0.300101, Final residual = 0.0178395, No Iterations 1
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 3.55422e-06, global = -2.44847e-06, cumulative = -3.11285e-06
DICPCG:  Solving for p, Initial residual = 0.0114237, Final residual = 1.12274e-07, No Iterations 5
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 3.48946e-11, global = -1.97666e-11, cumulative = -3.11287e-06
ExecutionTime = 0.099346 s  ClockTime = 0 s

Courant Number mean: 6.06108e-05 max: 0.000225148
[3] imbalance  0.000257 0.010254 0.010511 0.0105128
[1] imbalance  0.000277 0.010254 0.010531 0.0105128
[2] imbalance  0.000248 0.010254 0.010502 0.0105128
[0] imbalance  0.000253 0.010254 0.010507 0.0105128
Redistributing mesh with imbalance 0.00173599
deltaT = 1.52381e-06
Time = 3.90476e-06s

diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0

[3] #0  Foam::error::printStack(Foam::Ostream&)[0] #0  Foam::error::printStack(Foam::Ostream&)[1] #0  Foam::error::printStack(Foam::Ostream&)[2] #0  Foam::error::printStack(Foam::Ostream&) at ??:?
 at ??:?
 at ??:?
 at ??:?
[1] #1  Foam::sigFpe::sigHandler(int)[2] #1  Foam::sigFpe::sigHandler(int)[3] #1  Foam::sigFpe::sigHandler(int)[0] #1  Foam::sigFpe::sigHandler(int) at ??:?
[1] #2  ? at ??:?
[2] #2  ? at ??:?
[3] #2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
 in "/lib/x86_64-linux-gnu/libc.so.6"
[2] #3  [1] #3  Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&)Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
[0] #2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
[3] #3  Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) in "/lib/x86_64-linux-gnu/libc.so.6"
[0] #3  Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
 at ??:?
[3] #4  [2] #4  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&)Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?
[1] #4  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?
 at ??:?
 at ??:?
[3] #5  Foam::fluidThermo::nu() const[2] #5  Foam::fluidThermo::nu() const[1] #5  Foam::fluidThermo::nu() const at ??:?
[0] #4  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?
[0] #5  Foam::fluidThermo::nu() const at ??:?
[2] #6  Foam::laminarModels::Stokes<Foam::compressibleMomentumTransportModel>::nuEff() const at ??:?
 at ??:?
[1] #[3] #6  Foam::laminarModels::Stokes<Foam::compressibleMomentumTransportModel>::nuEff() const6  Foam::laminarModels::Stokes<Foam::compressibleMomentumTransportModel>::nuEff() const at ??:?
[0] #6  Foam::laminarModels::Stokes<Foam::compressibleMomentumTransportModel>::nuEff() const at ??:?
 at ??:?
 at ??:?
 at ??:?
[3] #7  [1] #7  [0] #7  Foam::linearViscousStress<Foam::laminarModel<Foam::compressibleMomentumTransportModel> >::divDevTau(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) constFoam::linearViscousStress<Foam::laminarModel<Foam::compressibleMomentumTransportModel> >::divDevTau(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) constFoam::linearViscousStress<Foam::laminarModel<Foam::compressibleMomentumTransportModel> >::divDevTau(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) const[2] #7  Foam::linearViscousStress<Foam::laminarModel<Foam::compressibleMomentumTransportModel> >::divDevTau(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) const at ??:?
[0] #8  Foam::solvers::isothermalFluid::momentumPredictor() at ??:?
[3] #8  Foam::solvers::isothermalFluid::momentumPredictor() at ??:?
[1] #8  Foam::solvers::isothermalFluid::momentumPredictor() at ??:?
[2] #8  Foam::solvers::isothermalFluid::momentumPredictor() at ??:?
[0] #9   at ??:?
[3] #9   at ??:?
[1] #9   at ??:?
...
Floating point exception
Floating point exception
Floating point exception
--------------------------------------------------------------------------
Primary job  terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
Floating point exception
--------------------------------------------------------------------------
mpirun detected that one or more processes exited with non-zero status, thus causing
the job to be terminated. The first process to do so was:

  Process name: [[8152,1],1]
  Exit code:    136