Skip to content

Commit 3882248

Browse files
AntoniaBergerAntoniaBerger
committed
Delete _dynRection and _dynReactionBulk from unit operations
1 parent 2cd0d4c commit 3882248

17 files changed

+43
-85
lines changed

‎src/libcadet/model/ColumnModel1D-InitialConditions.cpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -613,7 +613,8 @@ void ColumnModel1D::consistentInitialState(const SimulationTime& simTime, double
613613
}
614614

615615
// reset jacobian pattern
616-
setJacobianPattern(_globalJacDisc, _disc.curSection, _dynReactionBulk);
616+
bool hasBulkReaction = _reaction.hasReactions();
617+
setJacobianPattern(_globalJacDisc, _disc.curSection, hasBulkReaction);
617618

618619
} CADET_PARFOR_END;
619620
}

‎src/libcadet/model/ColumnModel1D.cpp‎

Lines changed: 1 addition & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -57,7 +57,7 @@ constexpr double SurfVolRatioSlab = 1.0;
5757

5858

5959
ColumnModel1D::ColumnModel1D(UnitOpIdx unitOpIdx) : UnitOperationBase(unitOpIdx),
60-
_globalJac(), _globalJacDisc(), _jacInlet(), _dynReactionBulk(nullptr),
60+
_globalJac(), _globalJacDisc(), _jacInlet(),
6161
_analyticJac(true), _jacobianAdDirs(0), _factorizeJacobian(false), _tempState(nullptr),
6262
_initC(0), _initCp(0), _initCs(0), _initState(0), _initStateDot(0)
6363
{
@@ -75,8 +75,6 @@ ColumnModel1D::~ColumnModel1D() CADET_NOEXCEPT
7575
_binding.clear(); // binding models are deleted in the respective particle model
7676
_dynReaction.clear(); // particle reaction models are deleted in the respective particle model
7777

78-
delete _dynReactionBulk;
79-
8078
_reaction.clearDynamicReactionModels();
8179
delete _linearSolver;
8280
}
@@ -337,7 +335,6 @@ bool ColumnModel1D::configureModelDiscretization(IParameterProvider& paramProvid
337335

338336
// ==== Construct and configure binding and particle reaction -> done in particle model, only pointers are copied here.
339337
_binding = std::vector<IBindingModel*>(_disc.nParType, nullptr);
340-
_dynReaction = std::vector<IDynamicReactionModel*>(_disc.nParType, nullptr);
341338
_reacParticle.resize(_disc.nParType, nullptr);
342339

343340
for (unsigned int parType = 0; parType < _disc.nParType; ++parType)

‎src/libcadet/model/ColumnModel1D.hpp‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -297,9 +297,8 @@ class ColumnModel1D : public UnitOperationBase
297297
std::vector<IParticleModel*> _particles; //!< Particle dispersion operator
298298

299299
parts::AxialConvectionDispersionOperatorBaseDG _convDispOp; //!< Convection dispersion operator base for interstitial volume transport
300-
IDynamicReactionModel* _dynReactionBulk; //!< Dynamic reactions in the bulk volume
301-
ReactionSystem _reaction;
302-
std::vector<const ReactionSystem*> _reacParticle;
300+
ReactionSystem _reaction; //!< Reaction system for bulk phase
301+
std::vector<const ReactionSystem*> _reacParticle; //!< Reaction systems for particle phases
303302

304303
cadet::linalg::EigenSolverBase* _linearSolver; //!< Linear solver
305304

‎src/libcadet/model/ColumnModel2D-InitialConditions.cpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -839,7 +839,8 @@ void ColumnModel2D::consistentInitialState(const SimulationTime& simTime, double
839839
}
840840

841841
// reset jacobian pattern //@todo can this be avoided?
842-
setJacobianPattern(_globalJacDisc, simTime.secIdx, _dynReactionBulk);
842+
bool hasBulkReactions = _reaction.hasReactions();
843+
setJacobianPattern(_globalJacDisc, simTime.secIdx, hasBulkReactions);
843844
}
844845

845846
/**

‎src/libcadet/model/ColumnModel2D.cpp‎

Lines changed: 1 addition & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -317,7 +317,7 @@ namespace model
317317
{
318318

319319
ColumnModel2D::ColumnModel2D(UnitOpIdx unitOpIdx) : UnitOperationBase(unitOpIdx),
320-
_dynReactionBulk(nullptr), _jacInlet(), _analyticJac(true), _jacobianAdDirs(0), _factorizeJacobian(false), _tempState(nullptr),
320+
_jacInlet(), _analyticJac(true), _jacobianAdDirs(0), _factorizeJacobian(false), _tempState(nullptr),
321321
_initC(0), _singleRadiusInitC(true), _initCp(0), _singleRadiusInitCp(true), _initCs(0), _singleRadiusInitCs(true), _initState(0), _initStateDot(0)
322322
{
323323
}
@@ -334,7 +334,6 @@ ColumnModel2D::~ColumnModel2D() CADET_NOEXCEPT
334334
_binding.clear(); // binding models are deleted in the respective particle model
335335
_dynReaction.clear(); // particle reaction models are deleted in the respective particle model
336336

337-
delete _dynReactionBulk;
338337

339338
delete[] _disc.parTypeOffset;
340339
delete[] _disc.nBound;
@@ -534,7 +533,6 @@ bool ColumnModel2D::configureModelDiscretization(IParameterProvider& paramProvid
534533

535534
// ==== Construct and configure binding and particle reaction -> done in particle model, only pointers are copied here.
536535
_binding = std::vector<IBindingModel*>(_disc.nParType, nullptr);
537-
_dynReaction = std::vector<IDynamicReactionModel*>(_disc.nParType, nullptr);
538536
_reacParticle = std::vector<const ReactionSystem*>(_disc.nParType, nullptr);
539537
for (unsigned int parType = 0; parType < _disc.nParType; ++parType)
540538
{
@@ -549,17 +547,10 @@ bool ColumnModel2D::configureModelDiscretization(IParameterProvider& paramProvid
549547
if (_singleBinding != !_particles[parType]->bindingParDep())
550548
throw InvalidParameterException("Binding particle type dependence must be the same for all particle types, check field BINDING_PARTYPE_DEPENDENT");
551549
}
552-
553-
if (_dynReaction[parType])
554-
{
555-
if (_singleDynReaction != !_particles[parType]->reactionParDep())
556-
throw InvalidParameterException("Reaction particle type dependence must be the same for all particle types, check field REACTION_PARTYPE_DEPENDENT");
557-
}
558550
}
559551
else // if no particle reaction or binding exists in first particle type, default to single mode
560552
{
561553
_singleBinding = _binding[parType] ? !_particles[parType]->bindingParDep() : true;
562-
_singleDynReaction = _dynReaction[parType] ? !_particles[parType]->reactionParDep() : true;
563554
}
564555
}
565556

‎src/libcadet/model/ColumnModel2D.hpp‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -303,9 +303,8 @@ class ColumnModel2D : public UnitOperationBase
303303
std::vector<IParticleModel*> _particles; //!< Particle dispersion operator
304304

305305
parts::TwoDimensionalConvectionDispersionOperatorDG _convDispOp; //!< Convection dispersion operator for interstitial volume transport
306-
IDynamicReactionModel* _dynReactionBulk; //!< Dynamic reactions in the bulk volume
307-
ReactionSystem _reaction;
308-
std::vector<const ReactionSystem*> _reacParticle;
306+
ReactionSystem _reaction; //!< Reaction system for bulk phase
307+
std::vector<const ReactionSystem*> _reacParticle; //!< Reaction systems for particle phases
309308

310309
cadet::linalg::EigenSolverBase* _linearSolver; //!< Linear solver
311310

‎src/libcadet/model/GeneralRateModel.cpp‎

Lines changed: 8 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -66,7 +66,7 @@ int schurComplementMultiplierGRM(void* userData, double const* x, double* z)
6666

6767
template <typename ConvDispOperator>
6868
GeneralRateModel<ConvDispOperator>::GeneralRateModel(UnitOpIdx unitOpIdx) : UnitOperationBase(unitOpIdx),
69-
_hasSurfaceDiffusion(0, false), _dynReactionBulk(nullptr),
69+
_hasSurfaceDiffusion(0, false),
7070
_jacP(nullptr), _jacPdisc(nullptr), _jacPF(nullptr), _jacFP(nullptr), _jacInlet(), _hasParDepSurfDiffusion(false),
7171
_analyticJac(true), _jacobianAdDirs(0), _factorizeJacobian(false), _tempState(nullptr),
7272
_initC(0), _initCp(0), _initCs(0), _initState(0), _initStateDot(0)
@@ -84,7 +84,8 @@ GeneralRateModel<ConvDispOperator>::~GeneralRateModel() CADET_NOEXCEPT
8484
delete[] _jacP;
8585
delete[] _jacPdisc;
8686

87-
delete _dynReactionBulk;
87+
_reaction.empty();
88+
_dynReaction.clear();
8889

8990
clearParDepSurfDiffusion();
9091
}
@@ -202,7 +203,6 @@ bool GeneralRateModel<ConvDispOperator>::configureModelDiscretization(IParameter
202203
_binding = std::vector<IBindingModel*>(_disc.nParType, nullptr);
203204
bool bindingConfSuccess = true;
204205
clearDynamicReactionModels();
205-
_dynReaction = std::vector<IDynamicReactionModel*>(_disc.nParType, nullptr);
206206
bool reactionConfSuccess = true;
207207

208208
for (int parType = 0; parType < _disc.nParType; parType++)
@@ -286,7 +286,6 @@ bool GeneralRateModel<ConvDispOperator>::configureModelDiscretization(IParameter
286286
_reaction.clearDynamicReactionModels();
287287
_reacParticle.clear();
288288
_reacParticle.resize(_disc.nParType);
289-
_dynReaction.resize(_disc.nParType, nullptr);
290289

291290
if (_disc.nParType > 0)
292291
{
@@ -2914,7 +2913,7 @@ bool GeneralRateModel<ConvDispOperator>::setParameter(const ParameterId& pId, do
29142913
if (_convDispOp.setParameter(pId, value))
29152914
return true;
29162915

2917-
if (model::setParameter(pId, value, std::vector<IDynamicReactionModel*>{ _dynReactionBulk }, true))
2916+
if (model::setParameter(pId, value, _reaction.getDynReactionVector("liquid"), true))
29182917
return true;
29192918
}
29202919

@@ -2938,7 +2937,7 @@ bool GeneralRateModel<ConvDispOperator>::setParameter(const ParameterId& pId, in
29382937

29392938
if (pId.unitOperation == _unitOpIdx)
29402939
{
2941-
if (model::setParameter(pId, value, std::vector<IDynamicReactionModel*>{ _dynReactionBulk }, true))
2940+
if (model::setParameter(pId, value, _reaction.getDynReactionVector("liquid"), true))
29422941
return true;
29432942
}
29442943

@@ -2956,7 +2955,7 @@ bool GeneralRateModel<ConvDispOperator>::setParameter(const ParameterId& pId, bo
29562955

29572956
if (pId.unitOperation == _unitOpIdx)
29582957
{
2959-
if (model::setParameter(pId, value, std::vector<IDynamicReactionModel*>{ _dynReactionBulk }, true))
2958+
if (model::setParameter(pId, value, _reaction.getDynReactionVector("liquid"), true))
29602959
return true;
29612960
}
29622961

@@ -3009,7 +3008,7 @@ void GeneralRateModel<ConvDispOperator>::setSensitiveParameterValue(const Parame
30093008
if (_convDispOp.setSensitiveParameterValue(_sensParams, pId, value))
30103009
return;
30113010

3012-
if (model::setSensitiveParameterValue(pId, value, _sensParams, std::vector<IDynamicReactionModel*>{ _dynReactionBulk }, true))
3011+
if (model::setSensitiveParameterValue(pId, value, _sensParams, _reaction.getDynReactionVector("liquid"), false))
30133012
return;
30143013
}
30153014

@@ -3109,7 +3108,7 @@ bool GeneralRateModel<ConvDispOperator>::setSensitiveParameter(const ParameterId
31093108
return true;
31103109
}
31113110

3112-
if (model::setSensitiveParameter(pId, adDirection, adValue, _sensParams, std::vector<IDynamicReactionModel*> { _dynReactionBulk }, true))
3111+
if (model::setSensitiveParameter(pId, adDirection, adValue, _sensParams, _reaction.getDynReactionVector("liquid"), false))
31133112
{
31143113
LOG(Debug) << "Found parameter " << pId << " in DynamicBulkReactionModel: Dir " << adDirection << " is set to " << adValue;
31153114
return true;

‎src/libcadet/model/GeneralRateModel.hpp‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -333,10 +333,9 @@ class GeneralRateModel : public UnitOperationBase
333333
// IExternalFunction* _extFun; //!< External function (owned by library user)
334334

335335
ConvDispOperator _convDispOp; //!< Convection dispersion operator for interstitial volume transport
336-
ReactionSystem _reaction;
337-
std::vector<ReactionSystem> _reacParticle;
336+
ReactionSystem _reaction; //!< Reaction system for bulk
337+
std::vector<ReactionSystem> _reacParticle; //!< Reaction systems for each particle type
338338

339-
IDynamicReactionModel* _dynReactionBulk; //!< Dynamic reactions in the bulk volume
340339

341340
linalg::BandMatrix* _jacP; //!< Particle jacobian diagonal blocks (all of them)
342341
linalg::FactorizableBandMatrix* _jacPdisc; //!< Particle jacobian diagonal blocks (all of them) with time derivatives from BDF method

‎src/libcadet/model/GeneralRateModel2D.cpp‎

Lines changed: 9 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -331,7 +331,7 @@ int schurComplementMultiplierGRM2D(void* userData, double const* x, double* z)
331331

332332

333333
GeneralRateModel2D::GeneralRateModel2D(UnitOpIdx unitOpIdx) : UnitOperationBase(unitOpIdx),
334-
_dynReactionBulk(nullptr), _jacP(nullptr), _jacPdisc(nullptr), _jacPF(nullptr), _jacFP(nullptr), _jacInlet(),
334+
_jacP(nullptr), _jacPdisc(nullptr), _jacPF(nullptr), _jacFP(nullptr), _jacInlet(),
335335
_analyticJac(true), _jacobianAdDirs(0), _factorizeJacobian(false), _tempState(nullptr),
336336
_initC(0), _singleRadiusInitC(true), _initCp(0), _singleRadiusInitCp(true), _initCs(0), _singleRadiusInitCs(true), _initState(0), _initStateDot(0)
337337
{
@@ -341,8 +341,6 @@ GeneralRateModel2D::~GeneralRateModel2D() CADET_NOEXCEPT
341341
{
342342
delete[] _tempState;
343343

344-
delete _dynReactionBulk;
345-
346344
delete[] _jacPF;
347345
delete[] _jacFP;
348346

@@ -419,7 +417,6 @@ bool GeneralRateModel2D::configureModelDiscretization(IParameterProvider& paramP
419417
_binding = std::vector<IBindingModel*>(_disc.nParType, nullptr);
420418
bool bindingConfSuccess = true;
421419
clearDynamicReactionModels();
422-
_dynReaction = std::vector<IDynamicReactionModel*>(_disc.nParType, nullptr);
423420
bool reactionConfSuccess = true;
424421

425422
for (int parType = 0; parType < _disc.nParType; parType++)
@@ -497,13 +494,6 @@ bool GeneralRateModel2D::configureModelDiscretization(IParameterProvider& paramP
497494
_singleBinding = !paramProvider.getInt("REACTION_PARTYPE_DEPENDENT");
498495
else
499496
_singleBinding = _disc.nParType == 1;
500-
501-
_dynReaction[parType] = helper.createDynamicReactionModel(dynReactModelName);
502-
if (!_dynReaction[parType])
503-
throw InvalidParameterException("Unknown dynamic reaction model " + dynReactModelName);
504-
505-
MultiplexedScopeSelector scopeGuard(paramProvider, "reaction", _dynReaction[parType]->usesParamProviderInDiscretizationConfig());
506-
reactionConfSuccess = reactionConfSuccess && _dynReaction[parType]->configureModelDiscretization(paramProvider, _disc.nComp, _disc.nBound + parType * _disc.nComp, _disc.boundOffset + parType * _disc.nComp) && reactionConfSuccess;
507497
}
508498

509499
// Set particle geometry
@@ -639,8 +629,6 @@ bool GeneralRateModel2D::configureModelDiscretization(IParameterProvider& paramP
639629

640630
// ==== Construct and configure dynamic bulk reaction model
641631

642-
_dynReaction.resize(_disc.nParType, nullptr);
643-
644632
_reaction.clearDynamicReactionModels();
645633

646634
_reacParticle.clear();
@@ -1419,7 +1407,7 @@ template <typename StateType, typename ResidualType, typename ParamType, bool wa
14191407
int GeneralRateModel2D::residualBulk(double t, unsigned int secIdx, StateType const* yBase, double const* yDotBase, ResidualType* resBase, util::ThreadLocalStorage& threadLocalMem)
14201408
{
14211409
_convDispOp.residual(*this, t, secIdx, yBase, yDotBase, resBase, wantJac, typename ParamSens<ParamType>::enabled());
1422-
if (!_dynReactionBulk || (_dynReactionBulk->numReactionsLiquid() == 0))
1410+
if (! _reaction.getDynReactionVector("liquid")[0])
14231411
return 0;
14241412

14251413
// Get offsets
@@ -1441,12 +1429,12 @@ int GeneralRateModel2D::residualBulk(double t, unsigned int secIdx, StateType co
14411429
if (!_reaction.getDynReactionVector("liquid")[i])
14421430
continue;
14431431

1444-
_dynReactionBulk->residualFluxAdd(t, secIdx, colPos, _disc.nComp, y, res, -1.0, tlmAlloc);
1432+
_reaction.getDynReactionVector("liquid")[i]->residualFluxAdd(t, secIdx, colPos, _disc.nComp, y, res, -1.0, tlmAlloc);
14451433

14461434
if (wantJac)
14471435
{
14481436
// static_cast should be sufficient here, but this statement is also analyzed when wantJac = false
1449-
_dynReactionBulk->analyticJacobianAdd(t, secIdx, colPos, _disc.nComp, reinterpret_cast<double const*>(y), -1.0, _convDispOp.jacobian().row(colCell * idxr.strideColRadialCell()), tlmAlloc);
1437+
_reaction.getDynReactionVector("liquid")[i]->analyticJacobianAdd(t, secIdx, colPos, _disc.nComp, reinterpret_cast<double const*>(y), -1.0, _convDispOp.jacobian().row(colCell * idxr.strideColRadialCell()), tlmAlloc);
14501438
}
14511439
}
14521440
}
@@ -2434,7 +2422,7 @@ bool GeneralRateModel2D::setParameter(const ParameterId& pId, double value)
24342422
if (_convDispOp.setParameter(pId, value))
24352423
return true;
24362424

2437-
if (model::setParameter(pId, value, std::vector<IDynamicReactionModel*>{ _dynReactionBulk }, true))
2425+
if (model::setParameter(pId, value, _reaction.getDynReactionVector("liquid"), true))
24382426
return true;
24392427
}
24402428

@@ -2454,7 +2442,7 @@ bool GeneralRateModel2D::setParameter(const ParameterId& pId, int value)
24542442

24552443
if (pId.unitOperation == _unitOpIdx)
24562444
{
2457-
if (model::setParameter(pId, value, std::vector<IDynamicReactionModel*>{ _dynReactionBulk }, true))
2445+
if (model::setParameter(pId, value, _reaction.getDynReactionVector("liquid"), true))
24582446
return true;
24592447
}
24602448

@@ -2468,7 +2456,7 @@ bool GeneralRateModel2D::setParameter(const ParameterId& pId, bool value)
24682456

24692457
if (pId.unitOperation == _unitOpIdx)
24702458
{
2471-
if (model::setParameter(pId, value, std::vector<IDynamicReactionModel*>{ _dynReactionBulk }, true))
2459+
if (model::setParameter(pId, value, _reaction.getDynReactionVector("liquid"), true))
24722460
return true;
24732461
}
24742462

@@ -2505,7 +2493,7 @@ void GeneralRateModel2D::setSensitiveParameterValue(const ParameterId& pId, doub
25052493
if (_convDispOp.setSensitiveParameterValue(_sensParams, pId, value))
25062494
return;
25072495

2508-
if (model::setSensitiveParameterValue(pId, value, _sensParams, std::vector<IDynamicReactionModel*>{ _dynReactionBulk }, true))
2496+
if (model::setSensitiveParameterValue(pId, value, _sensParams, _reaction.getDynReactionVector("liquid"), true))
25092497
return;
25102498
}
25112499

@@ -2586,7 +2574,7 @@ bool GeneralRateModel2D::setSensitiveParameter(const ParameterId& pId, unsigned
25862574
return true;
25872575
}
25882576

2589-
if (model::setSensitiveParameter(pId, adDirection, adValue, _sensParams, std::vector<IDynamicReactionModel*> { _dynReactionBulk }, true))
2577+
if (model::setSensitiveParameter(pId, adDirection, adValue, _sensParams, _reaction.getDynReactionVector("liquid"), true))
25902578
{
25912579
LOG(Debug) << "Found parameter " << pId << " in DynamicBulkReactionModel: Dir " << adDirection << " is set to " << adValue;
25922580
return true;

‎src/libcadet/model/GeneralRateModel2D.hpp‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -278,9 +278,8 @@ class GeneralRateModel2D : public UnitOperationBase
278278
// IExternalFunction* _extFun; //!< External function (owned by library user)
279279

280280
parts::TwoDimensionalConvectionDispersionOperator _convDispOp; //!< Convection dispersion operator for interstitial volume transport
281-
IDynamicReactionModel* _dynReactionBulk; //!< Dynamic reactions in the bulk volume
282-
ReactionSystem _reaction;
283-
std::vector<ReactionSystem> _reacParticle;
281+
ReactionSystem _reaction; //!< Reaction system for bulk
282+
std::vector<ReactionSystem> _reacParticle; //!< Reaction systems for particle bound states
284283

285284
linalg::BandMatrix* _jacP; //!< Particle jacobian diagonal blocks (all of them)
286285
linalg::FactorizableBandMatrix* _jacPdisc; //!< Particle jacobian diagonal blocks (all of them) with time derivatives from BDF method

0 commit comments

Comments
 (0)