-
Notifications
You must be signed in to change notification settings - Fork 26
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.
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
Conversation
Here is a draft for the 2D case: bool TwoDimensionalConvectionDispersionOperator::notifyDiscontinuousSectionTransition(double t, unsigned int secIdx)
{
bool hasChanged = false;
std::vector<double> _dir_old = _dir;
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
_dir = getSectionDependentSlice(_velocity, _nRad, secIdx);
for (unsigned int i = 0; i < _nRad; ++i)
{
if (_dir[i] * _dir_old[i] < 0.0)
{
hasChanged = true;
_curVelocity[i] *= -1.0;
}
}
}
// Change the sparsity pattern if necessary
if ((secIdx == 0) || hasChanged)
setSparsityPattern();
return hasChanged || (secIdx == 0);
} However, I still have two questions.
|
// If we don't have cross section area, velocity is given by parameter | ||
double _dir_old = _dir; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Don't prefix local variables with underscores.
if (_crossSection <= 0.0) | ||
_curVelocity = getSectionDependentScalar(_velocity, secIdx); | ||
_dir = static_cast<double>(getSectionDependentScalar(_velocity, secIdx)); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What happens if _velocity
is not just one of {-1.0, 1.0}
? I suggest to normalize (i.e., only assign -1.0
or 1.0
to _dir
) just to be on the safe side.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Good point! I will simply check the sign of the velocity that we read (also according to the documentation). Should I also change dir
to an int type?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Doesn't make a real difference. So it's your choice.
const double dir = static_cast<double>(getSectionDependentScalar(_velocity, secIdx)); | ||
if (dir < 0.0) | ||
_dir = static_cast<double>(getSectionDependentScalar(_velocity, secIdx)); | ||
if (_dir_old * static_cast<double>(_dir) < 0.0) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The static_cast
is not necessary as _dir
is already of type double
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Actually, if I don't cast, I get the following error:
193 | double dir_new = getSectionDependentScalar(_velocity, secIdx);
| ~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~
| |
| const sfad::Fwd<double>
/home/jo/software/cadet/CADET/src/libcadet/model/parts/ConvectionDispersionOperator.cpp:201:45: error: cannot convert ‘const sfad::Fwd<double>’ to ‘double’ in initialization
201 | double dir_new = getSectionDependentScalar(_velocity, secIdx);
| ~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~
| |
| const sfad::Fwd<double>
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sorry for the confusion:
- The cast is required here:
_dir = static_cast<double>(getSectionDependentScalar(_velocity, secIdx));
- It is not necessary here:
if (_dir_old * static_cast<double>(_dir) < 0.0)
Hm, that's more difficult as we're dealing with active* ptr = getSectionDependentSlice();
std::vector<double>(ptr, ptr + countOfItems); But since you may want to normalize the values first before putting them into the vector, you will need a loop anyway: active* ptr = getSectionDependentSlice();
for (int i = 0; i < _dir.size(); ++i)
{
if (ptr[i] >= 0.0)
vec[i] = 1.0;
else
vec[i] = -1.0;
} |
36938ee
to
c2c3481
Compare
@@ -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); |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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...
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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. 😄
@@ -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; |
There was a problem hiding this comment.
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.
_dir[i] = (dir_new[i] > 0); | ||
if (_dir[i] * dir_old[i] < 0.0) | ||
{ | ||
hasChanged = true; | ||
_curVelocity[i] *= -1.0; | ||
} |
There was a problem hiding this comment.
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;
There was a problem hiding this comment.
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?
} | ||
if (dir_old * _dir < 0.0) | ||
_curVelocity *= -1.0; |
There was a problem hiding this comment.
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;
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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.
_curVelocity = getSectionDependentScalar(_velocity, secIdx); | ||
{ | ||
double dir_new = static_cast<double>(getSectionDependentScalar(_velocity, secIdx)); | ||
_dir = (dir_new >= 0); |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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;
.
There was a problem hiding this comment.
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!
if (_crossSection <= 0.0) | ||
_curVelocity = getSectionDependentScalar(_velocity, secIdx); | ||
{ | ||
double dir_new = static_cast<double>(getSectionDependentScalar(_velocity, secIdx)); |
There was a problem hiding this comment.
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.
// If we don't have cross section area, velocity is given by parameter | ||
double dir_old = _dir; |
There was a problem hiding this comment.
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;
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)); |
There was a problem hiding this comment.
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.
A general remark: You're probably used to do a lot of Python programming. There, using lower case and underscores is common. In this codebase, we're doing |
Got it! Thanks for all your suggestions! |
} | ||
if (dir_old * _dir < 0.0) | ||
_curVelocity *= -1.0; |
There was a problem hiding this comment.
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 (dirOld * _dir < 0.0); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is a comparison int
vs double
. The result dirOld * _dir
will be converted to double
for the comparison. We can simplify this by comparing against 0
instead of 0.0
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I did that and just because I wanted to play around a bit more, I also tried the following:
const int dirNew = static_cast<int>(getSectionDependentScalar(_velocity, secIdx));
but I got this error:
/home/jo/software/cadet/CADET/src/libcadet/model/parts/ConvectionDispersionOperator.cpp:193:83: error: invalid static_cast from type ‘const sfad::Fwd<double>’ to type ‘int’
193 | const int dirNew = static_cast<int>(getSectionDependentScalar(_velocity, secIdx));
Why is that?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
getSectionDependentScalar()
returns an object of type active
. This object only implements an explicit conversion to double
. This is why you can convert this custom type to double
but not to anything else.
The conversion operator is defined here. Note that real_t
is double
in our case.
{ | ||
const double dirNew = static_cast<double>(getSectionDependentScalar(_velocity, secIdx)); | ||
_dir = (dirNew >= 0.0) ? 1 : -1; | ||
} |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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 (||
)?
There was a problem hiding this comment.
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;
}
@@ -729,6 +729,7 @@ bool TwoDimensionalConvectionDispersionOperator::configureModelDiscretization(IP | |||
// _radialCentroids.resize(nRad); | |||
_crossSections.resize(nRad); | |||
_curVelocity.resize(nRad); | |||
_dir.resize(nRad); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is not necessary as you overwrite the full vector later on. Doing
_dir = std::vector<int>(_nRad, 1);
means to
- throw away the current vector
_dir
(deallocate memory if it has capacity> 0
), - construct a new
std::vector<int>
that contains_nRad
elements all with value1
, and - assign the new vector to
_dir
.
When you do _dir.resize(nRad)
, you're telling the vector that you want to reserve memory for nRad
elements. The vector is still empty, but it doesn't need to allocate memory if you put in (up to) nRad
elements. So you're allocating memory here and throw it away later when you do _dir = std::vector<int>(...)
. We can avoid this by not resize()
ing the vector
.
Writing _dir(0)
in the constructor tells the vector
that you want it to hold 0
items initially. Hopefully, it does not allocate memory and sets its capacity to 0
also.
We call resize()
on the other vectors (e.g., _curVelocity
) since we are not replacing them with new vectors, but using indexed access. That is, we need the vectors to be non-empty.
bool ConvectionDispersionOperatorBase::notifyDiscontinuousSectionTransition(double t, unsigned int secIdx) | ||
{ | ||
double prevVelocity = static_cast<double>(_curVelocity); | ||
|
||
const int dirOld = _dir; | ||
// If we don't have cross section area, velocity is given by parameter | ||
if (_crossSection <= 0.0) | ||
_curVelocity = getSectionDependentScalar(_velocity, secIdx); | ||
else if (!_velocity.empty()) | ||
// 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 | ||
if (_crossSection <= 0.0 || !_velocity.empty()) | ||
{ | ||
if (secIdx > 0) | ||
{ | ||
const double dir = static_cast<double>(getSectionDependentScalar(_velocity, secIdx - 1)); | ||
if (dir < 0.0) | ||
prevVelocity *= -1.0; | ||
} | ||
const double dirNew = static_cast<double>(getSectionDependentScalar(_velocity, secIdx)); | ||
_dir = (dirNew >= 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 *= -1.0; | ||
if (dirOld * _dir < 0) | ||
_curVelocity *= -1; | ||
} | ||
|
||
return (prevVelocity * static_cast<double>(_curVelocity) < 0.0); | ||
return (dirOld * _dir < 0); | ||
} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Okay, this got really convoluted. Let's have a look at the original logic:
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
if (_crossSection <= 0.0)
_curVelocity = getSectionDependentScalar(_velocity, secIdx);
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;
}
return (prevVelocity * static_cast<double>(_curVelocity) < 0.0);
}
void ConvectionDispersionOperatorBase::setFlowRates(const active& in, const active& out, const active& colPorosity) CADET_NOEXCEPT
{
// If we have cross section area, interstitial velocity is given by network flow rates
if (_crossSection > 0.0)
_curVelocity = in / (_crossSection * colPorosity);
}
We have three cases:
_velocity
is non-empty and_crossSection <= 0.0
:_curVelocity
(including direction) is goverend by_velocity
._velocity
is non-empty and_crossSection > 0.0
: The magnitude of_curVelocity
is provided by the network insetFlowRates()
and its direction is given by the sign of_velocity
._velocity
is empty and_crossSection > 0.0
: Magnitude and direction are provided by the network insetFlowRates()
.
First, setFlowRates()
is called, then notifyDiscontinuousSectionTransition()
. In case of dynamic flow rates, we have subsequent calls to setFlowRates()
. These subsequent calls were resetting the direction of _curVelocity
, which caused the bug.
In order to preserve the direction of the flow in subsequent calls, you've introduced the _dir
variable and replaced setFlowRates()
:
void ConvectionDispersionOperatorBase::setFlowRates(const active& in, const active& out, const active& colPorosity) CADET_NOEXCEPT
{
// If we have cross section area, interstitial velocity is given by network flow rates
if (_crossSection > 0.0)
_curVelocity = _dir * in / (_crossSection * colPorosity);
}
Now we need to handle the three cases via notifyDiscontinuousSectionTransition()
:
bool ConvectionDispersionOperatorBase::notifyDiscontinuousSectionTransition(double t, unsigned int secIdx)
{
// setFlowRates() was called before, so _curVelocity has direction dirOld
const int dirOld = _dir;
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())
{
// Use network flow rate but take direction from _velocity
_dir = (getSectionDependentScalar(_velocity, secIdx) >= 0.0) ? 1 : -1;
// _curVelocity has correct magnitude but previous direction, so flip it if necessary
if (dirOld * _dir < 0)
_curVelocity *= -1.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);
}
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ah, I get it now! I did not fully understand how flow rates work if no cross section area is provided.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Looks good! We need some final testing before merging.
Thanks for your guidance and patience! :) |
Fixes #96
Still needs to be implemented for the TwoDimensionalConvectionDispersionOperator.