Fix test interdependencies

src/libcadet/linalg/BandedEigenSparseRowIterator.hpp

@@ -29,7 +29,7 @@ class BandedEigenSparseRowIterator
 	/**
 	 * @brief Creates a ConstBandedSparseRowIterator pointing nowhere
 	 */
-	BandedEigenSparseRowIterator() : _matrix(nullptr), _values(nullptr), _colIdx(nullptr), _row(-1), _numNonZero(0), _dummy(0.0) { }
+	BandedEigenSparseRowIterator() : _matrix(nullptr), _values(nullptr), _colIdx(nullptr), _row(-1), _numNonZero(0), _dummy(0.0), _size(0) { }
 
 	/**
 	 * @brief Creates a BandedSparseRowIterator of the given matrix row
@@ -38,8 +38,12 @@ class BandedEigenSparseRowIterator
 	 */
 	BandedEigenSparseRowIterator(Eigen::SparseMatrix<double, 0x1>& mat, int row)
 		: _matrix(&mat), _values(valuesOfRow(mat, row)), _colIdx(columnIndicesOfRow(mat, row)), _row(row),
-		_numNonZero(getInnerNumberOfNonZeros(mat, row)), _dummy(0.0)
+		_dummy(0.0), _size(mat.rows())
 	{
+		if (row < _size)
+			_numNonZero = getInnerNumberOfNonZeros(mat, row);
+		else
+			_numNonZero = 0;
 	}
 
 	~BandedEigenSparseRowIterator() CADET_NOEXCEPT { }
@@ -117,7 +121,7 @@ class BandedEigenSparseRowIterator
 	 */
 	inline double& native(int col)
 	{
-		cadet_assert((col >= 0) && (col < _matrix->rows()));
+		cadet_assert((col >= 0) && (col < _size));
 
 		// Try to find the element
 		// TODO: Use binary search
@@ -156,29 +160,41 @@ class BandedEigenSparseRowIterator
 
 	inline BandedEigenSparseRowIterator& operator++() CADET_NOEXCEPT
 	{
-		++_row;
-		updateOnRowChange();
+		if (_row + 1 < _size)
+		{
+			++_row;
+			updateOnRowChange();
+		}
 		return *this;
 	}
 
 	inline BandedEigenSparseRowIterator& operator--() CADET_NOEXCEPT
 	{
-		--_row;
-		updateOnRowChange();
+		if (_row - 1 > 0)
+		{
+			--_row;
+			updateOnRowChange();
+		}
 		return *this;
 	}
 
 	inline BandedEigenSparseRowIterator& operator+=(int idx) CADET_NOEXCEPT
 	{
-		_row += idx;
-		updateOnRowChange();
+		if (_row + idx < _size)
+		{
+			_row += idx;
+			updateOnRowChange();
+		}
 		return *this;
 	}
 
 	inline BandedEigenSparseRowIterator& operator-=(int idx) CADET_NOEXCEPT
 	{
-		_row -= idx;
-		updateOnRowChange();
+		if (_row - idx > 0)
+		{
+			_row -= idx;
+			updateOnRowChange();
+		}
 		return *this;
 	}
 
@@ -248,12 +264,13 @@ class BandedEigenSparseRowIterator
 		return mat.outerIndexPtr()[row + 1] - mat.outerIndexPtr()[row];
 	}
 
-	Eigen::SparseMatrix<double, 0x1>* _matrix;
-	double* _values;
-	int const* _colIdx;
-	int _row;
-	int _numNonZero;
-	double _dummy;
+	Eigen::SparseMatrix<double, 0x1>* _matrix; //<! Eigen library sparse matrix (should be square)
+	double* _values; //<! values of current row
+	int const* _colIdx; //<! column indices of entries in current row
+	int _row; //<! index of current row
+	int _numNonZero; //<! number of non-zeros in current row
+	double _dummy; //<! dummy value thats written to if entry is not found
+	int _size; //<! size (number of rows) of matrix
 };
 
 } // namespace linalg

src/libcadet/model/GeneralRateModel.cpp

@@ -87,14 +87,6 @@ GeneralRateModel<ConvDispOperator>::~GeneralRateModel() CADET_NOEXCEPT
 	delete _dynReactionBulk;
 
 	clearParDepSurfDiffusion();
-
-	delete[] _disc.nParCell;
-	delete[] _disc.parTypeOffset;
-	delete[] _disc.nParCellsBeforeType;
-	delete[] _disc.nBound;
-	delete[] _disc.boundOffset;
-	delete[] _disc.strideBound;
-	delete[] _disc.nBoundBeforeType;
 }
 
 template <typename ConvDispOperator>

src/libcadet/model/GeneralRateModel.hpp

@@ -290,6 +290,17 @@ class GeneralRateModel : public UnitOperationBase
 		unsigned int* boundOffset; //!< Array with offset to the first bound state of each component in the solid phase (particle type major ordering)
 		unsigned int* strideBound; //!< Total number of bound states for each particle type, additional last element contains total number of bound states for all types
 		unsigned int* nBoundBeforeType; //!< Array with number of bound states before a particle type (cumulative sum of strideBound)
+
+		~Discretization() // make sure this memory is freed correctly
+		{
+			delete[] nParCell;
+			delete[] parTypeOffset;
+			delete[] nParCellsBeforeType;
+			delete[] nBound;
+			delete[] boundOffset;
+			delete[] strideBound;
+			delete[] nBoundBeforeType;
+		}
 	};
 
 	enum class ParticleDiscretizationMode : int

src/libcadet/model/GeneralRateModelDG-InitialConditions.cpp

@@ -713,7 +713,7 @@ void GeneralRateModelDG::consistentInitialTimeDerivative(const SimulationTime& s
 
 			// Get iterators to beginning of solid phase
 			linalg::BandedEigenSparseRowIterator jacSolidOrig(_globalJac, idxr.offsetCp(ParticleTypeIndex{ type }, ParticleIndex{ par }) + j * static_cast<unsigned int>(idxr.strideParNode(type)) + static_cast<unsigned int>(idxr.strideParLiquid()));
-			linalg::BandedEigenSparseRowIterator jacSolid = jacPar - idxr.strideParBound(type);
+			linalg::BandedEigenSparseRowIterator jacSolid(_globalJacDisc, idxr.offsetCp(ParticleTypeIndex{ type }, ParticleIndex{ par }) + j * static_cast<unsigned int>(idxr.strideParNode(type)) + static_cast<unsigned int>(idxr.strideParLiquid()));
 
 			int const* const mask = _binding[type]->reactionQuasiStationarity();
 			double* const qShellDot = vecStateYdot + idxr.offsetCp(ParticleTypeIndex{ type }, ParticleIndex{ par }) + static_cast<int>(j) * idxr.strideParNode(type) + idxr.strideParLiquid();
@@ -1146,7 +1146,7 @@ void GeneralRateModelDG::consistentInitialSensitivity(const SimulationTime& simT
 				{
 					// Get iterators to beginning of solid phase
 					linalg::BandedEigenSparseRowIterator jacSolidOrig(_globalJac, idxr.offsetCp(ParticleTypeIndex{ static_cast<unsigned int>(type) }, ParticleIndex{ static_cast<unsigned int>(pblk) }) + j * idxr.strideParNode(type) + idxr.strideParLiquid());
-					linalg::BandedEigenSparseRowIterator jacSolid = jacPar - idxr.strideParBound(type);
+					linalg::BandedEigenSparseRowIterator jacSolid(_globalJacDisc, idxr.offsetCp(ParticleTypeIndex{ static_cast<unsigned int>(type) }, ParticleIndex{ static_cast<unsigned int>(pblk) }) + j * idxr.strideParNode(type) + idxr.strideParLiquid());
 
 					int const* const mask = _binding[type]->reactionQuasiStationarity();
 					double* const qShellDot = sensYdot + idxr.offsetCp(ParticleTypeIndex{ type }, ParticleIndex{ par }) + static_cast<int>(j) * idxr.strideParNode(type) + idxr.strideParLiquid();

src/libcadet/model/GeneralRateModelDG.cpp

@@ -55,7 +55,7 @@ constexpr double SurfVolRatioSlab = 1.0;
 
 
 GeneralRateModelDG::GeneralRateModelDG(UnitOpIdx unitOpIdx) : UnitOperationBase(unitOpIdx),
-	_hasSurfaceDiffusion(0, false), _dynReactionBulk(nullptr),
+	_hasSurfaceDiffusion(0, false), _singleParDepSurfDiffusion(false), _dynReactionBulk(nullptr),
 	_globalJac(), _globalJacDisc(), _jacInlet(), _hasParDepSurfDiffusion(false),
 	_analyticJac(true), _jacobianAdDirs(0), _factorizeJacobian(false), _tempState(nullptr),
 	_initC(0), _initCp(0), _initQ(0), _initState(0), _initStateDot(0)
@@ -70,37 +70,6 @@ GeneralRateModelDG::~GeneralRateModelDG() CADET_NOEXCEPT
 
 	clearParDepSurfDiffusion();
 
-	delete[] _disc.nParCell;
-	delete[] _disc.nParPointsBeforeType;
-	delete[] _disc.parPolyDeg;
-	delete[] _disc.nParNode;
-	delete[] _disc.nParPoints;
-	delete[] _disc.parExactInt;
-	delete[] _disc.parTypeOffset;
-	delete[] _disc.nBound;
-	delete[] _disc.boundOffset;
-	delete[] _disc.strideBound;
-	delete[] _disc.nBoundBeforeType;
-
-	delete[] _disc.offsetSurfDiff;
-
-	delete[] _disc.deltaR;
-	delete[] _disc.parNodes;
-	delete[] _disc.parPolyDerM;
-	delete[] _disc.minus_InvMM_ST;
-	delete[] _disc.parInvWeights;
-	delete[] _disc.parInvMM;
-	delete[] _disc.parInvMM_Leg;
-	delete[] _disc.Ir;
-	delete[] _disc.Dr;
-
-	delete[] _disc.DGjacParDispBlocks;
-
-	delete[] _disc.g_p;
-	delete[] _disc.g_pSum;
-	delete[] _disc.surfaceFluxParticle;
-	delete[] _disc.localFlux;
-
 	delete _linearSolver;
 }
 
@@ -160,7 +129,10 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP
 
 	paramProvider.pushScope("discretization");
 
-	_linearSolver = cadet::linalg::setLinearSolver(paramProvider.exists("LINEAR_SOLVER") ? paramProvider.getString("LINEAR_SOLVER") : "SparseLU");
+	const bool firstConfigCall = _tempState == nullptr; // used to not multiply allocate memory
+
+	if (firstConfigCall)
+		_linearSolver = cadet::linalg::setLinearSolver(paramProvider.exists("LINEAR_SOLVER") ? paramProvider.getString("LINEAR_SOLVER") : "SparseLU");
 
 	if (!newNBoundInterface && paramProvider.exists("NBOUND")) // done here and in this order for backwards compatibility
 		nBound = paramProvider.getIntArray("NBOUND");
@@ -215,10 +187,13 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP
 	std::vector<int> parPolyDeg(_disc.nParType);
 	std::vector<int> ParNelem(_disc.nParType);
 	std::vector<bool> parExactInt(_disc.nParType, true);
-	_disc.parPolyDeg = new unsigned int[_disc.nParType];
-	_disc.nParCell = new unsigned int[_disc.nParType];
-	_disc.parExactInt = new bool[_disc.nParType];
-	_disc.parGSM = new bool[_disc.nParType];
+	if (firstConfigCall)
+	{
+		_disc.parPolyDeg = new unsigned int[_disc.nParType];
+		_disc.nParCell = new unsigned int[_disc.nParType];
+		_disc.parExactInt = new bool[_disc.nParType];
+		_disc.parGSM = new bool[_disc.nParType];
+	}
 
 	if (paramProvider.exists("PAR_POLYDEG"))
 	{
@@ -300,7 +275,8 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP
 	// Compute discretization operators and initialize containers
 	_disc.initializeDG();
 
-	_disc.nBound = new unsigned int[_disc.nComp * _disc.nParType];
+	if (firstConfigCall)
+		_disc.nBound = new unsigned int[_disc.nComp * _disc.nParType];
 	if (nBound.size() < _disc.nComp * _disc.nParType)
 	{
 		// Multiplex number of bound states to all particle types
@@ -313,9 +289,12 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP
 	const unsigned int nTotalBound = std::accumulate(_disc.nBound, _disc.nBound + _disc.nComp * _disc.nParType, 0u);
 
 	// Precompute offsets and total number of bound states (DOFs in solid phase)
-	_disc.boundOffset = new unsigned int[_disc.nComp * _disc.nParType];
-	_disc.strideBound = new unsigned int[_disc.nParType + 1];
-	_disc.nBoundBeforeType = new unsigned int[_disc.nParType];
+	if (firstConfigCall)
+	{
+		_disc.boundOffset = new unsigned int[_disc.nComp * _disc.nParType];
+		_disc.strideBound = new unsigned int[_disc.nParType + 1];
+		_disc.nBoundBeforeType = new unsigned int[_disc.nParType];
+	}
 	_disc.strideBound[_disc.nParType] = nTotalBound;
 	_disc.nBoundBeforeType[0] = 0;
 	for (unsigned int j = 0; j < _disc.nParType; ++j)
@@ -335,8 +314,11 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP
 	}
 
 	// Precompute offsets of particle type DOFs
-	_disc.parTypeOffset = new unsigned int[_disc.nParType + 1];
-	_disc.nParPointsBeforeType = new unsigned int[_disc.nParType + 1];
+	if (firstConfigCall)
+	{
+		_disc.parTypeOffset = new unsigned int[_disc.nParType + 1];
+		_disc.nParPointsBeforeType = new unsigned int[_disc.nParType + 1];
+	}
 	_disc.parTypeOffset[0] = 0;
 	_disc.nParPointsBeforeType[0] = 0;
 	unsigned int nTotalParPoints = 0;
@@ -629,7 +611,8 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP
 	}
 
 	// Allocate memory
-	_tempState = new double[numDofs()];
+	if (firstConfigCall)
+		_tempState = new double[numDofs()];
 
 	if (_disc.exactInt)
 		_jacInlet.resize(_disc.nNodes, 1); // first cell depends on inlet concentration (same for every component)
@@ -651,6 +634,8 @@ bool GeneralRateModelDG::configure(IParameterProvider& paramProvider)
 {
 	_parameters.clear();
 
+	const bool firstConfigCall = _disc.offsetSurfDiff == nullptr; // used to not multiply allocate memory
+
 	const bool transportSuccess = _convDispOp.configure(_unitOpIdx, paramProvider, _parameters);
 
 	// Read geometry parameters
@@ -716,7 +701,8 @@ bool GeneralRateModelDG::configure(IParameterProvider& paramProvider)
 	_filmDiffusionMode = readAndRegisterMultiplexCompTypeSecParam(paramProvider, _parameters, _filmDiffusion, "FILM_DIFFUSION", _disc.nParType, _disc.nComp, _unitOpIdx);
 	_parDiffusionMode = readAndRegisterMultiplexCompTypeSecParam(paramProvider, _parameters, _parDiffusion, "PAR_DIFFUSION", _disc.nParType, _disc.nComp, _unitOpIdx);
 
-	_disc.offsetSurfDiff = new unsigned int[_disc.strideBound[_disc.nParType]];
+	if (firstConfigCall)
+		_disc.offsetSurfDiff = new unsigned int[_disc.strideBound[_disc.nParType]];
 	if (paramProvider.exists("PAR_SURFDIFFUSION"))
 		_parSurfDiffusionMode = readAndRegisterMultiplexBndCompTypeSecParam(paramProvider, _parameters, _parSurfDiffusion, "PAR_SURFDIFFUSION", _disc.nParType, _disc.nComp, _disc.strideBound, _disc.nBound, _unitOpIdx);
 	else
@@ -775,7 +761,8 @@ bool GeneralRateModelDG::configure(IParameterProvider& paramProvider)
 		registerParam2DArray(_parameters, _parTypeVolFrac, [=](bool multi, unsigned cell, unsigned int type) { return makeParamId(hashString("PAR_TYPE_VOLFRAC"), _unitOpIdx, CompIndep, type, BoundStateIndep, ReactionIndep, cell); }, _disc.nParType);
 
 	// Calculate the particle radial discretization variables (_parCellSize, _parCenterRadius, etc.)
-	_disc.deltaR = new active[_disc.offsetMetric[_disc.nParType]];
+	if (firstConfigCall)
+		_disc.deltaR = new active[_disc.offsetMetric[_disc.nParType]];
 	updateRadialDisc();
 
 	// Register initial conditions parameters
@@ -2166,22 +2153,17 @@ void GeneralRateModelDG::updateRadialDisc()
 	{
 		for (int cell = 0; cell < _disc.nParCell[parType]; cell++)
 		{
-			_disc.Ir[_disc.offsetMetric[parType] + cell] = Vector<active, Dynamic>::Zero(_disc.nParNode[parType]);
-
 			for (int node = 0; node < _disc.nParNode[parType]; node++)
 				_disc.Ir[_disc.offsetMetric[parType] + cell][node] = _disc.deltaR[_disc.offsetMetric[parType] + cell] / 2.0 * (_disc.parNodes[parType][node] + 1.0);
 
-			_disc.Dr[_disc.offsetMetric[parType] + cell].resize(_disc.nParNode[parType], _disc.nParNode[parType]);
-			_disc.Dr[_disc.offsetMetric[parType] + cell].setZero();
-
 			active r_L = _parCoreRadius[parType] + cell * _disc.deltaR[_disc.offsetMetric[parType] + cell]; // left boundary of current cell
 
 			_disc.Ir[_disc.offsetMetric[parType] + cell] = _disc.Ir[_disc.offsetMetric[parType] + cell] + VectorXd::Ones(_disc.nParNode[parType]) * r_L;
 
 			if (_parGeomSurfToVol[parType] == SurfVolRatioSphere)
 				_disc.Ir[_disc.offsetMetric[parType] + cell] = _disc.Ir[_disc.offsetMetric[parType] + cell].array().square();
 			else if (_parGeomSurfToVol[parType] == SurfVolRatioSlab)
-				_disc.Ir[_disc.offsetMetric[parType] + cell] = Vector<active, Dynamic>::Ones(_disc.nParNode[parType]); // no metrics for slab
+				_disc.Ir[_disc.offsetMetric[parType] + cell].setOnes(); // no metric terms for slab
 
 			// (D_r)_{i, j} = D_{i, j} * (r_j / r_i) [only needed for inexact integration]
 			_disc.Dr[_disc.offsetMetric[parType] + cell] = _disc.parPolyDerM[parType];