From 4e623c4d2e5b8a7da3b30c3495f46924c3753644 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Wed, 31 Jul 2024 09:59:44 +0200 Subject: [PATCH 01/17] Add PhaseTransitionModel.hpp and first changes for Linear Exchange --- src/libcadet/model/ModelUtils.hpp | 11 + src/libcadet/model/PhaseTransitionModel.hpp | 359 ++++++++++ .../model/exchange/LinearExchange.cpp | 642 ++++++++++++++++++ 3 files changed, 1012 insertions(+) create mode 100644 src/libcadet/model/PhaseTransitionModel.hpp create mode 100644 src/libcadet/model/exchange/LinearExchange.cpp diff --git a/src/libcadet/model/ModelUtils.hpp b/src/libcadet/model/ModelUtils.hpp index 5c5079f2e..fdad98b21 100644 --- a/src/libcadet/model/ModelUtils.hpp +++ b/src/libcadet/model/ModelUtils.hpp @@ -82,6 +82,17 @@ inline unsigned int numBoundStates(unsigned int const* const nBound, unsigned in return totalBound; } +/** + * @brief Returns the number of bound states + * @param [in] nBound Array with number of bound states for each component + * @param [in] nComp Number of components + * @return Number of bound states + */ +inline unsigned int numExchangeStates(unsigned int nChannel, unsigned int nComp) +{ + return nChannel * nComp; +} + /** * @brief Returns the number of binding components (i.e., components that have at least @c 1 bound state) * @param [in] nBound Array with number of bound states for each component diff --git a/src/libcadet/model/PhaseTransitionModel.hpp b/src/libcadet/model/PhaseTransitionModel.hpp new file mode 100644 index 000000000..e2981d525 --- /dev/null +++ b/src/libcadet/model/PhaseTransitionModel.hpp @@ -0,0 +1,359 @@ +// ============================================================================= +// CADET +// +// Copyright © 2008-2024: The CADET Authors +// Please see the AUTHORS and CONTRIBUTORS file. +// +// All rights reserved. This program and the accompanying materials +// are made available under the terms of the GNU Public License v3.0 (or, at +// your option, any later version) which accompanies this distribution, and +// is available at http://www.gnu.org/licenses/gpl.html +// ============================================================================= + +/** + * @file + * Defines the BindingModel interface. + */ + +#ifndef LIBCADET_PHASETRANSITIONMODELINTERFACE_HPP_ +#define LIBCADET_PHASETRANSITIONMODELINTERFACE_HPP_ + +#include + +#include "CompileTimeConfig.hpp" +#include "cadet/ParameterProvider.hpp" +#include "cadet/ParameterId.hpp" +#include "linalg/DenseMatrix.hpp" +#include "linalg/BandMatrix.hpp" + +#ifdef ENABLE_DG + #include "linalg/BandedEigenSparseRowIterator.hpp" +#endif + +#include "AutoDiff.hpp" +#include "SimulationTypes.hpp" +#include "Memory.hpp" + +namespace cadet +{ + +class IParameterProvider; +class IExternalFunction; + +struct ColumnPosition; + +namespace model +{ + +class IPhaseTransitionModel +{ +public: + + virtual ~IPhaseTransitionModel() CADET_NOEXCEPT { } + + /** + * @brief Returns the name of the binding model + * @details This name is also used to identify and create the binding model in the factory. + * @return Name of the binding model + */ + virtual const char* name() const CADET_NOEXCEPT = 0; + + /** + * @brief Returns whether the binding model requires additional parameters supplied by configure() + * @details After construction of an IPhaseTransitionModel object, configureModelDiscretization() is called. + * The binding model may require to read additional parameters from the adsorption group + * of the parameter provider. This opportunity is given by a call to configure(). + * However, a binding model may not require this. This function communicates to the outside + * whether a call to configure() is necessary. + * @return @c true if configure() has to be called, otherwise @c false + */ + virtual bool requiresConfiguration() const CADET_NOEXCEPT = 0; + + /** + * @brief Returns whether the dynamic reaction model uses the IParameterProvider in configureModelDiscretization() + * @details If the IParameterProvider is used in configureModelDiscretization(), it has to be in the correct scope. + * @return @c true if the IParameterProvider is used in configureModelDiscretization(), otherwise @c false + */ + virtual bool usesParamProviderInDiscretizationConfig() const CADET_NOEXCEPT = 0; + + /** + * @brief Sets fixed parameters of the binding model (e.g., the number of components and bound states) + * @details This function is called prior to configure() by the underlying model. + * It can only be called once. Whereas non-structural model parameters + * (e.g., rate constants) are configured by configure(), this function + * sets structural parameters (e.g., number of components and bound states). + * In particular, it determines whether each binding reaction is quasi-stationary. + * + * Before this function is called, usesParamProviderInDiscretizationConfig() + * is queried to determine whether the @p paramProvider is used by this + * function. If it is used, the @p paramProvider is put in the corresponding + * parameter scope. On exit, the @p paramProvider has to point to the same + * scope. + * @param [in] paramProvider Parameter provider + * @param [in] nComp Number of components + * @param [in] nBound Array of size @p nComp which contains the number of bound states for each component + * @param [in] boundOffset Array of size @p nComp with offsets to the first bound state of each component beginning from the solid phase + */ + virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBound, unsigned int const* boundOffset) = 0; + + /** + * @brief Configures the model by extracting all non-structural parameters (e.g., model parameters) from the given @p paramProvider + * @details The scope of the cadet::IParameterProvider is left unchanged on return. + * + * The structure of the model is left unchanged, that is, the number of degrees of + * freedom stays the same (e.g., number of bound states is left unchanged). Only + * true (non-structural) model parameters are read and changed. + * + * This function may only be called if configureModelDiscretization() has been called + * in the past. Contrary to configureModelDiscretization(), it can be called multiple + * times. + * @param [in] paramProvider Parameter provider + * @param [in] unitOpIdx Index of the unit operation this binding model belongs to + * @param [in] parTypeIdx Index of the particle type this binding model belongs to + * @return @c true if the configuration was successful, otherwise @c false + */ + virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) = 0; + + /** + * @brief Returns the ParameterId of all bound phase initial conditions (equations) + * @details The array has to be filled in the order of the equations. + * @param [out] params Array with ParameterId objects to fill + * @param [in] unitOpIdx Index of the unit operation this binding model belongs to + * @param [in] parTypeIdx Index of the particle type this binding model belongs to + */ + virtual void fillBoundPhaseInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) const CADET_NOEXCEPT = 0; + + /** + * @brief Returns the ParameterId of all bound phase initial conditions (equations) + * @details The array has to be filled in the order of the equations. + * @param [out] params Array with ParameterId objects to fill + * @param [in] unitOpIdx Index of the unit operation this binding model belongs to + * @param [in] parTypeIdx Index of the particle type this binding model belongs to + */ + virtual void fillChannelInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) const CADET_NOEXCEPT = 0; + + /** + * @brief Sets external functions for this binding model + * @details The external functions are not owned by this IPhaseTransitionModel. + * + * @param [in] extFuns Pointer to array of IExternalFunction objects of size @p size + * @param [in] size Number of elements in the IExternalFunction array @p extFuns + */ + virtual void setExternalFunctions(IExternalFunction** extFuns, unsigned int size) = 0; + + /** + * @brief Checks whether a given parameter exists + * @param [in] pId ParameterId that identifies the parameter uniquely + * @return @c true if the parameter exists, otherwise @c false + */ + virtual bool hasParameter(const ParameterId& pId) const = 0; + + /** + * @brief Returns all parameters with their current values that can be made sensitive + * @return Map with all parameters that can be made sensitive along with their current value + */ + virtual std::unordered_map getAllParameterValues() const = 0; + + /** + * @brief Sets a parameter value + * @details The parameter identified by its unique parameter is set to the given value. + * + * @param [in] pId ParameterId that identifies the parameter uniquely + * @param [in] value Value of the parameter + * + * @return @c true if the parameter has been successfully set to the given value, + * otherwise @c false (e.g., if the parameter is not available in this model) + */ + virtual bool setParameter(const ParameterId& pId, int value) = 0; + virtual bool setParameter(const ParameterId& pId, double value) = 0; + virtual bool setParameter(const ParameterId& pId, bool value) = 0; + + /** + * @brief Returns a pointer to the parameter identified by the given id + * @param [in] pId Parameter Id of the sensitive parameter + * @return a pointer to the parameter if the binding model contains the parameter, otherwise @c nullptr + */ + virtual active* getParameter(const ParameterId& pId) = 0; + + /** + * @brief Returns whether this binding model has a salt component + * @details The salt component is given by component 0, if present. + * @return @c true if salt is required, otherwise @c false + */ + virtual bool hasSalt() const CADET_NOEXCEPT = 0; + + + /** + * @brief Returns whether this binding model supports multi-state binding + * @return @c true if multi-state binding is supported, otherwise @c false + */ + virtual bool supportsMultistate() const CADET_NOEXCEPT = 0; + + /** + * @brief Returns whether this binding model supports non-binding components + * @details Non-binding components do not have an entry in the solid phase. + * @return @c true if non-binding components are supported, otherwise @c false + */ + virtual bool supportsNonExchange() const CADET_NOEXCEPT = 0; + + /** + * @brief Returns whether this binding model supports non-binding components + * @details Non-binding components do not have an entry in the solid phase. + * @return @c true if non-binding components are supported, otherwise @c false + */ + virtual bool supportsNonBinding() const CADET_NOEXCEPT = 0; + + /** + * @brief Returns whether this binding model has quasi-stationary reaction fluxes + * @return @c true if quasi-stationary fluxes are present, otherwise @c false + */ + virtual bool hasQuasiStationaryReactions() const CADET_NOEXCEPT = 0; + + /** + * @brief Returns whether this binding model has dynamic reaction fluxes + * @return @c true if dynamic fluxes are present, otherwise @c false + */ + virtual bool hasDynamicReactions() const CADET_NOEXCEPT = 0; + + /** + * @brief Returns whether this binding model depends on time + * @details Binding models may depend on time if external functions are used. + * @return @c true if the model is time-dependent, otherwise @c false + */ + virtual bool dependsOnTime() const CADET_NOEXCEPT = 0; + + /** + * @brief Returns whether this binding model requires workspace + * @details The workspace may be required for consistent initialization and / or evaluation + * of residual and Jacobian. A workspace is a memory buffer whose size is given by + * workspaceSize(). + * @return @c true if the model requires a workspace, otherwise @c false + */ + virtual bool requiresWorkspace() const CADET_NOEXCEPT = 0; + + /** + * @brief Returns the size of the required workspace in bytes + * @details The memory is required for externally dependent binding models. + * @param [in] nComp Number of components + * @param [in] totalNumBoundStates Total number of bound states + * @param [in] nBoundStates Array with bound states for each component + * @return Size of the workspace in bytes + */ + virtual unsigned int workspaceSize(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nBoundStates) const CADET_NOEXCEPT = 0; + + /** + * @brief Returns the amount of required AD seed vectors / directions + * @details Only internally required AD directions count (e.g., for Jacobian computation). + * Directions used for parameter sensitivities should not be included here. + * @return The number of required AD seed vectors / directions + */ + virtual unsigned int requiredADdirs() const CADET_NOEXCEPT = 0; + + /** + * @brief Evaluates the fluxes + * @details The binding model is responsible for calculating the flux from the mobile to the solid phase. + * + * This function is called simultaneously from multiple threads. + * It needs to overwrite all values of @p res as the result array @p res is not + * zeroed on entry. + * @param [in] t Current time point + * @param [in] secIdx Index of the current section + * @param [in] colPos Position in normalized coordinates (column inlet = 0, column outlet = 1; outer shell = 1, inner center = 0) + * @param [in] y Pointer to first bound state of the first component in the current particle shell + * @param [in] yCp Pointer to first component of the mobile phase in the current particle shell + * @param [out] res Pointer to result array that is filled with the fluxes + * @param [in,out] workSpace Memory work space + * @return @c 0 on success, @c -1 on non-recoverable error, and @c +1 on recoverable error + */ + virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, active const* yCp, active* res, LinearBufferAllocator workSpace, WithParamSensitivity) const = 0; + virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, active const* yCp, active* res, LinearBufferAllocator workSpace, WithoutParamSensitivity) const = 0; + virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double const* yCp, active* res, LinearBufferAllocator workSpace) const = 0; + virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double const* yCp, double* res, LinearBufferAllocator workSpace) const = 0; + + /** + * @brief Evaluates the Jacobian of the fluxes analytically + * @details This function is called simultaneously from multiple threads. + * It can be left out (empty implementation) if AD is used to evaluate the Jacobian. + * + * The Jacobian matrix is assumed to be zeroed out by the caller. + * @param [in] t Current time point + * @param [in] secIdx Index of the current section + * @param [in] colPos Position in normalized coordinates (column inlet = 0, column outlet = 1; outer shell = 1, inner center = 0) + * @param [in] y Pointer to first bound state of the first component in the current particle shell + * @param [in] offsetCp Offset from the first component of the mobile phase in the current particle shell to @p y + * @param [in,out] jac Row iterator pointing to the first bound states row of the underlying matrix in which the Jacobian is stored + * @param [in,out] workSpace Memory work space + */ + virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const = 0; + virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const = 0; +#ifdef ENABLE_DG + virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; +#endif + /** + * @brief Calculates the time derivative of the quasi-stationary bound state equations + * @details Calculates @f$ \frac{\partial \text{flux}_{\text{qs}}}{\partial t} @f$ for the quasi-stationary equations + * in the flux. + * + * This function is called simultaneously from multiple threads. + * It can be left out (empty implementation) which leads to slightly incorrect initial conditions + * when using externally dependent binding models. + * @param [in] t Current time point + * @param [in] secIdx Index of the current section + * @param [in] colPos Position in normalized coordinates (column inlet = 0, column outlet = 1; outer shell = 1, inner center = 0) + * @param [in] yCp Pointer to first component in the liquid phase of the current particle shell + * @param [in] y Pointer to first bound state of the first component in the current particle shell + * @param [out] dResDt Pointer to array that stores the time derivative + * @param [in,out] workSpace Memory work space + */ + virtual void timeDerivativeQuasiStationaryFluxes(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yCp, double const* y, double* dResDt, LinearBufferAllocator workSpace) const = 0; + + /** + * @brief Returns an array that determines whether each binding reaction is quasi-stationary + * @details Each entry represents a truth value (converts to @c bool) that determines whether + * the corresponding binding reaction is quasi-stationary (@c true) or not (@c false). + * @return Array that determines quasi-stationarity of binding reactions + */ + virtual int const* reactionQuasiStationarity() const CADET_NOEXCEPT = 0; + + /** + * @brief Returns whether a nonlinear solver is required to perform consistent initialization + * @details This function is called before a nonlinear solver attempts to solve the algebraic + * equations. If consistent initialization can be performed using analytic / direct + * formulas (i.e., without running a nonlinear solver), this function would be the + * correct place to do it. + * + * This function is called simultaneously from multiple threads. + * @param [in] t Current time point + * @param [in] secIdx Index of the current section + * @param [in] colPos Position in normalized coordinates (column inlet = 0, column outlet = 1; outer shell = 1, inner center = 0) + * @param [in,out] y Pointer to first bound state of the first component in the current particle shell + * @param [in] yCp Pointer to first component of the mobile phase in the current particle shell + * @param [in,out] workSpace Memory work space + * @return @c true if a nonlinear solver is required for consistent initialization, otherwise @c false + */ + virtual bool preConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const = 0; + + /** + * @brief Called after a nonlinear solver has attempted consistent initialization + * @details In consistent initialization, first preConsistentInitialState() is called. + * Depending on the return value, a nonlinear solver is applied to the algebraic + * equations and, afterwards, this function is called to clean up or refine the + * solution. + * + * This function is called simultaneously from multiple threads. + * @param [in] t Current time point + * @param [in] secIdx Index of the current section + * @param [in] colPos Position in normalized coordinates (column inlet = 0, column outlet = 1; outer shell = 1, inner center = 0) + * @param [in,out] y Pointer to first bound state of the first component in the current particle shell + * @param [in] yCp Pointer to first component of the mobile phase in the current particle shell + * @param [in,out] workSpace Memory work space + */ + virtual void postConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const = 0; + +protected: +}; + +} // namespace model +} // namespace cadet + +#endif // LIBCADET_PHASETRANSITIONMODELINTERFACE_HPP_ diff --git a/src/libcadet/model/exchange/LinearExchange.cpp b/src/libcadet/model/exchange/LinearExchange.cpp new file mode 100644 index 000000000..f05496843 --- /dev/null +++ b/src/libcadet/model/exchange/LinearExchange.cpp @@ -0,0 +1,642 @@ +// ============================================================================= +// CADET +// +// Copyright © 2008-2024: The CADET Authors +// Please see the AUTHORS and CONTRIBUTORS file. +// +// All rights reserved. This program and the accompanying materials +// are made available under the terms of the GNU Public License v3.0 (or, at +// your option, any later version) which accompanies this distribution, and +// is available at http://www.gnu.org/licenses/gpl.html +// ============================================================================= + + + +#include "model/BindingModel.hpp" +#include "model/ExternalFunctionSupport.hpp" +#include "ParamIdUtil.hpp" +#include "model/ModelUtils.hpp" +#include "cadet/Exceptions.hpp" +#include "model/Parameters.hpp" +#include "LocalVector.hpp" +#include "SimulationTypes.hpp" +#include "Memory.hpp" + +#include +#include +#include +#include +#include +#include + +namespace cadet +{ + +namespace model +{ + +/** + * @brief Handles linear exchange model parameters that do not depend on external functions + */ +class LinearParamHandler : public ConstParamHandlerBase +{ +public: + + /** + * @brief Holds actual parameter data + */ + struct ConstParams + { + std::vector kA; //!< Adsorption rate + std::vector kD; //!< Desorption rate + }; + + typedef ConstParams params_t; + typedef ConstParams const* ParamsHandle; + + /** + * @brief Returns name of the binding model + * @return Name of the binding model + */ + static const char* identifier() { return "LINEAR_EX"; } + + LinearParamHandler() CADET_NOEXCEPT : _kA(&_localParams.kA), _kD(&_localParams.kD) { } + + /** + * @brief Reads parameters and verifies them + * @details See IexchangeModel::configure() for details. + * @param [in] paramProvider IParameterProvider used for reading parameters + * @param [in] nComp Number of components + * @param [in] nChannelStates Array with number of bound states for each component + * @return @c true if the parameters were read and validated successfully, otherwise @c false + */ + inline bool configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nChannelStates) + { + _kA.configure("LIN_KA", paramProvider, nComp, nChannelStates); + _kD.configure("LIN_KD", paramProvider, nComp, nChannelStates); + return validateConfig(nComp, nChannelStates); + } + + /** + * @brief Registers all local parameters in a map for further use + * @param [in,out] parameters Map in which the parameters are stored + * @param [in] unitOpIdx Index of the unit operation used for registering the parameters + * @param [in] nComp Number of components + * @param [in] nChannelStates Array with number of bound states for each component + */ + inline void registerParameters(std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nChannelStates) + { + _kA.registerParam("LIN_KA", parameters, unitOpIdx, parTypeIdx, nComp, nChannelStates); + _kD.registerParam("LIN_KD", parameters, unitOpIdx, parTypeIdx, nComp, nChannelStates); + } + + /** + * @brief Reserves space in the storage of the parameters + * @param [in] numElem Number of elements (total) + * @param [in] numSlices Number of slices / exchange site types + * @param [in] nComp Number of components + * @param [in] nChannelStates Array with number of bound states for each component + */ + inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nChannelStates) \ + { + _kA.reserve(numElem, numSlices, nComp, nChannelStates); + _kD.reserve(numElem, numSlices, nComp, nChannelStates); + } + + /** + * @brief Updates the parameter cache in order to take the external profile into account + * @param [in] t Current time + * @param [in] z Axial coordinate in the column + * @param [in] r Radial coordinate in the bead + * @param [in] secIdx Index of the current section + * @param [in] nComp Number of components + * @param [in] nChannelStates Array with number of bound states for each component + * @param [in,out] workSpace Memory buffer for updated data + * @return Externally dependent parameter values + */ + inline ParamsHandle update(double t, unsigned int secIdx, const ColumnPosition& colPos, unsigned int nComp, unsigned int const* nChannelStates, LinearBufferAllocator& workSpace) const + { + return &_localParams; + } + + /** + * @brief Calculates time derivative in case of external dependence + * @param [in] t Current time + * @param [in] z Axial coordinate in the column + * @param [in] r Radial coordinate in the bead + * @param [in] secIdx Index of the current section + * @param [in] nComp Number of components + * @param [in] nChannelStates Array with number of bound states for each component + * @param [in,out] workSpace Memory buffer for updated data + * @return Time derivatives of externally dependent parameters + */ + inline std::tuple updateTimeDerivative(double t, unsigned int secIdx, const ColumnPosition& colPos, unsigned int nComp, unsigned int const* nChannelStates, LinearBufferAllocator& workSpace) const + { + return std::make_tuple(&_localParams, nullptr); + } + +protected: + + /** + * @brief Validates recently read parameters + * @param [in] nComp Number of components + * @param [in] nChannelStates Array with number of bound states for each component + * @return @c true if the parameters were validated successfully, otherwise @c false + */ + inline bool validateConfig(unsigned int nComp, unsigned int const* nChannelStates) + { + if ((_kA.size() != _kD.size()) || (_kA.size() < nComp)) + throw InvalidParameterException("LIN_KA and LIN_KD have to have the same size"); + + return true; + } + + ConstParams _localParams; //!< Actual parameter data + + // Handlers provide configure(), reserve(), and registerParam() for parameters + ScalarComponentDependentParameter _kA; //!< Handler for adsorption rate + ScalarComponentDependentParameter _kD; //!< Handler for desorption rate +}; + +/** + * @brief Handles linear exchange model parameters that depend on an external function + */ +class ExtLinearParamHandler : public ExternalParamHandlerBase +{ +public: + + /** + * @brief Holds actual parameter data + * @details The parameter data will be stored in a memory buffer. This requires + * control over the location of the data, which is not provided by + * std::vector and util::SlicedVector. There are "local" equivalents, + * util::LocalVector and util::LocalSlicedVector, for these types that + * allow to control the placement of the payload data. A single value (i.e., + * a single active or double) can simply be put in the struct. + */ + struct VariableParams + { + util::LocalVector kA; //!< Adsorption rate + util::LocalVector kD; //!< Desorption rate + }; + + typedef VariableParams params_t; + typedef ConstBufferedScalar ParamsHandle; + + ExtLinearParamHandler() CADET_NOEXCEPT { } + + /** + * @brief Returns name of the exchange model + * @return Name of the exchange model + */ + static const char* identifier() { return "EXT_LINEAR"; } + + /** + * @brief Reads parameters and verifies them + * @details See IBindingModel::configure() for details. + * @param [in] paramProvider IParameterProvider used for reading parameters + * @param [in] nComp Number of components + * @param [in] nChannelStates Array with number of bound states for each component + * @return @c true if the parameters were read and validated successfully, otherwise @c false + */ + inline bool configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nChannelStates) + { + _kA.configure("LIN_KA", paramProvider, nComp, nChannelStates); + _kD.configure("LIN_KD", paramProvider, nComp, nChannelStates); + + // Number of externally dependent parameters (2) needs to be given to ExternalParamHandlerBase::configure() + ExternalParamHandlerBase::configure(paramProvider, 2); + return validateConfig(nComp, nChannelStates); + } + + /** + * @brief Registers all local parameters in a map for further use + * @param [in,out] parameters Map in which the parameters are stored + * @param [in] unitOpIdx Index of the unit operation used for registering the parameters + * @param [in] nComp Number of components + * @param [in] nChannelStates Array with number of bound states for each component + */ + inline void registerParameters(std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nChannelStates) + { + _kA.registerParam("LIN_KA", parameters, unitOpIdx, parTypeIdx, nComp, nChannelStates); + _kD.registerParam("LIN_KD", parameters, unitOpIdx, parTypeIdx, nComp, nChannelStates); + } + + /** + * @brief Reserves space in the storage of the parameters + * @param [in] numElem Number of elements (total) + * @param [in] numSlices Number of slices / exchange site types + * @param [in] nComp Number of components + * @param [in] nChannelStates Array with number of bound states for each component + */ + inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nChannelStates) \ + { + _kA.reserve(numElem, numSlices, nComp, nChannelStates); + _kD.reserve(numElem, numSlices, nComp, nChannelStates); + } + + /** + * @brief Updates the parameter cache in order to take the external profile into account + * @details The data and the returned value are constructed in the given @p workSpace memory buffer. + * @param [in] t Current time + * @param [in] z Axial coordinate in the column + * @param [in] r Radial coordinate in the bead + * @param [in] secIdx Index of the current section + * @param [in] nComp Number of components + * @param [in] nChannelStates Array with number of bound states for each component + * @param [in,out] workSpace Memory buffer for updated data + * @return Externally dependent parameter values + */ + inline ParamsHandle update(double t, unsigned int secIdx, const ColumnPosition& colPos, unsigned int nComp, unsigned int const* nChannelStates, LinearBufferAllocator& workSpace) const + { + // Allocate params_t and buffer for function evaluation + BufferedScalar localParams = workSpace.scalar(); + BufferedArray extFunBuffer = workSpace.array(2); + + // Evaluate external functions in buffer + evaluateExternalFunctions(t, secIdx, colPos, 2, static_cast(extFunBuffer)); + + // Prepare the buffer for the data and update the data + _kA.prepareCache(localParams->kA, workSpace); + _kA.update(cadet::util::dataOfLocalVersion(localParams->kA), extFunBuffer[0], nComp, nChannelStates); + + _kD.prepareCache(localParams->kD, workSpace); + _kD.update(cadet::util::dataOfLocalVersion(localParams->kD), extFunBuffer[1], nComp, nChannelStates); + + return localParams; + } + + /** + * @brief Calculates time derivative in case of external dependence + * @details The time derivatives are constructed in the given @p workSpace memory buffer. + * @param [in] t Current time + * @param [in] z Axial coordinate in the column + * @param [in] r Radial coordinate in the bead + * @param [in] secIdx Index of the current section + * @param [in] nComp Number of components + * @param [in] nChannelStates Array with number of bound states for each component + * @param [in,out] workSpace Memory buffer for updated data + * @return Tuple with externally dependent parameter values and their time derivatives + */ + inline std::tuple updateTimeDerivative(double t, unsigned int secIdx, const ColumnPosition& colPos, unsigned int nComp, unsigned int const* nChannelStates, LinearBufferAllocator& workSpace) const + { + // Allocate params_t for parameters and their time derivatives + BufferedScalar localParams = workSpace.scalar(); + BufferedScalar p = workSpace.scalar(); + + // Allocate buffer for external function values and their time derivatives + BufferedArray extFunBuffer = workSpace.array(2); + BufferedArray extDerivBuffer = workSpace.array(2); + + // Evaluate external functions and their time derivatives + evaluateExternalFunctions(t, secIdx, colPos, 2, static_cast(extFunBuffer)); + evaluateTimeDerivativeExternalFunctions(t, secIdx, colPos, 2, static_cast(extDerivBuffer)); + + // Prepare the buffer for the data and update the data + _kA.prepareCache(localParams->kA, workSpace); + _kA.update(cadet::util::dataOfLocalVersion(localParams->kA), extFunBuffer[0], nComp, nChannelStates); + + _kA.prepareCache(p->kA, workSpace); + _kA.updateTimeDerivative(cadet::util::dataOfLocalVersion(p->kA), extFunBuffer[0], extDerivBuffer[0], nComp, nChannelStates); + + _kD.prepareCache(localParams->kD, workSpace); + _kD.update(cadet::util::dataOfLocalVersion(localParams->kD), extFunBuffer[1], nComp, nChannelStates); + + _kD.prepareCache(p->kD, workSpace); + _kD.updateTimeDerivative(cadet::util::dataOfLocalVersion(p->kD), extFunBuffer[1], extDerivBuffer[1], nComp, nChannelStates); + + return std::make_tuple(std::move(localParams), std::move(p));; + } + + /** + * @brief Returns how much memory is required for caching in bytes + * @details Memory size in bytes. + * @param [in] nComp Number of components + * @param [in] totalNumBoundStates Total number of bound states + * @param [in] nChannelStates Array with bound states for each component + * @return Memory size in bytes + */ + inline std::size_t cacheSize(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nChannelStates) const CADET_NOEXCEPT + { + // Required buffer memory: + // + params_t object + // + buffer for external function evaluations (2 parameters) + // + buffer for external function time derivative evaluations (2 parameters) + // + buffer for actual parameter data (memory for _kA data + memory for _kD data) + // + buffer for parameter time derivatives (memory for _kA data + memory for _kD data) + return 2 * sizeof(params_t) + alignof(params_t) + + 2 * 2 * sizeof(double) + alignof(double) + + 2 * (_kA.additionalDynamicMemory(nComp, totalNumBoundStates, nChannelStates) + _kD.additionalDynamicMemory(nComp, totalNumBoundStates, nChannelStates)); + } + +protected: + + /** + * @brief Validates recently read parameters + * @param [in] nComp Number of components + * @param [in] nChannelStates Array with number of bound states for each component + * @return @c true if the parameters were validated successfully, otherwise @c false + */ + inline bool validateConfig(unsigned int nComp, unsigned int const* nChannelStates) + { + if ((_kA.size() != _kD.size()) || (_kA.size() < nComp)) + throw InvalidParameterException("LIN_KA and LIN_KD have to have the same size"); + + return true; + } + + // Handlers provide configure(), reserve(), and registerParam() for parameters + ExternalScalarComponentDependentParameter _kA; //!< Handler for adsorption rate + ExternalScalarComponentDependentParameter _kD; //!< Handler for desorption rate +}; + + +/** + * @brief Defines the linear exchange model + * @details Implements the linear adsorption model \f$ \frac{\mathrm{d}q}{\mathrm{d}t} = k_a c_p - k_d q \f$. + * Multiple bound states are not supported. Components without bound state (i.e., non-exchange components) + * are supported. + * + * See @cite Guiochon2006 + * @tparam ParamHandler_t Type that can add support for external function dependence + */ +class LinearExchangeBase : public IPhaseTransitionModel +{ +public: + + LinearExchangeBase() : _nComp(0), _nChannel(0), _reactionQuasistationarity(0, false) { } + virtual ~LinearExchangeBase() CADET_NOEXCEPT { } + + static const char* identifier() { return ParamHandler_t::identifier(); } + virtual const char* name() const CADET_NOEXCEPT { return ParamHandler_t::identifier(); } + virtual bool requiresConfiguration() const CADET_NOEXCEPT { return true; } + virtual bool usesParamProviderInDiscretizationConfig() const CADET_NOEXCEPT { return true; } + + virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int nChannel, unsigned int const* boundOffset) + { + _nComp = nComp; + _nChannel = nChannel; + + // if (hasMultipleBoundStates(nChannel, nComp)) + // throw InvalidParameterException("Linear exchange model does not support multiple bound states"); + // Why do we throw this error? (see next line) + + _reactionQuasistationarity.resize(numExchangeStates(nChannel, nComp), false); + + return true; + } + + virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) + { + _parameters.clear(); + + // Read parameters (k_a and k_d) + _paramHandler.configure(paramProvider, _nComp, _nChannel); + + // Register parameters + _paramHandler.registerParameters(_parameters, unitOpIdx, parTypeIdx, _nComp, _nChannel); + + return true; + } + + virtual void fillBoundPhaseInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) const CADET_NOEXCEPT + { + unsigned int ctr = 0; + for (int c = 0; c < _nComp; ++c) + { + for (unsigned int bp = 0; bp < _nChannel[c]; ++bp, ++ctr) + params[ctr] = makeParamId(hashString("INIT_Q"), unitOpIdx, c, parTypeIdx, bp, ReactionIndep, SectionIndep); + } + } + + virtual std::unordered_map getAllParameterValues() const + { + std::unordered_map data; + std::transform(_parameters.begin(), _parameters.end(), std::inserter(data, data.end()), + [](const std::pair& p) { return std::make_pair(p.first, static_cast(*p.second)); }); + + return data; + } + + virtual bool hasParameter(const ParameterId& pId) const + { + return _parameters.find(pId) != _parameters.end(); + } + + virtual bool setParameter(const ParameterId& pId, int value) + { + return false; + } + + virtual bool setParameter(const ParameterId& pId, double value) + { + auto paramHandle = _parameters.find(pId); + if (paramHandle != _parameters.end()) + { + paramHandle->second->setValue(value); + return true; + } + + return false; + } + + virtual bool setParameter(const ParameterId& pId, bool value) + { + return false; + } + + virtual active* getParameter(const ParameterId& pId){ + auto paramHandle = _parameters.find(pId); + if (paramHandle != _parameters.end()) + { + return paramHandle->second; + } + + return nullptr; + } + + virtual unsigned int workspaceSize(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nChannelStates) const CADET_NOEXCEPT + { + return _paramHandler.cacheSize(nComp, totalNumBoundStates, nChannelStates); + } + + virtual void setExternalFunctions(IExternalFunction** extFuns, unsigned int size) { _paramHandler.setExternalFunctions(extFuns, size); } + + // The next three flux() function implementations and two analyticJacobian() function + // implementations are usually hidden behind + // CADET_BINDINGMODELBASE_BOILERPLATE + // which just expands to the six implementations below. + + virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, + active const* y, active const* yCp, active* res, LinearBufferAllocator workSpace, WithParamSensitivity) const + { + return fluxImpl(t, secIdx, colPos, y, yCp, res, workSpace); + } + + virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, + active const* y, active const* yCp, active* res, LinearBufferAllocator workSpace, WithoutParamSensitivity) const + { + return fluxImpl(t, secIdx, colPos, y, yCp, res, workSpace); + } + + virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, + double const* y, double const* yCp, active* res, LinearBufferAllocator workSpace) const + { + return fluxImpl(t, secIdx, colPos, y, yCp, res, workSpace); + } + + virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, + double const* y, double const* yCp, double* res, LinearBufferAllocator workSpace) const + { + return fluxImpl(t, secIdx, colPos, y, yCp, res, workSpace); + } + + virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const + { + jacobianImpl(t, secIdx, colPos, y, offsetCp, jac, workSpace); + } + + virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const + { + jacobianImpl(t, secIdx, colPos, y, offsetCp, jac, workSpace); + } + + virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const + { + jacobianImpl(t, secIdx, colPos, y, offsetCp, jac, workSpace); + } + + virtual void timeDerivativeQuasiStationaryFluxes(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yCp, double const* y, double* dResDt, LinearBufferAllocator workSpace) const + { + if (!hasQuasiStationaryReactions()) + return; + + if (!ParamHandler_t::dependsOnTime()) + return; + + // Update external function and compute time derivative of parameters + typename ParamHandler_t::ParamsHandle p; + typename ParamHandler_t::ParamsHandle dpDt; + std::tie(p, dpDt) = _paramHandler.updateTimeDerivative(t, secIdx, colPos, _nComp, _nChannel, workSpace); + + unsigned int bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nChannel[i] == 0) + continue; + + dResDt[bndIdx] = -static_cast(dpDt->kA[i]) * yCp[i] + static_cast(dpDt->kD[i]) * y[bndIdx]; + + // Next bound component + ++bndIdx; + } + } + + virtual bool preConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const + { + // Due to mass conservation, we need to run a nonlinear solver (although the model is simple). + return true; + } + + virtual void postConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const + { + // There's nothing to do here since the algebraic equations have already been solved in preConsistentInitialState() + } + + virtual bool hasSalt() const CADET_NOEXCEPT { return false; } + virtual bool supportsMultistate() const CADET_NOEXCEPT { return false; } + virtual bool supportsNonExchange() const CADET_NOEXCEPT { return true; } + + virtual bool hasQuasiStationaryReactions() const CADET_NOEXCEPT + { + return std::any_of(_reactionQuasistationarity.begin(), _reactionQuasistationarity.end(), [](int i) -> bool { return i; }); + } + + virtual bool hasDynamicReactions() const CADET_NOEXCEPT + { + return std::any_of(_reactionQuasistationarity.begin(), _reactionQuasistationarity.end(), [](int i) -> bool { return !static_cast(i); }); + } + + virtual bool dependsOnTime() const CADET_NOEXCEPT { return ParamHandler_t::dependsOnTime(); } + virtual bool requiresWorkspace() const CADET_NOEXCEPT { return ParamHandler_t::requiresWorkspace(); } + virtual unsigned int requiredADdirs() const CADET_NOEXCEPT { return 0; } + virtual int const* reactionQuasiStationarity() const CADET_NOEXCEPT { return _reactionQuasistationarity.data(); } + +protected: + int _nComp; //!< Number of components + unsigned int const* _nChannel; //!< Array with number of bound states for each component + std::vector _reactionQuasistationarity; //!< Determines whether each bound state is quasi-stationary (@c true) or not (@c false) + + ParamHandler_t _paramHandler; //!< Parameters + + std::unordered_map _parameters; //!< Map used to translate ParameterIds to actual variables + + virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; } + + template + + int fluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, + StateType const* y, StateType const* yCp, ResidualType* res, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nChannel, workSpace); + + // Implement -k_a * c_{p,i} + k_d * q_i + + unsigned int bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nChannel[i] == 0) + continue; + + res[bndIdx] = -static_cast(p->kA[i]) * yCp[i] + static_cast(p->kD[i]) * y[bndIdx]; + + // Next bound component + ++bndIdx; + } + + return 0; + } + + template + inline void jacobianImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, RowIterator jac, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nChannel, workSpace); + + int bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nChannel[i] == 0) + continue; + + jac[0] = static_cast(p->kD[i]); // dres / dq_i + jac[i - bndIdx - offsetCp] = -static_cast(p->kA[i]); // dres / dc_{p,i} + // The distance from liquid phase to solid phase is reduced for each non-exchange component + // since a bound state is neglected. The number of neglected bound states so far is i - bndIdx. + // Thus, by going back offsetCp - (i - bndIdx) = -[ i - bndIdx - offsetCp ] we get to the corresponding + // liquid phase component. + + ++bndIdx; + ++jac; + } + } + +}; + +typedef LinearExchangeBase LinearExchange; +typedef LinearExchangeBase ExternalLinearExchange; + +namespace exchange +{ + void registerLinearModel(std::unordered_map>& exchange) + { + exchange[LinearExchange::identifier()] = []() { return new LinearExchange(); }; + exchange[ExternalLinearExchange::identifier()] = []() { return new ExternalLinearExchange(); }; + } +} // namespace exchange + +} // namespace model + +} // namespace cadet From 965e8c74d3ea5e85cfbad637882e4c216d37a6f9 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Wed, 31 Jul 2024 14:35:30 +0200 Subject: [PATCH 02/17] Add comments for possible changes --- .../model/MultiChannelTransportModel.cpp | 37 ++++++++++++++++++- src/libcadet/model/binding/LinearBinding.cpp | 6 +-- .../model/exchange/LinearExchange.cpp | 22 ++++------- 3 files changed, 45 insertions(+), 20 deletions(-) diff --git a/src/libcadet/model/MultiChannelTransportModel.cpp b/src/libcadet/model/MultiChannelTransportModel.cpp index 1d211862d..c82626f3e 100644 --- a/src/libcadet/model/MultiChannelTransportModel.cpp +++ b/src/libcadet/model/MultiChannelTransportModel.cpp @@ -18,6 +18,7 @@ #include "cadet/SolutionRecorder.hpp" #include "ConfigurationHelper.hpp" #include "model/ReactionModel.hpp" +#include "model/PhaseTransitionModel.hpp" #include "SimulationTypes.hpp" #include "Stencil.hpp" @@ -407,13 +408,24 @@ bool MultiChannelTransportModel::configureModelDiscretization(IParameterProvider if (_dynReactionBulk->usesParamProviderInDiscretizationConfig()) paramProvider.popScope(); } + bool exchangeConfSuccess = true; + + //_exchange = nullptr + std::vector exchModelNames = { "NONE" }; + if (paramProvider.exists("EXCHANGE_MODEL")) + exchModelNames = paramProvider.getStringArray("EXCHANGE_MODEL"); + + + bool exchConfSuccess = true; + //exchConfSuccess = _exchange->configureModelDiscretization(paramProvider, _disc.nComp, _disc.nChannel , _disc.boundOffset) && bindingConfSuccess; + //Q: + i * _disc.nBound + i * _disc.nComp for LRMP const bool transportSuccess = _convDispOp.configureModelDiscretization(paramProvider, helper, _disc.nComp, _disc.nCol, _disc.nChannel, _dynReactionBulk); // Setup the memory for tempState based on state vector _tempState = new double[numDofs()]; - return transportSuccess && reactionConfSuccess; + return transportSuccess && reactionConfSuccess && exchConfSuccess; } bool MultiChannelTransportModel::configure(IParameterProvider& paramProvider) @@ -437,8 +449,29 @@ bool MultiChannelTransportModel::configure(IParameterProvider& paramProvider) dynReactionConfSuccess = _dynReactionBulk->configure(paramProvider, _unitOpIdx, ParTypeIndep); paramProvider.popScope(); } + /* + if (!_exchange.empty()){ + + const unsigned int maxChannelStates = *std::max_element(_disc.strideBound, _disc.strideBound + _disc.nParType); + std::vector initParams(maxBoundStates); //Q was macht initParams ? + + _exchange->fillChannelInitialParameters(initParams.data(), _unitOpIdx, type); - return transportSuccess && dynReactionConfSuccess; + active* const iq = _initQ.data() + _disc.nBoundBeforeType[type]; + for (unsigned int i = 0; i < _disc.strideBound[type]; ++i){ + _parameters[initParams[i]] = iq + i; + } + + } + */ + bool exchangeConfSuccess = true; + /* if (!_exchange.empty()){ + + exchangeConfSuccess = _exchange-> configure(paramProvider, _unitOpIdx, type) && bindingConfSuccess; + + } + */ + return transportSuccess && dynReactionConfSuccess && exchangeConfSuccess; } diff --git a/src/libcadet/model/binding/LinearBinding.cpp b/src/libcadet/model/binding/LinearBinding.cpp index 676a02e9d..0247d5cd4 100644 --- a/src/libcadet/model/binding/LinearBinding.cpp +++ b/src/libcadet/model/binding/LinearBinding.cpp @@ -90,7 +90,6 @@ class LinearParamHandler : public ConstParamHandlerBase static const char* identifier() { return "LINEAR"; } LinearParamHandler() CADET_NOEXCEPT : _kA(&_localParams.kA), _kD(&_localParams.kD) { } - /** * @brief Reads parameters and verifies them * @details See IBindingModel::configure() for details. @@ -126,8 +125,7 @@ class LinearParamHandler : public ConstParamHandlerBase * @param [in] nComp Number of components * @param [in] nBoundStates Array with number of bound states for each component */ - inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nBoundStates) \ - { + inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nBoundStates) { _kA.reserve(numElem, numSlices, nComp, nBoundStates); _kD.reserve(numElem, numSlices, nComp, nBoundStates); } @@ -389,7 +387,7 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase * See @cite Guiochon2006 * @tparam ParamHandler_t Type that can add support for external function dependence */ -template +template // nicht fuer LinearExchange class LinearBindingBase : public IBindingModel { public: diff --git a/src/libcadet/model/exchange/LinearExchange.cpp b/src/libcadet/model/exchange/LinearExchange.cpp index f05496843..c62559afa 100644 --- a/src/libcadet/model/exchange/LinearExchange.cpp +++ b/src/libcadet/model/exchange/LinearExchange.cpp @@ -364,7 +364,7 @@ class LinearExchangeBase : public IPhaseTransitionModel { public: - LinearExchangeBase() : _nComp(0), _nChannel(0), _reactionQuasistationarity(0, false) { } + LinearExchangeBase() : _nComp(0), _nChannel(0) { } virtual ~LinearExchangeBase() CADET_NOEXCEPT { } static const char* identifier() { return ParamHandler_t::identifier(); } @@ -375,13 +375,7 @@ class LinearExchangeBase : public IPhaseTransitionModel virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int nChannel, unsigned int const* boundOffset) { _nComp = nComp; - _nChannel = nChannel; - - // if (hasMultipleBoundStates(nChannel, nComp)) - // throw InvalidParameterException("Linear exchange model does not support multiple bound states"); - // Why do we throw this error? (see next line) - - _reactionQuasistationarity.resize(numExchangeStates(nChannel, nComp), false); + _nChannel = nChannel; // nChannel realy int ? -> by not nBoundStates not ? return true; } @@ -390,8 +384,8 @@ class LinearExchangeBase : public IPhaseTransitionModel { _parameters.clear(); - // Read parameters (k_a and k_d) - _paramHandler.configure(paramProvider, _nComp, _nChannel); + // Read parameters (k_a and k_d,corss section area) + _paramHandler.configure(paramProvider, _nComp, _nChannel); // add cross section area in parameter // Register parameters _paramHandler.registerParameters(_parameters, unitOpIdx, parTypeIdx, _nComp, _nChannel); @@ -399,16 +393,16 @@ class LinearExchangeBase : public IPhaseTransitionModel return true; } - virtual void fillBoundPhaseInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) const CADET_NOEXCEPT + virtual void fillExchangeInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) const CADET_NOEXCEPT { unsigned int ctr = 0; for (int c = 0; c < _nComp; ++c) { - for (unsigned int bp = 0; bp < _nChannel[c]; ++bp, ++ctr) - params[ctr] = makeParamId(hashString("INIT_Q"), unitOpIdx, c, parTypeIdx, bp, ReactionIndep, SectionIndep); + for (unsigned int bp = 0; bp < _nChannel; ++bp, ++ctr) + params[ctr] = makeParamId(hashString("INIT_C"), unitOpIdx, c, parTypeIdx, bp, ReactionIndep, SectionIndep); } } - + // ----------------------------------------------------------------------------------------------------------------------------// virtual std::unordered_map getAllParameterValues() const { std::unordered_map data; From 5c77e7ba91a50f17aa129dcdf1ad5e7583b7033e Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Thu, 1 Aug 2024 12:15:02 +0200 Subject: [PATCH 03/17] Add comments where to add exchange modul --- src/libcadet/model/MultiChannelTransportModel.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/libcadet/model/MultiChannelTransportModel.cpp b/src/libcadet/model/MultiChannelTransportModel.cpp index c82626f3e..96f36e35c 100644 --- a/src/libcadet/model/MultiChannelTransportModel.cpp +++ b/src/libcadet/model/MultiChannelTransportModel.cpp @@ -479,6 +479,8 @@ unsigned int MultiChannelTransportModel::threadLocalMemorySize() const CADET_NOE { LinearMemorySizer lms; + //Add exchange memory + // Memory for residualImpl() if (_dynReactionBulk && _dynReactionBulk->requiresWorkspace()) lms.fitBlock(_dynReactionBulk->workspaceSize(_disc.nComp, 0, nullptr)); @@ -611,7 +613,7 @@ int MultiChannelTransportModel::residual(const SimulationTime& simTime, const Co // Evaluate residual do not compute Jacobian or parameter sensitivities return residualImpl(simTime.t, simTime.secIdx, simState.vecStateY, simState.vecStateYdot, res, threadLocalMem); } - +//Q Do we need different types of residual calculation ? int MultiChannelTransportModel::jacobian(const SimulationTime& simTime, const ConstSimulationState& simState, double* const res, const AdJacobianParams& adJac, util::ThreadLocalStorage& threadLocalMem) { BENCH_SCOPE(_timerResidual); @@ -871,6 +873,7 @@ void MultiChannelTransportModel::multiplyWithJacobian(const SimulationTime& simT */ void MultiChannelTransportModel::multiplyWithDerivativeJacobian(const SimulationTime& simTime, const ConstSimulationState& simState, double const* sDot, double* ret) { + // Add Exchange here _convDispOp.multiplyWithDerivativeJacobian(simTime, sDot, ret); // Handle inlet DOFs (all algebraic) @@ -879,6 +882,7 @@ void MultiChannelTransportModel::multiplyWithDerivativeJacobian(const Simulation void MultiChannelTransportModel::setExternalFunctions(IExternalFunction** extFuns, unsigned int size) { + // Add exchange here } unsigned int MultiChannelTransportModel::localOutletComponentIndex(unsigned int port) const CADET_NOEXCEPT From 720a8c191b8d3a994e244b2f6bbc5be3c03e5361 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Mon, 5 Aug 2024 12:31:08 +0200 Subject: [PATCH 04/17] Create new exchange Modul Resolves #256 --- src/libcadet/CMakeLists.txt | 5 + .../model/MultiChannelTransportModel.cpp | 28 +++-- src/libcadet/model/PhaseTransitionModel.hpp | 4 +- src/libcadet/model/UnitOperationBase.hpp | 2 + .../model/exchange/LinearExchange.cpp | 102 +++++++++--------- 5 files changed, 73 insertions(+), 68 deletions(-) diff --git a/src/libcadet/CMakeLists.txt b/src/libcadet/CMakeLists.txt index 84aa4b9ae..d1f477105 100644 --- a/src/libcadet/CMakeLists.txt +++ b/src/libcadet/CMakeLists.txt @@ -101,6 +101,11 @@ set(LIBCADET_BINDINGMODEL_SOURCES ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/HICWaterOnHydrophobicSurfacesBinding.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/HICConstantWaterActivityBinding.cpp ) +set(LIBCADET_PHASETRANSITIONMODELINTERFACE_HPP_ + ${CMAKE_SOURCE_DIR}/src/libcadet/model/exchange/LinearExchange.cpp +) + + # LIBCADET_REACTIONMODEL_SOURCES holds all source files of reaction models set(LIBCADET_REACTIONMODEL_SOURCES diff --git a/src/libcadet/model/MultiChannelTransportModel.cpp b/src/libcadet/model/MultiChannelTransportModel.cpp index 96f36e35c..5f7f9a9ce 100644 --- a/src/libcadet/model/MultiChannelTransportModel.cpp +++ b/src/libcadet/model/MultiChannelTransportModel.cpp @@ -449,28 +449,24 @@ bool MultiChannelTransportModel::configure(IParameterProvider& paramProvider) dynReactionConfSuccess = _dynReactionBulk->configure(paramProvider, _unitOpIdx, ParTypeIndep); paramProvider.popScope(); } - /* - if (!_exchange.empty()){ - - const unsigned int maxChannelStates = *std::max_element(_disc.strideBound, _disc.strideBound + _disc.nParType); - std::vector initParams(maxBoundStates); //Q was macht initParams ? + - _exchange->fillChannelInitialParameters(initParams.data(), _unitOpIdx, type); + if (_exchange[0]){ - active* const iq = _initQ.data() + _disc.nBoundBeforeType[type]; - for (unsigned int i = 0; i < _disc.strideBound[type]; ++i){ - _parameters[initParams[i]] = iq + i; - } - + std::vector initParams(_disc.nChannel); //Q was macht initParams ? brauche ich maxBoundStates ? + _exchange[0]->fillChannelInitialParameters(initParams.data(), _unitOpIdx); } - */ + + + // Reconfigure exchange model bool exchangeConfSuccess = true; - /* if (!_exchange.empty()){ - - exchangeConfSuccess = _exchange-> configure(paramProvider, _unitOpIdx, type) && bindingConfSuccess; + if (_exchange[0] && paramProvider.exists("exchange") && _exchange[0]->requiresConfiguration()){ + paramProvider.pushScope("exchange"); + exchangeConfSuccess = _exchange[0]->configure(paramProvider, _unitOpIdx, cadet::ParTypeIndep); // Brauche ich PartypeIndep ? + paramProvider.popScope(); } - */ + return transportSuccess && dynReactionConfSuccess && exchangeConfSuccess; } diff --git a/src/libcadet/model/PhaseTransitionModel.hpp b/src/libcadet/model/PhaseTransitionModel.hpp index e2981d525..cc667684e 100644 --- a/src/libcadet/model/PhaseTransitionModel.hpp +++ b/src/libcadet/model/PhaseTransitionModel.hpp @@ -114,6 +114,8 @@ class IPhaseTransitionModel */ virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) = 0; + virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx) = 0; + /** * @brief Returns the ParameterId of all bound phase initial conditions (equations) * @details The array has to be filled in the order of the equations. @@ -130,7 +132,7 @@ class IPhaseTransitionModel * @param [in] unitOpIdx Index of the unit operation this binding model belongs to * @param [in] parTypeIdx Index of the particle type this binding model belongs to */ - virtual void fillChannelInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) const CADET_NOEXCEPT = 0; + virtual void fillChannelInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx) const CADET_NOEXCEPT = 0; /** * @brief Sets external functions for this binding model diff --git a/src/libcadet/model/UnitOperationBase.hpp b/src/libcadet/model/UnitOperationBase.hpp index ad4b569fd..8c4cc2515 100644 --- a/src/libcadet/model/UnitOperationBase.hpp +++ b/src/libcadet/model/UnitOperationBase.hpp @@ -34,6 +34,7 @@ namespace model { class IBindingModel; +class IPhaseTransitionModel; class IDynamicReactionModel; /** @@ -78,6 +79,7 @@ class UnitOperationBase : public IUnitOperation UnitOpIdx _unitOpIdx; //!< Unit operation index std::vector _binding; //!< Binding model + std::vector _exchange; //!< Exchange model bool _singleBinding; //!< Determines whether only a single binding model is present std::vector _dynReaction; //!< Dynamic reaction model bool _singleDynReaction; //!< Determines whether only a single dynamic reaction model is present diff --git a/src/libcadet/model/exchange/LinearExchange.cpp b/src/libcadet/model/exchange/LinearExchange.cpp index c62559afa..95d4f7cd3 100644 --- a/src/libcadet/model/exchange/LinearExchange.cpp +++ b/src/libcadet/model/exchange/LinearExchange.cpp @@ -67,14 +67,14 @@ class LinearParamHandler : public ConstParamHandlerBase * @details See IexchangeModel::configure() for details. * @param [in] paramProvider IParameterProvider used for reading parameters * @param [in] nComp Number of components - * @param [in] nChannelStates Array with number of bound states for each component + * @param [in] nChannel Array with number of bound states for each component * @return @c true if the parameters were read and validated successfully, otherwise @c false */ - inline bool configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nChannelStates) + inline bool configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nChannel) { - _kA.configure("LIN_KA", paramProvider, nComp, nChannelStates); - _kD.configure("LIN_KD", paramProvider, nComp, nChannelStates); - return validateConfig(nComp, nChannelStates); + _kA.configure("LIN_KA", paramProvider, nComp, nChannel); + _kD.configure("LIN_KD", paramProvider, nComp, nChannel); + return validateConfig(nComp, nChannel); } /** @@ -82,12 +82,12 @@ class LinearParamHandler : public ConstParamHandlerBase * @param [in,out] parameters Map in which the parameters are stored * @param [in] unitOpIdx Index of the unit operation used for registering the parameters * @param [in] nComp Number of components - * @param [in] nChannelStates Array with number of bound states for each component + * @param [in] nChannel Array with number of bound states for each component */ - inline void registerParameters(std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nChannelStates) + inline void registerParameters(std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nChannel) { - _kA.registerParam("LIN_KA", parameters, unitOpIdx, parTypeIdx, nComp, nChannelStates); - _kD.registerParam("LIN_KD", parameters, unitOpIdx, parTypeIdx, nComp, nChannelStates); + _kA.registerParam("LIN_KA", parameters, unitOpIdx, parTypeIdx, nComp, nChannel); + _kD.registerParam("LIN_KD", parameters, unitOpIdx, parTypeIdx, nComp, nChannel); } /** @@ -95,12 +95,12 @@ class LinearParamHandler : public ConstParamHandlerBase * @param [in] numElem Number of elements (total) * @param [in] numSlices Number of slices / exchange site types * @param [in] nComp Number of components - * @param [in] nChannelStates Array with number of bound states for each component + * @param [in] nChannel Array with number of bound states for each component */ - inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nChannelStates) \ + inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nChannel) \ { - _kA.reserve(numElem, numSlices, nComp, nChannelStates); - _kD.reserve(numElem, numSlices, nComp, nChannelStates); + _kA.reserve(numElem, numSlices, nComp, nChannel); + _kD.reserve(numElem, numSlices, nComp, nChannel); } /** @@ -110,11 +110,11 @@ class LinearParamHandler : public ConstParamHandlerBase * @param [in] r Radial coordinate in the bead * @param [in] secIdx Index of the current section * @param [in] nComp Number of components - * @param [in] nChannelStates Array with number of bound states for each component + * @param [in] nChannel Array with number of bound states for each component * @param [in,out] workSpace Memory buffer for updated data * @return Externally dependent parameter values */ - inline ParamsHandle update(double t, unsigned int secIdx, const ColumnPosition& colPos, unsigned int nComp, unsigned int const* nChannelStates, LinearBufferAllocator& workSpace) const + inline ParamsHandle update(double t, unsigned int secIdx, const ColumnPosition& colPos, unsigned int nComp, unsigned int const* nChannel, LinearBufferAllocator& workSpace) const { return &_localParams; } @@ -126,11 +126,11 @@ class LinearParamHandler : public ConstParamHandlerBase * @param [in] r Radial coordinate in the bead * @param [in] secIdx Index of the current section * @param [in] nComp Number of components - * @param [in] nChannelStates Array with number of bound states for each component + * @param [in] nChannel Array with number of bound states for each component * @param [in,out] workSpace Memory buffer for updated data * @return Time derivatives of externally dependent parameters */ - inline std::tuple updateTimeDerivative(double t, unsigned int secIdx, const ColumnPosition& colPos, unsigned int nComp, unsigned int const* nChannelStates, LinearBufferAllocator& workSpace) const + inline std::tuple updateTimeDerivative(double t, unsigned int secIdx, const ColumnPosition& colPos, unsigned int nComp, unsigned int const* nChannel, LinearBufferAllocator& workSpace) const { return std::make_tuple(&_localParams, nullptr); } @@ -140,10 +140,10 @@ class LinearParamHandler : public ConstParamHandlerBase /** * @brief Validates recently read parameters * @param [in] nComp Number of components - * @param [in] nChannelStates Array with number of bound states for each component + * @param [in] nChannel Array with number of bound states for each component * @return @c true if the parameters were validated successfully, otherwise @c false */ - inline bool validateConfig(unsigned int nComp, unsigned int const* nChannelStates) + inline bool validateConfig(unsigned int nComp, unsigned int const* nChannel) { if ((_kA.size() != _kD.size()) || (_kA.size() < nComp)) throw InvalidParameterException("LIN_KA and LIN_KD have to have the same size"); @@ -196,17 +196,17 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase * @details See IBindingModel::configure() for details. * @param [in] paramProvider IParameterProvider used for reading parameters * @param [in] nComp Number of components - * @param [in] nChannelStates Array with number of bound states for each component + * @param [in] nChannel Array with number of bound states for each component * @return @c true if the parameters were read and validated successfully, otherwise @c false */ - inline bool configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nChannelStates) + inline bool configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nChannel) { - _kA.configure("LIN_KA", paramProvider, nComp, nChannelStates); - _kD.configure("LIN_KD", paramProvider, nComp, nChannelStates); + _kA.configure("LIN_KA", paramProvider, nComp, nChannel); + _kD.configure("LIN_KD", paramProvider, nComp, nChannel); // Number of externally dependent parameters (2) needs to be given to ExternalParamHandlerBase::configure() ExternalParamHandlerBase::configure(paramProvider, 2); - return validateConfig(nComp, nChannelStates); + return validateConfig(nComp, nChannel); } /** @@ -214,12 +214,12 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase * @param [in,out] parameters Map in which the parameters are stored * @param [in] unitOpIdx Index of the unit operation used for registering the parameters * @param [in] nComp Number of components - * @param [in] nChannelStates Array with number of bound states for each component + * @param [in] nChannel Array with number of bound states for each component */ - inline void registerParameters(std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nChannelStates) + inline void registerParameters(std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nChannel) { - _kA.registerParam("LIN_KA", parameters, unitOpIdx, parTypeIdx, nComp, nChannelStates); - _kD.registerParam("LIN_KD", parameters, unitOpIdx, parTypeIdx, nComp, nChannelStates); + _kA.registerParam("LIN_KA", parameters, unitOpIdx, parTypeIdx, nComp, nChannel); + _kD.registerParam("LIN_KD", parameters, unitOpIdx, parTypeIdx, nComp, nChannel); } /** @@ -227,12 +227,12 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase * @param [in] numElem Number of elements (total) * @param [in] numSlices Number of slices / exchange site types * @param [in] nComp Number of components - * @param [in] nChannelStates Array with number of bound states for each component + * @param [in] nChannel Array with number of bound states for each component */ - inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nChannelStates) \ + inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nChannel) \ { - _kA.reserve(numElem, numSlices, nComp, nChannelStates); - _kD.reserve(numElem, numSlices, nComp, nChannelStates); + _kA.reserve(numElem, numSlices, nComp, nChannel); + _kD.reserve(numElem, numSlices, nComp, nChannel); } /** @@ -243,11 +243,11 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase * @param [in] r Radial coordinate in the bead * @param [in] secIdx Index of the current section * @param [in] nComp Number of components - * @param [in] nChannelStates Array with number of bound states for each component + * @param [in] nChannel Array with number of bound states for each component * @param [in,out] workSpace Memory buffer for updated data * @return Externally dependent parameter values */ - inline ParamsHandle update(double t, unsigned int secIdx, const ColumnPosition& colPos, unsigned int nComp, unsigned int const* nChannelStates, LinearBufferAllocator& workSpace) const + inline ParamsHandle update(double t, unsigned int secIdx, const ColumnPosition& colPos, unsigned int nComp, unsigned int const* nChannel, LinearBufferAllocator& workSpace) const { // Allocate params_t and buffer for function evaluation BufferedScalar localParams = workSpace.scalar(); @@ -258,10 +258,10 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase // Prepare the buffer for the data and update the data _kA.prepareCache(localParams->kA, workSpace); - _kA.update(cadet::util::dataOfLocalVersion(localParams->kA), extFunBuffer[0], nComp, nChannelStates); + _kA.update(cadet::util::dataOfLocalVersion(localParams->kA), extFunBuffer[0], nComp, nChannel); _kD.prepareCache(localParams->kD, workSpace); - _kD.update(cadet::util::dataOfLocalVersion(localParams->kD), extFunBuffer[1], nComp, nChannelStates); + _kD.update(cadet::util::dataOfLocalVersion(localParams->kD), extFunBuffer[1], nComp, nChannel); return localParams; } @@ -274,11 +274,11 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase * @param [in] r Radial coordinate in the bead * @param [in] secIdx Index of the current section * @param [in] nComp Number of components - * @param [in] nChannelStates Array with number of bound states for each component + * @param [in] nChannel Array with number of bound states for each component * @param [in,out] workSpace Memory buffer for updated data * @return Tuple with externally dependent parameter values and their time derivatives */ - inline std::tuple updateTimeDerivative(double t, unsigned int secIdx, const ColumnPosition& colPos, unsigned int nComp, unsigned int const* nChannelStates, LinearBufferAllocator& workSpace) const + inline std::tuple updateTimeDerivative(double t, unsigned int secIdx, const ColumnPosition& colPos, unsigned int nComp, unsigned int const* nChannel, LinearBufferAllocator& workSpace) const { // Allocate params_t for parameters and their time derivatives BufferedScalar localParams = workSpace.scalar(); @@ -294,16 +294,16 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase // Prepare the buffer for the data and update the data _kA.prepareCache(localParams->kA, workSpace); - _kA.update(cadet::util::dataOfLocalVersion(localParams->kA), extFunBuffer[0], nComp, nChannelStates); + _kA.update(cadet::util::dataOfLocalVersion(localParams->kA), extFunBuffer[0], nComp, nChannel); _kA.prepareCache(p->kA, workSpace); - _kA.updateTimeDerivative(cadet::util::dataOfLocalVersion(p->kA), extFunBuffer[0], extDerivBuffer[0], nComp, nChannelStates); + _kA.updateTimeDerivative(cadet::util::dataOfLocalVersion(p->kA), extFunBuffer[0], extDerivBuffer[0], nComp, nChannel); _kD.prepareCache(localParams->kD, workSpace); - _kD.update(cadet::util::dataOfLocalVersion(localParams->kD), extFunBuffer[1], nComp, nChannelStates); + _kD.update(cadet::util::dataOfLocalVersion(localParams->kD), extFunBuffer[1], nComp, nChannel); _kD.prepareCache(p->kD, workSpace); - _kD.updateTimeDerivative(cadet::util::dataOfLocalVersion(p->kD), extFunBuffer[1], extDerivBuffer[1], nComp, nChannelStates); + _kD.updateTimeDerivative(cadet::util::dataOfLocalVersion(p->kD), extFunBuffer[1], extDerivBuffer[1], nComp, nChannel); return std::make_tuple(std::move(localParams), std::move(p));; } @@ -313,10 +313,10 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase * @details Memory size in bytes. * @param [in] nComp Number of components * @param [in] totalNumBoundStates Total number of bound states - * @param [in] nChannelStates Array with bound states for each component + * @param [in] nChannel Array with bound states for each component * @return Memory size in bytes */ - inline std::size_t cacheSize(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nChannelStates) const CADET_NOEXCEPT + inline std::size_t cacheSize(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nChannel) const CADET_NOEXCEPT { // Required buffer memory: // + params_t object @@ -326,7 +326,7 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase // + buffer for parameter time derivatives (memory for _kA data + memory for _kD data) return 2 * sizeof(params_t) + alignof(params_t) + 2 * 2 * sizeof(double) + alignof(double) - + 2 * (_kA.additionalDynamicMemory(nComp, totalNumBoundStates, nChannelStates) + _kD.additionalDynamicMemory(nComp, totalNumBoundStates, nChannelStates)); + + 2 * (_kA.additionalDynamicMemory(nComp, totalNumBoundStates, nChannel) + _kD.additionalDynamicMemory(nComp, totalNumBoundStates, nChannel)); } protected: @@ -334,10 +334,10 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase /** * @brief Validates recently read parameters * @param [in] nComp Number of components - * @param [in] nChannelStates Array with number of bound states for each component + * @param [in] nChannel Array with number of bound states for each component * @return @c true if the parameters were validated successfully, otherwise @c false */ - inline bool validateConfig(unsigned int nComp, unsigned int const* nChannelStates) + inline bool validateConfig(unsigned int nComp, unsigned int const* nChannel) { if ((_kA.size() != _kD.size()) || (_kA.size() < nComp)) throw InvalidParameterException("LIN_KA and LIN_KD have to have the same size"); @@ -384,8 +384,8 @@ class LinearExchangeBase : public IPhaseTransitionModel { _parameters.clear(); - // Read parameters (k_a and k_d,corss section area) - _paramHandler.configure(paramProvider, _nComp, _nChannel); // add cross section area in parameter + // Read parameters (k_a and k_d) + _paramHandler.configure(paramProvider, _nComp, _nChannel); // Register parameters _paramHandler.registerParameters(_parameters, unitOpIdx, parTypeIdx, _nComp, _nChannel); @@ -449,9 +449,9 @@ class LinearExchangeBase : public IPhaseTransitionModel return nullptr; } - virtual unsigned int workspaceSize(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nChannelStates) const CADET_NOEXCEPT + virtual unsigned int workspaceSize(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nChannel) const CADET_NOEXCEPT { - return _paramHandler.cacheSize(nComp, totalNumBoundStates, nChannelStates); + return _paramHandler.cacheSize(nComp, totalNumBoundStates, nChannel); } virtual void setExternalFunctions(IExternalFunction** extFuns, unsigned int size) { _paramHandler.setExternalFunctions(extFuns, size); } From 6b475b7a8e51ed72214e93549a617c6c6cc4bf3b Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Wed, 7 Aug 2024 10:15:35 +0200 Subject: [PATCH 05/17] flux and jac in LinearExchange.cpp Without template configurations --- .../model/MultiChannelTransportModel.cpp | 7 +- src/libcadet/model/PhaseTransitionModel.hpp | 334 +++--------------- src/libcadet/model/binding/LinearBinding.cpp | 4 +- .../model/exchange/LinearExchange.cpp | 234 +++++++----- ...ltiChannelConvectionDispersionOperator.cpp | 67 +--- ...ltiChannelConvectionDispersionOperator.hpp | 4 + 6 files changed, 218 insertions(+), 432 deletions(-) diff --git a/src/libcadet/model/MultiChannelTransportModel.cpp b/src/libcadet/model/MultiChannelTransportModel.cpp index 5f7f9a9ce..11d28aa65 100644 --- a/src/libcadet/model/MultiChannelTransportModel.cpp +++ b/src/libcadet/model/MultiChannelTransportModel.cpp @@ -450,7 +450,7 @@ bool MultiChannelTransportModel::configure(IParameterProvider& paramProvider) paramProvider.popScope(); } - + /* Zunkufscode if (_exchange[0]){ std::vector initParams(_disc.nChannel); //Q was macht initParams ? brauche ich maxBoundStates ? @@ -466,8 +466,9 @@ bool MultiChannelTransportModel::configure(IParameterProvider& paramProvider) exchangeConfSuccess = _exchange[0]->configure(paramProvider, _unitOpIdx, cadet::ParTypeIndep); // Brauche ich PartypeIndep ? paramProvider.popScope(); } - - return transportSuccess && dynReactionConfSuccess && exchangeConfSuccess; + */ + + return transportSuccess && dynReactionConfSuccess; } diff --git a/src/libcadet/model/PhaseTransitionModel.hpp b/src/libcadet/model/PhaseTransitionModel.hpp index cc667684e..56563dd3e 100644 --- a/src/libcadet/model/PhaseTransitionModel.hpp +++ b/src/libcadet/model/PhaseTransitionModel.hpp @@ -51,305 +51,69 @@ class IPhaseTransitionModel virtual ~IPhaseTransitionModel() CADET_NOEXCEPT { } - /** - * @brief Returns the name of the binding model - * @details This name is also used to identify and create the binding model in the factory. - * @return Name of the binding model - */ virtual const char* name() const CADET_NOEXCEPT = 0; - /** - * @brief Returns whether the binding model requires additional parameters supplied by configure() - * @details After construction of an IPhaseTransitionModel object, configureModelDiscretization() is called. - * The binding model may require to read additional parameters from the adsorption group - * of the parameter provider. This opportunity is given by a call to configure(). - * However, a binding model may not require this. This function communicates to the outside - * whether a call to configure() is necessary. - * @return @c true if configure() has to be called, otherwise @c false - */ - virtual bool requiresConfiguration() const CADET_NOEXCEPT = 0; - - /** - * @brief Returns whether the dynamic reaction model uses the IParameterProvider in configureModelDiscretization() - * @details If the IParameterProvider is used in configureModelDiscretization(), it has to be in the correct scope. - * @return @c true if the IParameterProvider is used in configureModelDiscretization(), otherwise @c false - */ - virtual bool usesParamProviderInDiscretizationConfig() const CADET_NOEXCEPT = 0; - - /** - * @brief Sets fixed parameters of the binding model (e.g., the number of components and bound states) - * @details This function is called prior to configure() by the underlying model. - * It can only be called once. Whereas non-structural model parameters - * (e.g., rate constants) are configured by configure(), this function - * sets structural parameters (e.g., number of components and bound states). - * In particular, it determines whether each binding reaction is quasi-stationary. - * - * Before this function is called, usesParamProviderInDiscretizationConfig() - * is queried to determine whether the @p paramProvider is used by this - * function. If it is used, the @p paramProvider is put in the corresponding - * parameter scope. On exit, the @p paramProvider has to point to the same - * scope. - * @param [in] paramProvider Parameter provider - * @param [in] nComp Number of components - * @param [in] nBound Array of size @p nComp which contains the number of bound states for each component - * @param [in] boundOffset Array of size @p nComp with offsets to the first bound state of each component beginning from the solid phase - */ - virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBound, unsigned int const* boundOffset) = 0; - - /** - * @brief Configures the model by extracting all non-structural parameters (e.g., model parameters) from the given @p paramProvider - * @details The scope of the cadet::IParameterProvider is left unchanged on return. - * - * The structure of the model is left unchanged, that is, the number of degrees of - * freedom stays the same (e.g., number of bound states is left unchanged). Only - * true (non-structural) model parameters are read and changed. - * - * This function may only be called if configureModelDiscretization() has been called - * in the past. Contrary to configureModelDiscretization(), it can be called multiple - * times. - * @param [in] paramProvider Parameter provider - * @param [in] unitOpIdx Index of the unit operation this binding model belongs to - * @param [in] parTypeIdx Index of the particle type this binding model belongs to - * @return @c true if the configuration was successful, otherwise @c false - */ + //virtual bool requiresConfiguration() const CADET_NOEXCEPT = 0; + + //virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBound, unsigned int const* boundOffset) = 0; + virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) = 0; - virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx) = 0; - - /** - * @brief Returns the ParameterId of all bound phase initial conditions (equations) - * @details The array has to be filled in the order of the equations. - * @param [out] params Array with ParameterId objects to fill - * @param [in] unitOpIdx Index of the unit operation this binding model belongs to - * @param [in] parTypeIdx Index of the particle type this binding model belongs to - */ - virtual void fillBoundPhaseInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) const CADET_NOEXCEPT = 0; - - /** - * @brief Returns the ParameterId of all bound phase initial conditions (equations) - * @details The array has to be filled in the order of the equations. - * @param [out] params Array with ParameterId objects to fill - * @param [in] unitOpIdx Index of the unit operation this binding model belongs to - * @param [in] parTypeIdx Index of the particle type this binding model belongs to - */ - virtual void fillChannelInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx) const CADET_NOEXCEPT = 0; - - /** - * @brief Sets external functions for this binding model - * @details The external functions are not owned by this IPhaseTransitionModel. - * - * @param [in] extFuns Pointer to array of IExternalFunction objects of size @p size - * @param [in] size Number of elements in the IExternalFunction array @p extFuns - */ - virtual void setExternalFunctions(IExternalFunction** extFuns, unsigned int size) = 0; - - /** - * @brief Checks whether a given parameter exists - * @param [in] pId ParameterId that identifies the parameter uniquely - * @return @c true if the parameter exists, otherwise @c false - */ - virtual bool hasParameter(const ParameterId& pId) const = 0; - - /** - * @brief Returns all parameters with their current values that can be made sensitive - * @return Map with all parameters that can be made sensitive along with their current value - */ - virtual std::unordered_map getAllParameterValues() const = 0; - - /** - * @brief Sets a parameter value - * @details The parameter identified by its unique parameter is set to the given value. - * - * @param [in] pId ParameterId that identifies the parameter uniquely - * @param [in] value Value of the parameter - * - * @return @c true if the parameter has been successfully set to the given value, - * otherwise @c false (e.g., if the parameter is not available in this model) - */ - virtual bool setParameter(const ParameterId& pId, int value) = 0; - virtual bool setParameter(const ParameterId& pId, double value) = 0; - virtual bool setParameter(const ParameterId& pId, bool value) = 0; - - /** - * @brief Returns a pointer to the parameter identified by the given id - * @param [in] pId Parameter Id of the sensitive parameter - * @return a pointer to the parameter if the binding model contains the parameter, otherwise @c nullptr - */ - virtual active* getParameter(const ParameterId& pId) = 0; - - /** - * @brief Returns whether this binding model has a salt component - * @details The salt component is given by component 0, if present. - * @return @c true if salt is required, otherwise @c false - */ - virtual bool hasSalt() const CADET_NOEXCEPT = 0; + // virtual void fillBoundPhaseInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) const CADET_NOEXCEPT = 0; + + //virtual void fillChannelInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx) const CADET_NOEXCEPT = 0; + + + //virtual void setExternalFunctions(IExternalFunction** extFuns, unsigned int size) = 0; //X + + //virtual std::unordered_map getAllParameterValues() const = 0; + + //static const char* identifier() CADET_NOEXCEPT { return "IPhaseTransitionModel"; } + //virtual bool hasParameter(const ParameterId& pId) const = 0; + //virtual bool setParameter(const ParameterId& pId, int value) = 0; + //virtual bool setParameter(const ParameterId& pId, double value) = 0; + //virtual bool setParameter(const ParameterId& pId, bool value) = 0; + + //virtual active* getParameter(const ParameterId& pId) = 0; + + //virtual bool hasSalt() const CADET_NOEXCEPT = 0; + + //virtual bool supportsMultistate() const CADET_NOEXCEPT = 0;l bool supportsNonExchange() const CADET_NOEXCEPT = 0; + + + //virtual bool supportsNonBinding() const CADET_NOEXCEPT = 0; + + //virtual bool hasQuasiStationaryReactions() const CADET_NOEXCEPT = 0; + + //virtual bool hasDynamicReactions() const CADET_NOEXCEPT = 0; + + //virtual bool dependsOnTime() const CADET_NOEXCEPT = 0; + + //virtual bool requiresWorkspace() const CADET_NOEXCEPT = 0; + + //virtual unsigned int workspaceSize(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nBoundStates) const CADET_NOEXCEPT = 0; + + //virtual unsigned int requiredADdirs() const CADET_NOEXCEPT = 0; + + virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, active const* y, active* res) const = 0; + //virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, active const* y, double* res) const = 0; + //virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, double const* y, active* res) const = 0; + //virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, double const* y, double* res) const = 0; - /** - * @brief Returns whether this binding model supports multi-state binding - * @return @c true if multi-state binding is supported, otherwise @c false - */ - virtual bool supportsMultistate() const CADET_NOEXCEPT = 0; - - /** - * @brief Returns whether this binding model supports non-binding components - * @details Non-binding components do not have an entry in the solid phase. - * @return @c true if non-binding components are supported, otherwise @c false - */ - virtual bool supportsNonExchange() const CADET_NOEXCEPT = 0; - - /** - * @brief Returns whether this binding model supports non-binding components - * @details Non-binding components do not have an entry in the solid phase. - * @return @c true if non-binding components are supported, otherwise @c false - */ - virtual bool supportsNonBinding() const CADET_NOEXCEPT = 0; - - /** - * @brief Returns whether this binding model has quasi-stationary reaction fluxes - * @return @c true if quasi-stationary fluxes are present, otherwise @c false - */ - virtual bool hasQuasiStationaryReactions() const CADET_NOEXCEPT = 0; - - /** - * @brief Returns whether this binding model has dynamic reaction fluxes - * @return @c true if dynamic fluxes are present, otherwise @c false - */ - virtual bool hasDynamicReactions() const CADET_NOEXCEPT = 0; - - /** - * @brief Returns whether this binding model depends on time - * @details Binding models may depend on time if external functions are used. - * @return @c true if the model is time-dependent, otherwise @c false - */ - virtual bool dependsOnTime() const CADET_NOEXCEPT = 0; - - /** - * @brief Returns whether this binding model requires workspace - * @details The workspace may be required for consistent initialization and / or evaluation - * of residual and Jacobian. A workspace is a memory buffer whose size is given by - * workspaceSize(). - * @return @c true if the model requires a workspace, otherwise @c false - */ - virtual bool requiresWorkspace() const CADET_NOEXCEPT = 0; - - /** - * @brief Returns the size of the required workspace in bytes - * @details The memory is required for externally dependent binding models. - * @param [in] nComp Number of components - * @param [in] totalNumBoundStates Total number of bound states - * @param [in] nBoundStates Array with bound states for each component - * @return Size of the workspace in bytes - */ - virtual unsigned int workspaceSize(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nBoundStates) const CADET_NOEXCEPT = 0; - - /** - * @brief Returns the amount of required AD seed vectors / directions - * @details Only internally required AD directions count (e.g., for Jacobian computation). - * Directions used for parameter sensitivities should not be included here. - * @return The number of required AD seed vectors / directions - */ - virtual unsigned int requiredADdirs() const CADET_NOEXCEPT = 0; - - /** - * @brief Evaluates the fluxes - * @details The binding model is responsible for calculating the flux from the mobile to the solid phase. - * - * This function is called simultaneously from multiple threads. - * It needs to overwrite all values of @p res as the result array @p res is not - * zeroed on entry. - * @param [in] t Current time point - * @param [in] secIdx Index of the current section - * @param [in] colPos Position in normalized coordinates (column inlet = 0, column outlet = 1; outer shell = 1, inner center = 0) - * @param [in] y Pointer to first bound state of the first component in the current particle shell - * @param [in] yCp Pointer to first component of the mobile phase in the current particle shell - * @param [out] res Pointer to result array that is filled with the fluxes - * @param [in,out] workSpace Memory work space - * @return @c 0 on success, @c -1 on non-recoverable error, and @c +1 on recoverable error - */ - virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, active const* yCp, active* res, LinearBufferAllocator workSpace, WithParamSensitivity) const = 0; - virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, active const* yCp, active* res, LinearBufferAllocator workSpace, WithoutParamSensitivity) const = 0; - virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double const* yCp, active* res, LinearBufferAllocator workSpace) const = 0; - virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double const* yCp, double* res, LinearBufferAllocator workSpace) const = 0; - - /** - * @brief Evaluates the Jacobian of the fluxes analytically - * @details This function is called simultaneously from multiple threads. - * It can be left out (empty implementation) if AD is used to evaluate the Jacobian. - * - * The Jacobian matrix is assumed to be zeroed out by the caller. - * @param [in] t Current time point - * @param [in] secIdx Index of the current section - * @param [in] colPos Position in normalized coordinates (column inlet = 0, column outlet = 1; outer shell = 1, inner center = 0) - * @param [in] y Pointer to first bound state of the first component in the current particle shell - * @param [in] offsetCp Offset from the first component of the mobile phase in the current particle shell to @p y - * @param [in,out] jac Row iterator pointing to the first bound states row of the underlying matrix in which the Jacobian is stored - * @param [in,out] workSpace Memory work space - */ - virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const = 0; - virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const = 0; + virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, active const* y, active* res, linalg::CompressedSparseMatrix jac) const = 0; + //virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, active const* y, double* res, linalg::BandMatrix::RowIterator jac) const = 0; + //virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, active* res, linalg::BandMatrix::RowIterator jac) const = 0; + // virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, double* res, linalg::BandMatrix::RowIterator jac) const = 0; #ifdef ENABLE_DG - virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; + //virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; #endif - /** - * @brief Calculates the time derivative of the quasi-stationary bound state equations - * @details Calculates @f$ \frac{\partial \text{flux}_{\text{qs}}}{\partial t} @f$ for the quasi-stationary equations - * in the flux. - * - * This function is called simultaneously from multiple threads. - * It can be left out (empty implementation) which leads to slightly incorrect initial conditions - * when using externally dependent binding models. - * @param [in] t Current time point - * @param [in] secIdx Index of the current section - * @param [in] colPos Position in normalized coordinates (column inlet = 0, column outlet = 1; outer shell = 1, inner center = 0) - * @param [in] yCp Pointer to first component in the liquid phase of the current particle shell - * @param [in] y Pointer to first bound state of the first component in the current particle shell - * @param [out] dResDt Pointer to array that stores the time derivative - * @param [in,out] workSpace Memory work space - */ virtual void timeDerivativeQuasiStationaryFluxes(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yCp, double const* y, double* dResDt, LinearBufferAllocator workSpace) const = 0; - /** - * @brief Returns an array that determines whether each binding reaction is quasi-stationary - * @details Each entry represents a truth value (converts to @c bool) that determines whether - * the corresponding binding reaction is quasi-stationary (@c true) or not (@c false). - * @return Array that determines quasi-stationarity of binding reactions - */ virtual int const* reactionQuasiStationarity() const CADET_NOEXCEPT = 0; - /** - * @brief Returns whether a nonlinear solver is required to perform consistent initialization - * @details This function is called before a nonlinear solver attempts to solve the algebraic - * equations. If consistent initialization can be performed using analytic / direct - * formulas (i.e., without running a nonlinear solver), this function would be the - * correct place to do it. - * - * This function is called simultaneously from multiple threads. - * @param [in] t Current time point - * @param [in] secIdx Index of the current section - * @param [in] colPos Position in normalized coordinates (column inlet = 0, column outlet = 1; outer shell = 1, inner center = 0) - * @param [in,out] y Pointer to first bound state of the first component in the current particle shell - * @param [in] yCp Pointer to first component of the mobile phase in the current particle shell - * @param [in,out] workSpace Memory work space - * @return @c true if a nonlinear solver is required for consistent initialization, otherwise @c false - */ virtual bool preConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const = 0; - /** - * @brief Called after a nonlinear solver has attempted consistent initialization - * @details In consistent initialization, first preConsistentInitialState() is called. - * Depending on the return value, a nonlinear solver is applied to the algebraic - * equations and, afterwards, this function is called to clean up or refine the - * solution. - * - * This function is called simultaneously from multiple threads. - * @param [in] t Current time point - * @param [in] secIdx Index of the current section - * @param [in] colPos Position in normalized coordinates (column inlet = 0, column outlet = 1; outer shell = 1, inner center = 0) - * @param [in,out] y Pointer to first bound state of the first component in the current particle shell - * @param [in] yCp Pointer to first component of the mobile phase in the current particle shell - * @param [in,out] workSpace Memory work space - */ virtual void postConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const = 0; protected: diff --git a/src/libcadet/model/binding/LinearBinding.cpp b/src/libcadet/model/binding/LinearBinding.cpp index 0247d5cd4..d49f906dd 100644 --- a/src/libcadet/model/binding/LinearBinding.cpp +++ b/src/libcadet/model/binding/LinearBinding.cpp @@ -525,7 +525,7 @@ class LinearBindingBase : public IBindingModel { return fluxImpl(t, secIdx, colPos, y, yCp, res, workSpace); } - + /* virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, active const* yCp, active* res, LinearBufferAllocator workSpace, WithoutParamSensitivity) const { @@ -543,7 +543,7 @@ class LinearBindingBase : public IBindingModel { return fluxImpl(t, secIdx, colPos, y, yCp, res, workSpace); } - + */ virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const { jacobianImpl(t, secIdx, colPos, y, offsetCp, jac, workSpace); diff --git a/src/libcadet/model/exchange/LinearExchange.cpp b/src/libcadet/model/exchange/LinearExchange.cpp index 95d4f7cd3..6914086e3 100644 --- a/src/libcadet/model/exchange/LinearExchange.cpp +++ b/src/libcadet/model/exchange/LinearExchange.cpp @@ -21,6 +21,7 @@ #include "LocalVector.hpp" #include "SimulationTypes.hpp" #include "Memory.hpp" +#include "model/PhaseTransitionModel.hpp" #include #include @@ -360,15 +361,16 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase * See @cite Guiochon2006 * @tparam ParamHandler_t Type that can add support for external function dependence */ +// template class LinearExchangeBase : public IPhaseTransitionModel { public: - LinearExchangeBase() : _nComp(0), _nChannel(0) { } + LinearExchangeBase() : _nComp(0), _nChannel(0),_nBound(0) { } virtual ~LinearExchangeBase() CADET_NOEXCEPT { } - static const char* identifier() { return ParamHandler_t::identifier(); } - virtual const char* name() const CADET_NOEXCEPT { return ParamHandler_t::identifier(); } + static const char* identifier() { return "a"; } + virtual const char* name() const CADET_NOEXCEPT { return "a"; } virtual bool requiresConfiguration() const CADET_NOEXCEPT { return true; } virtual bool usesParamProviderInDiscretizationConfig() const CADET_NOEXCEPT { return true; } @@ -385,10 +387,10 @@ class LinearExchangeBase : public IPhaseTransitionModel _parameters.clear(); // Read parameters (k_a and k_d) - _paramHandler.configure(paramProvider, _nComp, _nChannel); + //_paramHandler.configure(paramProvider, _nComp, _nBound); // Register parameters - _paramHandler.registerParameters(_parameters, unitOpIdx, parTypeIdx, _nComp, _nChannel); + //_paramHandler.registerParameters(_parameters, unitOpIdx, parTypeIdx, _nComp, _nBound); return true; } @@ -407,7 +409,7 @@ class LinearExchangeBase : public IPhaseTransitionModel { std::unordered_map data; std::transform(_parameters.begin(), _parameters.end(), std::inserter(data, data.end()), - [](const std::pair& p) { return std::make_pair(p.first, static_cast(*p.second)); }); + [](const std::pair& p) { return std::make_pair(p.first, static_cast(*p.second)); }); return data; } @@ -439,7 +441,7 @@ class LinearExchangeBase : public IPhaseTransitionModel return false; } - virtual active* getParameter(const ParameterId& pId){ + virtual active* getParameter(const ParameterId& pId) { auto paramHandle = _parameters.find(pId); if (paramHandle != _parameters.end()) { @@ -451,80 +453,76 @@ class LinearExchangeBase : public IPhaseTransitionModel virtual unsigned int workspaceSize(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nChannel) const CADET_NOEXCEPT { - return _paramHandler.cacheSize(nComp, totalNumBoundStates, nChannel); + // return _paramHandler.cacheSize(nComp, totalNumBoundStates, nChannel); + return 0; } - virtual void setExternalFunctions(IExternalFunction** extFuns, unsigned int size) { _paramHandler.setExternalFunctions(extFuns, size); } + virtual void configure(IExternalFunction** extFuns, unsigned int size) {} // The next three flux() function implementations and two analyticJacobian() function // implementations are usually hidden behind // CADET_BINDINGMODELBASE_BOILERPLATE // which just expands to the six implementations below. - virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, - active const* y, active const* yCp, active* res, LinearBufferAllocator workSpace, WithParamSensitivity) const - { - return fluxImpl(t, secIdx, colPos, y, yCp, res, workSpace); - } - - virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, - active const* y, active const* yCp, active* res, LinearBufferAllocator workSpace, WithoutParamSensitivity) const + virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, + std::vector crossSections, active const* y, active* res) const { - return fluxImpl(t, secIdx, colPos, y, yCp, res, workSpace); + return fluxImpl( nChannel, nComp, nCol, exchangeMatrix, crossSections, y, res); } - - virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, - double const* y, double const* yCp, active* res, LinearBufferAllocator workSpace) const + /* + virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, + std::vector crossSections, double const* y, active* res) const { - return fluxImpl(t, secIdx, colPos, y, yCp, res, workSpace); + return fluxImpl(nChannel, nComp, nCol, exchangeMatrix, crossSections, y, res); } - virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, - double const* y, double const* yCp, double* res, LinearBufferAllocator workSpace) const + virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, + std::vector crossSections, double const* y, double* res) const { - return fluxImpl(t, secIdx, colPos, y, yCp, res, workSpace); + return fluxImpl(nChannel, nComp, nCol, exchangeMatrix, crossSections, y, res); } - - virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const + */ + virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, active const* y, active* res, linalg::CompressedSparseMatrix jac) const { - jacobianImpl(t, secIdx, colPos, y, offsetCp, jac, workSpace); + return jacobianImpl(nChannel, nComp, nCol, _exchangeMatrix, y, res, jac); } - - virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const + /* + virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, active* res, linalg::BandMatrix::RowIterator jac) const { - jacobianImpl(t, secIdx, colPos, y, offsetCp, jac, workSpace); + return jacobianImpl(nChannel, nComp, nCol, _exchangeMatrix, y, res, jac); } - virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const + virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, double* res, linalg::BandMatrix::RowIterator jac) const { - jacobianImpl(t, secIdx, colPos, y, offsetCp, jac, workSpace); + return jacobianImpl(nChannel, nComp, nCol, _exchangeMatrix, y, res, jac); } - + */ virtual void timeDerivativeQuasiStationaryFluxes(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yCp, double const* y, double* dResDt, LinearBufferAllocator workSpace) const { if (!hasQuasiStationaryReactions()) return; - if (!ParamHandler_t::dependsOnTime()) - return; + //if (!ParamHandler_t::dependsOnTime()) + // return; // Update external function and compute time derivative of parameters - typename ParamHandler_t::ParamsHandle p; - typename ParamHandler_t::ParamsHandle dpDt; - std::tie(p, dpDt) = _paramHandler.updateTimeDerivative(t, secIdx, colPos, _nComp, _nChannel, workSpace); + //typename ParamHandler_t::ParamsHandle p; + //typename ParamHandler_t::ParamsHandle dpDt; + /* + std::tie(p, dpDt) = _paramHandler.updateTimeDerivative(t, secIdx, colPos, _nComp, _nBound, workSpace); unsigned int bndIdx = 0; for (int i = 0; i < _nComp; ++i) { // Skip components without bound states (bound state index bndIdx is not advanced) - if (_nChannel[i] == 0) + if (_nBound[i] == 0) continue; dResDt[bndIdx] = -static_cast(dpDt->kA[i]) * yCp[i] + static_cast(dpDt->kD[i]) * y[bndIdx]; // Next bound component ++bndIdx; - } + }*/ } virtual bool preConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const @@ -552,85 +550,143 @@ class LinearExchangeBase : public IPhaseTransitionModel return std::any_of(_reactionQuasistationarity.begin(), _reactionQuasistationarity.end(), [](int i) -> bool { return !static_cast(i); }); } - virtual bool dependsOnTime() const CADET_NOEXCEPT { return ParamHandler_t::dependsOnTime(); } - virtual bool requiresWorkspace() const CADET_NOEXCEPT { return ParamHandler_t::requiresWorkspace(); } + virtual bool dependsOnTime() const CADET_NOEXCEPT { return false; } + virtual bool requiresWorkspace() const CADET_NOEXCEPT { return false; } virtual unsigned int requiredADdirs() const CADET_NOEXCEPT { return 0; } virtual int const* reactionQuasiStationarity() const CADET_NOEXCEPT { return _reactionQuasistationarity.data(); } protected: int _nComp; //!< Number of components - unsigned int const* _nChannel; //!< Array with number of bound states for each component + unsigned int const* _nBound; //!< Array with number of bound states for each component + unsigned int _nChannel; //!< Total number of bound states std::vector _reactionQuasistationarity; //!< Determines whether each bound state is quasi-stationary (@c true) or not (@c false) - ParamHandler_t _paramHandler; //!< Parameters + //ParamHandler_t _paramHandler; //!< Parameters std::unordered_map _parameters; //!< Map used to translate ParameterIds to actual variables virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; } - template - - int fluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, - StateType const* y, StateType const* yCp, ResidualType* res, LinearBufferAllocator workSpace) const - { - typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nChannel, workSpace); - - // Implement -k_a * c_{p,i} + k_d * q_i - unsigned int bndIdx = 0; - for (int i = 0; i < _nComp; ++i) - { - // Skip components without bound states (bound state index bndIdx is not advanced) - if (_nChannel[i] == 0) - continue; - - res[bndIdx] = -static_cast(p->kA[i]) * yCp[i] + static_cast(p->kD[i]) * y[bndIdx]; - - // Next bound component - ++bndIdx; + //template + int fluxImpl(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, active const* y, active* res) const + { + const unsigned int offsetC = nChannel * nComp; + for (unsigned int col = 0; col < nCol; ++col) { + + const unsigned int offsetColBlock = col * nChannel * nComp; + active* const resColBlock = res + offsetC + offsetColBlock; + active const* const yColBlock = y + offsetC + offsetColBlock; + + for (unsigned int rad_orig = 0; rad_orig < nChannel; ++rad_orig) + { + const unsigned int offsetToRadOrigBlock = rad_orig * nComp; + const unsigned int offsetColRadOrigBlock = offsetColBlock + offsetToRadOrigBlock; + active* const resColRadOrigBlock = resColBlock + offsetToRadOrigBlock; + active const* const yColRadOrigBlock = yColBlock + offsetToRadOrigBlock; + + for (unsigned int rad_dest = 0; rad_dest < nChannel; ++rad_dest) + { + if (rad_orig == rad_dest) + continue; + + const unsigned int offsetToRadDestBlock = rad_dest * nComp; + const unsigned int offsetColRadDestBlock = offsetColBlock + offsetToRadDestBlock; + active* const resColRadDestBlock = resColBlock + offsetToRadDestBlock; + // StateType const* const yColRadDestBlock = yColBlock + offsetToRadDestBlock; + + for (unsigned int comp = 0; comp < nComp; ++comp) + { + const unsigned int offsetCur_orig = offsetColRadOrigBlock + comp; + const unsigned int offsetCur_dest = offsetColRadDestBlock + comp; + active const* const yCur_orig = yColRadOrigBlock + comp; + // StateType const* const yCur_dest = yColRadDestBlock + comp; + active* const resCur_orig = resColRadOrigBlock + comp; + active* const resCur_dest = resColRadDestBlock + comp; + + const active exchange_orig_dest_comp = static_cast(exchangeMatrix[rad_orig * nChannel * nComp + rad_dest * nComp + comp]); + if (cadet_likely(exchange_orig_dest_comp > 0.0)) + { + *resCur_orig += exchange_orig_dest_comp * yCur_orig[0]; + *resCur_dest -= exchange_orig_dest_comp * yCur_orig[0] * static_cast(crossSections[rad_orig]) / static_cast(crossSections[rad_dest]); + } + } + + } + + } } return 0; } - template - inline void jacobianImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, RowIterator jac, LinearBufferAllocator workSpace) const + + //template + inline void jacobianImpl(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, active const* y, active* res, linalg::CompressedSparseMatrix jac) const { - typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nChannel, workSpace); - int bndIdx = 0; - for (int i = 0; i < _nComp; ++i) + const unsigned int offsetC = nChannel * nComp; + for (unsigned int col = 0; col < nCol; ++col) { - // Skip components without bound states (bound state index bndIdx is not advanced) - if (_nChannel[i] == 0) - continue; - - jac[0] = static_cast(p->kD[i]); // dres / dq_i - jac[i - bndIdx - offsetCp] = -static_cast(p->kA[i]); // dres / dc_{p,i} - // The distance from liquid phase to solid phase is reduced for each non-exchange component - // since a bound state is neglected. The number of neglected bound states so far is i - bndIdx. - // Thus, by going back offsetCp - (i - bndIdx) = -[ i - bndIdx - offsetCp ] we get to the corresponding - // liquid phase component. - - ++bndIdx; - ++jac; + const unsigned int offsetColBlock = col * nChannel * nComp; + active* const resColBlock = res + offsetC + offsetColBlock; + active const* const yColBlock = y + offsetC + offsetColBlock; + + for (unsigned int rad_orig = 0; rad_orig < nChannel; ++rad_orig) + { + const unsigned int offsetToRadOrigBlock = rad_orig * nComp; + const unsigned int offsetColRadOrigBlock = offsetColBlock + offsetToRadOrigBlock; + active* const resColRadOrigBlock = resColBlock + offsetToRadOrigBlock; + active const* const yColRadOrigBlock = yColBlock + offsetToRadOrigBlock; + + for (unsigned int rad_dest = 0; rad_dest < nChannel; ++rad_dest) + { + if (rad_orig == rad_dest) + continue; + + const unsigned int offsetToRadDestBlock = rad_dest * nComp; + const unsigned int offsetColRadDestBlock = offsetColBlock + offsetToRadDestBlock; + active* const resColRadDestBlock = resColBlock + offsetToRadDestBlock; + // StateType const* const yColRadDestBlock = yColBlock + offsetToRadDestBlock; + + for (unsigned int comp = 0; comp < nComp; ++comp) + { + const unsigned int offsetCur_orig = offsetColRadOrigBlock + comp; + const unsigned int offsetCur_dest = offsetColRadDestBlock + comp; + active const* const yCur_orig = yColRadOrigBlock + comp; + // StateType const* const yCur_dest = yColRadDestBlock + comp; + active* const resCur_orig = resColRadOrigBlock + comp; + active* const resCur_dest = resColRadDestBlock + comp; + + const active exchange_orig_dest_comp = static_cast(exchangeMatrix[rad_orig * nChannel * nComp + rad_dest * nComp + comp]); + if (cadet_likely(exchange_orig_dest_comp > 0.0)) + { + jac.centered(offsetCur_orig, 0) += static_cast(exchange_orig_dest_comp); + jac.centered(offsetCur_dest, static_cast(offsetCur_orig) - static_cast(offsetCur_dest)) -= static_cast(exchange_orig_dest_comp); + } + } + + } + + } } - } + } }; -typedef LinearExchangeBase LinearExchange; -typedef LinearExchangeBase ExternalLinearExchange; + +//typedef LinearExchangeBase LinearExchange; +// typedef LinearExchangeBase ExternalLinearExchange; namespace exchange { void registerLinearModel(std::unordered_map>& exchange) { - exchange[LinearExchange::identifier()] = []() { return new LinearExchange(); }; - exchange[ExternalLinearExchange::identifier()] = []() { return new ExternalLinearExchange(); }; + //exchange[LinearExchange::identifier()] = []() { return new LinearExchange(); }; + //exchange[ExternalLinearExchange::identifier()] = []() { return new ExternalLinearExchange(); }; } } // namespace exchange -} // namespace model +} // namespace model -} // namespace cadet +}// namespace cadet diff --git a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp index c010e44a0..82e675a3a 100644 --- a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp +++ b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp @@ -21,6 +21,8 @@ #include "linalg/CompressedSparseMatrix.hpp" #include "model/parts/AxialConvectionDispersionKernel.hpp" #include "model/ParameterDependence.hpp" +#include "model/exchange/LinearExchange.cpp" + #ifdef SUPERLU_FOUND #include "linalg/SuperLUSparseMatrix.hpp" @@ -830,6 +832,8 @@ bool MultiChannelConvectionDispersionOperator::configure(UnitOpIdx unitOpIdx, IP setSparsityPattern(); + _phaseTransitionModel = new LinearExchangeBase(); + return true; } @@ -936,6 +940,7 @@ const active& MultiChannelConvectionDispersionOperator::axialDispersion(unsigned * @param [in] wantJac Determines whether the Jacobian is computed or not * @return @c 0 on success, @c -1 on non-recoverable error, and @c +1 on recoverable error */ +/* int MultiChannelConvectionDispersionOperator::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, double* res, bool wantJac, WithoutParamSensitivity) { if (wantJac) @@ -943,15 +948,15 @@ int MultiChannelConvectionDispersionOperator::residual(const IModel& model, doub else return residualImpl(model, t, secIdx, y, yDot, res); } - +*/ int MultiChannelConvectionDispersionOperator::residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithoutParamSensitivity) { if (wantJac) - return residualImpl(model, t, secIdx, y, yDot, res); + return residualImpl(model, t, secIdx, y, yDot, res); else - return residualImpl(model, t, secIdx, y, yDot, res); + return residualImpl(model, t, secIdx, y, yDot, res); } - +/* int MultiChannelConvectionDispersionOperator::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity) { if (wantJac) @@ -967,7 +972,7 @@ int MultiChannelConvectionDispersionOperator::residual(const IModel& model, doub else return residualImpl(model, t, secIdx, y, yDot, res); } - +*/ template int MultiChannelConvectionDispersionOperator::residualImpl(const IModel& model, double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res) { @@ -1005,60 +1010,16 @@ int MultiChannelConvectionDispersionOperator::residualImpl(const IModel& model, convdisp::residualKernelAxial(SimulationTime{t, secIdx}, y, yDot, res, _jacC.row(i * _nComp), fp); } + // Handle inter-channel transport if (cadet_unlikely(_nChannel <= 1)) return 0; - const unsigned int offsetC = _nChannel * _nComp; - for (unsigned int col = 0; col < _nCol; ++col) - { - const unsigned int offsetColBlock = col * _nChannel * _nComp; - ResidualType* const resColBlock = res + offsetC + offsetColBlock; - StateType const* const yColBlock = y + offsetC + offsetColBlock; - for (unsigned int rad_orig = 0; rad_orig < _nChannel; ++rad_orig) - { - const unsigned int offsetToRadOrigBlock = rad_orig * _nComp; - const unsigned int offsetColRadOrigBlock = offsetColBlock + offsetToRadOrigBlock; - ResidualType* const resColRadOrigBlock = resColBlock + offsetToRadOrigBlock; - StateType const* const yColRadOrigBlock = yColBlock + offsetToRadOrigBlock; - - for (unsigned int rad_dest = 0; rad_dest < _nChannel; ++rad_dest) - { - if (rad_orig == rad_dest) - continue; - - const unsigned int offsetToRadDestBlock = rad_dest * _nComp; - const unsigned int offsetColRadDestBlock = offsetColBlock + offsetToRadDestBlock; - ResidualType* const resColRadDestBlock = resColBlock + offsetToRadDestBlock; - // StateType const* const yColRadDestBlock = yColBlock + offsetToRadDestBlock; - - for (unsigned int comp = 0; comp < _nComp; ++comp) - { - const unsigned int offsetCur_orig = offsetColRadOrigBlock + comp; - const unsigned int offsetCur_dest = offsetColRadDestBlock + comp; - StateType const* const yCur_orig = yColRadOrigBlock + comp; - // StateType const* const yCur_dest = yColRadDestBlock + comp; - ResidualType* const resCur_orig = resColRadOrigBlock + comp; - ResidualType* const resCur_dest = resColRadDestBlock + comp; - - const ParamType exchange_orig_dest_comp = static_cast(_exchangeMatrix[rad_orig * _nChannel * _nComp + rad_dest * _nComp + comp]); - if (cadet_likely(exchange_orig_dest_comp > 0.0)) - { - *resCur_orig += exchange_orig_dest_comp * yCur_orig[0]; - *resCur_dest -= exchange_orig_dest_comp * yCur_orig[0] * static_cast(_crossSections[rad_orig]) / static_cast(_crossSections[rad_dest]); - - if (wantJac) - { - _jacC.centered(offsetCur_orig, 0) += static_cast(exchange_orig_dest_comp); - _jacC.centered(offsetCur_dest, static_cast(offsetCur_orig) - static_cast(offsetCur_dest)) -= static_cast(exchange_orig_dest_comp); - } - } - } - - } + _phaseTransitionModel->flux(_nChannel, _nComp, _nCol, _exchangeMatrix, _crossSections, y, res); - } + if (wantJac){ + _phaseTransitionModel->analyticJacobian(_nChannel, _nComp, _nCol, _exchangeMatrix, y, res, _jacC); } return 0; diff --git a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.hpp b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.hpp index bf947429c..8fe8a0369 100644 --- a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.hpp +++ b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.hpp @@ -26,6 +26,7 @@ #include "model/ParameterMultiplexing.hpp" #include "SimulationTypes.hpp" #include "ConfigurationHelper.hpp" +#include "model/PhaseTransitionModel.hpp" #include #include @@ -146,6 +147,7 @@ class MultiChannelConvectionDispersionOperator unsigned int _nChannel; //!< Number of channels bool _hasDynamicReactions; //!< Determines whether the model has dynamic reactions (only relevant for sparsity pattern) + active _colLength; //!< Column length \f$ L \f$ std::vector _crossSections; //!< Cross section area of each compartment @@ -160,6 +162,8 @@ class MultiChannelConvectionDispersionOperator std::vector _exchangeMatrix; //!< Matrix of exchange coeffs for the linear inter-channel transport + IPhaseTransitionModel* _phaseTransitionModel; //!< Phase transition model + IParameterParameterDependence* _dispersionDep; ArrayPool _stencilMemory; //!< Provides memory for the stencil From ea5b80cb35d795901c1de432a56cb2de682c9ab9 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Wed, 7 Aug 2024 11:40:43 +0200 Subject: [PATCH 06/17] Fix template configuration --- src/libcadet/model/PhaseTransitionModel.hpp | 15 ++--- src/libcadet/model/binding/LinearBinding.cpp | 4 +- .../model/exchange/LinearExchange.cpp | 67 ++++++++++--------- ...ltiChannelConvectionDispersionOperator.cpp | 17 +++-- 4 files changed, 50 insertions(+), 53 deletions(-) diff --git a/src/libcadet/model/PhaseTransitionModel.hpp b/src/libcadet/model/PhaseTransitionModel.hpp index 56563dd3e..736fd3306 100644 --- a/src/libcadet/model/PhaseTransitionModel.hpp +++ b/src/libcadet/model/PhaseTransitionModel.hpp @@ -96,15 +96,12 @@ class IPhaseTransitionModel //virtual unsigned int requiredADdirs() const CADET_NOEXCEPT = 0; - virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, active const* y, active* res) const = 0; - //virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, active const* y, double* res) const = 0; - //virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, double const* y, active* res) const = 0; - //virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, double const* y, double* res) const = 0; - - virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, active const* y, active* res, linalg::CompressedSparseMatrix jac) const = 0; - //virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, active const* y, double* res, linalg::BandMatrix::RowIterator jac) const = 0; - //virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, active* res, linalg::BandMatrix::RowIterator jac) const = 0; - // virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, double* res, linalg::BandMatrix::RowIterator jac) const = 0; + virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, active const* y, active* res, WithParamSensitivity) const = 0; + virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, active const* y, active* res, WithoutParamSensitivity) const = 0; + virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, double const* y, active* res, WithParamSensitivity) const = 0; + virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, double const* y, double* res, WithoutParamSensitivity) const = 0; + + virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, linalg::CompressedSparseMatrix jac) const = 0; #ifdef ENABLE_DG //virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; #endif diff --git a/src/libcadet/model/binding/LinearBinding.cpp b/src/libcadet/model/binding/LinearBinding.cpp index d49f906dd..048cf0b2b 100644 --- a/src/libcadet/model/binding/LinearBinding.cpp +++ b/src/libcadet/model/binding/LinearBinding.cpp @@ -525,7 +525,7 @@ class LinearBindingBase : public IBindingModel { return fluxImpl(t, secIdx, colPos, y, yCp, res, workSpace); } - /* + virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, active const* yCp, active* res, LinearBufferAllocator workSpace, WithoutParamSensitivity) const { @@ -543,7 +543,7 @@ class LinearBindingBase : public IBindingModel { return fluxImpl(t, secIdx, colPos, y, yCp, res, workSpace); } - */ + virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const { jacobianImpl(t, secIdx, colPos, y, offsetCp, jac, workSpace); diff --git a/src/libcadet/model/exchange/LinearExchange.cpp b/src/libcadet/model/exchange/LinearExchange.cpp index 6914086e3..01af45d57 100644 --- a/src/libcadet/model/exchange/LinearExchange.cpp +++ b/src/libcadet/model/exchange/LinearExchange.cpp @@ -465,26 +465,32 @@ class LinearExchangeBase : public IPhaseTransitionModel // which just expands to the six implementations below. virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, - std::vector crossSections, active const* y, active* res) const + std::vector crossSections, active const* y, active* res, WithParamSensitivity) const { - return fluxImpl( nChannel, nComp, nCol, exchangeMatrix, crossSections, y, res); + return fluxImpl( nChannel, nComp, nCol, exchangeMatrix, crossSections, y, res); } - /* + virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, - std::vector crossSections, double const* y, active* res) const + std::vector crossSections, active const* y, active* res, WithoutParamSensitivity) const + { + return fluxImpl(nChannel, nComp, nCol, exchangeMatrix, crossSections, y, res); + } + + virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, + std::vector crossSections, double const* y, active* res, WithParamSensitivity) const { return fluxImpl(nChannel, nComp, nCol, exchangeMatrix, crossSections, y, res); } virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, - std::vector crossSections, double const* y, double* res) const + std::vector crossSections, double const* y, double* res, WithoutParamSensitivity) const { - return fluxImpl(nChannel, nComp, nCol, exchangeMatrix, crossSections, y, res); + return fluxImpl(nChannel, nComp, nCol, exchangeMatrix, crossSections, y, res); } - */ - virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, active const* y, active* res, linalg::CompressedSparseMatrix jac) const + + virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, linalg::CompressedSparseMatrix jac) const { - return jacobianImpl(nChannel, nComp, nCol, _exchangeMatrix, y, res, jac); + return jacobianImpl(nChannel, nComp, nCol, _exchangeMatrix, y, jac); } /* virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, active* res, linalg::BandMatrix::RowIterator jac) const @@ -568,22 +574,22 @@ class LinearExchangeBase : public IPhaseTransitionModel virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; } - //template - int fluxImpl(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, active const* y, active* res) const + template + int fluxImpl(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, StateType const* y, ResidualType* res) const { const unsigned int offsetC = nChannel * nComp; - for (unsigned int col = 0; col < nCol; ++col) { - + for (unsigned int col = 0; col < nCol; ++col) + { const unsigned int offsetColBlock = col * nChannel * nComp; - active* const resColBlock = res + offsetC + offsetColBlock; - active const* const yColBlock = y + offsetC + offsetColBlock; + ResidualType* const resColBlock = res + offsetC + offsetColBlock; + StateType const* const yColBlock = y + offsetC + offsetColBlock; for (unsigned int rad_orig = 0; rad_orig < nChannel; ++rad_orig) { - const unsigned int offsetToRadOrigBlock = rad_orig * nComp; + const unsigned int offsetToRadOrigBlock = rad_orig * _nComp; const unsigned int offsetColRadOrigBlock = offsetColBlock + offsetToRadOrigBlock; - active* const resColRadOrigBlock = resColBlock + offsetToRadOrigBlock; - active const* const yColRadOrigBlock = yColBlock + offsetToRadOrigBlock; + ResidualType* const resColRadOrigBlock = resColBlock + offsetToRadOrigBlock; + StateType const* const yColRadOrigBlock = yColBlock + offsetToRadOrigBlock; for (unsigned int rad_dest = 0; rad_dest < nChannel; ++rad_dest) { @@ -592,23 +598,23 @@ class LinearExchangeBase : public IPhaseTransitionModel const unsigned int offsetToRadDestBlock = rad_dest * nComp; const unsigned int offsetColRadDestBlock = offsetColBlock + offsetToRadDestBlock; - active* const resColRadDestBlock = resColBlock + offsetToRadDestBlock; + ResidualType* const resColRadDestBlock = resColBlock + offsetToRadDestBlock; // StateType const* const yColRadDestBlock = yColBlock + offsetToRadDestBlock; for (unsigned int comp = 0; comp < nComp; ++comp) { const unsigned int offsetCur_orig = offsetColRadOrigBlock + comp; const unsigned int offsetCur_dest = offsetColRadDestBlock + comp; - active const* const yCur_orig = yColRadOrigBlock + comp; + StateType const* const yCur_orig = yColRadOrigBlock + comp; // StateType const* const yCur_dest = yColRadDestBlock + comp; - active* const resCur_orig = resColRadOrigBlock + comp; - active* const resCur_dest = resColRadDestBlock + comp; + ResidualType* const resCur_orig = resColRadOrigBlock + comp; + ResidualType* const resCur_dest = resColRadDestBlock + comp; - const active exchange_orig_dest_comp = static_cast(exchangeMatrix[rad_orig * nChannel * nComp + rad_dest * nComp + comp]); + const ParamType exchange_orig_dest_comp = static_cast(exchangeMatrix[rad_orig * nChannel * nComp + rad_dest * nComp + comp]); if (cadet_likely(exchange_orig_dest_comp > 0.0)) { *resCur_orig += exchange_orig_dest_comp * yCur_orig[0]; - *resCur_dest -= exchange_orig_dest_comp * yCur_orig[0] * static_cast(crossSections[rad_orig]) / static_cast(crossSections[rad_dest]); + *resCur_dest -= exchange_orig_dest_comp * yCur_orig[0] * static_cast(crossSections[rad_orig]) / static_cast(crossSections[rad_dest]); } } @@ -622,22 +628,20 @@ class LinearExchangeBase : public IPhaseTransitionModel //template - inline void jacobianImpl(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, active const* y, active* res, linalg::CompressedSparseMatrix jac) const + inline void jacobianImpl(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, double const* y, linalg::CompressedSparseMatrix jac) const { const unsigned int offsetC = nChannel * nComp; for (unsigned int col = 0; col < nCol; ++col) { const unsigned int offsetColBlock = col * nChannel * nComp; - active* const resColBlock = res + offsetC + offsetColBlock; - active const* const yColBlock = y + offsetC + offsetColBlock; + double const* const yColBlock = y + offsetC + offsetColBlock; for (unsigned int rad_orig = 0; rad_orig < nChannel; ++rad_orig) { const unsigned int offsetToRadOrigBlock = rad_orig * nComp; const unsigned int offsetColRadOrigBlock = offsetColBlock + offsetToRadOrigBlock; - active* const resColRadOrigBlock = resColBlock + offsetToRadOrigBlock; - active const* const yColRadOrigBlock = yColBlock + offsetToRadOrigBlock; + double const* const yColRadOrigBlock = yColBlock + offsetToRadOrigBlock; for (unsigned int rad_dest = 0; rad_dest < nChannel; ++rad_dest) { @@ -646,17 +650,14 @@ class LinearExchangeBase : public IPhaseTransitionModel const unsigned int offsetToRadDestBlock = rad_dest * nComp; const unsigned int offsetColRadDestBlock = offsetColBlock + offsetToRadDestBlock; - active* const resColRadDestBlock = resColBlock + offsetToRadDestBlock; // StateType const* const yColRadDestBlock = yColBlock + offsetToRadDestBlock; for (unsigned int comp = 0; comp < nComp; ++comp) { const unsigned int offsetCur_orig = offsetColRadOrigBlock + comp; const unsigned int offsetCur_dest = offsetColRadDestBlock + comp; - active const* const yCur_orig = yColRadOrigBlock + comp; + double const* const yCur_orig = yColRadOrigBlock + comp; // StateType const* const yCur_dest = yColRadDestBlock + comp; - active* const resCur_orig = resColRadOrigBlock + comp; - active* const resCur_dest = resColRadDestBlock + comp; const active exchange_orig_dest_comp = static_cast(exchangeMatrix[rad_orig * nChannel * nComp + rad_dest * nComp + comp]); if (cadet_likely(exchange_orig_dest_comp > 0.0)) diff --git a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp index 82e675a3a..c24fd184d 100644 --- a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp +++ b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp @@ -940,7 +940,7 @@ const active& MultiChannelConvectionDispersionOperator::axialDispersion(unsigned * @param [in] wantJac Determines whether the Jacobian is computed or not * @return @c 0 on success, @c -1 on non-recoverable error, and @c +1 on recoverable error */ -/* + int MultiChannelConvectionDispersionOperator::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, double* res, bool wantJac, WithoutParamSensitivity) { if (wantJac) @@ -948,15 +948,15 @@ int MultiChannelConvectionDispersionOperator::residual(const IModel& model, doub else return residualImpl(model, t, secIdx, y, yDot, res); } -*/ + int MultiChannelConvectionDispersionOperator::residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithoutParamSensitivity) { if (wantJac) - return residualImpl(model, t, secIdx, y, yDot, res); + return residualImpl(model, t, secIdx, y, yDot, res); else - return residualImpl(model, t, secIdx, y, yDot, res); + return residualImpl(model, t, secIdx, y, yDot, res); } -/* + int MultiChannelConvectionDispersionOperator::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity) { if (wantJac) @@ -972,7 +972,7 @@ int MultiChannelConvectionDispersionOperator::residual(const IModel& model, doub else return residualImpl(model, t, secIdx, y, yDot, res); } -*/ + template int MultiChannelConvectionDispersionOperator::residualImpl(const IModel& model, double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res) { @@ -1015,11 +1015,10 @@ int MultiChannelConvectionDispersionOperator::residualImpl(const IModel& model, if (cadet_unlikely(_nChannel <= 1)) return 0; - - _phaseTransitionModel->flux(_nChannel, _nComp, _nCol, _exchangeMatrix, _crossSections, y, res); + _phaseTransitionModel->flux(_nChannel, _nComp, _nCol, _exchangeMatrix, _crossSections, y, res, typename ParamSens::enabled()); if (wantJac){ - _phaseTransitionModel->analyticJacobian(_nChannel, _nComp, _nCol, _exchangeMatrix, y, res, _jacC); + _phaseTransitionModel->analyticJacobian(_nChannel, _nComp, _nCol, _exchangeMatrix, reinterpret_cast(y), _jacC); } return 0; From 11ef2f3fba74cb7e5c2b5b2a449378257cc107b7 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Wed, 7 Aug 2024 16:05:00 +0200 Subject: [PATCH 07/17] changed parameter --- .../model/MultiChannelTransportModel.hpp | 6 ++++++ src/libcadet/model/PhaseTransitionModel.hpp | 2 +- src/libcadet/model/exchange/LinearExchange.cpp | 16 +++++++++++----- .../MultiChannelConvectionDispersionOperator.cpp | 2 +- 4 files changed, 19 insertions(+), 7 deletions(-) diff --git a/src/libcadet/model/MultiChannelTransportModel.hpp b/src/libcadet/model/MultiChannelTransportModel.hpp index aa163694c..d8ee9c531 100644 --- a/src/libcadet/model/MultiChannelTransportModel.hpp +++ b/src/libcadet/model/MultiChannelTransportModel.hpp @@ -204,8 +204,14 @@ class MultiChannelTransportModel : public UnitOperationBase parts::MultiChannelConvectionDispersionOperator _convDispOp; //!< Convection dispersion operator for interstitial volume transport IDynamicReactionModel* _dynReactionBulk; //!< Dynamic reactions in the bulk volume + IPhaseTransitionModel* _phaseTransModel; //!< Phase transition model linalg::DoubleSparseMatrix _jacInlet; //!< Jacobian inlet DOF block matrix connects inlet DOFs to first bulk cells + linalg::CompressedSparseMatrix _jacC; //!< Jacobian + + std::vector _exchangeMatrix; //!< Matrix of exchange coeffs for the linear inter-channel transport + std::vector _crossSections; //!< Cross section area of each compartment + bool _analyticJac; //!< Determines whether AD or analytic Jacobians are used unsigned int _jacobianAdDirs; //!< Number of AD seed vectors required for Jacobian computation diff --git a/src/libcadet/model/PhaseTransitionModel.hpp b/src/libcadet/model/PhaseTransitionModel.hpp index 736fd3306..ef849647f 100644 --- a/src/libcadet/model/PhaseTransitionModel.hpp +++ b/src/libcadet/model/PhaseTransitionModel.hpp @@ -101,7 +101,7 @@ class IPhaseTransitionModel virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, double const* y, active* res, WithParamSensitivity) const = 0; virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, double const* y, double* res, WithoutParamSensitivity) const = 0; - virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, linalg::CompressedSparseMatrix jac) const = 0; + virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, linalg::BandedSparseRowIterator jac) const = 0; #ifdef ENABLE_DG //virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; #endif diff --git a/src/libcadet/model/exchange/LinearExchange.cpp b/src/libcadet/model/exchange/LinearExchange.cpp index 01af45d57..57a510756 100644 --- a/src/libcadet/model/exchange/LinearExchange.cpp +++ b/src/libcadet/model/exchange/LinearExchange.cpp @@ -488,7 +488,7 @@ class LinearExchangeBase : public IPhaseTransitionModel return fluxImpl(nChannel, nComp, nCol, exchangeMatrix, crossSections, y, res); } - virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, linalg::CompressedSparseMatrix jac) const + virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, linalg::BandedSparseRowIterator jac) const { return jacobianImpl(nChannel, nComp, nCol, _exchangeMatrix, y, jac); } @@ -628,7 +628,7 @@ class LinearExchangeBase : public IPhaseTransitionModel //template - inline void jacobianImpl(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, double const* y, linalg::CompressedSparseMatrix jac) const + inline void jacobianImpl(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, double const* y, linalg::BandedSparseRowIterator jacBegin) const { const unsigned int offsetC = nChannel * nComp; @@ -660,10 +660,16 @@ class LinearExchangeBase : public IPhaseTransitionModel // StateType const* const yCur_dest = yColRadDestBlock + comp; const active exchange_orig_dest_comp = static_cast(exchangeMatrix[rad_orig * nChannel * nComp + rad_dest * nComp + comp]); + if (cadet_likely(exchange_orig_dest_comp > 0.0)) - { - jac.centered(offsetCur_orig, 0) += static_cast(exchange_orig_dest_comp); - jac.centered(offsetCur_dest, static_cast(offsetCur_orig) - static_cast(offsetCur_dest)) -= static_cast(exchange_orig_dest_comp); + { + linalg::BandedSparseRowIterator jacorig; + jacorig = jacBegin + offsetCur_orig; + jacorig[0] += static_cast(exchange_orig_dest_comp); + + linalg::BandedSparseRowIterator jacdest; + jacdest = jacBegin + offsetCur_dest; + jacdest[static_cast(offsetCur_orig) - static_cast(offsetCur_dest)] -= static_cast(exchange_orig_dest_comp); } } diff --git a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp index c24fd184d..0b1355e09 100644 --- a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp +++ b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp @@ -1018,7 +1018,7 @@ int MultiChannelConvectionDispersionOperator::residualImpl(const IModel& model, _phaseTransitionModel->flux(_nChannel, _nComp, _nCol, _exchangeMatrix, _crossSections, y, res, typename ParamSens::enabled()); if (wantJac){ - _phaseTransitionModel->analyticJacobian(_nChannel, _nComp, _nCol, _exchangeMatrix, reinterpret_cast(y), _jacC); + _phaseTransitionModel->analyticJacobian(_nChannel, _nComp, _nCol, _exchangeMatrix, reinterpret_cast(y), _jacC.row(0)); } return 0; From b138a934c5b8293c95dca70381c6bd3daa4bc505 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Thu, 8 Aug 2024 12:11:06 +0200 Subject: [PATCH 08/17] fix false parameter declaration --- src/libcadet/model/MultiChannelTransportModel.cpp | 4 ++-- src/libcadet/model/exchange/LinearExchange.cpp | 8 +++++--- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/libcadet/model/MultiChannelTransportModel.cpp b/src/libcadet/model/MultiChannelTransportModel.cpp index 11d28aa65..d997a0ac4 100644 --- a/src/libcadet/model/MultiChannelTransportModel.cpp +++ b/src/libcadet/model/MultiChannelTransportModel.cpp @@ -411,10 +411,10 @@ bool MultiChannelTransportModel::configureModelDiscretization(IParameterProvider bool exchangeConfSuccess = true; //_exchange = nullptr - std::vector exchModelNames = { "NONE" }; + /*std::vector exchModelNames = {"NONE"}; if (paramProvider.exists("EXCHANGE_MODEL")) exchModelNames = paramProvider.getStringArray("EXCHANGE_MODEL"); - + */ bool exchConfSuccess = true; //exchConfSuccess = _exchange->configureModelDiscretization(paramProvider, _disc.nComp, _disc.nChannel , _disc.boundOffset) && bindingConfSuccess; diff --git a/src/libcadet/model/exchange/LinearExchange.cpp b/src/libcadet/model/exchange/LinearExchange.cpp index 57a510756..12cdf02bc 100644 --- a/src/libcadet/model/exchange/LinearExchange.cpp +++ b/src/libcadet/model/exchange/LinearExchange.cpp @@ -366,7 +366,7 @@ class LinearExchangeBase : public IPhaseTransitionModel { public: - LinearExchangeBase() : _nComp(0), _nChannel(0),_nBound(0) { } + LinearExchangeBase() : _nComp(0), _nChannel(0),_nBound(0), _nCol(0) { } virtual ~LinearExchangeBase() CADET_NOEXCEPT { } static const char* identifier() { return "a"; } @@ -374,10 +374,11 @@ class LinearExchangeBase : public IPhaseTransitionModel virtual bool requiresConfiguration() const CADET_NOEXCEPT { return true; } virtual bool usesParamProviderInDiscretizationConfig() const CADET_NOEXCEPT { return true; } - virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int nChannel, unsigned int const* boundOffset) + virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int nChannel, unsigned int nCol) { _nComp = nComp; _nChannel = nChannel; // nChannel realy int ? -> by not nBoundStates not ? + _nCol = nCol; return true; } @@ -565,6 +566,7 @@ class LinearExchangeBase : public IPhaseTransitionModel int _nComp; //!< Number of components unsigned int const* _nBound; //!< Array with number of bound states for each component unsigned int _nChannel; //!< Total number of bound states + unsigned int _nCol; //!< Number of columns std::vector _reactionQuasistationarity; //!< Determines whether each bound state is quasi-stationary (@c true) or not (@c false) //ParamHandler_t _paramHandler; //!< Parameters @@ -586,7 +588,7 @@ class LinearExchangeBase : public IPhaseTransitionModel for (unsigned int rad_orig = 0; rad_orig < nChannel; ++rad_orig) { - const unsigned int offsetToRadOrigBlock = rad_orig * _nComp; + const unsigned int offsetToRadOrigBlock = rad_orig * nComp; const unsigned int offsetColRadOrigBlock = offsetColBlock + offsetToRadOrigBlock; ResidualType* const resColRadOrigBlock = resColBlock + offsetToRadOrigBlock; StateType const* const yColRadOrigBlock = yColBlock + offsetToRadOrigBlock; From 68324639c4f64d19719f0e6b1e8a7b7988abc7da Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Thu, 8 Aug 2024 18:12:53 +0200 Subject: [PATCH 09/17] Add linear flux residual to MCT --- src/libcadet/CMakeLists.txt | 1 + src/libcadet/ConfigurationHelper.hpp | 4 ++ src/libcadet/ExchangeModelFactory.cpp | 75 ++++++++++++++++++++++ src/libcadet/ExchangeModelFactory.hpp | 90 +++++++++++++++++++++++++++ src/libcadet/ModelBuilderImpl.cpp | 5 ++ src/libcadet/ModelBuilderImpl.hpp | 3 + test/Dummies.hpp | 1 + 7 files changed, 179 insertions(+) create mode 100644 src/libcadet/ExchangeModelFactory.cpp create mode 100644 src/libcadet/ExchangeModelFactory.hpp diff --git a/src/libcadet/CMakeLists.txt b/src/libcadet/CMakeLists.txt index d1f477105..31c0653d2 100644 --- a/src/libcadet/CMakeLists.txt +++ b/src/libcadet/CMakeLists.txt @@ -128,6 +128,7 @@ set(LIBCADET_SOURCES ${CMAKE_SOURCE_DIR}/src/libcadet/AdUtils.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/Weno.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/BindingModelFactory.cpp + ${CMAKE_SOURCE_DIR}/src/libcadet/ExchangeModelFactory.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/ReactionModelFactory.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/ParameterDependenceFactory.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/graph/GraphAlgos.cpp diff --git a/src/libcadet/ConfigurationHelper.hpp b/src/libcadet/ConfigurationHelper.hpp index f06d09a53..b26b3c297 100644 --- a/src/libcadet/ConfigurationHelper.hpp +++ b/src/libcadet/ConfigurationHelper.hpp @@ -30,6 +30,7 @@ class IExternalFunction; namespace model { class IBindingModel; + class IPhaseTransitionModel; class IDynamicReactionModel; class IParameterStateDependence; class IParameterParameterDependence; @@ -58,6 +59,9 @@ class IConfigHelper */ virtual model::IBindingModel* createBindingModel(const std::string& name) const = 0; + virtual model::IPhaseTransitionModel* createExchangeModel(const std::string& name) const = 0; + + /** * @brief Checks if there is an IBindingModel of the given @p name * @param [in] name Name of the IBindingModel object diff --git a/src/libcadet/ExchangeModelFactory.cpp b/src/libcadet/ExchangeModelFactory.cpp new file mode 100644 index 000000000..f7b521ca0 --- /dev/null +++ b/src/libcadet/ExchangeModelFactory.cpp @@ -0,0 +1,75 @@ +// ============================================================================= +// CADET +// +// Copyright © 2008-2022: The CADET Authors +// Please see the AUTHORS and CONTRIBUTORS file. +// +// All rights reserved. This program and the accompanying materials +// are made available under the terms of the GNU Public License v3.0 (or, at +// your option, any later version) which accompanies this distribution, and +// is available at http://www.gnu.org/licenses/gpl.html +// ============================================================================= + +#include "ExchangeModelFactory.hpp" +#include "cadet/Exceptions.hpp" + +#include "model/binding/SimplifiedMultiStateStericMassActionBinding.hpp" + +namespace cadet +{ + namespace model + { + namespace exchange + { + void registerLinearExModel(std::unordered_map>& exchange); + } + } + + ExchangeModelFactory::ExchangeModelFactory() + { + // Register all ExchangeModels here + model::exchange::registerLinearExModel(_exchangeModels); + + } + + ExchangeModelFactory::~ExchangeModelFactory() { } + + template + void ExchangeModelFactory::registerModel(const std::string& name) + { + _exchangeModels[name] = []() { return new ExchangeModel_t(); }; + } + + template + void ExchangeModelFactory::registerModel() + { + registerModel(ExchangeModel_t::identifier()); + } + + model::IPhaseTransitionModel* ExchangeModelFactory::create(const std::string& name) const + { + const auto it = _exchangeModels.find(name); + if (it == _exchangeModels.end()) + { + // BindingModel was not found + return nullptr; + } + + // Call factory function (thanks to type erasure of std::function we can store + // all factory functions in one container) + return it->second(); + } + + void ExchangeModelFactory::registerModel(const std::string& name, std::function factory) + { + if (_exchangeModels.find(name) == _exchangeModels.end()) + _exchangeModels[name] = factory; + else + throw InvalidParameterException("IPhaseTransition implementation with the name " + name + " is already registered and cannot be overwritten"); + } + + bool ExchangeModelFactory::exists(const std::string& name) const + { + return _exchangeModels.find(name) != _exchangeModels.end(); + } +} // namespace cadet diff --git a/src/libcadet/ExchangeModelFactory.hpp b/src/libcadet/ExchangeModelFactory.hpp new file mode 100644 index 000000000..16b2568a3 --- /dev/null +++ b/src/libcadet/ExchangeModelFactory.hpp @@ -0,0 +1,90 @@ +// ============================================================================= +// CADET +// +// Copyright © 2008-2024: The CADET Authors +// Please see the AUTHORS and CONTRIBUTORS file. +// +// All rights reserved. This program and the accompanying materials +// are made available under the terms of the GNU Public License v3.0 (or, at +// your option, any later version) which accompanies this distribution, and +// is available at http://www.gnu.org/licenses/gpl.html +// ============================================================================= + +/** + * @file + * Defines the ExchangeModelFactory + */ + +#ifndef LIBCADET_EXCHANGEMODELFACTORY_HPP_ +#define LIBCADET_EXCHANGEMODELFACTORY_HPP_ + +#include +#include +#include + +namespace cadet +{ + + namespace model + { + class IPhaseTransitionModel; + } + + /** + * @brief Creates binding models + */ + class ExchangeModelFactory + { + public: + /** + * @brief Construct the BindingModelFactory + * @details All internal binding models are registered here. + */ + ExchangeModelFactory(); + + ~ExchangeModelFactory(); + + /** + * @brief Creates binding models with the given @p name + * @param [in] name Name of the binding model + * @return The binding model or @c NULL if a binding model with this name does not exist + */ + model::IPhaseTransitionModel* create(const std::string& name) const; + + /** + * @brief Registers the given binding model implementation + * @param [in] name Name of the IBindingModel implementation + * @param [in] factory Function that creates an object of the IBindingModel class + */ + void registerModel(const std::string& name, std::function factory); + + /** + * @brief Returns whether a binding model of the given name @p name exists + * @param [in] name Name of the binding model + * @return @c true if a binding model of this name exists, otherwise @c false + */ + bool exists(const std::string& name) const; + protected: + + /** + * @brief Registers an IBindingModel + * @param [in] name Name of the binding model + * @tparam BindingModel_t Type of the binding model + */ + template + void registerModel(const std::string& name); + + /** + * @brief Registers an IBindingModel + * @details The name of the binding model is inferred from the static function IBindingModel::identifier(). + * @tparam BindingModel_t Type of the binding model + */ + template + void registerModel(); + + std::unordered_map> _exchangeModels; //!< Map with factory functions + }; + +} // namespace cadet + +#endif // LIBCADET_BINDINGMODELFACTORY_HPP_ diff --git a/src/libcadet/ModelBuilderImpl.cpp b/src/libcadet/ModelBuilderImpl.cpp index ca5a27e54..d7991d0f5 100644 --- a/src/libcadet/ModelBuilderImpl.cpp +++ b/src/libcadet/ModelBuilderImpl.cpp @@ -255,6 +255,11 @@ namespace cadet return _bindingModels.create(name); } + model::IPhaseTransitionModel* ModelBuilder::createExchangeModel(const std::string& name) const + { + return _exchangeModels.create(name); + } + bool ModelBuilder::isValidBindingModel(const std::string& name) const { return _bindingModels.exists(name); diff --git a/src/libcadet/ModelBuilderImpl.hpp b/src/libcadet/ModelBuilderImpl.hpp index 7630b7d50..efdeab30b 100644 --- a/src/libcadet/ModelBuilderImpl.hpp +++ b/src/libcadet/ModelBuilderImpl.hpp @@ -20,6 +20,7 @@ #include "cadet/ModelBuilder.hpp" #include "BindingModelFactory.hpp" +#include "ExchangeModelFactory.hpp" #include "ReactionModelFactory.hpp" #include "ParameterDependenceFactory.hpp" #include "ConfigurationHelper.hpp" @@ -60,6 +61,7 @@ class ModelBuilder : public IModelBuilder, public IConfigHelper virtual IInletProfile* createInletProfile(const std::string& type) const; virtual model::IBindingModel* createBindingModel(const std::string& name) const; + virtual model::IPhaseTransitionModel* createExchangeModel(const std::string& name) const; virtual bool isValidBindingModel(const std::string& name) const; virtual model::IDynamicReactionModel* createDynamicReactionModel(const std::string& name) const; virtual bool isValidDynamicReactionModel(const std::string& name) const; @@ -88,6 +90,7 @@ class ModelBuilder : public IModelBuilder, public IConfigHelper void registerModel(); BindingModelFactory _bindingModels; //!< Factory for IBindingModel implementations + ExchangeModelFactory _exchangeModels; //!< Factory for IPhaseTransitionModel implementations ReactionModelFactory _reactionModels; //!< Factory for IDynamicReactionModel implementations ParameterDependenceFactory _paramDeps; //!< Factory for IParameterStateDependence implementations diff --git a/test/Dummies.hpp b/test/Dummies.hpp index 309f5cf02..f07a57683 100644 --- a/test/Dummies.hpp +++ b/test/Dummies.hpp @@ -26,6 +26,7 @@ namespace virtual cadet::IInletProfile* createInletProfile(const std::string& type) const { return nullptr; } virtual cadet::model::IBindingModel* createBindingModel(const std::string& name) const { return nullptr; } + virtual cadet::model::IPhaseTransitionModel* createExchangeModel(const std::string& name) const { return nullptr; } virtual bool isValidBindingModel(const std::string& name) const { return false; } virtual cadet::IExternalFunction* createExternalFunction(const std::string& type) const { return nullptr; } virtual cadet::model::IDynamicReactionModel* createDynamicReactionModel(const std::string& name) const { return nullptr; } From 163dc3b65b129011d827e3f6114ee582f65c2656 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Fri, 9 Aug 2024 15:30:29 +0200 Subject: [PATCH 10/17] add jacobian to MCT modul --- .../model/MultiChannelTransportModel.cpp | 41 ++--- .../model/MultiChannelTransportModel.hpp | 8 +- src/libcadet/model/PhaseTransitionModel.hpp | 14 +- src/libcadet/model/UnitOperationBase.cpp | 10 ++ src/libcadet/model/UnitOperationBase.hpp | 1 + .../model/exchange/LinearExchange.cpp | 160 ++++++++---------- ...ltiChannelConvectionDispersionOperator.cpp | 71 +++++++- ...ltiChannelConvectionDispersionOperator.hpp | 2 + 8 files changed, 181 insertions(+), 126 deletions(-) diff --git a/src/libcadet/model/MultiChannelTransportModel.cpp b/src/libcadet/model/MultiChannelTransportModel.cpp index d997a0ac4..de0d06dd9 100644 --- a/src/libcadet/model/MultiChannelTransportModel.cpp +++ b/src/libcadet/model/MultiChannelTransportModel.cpp @@ -408,24 +408,23 @@ bool MultiChannelTransportModel::configureModelDiscretization(IParameterProvider if (_dynReactionBulk->usesParamProviderInDiscretizationConfig()) paramProvider.popScope(); } - bool exchangeConfSuccess = true; - //_exchange = nullptr - /*std::vector exchModelNames = {"NONE"}; - if (paramProvider.exists("EXCHANGE_MODEL")) - exchModelNames = paramProvider.getStringArray("EXCHANGE_MODEL"); - */ + + clearExchangeModels(); - bool exchConfSuccess = true; - //exchConfSuccess = _exchange->configureModelDiscretization(paramProvider, _disc.nComp, _disc.nChannel , _disc.boundOffset) && bindingConfSuccess; - //Q: + i * _disc.nBound + i * _disc.nComp for LRMP + _exchange.push_back(nullptr); + + _exchange[0] = helper.createExchangeModel("LINEAR_EX"); + + bool exchangeConfSuccess = true; + exchangeConfSuccess = _exchange[0]->configureModelDiscretization(paramProvider, _disc.nComp, _disc.nChannel, _disc.nCol); const bool transportSuccess = _convDispOp.configureModelDiscretization(paramProvider, helper, _disc.nComp, _disc.nCol, _disc.nChannel, _dynReactionBulk); // Setup the memory for tempState based on state vector _tempState = new double[numDofs()]; - return transportSuccess && reactionConfSuccess && exchConfSuccess; + return transportSuccess && reactionConfSuccess && exchangeConfSuccess; } bool MultiChannelTransportModel::configure(IParameterProvider& paramProvider) @@ -457,18 +456,15 @@ bool MultiChannelTransportModel::configure(IParameterProvider& paramProvider) _exchange[0]->fillChannelInitialParameters(initParams.data(), _unitOpIdx); } - + */ // Reconfigure exchange model bool exchangeConfSuccess = true; - if (_exchange[0] && paramProvider.exists("exchange") && _exchange[0]->requiresConfiguration()){ + + exchangeConfSuccess = _exchange[0]->configure(paramProvider, _unitOpIdx); + - paramProvider.pushScope("exchange"); - exchangeConfSuccess = _exchange[0]->configure(paramProvider, _unitOpIdx, cadet::ParTypeIndep); // Brauche ich PartypeIndep ? - paramProvider.popScope(); - } - */ - return transportSuccess && dynReactionConfSuccess; + return transportSuccess && dynReactionConfSuccess && exchangeConfSuccess; } @@ -736,13 +732,21 @@ int MultiChannelTransportModel::residual(const SimulationTime& simTime, const Co template int MultiChannelTransportModel::residualImpl(double t, unsigned int secIdx, StateType const* const y, double const* const yDot, ResidualType* const res, util::ThreadLocalStorage& threadLocalMem) { + // Handle inlet DOFs, which are simply copied to res for (unsigned int i = 0; i < _disc.nComp * _disc.nChannel; ++i) { res[i] = y[i]; } + _convDispOp.residual(*this, t, secIdx, y, yDot, res, wantJac, typename ParamSens::enabled()); + + linalg::BandedSparseRowIterator jacbegin = _convDispOp.getJacRowIterator(0); + if (_disc.nChannel >= 1) { + _exchange[0]->residual(y, res, typename ParamSens::enabled(), wantJac, jacbegin); + } + if (!_dynReactionBulk || (_dynReactionBulk->numReactionsLiquid() == 0)) return 0; @@ -870,7 +874,6 @@ void MultiChannelTransportModel::multiplyWithJacobian(const SimulationTime& simT */ void MultiChannelTransportModel::multiplyWithDerivativeJacobian(const SimulationTime& simTime, const ConstSimulationState& simState, double const* sDot, double* ret) { - // Add Exchange here _convDispOp.multiplyWithDerivativeJacobian(simTime, sDot, ret); // Handle inlet DOFs (all algebraic) diff --git a/src/libcadet/model/MultiChannelTransportModel.hpp b/src/libcadet/model/MultiChannelTransportModel.hpp index d8ee9c531..0cab9eee7 100644 --- a/src/libcadet/model/MultiChannelTransportModel.hpp +++ b/src/libcadet/model/MultiChannelTransportModel.hpp @@ -204,13 +204,11 @@ class MultiChannelTransportModel : public UnitOperationBase parts::MultiChannelConvectionDispersionOperator _convDispOp; //!< Convection dispersion operator for interstitial volume transport IDynamicReactionModel* _dynReactionBulk; //!< Dynamic reactions in the bulk volume - IPhaseTransitionModel* _phaseTransModel; //!< Phase transition model + + std::vector _exchange; //!< Exchange transition model - linalg::DoubleSparseMatrix _jacInlet; //!< Jacobian inlet DOF block matrix connects inlet DOFs to first bulk cells - linalg::CompressedSparseMatrix _jacC; //!< Jacobian - std::vector _exchangeMatrix; //!< Matrix of exchange coeffs for the linear inter-channel transport - std::vector _crossSections; //!< Cross section area of each compartment + linalg::DoubleSparseMatrix _jacInlet; //!< Jacobian inlet DOF block matrix connects inlet DOFs to first bulk cells bool _analyticJac; //!< Determines whether AD or analytic Jacobians are used diff --git a/src/libcadet/model/PhaseTransitionModel.hpp b/src/libcadet/model/PhaseTransitionModel.hpp index ef849647f..f14da70d1 100644 --- a/src/libcadet/model/PhaseTransitionModel.hpp +++ b/src/libcadet/model/PhaseTransitionModel.hpp @@ -55,9 +55,9 @@ class IPhaseTransitionModel //virtual bool requiresConfiguration() const CADET_NOEXCEPT = 0; - //virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBound, unsigned int const* boundOffset) = 0; + virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int nChannel, unsigned int nCol) = 0; - virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) = 0; + virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx) = 0; // virtual void fillBoundPhaseInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) const CADET_NOEXCEPT = 0; @@ -96,12 +96,12 @@ class IPhaseTransitionModel //virtual unsigned int requiredADdirs() const CADET_NOEXCEPT = 0; - virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, active const* y, active* res, WithParamSensitivity) const = 0; - virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, active const* y, active* res, WithoutParamSensitivity) const = 0; - virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, double const* y, active* res, WithParamSensitivity) const = 0; - virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, double const* y, double* res, WithoutParamSensitivity) const = 0; + virtual int residual(active const* y, active* res, WithParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const = 0; + virtual int residual(active const* y, active* res, WithoutParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const = 0; + virtual int residual(double const* y, active* res, WithParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const = 0; + virtual int residual(double const* y, double* res, WithoutParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const = 0; - virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, linalg::BandedSparseRowIterator jac) const = 0; + //virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, linalg::BandedSparseRowIterator jac) const = 0; #ifdef ENABLE_DG //virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; #endif diff --git a/src/libcadet/model/UnitOperationBase.cpp b/src/libcadet/model/UnitOperationBase.cpp index d54c54307..f18dd3999 100644 --- a/src/libcadet/model/UnitOperationBase.cpp +++ b/src/libcadet/model/UnitOperationBase.cpp @@ -75,6 +75,16 @@ void UnitOperationBase::clearBindingModels() CADET_NOEXCEPT _binding.clear(); } +void UnitOperationBase::clearExchangeModels() CADET_NOEXCEPT +{ + + for (IPhaseTransitionModel* bm : _exchange){ + delete bm; + } + + _exchange.clear(); +} + void UnitOperationBase::clearDynamicReactionModels() CADET_NOEXCEPT { if (_singleDynReaction) diff --git a/src/libcadet/model/UnitOperationBase.hpp b/src/libcadet/model/UnitOperationBase.hpp index 8c4cc2515..d60ae8841 100644 --- a/src/libcadet/model/UnitOperationBase.hpp +++ b/src/libcadet/model/UnitOperationBase.hpp @@ -72,6 +72,7 @@ class UnitOperationBase : public IUnitOperation void clearBindingModels() CADET_NOEXCEPT; void clearDynamicReactionModels() CADET_NOEXCEPT; + void clearExchangeModels() CADET_NOEXCEPT; void configureNonlinearSolver(IParameterProvider& paramProvider); void configureNonlinearSolver(); diff --git a/src/libcadet/model/exchange/LinearExchange.cpp b/src/libcadet/model/exchange/LinearExchange.cpp index 12cdf02bc..cec14d4bd 100644 --- a/src/libcadet/model/exchange/LinearExchange.cpp +++ b/src/libcadet/model/exchange/LinearExchange.cpp @@ -22,6 +22,8 @@ #include "SimulationTypes.hpp" #include "Memory.hpp" #include "model/PhaseTransitionModel.hpp" +#include "ParamReaderHelper.hpp" +#include "model/parts/MultiChannelConvectionDispersionOperator.hpp" #include #include @@ -369,8 +371,8 @@ class LinearExchangeBase : public IPhaseTransitionModel LinearExchangeBase() : _nComp(0), _nChannel(0),_nBound(0), _nCol(0) { } virtual ~LinearExchangeBase() CADET_NOEXCEPT { } - static const char* identifier() { return "a"; } - virtual const char* name() const CADET_NOEXCEPT { return "a"; } + static const char* identifier() { return "LINEAR_EX"; } + virtual const char* name() const CADET_NOEXCEPT { return "LINEAR_EX"; } virtual bool requiresConfiguration() const CADET_NOEXCEPT { return true; } virtual bool usesParamProviderInDiscretizationConfig() const CADET_NOEXCEPT { return true; } @@ -383,14 +385,21 @@ class LinearExchangeBase : public IPhaseTransitionModel return true; } - virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) + virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx) { _parameters.clear(); - // Read parameters (k_a and k_d) - //_paramHandler.configure(paramProvider, _nComp, _nBound); + // Read parameterreadParameterMatrixs (exchange matrix and cross section) + readParameterMatrix(_exchangeMatrix, paramProvider, "EXCHANGE_MATRIX", _nChannel * _nChannel * _nComp, 1); // include parameterPeaderHelp in exchange modul + + _crossSections = paramProvider.getDoubleArray("CHANNEL_CROSS_SECTION_AREAS"); + //ad::copyToAd(initCrossSections.data(), _crossSections.data(), _nChannel); + + //paramProvider.getDoubleArray("EXCHANGE_MATRIX"); + //paramProvider.getDoubleArray("CROSS_SECTION"); //push scrope and pop scope navigation - // Register parameters + // Spaeter + // Register parameters //_paramHandler.registerParameters(_parameters, unitOpIdx, parTypeIdx, _nComp, _nBound); return true; @@ -465,35 +474,45 @@ class LinearExchangeBase : public IPhaseTransitionModel // CADET_BINDINGMODELBASE_BOILERPLATE // which just expands to the six implementations below. - virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, - std::vector crossSections, active const* y, active* res, WithParamSensitivity) const + virtual int residual(active const* y, active* res, WithParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const { - return fluxImpl( nChannel, nComp, nCol, exchangeMatrix, crossSections, y, res); + if (wantJac) + return residualImpl(reinterpret_cast(y), res, jacBegin); + else + return residualImpl(y, res, jacBegin); } - virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, - std::vector crossSections, active const* y, active* res, WithoutParamSensitivity) const - { - return fluxImpl(nChannel, nComp, nCol, exchangeMatrix, crossSections, y, res); + virtual int residual(active const* y, active* res, WithoutParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const + { + if (wantJac) + return residualImpl(reinterpret_cast(y), res, jacBegin); + else + return residualImpl(y, res, jacBegin); } - - virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, - std::vector crossSections, double const* y, active* res, WithParamSensitivity) const + + virtual int residual(double const* y, active* res, WithParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const { - return fluxImpl(nChannel, nComp, nCol, exchangeMatrix, crossSections, y, res); + if(wantJac) + return residualImpl(y, res, jacBegin); + else + return residualImpl(y, res, jacBegin); } - virtual int flux(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, - std::vector crossSections, double const* y, double* res, WithoutParamSensitivity) const + virtual int residual(double const* y, double* res, WithoutParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const { - return fluxImpl(nChannel, nComp, nCol, exchangeMatrix, crossSections, y, res); + if (wantJac) + return residualImpl(y, res, jacBegin); + else + return residualImpl(y, res, jacBegin); } - + + + /* virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, linalg::BandedSparseRowIterator jac) const { return jacobianImpl(nChannel, nComp, nCol, _exchangeMatrix, y, jac); } - /* + virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, active* res, linalg::BandMatrix::RowIterator jac) const { return jacobianImpl(nChannel, nComp, nCol, _exchangeMatrix, y, res, jac); @@ -564,11 +583,14 @@ class LinearExchangeBase : public IPhaseTransitionModel protected: int _nComp; //!< Number of components - unsigned int const* _nBound; //!< Array with number of bound states for each component + unsigned int const* _nBound = nullptr; //!< Array with number of bound states for each component unsigned int _nChannel; //!< Total number of bound states unsigned int _nCol; //!< Number of columns std::vector _reactionQuasistationarity; //!< Determines whether each bound state is quasi-stationary (@c true) or not (@c false) + std::vector _exchangeMatrix; //!< Matrix of exchange coeffs for the linear inter-channel transport + std::vector _crossSections; + parts::MultiChannelConvectionDispersionOperator _conDis; //ParamHandler_t _paramHandler; //!< Parameters std::unordered_map _parameters; //!< Map used to translate ParameterIds to actual variables @@ -576,34 +598,35 @@ class LinearExchangeBase : public IPhaseTransitionModel virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; } - template - int fluxImpl(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, std::vector crossSections, StateType const* y, ResidualType* res) const - { - const unsigned int offsetC = nChannel * nComp; - for (unsigned int col = 0; col < nCol; ++col) + template + int residualImpl(StateType const* y, ResidualType* res, linalg::BandedSparseRowIterator jacBegin) const + { + + const unsigned int offsetC = _nChannel * _nComp; + for (unsigned int col = 0; col < _nCol; ++col) { - const unsigned int offsetColBlock = col * nChannel * nComp; + const unsigned int offsetColBlock = col * _nChannel * _nComp; ResidualType* const resColBlock = res + offsetC + offsetColBlock; StateType const* const yColBlock = y + offsetC + offsetColBlock; - for (unsigned int rad_orig = 0; rad_orig < nChannel; ++rad_orig) + for (unsigned int rad_orig = 0; rad_orig < _nChannel; ++rad_orig) { - const unsigned int offsetToRadOrigBlock = rad_orig * nComp; + const unsigned int offsetToRadOrigBlock = rad_orig * _nComp; const unsigned int offsetColRadOrigBlock = offsetColBlock + offsetToRadOrigBlock; ResidualType* const resColRadOrigBlock = resColBlock + offsetToRadOrigBlock; StateType const* const yColRadOrigBlock = yColBlock + offsetToRadOrigBlock; - for (unsigned int rad_dest = 0; rad_dest < nChannel; ++rad_dest) + for (unsigned int rad_dest = 0; rad_dest < _nChannel; ++rad_dest) { if (rad_orig == rad_dest) continue; - const unsigned int offsetToRadDestBlock = rad_dest * nComp; + const unsigned int offsetToRadDestBlock = rad_dest * _nComp; const unsigned int offsetColRadDestBlock = offsetColBlock + offsetToRadDestBlock; ResidualType* const resColRadDestBlock = resColBlock + offsetToRadDestBlock; // StateType const* const yColRadDestBlock = yColBlock + offsetToRadDestBlock; - for (unsigned int comp = 0; comp < nComp; ++comp) + for (unsigned int comp = 0; comp < _nComp; ++comp) { const unsigned int offsetCur_orig = offsetColRadOrigBlock + comp; const unsigned int offsetCur_dest = offsetColRadDestBlock + comp; @@ -612,86 +635,47 @@ class LinearExchangeBase : public IPhaseTransitionModel ResidualType* const resCur_orig = resColRadOrigBlock + comp; ResidualType* const resCur_dest = resColRadDestBlock + comp; - const ParamType exchange_orig_dest_comp = static_cast(exchangeMatrix[rad_orig * nChannel * nComp + rad_dest * nComp + comp]); + const ParamType exchange_orig_dest_comp = static_cast(_exchangeMatrix[rad_orig * _nChannel * _nComp + rad_dest * _nComp + comp]); if (cadet_likely(exchange_orig_dest_comp > 0.0)) { *resCur_orig += exchange_orig_dest_comp * yCur_orig[0]; - *resCur_dest -= exchange_orig_dest_comp * yCur_orig[0] * static_cast(crossSections[rad_orig]) / static_cast(crossSections[rad_dest]); - } - } + *resCur_dest -= exchange_orig_dest_comp * yCur_orig[0] * static_cast(_crossSections[rad_orig]) / static_cast(_crossSections[rad_dest]); - } + if (wantJac) { - } - } + linalg::BandedSparseRowIterator jacorig; + jacorig = jacBegin + offsetCur_orig; + jacorig[0] += static_cast(exchange_orig_dest_comp); - return 0; - } + linalg::BandedSparseRowIterator jacdest; + jacdest = jacBegin + offsetCur_dest; + jacdest[static_cast(offsetCur_orig) - static_cast(offsetCur_dest)] -= static_cast(exchange_orig_dest_comp); - - //template - inline void jacobianImpl(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector exchangeMatrix, double const* y, linalg::BandedSparseRowIterator jacBegin) const - { - - const unsigned int offsetC = nChannel * nComp; - for (unsigned int col = 0; col < nCol; ++col) - { - const unsigned int offsetColBlock = col * nChannel * nComp; - double const* const yColBlock = y + offsetC + offsetColBlock; + } - for (unsigned int rad_orig = 0; rad_orig < nChannel; ++rad_orig) - { - const unsigned int offsetToRadOrigBlock = rad_orig * nComp; - const unsigned int offsetColRadOrigBlock = offsetColBlock + offsetToRadOrigBlock; - double const* const yColRadOrigBlock = yColBlock + offsetToRadOrigBlock; - for (unsigned int rad_dest = 0; rad_dest < nChannel; ++rad_dest) - { - if (rad_orig == rad_dest) - continue; - - const unsigned int offsetToRadDestBlock = rad_dest * nComp; - const unsigned int offsetColRadDestBlock = offsetColBlock + offsetToRadDestBlock; - // StateType const* const yColRadDestBlock = yColBlock + offsetToRadDestBlock; - - for (unsigned int comp = 0; comp < nComp; ++comp) - { - const unsigned int offsetCur_orig = offsetColRadOrigBlock + comp; - const unsigned int offsetCur_dest = offsetColRadDestBlock + comp; - double const* const yCur_orig = yColRadOrigBlock + comp; - // StateType const* const yCur_dest = yColRadDestBlock + comp; + } - const active exchange_orig_dest_comp = static_cast(exchangeMatrix[rad_orig * nChannel * nComp + rad_dest * nComp + comp]); - if (cadet_likely(exchange_orig_dest_comp > 0.0)) - { - linalg::BandedSparseRowIterator jacorig; - jacorig = jacBegin + offsetCur_orig; - jacorig[0] += static_cast(exchange_orig_dest_comp); - - linalg::BandedSparseRowIterator jacdest; - jacdest = jacBegin + offsetCur_dest; - jacdest[static_cast(offsetCur_orig) - static_cast(offsetCur_dest)] -= static_cast(exchange_orig_dest_comp); - } } } } } - + return 0; } }; -//typedef LinearExchangeBase LinearExchange; // typedef LinearExchangeBase ExternalLinearExchange; +typedef LinearExchangeBase LinearExchange; namespace exchange { - void registerLinearModel(std::unordered_map>& exchange) + void registerLinearExModel(std::unordered_map>& exchange) { - //exchange[LinearExchange::identifier()] = []() { return new LinearExchange(); }; + exchange[LinearExchange::identifier()] = []() { return new LinearExchange(); }; //exchange[ExternalLinearExchange::identifier()] = []() { return new ExternalLinearExchange(); }; } } // namespace exchange diff --git a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp index 0b1355e09..0b31bb6b3 100644 --- a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp +++ b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp @@ -748,8 +748,7 @@ bool MultiChannelConvectionDispersionOperator::configure(UnitOpIdx unitOpIdx, IP _colLength = paramProvider.getDouble("COL_LENGTH"); const std::vector initCrossSections = paramProvider.getDoubleArray("CHANNEL_CROSS_SECTION_AREAS"); ad::copyToAd(initCrossSections.data(), _crossSections.data(), _nChannel); - readParameterMatrix(_exchangeMatrix, paramProvider, "EXCHANGE_MATRIX", _nChannel * _nChannel * _nComp, 1); - + readParameterMatrix(_exchangeMatrix, paramProvider, "EXCHANGE_MATRIX", _nChannel * _nChannel * _nComp, 1); // include parameterPeaderHelp in exchange modul // Read section dependent parameters (transport) // Read VELOCITY @@ -941,6 +940,8 @@ const active& MultiChannelConvectionDispersionOperator::axialDispersion(unsigned * @return @c 0 on success, @c -1 on non-recoverable error, and @c +1 on recoverable error */ + + int MultiChannelConvectionDispersionOperator::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, double* res, bool wantJac, WithoutParamSensitivity) { if (wantJac) @@ -1010,17 +1011,72 @@ int MultiChannelConvectionDispersionOperator::residualImpl(const IModel& model, convdisp::residualKernelAxial(SimulationTime{t, secIdx}, y, yDot, res, _jacC.row(i * _nComp), fp); } + // Handle inter-channel transport - if (cadet_unlikely(_nChannel <= 1)) + if (cadet_unlikely(_nChannel <= 1)) { return 0; + } + /* + linalg::BandedSparseRowIterator jacBegin = _jacC.row(0); - _phaseTransitionModel->flux(_nChannel, _nComp, _nCol, _exchangeMatrix, _crossSections, y, res, typename ParamSens::enabled()); + const unsigned int offsetC = _nChannel * _nComp; + for (unsigned int col = 0; col < _nCol; ++col) + { + const unsigned int offsetColBlock = col * _nChannel * _nComp; + ResidualType* const resColBlock = res + offsetC + offsetColBlock; + StateType const* const yColBlock = y + offsetC + offsetColBlock; - if (wantJac){ - _phaseTransitionModel->analyticJacobian(_nChannel, _nComp, _nCol, _exchangeMatrix, reinterpret_cast(y), _jacC.row(0)); - } + for (unsigned int rad_orig = 0; rad_orig < _nChannel; ++rad_orig) + { + const unsigned int offsetToRadOrigBlock = rad_orig * _nComp; + const unsigned int offsetColRadOrigBlock = offsetColBlock + offsetToRadOrigBlock; + ResidualType* const resColRadOrigBlock = resColBlock + offsetToRadOrigBlock; + StateType const* const yColRadOrigBlock = yColBlock + offsetToRadOrigBlock; + + for (unsigned int rad_dest = 0; rad_dest < _nChannel; ++rad_dest) + { + if (rad_orig == rad_dest) + continue; + const unsigned int offsetToRadDestBlock = rad_dest * _nComp; + const unsigned int offsetColRadDestBlock = offsetColBlock + offsetToRadDestBlock; + ResidualType* const resColRadDestBlock = resColBlock + offsetToRadDestBlock; + // StateType const* const yColRadDestBlock = yColBlock + offsetToRadDestBlock; + + for (unsigned int comp = 0; comp < _nComp; ++comp) + { + const unsigned int offsetCur_orig = offsetColRadOrigBlock + comp; + const unsigned int offsetCur_dest = offsetColRadDestBlock + comp; + StateType const* const yCur_orig = yColRadOrigBlock + comp; + // StateType const* const yCur_dest = yColRadDestBlock + comp; + ResidualType* const resCur_orig = resColRadOrigBlock + comp; + ResidualType* const resCur_dest = resColRadDestBlock + comp; + + const ParamType exchange_orig_dest_comp = static_cast(_exchangeMatrix[rad_orig * _nChannel * _nComp + rad_dest * _nComp + comp]); + if (cadet_likely(exchange_orig_dest_comp > 0.0)) + { + *resCur_orig += exchange_orig_dest_comp * yCur_orig[0]; + *resCur_dest -= exchange_orig_dest_comp * yCur_orig[0] * static_cast(_crossSections[rad_orig]) / static_cast(_crossSections[rad_dest]); + + if (wantJac) { + + linalg::BandedSparseRowIterator jacorig; + jacorig = jacBegin + offsetCur_orig; + jacorig[0] += static_cast(exchange_orig_dest_comp); + + linalg::BandedSparseRowIterator jacdest; + jacdest = jacBegin + offsetCur_dest; + jacdest[static_cast(offsetCur_orig) - static_cast(offsetCur_dest)] -= static_cast(exchange_orig_dest_comp); + } + } + + } + + } + } + } + */ return 0; } @@ -1039,6 +1095,7 @@ void MultiChannelConvectionDispersionOperator::setSparsityPattern() for (unsigned int i = 0; i < _nChannel; ++i) cadet::model::parts::convdisp::sparsityPatternAxial(pattern.row(i * _nComp), _nComp, _nCol, _nComp * _nChannel, static_cast(_curVelocity[i]), _weno); + // auslagern zum exchange if (_nChannel > 1) { for (unsigned int col = 0; col < _nCol; ++col) diff --git a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.hpp b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.hpp index 8fe8a0369..f2268a961 100644 --- a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.hpp +++ b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.hpp @@ -72,6 +72,8 @@ class MultiChannelConvectionDispersionOperator bool configure(UnitOpIdx unitOpIdx, IParameterProvider& paramProvider, std::unordered_map& parameters); bool notifyDiscontinuousSectionTransition(double t, unsigned int secIdx); + linalg::BandedSparseRowIterator getJacRowIterator(const int rowIdx) { return _jacC.row(rowIdx); } + int residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, double* res, bool wantJac, WithoutParamSensitivity); int residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithoutParamSensitivity); int residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity); From 4025d83a210fa2c500dc2e9b1e452e052146db4f Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Thu, 15 Aug 2024 12:33:00 +0200 Subject: [PATCH 11/17] Clean Up --- .../model/MultiChannelTransportModel.cpp | 9 - src/libcadet/model/PhaseTransitionModel.hpp | 49 -- .../model/exchange/LinearExchange.cpp | 432 +----------------- ...ltiChannelConvectionDispersionOperator.cpp | 62 +-- 4 files changed, 8 insertions(+), 544 deletions(-) diff --git a/src/libcadet/model/MultiChannelTransportModel.cpp b/src/libcadet/model/MultiChannelTransportModel.cpp index de0d06dd9..cbb3a6d57 100644 --- a/src/libcadet/model/MultiChannelTransportModel.cpp +++ b/src/libcadet/model/MultiChannelTransportModel.cpp @@ -449,17 +449,8 @@ bool MultiChannelTransportModel::configure(IParameterProvider& paramProvider) paramProvider.popScope(); } - /* Zunkufscode - if (_exchange[0]){ - - std::vector initParams(_disc.nChannel); //Q was macht initParams ? brauche ich maxBoundStates ? - _exchange[0]->fillChannelInitialParameters(initParams.data(), _unitOpIdx); - } - - */ // Reconfigure exchange model bool exchangeConfSuccess = true; - exchangeConfSuccess = _exchange[0]->configure(paramProvider, _unitOpIdx); diff --git a/src/libcadet/model/PhaseTransitionModel.hpp b/src/libcadet/model/PhaseTransitionModel.hpp index f14da70d1..865c83e9c 100644 --- a/src/libcadet/model/PhaseTransitionModel.hpp +++ b/src/libcadet/model/PhaseTransitionModel.hpp @@ -59,60 +59,11 @@ class IPhaseTransitionModel virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx) = 0; - // virtual void fillBoundPhaseInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) const CADET_NOEXCEPT = 0; - - //virtual void fillChannelInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx) const CADET_NOEXCEPT = 0; - - - //virtual void setExternalFunctions(IExternalFunction** extFuns, unsigned int size) = 0; //X - - //virtual std::unordered_map getAllParameterValues() const = 0; - - //static const char* identifier() CADET_NOEXCEPT { return "IPhaseTransitionModel"; } - - //virtual bool hasParameter(const ParameterId& pId) const = 0; - //virtual bool setParameter(const ParameterId& pId, int value) = 0; - //virtual bool setParameter(const ParameterId& pId, double value) = 0; - //virtual bool setParameter(const ParameterId& pId, bool value) = 0; - - //virtual active* getParameter(const ParameterId& pId) = 0; - - //virtual bool hasSalt() const CADET_NOEXCEPT = 0; - - //virtual bool supportsMultistate() const CADET_NOEXCEPT = 0;l bool supportsNonExchange() const CADET_NOEXCEPT = 0; - - - //virtual bool supportsNonBinding() const CADET_NOEXCEPT = 0; - - //virtual bool hasQuasiStationaryReactions() const CADET_NOEXCEPT = 0; - - //virtual bool hasDynamicReactions() const CADET_NOEXCEPT = 0; - - //virtual bool dependsOnTime() const CADET_NOEXCEPT = 0; - - //virtual bool requiresWorkspace() const CADET_NOEXCEPT = 0; - - //virtual unsigned int workspaceSize(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nBoundStates) const CADET_NOEXCEPT = 0; - - //virtual unsigned int requiredADdirs() const CADET_NOEXCEPT = 0; - virtual int residual(active const* y, active* res, WithParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const = 0; virtual int residual(active const* y, active* res, WithoutParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const = 0; virtual int residual(double const* y, active* res, WithParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const = 0; virtual int residual(double const* y, double* res, WithoutParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const = 0; - //virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, linalg::BandedSparseRowIterator jac) const = 0; -#ifdef ENABLE_DG - //virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; -#endif - virtual void timeDerivativeQuasiStationaryFluxes(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yCp, double const* y, double* dResDt, LinearBufferAllocator workSpace) const = 0; - - virtual int const* reactionQuasiStationarity() const CADET_NOEXCEPT = 0; - - virtual bool preConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const = 0; - - virtual void postConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const = 0; - protected: }; diff --git a/src/libcadet/model/exchange/LinearExchange.cpp b/src/libcadet/model/exchange/LinearExchange.cpp index cec14d4bd..733b907fd 100644 --- a/src/libcadet/model/exchange/LinearExchange.cpp +++ b/src/libcadet/model/exchange/LinearExchange.cpp @@ -37,338 +37,17 @@ namespace cadet namespace model { - -/** - * @brief Handles linear exchange model parameters that do not depend on external functions - */ -class LinearParamHandler : public ConstParamHandlerBase -{ -public: - - /** - * @brief Holds actual parameter data - */ - struct ConstParams - { - std::vector kA; //!< Adsorption rate - std::vector kD; //!< Desorption rate - }; - - typedef ConstParams params_t; - typedef ConstParams const* ParamsHandle; - - /** - * @brief Returns name of the binding model - * @return Name of the binding model - */ - static const char* identifier() { return "LINEAR_EX"; } - - LinearParamHandler() CADET_NOEXCEPT : _kA(&_localParams.kA), _kD(&_localParams.kD) { } - - /** - * @brief Reads parameters and verifies them - * @details See IexchangeModel::configure() for details. - * @param [in] paramProvider IParameterProvider used for reading parameters - * @param [in] nComp Number of components - * @param [in] nChannel Array with number of bound states for each component - * @return @c true if the parameters were read and validated successfully, otherwise @c false - */ - inline bool configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nChannel) - { - _kA.configure("LIN_KA", paramProvider, nComp, nChannel); - _kD.configure("LIN_KD", paramProvider, nComp, nChannel); - return validateConfig(nComp, nChannel); - } - - /** - * @brief Registers all local parameters in a map for further use - * @param [in,out] parameters Map in which the parameters are stored - * @param [in] unitOpIdx Index of the unit operation used for registering the parameters - * @param [in] nComp Number of components - * @param [in] nChannel Array with number of bound states for each component - */ - inline void registerParameters(std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nChannel) - { - _kA.registerParam("LIN_KA", parameters, unitOpIdx, parTypeIdx, nComp, nChannel); - _kD.registerParam("LIN_KD", parameters, unitOpIdx, parTypeIdx, nComp, nChannel); - } - - /** - * @brief Reserves space in the storage of the parameters - * @param [in] numElem Number of elements (total) - * @param [in] numSlices Number of slices / exchange site types - * @param [in] nComp Number of components - * @param [in] nChannel Array with number of bound states for each component - */ - inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nChannel) \ - { - _kA.reserve(numElem, numSlices, nComp, nChannel); - _kD.reserve(numElem, numSlices, nComp, nChannel); - } - - /** - * @brief Updates the parameter cache in order to take the external profile into account - * @param [in] t Current time - * @param [in] z Axial coordinate in the column - * @param [in] r Radial coordinate in the bead - * @param [in] secIdx Index of the current section - * @param [in] nComp Number of components - * @param [in] nChannel Array with number of bound states for each component - * @param [in,out] workSpace Memory buffer for updated data - * @return Externally dependent parameter values - */ - inline ParamsHandle update(double t, unsigned int secIdx, const ColumnPosition& colPos, unsigned int nComp, unsigned int const* nChannel, LinearBufferAllocator& workSpace) const - { - return &_localParams; - } - - /** - * @brief Calculates time derivative in case of external dependence - * @param [in] t Current time - * @param [in] z Axial coordinate in the column - * @param [in] r Radial coordinate in the bead - * @param [in] secIdx Index of the current section - * @param [in] nComp Number of components - * @param [in] nChannel Array with number of bound states for each component - * @param [in,out] workSpace Memory buffer for updated data - * @return Time derivatives of externally dependent parameters - */ - inline std::tuple updateTimeDerivative(double t, unsigned int secIdx, const ColumnPosition& colPos, unsigned int nComp, unsigned int const* nChannel, LinearBufferAllocator& workSpace) const - { - return std::make_tuple(&_localParams, nullptr); - } - -protected: - - /** - * @brief Validates recently read parameters - * @param [in] nComp Number of components - * @param [in] nChannel Array with number of bound states for each component - * @return @c true if the parameters were validated successfully, otherwise @c false - */ - inline bool validateConfig(unsigned int nComp, unsigned int const* nChannel) - { - if ((_kA.size() != _kD.size()) || (_kA.size() < nComp)) - throw InvalidParameterException("LIN_KA and LIN_KD have to have the same size"); - - return true; - } - - ConstParams _localParams; //!< Actual parameter data - - // Handlers provide configure(), reserve(), and registerParam() for parameters - ScalarComponentDependentParameter _kA; //!< Handler for adsorption rate - ScalarComponentDependentParameter _kD; //!< Handler for desorption rate -}; - -/** - * @brief Handles linear exchange model parameters that depend on an external function - */ -class ExtLinearParamHandler : public ExternalParamHandlerBase -{ -public: - - /** - * @brief Holds actual parameter data - * @details The parameter data will be stored in a memory buffer. This requires - * control over the location of the data, which is not provided by - * std::vector and util::SlicedVector. There are "local" equivalents, - * util::LocalVector and util::LocalSlicedVector, for these types that - * allow to control the placement of the payload data. A single value (i.e., - * a single active or double) can simply be put in the struct. - */ - struct VariableParams - { - util::LocalVector kA; //!< Adsorption rate - util::LocalVector kD; //!< Desorption rate - }; - - typedef VariableParams params_t; - typedef ConstBufferedScalar ParamsHandle; - - ExtLinearParamHandler() CADET_NOEXCEPT { } - - /** - * @brief Returns name of the exchange model - * @return Name of the exchange model - */ - static const char* identifier() { return "EXT_LINEAR"; } - - /** - * @brief Reads parameters and verifies them - * @details See IBindingModel::configure() for details. - * @param [in] paramProvider IParameterProvider used for reading parameters - * @param [in] nComp Number of components - * @param [in] nChannel Array with number of bound states for each component - * @return @c true if the parameters were read and validated successfully, otherwise @c false - */ - inline bool configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nChannel) - { - _kA.configure("LIN_KA", paramProvider, nComp, nChannel); - _kD.configure("LIN_KD", paramProvider, nComp, nChannel); - - // Number of externally dependent parameters (2) needs to be given to ExternalParamHandlerBase::configure() - ExternalParamHandlerBase::configure(paramProvider, 2); - return validateConfig(nComp, nChannel); - } - - /** - * @brief Registers all local parameters in a map for further use - * @param [in,out] parameters Map in which the parameters are stored - * @param [in] unitOpIdx Index of the unit operation used for registering the parameters - * @param [in] nComp Number of components - * @param [in] nChannel Array with number of bound states for each component - */ - inline void registerParameters(std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nChannel) - { - _kA.registerParam("LIN_KA", parameters, unitOpIdx, parTypeIdx, nComp, nChannel); - _kD.registerParam("LIN_KD", parameters, unitOpIdx, parTypeIdx, nComp, nChannel); - } - - /** - * @brief Reserves space in the storage of the parameters - * @param [in] numElem Number of elements (total) - * @param [in] numSlices Number of slices / exchange site types - * @param [in] nComp Number of components - * @param [in] nChannel Array with number of bound states for each component - */ - inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nChannel) \ - { - _kA.reserve(numElem, numSlices, nComp, nChannel); - _kD.reserve(numElem, numSlices, nComp, nChannel); - } - - /** - * @brief Updates the parameter cache in order to take the external profile into account - * @details The data and the returned value are constructed in the given @p workSpace memory buffer. - * @param [in] t Current time - * @param [in] z Axial coordinate in the column - * @param [in] r Radial coordinate in the bead - * @param [in] secIdx Index of the current section - * @param [in] nComp Number of components - * @param [in] nChannel Array with number of bound states for each component - * @param [in,out] workSpace Memory buffer for updated data - * @return Externally dependent parameter values - */ - inline ParamsHandle update(double t, unsigned int secIdx, const ColumnPosition& colPos, unsigned int nComp, unsigned int const* nChannel, LinearBufferAllocator& workSpace) const - { - // Allocate params_t and buffer for function evaluation - BufferedScalar localParams = workSpace.scalar(); - BufferedArray extFunBuffer = workSpace.array(2); - - // Evaluate external functions in buffer - evaluateExternalFunctions(t, secIdx, colPos, 2, static_cast(extFunBuffer)); - - // Prepare the buffer for the data and update the data - _kA.prepareCache(localParams->kA, workSpace); - _kA.update(cadet::util::dataOfLocalVersion(localParams->kA), extFunBuffer[0], nComp, nChannel); - - _kD.prepareCache(localParams->kD, workSpace); - _kD.update(cadet::util::dataOfLocalVersion(localParams->kD), extFunBuffer[1], nComp, nChannel); - - return localParams; - } - - /** - * @brief Calculates time derivative in case of external dependence - * @details The time derivatives are constructed in the given @p workSpace memory buffer. - * @param [in] t Current time - * @param [in] z Axial coordinate in the column - * @param [in] r Radial coordinate in the bead - * @param [in] secIdx Index of the current section - * @param [in] nComp Number of components - * @param [in] nChannel Array with number of bound states for each component - * @param [in,out] workSpace Memory buffer for updated data - * @return Tuple with externally dependent parameter values and their time derivatives - */ - inline std::tuple updateTimeDerivative(double t, unsigned int secIdx, const ColumnPosition& colPos, unsigned int nComp, unsigned int const* nChannel, LinearBufferAllocator& workSpace) const - { - // Allocate params_t for parameters and their time derivatives - BufferedScalar localParams = workSpace.scalar(); - BufferedScalar p = workSpace.scalar(); - - // Allocate buffer for external function values and their time derivatives - BufferedArray extFunBuffer = workSpace.array(2); - BufferedArray extDerivBuffer = workSpace.array(2); - - // Evaluate external functions and their time derivatives - evaluateExternalFunctions(t, secIdx, colPos, 2, static_cast(extFunBuffer)); - evaluateTimeDerivativeExternalFunctions(t, secIdx, colPos, 2, static_cast(extDerivBuffer)); - - // Prepare the buffer for the data and update the data - _kA.prepareCache(localParams->kA, workSpace); - _kA.update(cadet::util::dataOfLocalVersion(localParams->kA), extFunBuffer[0], nComp, nChannel); - - _kA.prepareCache(p->kA, workSpace); - _kA.updateTimeDerivative(cadet::util::dataOfLocalVersion(p->kA), extFunBuffer[0], extDerivBuffer[0], nComp, nChannel); - - _kD.prepareCache(localParams->kD, workSpace); - _kD.update(cadet::util::dataOfLocalVersion(localParams->kD), extFunBuffer[1], nComp, nChannel); - - _kD.prepareCache(p->kD, workSpace); - _kD.updateTimeDerivative(cadet::util::dataOfLocalVersion(p->kD), extFunBuffer[1], extDerivBuffer[1], nComp, nChannel); - - return std::make_tuple(std::move(localParams), std::move(p));; - } - - /** - * @brief Returns how much memory is required for caching in bytes - * @details Memory size in bytes. - * @param [in] nComp Number of components - * @param [in] totalNumBoundStates Total number of bound states - * @param [in] nChannel Array with bound states for each component - * @return Memory size in bytes - */ - inline std::size_t cacheSize(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nChannel) const CADET_NOEXCEPT - { - // Required buffer memory: - // + params_t object - // + buffer for external function evaluations (2 parameters) - // + buffer for external function time derivative evaluations (2 parameters) - // + buffer for actual parameter data (memory for _kA data + memory for _kD data) - // + buffer for parameter time derivatives (memory for _kA data + memory for _kD data) - return 2 * sizeof(params_t) + alignof(params_t) - + 2 * 2 * sizeof(double) + alignof(double) - + 2 * (_kA.additionalDynamicMemory(nComp, totalNumBoundStates, nChannel) + _kD.additionalDynamicMemory(nComp, totalNumBoundStates, nChannel)); - } - -protected: - - /** - * @brief Validates recently read parameters - * @param [in] nComp Number of components - * @param [in] nChannel Array with number of bound states for each component - * @return @c true if the parameters were validated successfully, otherwise @c false - */ - inline bool validateConfig(unsigned int nComp, unsigned int const* nChannel) - { - if ((_kA.size() != _kD.size()) || (_kA.size() < nComp)) - throw InvalidParameterException("LIN_KA and LIN_KD have to have the same size"); - - return true; - } - - // Handlers provide configure(), reserve(), and registerParam() for parameters - ExternalScalarComponentDependentParameter _kA; //!< Handler for adsorption rate - ExternalScalarComponentDependentParameter _kD; //!< Handler for desorption rate -}; - - /** * @brief Defines the linear exchange model - * @details Implements the linear adsorption model \f$ \frac{\mathrm{d}q}{\mathrm{d}t} = k_a c_p - k_d q \f$. - * Multiple bound states are not supported. Components without bound state (i.e., non-exchange components) - * are supported. - * - * See @cite Guiochon2006 - * @tparam ParamHandler_t Type that can add support for external function dependence + * @details Implements the linear exchange model for the MCT model + * The exchange is given by a matrix of exchange coefficients eij for each component i and j and the chross Ai section from channel i. + * The exchange from channel i to all other channel j is given by \f$ \frac{\mathrm{d}ci}{\mathrm{d}t} = \sum_j eij cj Aj/Ai - eji ci \f$. */ -// template class LinearExchangeBase : public IPhaseTransitionModel { public: - LinearExchangeBase() : _nComp(0), _nChannel(0),_nBound(0), _nCol(0) { } + LinearExchangeBase() : _nComp(0), _nChannel(0), _nCol(0) { } virtual ~LinearExchangeBase() CADET_NOEXCEPT { } static const char* identifier() { return "LINEAR_EX"; } @@ -388,19 +67,8 @@ class LinearExchangeBase : public IPhaseTransitionModel virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx) { _parameters.clear(); - - // Read parameterreadParameterMatrixs (exchange matrix and cross section) readParameterMatrix(_exchangeMatrix, paramProvider, "EXCHANGE_MATRIX", _nChannel * _nChannel * _nComp, 1); // include parameterPeaderHelp in exchange modul - _crossSections = paramProvider.getDoubleArray("CHANNEL_CROSS_SECTION_AREAS"); - //ad::copyToAd(initCrossSections.data(), _crossSections.data(), _nChannel); - - //paramProvider.getDoubleArray("EXCHANGE_MATRIX"); - //paramProvider.getDoubleArray("CROSS_SECTION"); //push scrope and pop scope navigation - - // Spaeter - // Register parameters - //_paramHandler.registerParameters(_parameters, unitOpIdx, parTypeIdx, _nComp, _nBound); return true; } @@ -469,11 +137,6 @@ class LinearExchangeBase : public IPhaseTransitionModel virtual void configure(IExternalFunction** extFuns, unsigned int size) {} - // The next three flux() function implementations and two analyticJacobian() function - // implementations are usually hidden behind - // CADET_BINDINGMODELBASE_BOILERPLATE - // which just expands to the six implementations below. - virtual int residual(active const* y, active* res, WithParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const { if (wantJac) @@ -507,91 +170,15 @@ class LinearExchangeBase : public IPhaseTransitionModel } - /* - virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, linalg::BandedSparseRowIterator jac) const - { - return jacobianImpl(nChannel, nComp, nCol, _exchangeMatrix, y, jac); - } - - virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, active* res, linalg::BandMatrix::RowIterator jac) const - { - return jacobianImpl(nChannel, nComp, nCol, _exchangeMatrix, y, res, jac); - } - - virtual void analyticJacobian(unsigned int nChannel, unsigned int nComp, unsigned int nCol, std::vector _exchangeMatrix, double const* y, double* res, linalg::BandMatrix::RowIterator jac) const - { - return jacobianImpl(nChannel, nComp, nCol, _exchangeMatrix, y, res, jac); - } - */ - virtual void timeDerivativeQuasiStationaryFluxes(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yCp, double const* y, double* dResDt, LinearBufferAllocator workSpace) const - { - if (!hasQuasiStationaryReactions()) - return; - - //if (!ParamHandler_t::dependsOnTime()) - // return; - - // Update external function and compute time derivative of parameters - //typename ParamHandler_t::ParamsHandle p; - //typename ParamHandler_t::ParamsHandle dpDt; - /* - std::tie(p, dpDt) = _paramHandler.updateTimeDerivative(t, secIdx, colPos, _nComp, _nBound, workSpace); - - unsigned int bndIdx = 0; - for (int i = 0; i < _nComp; ++i) - { - // Skip components without bound states (bound state index bndIdx is not advanced) - if (_nBound[i] == 0) - continue; - - dResDt[bndIdx] = -static_cast(dpDt->kA[i]) * yCp[i] + static_cast(dpDt->kD[i]) * y[bndIdx]; - - // Next bound component - ++bndIdx; - }*/ - } - - virtual bool preConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const - { - // Due to mass conservation, we need to run a nonlinear solver (although the model is simple). - return true; - } - - virtual void postConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const - { - // There's nothing to do here since the algebraic equations have already been solved in preConsistentInitialState() - } - - virtual bool hasSalt() const CADET_NOEXCEPT { return false; } - virtual bool supportsMultistate() const CADET_NOEXCEPT { return false; } - virtual bool supportsNonExchange() const CADET_NOEXCEPT { return true; } - - virtual bool hasQuasiStationaryReactions() const CADET_NOEXCEPT - { - return std::any_of(_reactionQuasistationarity.begin(), _reactionQuasistationarity.end(), [](int i) -> bool { return i; }); - } - - virtual bool hasDynamicReactions() const CADET_NOEXCEPT - { - return std::any_of(_reactionQuasistationarity.begin(), _reactionQuasistationarity.end(), [](int i) -> bool { return !static_cast(i); }); - } - - virtual bool dependsOnTime() const CADET_NOEXCEPT { return false; } - virtual bool requiresWorkspace() const CADET_NOEXCEPT { return false; } - virtual unsigned int requiredADdirs() const CADET_NOEXCEPT { return 0; } - virtual int const* reactionQuasiStationarity() const CADET_NOEXCEPT { return _reactionQuasistationarity.data(); } - protected: int _nComp; //!< Number of components - unsigned int const* _nBound = nullptr; //!< Array with number of bound states for each component unsigned int _nChannel; //!< Total number of bound states unsigned int _nCol; //!< Number of columns std::vector _reactionQuasistationarity; //!< Determines whether each bound state is quasi-stationary (@c true) or not (@c false) std::vector _exchangeMatrix; //!< Matrix of exchange coeffs for the linear inter-channel transport - std::vector _crossSections; - parts::MultiChannelConvectionDispersionOperator _conDis; - //ParamHandler_t _paramHandler; //!< Parameters + std::vector _crossSections; //!< Cross sections of the channels + parts::MultiChannelConvectionDispersionOperator _conDis; //!< Convection dispersion operator std::unordered_map _parameters; //!< Map used to translate ParameterIds to actual variables @@ -624,14 +211,12 @@ class LinearExchangeBase : public IPhaseTransitionModel const unsigned int offsetToRadDestBlock = rad_dest * _nComp; const unsigned int offsetColRadDestBlock = offsetColBlock + offsetToRadDestBlock; ResidualType* const resColRadDestBlock = resColBlock + offsetToRadDestBlock; - // StateType const* const yColRadDestBlock = yColBlock + offsetToRadDestBlock; for (unsigned int comp = 0; comp < _nComp; ++comp) { const unsigned int offsetCur_orig = offsetColRadOrigBlock + comp; const unsigned int offsetCur_dest = offsetColRadDestBlock + comp; StateType const* const yCur_orig = yColRadOrigBlock + comp; - // StateType const* const yCur_dest = yColRadDestBlock + comp; ResidualType* const resCur_orig = resColRadOrigBlock + comp; ResidualType* const resCur_dest = resColRadDestBlock + comp; @@ -667,16 +252,13 @@ class LinearExchangeBase : public IPhaseTransitionModel } }; - -// typedef LinearExchangeBase ExternalLinearExchange; - typedef LinearExchangeBase LinearExchange; + namespace exchange { void registerLinearExModel(std::unordered_map>& exchange) { exchange[LinearExchange::identifier()] = []() { return new LinearExchange(); }; - //exchange[ExternalLinearExchange::identifier()] = []() { return new ExternalLinearExchange(); }; } } // namespace exchange diff --git a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp index 0b31bb6b3..2e66eb345 100644 --- a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp +++ b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp @@ -1013,70 +1013,10 @@ int MultiChannelConvectionDispersionOperator::residualImpl(const IModel& model, - // Handle inter-channel transport if (cadet_unlikely(_nChannel <= 1)) { return 0; } - /* - linalg::BandedSparseRowIterator jacBegin = _jacC.row(0); - - const unsigned int offsetC = _nChannel * _nComp; - for (unsigned int col = 0; col < _nCol; ++col) - { - const unsigned int offsetColBlock = col * _nChannel * _nComp; - ResidualType* const resColBlock = res + offsetC + offsetColBlock; - StateType const* const yColBlock = y + offsetC + offsetColBlock; - - for (unsigned int rad_orig = 0; rad_orig < _nChannel; ++rad_orig) - { - const unsigned int offsetToRadOrigBlock = rad_orig * _nComp; - const unsigned int offsetColRadOrigBlock = offsetColBlock + offsetToRadOrigBlock; - ResidualType* const resColRadOrigBlock = resColBlock + offsetToRadOrigBlock; - StateType const* const yColRadOrigBlock = yColBlock + offsetToRadOrigBlock; - - for (unsigned int rad_dest = 0; rad_dest < _nChannel; ++rad_dest) - { - if (rad_orig == rad_dest) - continue; - - const unsigned int offsetToRadDestBlock = rad_dest * _nComp; - const unsigned int offsetColRadDestBlock = offsetColBlock + offsetToRadDestBlock; - ResidualType* const resColRadDestBlock = resColBlock + offsetToRadDestBlock; - // StateType const* const yColRadDestBlock = yColBlock + offsetToRadDestBlock; - - for (unsigned int comp = 0; comp < _nComp; ++comp) - { - const unsigned int offsetCur_orig = offsetColRadOrigBlock + comp; - const unsigned int offsetCur_dest = offsetColRadDestBlock + comp; - StateType const* const yCur_orig = yColRadOrigBlock + comp; - // StateType const* const yCur_dest = yColRadDestBlock + comp; - ResidualType* const resCur_orig = resColRadOrigBlock + comp; - ResidualType* const resCur_dest = resColRadDestBlock + comp; - - const ParamType exchange_orig_dest_comp = static_cast(_exchangeMatrix[rad_orig * _nChannel * _nComp + rad_dest * _nComp + comp]); - if (cadet_likely(exchange_orig_dest_comp > 0.0)) - { - *resCur_orig += exchange_orig_dest_comp * yCur_orig[0]; - *resCur_dest -= exchange_orig_dest_comp * yCur_orig[0] * static_cast(_crossSections[rad_orig]) / static_cast(_crossSections[rad_dest]); - - if (wantJac) { - - linalg::BandedSparseRowIterator jacorig; - jacorig = jacBegin + offsetCur_orig; - jacorig[0] += static_cast(exchange_orig_dest_comp); - - linalg::BandedSparseRowIterator jacdest; - jacdest = jacBegin + offsetCur_dest; - jacdest[static_cast(offsetCur_orig) - static_cast(offsetCur_dest)] -= static_cast(exchange_orig_dest_comp); - } - } - - } - - } - } - } - */ + return 0; } From c9715a53bb7293dbcf262f8ec6a6ba9778c2ca1b Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Tue, 3 Sep 2024 15:43:22 +0200 Subject: [PATCH 12/17] clean up changes according to review process --- src/libcadet/CMakeLists.txt | 2 +- src/libcadet/ConfigurationHelper.hpp | 4 ++-- src/libcadet/ExchangeModelFactory.cpp | 7 +++---- src/libcadet/ExchangeModelFactory.hpp | 8 ++++---- src/libcadet/ModelBuilderImpl.cpp | 2 +- src/libcadet/ModelBuilderImpl.hpp | 4 ++-- ...seTransitionModel.hpp => ExchangeModel.hpp} | 18 +++++++++--------- src/libcadet/model/ModelUtils.hpp | 6 +++--- .../model/MultiChannelTransportModel.cpp | 10 ++++++---- .../model/MultiChannelTransportModel.hpp | 2 +- src/libcadet/model/UnitOperationBase.cpp | 2 +- src/libcadet/model/UnitOperationBase.hpp | 4 ++-- src/libcadet/model/exchange/LinearExchange.cpp | 11 +++-------- ...ultiChannelConvectionDispersionOperator.cpp | 3 +-- ...ultiChannelConvectionDispersionOperator.hpp | 7 +++---- test/Dummies.hpp | 2 +- 16 files changed, 43 insertions(+), 49 deletions(-) rename src/libcadet/model/{PhaseTransitionModel.hpp => ExchangeModel.hpp} (89%) diff --git a/src/libcadet/CMakeLists.txt b/src/libcadet/CMakeLists.txt index 31c0653d2..b898398b5 100644 --- a/src/libcadet/CMakeLists.txt +++ b/src/libcadet/CMakeLists.txt @@ -101,7 +101,7 @@ set(LIBCADET_BINDINGMODEL_SOURCES ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/HICWaterOnHydrophobicSurfacesBinding.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/HICConstantWaterActivityBinding.cpp ) -set(LIBCADET_PHASETRANSITIONMODELINTERFACE_HPP_ +set(LIBCADET_EXCHANGEMODEL_SOURCES ${CMAKE_SOURCE_DIR}/src/libcadet/model/exchange/LinearExchange.cpp ) diff --git a/src/libcadet/ConfigurationHelper.hpp b/src/libcadet/ConfigurationHelper.hpp index b26b3c297..98bde769c 100644 --- a/src/libcadet/ConfigurationHelper.hpp +++ b/src/libcadet/ConfigurationHelper.hpp @@ -30,7 +30,7 @@ class IExternalFunction; namespace model { class IBindingModel; - class IPhaseTransitionModel; + class IExchangeModel; class IDynamicReactionModel; class IParameterStateDependence; class IParameterParameterDependence; @@ -59,7 +59,7 @@ class IConfigHelper */ virtual model::IBindingModel* createBindingModel(const std::string& name) const = 0; - virtual model::IPhaseTransitionModel* createExchangeModel(const std::string& name) const = 0; + virtual model::IExchangeModel* createExchangeModel(const std::string& name) const = 0; /** diff --git a/src/libcadet/ExchangeModelFactory.cpp b/src/libcadet/ExchangeModelFactory.cpp index f7b521ca0..8a38222b7 100644 --- a/src/libcadet/ExchangeModelFactory.cpp +++ b/src/libcadet/ExchangeModelFactory.cpp @@ -13,7 +13,6 @@ #include "ExchangeModelFactory.hpp" #include "cadet/Exceptions.hpp" -#include "model/binding/SimplifiedMultiStateStericMassActionBinding.hpp" namespace cadet { @@ -21,7 +20,7 @@ namespace cadet { namespace exchange { - void registerLinearExModel(std::unordered_map>& exchange); + void registerLinearExModel(std::unordered_map>& exchange); } } @@ -46,7 +45,7 @@ namespace cadet registerModel(ExchangeModel_t::identifier()); } - model::IPhaseTransitionModel* ExchangeModelFactory::create(const std::string& name) const + model::IExchangeModel* ExchangeModelFactory::create(const std::string& name) const { const auto it = _exchangeModels.find(name); if (it == _exchangeModels.end()) @@ -60,7 +59,7 @@ namespace cadet return it->second(); } - void ExchangeModelFactory::registerModel(const std::string& name, std::function factory) + void ExchangeModelFactory::registerModel(const std::string& name, std::function factory) { if (_exchangeModels.find(name) == _exchangeModels.end()) _exchangeModels[name] = factory; diff --git a/src/libcadet/ExchangeModelFactory.hpp b/src/libcadet/ExchangeModelFactory.hpp index 16b2568a3..d4b14731f 100644 --- a/src/libcadet/ExchangeModelFactory.hpp +++ b/src/libcadet/ExchangeModelFactory.hpp @@ -27,7 +27,7 @@ namespace cadet namespace model { - class IPhaseTransitionModel; + class IExchangeModel; } /** @@ -49,14 +49,14 @@ namespace cadet * @param [in] name Name of the binding model * @return The binding model or @c NULL if a binding model with this name does not exist */ - model::IPhaseTransitionModel* create(const std::string& name) const; + model::IExchangeModel* create(const std::string& name) const; /** * @brief Registers the given binding model implementation * @param [in] name Name of the IBindingModel implementation * @param [in] factory Function that creates an object of the IBindingModel class */ - void registerModel(const std::string& name, std::function factory); + void registerModel(const std::string& name, std::function factory); /** * @brief Returns whether a binding model of the given name @p name exists @@ -82,7 +82,7 @@ namespace cadet template void registerModel(); - std::unordered_map> _exchangeModels; //!< Map with factory functions + std::unordered_map> _exchangeModels; //!< Map with factory functions }; } // namespace cadet diff --git a/src/libcadet/ModelBuilderImpl.cpp b/src/libcadet/ModelBuilderImpl.cpp index d7991d0f5..6eb76911a 100644 --- a/src/libcadet/ModelBuilderImpl.cpp +++ b/src/libcadet/ModelBuilderImpl.cpp @@ -255,7 +255,7 @@ namespace cadet return _bindingModels.create(name); } - model::IPhaseTransitionModel* ModelBuilder::createExchangeModel(const std::string& name) const + model::IExchangeModel* ModelBuilder::createExchangeModel(const std::string& name) const { return _exchangeModels.create(name); } diff --git a/src/libcadet/ModelBuilderImpl.hpp b/src/libcadet/ModelBuilderImpl.hpp index efdeab30b..1136e757d 100644 --- a/src/libcadet/ModelBuilderImpl.hpp +++ b/src/libcadet/ModelBuilderImpl.hpp @@ -61,7 +61,7 @@ class ModelBuilder : public IModelBuilder, public IConfigHelper virtual IInletProfile* createInletProfile(const std::string& type) const; virtual model::IBindingModel* createBindingModel(const std::string& name) const; - virtual model::IPhaseTransitionModel* createExchangeModel(const std::string& name) const; + virtual model::IExchangeModel* createExchangeModel(const std::string& name) const; virtual bool isValidBindingModel(const std::string& name) const; virtual model::IDynamicReactionModel* createDynamicReactionModel(const std::string& name) const; virtual bool isValidDynamicReactionModel(const std::string& name) const; @@ -90,7 +90,7 @@ class ModelBuilder : public IModelBuilder, public IConfigHelper void registerModel(); BindingModelFactory _bindingModels; //!< Factory for IBindingModel implementations - ExchangeModelFactory _exchangeModels; //!< Factory for IPhaseTransitionModel implementations + ExchangeModelFactory _exchangeModels; //!< Factory for IExchangeModel implementations ReactionModelFactory _reactionModels; //!< Factory for IDynamicReactionModel implementations ParameterDependenceFactory _paramDeps; //!< Factory for IParameterStateDependence implementations diff --git a/src/libcadet/model/PhaseTransitionModel.hpp b/src/libcadet/model/ExchangeModel.hpp similarity index 89% rename from src/libcadet/model/PhaseTransitionModel.hpp rename to src/libcadet/model/ExchangeModel.hpp index 865c83e9c..70ef4a03e 100644 --- a/src/libcadet/model/PhaseTransitionModel.hpp +++ b/src/libcadet/model/ExchangeModel.hpp @@ -1,9 +1,9 @@ // ============================================================================= // CADET -// +// // Copyright © 2008-2024: The CADET Authors // Please see the AUTHORS and CONTRIBUTORS file. -// +// // All rights reserved. This program and the accompanying materials // are made available under the terms of the GNU Public License v3.0 (or, at // your option, any later version) which accompanies this distribution, and @@ -11,12 +11,12 @@ // ============================================================================= /** - * @file - * Defines the BindingModel interface. + * @file + * Defines the ExchangeModel interface. */ -#ifndef LIBCADET_PHASETRANSITIONMODELINTERFACE_HPP_ -#define LIBCADET_PHASETRANSITIONMODELINTERFACE_HPP_ +#ifndef LIBCADET_EXCHANGEMODELINTERFACE_HPP_ +#define LIBCADET_EXCHANGEMODELINTERFACE_HPP_ #include @@ -45,11 +45,11 @@ struct ColumnPosition; namespace model { -class IPhaseTransitionModel -{ +class IExchangeModel +{ public: - virtual ~IPhaseTransitionModel() CADET_NOEXCEPT { } + virtual ~IExchangeModel() CADET_NOEXCEPT { } virtual const char* name() const CADET_NOEXCEPT = 0; diff --git a/src/libcadet/model/ModelUtils.hpp b/src/libcadet/model/ModelUtils.hpp index fdad98b21..703d3b17f 100644 --- a/src/libcadet/model/ModelUtils.hpp +++ b/src/libcadet/model/ModelUtils.hpp @@ -83,10 +83,10 @@ inline unsigned int numBoundStates(unsigned int const* const nBound, unsigned in } /** - * @brief Returns the number of bound states - * @param [in] nBound Array with number of bound states for each component + * @brief Returns the number of channel states for each component + * @param [in] nChannel Number of channels * @param [in] nComp Number of components - * @return Number of bound states + * @return Number of channel states of each component */ inline unsigned int numExchangeStates(unsigned int nChannel, unsigned int nComp) { diff --git a/src/libcadet/model/MultiChannelTransportModel.cpp b/src/libcadet/model/MultiChannelTransportModel.cpp index cbb3a6d57..952f6d004 100644 --- a/src/libcadet/model/MultiChannelTransportModel.cpp +++ b/src/libcadet/model/MultiChannelTransportModel.cpp @@ -18,7 +18,7 @@ #include "cadet/SolutionRecorder.hpp" #include "ConfigurationHelper.hpp" #include "model/ReactionModel.hpp" -#include "model/PhaseTransitionModel.hpp" +#include "model/ExchangeModel.hpp" #include "SimulationTypes.hpp" #include "Stencil.hpp" @@ -414,7 +414,11 @@ bool MultiChannelTransportModel::configureModelDiscretization(IParameterProvider _exchange.push_back(nullptr); - _exchange[0] = helper.createExchangeModel("LINEAR_EX"); + if (paramProvider.exists("EXCHANGE_MODEL")) + _exchange[0] = helper.createExchangeModel(paramProvider.getString("EXCHANGE_MODEL")); + + if (!_exchange[0]) + _exchange[0] = helper.createExchangeModel("LINEAR_EX"); bool exchangeConfSuccess = true; exchangeConfSuccess = _exchange[0]->configureModelDiscretization(paramProvider, _disc.nComp, _disc.nChannel, _disc.nCol); @@ -463,8 +467,6 @@ unsigned int MultiChannelTransportModel::threadLocalMemorySize() const CADET_NOE { LinearMemorySizer lms; - //Add exchange memory - // Memory for residualImpl() if (_dynReactionBulk && _dynReactionBulk->requiresWorkspace()) lms.fitBlock(_dynReactionBulk->workspaceSize(_disc.nComp, 0, nullptr)); diff --git a/src/libcadet/model/MultiChannelTransportModel.hpp b/src/libcadet/model/MultiChannelTransportModel.hpp index 0cab9eee7..67f779241 100644 --- a/src/libcadet/model/MultiChannelTransportModel.hpp +++ b/src/libcadet/model/MultiChannelTransportModel.hpp @@ -205,7 +205,7 @@ class MultiChannelTransportModel : public UnitOperationBase parts::MultiChannelConvectionDispersionOperator _convDispOp; //!< Convection dispersion operator for interstitial volume transport IDynamicReactionModel* _dynReactionBulk; //!< Dynamic reactions in the bulk volume - std::vector _exchange; //!< Exchange transition model + std::vector _exchange; //!< Exchange transition model linalg::DoubleSparseMatrix _jacInlet; //!< Jacobian inlet DOF block matrix connects inlet DOFs to first bulk cells diff --git a/src/libcadet/model/UnitOperationBase.cpp b/src/libcadet/model/UnitOperationBase.cpp index f18dd3999..dde7cbb42 100644 --- a/src/libcadet/model/UnitOperationBase.cpp +++ b/src/libcadet/model/UnitOperationBase.cpp @@ -78,7 +78,7 @@ void UnitOperationBase::clearBindingModels() CADET_NOEXCEPT void UnitOperationBase::clearExchangeModels() CADET_NOEXCEPT { - for (IPhaseTransitionModel* bm : _exchange){ + for (IExchangeModel* bm : _exchange){ delete bm; } diff --git a/src/libcadet/model/UnitOperationBase.hpp b/src/libcadet/model/UnitOperationBase.hpp index d60ae8841..4e99915a7 100644 --- a/src/libcadet/model/UnitOperationBase.hpp +++ b/src/libcadet/model/UnitOperationBase.hpp @@ -34,7 +34,7 @@ namespace model { class IBindingModel; -class IPhaseTransitionModel; +class IExchangeModel; class IDynamicReactionModel; /** @@ -80,7 +80,7 @@ class UnitOperationBase : public IUnitOperation UnitOpIdx _unitOpIdx; //!< Unit operation index std::vector _binding; //!< Binding model - std::vector _exchange; //!< Exchange model + std::vector _exchange; //!< Exchange model bool _singleBinding; //!< Determines whether only a single binding model is present std::vector _dynReaction; //!< Dynamic reaction model bool _singleDynReaction; //!< Determines whether only a single dynamic reaction model is present diff --git a/src/libcadet/model/exchange/LinearExchange.cpp b/src/libcadet/model/exchange/LinearExchange.cpp index 733b907fd..d36df09f9 100644 --- a/src/libcadet/model/exchange/LinearExchange.cpp +++ b/src/libcadet/model/exchange/LinearExchange.cpp @@ -21,7 +21,7 @@ #include "LocalVector.hpp" #include "SimulationTypes.hpp" #include "Memory.hpp" -#include "model/PhaseTransitionModel.hpp" +#include "model/ExchangeModel.hpp" #include "ParamReaderHelper.hpp" #include "model/parts/MultiChannelConvectionDispersionOperator.hpp" @@ -43,7 +43,7 @@ namespace model * The exchange is given by a matrix of exchange coefficients eij for each component i and j and the chross Ai section from channel i. * The exchange from channel i to all other channel j is given by \f$ \frac{\mathrm{d}ci}{\mathrm{d}t} = \sum_j eij cj Aj/Ai - eji ci \f$. */ -class LinearExchangeBase : public IPhaseTransitionModel +class LinearExchangeBase : public IExchangeModel { public: @@ -129,11 +129,6 @@ class LinearExchangeBase : public IPhaseTransitionModel return nullptr; } - virtual unsigned int workspaceSize(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nChannel) const CADET_NOEXCEPT - { - // return _paramHandler.cacheSize(nComp, totalNumBoundStates, nChannel); - return 0; - } virtual void configure(IExternalFunction** extFuns, unsigned int size) {} @@ -256,7 +251,7 @@ typedef LinearExchangeBase LinearExchange; namespace exchange { - void registerLinearExModel(std::unordered_map>& exchange) + void registerLinearExModel(std::unordered_map>& exchange) { exchange[LinearExchange::identifier()] = []() { return new LinearExchange(); }; } diff --git a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp index 2e66eb345..ac43ff42d 100644 --- a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp +++ b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp @@ -831,7 +831,7 @@ bool MultiChannelConvectionDispersionOperator::configure(UnitOpIdx unitOpIdx, IP setSparsityPattern(); - _phaseTransitionModel = new LinearExchangeBase(); + _exchangeModel = new LinearExchangeBase(); return true; } @@ -1035,7 +1035,6 @@ void MultiChannelConvectionDispersionOperator::setSparsityPattern() for (unsigned int i = 0; i < _nChannel; ++i) cadet::model::parts::convdisp::sparsityPatternAxial(pattern.row(i * _nComp), _nComp, _nCol, _nComp * _nChannel, static_cast(_curVelocity[i]), _weno); - // auslagern zum exchange if (_nChannel > 1) { for (unsigned int col = 0; col < _nCol; ++col) diff --git a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.hpp b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.hpp index f2268a961..b417486b4 100644 --- a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.hpp +++ b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.hpp @@ -26,7 +26,7 @@ #include "model/ParameterMultiplexing.hpp" #include "SimulationTypes.hpp" #include "ConfigurationHelper.hpp" -#include "model/PhaseTransitionModel.hpp" +#include "model/ExchangeModel.hpp" #include #include @@ -149,7 +149,6 @@ class MultiChannelConvectionDispersionOperator unsigned int _nChannel; //!< Number of channels bool _hasDynamicReactions; //!< Determines whether the model has dynamic reactions (only relevant for sparsity pattern) - active _colLength; //!< Column length \f$ L \f$ std::vector _crossSections; //!< Cross section area of each compartment @@ -162,9 +161,9 @@ class MultiChannelConvectionDispersionOperator std::vector _dir; //!< Current flow direction bool _singleVelocity; //!< Determines whether only one velocity for all compartments is given - std::vector _exchangeMatrix; //!< Matrix of exchange coeffs for the linear inter-channel transport + std::vector _exchangeMatrix; //!< Matrix of exchange coeffs for the inter-channel transport - IPhaseTransitionModel* _phaseTransitionModel; //!< Phase transition model + IExchangeModel* _exchangeModel; //!< Phase transition model IParameterParameterDependence* _dispersionDep; diff --git a/test/Dummies.hpp b/test/Dummies.hpp index f07a57683..03a558160 100644 --- a/test/Dummies.hpp +++ b/test/Dummies.hpp @@ -26,7 +26,7 @@ namespace virtual cadet::IInletProfile* createInletProfile(const std::string& type) const { return nullptr; } virtual cadet::model::IBindingModel* createBindingModel(const std::string& name) const { return nullptr; } - virtual cadet::model::IPhaseTransitionModel* createExchangeModel(const std::string& name) const { return nullptr; } + virtual cadet::model::IExchangeModel* createExchangeModel(const std::string& name) const { return nullptr; } virtual bool isValidBindingModel(const std::string& name) const { return false; } virtual cadet::IExternalFunction* createExternalFunction(const std::string& type) const { return nullptr; } virtual cadet::model::IDynamicReactionModel* createDynamicReactionModel(const std::string& name) const { return nullptr; } From b65c242a4ad631c321567cac16c857e8f0caac93 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Tue, 1 Oct 2024 18:29:56 +0200 Subject: [PATCH 13/17] Little changes --- src/libcadet/ExchangeModelFactory.cpp | 2 +- src/libcadet/ExchangeModelFactory.hpp | 38 +++++++++---------- .../model/MultiChannelTransportModel.cpp | 2 +- src/libcadet/model/binding/LinearBinding.cpp | 9 +++-- .../model/exchange/LinearExchange.cpp | 4 +- 5 files changed, 28 insertions(+), 27 deletions(-) diff --git a/src/libcadet/ExchangeModelFactory.cpp b/src/libcadet/ExchangeModelFactory.cpp index 8a38222b7..a127e5065 100644 --- a/src/libcadet/ExchangeModelFactory.cpp +++ b/src/libcadet/ExchangeModelFactory.cpp @@ -50,7 +50,7 @@ namespace cadet const auto it = _exchangeModels.find(name); if (it == _exchangeModels.end()) { - // BindingModel was not found + // ExchangeModel was not found return nullptr; } diff --git a/src/libcadet/ExchangeModelFactory.hpp b/src/libcadet/ExchangeModelFactory.hpp index d4b14731f..e01deb854 100644 --- a/src/libcadet/ExchangeModelFactory.hpp +++ b/src/libcadet/ExchangeModelFactory.hpp @@ -31,53 +31,53 @@ namespace cadet } /** - * @brief Creates binding models + * @brief Creates Exchange models */ class ExchangeModelFactory { public: /** - * @brief Construct the BindingModelFactory - * @details All internal binding models are registered here. + * @brief Construct the ExchangeModelFactory + * @details All internal exchange models are registered here. */ ExchangeModelFactory(); ~ExchangeModelFactory(); /** - * @brief Creates binding models with the given @p name - * @param [in] name Name of the binding model - * @return The binding model or @c NULL if a binding model with this name does not exist + * @brief Creates exchange models with the given @p name + * @param [in] name Name of the exchange model + * @return The exchange model or @c NULL if a exchange model with this name does not exist */ model::IExchangeModel* create(const std::string& name) const; /** - * @brief Registers the given binding model implementation - * @param [in] name Name of the IBindingModel implementation - * @param [in] factory Function that creates an object of the IBindingModel class + * @brief Registers the given exchange model implementation + * @param [in] name Name of the IExchangeModel implementation + * @param [in] factory Function that creates an object of the IExchangeModel class */ void registerModel(const std::string& name, std::function factory); /** - * @brief Returns whether a binding model of the given name @p name exists - * @param [in] name Name of the binding model - * @return @c true if a binding model of this name exists, otherwise @c false + * @brief Returns whether a exchange model of the given name @p name exists + * @param [in] name Name of the exchange model + * @return @c true if a exchange model of this name exists, otherwise @c false */ bool exists(const std::string& name) const; protected: /** - * @brief Registers an IBindingModel - * @param [in] name Name of the binding model - * @tparam BindingModel_t Type of the binding model + * @brief Registers an IExchangeModel + * @param [in] name Name of the exchange model + * @tparam ExchangeModel_t Type of the exchange model */ template void registerModel(const std::string& name); /** - * @brief Registers an IBindingModel - * @details The name of the binding model is inferred from the static function IBindingModel::identifier(). - * @tparam BindingModel_t Type of the binding model + * @brief Registers an IExchangeModel + * @details The name of the exchange model is inferred from the static function IExchangeModel::identifier(). + * @tparam ExchangeModel_t Type of the exchange model */ template void registerModel(); @@ -87,4 +87,4 @@ namespace cadet } // namespace cadet -#endif // LIBCADET_BINDINGMODELFACTORY_HPP_ +#endif // LIBCADET_EXCHANGEMODELFACTORY_HPP_ diff --git a/src/libcadet/model/MultiChannelTransportModel.cpp b/src/libcadet/model/MultiChannelTransportModel.cpp index 952f6d004..14c48aea9 100644 --- a/src/libcadet/model/MultiChannelTransportModel.cpp +++ b/src/libcadet/model/MultiChannelTransportModel.cpp @@ -418,7 +418,7 @@ bool MultiChannelTransportModel::configureModelDiscretization(IParameterProvider _exchange[0] = helper.createExchangeModel(paramProvider.getString("EXCHANGE_MODEL")); if (!_exchange[0]) - _exchange[0] = helper.createExchangeModel("LINEAR_EX"); + _exchange[0] = helper.createExchangeModel("LINEAR"); bool exchangeConfSuccess = true; exchangeConfSuccess = _exchange[0]->configureModelDiscretization(paramProvider, _disc.nComp, _disc.nChannel, _disc.nCol); diff --git a/src/libcadet/model/binding/LinearBinding.cpp b/src/libcadet/model/binding/LinearBinding.cpp index 048cf0b2b..0c7a1a2db 100644 --- a/src/libcadet/model/binding/LinearBinding.cpp +++ b/src/libcadet/model/binding/LinearBinding.cpp @@ -125,7 +125,8 @@ class LinearParamHandler : public ConstParamHandlerBase * @param [in] nComp Number of components * @param [in] nBoundStates Array with number of bound states for each component */ - inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nBoundStates) { + inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nBoundStates) + { _kA.reserve(numElem, numSlices, nComp, nBoundStates); _kD.reserve(numElem, numSlices, nComp, nBoundStates); } @@ -387,7 +388,7 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase * See @cite Guiochon2006 * @tparam ParamHandler_t Type that can add support for external function dependence */ -template // nicht fuer LinearExchange +template class LinearBindingBase : public IBindingModel { public: @@ -525,7 +526,7 @@ class LinearBindingBase : public IBindingModel { return fluxImpl(t, secIdx, colPos, y, yCp, res, workSpace); } - + virtual int flux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, active const* yCp, active* res, LinearBufferAllocator workSpace, WithoutParamSensitivity) const { @@ -543,7 +544,7 @@ class LinearBindingBase : public IBindingModel { return fluxImpl(t, secIdx, colPos, y, yCp, res, workSpace); } - + virtual void analyticJacobian(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int offsetCp, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const { jacobianImpl(t, secIdx, colPos, y, offsetCp, jac, workSpace); diff --git a/src/libcadet/model/exchange/LinearExchange.cpp b/src/libcadet/model/exchange/LinearExchange.cpp index d36df09f9..373e2b806 100644 --- a/src/libcadet/model/exchange/LinearExchange.cpp +++ b/src/libcadet/model/exchange/LinearExchange.cpp @@ -50,8 +50,8 @@ class LinearExchangeBase : public IExchangeModel LinearExchangeBase() : _nComp(0), _nChannel(0), _nCol(0) { } virtual ~LinearExchangeBase() CADET_NOEXCEPT { } - static const char* identifier() { return "LINEAR_EX"; } - virtual const char* name() const CADET_NOEXCEPT { return "LINEAR_EX"; } + static const char* identifier() { return "LINEAR"; } + virtual const char* name() const CADET_NOEXCEPT { return "LINEAR"; } virtual bool requiresConfiguration() const CADET_NOEXCEPT { return true; } virtual bool usesParamProviderInDiscretizationConfig() const CADET_NOEXCEPT { return true; } From 6f200b0fb219b0c2677090c098a5e8262d971003 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Tue, 1 Oct 2024 18:40:59 +0200 Subject: [PATCH 14/17] More little changes --- src/libcadet/model/MultiChannelTransportModel.cpp | 3 +-- src/libcadet/model/MultiChannelTransportModel.hpp | 3 --- src/libcadet/model/exchange/LinearExchange.cpp | 10 ++++------ .../parts/MultiChannelConvectionDispersionOperator.cpp | 1 - .../parts/MultiChannelConvectionDispersionOperator.hpp | 1 - 5 files changed, 5 insertions(+), 13 deletions(-) diff --git a/src/libcadet/model/MultiChannelTransportModel.cpp b/src/libcadet/model/MultiChannelTransportModel.cpp index 14c48aea9..f53171cbe 100644 --- a/src/libcadet/model/MultiChannelTransportModel.cpp +++ b/src/libcadet/model/MultiChannelTransportModel.cpp @@ -599,7 +599,6 @@ int MultiChannelTransportModel::residual(const SimulationTime& simTime, const Co // Evaluate residual do not compute Jacobian or parameter sensitivities return residualImpl(simTime.t, simTime.secIdx, simState.vecStateY, simState.vecStateYdot, res, threadLocalMem); } -//Q Do we need different types of residual calculation ? int MultiChannelTransportModel::jacobian(const SimulationTime& simTime, const ConstSimulationState& simState, double* const res, const AdJacobianParams& adJac, util::ThreadLocalStorage& threadLocalMem) { BENCH_SCOPE(_timerResidual); @@ -875,7 +874,7 @@ void MultiChannelTransportModel::multiplyWithDerivativeJacobian(const Simulation void MultiChannelTransportModel::setExternalFunctions(IExternalFunction** extFuns, unsigned int size) { - // Add exchange here + // TODO: Add exchange here } unsigned int MultiChannelTransportModel::localOutletComponentIndex(unsigned int port) const CADET_NOEXCEPT diff --git a/src/libcadet/model/MultiChannelTransportModel.hpp b/src/libcadet/model/MultiChannelTransportModel.hpp index 67f779241..522aa3f3b 100644 --- a/src/libcadet/model/MultiChannelTransportModel.hpp +++ b/src/libcadet/model/MultiChannelTransportModel.hpp @@ -204,16 +204,13 @@ class MultiChannelTransportModel : public UnitOperationBase parts::MultiChannelConvectionDispersionOperator _convDispOp; //!< Convection dispersion operator for interstitial volume transport IDynamicReactionModel* _dynReactionBulk; //!< Dynamic reactions in the bulk volume - std::vector _exchange; //!< Exchange transition model - linalg::DoubleSparseMatrix _jacInlet; //!< Jacobian inlet DOF block matrix connects inlet DOFs to first bulk cells bool _analyticJac; //!< Determines whether AD or analytic Jacobians are used unsigned int _jacobianAdDirs; //!< Number of AD seed vectors required for Jacobian computation - bool _factorizeJacobian; //!< Determines whether the Jacobian needs to be factorized double* _tempState; //!< Temporary storage with the size of the state vector or larger if binding models require it diff --git a/src/libcadet/model/exchange/LinearExchange.cpp b/src/libcadet/model/exchange/LinearExchange.cpp index 373e2b806..494cbb424 100644 --- a/src/libcadet/model/exchange/LinearExchange.cpp +++ b/src/libcadet/model/exchange/LinearExchange.cpp @@ -12,7 +12,6 @@ -#include "model/BindingModel.hpp" #include "model/ExternalFunctionSupport.hpp" #include "ParamIdUtil.hpp" #include "model/ModelUtils.hpp" @@ -40,8 +39,8 @@ namespace model /** * @brief Defines the linear exchange model * @details Implements the linear exchange model for the MCT model - * The exchange is given by a matrix of exchange coefficients eij for each component i and j and the chross Ai section from channel i. - * The exchange from channel i to all other channel j is given by \f$ \frac{\mathrm{d}ci}{\mathrm{d}t} = \sum_j eij cj Aj/Ai - eji ci \f$. + * The exchange is given by a matrix of exchange coefficients e_{ij} for each component i and j and the chross A_{i} section from channel i. + * The exchange from channel i to all other channel j is given by \f$ \frac{\mathrm{d}ci}{\mathrm{d}t} = \sum_j e_{ij} c_j A_j/A_i - e_{ji} c_i \f$. */ class LinearExchangeBase : public IExchangeModel { @@ -58,7 +57,7 @@ class LinearExchangeBase : public IExchangeModel virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int nChannel, unsigned int nCol) { _nComp = nComp; - _nChannel = nChannel; // nChannel realy int ? -> by not nBoundStates not ? + _nChannel = nChannel; _nCol = nCol; return true; @@ -67,7 +66,7 @@ class LinearExchangeBase : public IExchangeModel virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx) { _parameters.clear(); - readParameterMatrix(_exchangeMatrix, paramProvider, "EXCHANGE_MATRIX", _nChannel * _nChannel * _nComp, 1); // include parameterPeaderHelp in exchange modul + readParameterMatrix(_exchangeMatrix, paramProvider, "EXCHANGE_MATRIX", _nChannel * _nChannel * _nComp, 1); // TODO include parameterPeaderHelp in exchange modul _crossSections = paramProvider.getDoubleArray("CHANNEL_CROSS_SECTION_AREAS"); return true; @@ -173,7 +172,6 @@ class LinearExchangeBase : public IExchangeModel std::vector _exchangeMatrix; //!< Matrix of exchange coeffs for the linear inter-channel transport std::vector _crossSections; //!< Cross sections of the channels - parts::MultiChannelConvectionDispersionOperator _conDis; //!< Convection dispersion operator std::unordered_map _parameters; //!< Map used to translate ParameterIds to actual variables diff --git a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp index ac43ff42d..896968be6 100644 --- a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp +++ b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp @@ -21,7 +21,6 @@ #include "linalg/CompressedSparseMatrix.hpp" #include "model/parts/AxialConvectionDispersionKernel.hpp" #include "model/ParameterDependence.hpp" -#include "model/exchange/LinearExchange.cpp" #ifdef SUPERLU_FOUND diff --git a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.hpp b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.hpp index b417486b4..a3bff35bc 100644 --- a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.hpp +++ b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.hpp @@ -162,7 +162,6 @@ class MultiChannelConvectionDispersionOperator bool _singleVelocity; //!< Determines whether only one velocity for all compartments is given std::vector _exchangeMatrix; //!< Matrix of exchange coeffs for the inter-channel transport - IExchangeModel* _exchangeModel; //!< Phase transition model IParameterParameterDependence* _dispersionDep; From e21b192618e2a6809d038a8b00d237880eac4f37 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Fri, 18 Oct 2024 10:19:56 +0200 Subject: [PATCH 15/17] more and more little changes --- src/libcadet/ExchangeModelFactory.cpp | 2 +- src/libcadet/model/MultiChannelTransportModel.cpp | 2 -- .../model/parts/MultiChannelConvectionDispersionOperator.cpp | 2 -- 3 files changed, 1 insertion(+), 5 deletions(-) diff --git a/src/libcadet/ExchangeModelFactory.cpp b/src/libcadet/ExchangeModelFactory.cpp index a127e5065..0f03f5c2d 100644 --- a/src/libcadet/ExchangeModelFactory.cpp +++ b/src/libcadet/ExchangeModelFactory.cpp @@ -64,7 +64,7 @@ namespace cadet if (_exchangeModels.find(name) == _exchangeModels.end()) _exchangeModels[name] = factory; else - throw InvalidParameterException("IPhaseTransition implementation with the name " + name + " is already registered and cannot be overwritten"); + throw InvalidParameterException("IExchange implementation with the name " + name + " is already registered and cannot be overwritten"); } bool ExchangeModelFactory::exists(const std::string& name) const diff --git a/src/libcadet/model/MultiChannelTransportModel.cpp b/src/libcadet/model/MultiChannelTransportModel.cpp index f53171cbe..2bc45f7e5 100644 --- a/src/libcadet/model/MultiChannelTransportModel.cpp +++ b/src/libcadet/model/MultiChannelTransportModel.cpp @@ -724,14 +724,12 @@ int MultiChannelTransportModel::residual(const SimulationTime& simTime, const Co template int MultiChannelTransportModel::residualImpl(double t, unsigned int secIdx, StateType const* const y, double const* const yDot, ResidualType* const res, util::ThreadLocalStorage& threadLocalMem) { - // Handle inlet DOFs, which are simply copied to res for (unsigned int i = 0; i < _disc.nComp * _disc.nChannel; ++i) { res[i] = y[i]; } - _convDispOp.residual(*this, t, secIdx, y, yDot, res, wantJac, typename ParamSens::enabled()); linalg::BandedSparseRowIterator jacbegin = _convDispOp.getJacRowIterator(0); diff --git a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp index 896968be6..357a67f67 100644 --- a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp +++ b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp @@ -830,8 +830,6 @@ bool MultiChannelConvectionDispersionOperator::configure(UnitOpIdx unitOpIdx, IP setSparsityPattern(); - _exchangeModel = new LinearExchangeBase(); - return true; } From 320ba88515c527214d8d078236f4c432bce0d117 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Fri, 18 Oct 2024 11:05:27 +0200 Subject: [PATCH 16/17] fix dependency --- src/libcadet/model/ExchangeModel.hpp | 2 +- src/libcadet/model/exchange/LinearExchange.cpp | 4 +--- .../model/parts/MultiChannelConvectionDispersionOperator.cpp | 1 + 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/src/libcadet/model/ExchangeModel.hpp b/src/libcadet/model/ExchangeModel.hpp index 70ef4a03e..8c9a62b53 100644 --- a/src/libcadet/model/ExchangeModel.hpp +++ b/src/libcadet/model/ExchangeModel.hpp @@ -70,4 +70,4 @@ class IExchangeModel } // namespace model } // namespace cadet -#endif // LIBCADET_PHASETRANSITIONMODELINTERFACE_HPP_ +#endif // LIBCADET_EXCHANGEMODELINTERFACE_HPP_ diff --git a/src/libcadet/model/exchange/LinearExchange.cpp b/src/libcadet/model/exchange/LinearExchange.cpp index 494cbb424..60f64c20f 100644 --- a/src/libcadet/model/exchange/LinearExchange.cpp +++ b/src/libcadet/model/exchange/LinearExchange.cpp @@ -10,8 +10,7 @@ // is available at http://www.gnu.org/licenses/gpl.html // ============================================================================= - - +#include "model/ExchangeModel.hpp" #include "model/ExternalFunctionSupport.hpp" #include "ParamIdUtil.hpp" #include "model/ModelUtils.hpp" @@ -20,7 +19,6 @@ #include "LocalVector.hpp" #include "SimulationTypes.hpp" #include "Memory.hpp" -#include "model/ExchangeModel.hpp" #include "ParamReaderHelper.hpp" #include "model/parts/MultiChannelConvectionDispersionOperator.hpp" diff --git a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp index 357a67f67..2ea3a78d0 100644 --- a/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp +++ b/src/libcadet/model/parts/MultiChannelConvectionDispersionOperator.cpp @@ -21,6 +21,7 @@ #include "linalg/CompressedSparseMatrix.hpp" #include "model/parts/AxialConvectionDispersionKernel.hpp" #include "model/ParameterDependence.hpp" +#include "model/exchange/LinearExchange.cpp" #ifdef SUPERLU_FOUND From f89eb6f133817ad4eeec129ba25578e2dda71cdf Mon Sep 17 00:00:00 2001 From: Jan Breuer Date: Wed, 20 Nov 2024 13:32:01 +0100 Subject: [PATCH 17/17] Update copyright statement --- include/cadet/cadet.hpp | 2 +- src/libcadet/ExchangeModelFactory.cpp | 4 ++-- src/libcadet/ExchangeModelFactory.hpp | 4 ++-- src/libcadet/HighResKoren.hpp | 4 ++-- src/libcadet/model/ExchangeModel.hpp | 4 ++-- src/libcadet/model/exchange/LinearExchange.cpp | 4 ++-- src/libcadet/model/reaction/CrystallizationReaction.cpp | 4 ++-- test/Crystallization.cpp | 4 ++-- 8 files changed, 15 insertions(+), 15 deletions(-) diff --git a/include/cadet/cadet.hpp b/include/cadet/cadet.hpp index 52ad85946..b70864c3d 100644 --- a/include/cadet/cadet.hpp +++ b/include/cadet/cadet.hpp @@ -16,7 +16,7 @@ * * @authors Please refer to CONTRIBUTING.md * @version 5.0.1 - * @date 2008-2024 + * @date 2008-present * @copyright GNU General Public License v3.0 (or, at your option, any later version). */ diff --git a/src/libcadet/ExchangeModelFactory.cpp b/src/libcadet/ExchangeModelFactory.cpp index 0f03f5c2d..fa480e3ab 100644 --- a/src/libcadet/ExchangeModelFactory.cpp +++ b/src/libcadet/ExchangeModelFactory.cpp @@ -1,8 +1,8 @@ // ============================================================================= // CADET // -// Copyright © 2008-2022: The CADET Authors -// Please see the AUTHORS and CONTRIBUTORS file. +// Copyright © 2008-present: The CADET-Core Authors +// Please see the AUTHORS.md file. // // All rights reserved. This program and the accompanying materials // are made available under the terms of the GNU Public License v3.0 (or, at diff --git a/src/libcadet/ExchangeModelFactory.hpp b/src/libcadet/ExchangeModelFactory.hpp index e01deb854..794384221 100644 --- a/src/libcadet/ExchangeModelFactory.hpp +++ b/src/libcadet/ExchangeModelFactory.hpp @@ -1,8 +1,8 @@ // ============================================================================= // CADET // -// Copyright © 2008-2024: The CADET Authors -// Please see the AUTHORS and CONTRIBUTORS file. +// Copyright © 2008-present: The CADET-Core Authors +// Please see the AUTHORS.md file. // // All rights reserved. This program and the accompanying materials // are made available under the terms of the GNU Public License v3.0 (or, at diff --git a/src/libcadet/HighResKoren.hpp b/src/libcadet/HighResKoren.hpp index 861f28597..4999a17e1 100644 --- a/src/libcadet/HighResKoren.hpp +++ b/src/libcadet/HighResKoren.hpp @@ -1,8 +1,8 @@ // ============================================================================= // CADET // -// Copyright © 2008-2022: 2008-present: The CADET-Core Authors -// Please see the AUTHORS and CONTRIBUTORS file. +// Copyright © 2008-present: 2008-present: The CADET-Core Authors +// Please see the AUTHORS.md file. // // All rights reserved. This program and the accompanying materials // are made available under the terms of the GNU Public License v3.0 (or, at diff --git a/src/libcadet/model/ExchangeModel.hpp b/src/libcadet/model/ExchangeModel.hpp index 8c9a62b53..61a786aee 100644 --- a/src/libcadet/model/ExchangeModel.hpp +++ b/src/libcadet/model/ExchangeModel.hpp @@ -1,8 +1,8 @@ // ============================================================================= // CADET // -// Copyright © 2008-2024: The CADET Authors -// Please see the AUTHORS and CONTRIBUTORS file. +// Copyright © 2008-present: The CADET-Core Authors +// Please see the AUTHORS.md file. // // All rights reserved. This program and the accompanying materials // are made available under the terms of the GNU Public License v3.0 (or, at diff --git a/src/libcadet/model/exchange/LinearExchange.cpp b/src/libcadet/model/exchange/LinearExchange.cpp index 60f64c20f..c65b13660 100644 --- a/src/libcadet/model/exchange/LinearExchange.cpp +++ b/src/libcadet/model/exchange/LinearExchange.cpp @@ -1,8 +1,8 @@ // ============================================================================= // CADET // -// Copyright © 2008-2024: The CADET Authors -// Please see the AUTHORS and CONTRIBUTORS file. +// Copyright © 2008-present: The CADET-Core Authors +// Please see the AUTHORS.md file. // // All rights reserved. This program and the accompanying materials // are made available under the terms of the GNU Public License v3.0 (or, at diff --git a/src/libcadet/model/reaction/CrystallizationReaction.cpp b/src/libcadet/model/reaction/CrystallizationReaction.cpp index d9df21641..77eea953b 100644 --- a/src/libcadet/model/reaction/CrystallizationReaction.cpp +++ b/src/libcadet/model/reaction/CrystallizationReaction.cpp @@ -1,8 +1,8 @@ // ============================================================================= // CADET // -// Copyright © 2008-2022: 2008-present: The CADET-Core Authors -// Please see the AUTHORS and CONTRIBUTORS file. +// Copyright © 2008-present: 2008-present: The CADET-Core Authors +// Please see the AUTHORS.md file. // // All rights reserved. This program and the accompanying materials // are made available under the terms of the GNU Public License v3.0 (or, at diff --git a/test/Crystallization.cpp b/test/Crystallization.cpp index 032200eee..e6b077c02 100644 --- a/test/Crystallization.cpp +++ b/test/Crystallization.cpp @@ -1,8 +1,8 @@ // ============================================================================= // CADET // -// Copyright © 2008-2024: 2008-present: The CADET-Core Authors -// Please see the AUTHORS and CONTRIBUTORS file. +// Copyright © 2008-present: The CADET-Core Authors +// Please see the AUTHORS.md file. // // All rights reserved. This program and the accompanying materials // are made available under the terms of the GNU Public License v3.0 (or, at