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

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class CompType, class ThermoType>
Foam::FiniteRateChemistryReductionMethods::PFA_FiniteRate<CompType, ThermoType>::PFA_FiniteRate
(
    const IOdictionary& dict,
    TDACFiniteRateChemistryModel<CompType, ThermoType>& chemistry
)
  :
FiniteRateChemistryReductionMethod<CompType, ThermoType>(dict, chemistry),
searchInitSet_(this->coeffsDict_.subDict("initialSet").size())
{
  label j=0;

  dictionary initSet = this->coeffsDict_.subDict("initialSet");

  for (label i=0; i<chemistry.nSpecie(); i++) {
    if (initSet.found(chemistry.Y()[i].member())) {
      searchInitSet_[j++] = i;
    }
  }

  if (j<searchInitSet_.size()) {
    FatalErrorInFunction
        << searchInitSet_.size()-j
        << " species in the initial set is not in the mechanism "
        << initSet
        << abort(FatalError);
  }
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

template<class CompType, class ThermoType>
Foam::FiniteRateChemistryReductionMethods::PFA_FiniteRate<CompType, ThermoType>::~PFA_FiniteRate()
{}


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

template<class CompType, class ThermoType>
void Foam::FiniteRateChemistryReductionMethods::PFA_FiniteRate<CompType, ThermoType>::reduceMechanism
(
    const scalarField &c,
    const scalar T,
    const scalar p
)
{
  scalarField& completeC(this->chemistry_.completeC());
  scalarField c1(this->chemistry_.nEqns(), 0.0);

  for (label i=0; i<this->nSpecie_; i++) {
    c1[i] = c[i];
    completeC[i] = c[i];
  }

  c1[this->nSpecie_] = T;
  c1[this->nSpecie_+1] = p;

  // Compute the rAB matrix
  RectangularMatrix<scalar> PAB(this->nSpecie_,this->nSpecie_,0.0);
  RectangularMatrix<scalar> CAB(this->nSpecie_,this->nSpecie_,0.0);
  scalarField PA(this->nSpecie_,0.0);
  scalarField CA(this->nSpecie_,0.0);

  // Number of initialized rAB for each lines
  Field<label> NbrABInit(this->nSpecie_,0);
  // Position of the initialized rAB, -1 when not initialized
  RectangularMatrix<label> rABPos(this->nSpecie_, this->nSpecie_, -1);
  // Index of the other species involved in the rABNum
  RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);

  scalar pf, cf, pr, cr;
  label lRef, rRef;
  forAll(this->chemistry_.reactions(), i) {
    const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
    // for each reaction compute omegai
    scalar omegai = this->chemistry_.omega
                    (
                        R, c1, T, p, pf, cf, lRef, pr, cr, rRef
                    );

    // then for each pair of species composing this reaction,
    // compute the rAB matrix (separate the numerator and
    // denominator)

    DynamicList<scalar> wA(R.lhs().size()+R.rhs().size());
    DynamicList<label> wAID(R.lhs().size()+R.rhs().size());

    forAll(R.lhs(), s) { // compute rAB for all species in the left hand side
      label ss = R.lhs()[s].index;
      scalar sl = -R.lhs()[s].stoichCoeff; // vAi = v''-v' => here -v'
      bool found(false);
      forAll(wAID, id) {
        if (ss==wAID[id]) {
          wA[id] += sl*omegai;
          found = true;
          break;
        }
      }
      if (!found) {
        wA.append(sl*omegai);
        wAID.append(ss);
      }
    }
    forAll(R.rhs(), s) { // compute rAB for all species in the right hand side
      label ss = R.rhs()[s].index;
      scalar sl = R.rhs()[s].stoichCoeff; // vAi = v''-v' => here v''
      bool found(false);
      forAll(wAID, id) {
        if (ss==wAID[id]) {
          wA[id] += sl*omegai;
          found = true;
          break;
        }
      }
      if (!found) {
        wA.append(sl*omegai);
        wAID.append(ss);
      }
    }

    wAID.shrink();
    forAll(wAID, id) {
      label curID = wAID[id];
      scalar curwA = wA[id];
      List<bool> deltaBi(this->nSpecie_, false);
      FIFOStack<label> usedIndex;
      forAll(R.lhs(),j) {
        label sj = R.lhs()[j].index;
        usedIndex.push(sj);
        deltaBi[sj] = true;
      }
      forAll(R.rhs(),j) {
        label sj = R.rhs()[j].index;
        usedIndex.push(sj);
        deltaBi[sj] = true;
      }

      deltaBi[curID] = false;

      while(!usedIndex.empty()) {
        label curIndex = usedIndex.pop();

        if (deltaBi[curIndex]) {
          deltaBi[curIndex] = false;
          if (rABPos(curID, curIndex)==-1) {
            rABPos(curID, curIndex) = NbrABInit[curID];
            rABOtherSpec(curID, NbrABInit[curID]) = curIndex;
            if (curwA > 0.0) {
              PAB(curID, NbrABInit[curID]) = curwA;
            } else {
              CAB(curID, NbrABInit[curID]) = -curwA;
            }
            NbrABInit[curID]++;
          } else {
            if (curwA > 0.0) {
              PAB(curID, rABPos(curID, curIndex)) += curwA;
            } else {
              CAB(curID, rABPos(curID, curIndex)) += -curwA;
            }
          }
        }
      }
      // Now that every species of the reactions has been visited, we can
      // compute the production and consumption rate. It avoids getting
      // wrong results when species are present in both lhs and rhs
      if (curwA > 0.0) {
        if (PA[curID] == 0.0) {
          PA[curID] = curwA;
        } else {
          PA[curID] += curwA;
        }
      } else {
        if (CA[curID] == 0.0) {
          CA[curID] = -curwA;
        } else {
          CA[curID] += -curwA;
        }
      }
    }
  }
  // rii = 0.0 by definition

  // Compute second generation link strength
  // For all species A look at all rAri of the connected species ri and
  // compute rriB with all the connected species of ri, B different of A.  If
  // a new species is connected, add it to the list of connected species.  It
  // is a connection of second generation and it will be aggregated in the
  // final step to evaluate the total connection strength (or path flux).
  // Compute rsecond=rAri*rriB with A!=ri!=B
  RectangularMatrix<scalar> PAB2nd(this->nSpecie_,this->nSpecie_,0.0);
  RectangularMatrix<scalar> CAB2nd(this->nSpecie_,this->nSpecie_,0.0);

  // Number of initialized rAB for each lines
  Field<label> NbrABInit2nd(this->nSpecie_, 0);

  // Position of the initialized rAB, -1 when not initialized
  RectangularMatrix<label> rABPos2nd(this->nSpecie_, this->nSpecie_, -1);

  // Index of the other species involved in the rABNum
  RectangularMatrix<label> rABOtherSpec2nd
  (
      this->nSpecie_, this->nSpecie_, -1
  );

  forAll(NbrABInit, A) {
    for (int i=0; i<NbrABInit[A]; i++) {
      label ri = rABOtherSpec(A, i);
      scalar maxPACA = max(PA[ri],CA[ri]);
      if (maxPACA > vSmall) {
        for (int j=0; j<NbrABInit[ri]; j++) {
          label B = rABOtherSpec(ri, j);
          if (B != A) { // if B!=A and also !=ri by definition
            if (rABPos2nd(A, B)==-1) {
              rABPos2nd(A, B) = NbrABInit2nd[A]++;
              rABOtherSpec2nd(A, rABPos2nd(A, B)) = B;
              PAB2nd(A, rABPos2nd(A, B)) =
                  PAB(A, i)*PAB(ri, j)/maxPACA;
              CAB2nd(A, rABPos2nd(A, B)) =
                  CAB(A, i)*CAB(ri, j)/maxPACA;
            } else {
              PAB2nd(A, rABPos2nd(A, B)) +=
                  PAB(A, i)*PAB(ri, j)/maxPACA;
              CAB2nd(A, rABPos2nd(A, B)) +=
                  CAB(A, i)*CAB(ri, j)/maxPACA;
            }
          }
        }
      }
    }
  }

  // Using the rAB matrix (numerator and denominator separated)
  label speciesNumber = 0;

  // set all species to inactive and activate them according
  // to rAB and initial set
  for (label i=0; i<this->nSpecie_; i++) {
    this->activeSpecies_[i] = false;
  }

  // Initialize the FIFOStack for search set
  const labelList& SIS(this->searchInitSet_);
  FIFOStack<label> Q;

  for (label i=0; i<SIS.size(); i++) {
    label q = SIS[i];
    this->activeSpecies_[q] = true;
    speciesNumber++;
    Q.push(q);
  }

  // Execute the main loop for R-value
  while (!Q.empty()) {
    label u = Q.pop();
    scalar Den = max(PA[u],CA[u]);

    if (Den!=0.0) {
      // first generation
      for (label v=0; v<NbrABInit[u]; v++) {
        label otherSpec = rABOtherSpec(u, v);
        scalar rAB = (PAB(u, v)+CAB(u, v))/Den; // first generation
        label id2nd = rABPos2nd(u, otherSpec);
        if (id2nd !=-1) { // if there is a second generation link
          rAB += (PAB2nd(u, id2nd)+CAB2nd(u, id2nd))/Den;
        }
        // the link is stronger than the user-defined tolerance
        if
        (
            rAB >= this->tolerance()
            && !this->activeSpecies_[otherSpec]
        ) {
          Q.push(otherSpec);
          this->activeSpecies_[otherSpec] = true;
          speciesNumber++;
        }

      }
      // second generation link only (for those without first link)
      for (label v=0; v<NbrABInit2nd[u]; v++) {
        label otherSpec = rABOtherSpec2nd(u, v);
        scalar rAB = (PAB2nd(u, v)+CAB2nd(u, v))/Den;
        // the link is stronger than the user-defined tolerance
        if
        (
            rAB >= this->tolerance()
            && !this->activeSpecies_[otherSpec]
        ) {
          Q.push(otherSpec);
          this->activeSpecies_[otherSpec] = true;
          speciesNumber++;
        }
      }
    }
  }

  // Put a flag on the reactions containing at least one removed species
  forAll(this->chemistry_.reactions(), i) {
    const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
    this->chemistry_.reactionsDisabled()[i] = false;

    forAll(R.lhs(), s) {
      label ss = R.lhs()[s].index;
      if (!this->activeSpecies_[ss]) {
        this->chemistry_.reactionsDisabled()[i] = true;
        break;
      }
    }
    if (!this->chemistry_.reactionsDisabled()[i]) {
      forAll(R.rhs(), s) {
        label ss = R.rhs()[s].index;
        if (!this->activeSpecies_[ss]) {
          this->chemistry_.reactionsDisabled()[i] = true;
          break;
        }
      }
    }
  }

  this->NsSimp_ = speciesNumber;
  scalarField& simplifiedC(this->chemistry_.simplifiedC());
  simplifiedC.setSize(this->NsSimp_+2);
  DynamicList<label>& s2c(this->chemistry_.simplifiedToCompleteIndex());
  s2c.setSize(this->NsSimp_);
  Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());

  label j = 0;
  for (label i=0; i<this->nSpecie_; i++) {
    if (this->activeSpecies_[i]) {
      s2c[j] = i;
      simplifiedC[j] = c[i];
      c2s[i] = j++;
      if (!this->chemistry_.active(i)) {
        this->chemistry_.setActive(i);
      }
    } else {
      c2s[i] = -1;
    }
  }
  simplifiedC[this->NsSimp_] = T;
  simplifiedC[this->NsSimp_+1] = p;
  this->chemistry_.setNsDAC(this->NsSimp_);
  // change temporary Ns in chemistryModel
  // to make the function nEqns working
  this->chemistry_.setNSpecie(this->NsSimp_);
}


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