diff --git a/src/libcadet/model/LumpedRateModelWithPoresDG2D.cpp b/src/libcadet/model/LumpedRateModelWithPoresDG2D.cpp index ecd522a26..d77b50904 100644 --- a/src/libcadet/model/LumpedRateModelWithPoresDG2D.cpp +++ b/src/libcadet/model/LumpedRateModelWithPoresDG2D.cpp @@ -449,6 +449,8 @@ bool LumpedRateModelWithPoresDG2D::configureModelDiscretization(IParameterProvid _disc.axNPoints = _convDispOp.axNPoints(); _disc.radNPoints = _convDispOp.radNPoints(); + _disc.radNNodes = _convDispOp.radNNodes(); + _disc.radNElem = _convDispOp.radNElem(); _disc.nBulkPoints = _disc.axNPoints * _disc.radNPoints; // Precompute offsets of particle type DOFs @@ -626,11 +628,11 @@ bool LumpedRateModelWithPoresDG2D::configure(IParameterProvider& paramProvider) // Let PAR_TYPE_VOLFRAC default to 1.0 for backwards compatibility if (paramProvider.exists("PAR_TYPE_VOLFRAC")) - _parTypeVolFracMode = readAndRegisterMultiplexParam(paramProvider, _parameters, _parTypeVolFrac, "PAR_TYPE_VOLFRAC", _disc.axNPoints, _disc.radNPoints, _disc.nParType, _unitOpIdx); + _parTypeVolFracMode = readAndRegisterMultiplexParam(paramProvider, _parameters, _parTypeVolFrac, "PAR_TYPE_VOLFRAC", _disc.axNPoints, _disc.radNElem, _disc.nParType, _unitOpIdx); else { // Only one particle type present - _parTypeVolFrac.resize(_disc.axNPoints * _disc.radNPoints, 1.0); + _parTypeVolFrac.resize(_disc.axNPoints * _disc.radNElem, 1.0); _parTypeVolFracMode = MultiplexMode::Independent; } @@ -638,18 +640,18 @@ bool LumpedRateModelWithPoresDG2D::configure(IParameterProvider& paramProvider) // Check whether all sizes are matched if (_disc.nParType != _parRadius.size()) throw InvalidParameterException("Number of elements in field PAR_RADIUS does not match number of particle types"); - if (_disc.nParType * _disc.nBulkPoints != _parTypeVolFrac.size()) - throw InvalidParameterException("Number of elements in field PAR_TYPE_VOLFRAC does not match number of particle types"); + if (_disc.nParType * _disc.axNPoints * _disc.radNElem != _parTypeVolFrac.size()) + throw InvalidParameterException("Number of elements in field PAR_TYPE_VOLFRAC does not match number of particle types times RAD_NELEM * (AX_POLYDEG+1)*AX_NELEM"); if (_disc.nParType != _parPorosity.size()) throw InvalidParameterException("Number of elements in field PAR_POROSITY does not match number of particle types"); // Check that particle volume fractions sum to 1.0 - for (unsigned int i = 0; i < _disc.axNPoints * _disc.radNPoints; ++i) + for (unsigned int i = 0; i < _disc.axNPoints * _disc.radNElem; ++i) { const double volFracSum = std::accumulate(_parTypeVolFrac.begin() + i * _disc.nParType, _parTypeVolFrac.begin() + (i+1) * _disc.nParType, 0.0, [](double a, const active& b) -> double { return a + static_cast(b); }); if (std::abs(1.0 - volFracSum) > 1e-10) - throw InvalidParameterException("Sum of field PAR_TYPE_VOLFRAC differs from 1.0 (is " + std::to_string(volFracSum) + ") in axial cell " + std::to_string(i / _disc.radNPoints) + " radial cell " + std::to_string(i % _disc.radNPoints)); + throw InvalidParameterException("Sum of field PAR_TYPE_VOLFRAC differs from 1.0 (is " + std::to_string(volFracSum) + ") in axial element " + std::to_string(i / _disc.radNElem) + " radial element " + std::to_string(i % _disc.radNElem)); } // Read vectorial parameters (which may also be section dependent; transport) @@ -671,7 +673,7 @@ bool LumpedRateModelWithPoresDG2D::configure(IParameterProvider& paramProvider) // Register initial conditions parameters registerParam1DArray(_parameters, _initC, [=](bool multi, unsigned int comp) { return makeParamId(hashString("INIT_C"), _unitOpIdx, comp, ParTypeIndep, BoundStateIndep, ReactionIndep, SectionIndep); }); - if (_disc.radNPoints > 1) + if (_disc.radNElem > 1) registerParam2DArray(_parameters, _initC, [=](bool multi, unsigned int rad, unsigned int comp) { return makeParamId(hashString("INIT_C"), _unitOpIdx, comp, ParTypeIndep, BoundStateIndep, rad, SectionIndep); }, _disc.nComp); if (_singleBinding) @@ -679,9 +681,9 @@ bool LumpedRateModelWithPoresDG2D::configure(IParameterProvider& paramProvider) for (unsigned int c = 0; c < _disc.nComp; ++c) _parameters[makeParamId(hashString("INIT_CP"), _unitOpIdx, c, ParTypeIndep, BoundStateIndep, ReactionIndep, SectionIndep)] = &_initCp[c]; - if (_disc.radNPoints > 1) + if (_disc.radNElem > 1) { - for (unsigned int r = 0; r < _disc.radNPoints; ++r) + for (unsigned int r = 0; r < _disc.radNElem; ++r) { for (unsigned int c = 0; c < _disc.nComp; ++c) _parameters[makeParamId(hashString("INIT_CP"), _unitOpIdx, c, ParTypeIndep, BoundStateIndep, r, SectionIndep)] = &_initCp[r * _disc.nComp * _disc.nParType + c]; @@ -691,7 +693,7 @@ bool LumpedRateModelWithPoresDG2D::configure(IParameterProvider& paramProvider) else { registerParam2DArray(_parameters, _initCp, [=](bool multi, unsigned int type, unsigned int comp) { return makeParamId(hashString("INIT_CP"), _unitOpIdx, comp, type, BoundStateIndep, ReactionIndep, SectionIndep); }, _disc.nComp); - if (_disc.radNPoints > 1) + if (_disc.radNElem > 1) registerParam3DArray(_parameters, _initCp, [=](bool multi, unsigned int rad, unsigned int type, unsigned int comp) { return makeParamId(hashString("INIT_CP"), _unitOpIdx, comp, type, BoundStateIndep, rad, SectionIndep); }, _disc.nComp, _disc.nParType); } @@ -722,13 +724,13 @@ bool LumpedRateModelWithPoresDG2D::configure(IParameterProvider& paramProvider) } // Register radially dependent - if (_disc.radNPoints > 1) + if (_disc.radNElem > 1) { if (_singleBinding) { _binding[0]->fillBoundPhaseInitialParameters(initParams.data(), _unitOpIdx, ParTypeIndep); - for (unsigned int r = 0; r < _disc.radNPoints; ++r) + for (unsigned int r = 0; r < _disc.radNElem; ++r) { for (ParameterId& pId : initParams) pId.reaction = r; @@ -740,7 +742,7 @@ bool LumpedRateModelWithPoresDG2D::configure(IParameterProvider& paramProvider) } else { - for (unsigned int r = 0; r < _disc.radNPoints; ++r) + for (unsigned int r = 0; r < _disc.radNElem; ++r) { for (unsigned int type = 0; type < _disc.nParType; ++type) { @@ -1303,7 +1305,8 @@ int LumpedRateModelWithPoresDG2D::residualFlux(double t, unsigned int secIdx, St for (unsigned int j = 0; j < _disc.radNPoints; ++j) { - invBetaC = 1.0 / static_cast(_convDispOp.columnPorosity(j)) - 1.0; + const int radZone = std::floor(j / _disc.radNNodes); + invBetaC = 1.0 / static_cast(_convDispOp.columnPorosity(radZone)) - 1.0; jacCF_val = invBetaC * _parGeomSurfToVol[type] / radius; jacPF_val = -_parGeomSurfToVol[type] / (epsP * radius); @@ -1319,7 +1322,7 @@ int LumpedRateModelWithPoresDG2D::residualFlux(double t, unsigned int secIdx, St if (wantRes) { // Add flux to column void / bulk volume equations - resCol[ClIdx] += jacCF_val * static_cast(filmDiff[comp]) * static_cast(_parTypeVolFrac[type + i * _disc.nParType * _disc.radNPoints + j * _disc.nParType]) * (yCol[ClIdx] - yCol[CpIdx]); + resCol[ClIdx] += jacCF_val * static_cast(filmDiff[comp]) * static_cast(_parTypeVolFrac[type + i * _disc.nParType * _disc.radNElem + radZone * _disc.nParType]) * (yCol[ClIdx] - yCol[CpIdx]); // Add flux to particle / bead volume equations resCol[CpIdx] += jacPF_val / static_cast(poreAccFactor[comp]) * static_cast(filmDiff[comp]) * (yCol[ClIdx] - yCol[CpIdx]); @@ -1330,11 +1333,11 @@ int LumpedRateModelWithPoresDG2D::residualFlux(double t, unsigned int secIdx, St // add Cl on Cl entries (added since entries already set in transport jacobian) // row: already at current bulk node and component // col: already at current bulk node and component - _globalJac.coeffRef(ClIdx, ClIdx) += static_cast(jacCF_val) * static_cast(filmDiff[comp]) * static_cast(_parTypeVolFrac[type + i * _disc.nParType * _disc.radNPoints + j * _disc.nParType]); + _globalJac.coeffRef(ClIdx, ClIdx) += static_cast(jacCF_val) * static_cast(filmDiff[comp]) * static_cast(_parTypeVolFrac[type + i * _disc.nParType * _disc.radNElem + radZone * _disc.nParType]); // add Cl on Cp entries // row: already at current bulk node and component // col: jump to particle phase - _globalJac.coeffRef(ClIdx, CpIdx) = -static_cast(jacCF_val) * static_cast(filmDiff[comp]) * static_cast(_parTypeVolFrac[type + i * _disc.nParType * _disc.radNPoints + j * _disc.nParType]); + _globalJac.coeffRef(ClIdx, CpIdx) = -static_cast(jacCF_val) * static_cast(filmDiff[comp]) * static_cast(_parTypeVolFrac[type + i * _disc.nParType * _disc.radNElem + radZone * _disc.nParType]); // add Cp on Cp entries (added since entries already set in particle jacobian) // row: already at particle. already at current node and liquid state @@ -1526,10 +1529,10 @@ unsigned int LumpedRateModelWithPoresDG2D::localOutletComponentIndex(unsigned in { // Inlets are duplicated so need to be accounted for if (static_cast(_convDispOp.currentVelocity(port)) >= 0.0) - // Forward Flow: outlet is last cell + // Forward Flow: outlet is last node return _disc.nComp * _disc.radNPoints + (_disc.axNPoints - 1) * _disc.nComp * _disc.radNPoints + port * _disc.nComp; else - // Backward flow: Outlet is first cell + // Backward flow: Outlet is first node return _disc.nComp * _disc.radNPoints + _disc.nComp * port; } @@ -1558,7 +1561,7 @@ bool LumpedRateModelWithPoresDG2D::setParameter(const ParameterId& pId, double v { if (pId.unitOperation == _unitOpIdx) { - if (multiplexParameterValue(pId, hashString("PAR_TYPE_VOLFRAC"), _parTypeVolFracMode, _parTypeVolFrac, _disc.axNPoints, _disc.radNPoints, _disc.nParType, value, nullptr)) + if (multiplexParameterValue(pId, hashString("PAR_TYPE_VOLFRAC"), _parTypeVolFracMode, _parTypeVolFrac, _disc.axNPoints, _disc.radNElem, _disc.nParType, value, nullptr)) return true; if (multiplexCompTypeSecParameterValue(pId, hashString("PORE_ACCESSIBILITY"), _poreAccessFactorMode, _poreAccessFactor, _disc.nParType, _disc.nComp, value, nullptr)) return true; @@ -1588,7 +1591,7 @@ void LumpedRateModelWithPoresDG2D::setSensitiveParameterValue(const ParameterId& { if (pId.unitOperation == _unitOpIdx) { - if (multiplexParameterValue(pId, hashString("PAR_TYPE_VOLFRAC"), _parTypeVolFracMode, _parTypeVolFrac, _disc.axNPoints, _disc.radNPoints, _disc.nParType, value, &_sensParams)) + if (multiplexParameterValue(pId, hashString("PAR_TYPE_VOLFRAC"), _parTypeVolFracMode, _parTypeVolFrac, _disc.axNPoints, _disc.radNElem, _disc.nParType, value, &_sensParams)) return; if (multiplexCompTypeSecParameterValue(pId, hashString("PORE_ACCESSIBILITY"), _poreAccessFactorMode, _poreAccessFactor, _disc.nParType, _disc.nComp, value, &_sensParams)) return; @@ -1613,7 +1616,7 @@ bool LumpedRateModelWithPoresDG2D::setSensitiveParameter(const ParameterId& pId, { if (pId.unitOperation == _unitOpIdx) { - if (multiplexParameterAD(pId, hashString("PAR_TYPE_VOLFRAC"), _parTypeVolFracMode, _parTypeVolFrac, _disc.axNPoints, _disc.radNPoints, _disc.nParType, adDirection, adValue, _sensParams)) + if (multiplexParameterAD(pId, hashString("PAR_TYPE_VOLFRAC"), _parTypeVolFracMode, _parTypeVolFrac, _disc.axNPoints, _disc.radNElem, _disc.nParType, adDirection, adValue, _sensParams)) { LOG(Debug) << "Found parameter " << pId << ": Dir " << adDirection << " is set to " << adValue; return true; @@ -1728,8 +1731,8 @@ int LumpedRateModelWithPoresDG2D::Exporter::writeParticleMobilePhase(unsigned in int LumpedRateModelWithPoresDG2D::Exporter::writeInlet(unsigned int port, double* buffer) const { - cadet_assert(port < _disc.radNPoints); - std::copy_n(_data + port * _disc.nComp, _disc.nComp, buffer); + cadet_assert(port < _disc.radNElem); + std::copy_n(_data + port * _disc.nComp * _disc.radNNodes, _disc.nComp * _disc.radNNodes, buffer); return _disc.nComp; } @@ -1741,12 +1744,12 @@ int LumpedRateModelWithPoresDG2D::Exporter::writeInlet(double* buffer) const int LumpedRateModelWithPoresDG2D::Exporter::writeOutlet(unsigned int port, double* buffer) const { - cadet_assert(port < _disc.radNPoints); + cadet_assert(port < _disc.radNElem); if (_model._convDispOp.currentVelocity(port) >= 0) - std::copy_n(&_idx.c(_data, _disc.axNPoints - 1, port, 0), _disc.nComp, buffer); + std::copy_n(&_idx.c(_data, _disc.axNPoints - 1, port * _disc.radNNodes, 0), _disc.nComp * _disc.radNNodes, buffer); else - std::copy_n(&_idx.c(_data, 0, port, 0), _disc.nComp, buffer); + std::copy_n(&_idx.c(_data, 0, port * _disc.radNNodes, 0), _disc.nComp * _disc.radNNodes, buffer); return _disc.nComp; } diff --git a/src/libcadet/model/LumpedRateModelWithPoresDG2D.hpp b/src/libcadet/model/LumpedRateModelWithPoresDG2D.hpp index 6775d3331..f3b7958be 100644 --- a/src/libcadet/model/LumpedRateModelWithPoresDG2D.hpp +++ b/src/libcadet/model/LumpedRateModelWithPoresDG2D.hpp @@ -78,8 +78,8 @@ class LumpedRateModelWithPoresDG2D : public UnitOperationBase virtual UnitOpIdx unitOperationId() const CADET_NOEXCEPT { return _unitOpIdx; } virtual unsigned int numComponents() const CADET_NOEXCEPT { return _disc.nComp; } virtual void setFlowRates(active const* in, active const* out) CADET_NOEXCEPT; - virtual unsigned int numInletPorts() const CADET_NOEXCEPT { return _disc.radNPoints; } - virtual unsigned int numOutletPorts() const CADET_NOEXCEPT { return _disc.radNPoints; } + virtual unsigned int numInletPorts() const CADET_NOEXCEPT { return _disc.radNElem; } + virtual unsigned int numOutletPorts() const CADET_NOEXCEPT { return _disc.radNElem; } virtual bool canAccumulate() const CADET_NOEXCEPT { return false; } static const char* identifier() { return "LUMPED_RATE_MODEL_WITH_PORES_2D"; } @@ -230,6 +230,8 @@ class LumpedRateModelWithPoresDG2D : public UnitOperationBase unsigned int nComp; //!< Number of components unsigned int axNPoints; //!< Number of axial discrete points unsigned int radNPoints; //!< Number of radial discrete points + unsigned int radNElem; //!< Number of radial DG elements, which can be used as radial zones for radial parameter inhomogeneities + unsigned int radNNodes; //!< Number of radial DG nodes per element unsigned int nBulkPoints; //!< Number of total bulk discrete points unsigned int nParType; //!< Number of particle types unsigned int* parTypeOffset; //!< Array with offsets (in particle block) to particle type, additional last element contains total number of particle DOFs @@ -375,8 +377,8 @@ class LumpedRateModelWithPoresDG2D : public UnitOperationBase virtual unsigned int numComponents() const CADET_NOEXCEPT { return _disc.nComp; } virtual unsigned int numPrimaryCoordinates() const CADET_NOEXCEPT { return _disc.axNPoints; } virtual unsigned int numSecondaryCoordinates() const CADET_NOEXCEPT { return _disc.radNPoints; } - virtual unsigned int numInletPorts() const CADET_NOEXCEPT { return _disc.radNPoints; } - virtual unsigned int numOutletPorts() const CADET_NOEXCEPT { return _disc.radNPoints; } + virtual unsigned int numInletPorts() const CADET_NOEXCEPT { return _disc.radNElem; } + virtual unsigned int numOutletPorts() const CADET_NOEXCEPT { return _disc.radNElem; } virtual unsigned int numParticleTypes() const CADET_NOEXCEPT { return _disc.nParType; } virtual unsigned int numParticleShells(unsigned int parType) const CADET_NOEXCEPT { return 1; } virtual unsigned int numBoundStates(unsigned int parType) const CADET_NOEXCEPT { return _disc.strideBound[parType]; } diff --git a/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperatorDG.cpp b/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperatorDG.cpp index 590ffed60..bf343e921 100644 --- a/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperatorDG.cpp +++ b/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperatorDG.cpp @@ -30,10 +30,10 @@ using namespace Eigen; namespace { -cadet::model::MultiplexMode readAndRegisterMultiplexParam(cadet::IParameterProvider& paramProvider, std::unordered_map& parameters, std::vector& values, const std::string& name, unsigned int nComp, unsigned int radNPoints, cadet::UnitOpIdx uoi) +cadet::model::MultiplexMode readAndRegisterMultiplexParam(cadet::IParameterProvider& paramProvider, std::unordered_map& parameters, std::vector& values, const std::string& name, unsigned int nComp, unsigned int radNElem, cadet::UnitOpIdx uoi) { cadet::model::MultiplexMode mode = cadet::model::MultiplexMode::Independent; - readParameterMatrix(values, paramProvider, name, nComp * radNPoints, 1); + readParameterMatrix(values, paramProvider, name, nComp * radNElem, 1); unsigned int nSec = 1; if (paramProvider.exists(name + "_MULTIPLEX")) { @@ -47,8 +47,8 @@ cadet::model::MultiplexMode readAndRegisterMultiplexParam(cadet::IParameterProvi else if (modeConfig == 1) { mode = cadet::model::MultiplexMode::Radial; - if (values.size() != radNPoints) - throw cadet::InvalidParameterException("Number of elements in field " + name + " inconsistent with " + name + "_MULTIPLEX (should be " + std::to_string(radNPoints) + ")"); + if (values.size() != radNElem) + throw cadet::InvalidParameterException("Number of elements in field " + name + " inconsistent with " + name + "_MULTIPLEX (should be " + std::to_string(radNElem) + ")"); } else if (modeConfig == 2) { @@ -59,8 +59,8 @@ cadet::model::MultiplexMode readAndRegisterMultiplexParam(cadet::IParameterProvi else if (modeConfig == 3) { mode = cadet::model::MultiplexMode::ComponentRadial; - if (values.size() != nComp * radNPoints) - throw cadet::InvalidParameterException("Number of elements in field " + name + " inconsistent with " + name + "_MULTIPLEX (should be " + std::to_string(nComp * radNPoints) + ")"); + if (values.size() != nComp * radNElem) + throw cadet::InvalidParameterException("Number of elements in field " + name + " inconsistent with " + name + "_MULTIPLEX (should be " + std::to_string(nComp * radNElem) + ")"); } else if (modeConfig == 4) { @@ -70,10 +70,10 @@ cadet::model::MultiplexMode readAndRegisterMultiplexParam(cadet::IParameterProvi else if (modeConfig == 5) { mode = cadet::model::MultiplexMode::RadialSection; - if (values.size() % radNPoints != 0) - throw cadet::InvalidParameterException("Number of elements in field " + name + " is not a positive multiple of radNPoints (" + std::to_string(radNPoints) + ")"); + if (values.size() % radNElem != 0) + throw cadet::InvalidParameterException("Number of elements in field " + name + " is not a positive multiple of radNElem (" + std::to_string(radNElem) + ")"); - nSec = values.size() / radNPoints; + nSec = values.size() / radNElem; } else if (modeConfig == 6) { @@ -86,10 +86,10 @@ cadet::model::MultiplexMode readAndRegisterMultiplexParam(cadet::IParameterProvi else if (modeConfig == 7) { mode = cadet::model::MultiplexMode::ComponentRadialSection; - if (values.size() % (nComp * radNPoints) != 0) - throw cadet::InvalidParameterException("Number of elements in field " + name + " is not a positive multiple of NCOMP * radNPoints (" + std::to_string(nComp * radNPoints) + ")"); + if (values.size() % (nComp * radNElem) != 0) + throw cadet::InvalidParameterException("Number of elements in field " + name + " is not a positive multiple of NCOMP * radNElem (" + std::to_string(nComp * radNElem) + ")"); - nSec = values.size() / (nComp * radNPoints); + nSec = values.size() / (nComp * radNElem); } } else @@ -98,24 +98,24 @@ cadet::model::MultiplexMode readAndRegisterMultiplexParam(cadet::IParameterProvi mode = cadet::model::MultiplexMode::Independent; else if (values.size() == nComp) mode = cadet::model::MultiplexMode::Component; - else if (values.size() == radNPoints) + else if (values.size() == radNElem) mode = cadet::model::MultiplexMode::Radial; - else if (values.size() == radNPoints * nComp) + else if (values.size() == radNElem * nComp) mode = cadet::model::MultiplexMode::ComponentRadial; else if (values.size() % nComp == 0) { mode = cadet::model::MultiplexMode::ComponentSection; nSec = values.size() / nComp; } - else if (values.size() % radNPoints == 0) + else if (values.size() % radNElem == 0) { mode = cadet::model::MultiplexMode::RadialSection; - nSec = values.size() / radNPoints; + nSec = values.size() / radNElem; } - else if (values.size() % (radNPoints * nComp) == 0) + else if (values.size() % (radNElem * nComp) == 0) { mode = cadet::model::MultiplexMode::ComponentRadialSection; - nSec = values.size() / (nComp * radNPoints); + nSec = values.size() / (nComp * radNElem); } else throw cadet::InvalidParameterException("Could not infer multiplex mode of field " + name + ", set " + name + "_MULTIPLEX or change number of elements"); @@ -129,24 +129,24 @@ cadet::model::MultiplexMode readAndRegisterMultiplexParam(cadet::IParameterProvi case cadet::model::MultiplexMode::Independent: case cadet::model::MultiplexMode::Section: { - std::vector p(nComp * radNPoints * nSec); + std::vector p(nComp * radNElem * nSec); for (unsigned int s = 0; s < nSec; ++s) - std::fill(p.begin() + s * radNPoints * nComp, p.begin() + (s+1) * radNPoints * nComp, values[s]); + std::fill(p.begin() + s * radNElem * nComp, p.begin() + (s+1) * radNElem * nComp, values[s]); values = std::move(p); for (unsigned int s = 0; s < nSec; ++s) - parameters[cadet::makeParamId(nameHash, uoi, cadet::CompIndep, cadet::ParTypeIndep, cadet::BoundStateIndep, cadet::ReactionIndep, (mode == cadet::model::MultiplexMode::Independent) ? cadet::SectionIndep : s)] = &values[s * radNPoints * nComp]; + parameters[cadet::makeParamId(nameHash, uoi, cadet::CompIndep, cadet::ParTypeIndep, cadet::BoundStateIndep, cadet::ReactionIndep, (mode == cadet::model::MultiplexMode::Independent) ? cadet::SectionIndep : s)] = &values[s * radNElem * nComp]; } break; case cadet::model::MultiplexMode::Component: case cadet::model::MultiplexMode::ComponentSection: { - std::vector p(nComp * radNPoints * nSec); + std::vector p(nComp * radNElem * nSec); for (unsigned int s = 0; s < nSec; ++s) { for (unsigned int i = 0; i < nComp; ++i) - std::copy(values.begin() + s * nComp, values.begin() + (s+1) * nComp, p.begin() + i * nComp + s * nComp * radNPoints); + std::copy(values.begin() + s * nComp, values.begin() + (s+1) * nComp, p.begin() + i * nComp + s * nComp * radNElem); } values = std::move(p); @@ -154,29 +154,29 @@ cadet::model::MultiplexMode readAndRegisterMultiplexParam(cadet::IParameterProvi for (unsigned int s = 0; s < nSec; ++s) { for (unsigned int i = 0; i < nComp; ++i) - parameters[cadet::makeParamId(nameHash, uoi, i, cadet::ParTypeIndep, cadet::BoundStateIndep, cadet::ReactionIndep, (mode == cadet::model::MultiplexMode::Component) ? cadet::SectionIndep : s)] = &values[s * radNPoints * nComp + i]; + parameters[cadet::makeParamId(nameHash, uoi, i, cadet::ParTypeIndep, cadet::BoundStateIndep, cadet::ReactionIndep, (mode == cadet::model::MultiplexMode::Component) ? cadet::SectionIndep : s)] = &values[s * radNElem * nComp + i]; } } break; case cadet::model::MultiplexMode::Radial: case cadet::model::MultiplexMode::RadialSection: { - std::vector p(nComp * radNPoints * nSec); - for (unsigned int i = 0; i < radNPoints * nSec; ++i) + std::vector p(nComp * radNElem * nSec); + for (unsigned int i = 0; i < radNElem * nSec; ++i) std::fill(p.begin() + i * nComp, p.begin() + (i+1) * nComp, values[i]); values = std::move(p); for (unsigned int s = 0; s < nSec; ++s) { - for (unsigned int i = 0; i < radNPoints; ++i) - parameters[cadet::makeParamId(nameHash, uoi, cadet::CompIndep, i, cadet::BoundStateIndep, cadet::ReactionIndep, (mode == cadet::model::MultiplexMode::Radial) ? cadet::SectionIndep : s)] = &values[s * radNPoints * nComp + i * nComp]; + for (unsigned int i = 0; i < radNElem; ++i) + parameters[cadet::makeParamId(nameHash, uoi, cadet::CompIndep, i, cadet::BoundStateIndep, cadet::ReactionIndep, (mode == cadet::model::MultiplexMode::Radial) ? cadet::SectionIndep : s)] = &values[s * radNElem * nComp + i * nComp]; } } break; case cadet::model::MultiplexMode::ComponentRadial: case cadet::model::MultiplexMode::ComponentRadialSection: - cadet::registerParam3DArray(parameters, values, [=](bool multi, unsigned int sec, unsigned int compartment, unsigned int comp) { return cadet::makeParamId(nameHash, uoi, comp, compartment, cadet::BoundStateIndep, cadet::ReactionIndep, multi ? sec : cadet::SectionIndep); }, nComp, radNPoints); + cadet::registerParam3DArray(parameters, values, [=](bool multi, unsigned int sec, unsigned int compartment, unsigned int comp) { return cadet::makeParamId(nameHash, uoi, comp, compartment, cadet::BoundStateIndep, cadet::ReactionIndep, multi ? sec : cadet::SectionIndep); }, nComp, radNElem); break; case cadet::model::MultiplexMode::Axial: case cadet::model::MultiplexMode::AxialRadial: @@ -190,7 +190,7 @@ cadet::model::MultiplexMode readAndRegisterMultiplexParam(cadet::IParameterProvi return mode; } -bool multiplexParameterValue(const cadet::ParameterId& pId, cadet::StringHash nameHash, cadet::model::MultiplexMode mode, std::vector& data, unsigned int nComp, unsigned int radNPoints, double value, std::unordered_set const* sensParams) +bool multiplexParameterValue(const cadet::ParameterId& pId, cadet::StringHash nameHash, cadet::model::MultiplexMode mode, std::vector& data, unsigned int nComp, unsigned int radNElem, double value, std::unordered_set const* sensParams) { if (pId.name != nameHash) return false; @@ -217,11 +217,11 @@ bool multiplexParameterValue(const cadet::ParameterId& pId, cadet::StringHash na || (pId.reaction != cadet::ReactionIndep) || (pId.section == cadet::SectionIndep)) return false; - if (sensParams && !cadet::contains(*sensParams, &data[pId.section * nComp * radNPoints])) + if (sensParams && !cadet::contains(*sensParams, &data[pId.section * nComp * radNElem])) return false; - for (unsigned int i = 0; i < nComp * radNPoints; ++i) - data[i + pId.section * nComp * radNPoints].setValue(value); + for (unsigned int i = 0; i < nComp * radNElem; ++i) + data[i + pId.section * nComp * radNElem].setValue(value); return true; } @@ -234,7 +234,7 @@ bool multiplexParameterValue(const cadet::ParameterId& pId, cadet::StringHash na if (sensParams && !cadet::contains(*sensParams, &data[pId.component])) return false; - for (unsigned int i = 0; i < radNPoints; ++i) + for (unsigned int i = 0; i < radNElem; ++i) data[i * nComp + pId.component].setValue(value); return true; @@ -245,11 +245,11 @@ bool multiplexParameterValue(const cadet::ParameterId& pId, cadet::StringHash na || (pId.reaction != cadet::ReactionIndep) || (pId.section == cadet::SectionIndep)) return false; - if (sensParams && !cadet::contains(*sensParams, &data[pId.component + pId.section * nComp * radNPoints])) + if (sensParams && !cadet::contains(*sensParams, &data[pId.component + pId.section * nComp * radNElem])) return false; - for (unsigned int i = 0; i < radNPoints; ++i) - data[i * nComp + pId.component + pId.section * nComp * radNPoints].setValue(value); + for (unsigned int i = 0; i < radNElem; ++i) + data[i * nComp + pId.component + pId.section * nComp * radNElem].setValue(value); return true; } @@ -273,11 +273,11 @@ bool multiplexParameterValue(const cadet::ParameterId& pId, cadet::StringHash na || (pId.reaction != cadet::ReactionIndep) || (pId.section == cadet::SectionIndep)) return false; - if (sensParams && !cadet::contains(*sensParams, &data[pId.particleType * nComp + pId.section * nComp * radNPoints])) + if (sensParams && !cadet::contains(*sensParams, &data[pId.particleType * nComp + pId.section * nComp * radNElem])) return false; for (unsigned int i = 0; i < nComp; ++i) - data[i + pId.particleType * nComp + pId.section * nComp * radNPoints].setValue(value); + data[i + pId.particleType * nComp + pId.section * nComp * radNElem].setValue(value); return true; } @@ -300,10 +300,10 @@ bool multiplexParameterValue(const cadet::ParameterId& pId, cadet::StringHash na || (pId.reaction != cadet::ReactionIndep) || (pId.section == cadet::SectionIndep)) return false; - if (sensParams && !cadet::contains(*sensParams, &data[pId.component + pId.particleType * nComp + pId.section * nComp * radNPoints])) + if (sensParams && !cadet::contains(*sensParams, &data[pId.component + pId.particleType * nComp + pId.section * nComp * radNElem])) return false; - data[pId.component + pId.particleType * nComp + pId.section * nComp * radNPoints].setValue(value); + data[pId.component + pId.particleType * nComp + pId.section * nComp * radNElem].setValue(value); return true; } @@ -319,7 +319,7 @@ bool multiplexParameterValue(const cadet::ParameterId& pId, cadet::StringHash na return false; } -bool multiplexParameterAD(const cadet::ParameterId& pId, cadet::StringHash nameHash, cadet::model::MultiplexMode mode, std::vector& data, unsigned int nComp, unsigned int radNPoints, unsigned int adDirection, double adValue, std::unordered_set& sensParams) +bool multiplexParameterAD(const cadet::ParameterId& pId, cadet::StringHash nameHash, cadet::model::MultiplexMode mode, std::vector& data, unsigned int nComp, unsigned int radNElem, unsigned int adDirection, double adValue, std::unordered_set& sensParams) { if (pId.name != nameHash) return false; @@ -345,10 +345,10 @@ bool multiplexParameterAD(const cadet::ParameterId& pId, cadet::StringHash nameH || (pId.reaction != cadet::ReactionIndep) || (pId.section == cadet::SectionIndep)) return false; - sensParams.insert(&data[pId.section * nComp * radNPoints]); + sensParams.insert(&data[pId.section * nComp * radNElem]); - for (unsigned int i = 0; i < nComp * radNPoints; ++i) - data[i + pId.section * nComp * radNPoints].setADValue(adDirection, adValue); + for (unsigned int i = 0; i < nComp * radNElem; ++i) + data[i + pId.section * nComp * radNElem].setADValue(adDirection, adValue); return true; } @@ -360,7 +360,7 @@ bool multiplexParameterAD(const cadet::ParameterId& pId, cadet::StringHash nameH sensParams.insert(&data[pId.component]); - for (unsigned int i = 0; i < radNPoints; ++i) + for (unsigned int i = 0; i < radNElem; ++i) data[i * nComp + pId.component].setADValue(adDirection, adValue); return true; @@ -371,10 +371,10 @@ bool multiplexParameterAD(const cadet::ParameterId& pId, cadet::StringHash nameH || (pId.reaction != cadet::ReactionIndep) || (pId.section == cadet::SectionIndep)) return false; - sensParams.insert(&data[pId.component + pId.section * nComp * radNPoints]); + sensParams.insert(&data[pId.component + pId.section * nComp * radNElem]); - for (unsigned int i = 0; i < radNPoints; ++i) - data[i * nComp + pId.component + pId.section * nComp * radNPoints].setADValue(adDirection, adValue); + for (unsigned int i = 0; i < radNElem; ++i) + data[i * nComp + pId.component + pId.section * nComp * radNElem].setADValue(adDirection, adValue); return true; } @@ -397,10 +397,10 @@ bool multiplexParameterAD(const cadet::ParameterId& pId, cadet::StringHash nameH || (pId.reaction != cadet::ReactionIndep) || (pId.section == cadet::SectionIndep)) return false; - sensParams.insert(&data[pId.particleType * nComp + pId.section * nComp * radNPoints]); + sensParams.insert(&data[pId.particleType * nComp + pId.section * nComp * radNElem]); for (unsigned int i = 0; i < nComp; ++i) - data[i + pId.particleType * nComp + pId.section * nComp * radNPoints].setADValue(adDirection, adValue); + data[i + pId.particleType * nComp + pId.section * nComp * radNElem].setADValue(adDirection, adValue); return true; } @@ -422,9 +422,9 @@ bool multiplexParameterAD(const cadet::ParameterId& pId, cadet::StringHash nameH || (pId.reaction != cadet::ReactionIndep) || (pId.section == cadet::SectionIndep)) return false; - sensParams.insert(&data[pId.component + pId.particleType * nComp + pId.section * nComp * radNPoints]); + sensParams.insert(&data[pId.component + pId.particleType * nComp + pId.section * nComp * radNElem]); - data[pId.component + pId.particleType * nComp + pId.section * nComp * radNPoints].setADValue(adDirection, adValue); + data[pId.component + pId.particleType * nComp + pId.section * nComp * radNElem].setADValue(adDirection, adValue); return true; } @@ -455,8 +455,7 @@ namespace parts * @brief Creates a TwoDimensionalConvectionDispersionOperatorDG */ TwoDimensionalConvectionDispersionOperatorDG::TwoDimensionalConvectionDispersionOperatorDG() : _colPorosities(0), _dir(0), _dispersionDep(nullptr), -_radLiftMCyl(nullptr), _transMrCyl(nullptr), _invTransMrCyl(nullptr), _transTildeMr(nullptr), _transTildeMrDash(nullptr), _transTildeSrDash(nullptr), -_SrCyl(nullptr), _jacConvection(nullptr), _jacAxDispersion(nullptr), _jacRadDispersion(nullptr) +_radLiftMCyl(nullptr), _transMrCyl(nullptr), _invTransMrCyl(nullptr), _SrCyl(nullptr), _jacConvection(nullptr), _jacAxDispersion(nullptr), _jacRadDispersion(nullptr) { } @@ -466,9 +465,6 @@ TwoDimensionalConvectionDispersionOperatorDG::~TwoDimensionalConvectionDispersio delete[] _radLiftMCyl; delete[] _transMrCyl; delete[] _invTransMrCyl; - delete[] _transTildeMr; - delete[] _transTildeMrDash; - delete[] _transTildeSrDash; delete[] _SrCyl; delete[] _jacConvection; delete[] _jacAxDispersion; @@ -605,9 +601,6 @@ void TwoDimensionalConvectionDispersionOperatorDG::initializeDG() { _transMrCyl = new MatrixXd[_radNElem]; _invTransMrCyl = new MatrixXd[_radNElem]; - _transTildeMr = new MatrixXd[_radNElem]; - _transTildeMrDash = new MatrixXd[_radNElem]; - _transTildeSrDash = new MatrixXd[_radNElem]; _SrCyl = new MatrixXd[_radNElem]; _radLiftMCyl = new MatrixXd[_radNElem]; // Jacobian blocks @@ -716,11 +709,10 @@ bool TwoDimensionalConvectionDispersionOperatorDG::configureModelDiscretization( _radNodeStride = radNodeStride; _radElemStride = _radNNodes * _radNodeStride; - //_radialCoordinates.resize(_radNPoints + 1); // todo not needed, delete? _radialElemInterfaces.resize(_radNElem + 1); _radDelta.resize(_radNElem); - _nodalCrossSections.resize(_radNPoints); - _curVelocity.resize(_radNPoints); + _elementCrossSections.resize(_radNElem); + _curVelocity.resize(_radNElem); return true; } @@ -741,12 +733,12 @@ bool TwoDimensionalConvectionDispersionOperatorDG::configure(UnitOpIdx unitOpIdx readScalarParameterOrArray(_colPorosities, paramProvider, "COL_POROSITY", 1); _axDelta = _colLength / _axNElem; // deltaR is treated in updateRadialDisc - if ((_colPorosities.size() != 1) && (_colPorosities.size() != _radNPoints)) - throw InvalidParameterException("Number of elements in field COL_POROSITY is neither 1 nor radNPoints (" + std::to_string(_radNPoints) + ")"); + if ((_colPorosities.size() != 1) && (_colPorosities.size() != _radNElem)) + throw InvalidParameterException("Number of elements in field COL_POROSITY is neither 1 nor radNElem (" + std::to_string(_radNElem) + ")"); _singlePorosity = (_colPorosities.size() == 1); if (_singlePorosity) - _colPorosities = std::vector(_radNPoints, _colPorosities[0]); + _colPorosities = std::vector(_radNElem, _colPorosities[0]); // Read radial discretization mode and default to "EQUIDISTANT" paramProvider.pushScope("discretization"); @@ -791,26 +783,26 @@ bool TwoDimensionalConvectionDispersionOperatorDG::configure(UnitOpIdx unitOpIdx // Rad-dep, sec-dep _singleVelocity = false; - if (!_singleVelocity && (_velocity.size() % _radNPoints != 0)) - throw InvalidParameterException("Number of elements in field VELOCITY is not a positive multiple of radNPoints (" + std::to_string(_radNPoints) + ")"); + if (!_singleVelocity && (_velocity.size() % _radNElem != 0)) + throw InvalidParameterException("Number of elements in field VELOCITY is not a positive multiple of radNElem (" + std::to_string(_radNElem) + ")"); if ((mode == 0) && (_velocity.size() != 1)) throw InvalidParameterException("Number of elements in field VELOCITY inconsistent with VELOCITY_MULTIPLEX (should be 1)"); - if ((mode == 1) && (_velocity.size() != _radNPoints)) - throw InvalidParameterException("Number of elements in field VELOCITY inconsistent with VELOCITY_MULTIPLEX (should be " + std::to_string(_radNPoints) + ")"); + if ((mode == 1) && (_velocity.size() != _radNElem)) + throw InvalidParameterException("Number of elements in field VELOCITY inconsistent with VELOCITY_MULTIPLEX (should be " + std::to_string(_radNElem) + ")"); } else { // Infer radial dependence of VELOCITY: - // size not divisible by radNPoints -> radial independent - _singleVelocity = ((_velocity.size() % _radNPoints) != 0); + // size not divisible by _radNElem -> radial independent + _singleVelocity = ((_velocity.size() % _radNElem) != 0); } // Expand _velocity to make it component dependent if (_singleVelocity) { - std::vector expanded(_velocity.size() * _radNPoints); + std::vector expanded(_velocity.size() * _radNElem); for (std::size_t i = 0; i < _velocity.size(); ++i) - std::fill(expanded.begin() + i * _radNPoints, expanded.begin() + (i + 1) * _radNPoints, _velocity[i]); + std::fill(expanded.begin() + i * _radNElem, expanded.begin() + (i + 1) * _radNElem, _velocity[i]); _velocity = std::move(expanded); } @@ -818,17 +810,17 @@ bool TwoDimensionalConvectionDispersionOperatorDG::configure(UnitOpIdx unitOpIdx else { _singleVelocity = false; - _velocity.resize(_radNPoints, 1.0); + _velocity.resize(_radNElem, 1.0); } // Register VELOCITY if (_singleVelocity) { - if (_velocity.size() > _radNPoints) + if (_velocity.size() > _radNElem) { // Register only the first item in each section - for (std::size_t i = 0; i < _velocity.size() / _radNPoints; ++i) - parameters[makeParamId(hashString("VELOCITY"), unitOpIdx, CompIndep, ParTypeIndep, BoundStateIndep, ReactionIndep, i)] = &_velocity[i * _radNPoints]; + for (std::size_t i = 0; i < _velocity.size() / _radNElem; ++i) + parameters[makeParamId(hashString("VELOCITY"), unitOpIdx, CompIndep, ParTypeIndep, BoundStateIndep, ReactionIndep, i)] = &_velocity[i * _radNElem]; } else { @@ -837,12 +829,12 @@ bool TwoDimensionalConvectionDispersionOperatorDG::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); }, _radNPoints); + registerParam2DArray(parameters, _velocity, [=](bool multi, unsigned int sec, unsigned int compartment) { return makeParamId(hashString("VELOCITY"), unitOpIdx, CompIndep, compartment, BoundStateIndep, ReactionIndep, multi ? sec : SectionIndep); }, _radNElem); - _dir = std::vector(_radNPoints, 1); + _dir = std::vector(_radNElem, 1); - _axialDispersionMode = readAndRegisterMultiplexParam(paramProvider, parameters, _axialDispersion, "COL_DISPERSION", _nComp, _radNPoints, unitOpIdx); - _radialDispersionMode = readAndRegisterMultiplexParam(paramProvider, parameters, _radialDispersion, "COL_DISPERSION_RADIAL", _nComp, _radNPoints, unitOpIdx); + _axialDispersionMode = readAndRegisterMultiplexParam(paramProvider, parameters, _axialDispersion, "COL_DISPERSION", _nComp, _radNElem, unitOpIdx); + _radialDispersionMode = readAndRegisterMultiplexParam(paramProvider, parameters, _radialDispersion, "COL_DISPERSION_RADIAL", _nComp, _radNElem, unitOpIdx); // Add parameters to map parameters[makeParamId(hashString("COL_LENGTH"), unitOpIdx, CompIndep, ParTypeIndep, BoundStateIndep, ReactionIndep, SectionIndep)] = &_colLength; @@ -852,24 +844,11 @@ bool TwoDimensionalConvectionDispersionOperatorDG::configure(UnitOpIdx unitOpIdx // compute standard DG operators initializeDG(); - // interpolate radial dispersion to quadrature nodes - _curRadialDispersionTilde = std::vector(_radNElem * _qNNodes, 0.0); - _curAxialDispersionTilde = std::vector(_radNElem * _qNNodes, 0.0); - // todo component dependence! getSectionDependentSlice gives dispersion parameter in radial position major - const active* const curRadialDispersion = getSectionDependentSlice(_radialDispersion, _radNPoints * _nComp, 0); - const active* const curAxialDispersion = getSectionDependentSlice(_axialDispersion, _radNPoints * _nComp, 0); - for (unsigned int i = 0; i < _radInterpolationM.rows(); i++) { - for (unsigned int j = 0; j < _radInterpolationM.cols(); j++) { - _curRadialDispersionTilde[i] += _radInterpolationM(i, j) * curRadialDispersion[j]; - _curAxialDispersionTilde[i] += _radInterpolationM(i, j) * curAxialDispersion[j]; - } - } - // compute nodal cross section areas and radial element interfaces updateRadialDisc(); // compute DG operators that depend on radial geometry and dispersion after updateRadialDisc - const active* const d_rad = getSectionDependentSlice(_radialDispersion, _radNPoints * _nComp, 0); + const active* const d_rad = getSectionDependentSlice(_radialDispersion, _radNElem * _nComp, 0); // todo component dependence of radial dispersion const int comp = 0; for (unsigned int rElem = 0; rElem < _radNElem; rElem++) @@ -882,10 +861,6 @@ bool TwoDimensionalConvectionDispersionOperatorDG::configure(UnitOpIdx unitOpIdx _transMrCyl[rElem] = (static_cast(_radialElemInterfaces[rElem]) * dgtoolbox::mMatrix(_radPolyDeg, _radNodes, 0.0, 0.0) + static_cast(_radDelta[rElem]) / 2.0 * dgtoolbox::mMatrix(_radPolyDeg, _radNodes, 0.0, 1.0)).transpose(); _invTransMrCyl[rElem] = _transMrCyl[rElem].inverse(); - _transTildeMr[rElem] = calcTildeMr(rElem, &_curAxialDispersionTilde[0]).transpose(); // todo still needed in new derivation? - _transTildeMrDash[rElem] = calcTildeMrDash(rElem, &_curAxialDispersionTilde[0]).transpose(); // todo still needed in new derivation? - _transTildeSrDash[rElem] = calcTildeSrDash(rElem, &_curRadialDispersionTilde[0]).transpose(); // todo still needed in new derivation? - _SrCyl[rElem] = _transMrCyl[rElem].transpose() * dgtoolbox::derivativeMatrix(_radPolyDeg, _radNodes); } @@ -907,20 +882,14 @@ bool TwoDimensionalConvectionDispersionOperatorDG::notifyDiscontinuousSectionTra bool hasChanged = false; // todo update operators for section dependent parameters - const active* const curRadialDispersion = getSectionDependentSlice(_radialDispersion, _radNPoints * _nComp, 0); - for (unsigned int i = 0; i < _radInterpolationM.rows(); i++) { - for (unsigned int j = 0; j < _radInterpolationM.cols(); j++) { - _curRadialDispersionTilde[i] += _radInterpolationM(i, j) * curRadialDispersion[j]; - } - } 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 - active const* const dirNew = getSectionDependentSlice(_velocity, _radNPoints, secIdx); + active const* const dirNew = getSectionDependentSlice(_velocity, _radNElem, secIdx); - for (unsigned int i = 0; i < _radNPoints; ++i) + for (unsigned int i = 0; i < _radNElem; ++i) { const int newDir = (dirNew[i] >= 0) ? 1 : -1; if (_dir[i] * newDir < 0) @@ -950,13 +919,13 @@ bool TwoDimensionalConvectionDispersionOperatorDG::notifyDiscontinuousSectionTra */ void TwoDimensionalConvectionDispersionOperatorDG::setFlowRates(int compartment, const active& in, const active& out) CADET_NOEXCEPT { - _curVelocity[compartment] = _dir[compartment] * in / (_nodalCrossSections[compartment] * _colPorosities[compartment]); + _curVelocity[compartment] = _dir[compartment] * in / (_elementCrossSections[compartment] * _colPorosities[compartment]); } void TwoDimensionalConvectionDispersionOperatorDG::setFlowRates(active const* in, active const* out) CADET_NOEXCEPT { - for (unsigned int compartment = 0; compartment < _radNPoints; ++compartment) - _curVelocity[compartment] = in[compartment] / (_nodalCrossSections[compartment] * _colPorosities[compartment]); + for (unsigned int compartment = 0; compartment < _radNElem; ++compartment) + _curVelocity[compartment] = in[compartment] / (_elementCrossSections[compartment] * _colPorosities[compartment]); } double TwoDimensionalConvectionDispersionOperatorDG::inletFactor(unsigned int idxSec, int idxRad) const CADET_NOEXCEPT // todo is this function needed? @@ -967,12 +936,12 @@ double TwoDimensionalConvectionDispersionOperatorDG::inletFactor(unsigned int id const active& TwoDimensionalConvectionDispersionOperatorDG::axialDispersion(unsigned int idxSec, int idxRad, int idxComp) const CADET_NOEXCEPT { - return *(getSectionDependentSlice(_axialDispersion, _radNPoints * _nComp, idxSec) + idxRad * _nComp + idxComp); + return *(getSectionDependentSlice(_axialDispersion, _radNElem * _nComp, idxSec) + idxRad * _nComp + idxComp); } const active& TwoDimensionalConvectionDispersionOperatorDG::radialDispersion(unsigned int idxSec, int idxRad, int idxComp) const CADET_NOEXCEPT { - return *(getSectionDependentSlice(_radialDispersion, _radNPoints * _nComp, idxSec) + idxRad * _nComp + idxComp); + return *(getSectionDependentSlice(_radialDispersion, _radNElem * _nComp, idxSec) + idxRad * _nComp + idxComp); } /** @@ -1015,8 +984,8 @@ int TwoDimensionalConvectionDispersionOperatorDG::residualImpl(const IModel& mod const int auxTildeAxNodeStride = _radNElem * _qNNodes; const int auxTildeAxElemStride = _axNNodes * auxTildeAxNodeStride; - const active* const curAxialDispersion = getSectionDependentSlice(_axialDispersion, _radNPoints * _nComp, 0); - const active* const curRadialDispersion = getSectionDependentSlice(_radialDispersion, _radNPoints * _nComp, 0); + const active* const curAxialDispersion = getSectionDependentSlice(_axialDispersion, _radNElem * _nComp, secIdx); + const active* const curRadialDispersion = getSectionDependentSlice(_radialDispersion, _radNElem * _nComp, secIdx); for (unsigned int comp = 0; comp < _nComp; comp++) { @@ -1200,18 +1169,6 @@ int TwoDimensionalConvectionDispersionOperatorDG::residualImpl(const IModel& mod ); // Axial dispersion - - // TODO old code for radial position dependent parameters : delete or fix - //// Step 1: To get _GzTilde * _transTildeMrDash * _invTransMrCyl, we need to multiply a matrix A of dimensionality N x M with M x M blocks of a matrix B of dimensionality M*K x M. - //// Additionally Matrix A might have actives while B has doubles and the resulting matrix is immediately multiplied by another double matrix C with dimensionality N x M. - //// We thus first multiply B and C, then cast the result to A's scalar type and block-wise multiply by the result of B * C which again has dimensionality M*K x M. - //MatrixMap matrixCache(reinterpret_cast(&_matrixProductCache[0]), _axNNodes, _radNNodes, Stride(1, _radNNodes)); - //matrixCache.setZero(); - // - //for (int qNode = 0; qNode < _qNNodes; qNode++) - //{ - // matrixCache += _GzTilde.template cast() * (_transTildeMrDash[rEidx].block(qNode * _radNNodes, 0, _radNNodes, _radNNodes) * _invTransMrCyl[rEidx]).template cast(); - //} _Res -= 2.0 / static_cast(_axDelta) * _axInvMM.template cast() * ( // surface integral @@ -1443,7 +1400,7 @@ bool TwoDimensionalConvectionDispersionOperatorDG::computeConvDispJacobianBlocks }); /* axial dispersion block */ - const active* const curAxialDispersion = getSectionDependentSlice(_axialDispersion, _radNPoints * _nComp, 0); + const active* const curAxialDispersion = getSectionDependentSlice(_axialDispersion, _radNElem * _nComp, 0); const double axDisp = static_cast(curAxialDispersion[0]); MatrixXd gStarZ = MatrixXd::Zero(Np, 5 * Np); @@ -1487,7 +1444,7 @@ bool TwoDimensionalConvectionDispersionOperatorDG::computeConvDispJacobianBlocks } /* radial dispersion block */ - const active* const curRadialDispersion = getSectionDependentSlice(_radialDispersion, _radNPoints * _nComp, 0); + const active* const curRadialDispersion = getSectionDependentSlice(_radialDispersion, _radNElem * _nComp, 0); const double radDisp = static_cast(curRadialDispersion[0]); MatrixXd gStarR = MatrixXd::Zero(Np, 5 * Np); @@ -1855,14 +1812,7 @@ void TwoDimensionalConvectionDispersionOperatorDG::setEquidistantRadialDisc() else _radialElemInterfaces[r + 1] = h * (r + 1); - subcellLeftEnd = dgtoolbox::mapRefToPhys(_radDelta, r, -1.0); - VectorXd radWeights = _radInvWeights.cwiseInverse(); - for (unsigned int node = 0; node < _radNNodes; ++node) - { - subcellRightEnd = subcellLeftEnd + radWeights[node] / 2.0 * _radDelta[r]; - _nodalCrossSections[r * _radNNodes + node] = pi * (pow(subcellRightEnd, 2.0) - pow(subcellLeftEnd, 2.0)); - subcellLeftEnd = subcellRightEnd; - } + _elementCrossSections[r] = pi * (pow(_radialElemInterfaces[r + 1], 2.0) - pow(_radialElemInterfaces[r], 2.0)); } } @@ -1884,14 +1834,7 @@ void TwoDimensionalConvectionDispersionOperatorDG::setEquivolumeRadialDisc() _radDelta[r] = _radialElemInterfaces[r + 1] - _radialElemInterfaces[r]; - subcellLeftEnd = dgtoolbox::mapRefToPhys(_radDelta, r, -1.0); - VectorXd radWeights = _radInvWeights.cwiseInverse(); - for (unsigned int node = 0; node < _radNNodes; ++node) - { - subcellRightEnd = subcellLeftEnd + radWeights[node] / 2.0 * _radDelta[r]; - _nodalCrossSections[r * _radNNodes + node] = pi * (pow(subcellRightEnd, 2.0) - pow(subcellLeftEnd, 2.0)); - subcellLeftEnd = subcellRightEnd; - } + _elementCrossSections[r] = pi * (pow(_radialElemInterfaces[r + 1], 2.0) - pow(_radialElemInterfaces[r], 2.0)); } } @@ -1905,14 +1848,7 @@ void TwoDimensionalConvectionDispersionOperatorDG::setUserdefinedRadialDisc() { _radDelta[r] = _radialElemInterfaces[r + 1] - _radialElemInterfaces[r]; - subcellLeftEnd = dgtoolbox::mapRefToPhys(_radDelta, r, -1.0); - VectorXd radWeights = _radInvWeights.cwiseInverse(); - for (unsigned int node = 0; node < _radNNodes; ++node) - { - subcellRightEnd = subcellLeftEnd + radWeights[node] / 2.0 * _radDelta[r]; - _nodalCrossSections[r * _radNNodes + node] = pi * (pow(subcellRightEnd, 2.0) - pow(subcellLeftEnd, 2.0)); - subcellLeftEnd = subcellRightEnd; - } + _elementCrossSections[r] = pi * (pow(_radialElemInterfaces[r + 1], 2.0) - pow(_radialElemInterfaces[r], 2.0)); } } @@ -1931,21 +1867,21 @@ bool TwoDimensionalConvectionDispersionOperatorDG::setParameter(const ParameterI if (_singlePorosity && (pId.name == hashString("COL_POROSITY")) && (pId.component == CompIndep) && (pId.boundState == BoundStateIndep) && (pId.reaction == ReactionIndep) && (pId.section == SectionIndep) && (pId.particleType == ParTypeIndep)) { - for (unsigned int i = 0; i < _radNPoints; ++i) + for (unsigned int i = 0; i < _radNElem; ++i) _colPorosities[i].setValue(value); return true; } if (_singleVelocity && (pId.name == hashString("VELOCITY")) && (pId.component == CompIndep) && (pId.boundState == BoundStateIndep) && (pId.reaction == ReactionIndep)) { - if (_velocity.size() > _radNPoints) + if (_velocity.size() > _radNElem) { // Section dependent if (pId.section == SectionIndep) return false; - for (unsigned int i = 0; i < _radNPoints; ++i) - _velocity[pId.section * _radNPoints + i].setValue(value); + for (unsigned int i = 0; i < _radNElem; ++i) + _velocity[pId.section * _radNElem + i].setValue(value); } else { @@ -1953,16 +1889,16 @@ bool TwoDimensionalConvectionDispersionOperatorDG::setParameter(const ParameterI if (pId.section != SectionIndep) return false; - for (unsigned int i = 0; i < _radNPoints; ++i) + for (unsigned int i = 0; i < _radNElem; ++i) _velocity[i].setValue(value); } } - const bool ad = multiplexParameterValue(pId, hashString("COL_DISPERSION"), _axialDispersionMode, _axialDispersion, _nComp, _radNPoints, value, nullptr); + const bool ad = multiplexParameterValue(pId, hashString("COL_DISPERSION"), _axialDispersionMode, _axialDispersion, _nComp, _radNElem, value, nullptr); if (ad) return true; - const bool adr = multiplexParameterValue(pId, hashString("COL_DISPERSION_RADIAL"), _radialDispersionMode, _radialDispersion, _nComp, _radNPoints, value, nullptr); + const bool adr = multiplexParameterValue(pId, hashString("COL_DISPERSION_RADIAL"), _radialDispersionMode, _radialDispersion, _nComp, _radNElem, value, nullptr); if (adr) return true; @@ -1976,7 +1912,7 @@ bool TwoDimensionalConvectionDispersionOperatorDG::setSensitiveParameterValue(co { if (contains(sensParams, &_colPorosities[0])) { - for (unsigned int i = 0; i < _radNPoints; ++i) + for (unsigned int i = 0; i < _radNElem; ++i) _colPorosities[i].setValue(value); return true; } @@ -1984,14 +1920,14 @@ bool TwoDimensionalConvectionDispersionOperatorDG::setSensitiveParameterValue(co if (_singleVelocity && (pId.name == hashString("VELOCITY")) && (pId.component == CompIndep) && (pId.boundState == BoundStateIndep) && (pId.reaction == ReactionIndep)) { - if (_velocity.size() > _radNPoints) + if (_velocity.size() > _radNElem) { // Section dependent - if ((pId.section == SectionIndep) || !contains(sensParams, &_velocity[pId.section * _radNPoints])) + if ((pId.section == SectionIndep) || !contains(sensParams, &_velocity[pId.section * _radNElem])) return false; - for (unsigned int i = 0; i < _radNPoints; ++i) - _velocity[pId.section * _radNPoints + i].setValue(value); + for (unsigned int i = 0; i < _radNElem; ++i) + _velocity[pId.section * _radNElem + i].setValue(value); } else { @@ -1999,16 +1935,16 @@ bool TwoDimensionalConvectionDispersionOperatorDG::setSensitiveParameterValue(co if ((pId.section != SectionIndep) || !contains(sensParams, &_velocity[0])) return false; - for (unsigned int i = 0; i < _radNPoints; ++i) + for (unsigned int i = 0; i < _radNElem; ++i) _velocity[i].setValue(value); } } - const bool ad = multiplexParameterValue(pId, hashString("COL_DISPERSION"), _axialDispersionMode, _axialDispersion, _nComp, _radNPoints, value, &sensParams); + const bool ad = multiplexParameterValue(pId, hashString("COL_DISPERSION"), _axialDispersionMode, _axialDispersion, _nComp, _radNElem, value, &sensParams); if (ad) return true; - const bool adr = multiplexParameterValue(pId, hashString("COL_DISPERSION_RADIAL"), _radialDispersionMode, _radialDispersion, _nComp, _radNPoints, value, &sensParams); + const bool adr = multiplexParameterValue(pId, hashString("COL_DISPERSION_RADIAL"), _radialDispersionMode, _radialDispersion, _nComp, _radNElem, value, &sensParams); if (adr) return true; @@ -2021,7 +1957,7 @@ bool TwoDimensionalConvectionDispersionOperatorDG::setSensitiveParameter(std::un && (pId.reaction == ReactionIndep) && (pId.section == SectionIndep) && (pId.particleType == ParTypeIndep)) { sensParams.insert(&_colPorosities[0]); - for (unsigned int i = 0; i < _radNPoints; ++i) + for (unsigned int i = 0; i < _radNElem; ++i) _colPorosities[i].setADValue(adDirection, adValue); return true; @@ -2029,15 +1965,15 @@ bool TwoDimensionalConvectionDispersionOperatorDG::setSensitiveParameter(std::un if (_singleVelocity && (pId.name == hashString("VELOCITY")) && (pId.component == CompIndep) && (pId.boundState == BoundStateIndep) && (pId.reaction == ReactionIndep)) { - if (_velocity.size() > _radNPoints) + if (_velocity.size() > _radNElem) { // Section dependent if (pId.section == SectionIndep) return false; - sensParams.insert(&_velocity[pId.section * _radNPoints]); - for (unsigned int i = 0; i < _radNPoints; ++i) - _velocity[pId.section * _radNPoints + i].setADValue(adDirection, adValue); + sensParams.insert(&_velocity[pId.section * _radNElem]); + for (unsigned int i = 0; i < _radNElem; ++i) + _velocity[pId.section * _radNElem + i].setADValue(adDirection, adValue); } else { @@ -2046,16 +1982,16 @@ bool TwoDimensionalConvectionDispersionOperatorDG::setSensitiveParameter(std::un return false; sensParams.insert(&_velocity[0]); - for (unsigned int i = 0; i < _radNPoints; ++i) + for (unsigned int i = 0; i < _radNElem; ++i) _velocity[i].setADValue(adDirection, adValue); } } - const bool ad = multiplexParameterAD(pId, hashString("COL_DISPERSION"), _axialDispersionMode, _axialDispersion, _nComp, _radNPoints, adDirection, adValue, sensParams); + const bool ad = multiplexParameterAD(pId, hashString("COL_DISPERSION"), _axialDispersionMode, _axialDispersion, _nComp, _radNElem, adDirection, adValue, sensParams); if (ad) return true; - const bool adr = multiplexParameterAD(pId, hashString("COL_DISPERSION_RADIAL"), _radialDispersionMode, _radialDispersion, _nComp, _radNPoints, adDirection, adValue, sensParams); + const bool adr = multiplexParameterAD(pId, hashString("COL_DISPERSION_RADIAL"), _radialDispersionMode, _radialDispersion, _nComp, _radNElem, adDirection, adValue, sensParams); if (adr) return true; diff --git a/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperatorDG.hpp b/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperatorDG.hpp index cc3ab671f..77c5aa9ad 100644 --- a/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperatorDG.hpp +++ b/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperatorDG.hpp @@ -114,6 +114,8 @@ class TwoDimensionalConvectionDispersionOperatorDG inline const int radNNodes() const { return _radNNodes; } inline const int radNElem() const { return _radNElem; } inline const int elemNPoints() const { return _elemNPoints; } + inline const unsigned int axNPoints() const { return _axNPoints; } + inline const unsigned int radNPoints() const { return _radNPoints; } double relativeAxialCoordinate(unsigned int idx) const { @@ -128,8 +130,6 @@ class TwoDimensionalConvectionDispersionOperatorDG return (elem * static_cast(_radDelta[elem]) + 0.5 * static_cast(_radDelta[elem]) * (1.0 + _radNodes[idx % _radNNodes])) / static_cast(_colRadius); } - inline const unsigned int axNPoints() const CADET_NOEXCEPT { return _axNPoints; } - inline const unsigned int radNPoints() const CADET_NOEXCEPT { return _radNPoints; } inline bool isCurrentFlowForward(int idx) const CADET_NOEXCEPT { return _curVelocity[idx] >= 0.0; } const active& axialDispersion(unsigned int idxSec, int idxRad, int idxComp) const CADET_NOEXCEPT; const active& radialDispersion(unsigned int idxSec, int idxRad, int idxComp) const CADET_NOEXCEPT; @@ -197,9 +197,9 @@ class TwoDimensionalConvectionDispersionOperatorDG unsigned int _radNElem; //!< Number of radial elements unsigned int _radPolyDeg; //!< Polynomial degree of radial discretization unsigned int _radNNodes; //!< Number of radial discrete points - unsigned int _quadratureRule; //!< Numerical quadrature rule - unsigned int _quadratureOrder; //!< Order of the numerical quadrature - unsigned int _qNNodes; //!< Number of quadrature nodes + unsigned int _quadratureRule; //!< Numerical quadrature rule // todo needed? + unsigned int _quadratureOrder; //!< Order of the numerical quadrature // todo needed? + unsigned int _qNNodes; //!< Number of quadrature nodes // todo needed? unsigned int _elemNPoints; //!< Number of discrete points per 2D element unsigned int _bulkNPoints; //!< Number of total 2D grid points (bulk grid) // strides @@ -229,9 +229,6 @@ class TwoDimensionalConvectionDispersionOperatorDG Eigen::MatrixXd* _transMrCyl; //!< Radial transposed mass matrix with cylinder metrics Eigen::MatrixXd* _invTransMrCyl; //!< Radial inverted transposed mass matrix with cylinder metrics Eigen::MatrixXd _radInterpolationM; //!< Polynomial interpolation matrix from (radial) LGL nodes to quadrature nodes - Eigen::MatrixXd* _transTildeMr; //!< Main eq. transposed mass matrix adjusted for cylindrical metrics and dispersion - Eigen::MatrixXd* _transTildeMrDash; //!< Main eq. transposed mass matrix on quadrature nodes adjusted for cylindrical metrics and dispersion - Eigen::MatrixXd* _transTildeSrDash; //!< Main eq. transposed stiffness matrix on quadrature nodes adjusted for cylindrical metrics and dispersion Eigen::MatrixXd* _SrCyl; //!< Main eq. stiffness matrix with cylindrical metrics // Jacobian blocks Eigen::MatrixXd* _jacConvection; //!< Convection Jacobian blocks for each radial element @@ -260,7 +257,7 @@ class TwoDimensionalConvectionDispersionOperatorDG active _axDelta; //!< Axial equidistant element spacing std::vector _radDelta; //!< Radial element spacing std::vector _radialElemInterfaces; //!< Coordinates of the element interfaces - std::vector _nodalCrossSections; //!< cross section area for each node + std::vector _elementCrossSections; //!< cross section area for each radial element RadialDiscretizationMode _radialDiscretizationMode; @@ -270,8 +267,6 @@ class TwoDimensionalConvectionDispersionOperatorDG std::vector _axialDispersion; //!< Axial dispersion coefficient \f$ D_{\text{ax}} \f$ MultiplexMode _axialDispersionMode; //!< Multiplex mode of the axial dispersion std::vector _radialDispersion; //!< Radial dispersion coefficient \f$ D_{\rho} \f$ at interpolation nodes - std::vector _curRadialDispersionTilde; //!< Radial dispersion coefficient \f$ D_{\rho} \f$ at radial quadrature nodes - std::vector _curAxialDispersionTilde; //!< Axial dispersion coefficient \f$ D_{\rho} \f$ at radial quadrature nodes MultiplexMode _radialDispersionMode; //!< Multiplex mode of the radial dispersion std::vector _velocity; //!< Interstitial velocity parameter std::vector _curVelocity; //!< Current interstitial velocity \f$ u \f$