Skip to content

Commit

Permalink
[oneD] Simplify/improve boundary logic
Browse files Browse the repository at this point in the history
  • Loading branch information
@ischoegl @speth
ischoegl authored and speth committed Aug 13, 2023
1 parent a997a31 commit 5152b22
Show file tree
Hide file tree
Showing 2 changed files with 36 additions and 39 deletions.
30 changes: 20 additions & 10 deletions src/oneD/Boundary1D.cpp
View file
Expand Up @@ -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()) {
Expand Down Expand Up @@ -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++) {
Expand Down Expand Up @@ -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) {
Expand All @@ -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++) {
Expand Down Expand Up @@ -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
Expand All @@ -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++) {
Expand Down
45 changes: 16 additions & 29 deletions src/oneD/StFlow.cpp
View file
Expand Up @@ -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);
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_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_L, 0)] = -rho_u(x, 0);
diag[index(c_offset_L, 0)] = 0;
}

// The default boundary condition for species is zero flux. However,
Expand Down Expand Up @@ -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;
}
Expand Down Expand Up @@ -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++) {
Expand All @@ -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)
Expand Down Expand Up @@ -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);
}
}
Expand Down

0 comments on commit 5152b22

Please sign in to comment.