/** @file Qmicroturbine.cc
    @brief Main file

    This file is part of LIBPF
    All rights reserved; do not distribute without permission.
    @author (C) Copyright 2009-2010 Paolo Greppi 3iP

    Repository revision: $Rev: 551 $
    File revision date: $Date: 2010-05-15 12:53:03 +0200 (sab, 15 mag 2010) $
 */

static const int verbositySubSystem = 0;

#include "models.h"
#include "stream.h"
#include "graph.h"
#include "usermodels.h"
#include "purecomps.h"

#include <boost/assign/list_of.hpp>
#include <boost/assign/list_inserter.hpp>

void setComponents(void) {
  components.addcomp(new purecomps::water);
  components.addcomp(new purecomps::N2);
  components.addcomp(new purecomps::O2);
  components.addcomp(new purecomps::methane);
  components.addcomp(new purecomps::CO);
  components.addcomp(new purecomps::CO2);
  components.addcomp(new purecomps::H2);
} // void setComponents(void)

// /////////////////////////////////////////////////////////////////////////////
// ///////////////////////////////// microturbo ////////////////////////////////
// /////////////////////////////////////////////////////////////////////////////

#define standalone

#ifdef standalone
#define AAA flowsheet<zero_zero>
#else
#define AAA flowsheet<vertex<3,2> >
#endif

class microturbo : public modelBase, public AAA {
private:
  static std::string type_;
  void maketables(long);
public:
  microturbo(const std::string &t, const std::string &d) : persistent_item_interface(-1), persistent_item(-1, t, d), modelBaseInterface(-1), modelBase(t, d) {
    maketables(-1);
  };
  microturbo(long cid) : persistent_item_interface(cid), persistent_item(cid, "", ""), modelBaseInterface(cid), modelBase(cid), AAA(cid) {
    maketables(cid);
    readtables(cid);
  }

  Qdouble LHV_CH4;      ///< methane Lower Heating Value
  Qdouble LHV_H2;       ///< hydrogen Lower Heating Value
  Qdouble LHV_CO;       ///< carbon monoxide Lower Heating Value

  Qdouble powerIn;      ///< Lower Heating Value input
  Qdouble power;        ///< Gas Turbine Subsystem power output
  Qdouble powerEl;      ///< Total net electric output
  Qdouble etaEl;        ///< Overall electric efficiency

  Qdouble beta0;        ///< compressor compression ratio (cr) at the nominal point
  Qdouble mC0;          ///< compressor mass flow at the nominal point (S01)
  Qdouble TC0;          ///< compressor inlet temperature at the nominal point (S01)
  Qdouble PC0;          ///< compressor inlet pressure at the nominal point (S01)
  Qdouble mrC0;         ///< compressor corrected mass flow at the nominal point
  Qdouble mT0;          ///< turbine mass flow at the nominal point (S05)
  Qdouble PT0;          ///< turbine inlet pressure at the nominal point (S05)
  Qdouble mrT0;         ///< turbine corrected mass flow at the nominal point

  Qdouble TIT;          ///< should be Turbine Inlet Temperature
  Qdouble mC;           ///< should-be compressor uncorrected mass flow
  Qdouble mT;           ///< should-be turbine uncorrected mass flow
  Qdouble maC;          ///< actual compressor uncorrected mass flow
  Qdouble maT;          ///< actual turbine uncorrected mass flow

  Qdouble nr;           ///< compressor/turbine reduced frequency

  void makeuserassembly(std::list<assignment *>::iterator &p);
  void setup(void);
  void calculate(int MAXITS = 100, int level = 0);
  const std::string &type(void) const {
    return type_;
  }

  int smiter_precalc(void) { return 5; }
  bool supports_simultaneous(void) { return true; };
  int maxit(void) { return 200; }
}; // microturbo

std::string microturbo::type_("microturbo");

void microturbo::maketables(long cid) {
  static const int verbosityLocal = 0;
  try {
    diagnostic(2, " entered for " << tag());

    // Non embedded objects
    MAKEQDOUBLE(LHV_CH4, 802618.0, "kJ/kmol", "methane Lower Heating Value");
    MAKEQDOUBLE(LHV_H2, 241814.0, "kJ/kmol", "hydrogen Lower Heating Value");
    MAKEQDOUBLE(LHV_CO, 282980.0, "kJ/kmol", "carbon monoxide Lower Heating Value");

    MAKEQDOUBLE(powerIn, 0.0, "W", "Lower Heating Value input");

    MAKEQDOUBLE(power, 0.0, "W", "Gas Turbine Subsystem power output");
    MAKEQDOUBLE(powerEl, 0.0, "W", "Total net electric output");
    MAKEQDOUBLE(etaEl, 0.0, "", "Overall electric efficiency");

    MAKEQDOUBLE(beta0, 3.6, "", "compressor compression ratio at the nominal point");
    MAKEQDOUBLE(mC0, 1.66667, "kg/s", "compressor mass flow at the nominal point");
    MAKEQDOUBLE(TC0, 296.15, "K", "compressor inlet temperature at the nominal point");
    MAKEQDOUBLE(PC0, 101325., "Pa", "compressor inlet pressure at the nominal point");
    MAKEQDOUBLE(mrC0, 0.0, "m*s", "compressor corrected mass flow at the nominal point");
    MAKEQDOUBLE(mT0, 1.87407, "kg/s", "turbine mass flow at the nominal point");
    MAKEQDOUBLE(PT0, 347676., "Pa", "turbine inlet pressure at the nominal point");
    MAKEQDOUBLE(mrT0, 0.0, "m*s", "turbine corrected mass flow at the nominal point");

    MAKEQDOUBLE(TIT, t0, "K", "should be Turbine Inlet Temperature");
    MAKEQDOUBLE(mC, 0.0, "kg/s", "should-be compressor uncorrected mass flow");
    MAKEQDOUBLE(mT, 0.0, "kg/s", "should-be turbine uncorrected mass flow");
    MAKEQDOUBLE(maC, 0.0, "kg/s", "actual compressor uncorrected mass flow");
    MAKEQDOUBLE(maT, 0.0, "kg/s", "actual turbine uncorrected mass flow");

    MAKEQDOUBLE(nr, 0.0, "", "compressor and turbine reduced frequency");

    if (cid == -1) {
      diagnostic(2, "Define unit operations");
      makeVertex<mixer>("MIX", "Mixes fuel preheated air and FCS purge", cid);
      makeVertex<splitter>("SPLIT", "Compressor outlet splitter", cid, boost::assign::map_list_of<std::string, long>("nStreams", 2));

      makeVertex<compressor>("T", "Turbine", cid);
      makeVertex<compressor>("C", "Compressor", cid);

      std::list<std::string> rHC = boost::assign::list_of("reactionCH4_TC")
                                               ("reactionCO_TC")
                                               ("reactionH2_TC");
      makeVertex<genflashN1>("RX", "Burner", rHC, cid);
      makeVertex<hx>("HX", "Turbine recuperator", cid);

      diagnostic(2, "Define stream objects and connect");
      makeEdge<stream_V>("S01", "Air to GTS", "source", "out", "C", "in", cid);
      makeEdge<stream_V>("S02", "NG to GTS burner", "source", "out", "RX", "in", cid);
      makeEdge<stream_V>("S03", "Hot air recovery", "source", "out", "MIX", "in", cid);

      makeEdge<stream_V>("S04", "RX burner inlet", "MIX", "out", "RX", "in", cid);
      makeEdge<stream_V>("S05", "Turbine Inlet", "RX", "out", "T", "in", cid);
      makeEdge<stream_V>("S06", "Turbine Outlet", "T", "out", "HX", "hotin", cid);
      makeEdge<stream_V>("S07", "Compressed Air", "C", "out", "SPLIT", "in", cid);
      makeEdge<stream_V>("S08", "HX cold inlet", "SPLIT", "out2", "HX", "coldin", cid);
      makeEdge<stream_V>("S09", "HX cold outlet", "HX", "coldout", "MIX", "in", cid);

      makeEdge<stream_V>("S10", "Air to FCS", "SPLIT", "out1", "sink", "in", cid);
      makeEdge<stream_V>("S11", "Recuperator exhaust", "HX", "hotout", "sink", "in", cid);
    }
  } // try
  catch (error &e) {
    e.append(CURRENT_FUNCTION);
    throw;
  } // catch
} // microturbo::maketables

void microturbo::setup(void) {
  static const int verbosityLocal = 0;
  try {
    diagnostic(2, " entered for " << tag());

    double Pset = 4.0; // max pressure in bar

    diagnostic(3, "Calling flowsheet::setup to initialize any embedded flowsheet");
    AAA::setup();

    diagnostic(3, "Setting input variables");
    beta0 = 3.94769306686405;
    mC0 = Qdouble(0.306891411, "kg/s");
    TC0 = Qdouble(288.15, "K");
    PC0 = Qdouble(101325., "Pa");
    mT0 = Qdouble(0.309236235, "kg/s");
    PT0 = Qdouble(385000., "Pa");

    TIT = Qdouble(1144.0, "K");

    mrC0 = mC0 * sqrt(TC0 / T0) / PC0;
    mrT0 = mT0 * sqrt(TIT / T0) / PT0;

#ifdef standalone
    // fresh air
    O("S01")->Q("T")->set(15.0 + 273.15, "K");
    O("S01")->Q("P")->set(1.0, "atm");
    O("S01")->S("flowoption")->set("Nx");
    my_cast<stream *>(O("S01"), CURRENT_FUNCTION)->clearcomposition();
    O("S01")->O("Tphase")->Q("ndot")->set(38.49159293, "kmol/h");
    O("S01")->O("Tphase")->Q("x", "N2")->set(0.792);
    O("S01")->O("Tphase")->Q("x", "O2")->set(0.198);
    O("S01")->O("Tphase")->Q("x", "H2O")->set(0.01);

    // feed natural gas to turbine
    O("S02")->Q("T")->set(15.0 + 273.15, "K");
    O("S02")->Q("P")->set(Pset, "bar");
    O("S02")->S("flowoption")->set("Nx");
    my_cast<stream *>(O("S02"), CURRENT_FUNCTION)->clearcomposition();
    O("S02")->O("Tphase")->Q("ndot")->set(0.522280634, "kmol/h");
    O("S02")->O("Tphase")->Q("x", "CH4")->set(0.99);
    O("S02")->O("Tphase")->Q("x", "N2")->set(0.01);

    // hot recovered air
    O("S03")->Q("T")->set(15.0 + 273.15, "K");
    O("S03")->Q("P")->set(1.0, "atm");
    O("S03")->S("flowoption")->set("Nx");
    my_cast<stream *>(O("S03"), CURRENT_FUNCTION)->clearcomposition();
    O("S03")->O("Tphase")->Q("ndot")->set(0.0, "kmol/h");
    O("S03")->O("Tphase")->Q("x", "N2")->set(0.792);
    O("S03")->O("Tphase")->Q("x", "O2")->set(0.198);
    O("S03")->O("Tphase")->Q("x", "H2O")->set(0.01);
#endif // standalone

    // Burner
    O("RX")->Q("deltaP")->set(100.0, "mbar");
    O("RX")->Q("duty")->set(0.0, "W");
    O("RX")->S("option")->set("DH");
    // total combustion of methane
    O("RX")->O("reactions", 0)->Q("z")->set(1.0);
    // total combustion of CO
    O("RX")->O("reactions", 1)->Q("z")->set(1.0);
    // total combustion of H2
    O("RX")->O("reactions", 2)->Q("z")->set(1.0);

    // Compressor outlet splitter, 1-(stream16/stream15)
    O("SPLIT")->Q("outSplit", 0)->set(0.0);
    O("SPLIT")->Q("outSplit", 1)->set(1.0);

    // Turbine recuperator
    O("HX")->Q("deltaPhot")->set(50.0, "mbar");
    O("HX")->S("hotoption")->set("D");
    O("HX")->Q("deltaPcold")->set(50.0, "mbar");
    O("HX")->S("coldoption")->set("D");
    // O("HX")->S("option")->set("Thot");
    // O("HX")->Q("hotT")->set(324.9144+273.15, "K");
    O("HX")->S("option")->set("UA");
    O("HX")->Q("U")->set(100.0, "W/(m2*K)");
    O("HX")->Q("A")->set(14.2866994561544, "m2");

    O("C")->S("option")->set("P");
    O("C")->Q("P")->set(Pset, "bar");
    O("C")->Q("theta")->set(0.85);
    O("C")->Q("etaM")->set(0.99);
    O("C")->Q("etaE")->set(1.0);

    O("T")->S("option")->set("P");
    O("T")->Q("P")->set(1013.25 + 50.0, "mbar");
    O("T")->Q("theta")->set(0.719017307023255);
    O("T")->Q("etaM")->set(0.99);
    O("T")->Q("etaE")->set(1.0);

    diagnostic(3, "Initializing cut streams");
    // Turbine recuperator (HX) cold outlet, streamB52.17
    O("S09")->Q("T")->set(540.2536 + 273.15, "K");
    O("S09")->Q("P")->set(Pset, "bar");
    O("S09")->S("flowoption")->set("Nx");
    my_cast<stream *>(O("S09"), CURRENT_FUNCTION)->clearcomposition();
    O("S09")->O("Tphase")->Q("ndot")->set(63.69577 / 2.0, "kmol/h");
    O("S09")->O("Tphase")->Q("x", "N2")->set(0.79);
    O("S09")->O("Tphase")->Q("x", "O2")->set(0.21);

    diagnostic(3, "Defining cut streams");
    flowsheetBase::addcut("S09");

    diagnostic(3, "Making selected outputs visible in GUI");
    powerIn.setOutput();
    powerEl.setOutput();
    etaEl.setOutput();

    TIT.setOutput();
    mC.setOutput();
    mT.setOutput();
    maC.setOutput();
    maT.setOutput();
    nr.setOutput();

    O("RX")->Q("T")->setOutput();
    O("RX")->Q("P")->setOutput();
    O("RX")->Q("duty")->setOutput();
    O("RX")->Q("deltaP")->setOutput();
    O("RX")->O("reactions", 0)->Q("conv")->setOutput();
    O("RX")->O("reactions", 0)->Q("z")->setOutput();
    O("RX")->O("reactions", 1)->Q("conv")->setOutput();
    O("RX")->O("reactions", 1)->Q("z")->setOutput();
    O("RX")->O("reactions", 2)->Q("conv")->setOutput();
    O("RX")->O("reactions", 2)->Q("z")->setOutput();
  } // try
  catch (error &e) {
    e.append(CURRENT_FUNCTION);
    throw;
  } // catch
} // microturbo::setup

void microturbo::calculate(int MAXITS, int level) {
  static const int verbosityLocal = 0;
  try {
    diagnostic(2, " entered for " << tag());
    reset_errors();

    AAA::calculate(MAXITS, level);

    powerIn = *O("S02:Tphase")->Q("ndotcomps", "CH4") * LHV_CH4 +
              *O("S02:Tphase")->Q("ndotcomps", "H2") * LHV_H2 +
              *O("S02:Tphase")->Q("ndotcomps", "CO") * LHV_CO;

    powerEl = -*O("C")->Q("We") - *O("T")->Q("We");

    etaEl = powerEl / powerIn;
  } // try
  catch (error &e) {
    e.append(CURRENT_FUNCTION);
    throw;
  } // catch
} // microturbo::calculate

void microturbo::makeuserassembly(std::list<assignment *>::iterator &p) {
  static const int verbosityLocal = -1;
  try {
    diagnostic(2, " entered for " << tag());
    long i(0);

    diagnostic(3, "microturbo flowsheet specific equations");
    // turbine matching
    nr = sqrt((*O("C")->Q("cr") - One) / (beta0 - One));
    mC = mrC0 * nr * *O("S01")->Q("P") / sqrt(*O("S01")->Q("T") / T0);
    // mrT0 = mT0*sqrt(TIT/T0)/PT0;
    mT = mrT0 * nr * *O("S05")->Q("P") / sqrt(*O("S05")->Q("T") / T0);
    maC = *O("S01:Tphase")->Q("mdot");
    maT = *O("S05:Tphase")->Q("mdot");

    MakeAssignment(*Q("T.P"), *Q("T.P") - (*Q("S11.P") - Qdouble(101325.0, "Pa")), "set discharge pressure");

    MakeAssignment(*O("S01:Tphase")->Q("ndot"), mC / *O("S01:Tphase")->Q("AMW"),
                   "Adjust inlet air to balance GTS");
    MakeAssignment(*O("S02:Tphase")->Q("ndot"), *O("S02:Tphase")->Q("ndot") +
                   (TIT - *O("S05")->Q("T")) * *O("S05:Tphase")->Q("cp") * *O("S05:Tphase")->Q("mdot") / LHV_CH4,
                   "Adjust NG to GTS burner to fix TIT");

    // linear
    // Qdouble factor = One - (One - nr)/3.;
    // quadratic
    Qdouble factor = One - (nr - One) * (nr - One) / 1.5;
    // cubic
    // Qdouble factor = One + pow(nr-One, 3.)/0.27;
    MakeAssignment(*Q("C.theta"), 0.85 * factor, "compressor efficiency");
    // MakeAssignment(*Q("T.theta"), 0.719017307023255*pow(factor, 0.1), "turbine efficiency");
    /*
        // Enable these to match vendor data
        MakeAssignment(*Q("HX.A"), *Q("HX.A") * Qdouble(275.+273.15, "K") / *Q("S11.T"),
          "Match vendor claimed exhaust temperature");
        MakeAssignment(*Q("T.theta"), *Q("T.theta") * sqrt(*Q("S05.T") / Qdouble(1144., "K")),
          "Match vendor claimed TIT");
     */
  } // try
  catch (error &e) {
    e.append(CURRENT_FUNCTION);
    throw;
  } // catch
} // microturbo::makeuserassembly

const char *default_type = "microturbo";

// /////////////////////////////////////////////////////////////////////////////
// //////////////////////////// registrar_mymodels /////////////////////////////
// /////////////////////////////////////////////////////////////////////////////

class registrar_mymodels : public registrar {
private:
  static bool done_;
public:
  registrar_mymodels(void);
}; // class registrar_mymodels

registrar_mymodels::registrar_mymodels(void) {
  static const int verbosityLocal = 0;
  diagnostic(2, "Entered");
  if (!done_) {
    diagnostic(3, "Registering");
    model_factory().Register<microturbo>("microturbo");

    done_ = true;
  }
} // registrar_mymodels::registrar_mymodels

bool registrar_mymodels::done_(false);

static registrar_mymodels r_;

#include "Q.cc"