Skip to content

Commit 38e0609

Browse files
jbreue16jbreue16
committed
Fix/improve reaction config scope in particles
1 parent 1196825 commit 38e0609

File tree

8 files changed

+55
-118
lines changed

8 files changed

+55
-118
lines changed

‎src/libcadet/model/ColumnModel1D.cpp‎

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

310-
if (_dynReactionBulk->usesParamProviderInDiscretizationConfig())
311-
paramProvider.pushScope("reaction_bulk");
312-
310+
MultiplexedScopeSelector scopeGuard(paramProvider, "reaction_bulk", _dynReactionBulk->usesParamProviderInDiscretizationConfig());
313311
reactionConfSuccess = _dynReactionBulk->configureModelDiscretization(paramProvider, _disc.nComp, nullptr, nullptr);
314-
315-
if (_dynReactionBulk->usesParamProviderInDiscretizationConfig())
316-
paramProvider.popScope();
317312
}
318313

319314
// ==== Construct and configure binding and particle reaction -> done in particle model, only pointers are copied here.
@@ -474,9 +469,8 @@ bool ColumnModel1D::configure(IParameterProvider& paramProvider)
474469
bool dynReactionConfSuccess = true;
475470
if (_dynReactionBulk && _dynReactionBulk->requiresConfiguration())
476471
{
477-
paramProvider.pushScope("reaction_bulk");
472+
MultiplexedScopeSelector scopeGuard(paramProvider, "reaction_bulk", _dynReactionBulk->requiresConfiguration());
478473
dynReactionConfSuccess = _dynReactionBulk->configure(paramProvider, _unitOpIdx, ParTypeIndep);
479-
paramProvider.popScope();
480474
}
481475

482476
// 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
@@ -511,13 +511,8 @@ bool ColumnModel2D::configureModelDiscretization(IParameterProvider& paramProvid
511511
if (!_dynReactionBulk)
512512
throw InvalidParameterException("Unknown dynamic reaction model " + dynReactName);
513513

514-
if (_dynReactionBulk->usesParamProviderInDiscretizationConfig())
515-
paramProvider.pushScope("reaction_bulk");
516-
514+
MultiplexedScopeSelector scopeGuard(paramProvider, "reaction_bulk", _dynReactionBulk->usesParamProviderInDiscretizationConfig());
517515
reactionConfSuccess = _dynReactionBulk->configureModelDiscretization(paramProvider, _disc.nComp, nullptr, nullptr);
518-
519-
if (_dynReactionBulk->usesParamProviderInDiscretizationConfig())
520-
paramProvider.popScope();
521516
}
522517

523518
// ==== Construct and configure binding and particle reaction -> done in particle model, only pointers are copied here.
@@ -699,9 +694,8 @@ bool ColumnModel2D::configure(IParameterProvider& paramProvider)
699694
bool dynReactionConfSuccess = true;
700695
if (_dynReactionBulk && _dynReactionBulk->requiresConfiguration())
701696
{
702-
paramProvider.pushScope("reaction_bulk");
697+
MultiplexedScopeSelector scopeGuard(paramProvider, "reaction_bulk", _dynReactionBulk->requiresConfiguration());
703698
dynReactionConfSuccess = _dynReactionBulk->configure(paramProvider, _unitOpIdx, ParTypeIndep);
704-
paramProvider.popScope();
705699
}
706700

707701
setJacobianPattern(_globalJac, 0, _dynReactionBulk);

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

Lines changed: 21 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -62,54 +62,31 @@ namespace model
6262
parTypeIdxString << std::setfill('0') << std::setw(3) << std::setprecision(0) << _parTypeIdx;
6363
paramProvider.pushScope("particle_type_" + parTypeIdxString.str());
6464

65-
// ==== Construct binding model
65+
// ==== Construct and configure 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]);
89-
90-
std::vector<int> nBound = paramProvider.getIntArray("NBOUND");
91-
if (nBound.size() != _nComp)
92-
throw InvalidParameterException("Field NBOUND does not contain NCOMP = " + std::to_string(_nComp) + " entries for particle type " + std::to_string(_parTypeIdx));
93-
94-
if (!_nBound)
95-
_nBound = std::make_shared<unsigned int[]>(_nComp);
96-
std::copy_n(nBound.begin(), _nComp, _nBound.get());
76+
throw InvalidParameterException("Unknown binding model " + bindModelName);
9777

98-
_bindingParDep = true;
99-
100-
if (_binding->usesParamProviderInDiscretizationConfig())
78+
if (bindModelName == "NONE")
79+
_nBound = std::make_shared<unsigned int[]>(_nComp, 0);
80+
else
10181
{
102-
paramProvider.pushScope("adsorption");
103-
104-
if (paramProvider.exists("BINDING_PARTYPE_DEPENDENT"))
105-
_bindingParDep = paramProvider.getBool("BINDING_PARTYPE_DEPENDENT");
82+
std::vector<int> nBound = paramProvider.getIntArray("NBOUND");
83+
if (nBound.size() != _nComp)
84+
throw InvalidParameterException("Field NBOUND does not contain NCOMP = " + std::to_string(_nComp) + " entries for particle type " + std::to_string(_parTypeIdx));
10685

107-
paramProvider.popScope();
86+
if (!_nBound)
87+
_nBound = std::make_shared<unsigned int[]>(_nComp);
88+
std::copy_n(nBound.begin(), _nComp, _nBound.get());
10889
}
109-
else if (bindModelNames[0] == "NONE")
110-
_nBound = std::make_shared<unsigned int[]>(_nComp, 0);
111-
else
112-
throw InvalidParameterException("Binding model " + bindModelNames[0] + " 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,14 +110,11 @@ 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();
113+
{
114+
MultiplexedScopeSelector scopeGuard(paramProvider, "adsorption", _binding->usesParamProviderInDiscretizationConfig());
115+
bindingConfSuccess = _binding->configureModelDiscretization(paramProvider, _nComp, _nBound.get(), _parDiffOp->offsetBoundComp());
116+
_bindingParDep = paramProvider.exists("BINDING_PARTYPE_DEPENDENT") ? paramProvider.getBool("BINDING_PARTYPE_DEPENDENT") : true;
117+
}
144118

145119
// ==== Construct and configure dynamic reaction model
146120

@@ -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: 22 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -105,51 +105,34 @@ namespace model
105105
// ==== Construct and configure binding model
106106

107107
_binding = nullptr;
108-
std::vector<std::string> bindModelNames = { "NONE" };
108+
std::string bindModelName = "NONE";
109109
bool bindingConfSuccess = true;
110110

111111
if (paramProvider.exists("ADSORPTION_MODEL"))
112-
{
113-
bindModelNames = paramProvider.getStringArray("ADSORPTION_MODEL");
114-
115-
if (_bindingParDep && (bindModelNames.size() != 1))
116-
throw InvalidParameterException("Field ADSORPTION_MODEL requires (only) 1 element");
117-
118-
if (paramProvider.exists("adsorption"))
119-
{
120-
paramProvider.pushScope("adsorption");
121-
_bindingParDep = paramProvider.exists("BINDING_PARTYPE_DEPENDENT") ? paramProvider.getInt("BINDING_PARTYPE_DEPENDENT") : true;
122-
paramProvider.popScope();
123-
}
124-
}
112+
bindModelName = paramProvider.getString("ADSORPTION_MODEL");
125113

126-
_binding = helper.createBindingModel(bindModelNames[0]);
114+
_binding = helper.createBindingModel(bindModelName);
127115
if (!_binding)
128-
throw InvalidParameterException("Unknown binding model " + bindModelNames[0]);
129-
130-
_bindingParDep = true;
131-
132-
if (_binding->usesParamProviderInDiscretizationConfig())
133-
{
134-
paramProvider.pushScope("adsorption");
135-
136-
if (paramProvider.exists("BINDING_PARTYPE_DEPENDENT"))
137-
_bindingParDep = paramProvider.getBool("BINDING_PARTYPE_DEPENDENT");
116+
throw InvalidParameterException("Unknown binding model " + bindModelName);
138117

139-
paramProvider.popScope();
140-
}
141-
else if (bindModelNames[0] == "NONE")
118+
if (bindModelName == "NONE")
142119
_nBound = std::make_shared<unsigned int[]>(_nComp, 0);
143120
else
144-
throw InvalidParameterException("Binding model " + bindModelNames[0] + " was specified, but group \"adsorption\" is missing for particle type " + std::to_string(_parTypeIdx));
145-
146-
if (_binding->usesParamProviderInDiscretizationConfig())
147-
paramProvider.pushScope("adsorption");
121+
{
122+
std::vector<int> nBound = paramProvider.getIntArray("NBOUND");
123+
if (nBound.size() != _nComp)
124+
throw InvalidParameterException("Field NBOUND does not contain NCOMP = " + std::to_string(_nComp) + " entries for particle type " + std::to_string(_parTypeIdx));
148125

149-
bindingConfSuccess = _binding->configureModelDiscretization(paramProvider, _nComp, _nBound.get(), _boundOffset);
126+
if (!_nBound)
127+
_nBound = std::make_shared<unsigned int[]>(_nComp);
128+
std::copy_n(nBound.begin(), _nComp, _nBound.get());
129+
}
150130

151-
if (_binding->usesParamProviderInDiscretizationConfig())
152-
paramProvider.popScope();
131+
{
132+
MultiplexedScopeSelector scopeGuard(paramProvider, "adsorption", _binding->usesParamProviderInDiscretizationConfig());
133+
bindingConfSuccess = _binding->configureModelDiscretization(paramProvider, _nComp, _nBound.get(), _boundOffset);
134+
_bindingParDep = paramProvider.exists("BINDING_PARTYPE_DEPENDENT") ? paramProvider.getBool("BINDING_PARTYPE_DEPENDENT") : true;
135+
}
153136

154137
// ==== Construct and configure dynamic reaction model
155138

@@ -169,14 +152,11 @@ namespace model
169152
if (!_dynReaction)
170153
throw InvalidParameterException("Unknown dynamic reaction model " + dynReactModelNames[0]);
171154

172-
if (paramProvider.exists("reaction"))
173155
{
174-
paramProvider.pushScope("reaction");
156+
MultiplexedScopeSelector scopeGuard(paramProvider, "reaction", _dynReaction->usesParamProviderInDiscretizationConfig());
157+
reactionConfSuccess = _dynReaction->configureModelDiscretization(paramProvider, _nComp, _nBound.get(), _boundOffset) && reactionConfSuccess;
175158
_reactionParDep = paramProvider.exists("REACTION_PARTYPE_DEPENDENT") ? paramProvider.getInt("REACTION_PARTYPE_DEPENDENT") : true;
176-
paramProvider.popScope();
177159
}
178-
179-
reactionConfSuccess = _dynReaction->configureModelDiscretization(paramProvider, _nComp, _nBound.get(), _boundOffset) && reactionConfSuccess;
180160
}
181161

182162
paramProvider.popScope(); // particle_type_{:03}
@@ -420,24 +400,22 @@ namespace model
420400

421401
void HomogeneousParticle::setParJacPattern(std::vector<Eigen::Triplet<double>>& tripletList, const unsigned int offsetPar, const unsigned int offsetBulk, unsigned int colNode, unsigned int secIdx) const
422402
{
423-
// binding Jacobian pattern
424-
// add dense nComp * nBound blocks, since all solid and liquid entries can be coupled through binding.
403+
// binding and reaction Jacobian pattern
404+
// add dense nComp * nBound blocks, since all solid and liquid entries can be coupled through reaction and or binding.
425405
for (unsigned int parState = 0; parState < _nComp + _strideBound; parState++) {
426406
for (unsigned int toParState = 0; toParState < _nComp + _strideBound; toParState++) {
427407
tripletList.push_back(Eigen::Triplet<double>(offsetPar + parState, offsetPar + toParState, 0.0));
428408
}
429409
}
430410

431411
// flux Jacobian pattern
432-
433412
for (unsigned int comp = 0; comp < _nComp; comp++)
434413
{
435414
tripletList.push_back(Eigen::Triplet<double>(offsetPar + comp, offsetBulk + comp, 0.0));
436415
tripletList.push_back(Eigen::Triplet<double>(offsetBulk + comp, offsetPar + comp, 0.0));
437416
}
438417
}
439418

440-
441419
unsigned int HomogeneousParticle::jacobianNNZperParticle() const
442420
{
443421
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
@@ -915,11 +915,12 @@ namespace parts
915915
// ======================================== DG particle Jacobian ============================================= //
916916
// ========================================================================================================================================================== //
917917

918-
void ParticleDiffusionOperatorDG::parBindingPattern(std::vector<Eigen::Triplet<double>>& tripletList, const int offset, const unsigned int colNode)
918+
void ParticleDiffusionOperatorDG::parBindingAndReactionPattern(std::vector<Eigen::Triplet<double>>& tripletList, const int offset, const unsigned int colNode)
919919
{
920-
// every bound state might depend on every bound and liquid state
921920
for (int parNode = 0; parNode < _nParPoints; parNode++)
922921
{
922+
// binding: every bound state might depend on every bound and liquid state
923+
// reaction: every liquid and bound state might depend on every liquid and bound state
923924
for (int state = 0; state < strideParPoint(); state++)
924925
{
925926
for (int stateDep = 0; stateDep < strideParPoint(); stateDep++) {
@@ -931,14 +932,14 @@ namespace parts
931932
}
932933
}
933934
}
935+
934936
/**
935937
* @brief calculates the particle dispersion jacobian Pattern, including entries for the dependence of particle entries on bulk entries through film diffusion boundary condition
936938
* @detail Does NOT add film diffusion entries for the dependence of bulk conc. on particle conc.
937939
*/
938940
void ParticleDiffusionOperatorDG::setParticleJacobianPattern(std::vector<ParticleDiffusionOperatorDG::T>& tripletList, unsigned int offsetPar, unsigned int offsetBulk, unsigned int colNode, unsigned int secIdx)
939941
{
940-
941-
parBindingPattern(tripletList, offsetPar, colNode);
942+
parBindingAndReactionPattern(tripletList, offsetPar, colNode);
942943

943944
// Ordering of particle surface diffusion:
944945
// 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);

‎test/data/model_COL1D_GRM_dynLin_1comp_sensbenchmark2.json‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -79,7 +79,6 @@
7979
0.0
8080
],
8181
"ADSORPTION_MODEL": "LINEAR",
82-
"ADSORPTION_MODEL_MULTIPLEX": 1,
8382
"NBOUND": [
8483
1
8584
],

0 commit comments

Comments
 (0)