diff --git a/src/libcadet/model/parts/ConvectionDispersionOperator.cpp b/src/libcadet/model/parts/ConvectionDispersionOperator.cpp index 68c98f609..d1ccbea28 100644 --- a/src/libcadet/model/parts/ConvectionDispersionOperator.cpp +++ b/src/libcadet/model/parts/ConvectionDispersionOperator.cpp @@ -103,6 +103,7 @@ bool ConvectionDispersionOperatorBase::configure(UnitOpIdx unitOpIdx, IParameter { readScalarParameterOrArray(_velocity, paramProvider, "VELOCITY", 1); } + _dir = 1; readScalarParameterOrArray(_colDispersion, paramProvider, "COL_DISPERSION", 1); if (paramProvider.exists("COL_DISPERSION_MULTIPLEX")) @@ -185,29 +186,32 @@ bool ConvectionDispersionOperatorBase::configure(UnitOpIdx unitOpIdx, IParameter */ bool ConvectionDispersionOperatorBase::notifyDiscontinuousSectionTransition(double t, unsigned int secIdx) { - double prevVelocity = static_cast(_curVelocity); + // setFlowRates() was called before, so _curVelocity has direction dirOld + const int dirOld = _dir; - // If we don't have cross section area, velocity is given by parameter if (_crossSection <= 0.0) + { + // Use the provided _velocity (direction is also set), only update _dir _curVelocity = getSectionDependentScalar(_velocity, secIdx); + _dir = (_curVelocity >= 0.0) ? 1 : -1; + } else if (!_velocity.empty()) { - if (secIdx > 0) - { - const double dir = static_cast(getSectionDependentScalar(_velocity, secIdx - 1)); - if (dir < 0.0) - prevVelocity *= -1.0; - } + // Use network flow rate but take direction from _velocity + _dir = (getSectionDependentScalar(_velocity, secIdx) >= 0.0) ? 1 : -1; - // We have both cross section area and interstitial flow rate - // _curVelocity has already been set to the network flow rate in setFlowRates() - // the direction of the flow (i.e., sign of _curVelocity) is given by _velocity - const double dir = static_cast(getSectionDependentScalar(_velocity, secIdx)); - if (dir < 0.0) + // _curVelocity has correct magnitude but previous direction, so flip it if necessary + if (dirOld * _dir < 0) _curVelocity *= -1.0; } - return (prevVelocity * static_cast(_curVelocity) < 0.0); + // Remaining case: _velocity is empty and _crossSection <= 0.0 + // _curVelocity is goverend by network flow rate provided in setFlowRates(). + // Direction never changes (always forward, that is, _dir = 1)- + // No action required. + + // Detect change in flow direction + return (dirOld * _dir < 0); } /** @@ -221,7 +225,7 @@ void ConvectionDispersionOperatorBase::setFlowRates(const active& in, const acti { // If we have cross section area, interstitial velocity is given by network flow rates if (_crossSection > 0.0) - _curVelocity = in / (_crossSection * colPorosity); + _curVelocity = _dir * in / (_crossSection * colPorosity); } /** diff --git a/src/libcadet/model/parts/ConvectionDispersionOperator.hpp b/src/libcadet/model/parts/ConvectionDispersionOperator.hpp index f5d17cd28..5970f7292 100644 --- a/src/libcadet/model/parts/ConvectionDispersionOperator.hpp +++ b/src/libcadet/model/parts/ConvectionDispersionOperator.hpp @@ -121,6 +121,7 @@ class ConvectionDispersionOperatorBase std::vector _colDispersion; //!< Column dispersion (may be section dependent) \f$ D_{\text{ax}} \f$ std::vector _velocity; //!< Interstitial velocity (may be section dependent) \f$ u \f$ active _curVelocity; //!< Current interstitial velocity \f$ u \f$ in this time section + int _dir; //!< Current flow direction in this time section ArrayPool _stencilMemory; //!< Provides memory for the stencil double* _wenoDerivatives; //!< Holds derivatives of the WENO scheme diff --git a/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.cpp b/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.cpp index b2bceae15..bd03ac5b4 100644 --- a/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.cpp +++ b/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.cpp @@ -653,7 +653,7 @@ class TwoDimensionalConvectionDispersionOperator::DenseDirectSolver : public Two /** * @brief Creates a TwoDimensionalConvectionDispersionOperator */ -TwoDimensionalConvectionDispersionOperator::TwoDimensionalConvectionDispersionOperator() : _colPorosities(0), _stencilMemory(sizeof(active) * Weno::maxStencilSize()), +TwoDimensionalConvectionDispersionOperator::TwoDimensionalConvectionDispersionOperator() : _colPorosities(0), _dir(0), _stencilMemory(sizeof(active) * Weno::maxStencilSize()), _wenoDerivatives(new double[Weno::maxStencilSize()]), _weno(), _linearSolver(nullptr) { } @@ -850,6 +850,8 @@ bool TwoDimensionalConvectionDispersionOperator::configure(UnitOpIdx unitOpIdx, else registerParam2DArray(parameters, _velocity, [=](bool multi, unsigned int sec, unsigned int compartment) { return makeParamId(hashString("VELOCITY"), unitOpIdx, CompIndep, compartment, BoundStateIndep, ReactionIndep, multi ? sec : SectionIndep); }, _nRad); + _dir = std::vector(_nRad, 1); + _axialDispersionMode = readAndRegisterMultiplexParam(paramProvider, parameters, _axialDispersion, "COL_DISPERSION", _nComp, _nRad, unitOpIdx); _radialDispersionMode = readAndRegisterMultiplexParam(paramProvider, parameters, _radialDispersion, "COL_DISPERSION_RADIAL", _nComp, _nRad, unitOpIdx); @@ -877,25 +879,17 @@ bool TwoDimensionalConvectionDispersionOperator::notifyDiscontinuousSectionTrans { // _curVelocity has already been set to the network flow rate in setFlowRates() // the direction of the flow (i.e., sign of _curVelocity) is given by _velocity + active const* const dirNew = getSectionDependentSlice(_velocity, _nRad, secIdx); - active const* const dirs = getSectionDependentSlice(_velocity, _nRad, secIdx); - if (secIdx > 0) + for (unsigned int i = 0; i < _nRad; ++i) { - active const* const dirsPrev = getSectionDependentSlice(_velocity, _nRad, secIdx - 1); - for (unsigned int i = 0; i < _nRad; ++i) + const int newDir = (dirNew[i] >= 0) ? 1 : -1; + if (_dir[i] * newDir < 0) { - if (dirs[i] * dirsPrev[i] < 0.0) - { - hasChanged = true; - break; - } + hasChanged = true; + _curVelocity[i] *= -1; } - } - - for (unsigned int i = 0; i < _nRad; ++i) - { - if (dirs[i] < 0.0) - _curVelocity[i] *= -1.0; + _dir[i] = newDir; } } @@ -915,7 +909,7 @@ bool TwoDimensionalConvectionDispersionOperator::notifyDiscontinuousSectionTrans */ void TwoDimensionalConvectionDispersionOperator::setFlowRates(int compartment, const active& in, const active& out) CADET_NOEXCEPT { - _curVelocity[compartment] = in / (_crossSections[compartment] * _colPorosities[compartment]); + _curVelocity[compartment] = _dir[compartment] * in / (_crossSections[compartment] * _colPorosities[compartment]); } void TwoDimensionalConvectionDispersionOperator::setFlowRates(active const* in, active const* out) CADET_NOEXCEPT diff --git a/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.hpp b/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.hpp index f3f1f528f..c895c7b23 100644 --- a/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.hpp +++ b/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.hpp @@ -167,6 +167,7 @@ class TwoDimensionalConvectionDispersionOperator MultiplexMode _radialDispersionMode; //!< Multiplex mode of the radial dispersion std::vector _velocity; //!< Interstitial velocity parameter std::vector _curVelocity; //!< Current interstitial velocity \f$ u \f$ + std::vector _dir; //!< Current flow direction bool _singleVelocity; //!< Determines whether only one velocity for all compartments is given ArrayPool _stencilMemory; //!< Provides memory for the stencil