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.

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<double>(_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<double>(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<double>(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<double>(_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);
 }
 
 /**

src/libcadet/model/parts/ConvectionDispersionOperator.hpp

@@ -121,6 +121,7 @@ class ConvectionDispersionOperatorBase
 	std::vector<active> _colDispersion; //!< Column dispersion (may be section dependent) \f$ D_{\text{ax}} \f$
 	std::vector<active> _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

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<int>(_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

src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.hpp

@@ -167,6 +167,7 @@ class TwoDimensionalConvectionDispersionOperator
 	MultiplexMode _radialDispersionMode; //!< Multiplex mode of the radial dispersion
 	std::vector<active> _velocity; //!< Interstitial velocity parameter
 	std::vector<active> _curVelocity; //!< Current interstitial velocity \f$ u \f$
+	std::vector<int> _dir; //!< Current flow direction 
 	bool _singleVelocity; //!< Determines whether only one velocity for all compartments is given
 
 	ArrayPool _stencilMemory; //!< Provides memory for the stencil