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

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

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

// turn on to enable Quality Control mode
// #define QC

static const int verbositySubSystem = 0;

#include "graph.h"
#include "fuelcells.h"
#include "purecomps.h"

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

// /////////////////////////////////////////////////////////////////////////////
// //////////////////////////////// mcfcsystem /////////////////////////////////
// /////////////////////////////////////////////////////////////////////////////

/// AFC mcfc fuel cell system with concentrated parameter model of cell
class mcfcsystem : public modelBase, public flowsheet<zero_zero> {
private:
  static std::string type_;
  void maketables(long);
public:
  mcfcsystem(const std::string &t, const std::string &d) : persistent_item_interface(-1), persistent_item(-1, t, d), modelBaseInterface(-1), modelBase(t, d) {
    maketables(-1);
  };
  mcfcsystem(long cid) : persistent_item_interface(cid), persistent_item(cid, "", ""), modelBaseInterface(cid), modelBase(cid), flowsheet<zero_zero>(cid) {
    maketables(cid);
    readtables(cid);
  }

  Qdouble maxPower; ///< Inlet fuel LHV
  Qdouble etaStack; ///< stack electric yield
  Qdouble etaEl;    ///< stack + blower electric yield

  void calculate(int MAXITS = 100, int level = 0);
  void makeuserassembly(std::list<assignment *>::iterator &p);
  virtual void setup(void);

  virtual const std::string &type(void) const {
    return type_;
  }
  int smiter_precalc(void) { return 4; }
  bool supports_simultaneous(void) { return true; };
  int maxit(void) { return 200; }  
}; // mcfcsystem

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

void mcfcsystem::calculate(int MAXITS, int level) {
  diagnostic(2, "Entered for " << fulltag());

  flowsheet<zero_zero>::calculate(MAXITS, level);
  maxPower = *O("ANIN:Tphase")->Q("ndotcomps", "CH4") * Qdouble(802600.0, "kJ/kmol") +
             *O("ANIN:Tphase")->Q("ndotcomps", "H2") * Qdouble(241000.0, "kJ/kmol") +
             *O("ANIN:Tphase")->Q("ndotcomps", "CO") * Qdouble(283000.0, "kJ/kmol");
  etaStack = *O("FC:multiReactions", 0)->Q("Power") / maxPower;
  etaEl = (*O("FC:multiReactions", 0)->Q("Power") + *O("BLOWER")->Q("duty")) / maxPower;
} // mcfcsystem::calculate

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

    // Non embedded objects
    MAKEQDOUBLE(maxPower, 0.0, "W", "Inlet fuel LHV");
    MAKEQDOUBLE(etaStack, 0.0, "", "stack electric yield");
    MAKEQDOUBLE(etaEl, 0.0, "", "stack and blower electric yield");

    if (cid == -1) {
      diagnostic(2, "Define unit operations");
      makeVertex<mixer>("MIX1", "Mixes inlet air with recompressed reformer outlet", cid);
      makeVertex<mixer>("MIX2", "Mixes anode and cathode outlets", cid);
      makeVertex<compressor>("BLOWER", "Compresses reformer outlet", cid);
      makeVertex<hx>("REGHEX", "Regenerative heat exchanger", cid);
      makeVertex<splitter>("SPLIT", "Cathode outlet splitter", cid, boost::assign::map_list_of<std::string, long>("nStreams", 2));
      makeVertex<hx>("REFORMER", "Reformer", cid,
                     boost::assign::map_list_of<std::string, long>("nReactions", 2)
                            ("nMultiReactions", 0),
                     boost::assign::map_list_of<std::string, std::string>("embeddedTypeReactions[0]", "reactionCH4_Ref_Eq")
                            ("embeddedTypeReactions[1]", "reactionWGS_Eq"));
      makeVertex<genflashN1>("CB", "Catalytic burner", cid,
                             boost::assign::map_list_of<std::string, long>("nReactions", 3),
                             boost::assign::map_list_of<std::string, std::string>("embeddedTypeReactions[0]", "reactionCH4_TC")
                                    ("embeddedTypeReactions[1]", "reactionCO_TC")
                                    ("embeddedTypeReactions[2]", "reactionH2_TC"));
      makeVertex<multihx>("FC", "MCFC stack", cid,
                          boost::assign::map_list_of<std::string, long>("nReactions", 1)
                                 ("nMultiReactions", 1)
                                 ("nStreams", 2),
                          boost::assign::map_list_of<std::string, std::string>("embeddedTypeReactions[0]", "reactionWGS_Eq")
                                 ("embeddedTypeMultiReactions[0]", "multiReaction_MCFC"));

      O("FC")->setdot("width=\"1.0\",height=\"1.0\",shapefile = \"FC.png\",shape=plaintext,color=lightgrey, fixedsize=true");

      diagnostic(2, "Define stream objects and connect");
      makeEdge<stream_V>("AIRIN", "Inlet fresh air", "source", "out", "MIX1", "in", cid);
      makeEdge<stream_V>("FUELIN", "Humidified fuel", "source", "out", "REGHEX", "coldin", cid);
      makeEdge<stream_V>("REFIN", "Preheated humidified fuel", "REGHEX", "coldout", "REFORMER", "in2", cid);

      makeEdge<stream_V>("ANIN", "Reformed cooled fuel to anode", "REGHEX", "hotout", "FC", "in1", cid);
      makeEdge<stream_V>("ANOUT", "Anode outlet", "FC", "out1", "MIX2", "in", cid);
      makeEdge<stream_V>("CATIN", "Cathode inlet", "MIX1", "out", "FC", "in2", cid);
      makeEdge<stream_V>("CATOUT", "Cathode outlet", "FC", "out2", "SPLIT", "in", cid);

      makeEdge<stream_V>("MIXIN", "Recompressed flue gas", "BLOWER", "out", "MIX1", "in", cid);
      makeEdge<stream_V>("REFOUT", "Reformed fuel", "REFORMER", "out2", "REGHEX", "hotin", cid);
      makeEdge<stream_V>("SPLITOUT", "Purged cathode outlet", "SPLIT", "out1", "sink", "in", cid);
      makeEdge<stream_V>("RECYCLE", "Recycled cathode outlet", "SPLIT", "out2", "MIX2", "in", cid);
      makeEdge<stream_V>("RECYCLE1", "Recycled cathode and anode outlet to CB", "MIX2", "out", "CB", "in", cid);
      makeEdge<stream_V>("RECYCLE2", "CB outlet", "CB", "out", "REFORMER", "in1", cid);
      makeEdge<stream_V>("RECYCLE3", "Recycle", "REFORMER", "out1", "BLOWER", "in", cid);
    }
  } // try
  catch (error &e) {
    e.append(CURRENT_FUNCTION);
    throw;
  } // catch
} // mcfcsystem::maketables

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

    diagnostic(3, "Calling flowsheet::setup to initialize any embedded flowsheet");
    flowsheet<zero_zero>::setup();

    diagnostic(3, "Setting input variables");
    O("FUELIN")->Q("T")->set(190.5 + 273.15, "K");
    O("FUELIN")->Q("P")->set(3.6, "atm");
    O("FUELIN")->S("flowoption")->set("Nx");
    O("FUELIN")->O("Tphase")->Q("ndot")->set(19.751, "kmol/h");
    O("FUELIN")->O("Tphase")->Q("x", "H2O")->set(0.778);
    O("FUELIN")->O("Tphase")->Q("x", "N2")->set(0.0);
    O("FUELIN")->O("Tphase")->Q("x", "O2")->set(0.0);
    O("FUELIN")->O("Tphase")->Q("x", "CH4")->set(0.222);
    O("FUELIN")->O("Tphase")->Q("x", "CO")->set(0.0);
    O("FUELIN")->O("Tphase")->Q("x", "CO2")->set(0.0);
    O("FUELIN")->O("Tphase")->Q("x", "H2")->set(0.0);

    O("AIRIN")->Q("T")->set(202.7 + 273.15, "K");
    O("AIRIN")->Q("P")->set(3.6, "atm");
    O("AIRIN")->S("flowoption")->set("Nx");
    O("AIRIN")->O("Tphase")->Q("ndot")->set(95.319, "kmol/h");
    O("AIRIN")->O("Tphase")->Q("x", "H2O")->set(0.00);
    O("AIRIN")->O("Tphase")->Q("x", "N2")->set(0.79);
    O("AIRIN")->O("Tphase")->Q("x", "O2")->set(0.21);
    for (int j = 3; j < NCOMPONENTS; ++j)
      O("AIRIN")->O("Tphase")->Q("x", j)->set(0.0);
    O("CB")->Q("deltaP")->set(40.0, "mbar");
    O("CB")->Q("duty")->set(0.0, "W");
    O("CB")->S("option")->set("DH");
    // total combustion of methane
    O("CB")->O("reactions", 0)->Q("z")->set(1.0);
    // total combustion of CO
    O("CB")->O("reactions", 1)->Q("z")->set(1.0);
    // total combustion of H2
    O("CB")->O("reactions", 2)->Q("z")->set(1.0);

    O("REGHEX")->Q("U")->set(100.0, "W/(m2*K)");
    O("REGHEX")->Q("A")->set(1.0, "m2");

    O("SPLIT")->Q("outSplit", 0)->set(0.205);

    // if hx
    O("REFORMER")->Q("deltaPhot")->set(30.0, "mbar");
    O("REFORMER")->S("hotoption")->set("D");
    O("REFORMER")->Q("deltaPcold")->set(0.0, "Pa");
    O("REFORMER")->S("coldoption")->set("D");
    O("REFORMER")->Q("U")->set(850.0, "W/(m2*K)");
    O("REFORMER")->Q("A")->set(7.0, "m2");
    O("REFORMER")->I("reactionSides", 0)->set_val(1); // reactionCH4_Ref_Eq
    O("REFORMER")->I("reactionSides", 1)->set_val(1); // reactionWGS_Eq
    // reforming of methane, reaction conversion estimate
    O("REFORMER")->O("reactions", 0)->Q("z")->set(0.8);
    O("REFORMER")->O("reactions", 0)->Q("deltaT")->set(0.0, "K");
    // water gas shift reaction, reaction conversion estimate
    O("REFORMER")->O("reactions", 1)->Q("z")->set(0.2);

    //  O("FC")->Q("A")->set(2.0*0.711025*300, "m2");
    O("FC")->S("MTDoption")->set("CSTR");
    O("FC")->Q("A")->set(100.0, "m2");
    O("FC")->Q("netDuty")->set(0.0, "kW"); // heat losses
    // anodic side
    O("FC")->Q("deltaP", 0)->set(4.0, "mbar");
    O("FC")->S("pressureOption", 0)->set("D");
    O("FC")->I("reactionSides", 0)->set_val(0);
    // cathodic side
    O("FC")->Q("deltaP", 1)->set(23.0, "mbar");
    O("FC")->S("pressureOption", 1)->set("D");
    // electrochemical reaction
    O("FC:multiReactions", 0)->S("kinetic")->set("literature");
    O("FC:multiReactions", 0)->Q("H")->set(2.0 * 0.595, "m");
    O("FC:multiReactions", 0)->Q("W")->set(1.195, "m");
    O("FC:multiReactions", 0)->Q("Voffset")->set(-0.01663, "V");
    O("FC:multiReactions", 0)->Q("z")->set(0.724412);
    O("FC:multiReactions", 0)->Q("nCells")->set(300.0);
    // heat transfer coefficients
    O("FC")->Q("U", 0)->set(390.0, "W/(m2*K)");
    O("FC")->Q("U", 1)->set(220.0, "W/(m2*K)");

    O("BLOWER")->Q("P")->set(3.6, "atm");
    O("BLOWER")->Q("deltas")->set(0.0, "W/K");
    O("BLOWER")->S("option")->set("P");

    diagnostic(3, "Initializing cut streams");
    O("REFOUT")->Q("T")->set(754.0 + 273.15, "K");
    O("REFOUT")->Q("P")->set(3.6e5, "Pa");
    O("REFOUT")->S("flowoption")->set("Nx");
    O("REFOUT")->O("Tphase")->Q("ndot")->set(25.4, "kmol/h");
    O("REFOUT")->O("Tphase")->Q("x", "H2O")->set(0.33);
    O("REFOUT")->O("Tphase")->Q("x", "N2")->set(0.0);
    O("REFOUT")->O("Tphase")->Q("x", "O2")->set(0.0);
    O("REFOUT")->O("Tphase")->Q("x", "CH4")->set(0.0001);
    O("REFOUT")->O("Tphase")->Q("x", "CO")->set(0.10);
    O("REFOUT")->O("Tphase")->Q("x", "CO2")->set(0.07);
    O("REFOUT")->O("Tphase")->Q("x", "H2")->set(0.50);

    O("MIXIN")->Q("T")->set(693.0 + 273.15, "K");
    O("MIXIN")->Q("P")->set(3.6e5, "Pa");
    O("MIXIN")->S("flowoption")->set("Nx");
    O("MIXIN")->O("Tphase")->Q("ndot")->set(400.0, "kmol/h");
    O("MIXIN")->O("Tphase")->Q("x", "H2O")->set(0.24);
    O("MIXIN")->O("Tphase")->Q("x", "N2")->set(0.60);
    O("MIXIN")->O("Tphase")->Q("x", "O2")->set(0.08);
    O("MIXIN")->O("Tphase")->Q("x", "CH4")->set(0.0);
    O("MIXIN")->O("Tphase")->Q("x", "CO")->set(0.00);
    O("MIXIN")->O("Tphase")->Q("x", "CO2")->set(0.08);
    O("MIXIN")->O("Tphase")->Q("x", "H2")->set(0.0);

    diagnostic(3, "Defining cut streams");
    addcut("MIXIN");
    addcut("REFOUT");

    diagnostic(3, "Making selected outputs visible in GUI");
    O("CB")->Q("T")->setOutput();
    O("CB")->Q("P")->setOutput();
    O("CB")->Q("duty")->setOutput();
    O("CB")->Q("deltaP")->setOutput();
    O("CB")->O("reactions", 0)->Q("conv")->setOutput();
    O("CB")->O("reactions", 0)->Q("z")->setOutput();
    O("CB")->O("reactions", 1)->Q("conv")->setOutput();
    O("CB")->O("reactions", 1)->Q("z")->setOutput();
    O("CB")->O("reactions", 2)->Q("conv")->setOutput();
    O("CB")->O("reactions", 2)->Q("z")->setOutput();

    O("REFORMER")->Q("coldT")->setOutput();
    O("REFORMER")->Q("hotT")->setOutput();
    O("REFORMER")->Q("coldP")->setOutput();
    O("REFORMER")->Q("hotP")->setOutput();
    O("REFORMER")->Q("dutyhot")->setOutput();
    O("REFORMER")->Q("deltaPcold")->setOutput();
    O("REFORMER")->Q("deltaPhot")->setOutput();

    O("BLOWER")->Q("T")->setOutput();
    O("BLOWER")->Q("P")->setOutput();
    O("BLOWER")->Q("duty")->setOutput();
    O("BLOWER")->Q("deltas")->setOutput();
    O("REGHEX")->Q("A")->setOutput();

    O("REGHEX")->Q("coldT")->setOutput();
    O("REGHEX")->Q("hotT")->setOutput();
    O("REGHEX")->Q("coldP")->setOutput();
    O("REGHEX")->Q("hotP")->setOutput();
    O("REGHEX")->Q("deltaPcold")->setOutput();
    O("REGHEX")->Q("deltaPhot")->setOutput();
    O("REGHEX")->Q("dutyhot")->setOutput();

    O("FC")->Q("Ts")->setOutput();
    O("FC")->Q("T", 0)->setOutput();
    O("FC")->Q("T", 1)->setOutput();
    O("FC")->Q("P", 0)->setOutput();
    O("FC")->Q("P", 1)->setOutput();
    O("FC")->O("multiReactions", 0)->Q("Power")->setOutput();
    O("FC")->O("multiReactions", 0)->Q("I")->setOutput();
    O("FC")->O("multiReactions", 0)->Q("J")->setOutput();
    O("FC")->O("multiReactions", 0)->Q("V")->setOutput();
    O("FC")->O("multiReactions", 0)->Q("conv")->setOutput();
    O("FC")->O("multiReactions", 0)->Q("z")->setOutput();
    O("FC")->O("reactions", 0)->Q("conv")->setOutput();
    O("FC")->O("reactions", 0)->Q("z")->setOutput();
  } // try
  catch (error &e) {
    e.append(CURRENT_FUNCTION);
    throw;
  } // catch
} // mcfcsystem::setup

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

    long i(0);

    diagnostic(3, "MCFC flowsheet specific equations");
    MakeAssignment(*O("FC:multiReactions", 0)->Q("z"),
                   *O("FC:multiReactions", 0)->Q("z") *
                   (0.75 / (1.0 - (*O("ANOUT:Tphase")->Q("ndotcomps", "H2") + *O("ANOUT:Tphase")->Q("ndotcomps", "CO")) / (*O("ANIN:Tphase")->Q("ndotcomps", "H2") + *O("ANIN:Tphase")->Q("ndotcomps", "CO")))),
                   "Set fuel utilization");
    /*
       diagnostic(0, "Voffset = " << *O("FC:multiReactions", 0)->Q("Voffset") << " I = " << *O("FC:multiReactions", 0)->Q("I") << " Voffset new = " << (*O("FC:multiReactions", 0)->Q("Voffset"))*(Qdouble(2200.0, "A")/(*O("FC:multiReactions", 0)->Q("I"))));
       MakeAssignment(*O("FC:multiReactions", 0)->Q("Voffset"),
        (*O("FC:multiReactions", 0)->Q("Voffset"))*pow(Qdouble(2200.0, "A")/(*O("FC:multiReactions", 0)->Q("I")), 2.0),
        "Set current");
       MakeAssignment(*O("FC:multiReactions", 0)->Q("Voffset"),
        (*O("FC:multiReactions", 0)->Q("Voffset"))*pow((*O("FC:multiReactions", 0)->Q("I"))/Qdouble(2200.0, "A"), 2.0),
        "Set current");
     */
  } // try
  catch (error &e) {
    e.append(CURRENT_FUNCTION);
    throw;
  } // catch
} // mcfcsystem::makeuserassembly

const char *default_type = "mcfcsystem";

// /////////////////////////////////////////////////////////////////////////////
// //////////////////////////// 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<mcfcsystem>("mcfcsystem");

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

bool registrar_mymodels::done_(false);

static registrar_mymodels r_;

#include "Q.cc"