Skip to content

Commit f01feea

Browse files
jbreue16jbreue16
committed
Fix/improve reaction config scope in particles
1 parent 85b4362 commit f01feea

File tree

7 files changed

+36
-93
lines changed

7 files changed

+36
-93
lines changed

‎src/libcadet/model/ColumnModel1D.cpp‎

Lines changed: 2 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -296,13 +296,8 @@ bool ColumnModel1D::configureModelDiscretization(IParameterProvider& paramProvid
296296
if (!_dynReactionBulk)
297297
throw InvalidParameterException("Unknown dynamic reaction model " + dynReactName);
298298

299-
if (_dynReactionBulk->usesParamProviderInDiscretizationConfig())
300-
paramProvider.pushScope("reaction_bulk");
301-
299+
MultiplexedScopeSelector scopeGuard(paramProvider, "reaction_bulk", _dynReactionBulk->usesParamProviderInDiscretizationConfig());
302300
reactionConfSuccess = _dynReactionBulk->configureModelDiscretization(paramProvider, _disc.nComp, nullptr, nullptr);
303-
304-
if (_dynReactionBulk->usesParamProviderInDiscretizationConfig())
305-
paramProvider.popScope();
306301
}
307302

308303
// ==== Construct and configure binding and particle reaction -> done in particle model, only pointers are copied here.
@@ -463,9 +458,8 @@ bool ColumnModel1D::configure(IParameterProvider& paramProvider)
463458
bool dynReactionConfSuccess = true;
464459
if (_dynReactionBulk && _dynReactionBulk->requiresConfiguration())
465460
{
466-
paramProvider.pushScope("reaction_bulk");
461+
MultiplexedScopeSelector scopeGuard(paramProvider, "reaction_bulk", _dynReactionBulk->requiresConfiguration());
467462
dynReactionConfSuccess = _dynReactionBulk->configure(paramProvider, _unitOpIdx, ParTypeIndep);
468-
paramProvider.popScope();
469463
}
470464

471465
// jaobian pattern set after binding and particle surface diffusion are configured

‎src/libcadet/model/ColumnModel1D.hpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -501,7 +501,7 @@ class ColumnModel1D : public UnitOperationBase
501501
}
502502
}
503503

504-
// particle jacobian (including film diffusion, isotherm and time derivative)
504+
// particle jacobian (including film diffusion, isotherm, reaction and time derivative)
505505
for (int colNode = 0; colNode < _disc.nPoints; colNode++)
506506
{
507507
for (int type = 0; type < _disc.nParType; type++)
@@ -519,7 +519,7 @@ class ColumnModel1D : public UnitOperationBase
519519
// inlet and bulk jacobian
520520
_convDispOp.calcTransportJacobian(_globalJac, _jacInlet, idxr.offsetC());
521521

522-
// particle jacobian (without isotherm, which is handled in residualKernel)
522+
// particle transport jacobian
523523
for (int colNode = 0; colNode < _disc.nPoints; colNode++)
524524
{
525525
for (int parType = 0; parType < _disc.nParType; parType++)

‎src/libcadet/model/ColumnModel2D.cpp‎

Lines changed: 2 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -508,13 +508,8 @@ bool ColumnModel2D::configureModelDiscretization(IParameterProvider& paramProvid
508508
if (!_dynReactionBulk)
509509
throw InvalidParameterException("Unknown dynamic reaction model " + dynReactName);
510510

511-
if (_dynReactionBulk->usesParamProviderInDiscretizationConfig())
512-
paramProvider.pushScope("reaction_bulk");
513-
511+
MultiplexedScopeSelector scopeGuard(paramProvider, "reaction_bulk", _dynReactionBulk->usesParamProviderInDiscretizationConfig());
514512
reactionConfSuccess = _dynReactionBulk->configureModelDiscretization(paramProvider, _disc.nComp, nullptr, nullptr);
515-
516-
if (_dynReactionBulk->usesParamProviderInDiscretizationConfig())
517-
paramProvider.popScope();
518513
}
519514

520515
// ==== Construct and configure binding and particle reaction -> done in particle model, only pointers are copied here.
@@ -696,9 +691,8 @@ bool ColumnModel2D::configure(IParameterProvider& paramProvider)
696691
bool dynReactionConfSuccess = true;
697692
if (_dynReactionBulk && _dynReactionBulk->requiresConfiguration())
698693
{
699-
paramProvider.pushScope("reaction_bulk");
694+
MultiplexedScopeSelector scopeGuard(paramProvider, "reaction_bulk", _dynReactionBulk->requiresConfiguration());
700695
dynReactionConfSuccess = _dynReactionBulk->configure(paramProvider, _unitOpIdx, ParTypeIndep);
701-
paramProvider.popScope();
702696
}
703697

704698
setJacobianPattern(_globalJac, 0, _dynReactionBulk);

‎src/libcadet/model/particle/GeneralRateParticle.cpp‎

Lines changed: 14 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -65,27 +65,15 @@ namespace model
6565
// ==== Construct binding model
6666

6767
_binding = nullptr;
68-
std::vector<std::string> bindModelNames = { "NONE" };
68+
std::string bindModelName = "NONE";
6969
bool bindingConfSuccess = true;
7070

7171
if (paramProvider.exists("ADSORPTION_MODEL"))
72-
{
73-
bindModelNames = paramProvider.getStringArray("ADSORPTION_MODEL");
74-
75-
if (_bindingParDep && (bindModelNames.size() != 1))
76-
throw InvalidParameterException("Field ADSORPTION_MODEL requires (only) 1 element");
72+
bindModelName = paramProvider.getString("ADSORPTION_MODEL");
7773

78-
if (paramProvider.exists("adsorption"))
79-
{
80-
paramProvider.pushScope("adsorption");
81-
_bindingParDep = paramProvider.exists("BINDING_PARTYPE_DEPENDENT") ? paramProvider.getInt("BINDING_PARTYPE_DEPENDENT") : true;
82-
paramProvider.popScope();
83-
}
84-
}
85-
86-
_binding = helper.createBindingModel(bindModelNames[0]);
74+
_binding = helper.createBindingModel(bindModelName);
8775
if (!_binding)
88-
throw InvalidParameterException("Unknown binding model " + bindModelNames[0]);
76+
throw InvalidParameterException("Unknown binding model " + bindModelName);
8977

9078
std::vector<int> nBound = paramProvider.getIntArray("NBOUND");
9179
if (nBound.size() != _nComp)
@@ -95,21 +83,10 @@ namespace model
9583
_nBound = std::make_shared<unsigned int[]>(_nComp);
9684
std::copy_n(nBound.begin(), _nComp, _nBound.get());
9785

98-
_bindingParDep = true;
99-
100-
if (_binding->usesParamProviderInDiscretizationConfig())
101-
{
102-
paramProvider.pushScope("adsorption");
103-
104-
if (paramProvider.exists("BINDING_PARTYPE_DEPENDENT"))
105-
_bindingParDep = paramProvider.getBool("BINDING_PARTYPE_DEPENDENT");
106-
107-
paramProvider.popScope();
108-
}
109-
else if (bindModelNames[0] == "NONE")
86+
if (bindModelName == "NONE")
11087
_nBound = std::make_shared<unsigned int[]>(_nComp, 0);
11188
else
112-
throw InvalidParameterException("Binding model " + bindModelNames[0] + " was specified, but group \"adsorption\" is missing for particle type " + std::to_string(_parTypeIdx));
89+
throw InvalidParameterException("Binding model " + bindModelName + " was specified, but group \"adsorption\" is missing for particle type " + std::to_string(_parTypeIdx));
11390

11491
// ==== Construct and configure particle transport and discretization
11592

@@ -133,15 +110,12 @@ namespace model
133110
const bool particleTransportConfigSuccess = _parDiffOp->configureModelDiscretization(paramProvider, helper, nComp, parTypeIdx, nParType, strideBulkComp);
134111

135112
// ==== Configure binding model
136-
137-
if (_binding->usesParamProviderInDiscretizationConfig())
138-
paramProvider.pushScope("adsorption");
139-
140-
bindingConfSuccess = _binding->configureModelDiscretization(paramProvider, _nComp, _nBound.get(), _parDiffOp->offsetBoundComp());
141-
142-
if (_binding->usesParamProviderInDiscretizationConfig())
143-
paramProvider.popScope();
144-
113+
_bindingParDep = true;
114+
{
115+
MultiplexedScopeSelector scopeGuard(paramProvider, "adsorption", _binding->usesParamProviderInDiscretizationConfig());
116+
bindingConfSuccess = _binding->configureModelDiscretization(paramProvider, _nComp, _nBound.get(), _parDiffOp->offsetBoundComp());
117+
_bindingParDep = paramProvider.exists("BINDING_PARTYPE_DEPENDENT") ? paramProvider.getInt("BINDING_PARTYPE_DEPENDENT") : true;
118+
}
145119
// ==== Construct and configure dynamic reaction model
146120

147121
_dynReaction = nullptr;
@@ -160,14 +134,11 @@ namespace model
160134
if (!_dynReaction)
161135
throw InvalidParameterException("Unknown dynamic reaction model " + dynReactModelNames[0]);
162136

163-
if (paramProvider.exists("reaction"))
164137
{
165-
paramProvider.pushScope("reaction");
138+
MultiplexedScopeSelector scopeGuard(paramProvider, "reaction", _dynReaction->usesParamProviderInDiscretizationConfig());
139+
reactionConfSuccess = _dynReaction->configureModelDiscretization(paramProvider, _nComp, _nBound.get(), _parDiffOp->offsetBoundComp()) && reactionConfSuccess;
166140
_reactionParDep = paramProvider.exists("REACTION_PARTYPE_DEPENDENT") ? paramProvider.getInt("REACTION_PARTYPE_DEPENDENT") : true;
167-
paramProvider.popScope();
168141
}
169-
170-
reactionConfSuccess = _dynReaction->configureModelDiscretization(paramProvider, _nComp, _nBound.get(), _parDiffOp->offsetBoundComp()) && reactionConfSuccess;
171142
}
172143

173144
paramProvider.popScope(); // particle_type_{:03}

‎src/libcadet/model/particle/HomogeneousParticle.cpp‎

Lines changed: 10 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -129,28 +129,16 @@ namespace model
129129

130130
_bindingParDep = true;
131131

132-
if (_binding->usesParamProviderInDiscretizationConfig())
133-
{
134-
paramProvider.pushScope("adsorption");
135-
136-
if (paramProvider.exists("BINDING_PARTYPE_DEPENDENT"))
137-
_bindingParDep = paramProvider.getBool("BINDING_PARTYPE_DEPENDENT");
138-
139-
paramProvider.popScope();
140-
}
141-
else if (bindModelNames[0] == "NONE")
132+
if (bindModelNames[0] == "NONE")
142133
_nBound = std::make_shared<unsigned int[]>(_nComp, 0);
143134
else
144135
throw InvalidParameterException("Binding model " + bindModelNames[0] + " was specified, but group \"adsorption\" is missing for particle type " + std::to_string(_parTypeIdx));
145136

146-
if (_binding->usesParamProviderInDiscretizationConfig())
147-
paramProvider.pushScope("adsorption");
148-
149-
bindingConfSuccess = _binding->configureModelDiscretization(paramProvider, _nComp, _nBound.get(), _boundOffset);
150-
151-
if (_binding->usesParamProviderInDiscretizationConfig())
152-
paramProvider.popScope();
153-
137+
{
138+
MultiplexedScopeSelector scopeGuard(paramProvider, "adsorption", _binding->usesParamProviderInDiscretizationConfig());
139+
bindingConfSuccess = _binding->configureModelDiscretization(paramProvider, _nComp, _nBound.get(), _boundOffset);
140+
_bindingParDep = paramProvider.getBool("BINDING_PARTYPE_DEPENDENT");
141+
}
154142
// ==== Construct and configure dynamic reaction model
155143

156144
_dynReaction = nullptr;
@@ -169,14 +157,11 @@ namespace model
169157
if (!_dynReaction)
170158
throw InvalidParameterException("Unknown dynamic reaction model " + dynReactModelNames[0]);
171159

172-
if (paramProvider.exists("reaction"))
173160
{
174-
paramProvider.pushScope("reaction");
161+
MultiplexedScopeSelector scopeGuard(paramProvider, "reaction", _dynReaction->usesParamProviderInDiscretizationConfig());
162+
reactionConfSuccess = _dynReaction->configureModelDiscretization(paramProvider, _nComp, _nBound.get(), _boundOffset) && reactionConfSuccess;
175163
_reactionParDep = paramProvider.exists("REACTION_PARTYPE_DEPENDENT") ? paramProvider.getInt("REACTION_PARTYPE_DEPENDENT") : true;
176-
paramProvider.popScope();
177164
}
178-
179-
reactionConfSuccess = _dynReaction->configureModelDiscretization(paramProvider, _nComp, _nBound.get(), _boundOffset) && reactionConfSuccess;
180165
}
181166

182167
paramProvider.popScope(); // particle_type_{:03}
@@ -420,24 +405,22 @@ namespace model
420405

421406
void HomogeneousParticle::setParJacPattern(std::vector<Eigen::Triplet<double>>& tripletList, const unsigned int offsetPar, const unsigned int offsetBulk, unsigned int colNode, unsigned int secIdx) const
422407
{
423-
// binding Jacobian pattern
424-
// add dense nComp * nBound blocks, since all solid and liquid entries can be coupled through binding.
408+
// binding and reaction Jacobian pattern
409+
// add dense nComp * nBound blocks, since all solid and liquid entries can be coupled through reaction and or binding.
425410
for (unsigned int parState = 0; parState < _nComp + _strideBound; parState++) {
426411
for (unsigned int toParState = 0; toParState < _nComp + _strideBound; toParState++) {
427412
tripletList.push_back(Eigen::Triplet<double>(offsetPar + parState, offsetPar + toParState, 0.0));
428413
}
429414
}
430415

431416
// flux Jacobian pattern
432-
433417
for (unsigned int comp = 0; comp < _nComp; comp++)
434418
{
435419
tripletList.push_back(Eigen::Triplet<double>(offsetPar + comp, offsetBulk + comp, 0.0));
436420
tripletList.push_back(Eigen::Triplet<double>(offsetBulk + comp, offsetPar + comp, 0.0));
437421
}
438422
}
439423

440-
441424
unsigned int HomogeneousParticle::jacobianNNZperParticle() const
442425
{
443426
return (_nComp + _strideBound) * (_nComp + _strideBound) + _nComp * 4; // reaction, binding patter + film diffusion pattern for one particle

‎src/libcadet/model/parts/ParticleDiffusionOperatorDG.cpp‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -912,11 +912,12 @@ namespace parts
912912
// ======================================== DG particle Jacobian ============================================= //
913913
// ========================================================================================================================================================== //
914914

915-
void ParticleDiffusionOperatorDG::parBindingPattern(std::vector<Eigen::Triplet<double>>& tripletList, const int offset, const unsigned int colNode)
915+
void ParticleDiffusionOperatorDG::parBindingAndReactionPattern(std::vector<Eigen::Triplet<double>>& tripletList, const int offset, const unsigned int colNode)
916916
{
917-
// every bound state might depend on every bound and liquid state
918917
for (int parNode = 0; parNode < _nParPoints; parNode++)
919918
{
919+
// binding: every bound state might depend on every bound and liquid state
920+
// reaction: every liquid and bound state might depend on every liquid and bound state
920921
for (int state = 0; state < strideParPoint(); state++)
921922
{
922923
for (int stateDep = 0; stateDep < strideParPoint(); stateDep++) {
@@ -928,14 +929,14 @@ namespace parts
928929
}
929930
}
930931
}
932+
931933
/**
932934
* @brief calculates the particle dispersion jacobian Pattern, including entries for the dependence of particle entries on bulk entries through film diffusion boundary condition
933935
* @detail Does NOT add film diffusion entries for the dependence of bulk conc. on particle conc.
934936
*/
935937
void ParticleDiffusionOperatorDG::setParticleJacobianPattern(std::vector<ParticleDiffusionOperatorDG::T>& tripletList, unsigned int offsetPar, unsigned int offsetBulk, unsigned int colNode, unsigned int secIdx)
936938
{
937-
938-
parBindingPattern(tripletList, offsetPar, colNode);
939+
parBindingAndReactionPattern(tripletList, offsetPar, colNode);
939940

940941
// Ordering of particle surface diffusion:
941942
// bnd0comp0, bnd1comp0, bnd0comp1, bnd1comp1, bnd0comp2, bnd1comp2

‎src/libcadet/model/parts/ParticleDiffusionOperatorDG.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -143,7 +143,7 @@ namespace parts
143143

144144
protected:
145145

146-
void parBindingPattern(std::vector<Eigen::Triplet<double>>& tripletList, const int offset, const unsigned int colNode);
146+
void parBindingAndReactionPattern(std::vector<Eigen::Triplet<double>>& tripletList, const int offset, const unsigned int colNode);
147147

148148
template <typename StateType, typename ResidualType, typename ParamType, bool wantJac, bool wantRes>
149149
int residualImpl(double t, unsigned int secIdx, StateType const* yPar, StateType const* yBulk, double const* yDotPar, ResidualType* resPar, linalg::BandedEigenSparseRowIterator& jacBase);

0 commit comments

Comments
 (0)