initial import
This commit is contained in:
Dave Goodwin
parent
c7f6cd18be
commit
1d742bf984
1 changed files with 785 additions and 0 deletions
785
Cantera/src/oneD/boundaries1D.cpp
Normal file
785
Cantera/src/oneD/boundaries1D.cpp
Normal file
|
|
@ -0,0 +1,785 @@
|
|||
/**
|
||||
* @file boundaries1D.cpp
|
||||
*/
|
||||
|
||||
/*
|
||||
* $Author$
|
||||
* $Revision$
|
||||
* $Date$
|
||||
*/
|
||||
|
||||
// Copyright 2002-3 California Institute of Technology
|
||||
|
||||
|
||||
|
||||
|
||||
#include "Inlet1D.h"
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
Bdry1D::Bdry1D() : Domain1D(1, 1, 0.0),
|
||||
m_flow_left(0), m_flow_right(0),
|
||||
m_ilr(0), m_left_nv(0), m_right_nv(0),
|
||||
m_left_loc(0), m_right_loc(0),
|
||||
m_left_points(0), m_nv(0),
|
||||
m_left_nsp(0), m_right_nsp(0),
|
||||
m_start_left(0), m_start_right(0),
|
||||
m_phase_left(0), m_phase_right(0), m_mdot(0.0) {
|
||||
m_type = cConnectorType;
|
||||
}
|
||||
|
||||
|
||||
void Bdry1D::
|
||||
_init(int n) {
|
||||
if (m_index < 0) {
|
||||
throw CanteraError("Bdry1D",
|
||||
"install in container before calling init.");
|
||||
}
|
||||
resize(n,1);
|
||||
|
||||
m_left_nsp = 0;
|
||||
m_right_nsp = 0;
|
||||
|
||||
// check for left and right flow objects
|
||||
if (m_index > 0) {
|
||||
Domain1D& r = container().domain(m_index-1);
|
||||
if (r.domainType() == cFlowType) {
|
||||
m_flow_left = (StFlow*)&r;
|
||||
m_left_nv = m_flow_left->nComponents();
|
||||
m_left_points = m_flow_left->nPoints();
|
||||
m_left_loc = container().start(m_index-1);
|
||||
m_left_nsp = m_left_nv - 4;
|
||||
m_phase_left = &m_flow_left->phase();
|
||||
}
|
||||
else
|
||||
throw CanteraError("Bdry1D::init",
|
||||
"Boundary domains can only be "
|
||||
"connected on the left to flow domains, not type "+int2str(r.domainType())
|
||||
+ " domains.");
|
||||
}
|
||||
|
||||
if (m_index < container().nDomains() - 1) {
|
||||
Domain1D& r = container().domain(m_index+1);
|
||||
if (r.domainType() == cFlowType) {
|
||||
m_flow_right = (StFlow*)&r;
|
||||
m_right_nv = m_flow_right->nComponents();
|
||||
m_right_loc = container().start(m_index+1);
|
||||
m_right_nsp = m_right_nv - 4;
|
||||
m_phase_right = &m_flow_right->phase();
|
||||
}
|
||||
else
|
||||
throw CanteraError("Bdry1D::init",
|
||||
"Boundary domains can only be "
|
||||
"connected on the right to flow domains, not type "+int2str(r.domainType())
|
||||
+ " domains.");
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
//----------------------------------------------------------
|
||||
//
|
||||
// Inlet1D methods
|
||||
//
|
||||
//----------------------------------------------------------
|
||||
|
||||
|
||||
void Inlet1D::
|
||||
setMoleFractions(string xin) {
|
||||
m_xstr = xin;
|
||||
if (m_flow) {
|
||||
m_flow->phase().setMoleFractionsByName(xin);
|
||||
m_flow->phase().getMassFractions(m_yin.begin());
|
||||
needJacUpdate();
|
||||
}
|
||||
}
|
||||
|
||||
void Inlet1D::
|
||||
setMoleFractions(doublereal* xin) {
|
||||
if (m_flow) {
|
||||
m_flow->phase().setMoleFractions(xin);
|
||||
m_flow->phase().getMassFractions(m_yin.begin());
|
||||
needJacUpdate();
|
||||
}
|
||||
}
|
||||
|
||||
string Inlet1D::
|
||||
componentName(int n) const {
|
||||
switch (n) {
|
||||
case 0: return "mdot"; break;
|
||||
case 1: return "temperature"; break;
|
||||
default: return "unknown";
|
||||
}
|
||||
}
|
||||
|
||||
void Inlet1D::
|
||||
init() {
|
||||
|
||||
_init(2);
|
||||
|
||||
// set bounds (mdot, T)
|
||||
const doublereal lower[2] = {-1.0e5, 200.0};
|
||||
const doublereal upper[2] = {1.0e5, 1.e5};
|
||||
setBounds(2, lower, 2, upper);
|
||||
|
||||
// set tolerances
|
||||
vector_fp rtol(2, 1e-4);
|
||||
vector_fp atol(2, 1.e-5);
|
||||
setTolerances(2, rtol.begin(), 2, atol.begin());
|
||||
|
||||
// if a flow domain is present on the left, then this must be
|
||||
// a right inlet. Note that an inlet object can only be a
|
||||
// terminal object.
|
||||
if (m_flow_left) {
|
||||
m_ilr = RightInlet;
|
||||
m_flow = m_flow_left;
|
||||
}
|
||||
else if (m_flow_right) {
|
||||
m_ilr = LeftInlet;
|
||||
m_flow = m_flow_right;
|
||||
}
|
||||
else {
|
||||
throw CanteraError("Inlet1D::init","no flow!");
|
||||
}
|
||||
|
||||
// components = u, V, T, lambda, + mass fractions
|
||||
m_nsp = m_flow->nComponents() - 4;
|
||||
m_yin.resize(m_nsp, 0.0);
|
||||
if (m_xstr != "")
|
||||
setMoleFractions(m_xstr);
|
||||
else
|
||||
m_yin[0] = 1.0;
|
||||
}
|
||||
|
||||
|
||||
void Inlet1D::
|
||||
eval(int jg, doublereal* xg, doublereal* rg,
|
||||
integer* diagg, doublereal rdt) {
|
||||
int k;
|
||||
if (jg >= 0 && (jg < firstPoint() - 2 || jg > lastPoint() + 2)) return;
|
||||
|
||||
// start of local part of global arrays
|
||||
doublereal* x = xg + loc();
|
||||
doublereal* r = rg + loc();
|
||||
integer* diag = diagg + loc();
|
||||
doublereal *xb, *rb;
|
||||
|
||||
// residual equations for the two local variables
|
||||
r[0] = m_mdot - x[0];
|
||||
r[1] = m_temp - x[1];
|
||||
|
||||
// both are algebraic constraints
|
||||
diag[0] = 0;
|
||||
diag[1] = 0;
|
||||
|
||||
// if it is a left inlet, then the flow solution vector
|
||||
// starts 2 to the right in the global solution vector
|
||||
if (m_ilr == LeftInlet) {
|
||||
xb = x + 2;
|
||||
rb = r + 2;
|
||||
|
||||
rb[2] = xb[2] - x[1]; // T
|
||||
|
||||
// spreading rate. Flow domain sets this to V(0),
|
||||
// so for finite spreading rate subtract m_V0.
|
||||
rb[1] -= m_V0;
|
||||
|
||||
rb[3] += x[0]; // lambda
|
||||
for (k = 0; k < m_nsp; k++) {
|
||||
if (m_flow->doSpecies(k)) {
|
||||
rb[4+k] += x[0]*m_yin[k];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// right inlet.
|
||||
else {
|
||||
int boffset = m_flow->nComponents();
|
||||
xb = x - boffset;
|
||||
rb = r - boffset;
|
||||
rb[1] -= m_V0;
|
||||
rb[2] = xb[2] - x[1]; // T
|
||||
rb[0] += x[0]; // u
|
||||
for (k = 1; k < m_nsp; k++) {
|
||||
if (m_flow->doSpecies(k)) {
|
||||
rb[4+k] += x[0]*(-xb[4+k] + m_yin[k]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void Inlet1D::
|
||||
save(XML_Node& o, doublereal* soln) {
|
||||
doublereal* s = soln + loc();
|
||||
//XML_Node& inlt = o.addChild("inlet");
|
||||
XML_Node& inlt = o.addChild("domain");
|
||||
inlt.addAttribute("id",id());
|
||||
inlt.addAttribute("points",1);
|
||||
inlt.addAttribute("type","inlet");
|
||||
inlt.addAttribute("components",nComponents());
|
||||
for (int k = 0; k < nComponents(); k++) {
|
||||
ctml::addFloat(inlt, componentName(k), s[k], "", "",0.0, 1.0);
|
||||
}
|
||||
}
|
||||
|
||||
void Inlet1D::
|
||||
restore(XML_Node& dom, doublereal* soln) {
|
||||
map<string, double> x;
|
||||
getFloats(dom, x);
|
||||
soln[0] = x["mdot"];
|
||||
soln[1] = x["temperature"];
|
||||
resize(2,1);
|
||||
}
|
||||
|
||||
|
||||
//--------------------------------------------------
|
||||
// Symm1D
|
||||
//--------------------------------------------------
|
||||
|
||||
string Symm1D::componentName(int n) const {
|
||||
switch (n) {
|
||||
case 0: return "dummy"; break;
|
||||
default: return "<unknown>";
|
||||
}
|
||||
}
|
||||
|
||||
void Symm1D::
|
||||
init() { _init(1);
|
||||
// set bounds (T)
|
||||
const doublereal lower = -1.0;
|
||||
const doublereal upper = 1.0;
|
||||
setBounds(1, &lower, 1, &upper);
|
||||
|
||||
// set tolerances
|
||||
const doublereal rtol = 1e-4;
|
||||
const doublereal atol = 1.e-4;
|
||||
setTolerances(1, &rtol, 1, &atol);
|
||||
}
|
||||
|
||||
void Symm1D::
|
||||
eval(int jg, doublereal* xg, doublereal* rg,
|
||||
integer* diagg, doublereal rdt) {
|
||||
if (jg >= 0 && (jg < firstPoint() - 2 || jg > lastPoint() + 2)) return;
|
||||
|
||||
// start of local part of global arrays
|
||||
doublereal* x = xg + loc();
|
||||
doublereal* r = rg + loc();
|
||||
integer* diag = diagg + loc();
|
||||
doublereal *xb, *rb;
|
||||
integer *db;
|
||||
|
||||
r[0] = x[0];
|
||||
diag[0] = 0;
|
||||
int nc;
|
||||
|
||||
if (m_flow_right) {
|
||||
nc = m_flow_right->nComponents();
|
||||
xb = x + 1;
|
||||
rb = r + 1;
|
||||
db = diag + 1;
|
||||
db[1] = 0;
|
||||
db[2] = 0;
|
||||
rb[1] = xb[1] - xb[1 + nc]; // zero dV/dz
|
||||
rb[2] = xb[2] - xb[2 + nc]; // zero dT/dz
|
||||
}
|
||||
|
||||
if (m_flow_left) {
|
||||
nc = m_flow_left->nComponents();
|
||||
xb = x - nc;
|
||||
rb = r - nc;
|
||||
db = diag - nc;
|
||||
db[1] = 0;
|
||||
db[2] = 0;
|
||||
rb[1] = xb[1] - xb[1 - nc]; // zero dV/dz
|
||||
rb[2] = xb[2] - xb[2 - nc]; // zero dT/dz
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Symm1D::
|
||||
save(XML_Node& o, doublereal* soln) {
|
||||
XML_Node& symm = o.addChild("domain");
|
||||
symm.addAttribute("id",id());
|
||||
symm.addAttribute("points",1);
|
||||
symm.addAttribute("type","outlet");
|
||||
symm.addAttribute("components",nComponents());
|
||||
}
|
||||
|
||||
void Symm1D::
|
||||
restore(XML_Node& dom, doublereal* soln) {
|
||||
resize(1,1);
|
||||
}
|
||||
|
||||
|
||||
//--------------------------------------------------
|
||||
// Outlet1D
|
||||
//--------------------------------------------------
|
||||
|
||||
string Outlet1D::componentName(int n) const {
|
||||
switch (n) {
|
||||
case 0: return "dummy"; break;
|
||||
default: return "<unknown>";
|
||||
}
|
||||
}
|
||||
|
||||
void Outlet1D::
|
||||
init() {
|
||||
_init(1);
|
||||
// set bounds (T)
|
||||
const doublereal lower = -1.0;
|
||||
const doublereal upper = 1.0;
|
||||
setBounds(1, &lower, 1, &upper);
|
||||
|
||||
// set tolerances
|
||||
const doublereal rtol = 1e-4;
|
||||
const doublereal atol = 1.e-4;
|
||||
setTolerances(1, &rtol, 1, &atol);
|
||||
}
|
||||
|
||||
|
||||
void Outlet1D::
|
||||
eval(int jg, doublereal* xg, doublereal* rg,
|
||||
integer* diagg, doublereal rdt) {
|
||||
if (jg >= 0 && (jg < firstPoint() - 2 || jg > lastPoint() + 2)) return;
|
||||
|
||||
// start of local part of global arrays
|
||||
doublereal* x = xg + loc();
|
||||
doublereal* r = rg + loc();
|
||||
integer* diag = diagg + loc();
|
||||
doublereal *xb, *rb;
|
||||
integer *db;
|
||||
|
||||
r[0] = x[0];
|
||||
diag[0] = 0;
|
||||
int nc, n;
|
||||
|
||||
if (m_flow_right) {
|
||||
nc = m_flow_right->nComponents();
|
||||
xb = x + 1;
|
||||
rb = r + 1;
|
||||
db = diag + 1;
|
||||
rb[0] = xb[3];
|
||||
rb[2] = xb[2] - xb[2 + nc];
|
||||
}
|
||||
|
||||
if (m_flow_left) {
|
||||
nc = m_flow_left->nComponents();
|
||||
xb = x - nc;
|
||||
rb = r - nc;
|
||||
db = diag - nc;
|
||||
rb[0] = xb[3];
|
||||
rb[2] = xb[2] - xb[2 - nc];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Outlet1D::
|
||||
save(XML_Node& o, doublereal* soln) {
|
||||
XML_Node& outlt = o.addChild("domain");
|
||||
outlt.addAttribute("id",id());
|
||||
outlt.addAttribute("points",1);
|
||||
outlt.addAttribute("type","outlet");
|
||||
outlt.addAttribute("components",nComponents());
|
||||
}
|
||||
|
||||
void Outlet1D::
|
||||
restore(XML_Node& dom, doublereal* soln) {
|
||||
resize(1,1);
|
||||
}
|
||||
|
||||
|
||||
//-----------------------------------------------------------
|
||||
//
|
||||
// Surf1D
|
||||
//
|
||||
//-----------------------------------------------------------
|
||||
|
||||
|
||||
|
||||
string Surf1D::componentName(int n) const {
|
||||
switch (n) {
|
||||
case 0: return "temperature"; break;
|
||||
default: return "<unknown>";
|
||||
}
|
||||
}
|
||||
|
||||
void Surf1D::
|
||||
init() {
|
||||
_init(1);
|
||||
// set bounds (T)
|
||||
const doublereal lower = 200.0;
|
||||
const doublereal upper = 1.e5;
|
||||
setBounds(1, &lower, 1, &upper);
|
||||
|
||||
// set tolerances
|
||||
const doublereal rtol = 1e-4;
|
||||
const doublereal atol = 1.e-4;
|
||||
setTolerances(1, &rtol, 1, &atol);
|
||||
}
|
||||
|
||||
|
||||
void Surf1D::
|
||||
eval(int jg, doublereal* xg, doublereal* rg,
|
||||
integer* diagg, doublereal rdt) {
|
||||
if (jg >= 0 && (jg < firstPoint() - 2 || jg > lastPoint() + 2)) return;
|
||||
|
||||
// start of local part of global arrays
|
||||
doublereal* x = xg + loc();
|
||||
doublereal* r = rg + loc();
|
||||
integer* diag = diagg + loc();
|
||||
doublereal *xb, *rb;
|
||||
|
||||
r[0] = x[0] - m_temp;
|
||||
diag[0] = 0;
|
||||
int nc;
|
||||
|
||||
if (m_flow_right) {
|
||||
rb = r + 1;
|
||||
xb = x + 1;
|
||||
rb[2] = xb[2] - x[0]; // specified T
|
||||
}
|
||||
|
||||
if (m_flow_left) {
|
||||
nc = m_flow_left->nComponents();
|
||||
rb = r - nc;
|
||||
xb = x - nc;
|
||||
rb[2] = xb[2] - x[0]; // specified T
|
||||
}
|
||||
}
|
||||
|
||||
void Surf1D::
|
||||
save(XML_Node& o, doublereal* soln) {
|
||||
doublereal* s = soln + loc();
|
||||
//XML_Node& inlt = o.addChild("inlet");
|
||||
XML_Node& inlt = o.addChild("domain");
|
||||
inlt.addAttribute("id",id());
|
||||
inlt.addAttribute("points",1);
|
||||
inlt.addAttribute("type","surface");
|
||||
inlt.addAttribute("components",nComponents());
|
||||
for (int k = 0; k < nComponents(); k++) {
|
||||
ctml::addFloat(inlt, componentName(k), s[k], "", "",0.0, 1.0);
|
||||
}
|
||||
}
|
||||
|
||||
void Surf1D::
|
||||
restore(XML_Node& dom, doublereal* soln) {
|
||||
map<string, double> x;
|
||||
getFloats(dom, x);
|
||||
soln[0] = x["temperature"];
|
||||
resize(1,1);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
/////////////////////////////////////////////////////////////
|
||||
//
|
||||
// surf1D
|
||||
//
|
||||
////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
// string ChemSurf1D::
|
||||
// componentName(int n) const {
|
||||
// /// @todo why dummy?
|
||||
// switch (n) {
|
||||
// case 0: return "dummy"; break;
|
||||
// case 1: return "temperature"; break;
|
||||
// default: return "<unknown>";
|
||||
// }
|
||||
// }
|
||||
|
||||
// // Set the kinetics manager for the surface.
|
||||
// void ChemSurf1D::
|
||||
// setKinetics(InterfaceKinetics* kin) {
|
||||
// m_kin = kin;
|
||||
// int np = kin->nPhases();
|
||||
// m_sphase = 0;
|
||||
// for (int n = 0; n < np; n++) {
|
||||
// if (kin->phase(n).eosType() == cSurf) {
|
||||
// m_sphase = (SurfPhase*)&m_kin->phase(n);
|
||||
// m_nsurf = n;
|
||||
// }
|
||||
// else {
|
||||
// m_bulk.push_back(&kin->phase(n));
|
||||
// m_nbulk.push_back(n);
|
||||
// }
|
||||
// }
|
||||
// if (!m_sphase)
|
||||
// throw CanteraError("setKinetics","no surface phase defined");
|
||||
|
||||
// m_nsp = m_sphase->nSpecies();
|
||||
// resize(m_nsp,1);
|
||||
// if (m_bulk.size() == 1) {
|
||||
// m_bulk.push_back(0);
|
||||
// }
|
||||
// }
|
||||
|
||||
// void ChemSurf1D::
|
||||
// init() {
|
||||
// cout << "ChemSurf1D::init" << endl;
|
||||
// if (m_index < 0) {
|
||||
// throw CanteraError("Surf1D",
|
||||
// "install in container before calling init.");
|
||||
// }
|
||||
// resize(m_nsp,1);
|
||||
// m_mult.resize(m_nsp, 1.0);
|
||||
// m_do_surf_species.resize(m_nsp, true);
|
||||
// if (!m_sphase) m_do_surf_species[0] = false;
|
||||
// m_fixed_cov.resize(m_nsp, 1.0/m_nsp);
|
||||
|
||||
// // set bounds
|
||||
// vector_fp lower(m_nsp, -1.e-3);
|
||||
// vector_fp upper(m_nsp, 1.0);
|
||||
// setBounds(m_nsp, lower.begin(), m_nsp, upper.begin());
|
||||
|
||||
// // set tolerances
|
||||
// vector_fp rtol(m_nsp, 1e-4);
|
||||
// vector_fp atol(m_nsp, 1.e-10);
|
||||
// setTolerances(m_nsp, rtol.begin(), m_nsp, atol.begin());
|
||||
|
||||
// m_left_nsp = 0;
|
||||
// m_right_nsp = 0;
|
||||
|
||||
// // check for left and right flow objects
|
||||
// if (m_index > 0) {
|
||||
// Domain1D& r = container().domain(m_index-1);
|
||||
// if (r.domainType() == cFlowType) {
|
||||
// m_flow_left = (StFlow*)&r;
|
||||
// m_left_nv = m_flow_left->nComponents();
|
||||
// m_left_points = m_flow_left->nPoints();
|
||||
// m_left_loc = container().start(m_index-1);
|
||||
// m_left_nsp = m_left_nv - 4;
|
||||
// m_phase_left = &m_flow_left->phase();
|
||||
// m_molwt_left = m_phase_left->molecularWeights().begin();
|
||||
// if (m_phase_left == m_bulk[0])
|
||||
// m_start_left = m_kin->start(m_nbulk[0]);
|
||||
// else if (m_phase_left == m_bulk[1])
|
||||
// m_start_left = m_kin->start(m_nbulk[1]);
|
||||
// else
|
||||
// throw CanteraError("ChemSurf1D::init",
|
||||
// "left gas does not match one in surface mechanism");
|
||||
// }
|
||||
// else
|
||||
// throw CanteraError("ChemSurf1D::init",
|
||||
// "Surface domains can only be "
|
||||
// "connected to flow domains.");
|
||||
// }
|
||||
|
||||
// if (m_index < container().nDomains() - 1) {
|
||||
// Domain1D& r = container().domain(m_index+1);
|
||||
// if (r.domainType() == cFlowType) {
|
||||
// m_flow_right = (StFlow*)&r;
|
||||
// m_right_nv = m_flow_right->nComponents();
|
||||
// m_right_loc = container().start(m_index+1);
|
||||
// m_right_nsp = m_right_nv - 4;
|
||||
// m_phase_right = &m_flow_right->phase();
|
||||
// m_molwt_right = m_phase_right->molecularWeights().begin();
|
||||
// if (m_phase_right == m_bulk[0])
|
||||
// m_start_right = m_kin->start(m_nbulk[0]);
|
||||
// else if (m_phase_right == m_bulk[1])
|
||||
// m_start_right = m_kin->start(m_nbulk[1]);
|
||||
// else
|
||||
// throw CanteraError("ChemSurf1D::init",
|
||||
// "right gas does not match one in surface mechanism");
|
||||
// }
|
||||
// else
|
||||
// throw CanteraError("ChemSurf1D::init",
|
||||
// "Surface domains can only be "
|
||||
// "connected to flow domains.");
|
||||
// }
|
||||
// m_work.resize(m_kin->nTotalSpecies());
|
||||
// }
|
||||
|
||||
|
||||
// void ChemSurf1D::eval(int jg, doublereal* xg, doublereal* rg,
|
||||
// integer* diagg, doublereal rdt) {
|
||||
// int k;
|
||||
|
||||
// // if computing a Jacobian (jg > 0), and the global point is
|
||||
// // outside the points the surface can influence, then skip
|
||||
// // evaluating the residual
|
||||
// if (jg >= 0 && (jg < firstPoint() - 2
|
||||
// || jg > lastPoint() + 2)) return;
|
||||
|
||||
// // start of local part of global arrays
|
||||
// doublereal* x = xg + loc();
|
||||
// doublereal* r = rg + loc();
|
||||
// integer* diag = diagg + loc();
|
||||
|
||||
// // set the coverages
|
||||
// doublereal sum = 0.0;
|
||||
// for (k = 0; k < m_nsp; k++) {
|
||||
// m_work[k] = x[k];
|
||||
// sum += x[k];
|
||||
// }
|
||||
// m_sphase->setCoverages(m_work.begin());
|
||||
|
||||
// // set the left gas state to the adjacent point
|
||||
|
||||
// int leftloc = 0, rightloc = 0;
|
||||
// int pnt = 0;
|
||||
|
||||
// if (m_flow_left) {
|
||||
// leftloc = m_flow_left->loc();
|
||||
// pnt = m_flow_left->nPoints() - 1;
|
||||
// m_flow_left->setGas(xg + leftloc, pnt);
|
||||
// }
|
||||
|
||||
// if (m_flow_right) {
|
||||
// rightloc = m_flow_right->loc();
|
||||
// m_flow_right->setGas(xg + rightloc, 0);
|
||||
// }
|
||||
|
||||
// m_kin->getNetProductionRates(m_work.begin());
|
||||
// doublereal rs0 = 1.0/m_sphase->siteDensity();
|
||||
|
||||
// scale(m_work.begin(), m_work.end(), m_work.begin(), m_mult[0]);
|
||||
|
||||
// bool enabled = true;
|
||||
// int ioffset = m_kin->start(m_nsurf); // m_left_nsp + m_right_nsp;
|
||||
// doublereal maxx = -1.0;
|
||||
// int imx = -1;
|
||||
// for (k = 0; k < m_nsp; k++) {
|
||||
// r[k] = m_work[k + ioffset] * m_sphase->size(k) * rs0;
|
||||
// r[k] -= rdt*(x[k] - prevSoln(k,0));
|
||||
// diag[k] = 1;
|
||||
// if (x[k] > maxx) {
|
||||
// maxx = x[k];
|
||||
// imx = k;
|
||||
// }
|
||||
// if (!m_do_surf_species[k]) {
|
||||
// r[k] = x[k] - m_fixed_cov[k];
|
||||
// diag[k] = 0;
|
||||
// enabled = false;
|
||||
// }
|
||||
// }
|
||||
// if (enabled) {
|
||||
// r[imx] = 1.0 - sum;
|
||||
// diag[imx] = 0;
|
||||
// }
|
||||
|
||||
// // gas-phase residuals
|
||||
// doublereal rho;
|
||||
// if (m_flow_left) {
|
||||
// rho = m_phase_left->density();
|
||||
// // doublereal rdz = 2.0/
|
||||
// // (m_flow_left->z(m_left_points-1) -
|
||||
// // m_flow_left->z(m_left_points - 2));
|
||||
|
||||
// for (k = 0; k < m_left_nsp; k++)
|
||||
// m_work[k + m_start_left] *= m_molwt_left[k];
|
||||
|
||||
// int ileft = loc() - m_left_nv;
|
||||
|
||||
// // if the energy equation is enabled at this point,
|
||||
// // set the gas temperature to the surface temperature
|
||||
// if (m_flow_left->doEnergy(pnt)) {
|
||||
// rg[ileft + 2] = xg[ileft + 2] - m_sphase->temperature();
|
||||
// }
|
||||
|
||||
// for (k = 1; k < m_left_nsp; k++) {
|
||||
// if (enabled && m_flow_left->doSpecies(k)) {
|
||||
// rg[ileft + 4 + k] += m_work[k + m_start_left];
|
||||
// //+= rdz*m_work[k + m_sp_left]/rho;
|
||||
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
|
||||
// if (m_flow_right) {
|
||||
// for (k = 0; k < m_right_nsp; k++)
|
||||
// m_work[k + m_start_right] *= m_molwt_right[k];
|
||||
|
||||
// int iright = loc() + m_nsp;
|
||||
// rg[iright + 2] -= m_sphase->temperature();
|
||||
// //r[iright + 3] = x[iright];
|
||||
// for (k = 0; k < m_right_nsp; k++) {
|
||||
// rg[iright + 4 + k] -= m_work[k + m_start_right];
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
|
||||
// void ChemSurf1D::
|
||||
// save(XML_Node& o, doublereal* soln) {
|
||||
// doublereal* s = soln + loc();
|
||||
// XML_Node& srf = o.addChild("surface");
|
||||
// for (int k = 0; k < m_nsp; k++) {
|
||||
// ctml::addFloat(srf, componentName(k), s[k], "", "coverage",
|
||||
// 0.0, 1.0);
|
||||
// }
|
||||
// }
|
||||
|
||||
|
||||
// /////////////////////////////////////////////////////////////
|
||||
// //
|
||||
// // surf1D
|
||||
// //
|
||||
// ////////////////////////////////////////////////////////////
|
||||
|
||||
// string ChemSurf1D::
|
||||
// componentName(int n) const {
|
||||
// /// @todo why dummy?
|
||||
// switch (n) {
|
||||
// case 0: return "temperature"; break;
|
||||
// default: return "<unknown>";
|
||||
// }
|
||||
// }
|
||||
|
||||
// void ChemSurf1D::
|
||||
// init() {
|
||||
// if (m_index < 0) {
|
||||
// throw CanteraError("Surf1D",
|
||||
// "install in container before calling init.");
|
||||
// }
|
||||
// if (m_index > 0) m_left_flow = true;
|
||||
// resize(1,1);
|
||||
|
||||
// // set bounds
|
||||
// vector_fp lower(1, 200.0);
|
||||
// vector_fp upper(1, 10000.0);
|
||||
// setBounds(1, lower.begin(), 1, upper.begin());
|
||||
|
||||
// // set tolerances
|
||||
// vector_fp rtol(1, 1e-4);
|
||||
// vector_fp atol(1, 1.e-5);
|
||||
// setTolerances(1, rtol.begin(), 1, atol.begin());
|
||||
// }
|
||||
|
||||
|
||||
// void ChemSurf1D::eval(int jg, doublereal* xg, doublereal* rg,
|
||||
// integer* diagg, doublereal rdt) {
|
||||
// int k;
|
||||
|
||||
// // if computing a Jacobian (jg > 0), and the global point is
|
||||
// // outside the points the surface can influence, then skip
|
||||
// // evaluating the residual
|
||||
// if (jg >= 0 && (jg < firstPoint() - 2
|
||||
// || jg > lastPoint() + 2)) return;
|
||||
|
||||
// // start of local part of global arrays
|
||||
// doublereal* x = xg + loc();
|
||||
// doublereal* r = rg + loc();
|
||||
// integer* diag = diagg + loc();
|
||||
|
||||
// // set the left gas state to the adjacent point
|
||||
|
||||
// // gas-phase residuals
|
||||
// doublereal rho;
|
||||
// }
|
||||
|
||||
// void ChemSurf1D::
|
||||
// save(XML_Node& o, doublereal* soln) {
|
||||
// doublereal* s = soln + loc();
|
||||
// XML_Node& srf = o.addChild("surface");
|
||||
// for (int k = 0; k < m_nsp; k++) {
|
||||
// ctml::addFloat(srf, componentName(k), s[k], "", "coverage",
|
||||
// 0.0, 1.0);
|
||||
// }
|
||||
// }
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
Loading…
Add table
Reference in a new issue