Skip to content

Commit

Permalink
fix residual
Browse files Browse the repository at this point in the history
  • Loading branch information
@schmoelder
schmoelder committed Jul 21, 2021
1 parent 4f69811 commit 74059fd
Showing 1 changed file with 15 additions and 22 deletions.
37 changes: 15 additions & 22 deletions src/libcadet/model/parts/TractorConvectionDispersionOperator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1118,16 +1118,14 @@ int TractorConvectionDispersionOperator::residualImpl(double t, unsigned int sec

for (unsigned int rad_orig = 0; rad_orig < _nRad; ++rad_orig)
{

const unsigned int offsetToRadOrigBlock = rad_orig * _nComp;
const unsigned int offsetColRadOrigBlock = offsetColBlock + offsetToRadOrigBlock;
ResidualType* const resColRadOrigBlock = resColBlock + offsetToRadOrigBlock;
StateType const* const yColRadOrigBlock = yColBlock + offsetToRadOrigBlock;


for (unsigned int rad_dest = 0; rad_dest < rad_orig; ++rad_dest)
/* for (unsigned int rad_dest = 0; rad_dest < (_nRad - rad_orig); ++rad_dest) */
for (unsigned int rad_dest = 0; rad_dest < _nRad; ++rad_dest)
{

const unsigned int offsetToRadDestBlock = rad_dest * _nComp;
const unsigned int offsetColRadDestBlock = offsetColBlock + offsetToRadDestBlock;
ResidualType* const resColRadDestBlock = resColBlock + offsetToRadDestBlock;
Expand All @@ -1143,25 +1141,21 @@ int TractorConvectionDispersionOperator::residualImpl(double t, unsigned int sec
ResidualType* const resCur_dest = resColRadDestBlock + comp;

// NOTE: Tractor implementation follows
// TODO: Look at actual flux equations
const ParamType exchange_orig_dest_comp = static_cast<ParamType>(_exchangeMatrix[rad_orig * _nRad * _nComp + rad_dest * _nComp + comp]);
*resCur_orig += exchange_orig_dest_comp * (yCur_orig[0] - yCur_dest[0]);
*resCur_dest -= exchange_orig_dest_comp * (yCur_orig[0] - yCur_dest[0]);

// TODO: We can handle forward and backward rates in one step and only iterate over the upper part of the matrix
const ParamType exchange_dest_orig_comp = static_cast<ParamType>(_exchangeMatrix[rad_dest * _nRad * _nComp + rad_orig * _nComp + comp]);
*resCur_orig += exchange_dest_orig_comp * (yCur_dest[0] - yCur_orig[0]);
*resCur_dest += exchange_dest_orig_comp * (yCur_orig[0] - yCur_dest[0]);

// TODO: @Sam, please check this
if (wantJac)
{
_jacC.centered(offsetCur_orig, 0) += static_cast<double>(exchange_orig_dest_comp);
_jacC.centered(offsetCur_orig, offsetCur_dest) -= static_cast<double>(exchange_orig_dest_comp);
_jacC.centered(offsetCur_dest, 0) -= static_cast<double>(exchange_orig_dest_comp);
_jacC.centered(offsetCur_dest, offsetCur_orig) += static_cast<double>(exchange_orig_dest_comp);
}
if (cadet_likely(exchange_orig_dest_comp > 0.0))
{
*resCur_orig += exchange_orig_dest_comp * yCur_orig[0];
*resCur_dest -= exchange_orig_dest_comp * yCur_orig[0];

// TODO: @Sam, please check this
if (wantJac)
{
_jacC.centered(offsetCur_orig, 0) += static_cast<double>(exchange_orig_dest_comp);
_jacC.centered(offsetCur_orig, offsetCur_dest) -= static_cast<double>(exchange_orig_dest_comp);
_jacC.centered(offsetCur_dest, 0) -= static_cast<double>(exchange_orig_dest_comp);
_jacC.centered(offsetCur_dest, offsetCur_orig) += static_cast<double>(exchange_orig_dest_comp);
}
}
}

}
Expand Down Expand Up @@ -1212,7 +1206,6 @@ void TractorConvectionDispersionOperator::setSparsityPattern()

for (unsigned int comp = 0; comp < _nComp; ++comp)
{
// TODO: Only add pattern if element in exchange matrix is actually non-zero
const ParamType exchange_orig_dest_comp = static_cast<ParamType>(_exchangeMatrix[rad_orig * _nRad * _nComp + rad_dest * _nComp + comp]);
if (exchange_orig_dest_comp == 0)
continue;
Expand Down

0 comments on commit 74059fd

Please sign in to comment.