Skip to content

Commit

Permalink
Add first numerical reference test
Browse files Browse the repository at this point in the history
  • Loading branch information
@jbreue16
jbreue16 committed Mar 6, 2024
1 parent eabf3f3 commit c1076bd
Show file tree
Hide file tree
Showing 5 changed files with 490 additions and 6 deletions.
163 changes: 160 additions & 3 deletions test/ColumnTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,9 @@
#include "JsonTestModels.hpp"
#include "JacobianHelper.hpp"
#include "UnitOperationTests.hpp"
#include "LoggingUtils.hpp"
#include "../include/io/hdf5/HDF5Reader.hpp"
#include "common/ParameterProviderImpl.hpp"

#include <cmath>
#include <functional>
Expand Down Expand Up @@ -175,7 +178,7 @@ namespace column
uint32_t _numElements;
};

void setNumAxialCells(cadet::JsonParameterProvider& jpp, unsigned int nCol)
void setNumAxialCells(cadet::JsonParameterProvider& jpp, unsigned int nCol, std::string unitID)
{
int level = 0;

Expand All @@ -184,9 +187,9 @@ namespace column
jpp.pushScope("model");
++level;
}
if (jpp.exists("unit_000"))
if (jpp.exists("unit_"+ unitID))
{
jpp.pushScope("unit_000");
jpp.pushScope("unit_" + unitID);
++level;
}

Expand All @@ -200,6 +203,31 @@ namespace column
jpp.popScope();
}

void setNumParCells(cadet::JsonParameterProvider& jpp, unsigned int nPar, std::string unitID)
{
int level = 0;

if (jpp.exists("model"))
{
jpp.pushScope("model");
++level;
}
if (jpp.exists("unit_" + unitID))
{
jpp.pushScope("unit_" + unitID);
++level;
}

jpp.pushScope("discretization");

jpp.set("NPAR", static_cast<int>(nPar));

jpp.popScope();

for (int l = 0; l < level; ++l)
jpp.popScope();
}

void setWenoOrder(cadet::JsonParameterProvider& jpp, int order)
{
int level = 0;
Expand Down Expand Up @@ -227,6 +255,55 @@ namespace column
jpp.popScope();
}

void setTimeSteppingResolution(nlohmann::json& setupJson)
{
nlohmann::json timeIntegrator;
timeIntegrator["ABSTOL"] = 1e-8;
timeIntegrator["ALGTOL"] = 1e-8;
timeIntegrator["INIT_STEP_SIZE"] = 1e-10;
timeIntegrator["MAX_STEPS"] = 1000000;
timeIntegrator["RELTOL"] = 1e-6;
setupJson["solver"]["time_integrator"] = timeIntegrator;
}

void copyNumericalMethod(cadet::ParameterProviderImpl<cadet::io::HDF5Reader>& pp, nlohmann::json& setupJson, const std::string unitID)
{
// copy over numerical methods from reference file. Note that we leave out spatial and time step resolution parameters
pp.pushScope("model");
pp.pushScope("solver");
nlohmann::json solver = setupJson["model"]["solver"];
solver["GS_TYPE"] = pp.getInt("GS_TYPE");
solver["MAX_KRYLOV"] = pp.getInt("MAX_KRYLOV");
solver["MAX_RESTARTS"] = pp.getInt("MAX_RESTARTS");
solver["SCHUR_SAFETY"] = pp.getDouble("SCHUR_SAFETY");
pp.popScope();
setupJson["model"]["solver"] = solver;

pp.pushScope("unit_" + unitID);
pp.pushScope("discretization");
nlohmann::json discretization = setupJson["model"]["unit_" + unitID]["discretization"];
discretization["NBOUND"] = pp.getIntArray("NBOUND");
discretization["RECONSTRUCTION"] = pp.getInt("RECONSTRUCTION");
discretization["USE_ANALYTIC_JACOBIAN"] = pp.getInt("USE_ANALYTIC_JACOBIAN");
pp.pushScope("weno");
nlohmann::json weno;
weno["WENO_ORDER"] = pp.getInt("WENO_ORDER");
weno["WENO_EPS"] = pp.getDouble("WENO_EPS");
weno["BOUNDARY_MODEL"] = pp.getInt("BOUNDARY_MODEL");
discretization["weno"] = weno;
pp.popScope();
setupJson["model"]["unit_" + unitID]["discretization"] = discretization;
pp.popScope();
pp.popScope();
pp.popScope();

pp.pushScope("solver");
setupJson["solver"]["CONSISTENT_INIT_MODE"] = pp.getInt("CONSISTENT_INIT_MODE");
setupJson["solver"]["CONSISTENT_INIT_MODE_SENS"] = pp.getInt("CONSISTENT_INIT_MODE_SENS");
setupJson["solver"]["NTHREADS"] = pp.getInt("NTHREADS");
pp.popScope();
}

void reverseFlow(cadet::JsonParameterProvider& jpp)
{
int level = 0;
Expand Down Expand Up @@ -1139,6 +1216,86 @@ namespace column
destroyModelBuilder(mb);
}

void testReferenceBenchmark(const std::string& setupFileRelPath, const std::string& refFileRelPath, const std::string& unitID, const std::vector<double> absTol, const std::vector<double> relTol, const unsigned int nCol, const unsigned int nPar, const bool compare_sens)
{

// read json model setup file
const std::string setupFile = std::string(getTestDirectory()) + setupFileRelPath;
JsonParameterProvider pp_setup(JsonParameterProvider::fromFile(setupFile));

// adjust numerical parameters
nlohmann::json* setupJson = pp_setup.data();
// copy over some numerical parameters from reference file, e.g. consistent initialization
cadet::io::HDF5Reader rd;
const std::string refFile = std::string(getTestDirectory()) + refFileRelPath;
rd.openFile(refFile, "r");
ParameterProviderImpl<cadet::io::HDF5Reader> pp_ref(rd);
copyNumericalMethod(pp_ref, *setupJson, unitID);
rd.closeFile();

// set remaining numerical parameters
setTimeSteppingResolution(*setupJson);
setNumAxialCells(pp_setup, nCol, unitID);
if (nPar > 0)
setNumParCells(pp_setup, nPar, unitID);

// run simulation
Driver drv;
drv.configure(pp_setup);
drv.run();

// get simulation result
InternalStorageUnitOpRecorder const* const simData = drv.solution()->unitOperation(1);
double const* sim_outlet = simData->outlet();

// read h5 reference data
rd.openFile(refFile, "r");

// get outlet and sensitivity reference
pp_ref.popScope();
pp_ref.pushScope("output");
pp_ref.pushScope("solution");
pp_ref.pushScope("unit_" + unitID);
const std::vector<double> ref_outlet = pp_ref.getDoubleArray("SOLUTION_OUTLET");
pp_ref.popScope();
pp_ref.popScope();

// compare the simulation results with the reference data

for (unsigned int i = 0; i < ref_outlet.size(); ++i)
CHECK((sim_outlet[i]) == cadet::test::makeApprox(ref_outlet[i], relTol[0], absTol[0]));


if (pp_ref.exists("sensitivity") && compare_sens)
{
pp_ref.pushScope("sensitivity");

unsigned int sensID = 0;
std::string sensParam = std::to_string(sensID);
sensParam = "param_" + std::string(3 - sensParam.length(), '0') + sensParam;

while (pp_ref.exists(sensParam))
{
pp_ref.pushScope(sensParam);
pp_ref.pushScope("unit_" + unitID);
const std::vector<double> ref_sens = pp_ref.getDoubleArray("SENS_OUTLET");
pp_ref.popScope();
pp_ref.popScope();

double const* sim_sens = simData->sensOutlet(sensID);

for (unsigned int i = 0; i < ref_sens.size(); ++i)
CHECK((sim_sens[i]) == cadet::test::makeApprox(ref_sens[i], relTol[sensID + 1], absTol[sensID + 1]));

sensID++;
sensParam = std::to_string(sensID);
sensParam = "param_" + std::string(3 - sensParam.length(), '0') + sensParam;
}
}

rd.closeFile();
}

} // namespace column
} // namespace test
} // namespace cadet
15 changes: 14 additions & 1 deletion test/ColumnTests.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,8 +38,9 @@ namespace column
* @details Overwrites the NCOL field in the discretization group of the given ParameterProvider.
* @param [in,out] jpp ParameterProvider to change the number of axial cells in
* @param [in] nCol Number of axial cells
* @param [in] unitID unit operation ID
*/
void setNumAxialCells(cadet::JsonParameterProvider& jpp, unsigned int nCol);
void setNumAxialCells(cadet::JsonParameterProvider& jpp, unsigned int nCol, std::string unitID="000");

/**
* @brief Sets the WENO order in a configuration of a column-like unit operation
Expand Down Expand Up @@ -279,6 +280,18 @@ namespace column
*/
void testInletDofJacobian(const std::string& uoType);

/**
* @brief Runs a simulation test comparing against numerical reference data (outlet data)
* @param [in] setupFileRelPath Path to the setup data file from the directory of this file
* @param [in] refFileRelPath Path to the reference data file from the directory of this file
* @param [in] unitID ID of the unit of interest
* @param [in] absTol Absolute error tolerance
* @param [in] relTol Relative error tolerance
* @param [in] nCol Number of axial cells
* @param [in] nPar Number of particle cells (if required)
*/
void testReferenceBenchmark(const std::string& setupFileRelPath, const std::string& refFileRelPath, const std::string& unitID, const std::vector<double> absTol, const std::vector<double> relTol, const unsigned int nCol, const unsigned int nPar=0, const bool compare_sens=false);

} // namespace column
} // namespace test
} // namespace cadet
Expand Down
13 changes: 11 additions & 2 deletions test/LumpedRateModelWithoutPores.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,15 +31,15 @@ TEST_CASE("LRM LWE forward vs backward flow", "[LRM],[Simulation],[CIlrm]")
cadet::test::column::testWenoForwardBackward("LUMPED_RATE_MODEL_WITHOUT_PORES", i, 6e-9, 6e-4);
}

TEST_CASE("LRM linear pulse vs analytic solution", "[LRM],[Simulation],[Analytic],[CIlrm]")
TEST_CASE("LRM linear pulse vs analytic solution", "[LRM],[Simulation],[Reference],[Analytic],[CIlrm]")
{
cadet::test::column::testAnalyticBenchmark("LUMPED_RATE_MODEL_WITHOUT_PORES", "/data/lrm-pulseBenchmark.data", true, true, 1024, 2e-5, 1e-7);
cadet::test::column::testAnalyticBenchmark("LUMPED_RATE_MODEL_WITHOUT_PORES", "/data/lrm-pulseBenchmark.data", true, false, 1024, 2e-5, 1e-7);
cadet::test::column::testAnalyticBenchmark("LUMPED_RATE_MODEL_WITHOUT_PORES", "/data/lrm-pulseBenchmark.data", false, true, 1024, 2e-5, 1e-7);
cadet::test::column::testAnalyticBenchmark("LUMPED_RATE_MODEL_WITHOUT_PORES", "/data/lrm-pulseBenchmark.data", false, false, 1024, 2e-5, 1e-7);
}

TEST_CASE("LRM non-binding linear pulse vs analytic solution", "[LRM],[Simulation],[Analytic],[NonBinding],[CIlrm]")
TEST_CASE("LRM non-binding linear pulse vs analytic solution", "[LRM],[Simulation],[Reference],[Analytic],[NonBinding],[CIlrm]")
{
cadet::test::column::testAnalyticNonBindingBenchmark("LUMPED_RATE_MODEL_WITHOUT_PORES", "/data/lrm-nonBinding.data", true, 1024, 2e-5, 1e-7);
cadet::test::column::testAnalyticNonBindingBenchmark("LUMPED_RATE_MODEL_WITHOUT_PORES", "/data/lrm-nonBinding.data", false, 1024, 2e-5, 1e-7);
Expand All @@ -52,6 +52,15 @@ TEST_CASE("LRM Jacobian forward vs backward flow", "[LRM],[UnitOp],[Residual],[J
cadet::test::column::testJacobianWenoForwardBackward("LUMPED_RATE_MODEL_WITHOUT_PORES", i);
}

TEST_CASE("LRM numerical Benchmark with parameter sensitivities for linear case", "[LRM],[Simulation],[Reference],[Sensitivity],[testHereLRM1]")
{
const std::string& setup = std::string("/data/setup_LRM_dynLin_1comp.json");
const std::string& reference = std::string("/data/ref_LRM_dynLin_1comp.h5");
const std::vector<double> absTol = { 1e-8, 1e-8, 1e-8, 1e-8 };
const std::vector<double> relTol = { 1.0e-3, 1.0e-1, 1.0e-3, 1.0e-3 };
cadet::test::column::testReferenceBenchmark(setup, reference, "001", absTol, relTol, 256, 0, true);
}

TEST_CASE("LRM time derivative Jacobian vs FD", "[LRM],[UnitOp],[Residual],[Jacobian],[FDtestLRM]")
{
cadet::test::column::testTimeDerivativeJacobianFD("LUMPED_RATE_MODEL_WITHOUT_PORES");
Expand Down
Binary file added test/data/ref_LRM_dynLin_1comp.h5
View file
Binary file not shown.
Loading

0 comments on commit c1076bd

Please sign in to comment.