From 2e0af622095ba17628ddc810a993bbaec3ad4315 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Schm=C3=B6lder?= <20299934+schmoelder@users.noreply.github.com> Date: Thu, 23 Sep 2021 10:49:16 +0200 Subject: [PATCH] Fix handling of flow direction for dynamic flow rates If dynamic flow rates were enabled, the column velocity (which is used to determine the flow direction if the cross section area is also provided) was not taken into account. `_curVelocity` was correct in the call to `residual()` after `notifyDiscontinuousSectionTransition()`. However, it was overwritten in subsequent calls to `residual()` which happened almost instantly (in terms of time) if dynamic flow rates are configured. --- .../parts/ConvectionDispersionOperator.cpp | 34 +++++++++++-------- .../parts/ConvectionDispersionOperator.hpp | 1 + ...imensionalConvectionDispersionOperator.cpp | 28 ++++++--------- ...imensionalConvectionDispersionOperator.hpp | 1 + 4 files changed, 32 insertions(+), 32 deletions(-) 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