Skip to content

Commit

Permalink
Add tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@AntoniaBerger
AntoniaBerger committed Oct 23, 2024
1 parent e85e10f commit c38a49b
Show file tree
Hide file tree
Showing 2 changed files with 39 additions and 24 deletions.
32 changes: 13 additions & 19 deletions src/libcadet/model/exchange/LangumirExchange.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ class LangumirExchangeBase : public IPhaseTransitionModel
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;
Expand All @@ -68,7 +68,7 @@ class LangumirExchangeBase : public IPhaseTransitionModel
_parameters.clear();
readParameterMatrix(_exchangeMatrix, paramProvider, "EXCHANGE_MATRIX", _nChannel * _nChannel * _nComp, 1); // include parameterPeaderHelp in exchange modul
_crossSections = paramProvider.getDoubleArray("CHANNEL_CROSS_SECTION_AREAS");
readParameterMatrix(_capacityMatrix, paramProvider, "CAPACITY_MATRIX", _nComp * _nChannel, 1); // include parameterPeaderHelp in exchange modul
readParameterMatrix(_capacityMatrix, paramProvider, "CAPACITY_MATRIX", _nComp * _nChannel, 1);

return true;
}
Expand Down Expand Up @@ -228,7 +228,8 @@ class LangumirExchangeBase : public IPhaseTransitionModel
if (cadet_likely(exchange_orig_dest_comp > 0.0))
{
double compSum_orig_dest = 0.0;
for (unsigned int component = 0; component < _nComp; component++) {
for (unsigned int component = 0; component < _nComp; component++)
{

StateType const* const compSum = yColRadDestBlock + component;
double compSum_dest = static_cast<double>(compSum[0]);
Expand All @@ -241,24 +242,17 @@ class LangumirExchangeBase : public IPhaseTransitionModel
}

compSum_orig_dest += compSum_dest / cMax_destComp;

//double chanSum_orig = 0.0;
//double chanSum_dest = 0.0;
//for (unsigned int channel = 0; channel < _nChannel; channel++) {
//
// chanSum_orig += static_cast<double>(yColBlock[component + channel * _nChannel]) / cMax_origComp;
//
// chanSum_dest += static_cast<double>(yColBlock[component + channel * _nChannel]) / cMax_destComp;
//}

}

*resCur_orig += exchange_orig_dest_comp * yCur_orig[0] * capacity_orig_comp * (1 - compSum_orig_dest);
if (capacity_orig_comp < 1e-16) {
if (capacity_orig_comp < 1e-16)
{
*resCur_orig += exchange_orig_dest_comp * yCur_orig[0];
}

*resCur_dest -= exchange_orig_dest_comp * yCur_orig[0] * static_cast<ParamType>(_crossSections[rad_orig]) / static_cast<ParamType>(_crossSections[rad_dest]) * capacity_orig_comp * (1 - compSum_orig_dest);
if (capacity_orig_comp < 1e-16) {
if (capacity_orig_comp < 1e-16)
{
*resCur_dest -= exchange_orig_dest_comp * yCur_orig[0] * static_cast<ParamType>(_crossSections[rad_orig]) / static_cast<ParamType>(_crossSections[rad_dest]);
}

Expand All @@ -268,22 +262,22 @@ class LangumirExchangeBase : public IPhaseTransitionModel
linalg::BandedSparseRowIterator jacorig;
jacorig = jacBegin + offsetCur_orig;
jacorig[0] += static_cast<double>(exchange_orig_dest_comp) * static_cast<double>(capacity_orig_comp * (1 - compSum_orig_dest));
if (capacity_orig_comp < 1e-16) {
if (capacity_orig_comp < 1e-16)
{
jacorig[0] -= static_cast<double>(exchange_orig_dest_comp);
}

linalg::BandedSparseRowIterator jacdest;
jacdest = jacBegin + offsetCur_dest;
jacdest[static_cast<int>(offsetCur_orig) - static_cast<int>(offsetCur_dest)] -= static_cast<double>(exchange_orig_dest_comp) * static_cast<double>(capacity_orig_comp) * (1 - compSum_orig_dest);
if (capacity_orig_comp < 1e-16) {
if (capacity_orig_comp < 1e-16)
{
jacdest[static_cast<int>(offsetCur_orig) - static_cast<int>(offsetCur_dest)] += static_cast<double>(exchange_orig_dest_comp);
}
}


}


}

}
Expand Down
31 changes: 26 additions & 5 deletions test/MultiChannelTransportModel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -469,7 +469,7 @@ TEST_CASE("MCT numerical Benchmark comparison with LRM (1 channel no exchange, n
cadet::test::column::testReferenceBenchmark(modelFilePath, refFilePath, "001", absTol, relTol, disc, false);
}

TEST_CASE("MCT numerical Benchmark comparison with linear binding LRM (2 channel with exchange, no reaction case)", "[MCT],[Simulation],[Reference],[mctReference],[CI]")
TEST_CASE("MCT numerical Benchmark comparison with linear binding LRM (2 channel with linear exchange, no reaction case)", "[MCT],[Simulation],[Reference],[mctReference],[CI]")
{
const std::string& modelFilePath = std::string("/data/model_MCT2ch_1comp_benchmark1.json");
const std::string& refFilePath = std::string("/data/ref_LRM_dynLin_1comp_benchmark2_FV_Z357.h5");
Expand All @@ -492,7 +492,7 @@ TEST_CASE("MCT numerical Benchmark for 1 channel no exchange, with reaction case
cadet::test::column::testReferenceBenchmark(modelFilePath, refFilePath, "001", absTol, relTol, disc, false);
}

TEST_CASE("MCT numerical Benchmark for 2 channels with one-way-exchange and reaction case", "[MCT],[Simulation],[Reference],[mctReference]") // todo CI flag: currently only runs locally but fails on server
TEST_CASE("MCT numerical Benchmark for 2 channels with one-way-linear exchange and reaction case", "[MCT],[Simulation],[Reference],[mctReference]") // todo CI flag: currently only runs locally but fails on server
{
const std::string& modelFilePath = std::string("/data/model_MCT2ch_oneWayEx_reac_benchmark1.json");
const std::string& refFilePath = std::string("/data/ref_MCT2ch_oneWayEx_reac_benchmark1_FV_Z256.h5");
Expand All @@ -503,7 +503,7 @@ TEST_CASE("MCT numerical Benchmark for 2 channels with one-way-exchange and reac
cadet::test::column::testReferenceBenchmark(modelFilePath, refFilePath, "001", absTol, relTol, disc, false);
}

TEST_CASE("MCT numerical Benchmark for 3 channels with two-way-exchange and reaction case", "[MCT],[Simulation],[Reference],[mctReference]") // todo CI flag: currently only runs locally but fails on server
TEST_CASE("MCT numerical Benchmark for 3 channels with two-way-linear exchange and reaction case", "[MCT],[Simulation],[Reference],[mctReference]") // todo CI flag: currently only runs locally but fails on server
{
const std::string& modelFilePath = std::string("/data/model_MCT3ch_twoWayExc_reac_benchmark1.json");
const std::string& refFilePath = std::string("/data/ref_MCT3ch_twoWayExc_reac_benchmark1_FV_Z256.h5");
Expand All @@ -522,7 +522,7 @@ TEST_CASE("MCT compare AD with analytical Jacobian for 1 channel without exchang
cadet::test::column::testJacobianAD(jpp, std::numeric_limits<float>::epsilon());
}

TEST_CASE("MCT compare AD with analytical Jacobian for 2 channels and exchange", "[MCT],[UnitOp],[Jacobian],[CI]")
TEST_CASE("MCT compare AD with analytical Jacobian for 2 channels and linear exchange", "[MCT],[UnitOp],[Jacobian],[CI]")
{
cadet::JsonParameterProvider jpp = createMCT({ 1.0, 1.0 }, { 1.0, 1.0 }, { 1.0, 0.2 }, { 1.0, 1.0 }, { 0.0, 0.01, 0.0, 0.0 }, 1e-4); // increased col dispersion so that jacobian entries are above tolerances
jpp.pushScope("model");
Expand All @@ -531,7 +531,7 @@ TEST_CASE("MCT compare AD with analytical Jacobian for 2 channels and exchange",
cadet::test::column::testJacobianAD(jpp, FDtolerance);
}

TEST_CASE("MCT compare AD with analytical Jacobian for 2 channels with opposing flow directions and exchange", "[MCT],[UnitOp],[Jacobian],[CI]")
TEST_CASE("MCT compare AD with analytical Jacobian for 2 channels with opposing flow directions and linear exchange", "[MCT],[UnitOp],[Jacobian],[CI]")
{
cadet::JsonParameterProvider jpp = createMCT({ 1.0, 1.0 }, { 1.0, -1.0 }, { 1.0, 0.2 }, { 1.0, 1.0 }, { 0.0, 0.01, 0.0, 0.0 }, 1e-4); // increased col dispersion so that jacobian entries are above tolerances
jpp.pushScope("model");
Expand Down Expand Up @@ -604,4 +604,25 @@ TEST_CASE("MCT dynamic reactions time derivative Jacobian vs FD modified bulk",
TEST_CASE("MCT dynamic reactions Jacobian vs AD bulk", "[MCT],[Jacobian],[AD],[ReactionModel],[CI]")
{
cadet::test::reaction::testUnitJacobianDynamicReactionsAD("MULTI_CHANNEL_TRANSPORT", "FV", true, false, false, std::numeric_limits<float>::epsilon() * 100.0);
}

TEST_CASE("MCT compare AD with analytical Jacobian for 2 channels and langmuir exchange", "[testHere]")
{
cadet::JsonParameterProvider jpp = createMCT({ 1.0, 1.0 }, { 1.0, 1.0 }, { 1.0, 0.2 }, { 1.0, 1.0 }, { 0.0, 0.01, 0.0, 0.0 }, 1e-4); // increased col dispersion so that jacobian entries are above tolerances
jpp.pushScope("model");
jpp.pushScope("unit_000");
const double FDtolerance = 0.02; // large tolerance to effectively disable FD pattern check, which fails with 0.0 != -0.01
cadet::test::column::testJacobianAD(jpp, FDtolerance);
}

TEST_CASE("MCT numerical Benchmark comparison with linear binding LRM (2 channel with langmuir exchange, no reaction case)", "[testHere]")
{
const std::string& modelFilePath = std::string("/data/model_MCT2ch_1comp_benchmark1.json");
const std::string& refFilePath = std::string("/data/ref_LRM_dynLin_1comp_benchmark2_FV_Z357.h5");
const std::vector<double> absTol = { RelApprox::defaultEpsilon() };
const std::vector<double> relTol = { RelApprox::defaultMargin() };
cadet::test::column::FVparams disc(357);
disc.setNRad(2); // will be used as NCHANNEL

cadet::test::column::testReferenceBenchmark(modelFilePath, refFilePath, "001", absTol, relTol, disc, false, 2);
}

0 comments on commit c38a49b

Please sign in to comment.