From 180d47e3b033d39c675d6ea87cb936895cbe77c3 Mon Sep 17 00:00:00 2001 From: Jan Breuer Date: Tue, 6 Aug 2024 17:40:04 +0200 Subject: [PATCH] Fix optional build of DG units --- src/libcadet/AdUtils.cpp | 204 +++++++++--------- src/libcadet/AdUtils.hpp | 67 +++--- src/libcadet/AutoDiff.hpp | 70 +++--- src/libcadet/linalg/Subset.hpp | 8 +- .../model/GeneralRateModelBuilder.cpp | 11 +- .../model/LumpedRateModelWithPoresBuilder.cpp | 10 +- .../LumpedRateModelWithoutPoresBuilder.cpp | 11 +- src/libcadet/model/ReactionModel.hpp | 4 +- src/libcadet/model/binding/DummyBinding.cpp | 2 + src/libcadet/model/binding/LinearBinding.cpp | 2 + test/CMakeLists.txt | 5 +- 11 files changed, 220 insertions(+), 174 deletions(-) diff --git a/src/libcadet/AdUtils.cpp b/src/libcadet/AdUtils.cpp index bf6b78efa..edd8986d2 100644 --- a/src/libcadet/AdUtils.cpp +++ b/src/libcadet/AdUtils.cpp @@ -73,65 +73,6 @@ void extractBandedJacobianFromAd(active const* const adVec, int adDirOffset, int } } -void extractBandedEigenJacobianFromAd(active const* const adVec, int adDirOffset, int diagDir, const int lowerBandwidth, const int upperBandwidth, Eigen::SparseMatrix& mat) -{ - const int stride = lowerBandwidth + 1 + upperBandwidth; - for (int eq = 0; eq < mat.rows(); ++eq) - { - // Start with lowest subdiagonal and stay in the range of the columns: - // diagDir might be too big for the matrix and, hence, dir ranges between - // diagDir - lowerBandwidth <= dir <= diagDir + upperBandwidth - int dir = diagDir - lowerBandwidth + eq % stride; - - // Loop over diagonals - for (int diag = 0; diag < stride; ++diag) - { - if (eq - lowerBandwidth + diag >= 0 && // left block boundary - eq - lowerBandwidth + diag < mat.rows() && // right block boundary - adVec[eq].getADValue(adDirOffset + dir) != 0.0 // keep pattern - ) - mat.coeffRef(eq, eq - lowerBandwidth + diag) = adVec[eq].getADValue(adDirOffset + dir); - - // Wrap around at end of row and jump to lowest subdiagonal - if (dir == diagDir + upperBandwidth) - dir = diagDir - lowerBandwidth; - else - ++dir; - } - } -} - -#ifdef ENABLE_DG - void extractBandedBlockEigenJacobianFromAd(active const* const adVec, int adDirOffset, int diagDir, const int lowerBandwidth, const int upperBandwidth, - const int blockOffset, const int nRows, Eigen::SparseMatrix& mat, const int matrixOffset) - { - const int stride = lowerBandwidth + 1 + upperBandwidth; - for (int eq = blockOffset; eq < blockOffset + nRows; ++eq) - { - // Start with lowest subdiagonal and stay in the range of the columns: - // diagDir might be too big for the matrix and, hence, dir ranges between - // diagDir - lowerBandwidth <= dir <= diagDir + upperBandwidth - int dir = diagDir - lowerBandwidth + eq % stride; - - // Loop over diagonals - for (int diag = 0; diag < stride; ++diag) - { - if (eq - lowerBandwidth + diag >= blockOffset && // left block boundary - eq - lowerBandwidth + diag < blockOffset + nRows && // right block boundary - adVec[eq].getADValue(adDirOffset + dir) != 0.0 // keep pattern - ) - mat.coeffRef(matrixOffset + eq, matrixOffset + eq - lowerBandwidth + diag) = adVec[eq].getADValue(adDirOffset + dir); - - // Wrap around at end of row and jump to lowest subdiagonal - if (dir == diagDir + upperBandwidth) - dir = diagDir - lowerBandwidth; - else - ++dir; - } - } - } -#endif - void prepareAdVectorSeedsForDenseMatrix(active* const adVec, int adDirOffset, int cols) { for (int col = 0; col < cols; ++col) @@ -285,63 +226,124 @@ double compareDenseJacobianWithBandedAd(active const* const adVec, int row, int return maxDiff; } -double compareBandedEigenJacobianWithAd(active const* const adVec, const int adDirOffset, const int diagDir, const int lowerBandwidth, const int upperBandwidth, - const int blockOffset, const int nRows, const Eigen::SparseMatrix& mat, const int matrixOffset) +void adMatrixVectorMultiply(const linalg::SparseMatrix& mat, double const* x, double* y, double alpha, double beta, int adDir) { - double maxDiff = 0.0; + const std::vector& rows = mat.rows(); + const std::vector& cols = mat.cols(); + const std::vector& values = mat.values(); + const int numNonZero = mat.numNonZero(); - const int stride = lowerBandwidth + 1 + upperBandwidth; - for (int eq = blockOffset; eq < blockOffset + nRows; ++eq) + for (int i = 0; i < numNonZero; ++i) { - // Start with lowest subdiagonal and stay in the range of the columns: - // diagDir might be too big for the matrix and, hence, dir ranges between - // diagDir - lowerBandwidth <= dir <= diagDir + upperBandwidth - int dir = diagDir - lowerBandwidth + eq % stride; + y[rows[i]] = alpha * values[i].getADValue(adDir) * x[cols[i]] + beta * y[rows[i]]; + } +} - // Loop over diagonals - for (int diag = 0; diag < stride; ++diag) +#ifdef ENABLE_DG + + void extractBandedEigenJacobianFromAd(active const* const adVec, int adDirOffset, int diagDir, const int lowerBandwidth, const int upperBandwidth, Eigen::SparseMatrix& mat) + { + const int stride = lowerBandwidth + 1 + upperBandwidth; + for (int eq = 0; eq < mat.rows(); ++eq) { - if (eq - lowerBandwidth + diag >= blockOffset && // left block boundary - eq - lowerBandwidth + diag < blockOffset + nRows && // right block boundary - adVec[eq].getADValue(adDirOffset + dir) != 0.0 // keep pattern - ) - { - double baseVal = adVec[eq].getADValue(adDirOffset + dir); - double matVal = mat.coeff(matrixOffset + eq, matrixOffset + eq - lowerBandwidth + diag); + // Start with lowest subdiagonal and stay in the range of the columns: + // diagDir might be too big for the matrix and, hence, dir ranges between + // diagDir - lowerBandwidth <= dir <= diagDir + upperBandwidth + int dir = diagDir - lowerBandwidth + eq % stride; - if (std::isnan(matVal) || std::isnan(baseVal)) - return std::numeric_limits::quiet_NaN(); - const double diff = std::abs(matVal - baseVal); + // Loop over diagonals + for (int diag = 0; diag < stride; ++diag) + { + if (eq - lowerBandwidth + diag >= 0 && // left block boundary + eq - lowerBandwidth + diag < mat.rows() && // right block boundary + adVec[eq].getADValue(adDirOffset + dir) != 0.0 // keep pattern + ) + mat.coeffRef(eq, eq - lowerBandwidth + diag) = adVec[eq].getADValue(adDirOffset + dir); - baseVal = std::abs(baseVal); - if (baseVal > 0.0) - maxDiff = std::max(maxDiff, diff / baseVal); + // Wrap around at end of row and jump to lowest subdiagonal + if (dir == diagDir + upperBandwidth) + dir = diagDir - lowerBandwidth; else - maxDiff = std::max(maxDiff, diff); + ++dir; } - - // Wrap around at end of row and jump to lowest subdiagonal - if (dir == diagDir + upperBandwidth) - dir = diagDir - lowerBandwidth; - else - ++dir; } } - return maxDiff; -} -void adMatrixVectorMultiply(const linalg::SparseMatrix& mat, double const* x, double* y, double alpha, double beta, int adDir) -{ - const std::vector& rows = mat.rows(); - const std::vector& cols = mat.cols(); - const std::vector& values = mat.values(); - const int numNonZero = mat.numNonZero(); + void extractBandedBlockEigenJacobianFromAd(active const* const adVec, int adDirOffset, int diagDir, const int lowerBandwidth, const int upperBandwidth, + const int blockOffset, const int nRows, Eigen::SparseMatrix& mat, const int matrixOffset) + { + const int stride = lowerBandwidth + 1 + upperBandwidth; + for (int eq = blockOffset; eq < blockOffset + nRows; ++eq) + { + // Start with lowest subdiagonal and stay in the range of the columns: + // diagDir might be too big for the matrix and, hence, dir ranges between + // diagDir - lowerBandwidth <= dir <= diagDir + upperBandwidth + int dir = diagDir - lowerBandwidth + eq % stride; - for (int i = 0; i < numNonZero; ++i) + // Loop over diagonals + for (int diag = 0; diag < stride; ++diag) + { + if (eq - lowerBandwidth + diag >= blockOffset && // left block boundary + eq - lowerBandwidth + diag < blockOffset + nRows && // right block boundary + adVec[eq].getADValue(adDirOffset + dir) != 0.0 // keep pattern + ) + mat.coeffRef(matrixOffset + eq, matrixOffset + eq - lowerBandwidth + diag) = adVec[eq].getADValue(adDirOffset + dir); + + // Wrap around at end of row and jump to lowest subdiagonal + if (dir == diagDir + upperBandwidth) + dir = diagDir - lowerBandwidth; + else + ++dir; + } + } + } + + double compareBandedEigenJacobianWithAd(active const* const adVec, const int adDirOffset, const int diagDir, const int lowerBandwidth, const int upperBandwidth, + const int blockOffset, const int nRows, const Eigen::SparseMatrix& mat, const int matrixOffset) { - y[rows[i]] = alpha * values[i].getADValue(adDir) * x[cols[i]] + beta * y[rows[i]]; + double maxDiff = 0.0; + + const int stride = lowerBandwidth + 1 + upperBandwidth; + for (int eq = blockOffset; eq < blockOffset + nRows; ++eq) + { + // Start with lowest subdiagonal and stay in the range of the columns: + // diagDir might be too big for the matrix and, hence, dir ranges between + // diagDir - lowerBandwidth <= dir <= diagDir + upperBandwidth + int dir = diagDir - lowerBandwidth + eq % stride; + + // Loop over diagonals + for (int diag = 0; diag < stride; ++diag) + { + if (eq - lowerBandwidth + diag >= blockOffset && // left block boundary + eq - lowerBandwidth + diag < blockOffset + nRows && // right block boundary + adVec[eq].getADValue(adDirOffset + dir) != 0.0 // keep pattern + ) + { + double baseVal = adVec[eq].getADValue(adDirOffset + dir); + double matVal = mat.coeff(matrixOffset + eq, matrixOffset + eq - lowerBandwidth + diag); + + if (std::isnan(matVal) || std::isnan(baseVal)) + return std::numeric_limits::quiet_NaN(); + const double diff = std::abs(matVal - baseVal); + + baseVal = std::abs(baseVal); + if (baseVal > 0.0) + maxDiff = std::max(maxDiff, diff / baseVal); + else + maxDiff = std::max(maxDiff, diff); + } + + // Wrap around at end of row and jump to lowest subdiagonal + if (dir == diagDir + upperBandwidth) + dir = diagDir - lowerBandwidth; + else + ++dir; + } + } + return maxDiff; } -} + +#endif } // namespace ad diff --git a/src/libcadet/AdUtils.hpp b/src/libcadet/AdUtils.hpp index 4db6c1cd6..666e14a00 100644 --- a/src/libcadet/AdUtils.hpp +++ b/src/libcadet/AdUtils.hpp @@ -85,20 +85,41 @@ void extractBandedJacobianFromAd(active const* const adVec, int adDirOffset, int */ void extractBandedBlockEigenJacobianFromAd(active const* const adVec, int adDirOffset, int diagDir, const int lowerBandwidth, const int upperBandwidth, const int blockOffset, const int nCols, Eigen::SparseMatrix& mat, const int matrixOffset = 0); -#endif -/** - * @brief Extracts a band matrix (Eigen lib) from band compressed AD seed vectors - * @details Uses the results of an AD computation with seed vectors set by prepareAdVectorSeedsForBandMatrix() to - assemble the Jacobian block which is a band matrix. The block must be on the main diagonal. - * @param [in] adVec Vector of AD datatypes with band compressed seed vectors - * @param [in] adDirOffset Offset in the AD directions (can be used to move past parameter sensitivity directions) - * @param [in] diagDir Diagonal direction index - * @param [in] lowerBandwidth lower band width - * @param [in] upperBandwidth upper band width - * @param [out] mat Eigen matrix to be populated with the Jacobian block - */ -void extractBandedEigenJacobianFromAd(active const* const adVec, int adDirOffset, int diagDir, const int lowerBandwidth, const int upperBandwidth, Eigen::SparseMatrix& mat); + /** + * @brief Extracts a band matrix (Eigen lib) from band compressed AD seed vectors + * @details Uses the results of an AD computation with seed vectors set by prepareAdVectorSeedsForBandMatrix() to + assemble the Jacobian block which is a band matrix. The block must be on the main diagonal. + * @param [in] adVec Vector of AD datatypes with band compressed seed vectors + * @param [in] adDirOffset Offset in the AD directions (can be used to move past parameter sensitivity directions) + * @param [in] diagDir Diagonal direction index + * @param [in] lowerBandwidth lower band width + * @param [in] upperBandwidth upper band width + * @param [out] mat Eigen matrix to be populated with the Jacobian block + */ + void extractBandedEigenJacobianFromAd(active const* const adVec, int adDirOffset, int diagDir, const int lowerBandwidth, const int upperBandwidth, Eigen::SparseMatrix& mat); + + /** + * @brief Compares a (block-) banded Jacobian in Eigen library row-major format with an AD version derived by band compressed AD seed vectors + * @details Uses the results of an AD computation with seed vectors set by prepareAdVectorSeedsForBandMatrix() to + compare the results with a given banded Jacobian. The AD Jacobian is treated as base and the analytic + Jacobian is compared against it. The relative difference + @f[ \Delta_{ij} = \begin{cases} \left\lvert \frac{ J_{\text{ana},ij} - J_{\text{ad},ij} }{ J_{\text{ad},ij} }\right\rvert, & J_{\text{ad},ij} \neq 0 \\ + \left\lvert J_{\text{ana},ij} - J_{\text{ad},ij} \right\rvert, & J_{\text{ad},ij} = 0 \end{cases} @f] + is computed for each matrix entry. The maximum of all @f$ \Delta_{ij} @f$ is returned. + * @param [in] adVec Vector of AD datatypes with band compressed seed vectors + * @param [in] adDirOffset Offset in the AD directions (can be used to move past parameter sensitivity directions) + * @param [in] diagDir Diagonal direction index + * @param [in] lowerBandwidth of the block + * @param [in] upperBandwidth of the block + * @param [in] blockOffset of the currently considered block + * @param [in] nRows of the matrix or number of equations to be compared + * @param [in] mat BandMatrix populated with the analytic Jacobian + * @return The maximum absolute relative difference between the matrix elements + */ + double compareBandedEigenJacobianWithAd(active const* const adVec, const int adDirOffset, const int diagDir, const int lowerBandwidth, const int upperBandwidth, const int blockOffset, const int nRows, const Eigen::SparseMatrix& mat, const int matrixOffset); + +#endif /** * @brief Sets seed vectors on an AD vector for computing a dense Jacobian @@ -151,26 +172,6 @@ void extractDenseJacobianFromBandedAd(active const* const adVec, int row, int ad */ double compareBandedJacobianWithAd(active const* const adVec, int adDirOffset, int diagDir, const linalg::BandMatrix& mat); -/** - * @brief Compares a (block-) banded Jacobian in Eigen library row-major format with an AD version derived by band compressed AD seed vectors - * @details Uses the results of an AD computation with seed vectors set by prepareAdVectorSeedsForBandMatrix() to - compare the results with a given banded Jacobian. The AD Jacobian is treated as base and the analytic - Jacobian is compared against it. The relative difference - @f[ \Delta_{ij} = \begin{cases} \left\lvert \frac{ J_{\text{ana},ij} - J_{\text{ad},ij} }{ J_{\text{ad},ij} }\right\rvert, & J_{\text{ad},ij} \neq 0 \\ - \left\lvert J_{\text{ana},ij} - J_{\text{ad},ij} \right\rvert, & J_{\text{ad},ij} = 0 \end{cases} @f] - is computed for each matrix entry. The maximum of all @f$ \Delta_{ij} @f$ is returned. - * @param [in] adVec Vector of AD datatypes with band compressed seed vectors - * @param [in] adDirOffset Offset in the AD directions (can be used to move past parameter sensitivity directions) - * @param [in] diagDir Diagonal direction index - * @param [in] lowerBandwidth of the block - * @param [in] upperBandwidth of the block - * @param [in] blockOffset of the currently considered block - * @param [in] nRows of the matrix or number of equations to be compared - * @param [in] mat BandMatrix populated with the analytic Jacobian - * @return The maximum absolute relative difference between the matrix elements - */ -double compareBandedEigenJacobianWithAd(active const* const adVec, const int adDirOffset, const int diagDir, const int lowerBandwidth, const int upperBandwidth, const int blockOffset, const int nRows, const Eigen::SparseMatrix& mat, const int matrixOffset); - /** * @brief Compares a dense Jacobian with an AD version derived by AD seed vectors * @details Uses the results of an AD computation with seed vectors set by prepareAdVectorSeedsForDenseMatrix() to diff --git a/src/libcadet/AutoDiff.hpp b/src/libcadet/AutoDiff.hpp index 437dbaace..be4803dbc 100644 --- a/src/libcadet/AutoDiff.hpp +++ b/src/libcadet/AutoDiff.hpp @@ -17,10 +17,14 @@ #ifndef LIBCADET_AUTODIFF_HPP_ #define LIBCADET_AUTODIFF_HPP_ +#include "CompileTimeConfig.hpp" #include "cadet/cadetCompilerInfo.hpp" #include "common/CompilerSpecific.hpp" -#include +#ifdef ENABLE_DG + #include +#endif + #include #if defined(ACTIVE_SFAD) || defined(ACTIVE_SETFAD) @@ -163,39 +167,41 @@ namespace cadet struct ActiveRefOrDouble { typedef const double type; }; } -namespace Eigen { +#ifdef ENABLE_DG + namespace Eigen { - template<> struct NumTraits - : NumTraits // permits to get the epsilon, dummy_precision, lowest, highest functions - { - typedef cadet::active Real; - typedef cadet::active NonInteger; - typedef cadet::active Nested; - - enum { - IsComplex = 0, - IsInteger = 0, - IsSigned = 1, - RequireInitialization = 1, - ReadCost = 1, - AddCost = 3, - MulCost = 3 + template<> struct NumTraits + : NumTraits // permits to get the epsilon, dummy_precision, lowest, highest functions + { + typedef cadet::active Real; + typedef cadet::active NonInteger; + typedef cadet::active Nested; + + enum { + IsComplex = 0, + IsInteger = 0, + IsSigned = 1, + RequireInitialization = 1, + ReadCost = 1, + AddCost = 3, + MulCost = 3 + }; }; - }; - // specify return types concerning active double scalar operations for eigen - template<> - struct ScalarBinaryOpTraits> { - typedef cadet::active ReturnType; - }; - template<> - struct ScalarBinaryOpTraits> { - typedef cadet::active ReturnType; - }; - template<> - struct ScalarBinaryOpTraits> { - typedef cadet::active ReturnType; - }; -} + // specify return types concerning active double scalar operations for eigen + template<> + struct ScalarBinaryOpTraits> { + typedef cadet::active ReturnType; + }; + template<> + struct ScalarBinaryOpTraits> { + typedef cadet::active ReturnType; + }; + template<> + struct ScalarBinaryOpTraits> { + typedef cadet::active ReturnType; + }; + } +#endif #endif // LIBCADET_AUTODIFF_HPP_ diff --git a/src/libcadet/linalg/Subset.hpp b/src/libcadet/linalg/Subset.hpp index de937cec1..3dfa772f2 100644 --- a/src/libcadet/linalg/Subset.hpp +++ b/src/libcadet/linalg/Subset.hpp @@ -25,8 +25,10 @@ #include "linalg/BandMatrix.hpp" #include -#include -#include +#ifdef ENABLE_DG + #include + #include +#endif namespace cadet { @@ -515,6 +517,7 @@ namespace linalg } } +#ifdef ENABLE_DG /** * @brief Copies a subset of a given Eigen matrix into this matrix * @details Copies a submatrix indentified by row and column masks from a given source in banded @@ -582,6 +585,7 @@ namespace linalg ptrDest += dest.stride(); } } +#endif } // namespace linalg diff --git a/src/libcadet/model/GeneralRateModelBuilder.cpp b/src/libcadet/model/GeneralRateModelBuilder.cpp index 9f55214ca..6a85e7571 100644 --- a/src/libcadet/model/GeneralRateModelBuilder.cpp +++ b/src/libcadet/model/GeneralRateModelBuilder.cpp @@ -1,5 +1,8 @@ #include "model/GeneralRateModel.hpp" +#include "CompileTimeConfig.hpp" +#ifdef ENABLE_DG #include "model/GeneralRateModelDG.hpp" +#endif #include "LoggingUtils.hpp" #include "Logging.hpp" @@ -18,10 +21,13 @@ namespace cadet if (paramProvider.exists("SPATIAL_METHOD")) { const std::string discName = paramProvider.getString("SPATIAL_METHOD"); - +#ifdef ENABLE_DG if (discName == "DG") model = new GeneralRateModelDG(uoId); else if (discName == "FV") +#else + if (discName == "FV") +#endif model = createAxialFVGRM(uoId); else { @@ -72,9 +78,10 @@ namespace cadet typedef GeneralRateModel AxialGRM; typedef GeneralRateModel RadialGRM; +#ifdef ENABLE_DG models[GeneralRateModelDG::identifier()] = selectAxialFlowDiscretizationGRM; models["GRM_DG"] = selectAxialFlowDiscretizationGRM; - +#endif models[AxialGRM::identifier()] = selectAxialFlowDiscretizationGRM; models["GRM"] = selectAxialFlowDiscretizationGRM; diff --git a/src/libcadet/model/LumpedRateModelWithPoresBuilder.cpp b/src/libcadet/model/LumpedRateModelWithPoresBuilder.cpp index 2b9c459c4..c713f85f9 100644 --- a/src/libcadet/model/LumpedRateModelWithPoresBuilder.cpp +++ b/src/libcadet/model/LumpedRateModelWithPoresBuilder.cpp @@ -1,5 +1,8 @@ #include "model/LumpedRateModelWithPores.hpp" +#include "CompileTimeConfig.hpp" +#ifdef ENABLE_DG #include "model/LumpedRateModelWithPoresDG.hpp" +#endif #include "LoggingUtils.hpp" #include "Logging.hpp" @@ -19,9 +22,13 @@ namespace cadet const std::string discName = paramProvider.getString("SPATIAL_METHOD"); +#ifdef ENABLE_DG if (discName == "DG") model = new LumpedRateModelWithPoresDG(uoId); else if (discName == "FV") +#else + if (discName == "FV") +#endif model = createAxialFVLRMP(uoId); else { @@ -72,9 +79,10 @@ namespace cadet typedef LumpedRateModelWithPores AxialLRMP; typedef LumpedRateModelWithPores RadialLRMP; +#ifdef ENABLE_DG models[LumpedRateModelWithPoresDG::identifier()] = selectAxialFlowDiscretizationLRMP; models["LRMP_DG"] = selectAxialFlowDiscretizationLRMP; - +#endif models[AxialLRMP::identifier()] = selectAxialFlowDiscretizationLRMP; models["LRMP"] = selectAxialFlowDiscretizationLRMP; diff --git a/src/libcadet/model/LumpedRateModelWithoutPoresBuilder.cpp b/src/libcadet/model/LumpedRateModelWithoutPoresBuilder.cpp index 0b5a9ad8c..5cf23237b 100644 --- a/src/libcadet/model/LumpedRateModelWithoutPoresBuilder.cpp +++ b/src/libcadet/model/LumpedRateModelWithoutPoresBuilder.cpp @@ -1,5 +1,8 @@ #include "model/LumpedRateModelWithoutPores.hpp" -#include "model/LumpedRateModelWithoutPoresDG.hpp" +#include "CompileTimeConfig.hpp" +#ifdef ENABLE_DG + #include "model/LumpedRateModelWithoutPoresDG.hpp" +#endif #include "LoggingUtils.hpp" #include "Logging.hpp" @@ -19,9 +22,13 @@ namespace model const std::string discName = paramProvider.getString("SPATIAL_METHOD"); +#ifdef ENABLE_DG if(discName == "DG") model = new LumpedRateModelWithoutPoresDG(uoId); else if (discName == "FV") +#else + if (discName == "FV") +#endif model = createAxialFVLRM(uoId); else { @@ -72,8 +79,10 @@ void registerLumpedRateModelWithoutPores(std::unordered_map AxialLRM; typedef LumpedRateModelWithoutPores RadialLRM; +#ifdef ENABLE_DG models[LumpedRateModelWithoutPoresDG::identifier()] = selectAxialFlowDiscretizationLRM; models["LRM_DG"] = selectAxialFlowDiscretizationLRM; +#endif models[AxialLRM::identifier()] = selectAxialFlowDiscretizationLRM; models["LRM"] = selectAxialFlowDiscretizationLRM; diff --git a/src/libcadet/model/ReactionModel.hpp b/src/libcadet/model/ReactionModel.hpp index eeadf9442..3555ac64f 100644 --- a/src/libcadet/model/ReactionModel.hpp +++ b/src/libcadet/model/ReactionModel.hpp @@ -25,7 +25,9 @@ #include "cadet/ParameterId.hpp" #include "linalg/DenseMatrix.hpp" #include "linalg/BandMatrix.hpp" -#include "linalg/BandedEigenSparseRowIterator.hpp" +#ifdef ENABLE_DG + #include "linalg/BandedEigenSparseRowIterator.hpp" +#endif #include "linalg/CompressedSparseMatrix.hpp" #include "AutoDiff.hpp" #include "Memory.hpp" diff --git a/src/libcadet/model/binding/DummyBinding.cpp b/src/libcadet/model/binding/DummyBinding.cpp index f971a955d..34b81637b 100644 --- a/src/libcadet/model/binding/DummyBinding.cpp +++ b/src/libcadet/model/binding/DummyBinding.cpp @@ -131,9 +131,11 @@ class DummyBinding : public IBindingModel { } +#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 { } +#endif virtual void timeDerivativeQuasiStationaryFluxes(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yCp, double const* y, double* dResDt, LinearBufferAllocator workSpace) const { } diff --git a/src/libcadet/model/binding/LinearBinding.cpp b/src/libcadet/model/binding/LinearBinding.cpp index 7f494af99..4d481da5b 100644 --- a/src/libcadet/model/binding/LinearBinding.cpp +++ b/src/libcadet/model/binding/LinearBinding.cpp @@ -556,10 +556,12 @@ class LinearBindingBase : public IBindingModel jacobianImpl(t, secIdx, colPos, y, offsetCp, jac, workSpace); } +#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 { jacobianImpl(t, secIdx, colPos, y, offsetCp, jac, workSpace); } +#endif virtual void timeDerivativeQuasiStationaryFluxes(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yCp, double const* y, double* dResDt, LinearBufferAllocator workSpace) const { diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 71f282549..eb89e6c74 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -61,9 +61,12 @@ if (ENABLE_2D_MODELS) list(APPEND TEST_ADDITIONAL_SOURCES SparseFactorizableMatrix.cpp TwoDimConvectionDispersionOperator.cpp) endif() +if (ENABLE_DG) + list(APPEND TEST_ADDITIONAL_SOURCES GeneralRateModelDG.cpp LumpedRateModelWithPoresDG.cpp LumpedRateModelWithoutPoresDG.cpp) +endif() + add_executable(testRunner testRunner.cpp JsonTestModels.cpp ColumnTests.cpp UnitOperationTests.cpp SimHelper.cpp ParticleHelper.cpp GeneralRateModel.cpp GeneralRateModel2D.cpp LumpedRateModelWithPores.cpp LumpedRateModelWithoutPores.cpp - GeneralRateModelDG.cpp LumpedRateModelWithPoresDG.cpp LumpedRateModelWithoutPoresDG.cpp RadialGeneralRateModel.cpp RadialLumpedRateModelWithPores.cpp RadialLumpedRateModelWithoutPores.cpp MultiChannelTransportModel.cpp CSTR-Residual.cpp CSTR-Simulation.cpp