Skip to content

Commit

Permalink
Add DG unit tests (#157)
Browse files Browse the repository at this point in the history 
Add tolerance option to AD jacobian comparison tests

Add relative tolerances to FV numerical reference tests
  • Loading branch information
@jbreue16
jbreue16 committed Jun 18, 2024
1 parent 4d044c8 commit 5a696d5
Show file tree
Hide file tree
Showing 30 changed files with 1,257 additions and 227 deletions.
51 changes: 48 additions & 3 deletions src/libcadet/AdUtils.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,10 +100,10 @@ void extractBandedEigenJacobianFromAd(active const* const adVec, int adDirOffset
}

void extractBandedBlockEigenJacobianFromAd(active const* const adVec, int adDirOffset, int diagDir, const int lowerBandwidth, const int upperBandwidth,
const int blockOffset, const int nCols, Eigen::SparseMatrix<double, Eigen::RowMajor>& mat, const int matrixOffset)
const int blockOffset, const int nRows, Eigen::SparseMatrix<double, Eigen::RowMajor>& mat, const int matrixOffset)
{
const int stride = lowerBandwidth + 1 + upperBandwidth;
for (int eq = blockOffset; eq < blockOffset + nCols; ++eq)
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
Expand All @@ -114,7 +114,7 @@ void extractBandedBlockEigenJacobianFromAd(active const* const adVec, int adDirO
for (int diag = 0; diag < stride; ++diag)
{
if (eq - lowerBandwidth + diag >= blockOffset && // left block boundary
eq - lowerBandwidth + diag < blockOffset + nCols && // right 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);
Expand Down Expand Up @@ -281,6 +281,51 @@ 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<double, Eigen::RowMajor>& mat, const int matrixOffset)
{
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<double>::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;
}

void adMatrixVectorMultiply(const linalg::SparseMatrix<active>& mat, double const* x, double* y, double alpha, double beta, int adDir)
{
const std::vector<int>& rows = mat.rows();
Expand Down
20 changes: 20 additions & 0 deletions src/libcadet/AdUtils.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -145,6 +145,26 @@ 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<double, Eigen::RowMajor>& 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
Expand Down
67 changes: 50 additions & 17 deletions src/libcadet/model/LumpedRateModelWithPoresDG.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -805,23 +805,56 @@ void LumpedRateModelWithPoresDG::extractJacobianFromAD(active const* const adRes
*/
void LumpedRateModelWithPoresDG::checkAnalyticJacobianAgainstAd(active const* const adRes, unsigned int adDirOffset) const
{
// todo write this function
//Indexer idxr(_disc);

//LOG(Debug) << "AD dir offset: " << adDirOffset << " DiagDirCol: " << _convDispOp.jacobian().lowerBandwidth();

//// Column
//const double maxDiffCol = _convDispOp.checkAnalyticJacobianAgainstAd(adRes, adDirOffset);

//// Particles
//double maxDiffPar = 0.0;
//for (unsigned int type = 0; type < _disc.nParType; ++type)
//{
// const linalg::BandMatrix& jacMat = _jacP[type];
// const double localDiff = ad::compareBandedJacobianWithAd(adRes + idxr.offsetCp(ParticleTypeIndex{ type }), adDirOffset, jacMat.lowerBandwidth(), jacMat);
// LOG(Debug) << "-> Par type " << type << " diff: " << localDiff;
// maxDiffPar = std::max(maxDiffPar, localDiff);
//}
Indexer idxr(_disc);

const int lowerBandwidth = (_disc.exactInt) ? 2 * _disc.nNodes * idxr.strideColNode() : _disc.nNodes * idxr.strideColNode();
const int upperBandwidth = lowerBandwidth;
const int stride = lowerBandwidth + 1 + upperBandwidth;

const double maxDiff = ad::compareBandedEigenJacobianWithAd(adRes + idxr.offsetC(), adDirOffset, lowerBandwidth, lowerBandwidth, upperBandwidth, 0, numPureDofs(), _globalJac, 0);

if (maxDiff > 1e-6)
int jojo = 0;

double maxDiffPar = 0.0;

for (unsigned int type = 0; type < _disc.nParType; type++)
{
int offsetParticleTypeDirs = adDirOffset + idxr.offsetCp(ParticleTypeIndex{ type }) - idxr.offsetC();

for (unsigned int par = 0; par < _disc.nPoints; par++)
{
const int eqOffset_res = idxr.offsetCp(ParticleTypeIndex{ type }, ParticleIndex{ par });
const int eqOffset_mat = idxr.offsetCp(ParticleTypeIndex{ type }, ParticleIndex{ par }) - idxr.offsetC();
for (unsigned int phase = 0; phase < idxr.strideParBlock(type); phase++)
{
for (unsigned int phaseTo = 0; phaseTo < idxr.strideParBlock(type); phaseTo++)
{

double baseVal = adRes[eqOffset_res + phase].getADValue(offsetParticleTypeDirs + phaseTo);
double matVal = _globalJac.coeff(eqOffset_mat + phase, eqOffset_mat + phaseTo);

if (std::isnan(matVal) || std::isnan(baseVal))
continue;

const double diff = std::abs(matVal - baseVal);

baseVal = std::abs(baseVal);
if (baseVal > 0.0)
maxDiffPar = std::max(maxDiffPar, diff / baseVal);
else
maxDiffPar = std::max(maxDiffPar, diff);

}
}
}
offsetParticleTypeDirs += idxr.strideParBlock(type);
}

if (maxDiffPar > 1e-6)
int jojo2 = 0;

LOG(Debug) << "AD dir offset: " << adDirOffset << " DiagBlockSize: " << stride << " MaxDiffBulk: " << maxDiff << " MaxDiffParticle: " << maxDiffPar;
}

#endif
Expand Down
20 changes: 13 additions & 7 deletions src/libcadet/model/LumpedRateModelWithoutPoresDG.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -387,8 +387,7 @@ namespace cadet
#ifndef CADET_CHECK_ANALYTIC_JACOBIAN
return _jacobianAdDirs;
#else
// If CADET_CHECK_ANALYTIC_JACOBIAN is active, we always need the AD directions for the Jacobian
return _jac.cols();// _jac.stride();@SAM?
return _convDispOp.requiredADdirs();
#endif
}

Expand Down Expand Up @@ -427,10 +426,10 @@ namespace cadet
const int stride = lowerBandwidth + 1 + upperBandwidth;

int diagDir = lowerBandwidth;
const int nCols = _jac.cols();
const int nRows = _jac.rows();
const int colOffset = 0;

ad::extractBandedBlockEigenJacobianFromAd(adVec, adDirOffset, diagDir, lowerBandwidth, upperBandwidth, colOffset, nCols, _jac);
ad::extractBandedBlockEigenJacobianFromAd(adVec, adDirOffset, diagDir, lowerBandwidth, upperBandwidth, colOffset, nRows, _jac);

}

Expand All @@ -446,9 +445,16 @@ namespace cadet
{
Indexer idxr(_disc);

// @TODO
//const double maxDiff = ad::compareBandedJacobianWithAd(adRes + idxr.offsetC(), adDirOffset, _jac.lowerBandwidth(), _jac);
//LOG(Debug) << "AD dir offset: " << adDirOffset << " DiagDirCol: " << _jac.lowerBandwidth() << " MaxDiff: " << maxDiff;
const int lowerBandwidth = (_disc.exactInt) ? 2 * _disc.nNodes * idxr.strideColNode() : _disc.nNodes * idxr.strideColNode();
const int upperBandwidth = lowerBandwidth;
const int stride = lowerBandwidth + 1 + upperBandwidth;

const double maxDiff = ad::compareBandedEigenJacobianWithAd(adRes + idxr.offsetC(), adDirOffset, lowerBandwidth, lowerBandwidth, upperBandwidth, 0, _jac.rows(), _jac, 0);

if (maxDiff > 1e-6)
int jojo = 0;

LOG(Debug) << "AD dir offset: " << adDirOffset << " DiagBlockSize: " << stride << " MaxDiff: " << maxDiff;
}

#endif
Expand Down
3 changes: 2 additions & 1 deletion test/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,8 @@ if (ENABLE_GRM_2D)
endif()

add_executable(testRunner testRunner.cpp JsonTestModels.cpp ColumnTests.cpp UnitOperationTests.cpp SimHelper.cpp ParticleHelper.cpp
GeneralRateModel.cpp GeneralRateModel2D.cpp LumpedRateModelWithPores.cpp LumpedRateModelWithoutPores.cpp LumpedRateModelWithoutPoresDG.cpp
GeneralRateModel.cpp GeneralRateModel2D.cpp LumpedRateModelWithPores.cpp LumpedRateModelWithoutPores.cpp
GeneralRateModelDG.cpp LumpedRateModelWithPoresDG.cpp LumpedRateModelWithoutPoresDG.cpp
RadialGeneralRateModel.cpp RadialLumpedRateModelWithPores.cpp RadialLumpedRateModelWithoutPores.cpp
CSTR-Residual.cpp CSTR-Simulation.cpp
ConvectionDispersionOperator.cpp
Expand Down
Loading

0 comments on commit 5a696d5

Please sign in to comment.