Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix handling of flow direction for dynamic flow rates #101

Merged
merged 7 commits into from
Sep 23, 2021
Merged

Fix handling of flow direction for dynamic flow rates #101

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
c37fac5
Fix handling of flow direction for dynamic flow rates
schmoelder Aug 16, 2021
971aaa5
added fix for 2D case
schmoelder Aug 17, 2021
c2c3481
included review
schmoelder Aug 17, 2021
0ecc576
review part 2
schmoelder Aug 19, 2021
5fd603c
review part 3
schmoelder Aug 19, 2021
df438f5
initialize _dir in constructor
schmoelder Aug 20, 2021
140136d
review part 4
schmoelder Aug 20, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 12 additions & 15 deletions src/libcadet/model/parts/ConvectionDispersionOperator.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"))
Expand Down Expand Up @@ -185,29 +186,25 @@ bool ConvectionDispersionOperatorBase::configure(UnitOpIdx unitOpIdx, IParameter
*/
bool ConvectionDispersionOperatorBase::notifyDiscontinuousSectionTransition(double t, unsigned int secIdx)
{
double prevVelocity = static_cast<double>(_curVelocity);

// If we don't have cross section area, velocity is given by parameter
double dir_old = _dir;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

_dir is an int now:

const int dir_old = _dir;

if (_crossSection <= 0.0)
_curVelocity = getSectionDependentScalar(_velocity, secIdx);
{
double dir_new = static_cast<double>(getSectionDependentScalar(_velocity, secIdx));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We don't need to change this value, so make it const, please.

_dir = (dir_new >= 0);
Copy link
Contributor

@sleweke-bayer sleweke-bayer Aug 19, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This will result in _dir having the value either 0 (test fails) or 1 (test succeeds). We need -1 or 1, so I'd suggest

_dir = (dir_new >= 0.0) ? 1 : -1;

This occurs some more times in the remaining code.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see my mistake. Can you please elaborate on the syntax of your suggestion?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This ? is called a ternary operator. It's similar to Python's

a = 2 if b == 3 else 6

which would translate to

int a = (b == 3) ? 2 : 6;

So it's (condition) ? trueValue : falseValue;.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's what I expected but I really was not sure about the operator. Thanks for explaining!

}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We also need to set _curVelocity here. In setFlowRates(), _curVelocity is only set if _crossSection > 0.0. Thus, we need to set _curVelocity here (inside the curVelocity <= 0.0 branch) in this case. Otherwise, it will be undefined.

Copy link
Contributor Author

@schmoelder schmoelder Aug 19, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think you're right, but I'm also a bit confused. ^^

(inside the curVelocity <= 0.0 branch)

Do you mean inside the _crossSection <= 0.0 branch?

However, a bit up, you also write:

However, the flipping must occur only in the else if branch above. Otherwise, this could happen:
Constant flow rate per section, transition from positive to negative.

With that you simply meant that it was wrong to put the statement that overwrites _dir outside the conditions but it needs to be set in both cases, right?

If so, doesn't this mean, I can combine both branches with a simple OR (||)?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you mean inside the _crossSection <= 0.0 branch?

Yes, sorry. I meant the _crossSection <= 0.0 branch.

With that you simply meant that it was wrong to put the statement that overwrites _dir outside the conditions but it needs to be set in both cases, right?

_dir needs to be set in all cases. But the flip of _curVelocity must only happen in the else if branch:

else if (!_velocity.empty())
{
	// ...

	const double dirNew = static_cast<double>(getSectionDependentScalar(_velocity, secIdx));
	_dir = (dirNew >= 0.0) ? 1 : -1;

	if (dirOld * _dir < 0.0)
		_curVelocity *= -1.0;
}

else if (!_velocity.empty())
{
if (secIdx > 0)
{
const double dir = static_cast<double>(getSectionDependentScalar(_velocity, secIdx - 1));
if (dir < 0.0)
prevVelocity *= -1.0;
}

// 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 *= -1.0;
double dir_new = static_cast<double>(getSectionDependentScalar(_velocity, secIdx));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We don't need to change this value, so make it const, please.

_dir = (dir_new >= 0);
}
if (dir_old * _dir < 0.0)
_curVelocity *= -1.0;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this working or should we do this instead:

_curVelocity = abs(_curVelocity) * _dir;

Copy link
Contributor Author

@schmoelder schmoelder Aug 19, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it works; it only changes the sign of the velocity iff dir_old != _dir; but I don't mind changing it.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're probably right. When a discontinuous section transition occurs, the ModelSystem first calls setFlowRates(). Here, the velocity gets the sign of the previous _dir. Then, notifyDiscontinuousSectionTransition() is called. Here, we reverse the direction if necessary. So I agree, this should be correct.

However, the flipping must occur only in the else if branch above. Otherwise, this could happen:
Constant flow rate per section, transition from positive to negative. _curVelocity is set in the if branch and it is negative as it should be. We detect flow reversal and flip it again. _curVelocity ends up positive, which is wrong.


return (prevVelocity * static_cast<double>(_curVelocity) < 0.0);
return (dir_old * _dir < 0.0);
}

/**
Expand All @@ -221,7 +218,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);
}

/**
Expand Down
1 change: 1 addition & 0 deletions src/libcadet/model/parts/ConvectionDispersionOperator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
26 changes: 10 additions & 16 deletions src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);

std::vector<int> _dir(_nRad, 1);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This creates a local variable that shadows the member variable with the same name.
You probably wanted to do

_dir = std::vector<int>(_nRad, 1);

You should also initialize the vector in the constructor with size 0 to prevent unnecessary memory allocation.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ah yes, I was a bit confused there...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You should also initialize the vector in the constructor with size 0 to prevent unnecessary memory allocation.

Sorry, can you tell me where would that be?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Around line 656 in TwoDimensionalConvectionDispersionOperator::TwoDimensionalConvectionDispersionOperator().
There's the syntax

SomeClass::SomeClass() : _var1(value), _var2(ctorArg1, ctorArg2) { /* ... /* }

The part between : and { is an initializer list. Here, we can initialize (or call the constructor) of all our member variables. But we can also leave out the ones we don't want to initialize here. This happens before the block { } is executed. Note that the variables that do appear have to be in the order of declaration in the class.

My suggestion is to add _dir(0) to this list. This creates a vector of size 0, which should prevent memory allocation.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks; How do you decide what to initialize? E.g. _curVelocity is not initialized there.

Also, should I resize it like other radially dependent parameter vectors are?
https://github.com/modsim/CADET/blob/0ecc57660ae38712f0715a82bdb15ecb1b7ef587/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.cpp#L731

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You don't need to resize() since you replace the whole vector later on (when you assign a new vector like in _dir = std::vector<int>(...).

It would be consistent to also initialize the other vectors with size 0. I was a little lazy, obviously. It doesn't matter much, though. Memory allocation only happens at the beginning of a simulation and we're talking about small sizes. I just wanted to raise awareness for memory stuff. 😄


_axialDispersionMode = readAndRegisterMultiplexParam(paramProvider, parameters, _axialDispersion, "COL_DISPERSION", _nComp, _nRad, unitOpIdx);
_radialDispersionMode = readAndRegisterMultiplexParam(paramProvider, parameters, _radialDispersion, "COL_DISPERSION_RADIAL", _nComp, _nRad, unitOpIdx);

Expand All @@ -872,30 +874,22 @@ bool TwoDimensionalConvectionDispersionOperator::configure(UnitOpIdx unitOpIdx,
bool TwoDimensionalConvectionDispersionOperator::notifyDiscontinuousSectionTransition(double t, unsigned int secIdx)
{
bool hasChanged = false;
std::vector<int> dir_old = _dir;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we can avoid this local (and memory allocating) variable. See below for a suggestion.


if (!_velocity.empty())
{
// _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 dirs = getSectionDependentSlice(_velocity, _nRad, secIdx);
if (secIdx > 0)
{
active const* const dirsPrev = getSectionDependentSlice(_velocity, _nRad, secIdx - 1);
for (unsigned int i = 0; i < _nRad; ++i)
{
if (dirs[i] * dirsPrev[i] < 0.0)
{
hasChanged = true;
break;
}
}
}
active const* const dir_new = getSectionDependentSlice(_velocity, _nRad, secIdx);

for (unsigned int i = 0; i < _nRad; ++i)
{
if (dirs[i] < 0.0)
_dir[i] = (dir_new[i] > 0);
if (_dir[i] * dir_old[i] < 0.0)
{
hasChanged = true;
_curVelocity[i] *= -1.0;
}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doing it this way allows to get rid of the local dir_new variable:

const int newDir = (dir_new[i] >= 0.0) ? 1 : -1;
if (_dir[i] * newDir < 0)
{
	hasChanged = true;
	_curVelocity[i] *= -1.0;
}
_dir[i] = newDir;

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doing it this way allows to get rid of the local dir_new variable:

here you mean dir_old, right?

}
}

Expand All @@ -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
Expand Down
1 change: 1 addition & 0 deletions src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down