GreenPolyProducts/cantera
5
0
Fork 0
mirror of https://github.com/Cantera/cantera.git synced 2025-02-25 18:55:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

[oneD] Simplify/improve boundary logic

Browse Source
This commit is contained in:
Ingmar Schoegl 2023-07-24 11:59:09 -06:00 committed by Ray Speth
parent a997a317ce
commit 5152b226cc
2 changed files with 37 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
src/oneD/Boundary1D.cpp
View File

@ -191,6 +191,8 @@ void Inlet1D::eval(size_t jg, double* xg, double* rg,
// The third flow residual is for T, where it is set to T(0). Subtract
// the local temperature to hold the flow T to the inlet T.
rb[c_offset_T] -= m_temp;
} else {
rb[c_offset_T] -= m_flow->T_fixed(0);
}
if (m_flow->isFree()) {
@ -221,6 +223,8 @@ void Inlet1D::eval(size_t jg, double* xg, double* rg,
rb[c_offset_V] -= m_V0;
if (m_flow->doEnergy(m_flow->nPoints() - 1)) {
rb[c_offset_T] -= m_temp; // T
} else {
rb[c_offset_T] -= m_flow->T_fixed(m_flow->nPoints() - 1);
}
rb[c_offset_U] += m_mdot; // u
for (size_t k = 0; k < m_nsp; k++) {
@ -380,6 +384,11 @@ void Outlet1D::eval(size_t jg, double* xg, double* rg, integer* diagg,
for (size_t k = c_offset_Y; k < nc; k++) {
rb[k] = xb[k] - xb[k + nc];
}
if (m_flow_left->doEnergy(0)) {
rb[c_offset_T] = xb[c_offset_T + nc] - xb[c_offset_T]; // zero T gradient
} else {
rb[c_offset_T] = xb[c_offset_T] - m_flow_left->T_fixed(0);
}
}
if (m_flow_left) {
@ -389,8 +398,11 @@ void Outlet1D::eval(size_t jg, double* xg, double* rg, integer* diagg,
double* rb = r - nc;
int* db = diag - nc;
if (m_flow_left->doEnergy(m_flow_left->nPoints()-1)) {
size_t last = m_flow_left->nPoints() - 1;
if (m_flow_left->doEnergy(last)) {
rb[c_offset_T] = xb[c_offset_T] - xb[c_offset_T - nc]; // zero T gradient
} else {
rb[c_offset_T] = xb[c_offset_T] - m_flow_left->T_fixed(last);
}
size_t kSkip = c_offset_Y + m_flow_left->rightExcessSpecies();
for (size_t k = c_offset_Y; k < nc; k++) {
@ -467,13 +479,11 @@ void OutletRes1D::eval(size_t jg, double* xg, double* rg,
double* xb = x;
double* rb = r;
// this seems wrong...
// zero Lambda
rb[c_offset_U] = xb[c_offset_L];
if (m_flow_right->doEnergy(0)) {
// zero gradient for T
rb[c_offset_T] = xb[c_offset_T] - xb[c_offset_T + nc];
} else {
rb[c_offset_T] = xb[c_offset_T] - m_flow_left->T_fixed(0);
}
// specified mass fractions
@ -488,11 +498,11 @@ void OutletRes1D::eval(size_t jg, double* xg, double* rg,
double* rb = r - nc;
int* db = diag - nc;
if (!m_flow_left->usesLambda()) {
rb[c_offset_L] = xb[c_offset_L]; // zero Lambda
}
if (m_flow_left->doEnergy(m_flow_left->nPoints()-1)) {
rb[c_offset_T] = xb[c_offset_T] - m_temp; // zero dT/dz
size_t last = m_flow_left->nPoints() - 1;
if (m_flow_left->doEnergy(last)) {
rb[c_offset_T] = xb[c_offset_T] - xb[c_offset_T - nc]; // zero T gradient
} else {
rb[c_offset_T] = xb[c_offset_T] - m_flow_left->T_fixed(last);
}
size_t kSkip = m_flow_left->rightExcessSpecies();
for (size_t k = c_offset_Y; k < nc; k++) {

47
src/oneD/StFlow.cpp
View File

@ -480,16 +480,12 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag,
// a result, these residual equations will force the solution
// variables to the values for the boundary object
rsd[index(c_offset_V,0)] = V(x,0);
if (doEnergy(0)) {
rsd[index(c_offset_T,0)] = T(x,0);
} else {
rsd[index(c_offset_T,0)] = T(x,0) - T_fixed(0);
}
if (m_isFree) {
rsd[index(c_offset_L, 0)] = lambda(x, 0);
} else {
// used to set mass flow rate (both axisymmetric and unstrained)
rsd[index(c_offset_T,0)] = T(x,0);
if (m_usesLambda) {
rsd[index(c_offset_L, 0)] = -rho_u(x, 0);
} else {
rsd[index(c_offset_L, 0)] = lambda(x, 0);
diag[index(c_offset_L, 0)] = 0;
}
// The default boundary condition for species is zero flux. However,
@ -574,10 +570,10 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag,
diag[index(c_offset_T, j)] = 0;
}
if (m_isFree) {
rsd[index(c_offset_L, j)] = lambda(x, j);
} else {
if (m_usesLambda) {
rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j - 1);
} else {
rsd[index(c_offset_L, j)] = lambda(x, j);
}
diag[index(c_offset_L, j)] = 0;
}
@ -1015,8 +1011,6 @@ void StFlow::evalRightBoundary(double* x, double* rsd, int* diag, double rdt)
rsd[index(c_offset_V,j)] = V(x,j);
double sum = 0.0;
rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
diag[index(c_offset_L, j)] = 0;
// set residual of poisson's equ to zero
rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)];
for (size_t k = 0; k < m_nsp; k++) {
@ -1025,24 +1019,16 @@ void StFlow::evalRightBoundary(double* x, double* rsd, int* diag, double rdt)
}
rsd[index(c_offset_Y + rightExcessSpecies(), j)] = 1.0 - sum;
diag[index(c_offset_Y + rightExcessSpecies(), j)] = 0;
if (m_isFree) {
rsd[index(c_offset_U,j)] = rho_u(x,j) - rho_u(x,j-1);
rsd[index(c_offset_T,j)] = T(x,j) - T(x,j-1);
} else if (m_usesLambda) {
rsd[index(c_offset_U,j)] = rho_u(x,j);
if (m_do_energy[j]) {
rsd[index(c_offset_T,j)] = T(x,j);
} else {
rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j);
}
if (m_usesLambda) {
rsd[index(c_offset_U, j)] = rho_u(x, j);
rsd[index(c_offset_L, j)] = lambda(x, j);
} else {
rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j - 1);
if (m_do_energy[j]) {
rsd[index(c_offset_T, j)] = T(x, j);
} else {
rsd[index(c_offset_T, j)] = T(x, j) - T_fixed(j);
}
rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j-1);
rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j-1);
}
diag[index(c_offset_L, j)] = 0;
rsd[index(c_offset_T, j)] = T(x, j);
}
void StFlow::evalContinuity(size_t j, double* x, double* rsd, int* diag, double rdt)
@ -1079,6 +1065,7 @@ void StFlow::evalContinuity(size_t j, double* x, double* rsd, int* diag, double
- (density(j+1)*V(x,j+1) + density(j)*V(x,j));
}
} else {
// unstrained
rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j - 1);
}
}