Skip to content

Commit

Permalink
Fix multiple particle type bug
Browse files Browse the repository at this point in the history
  • Loading branch information
@jbreue16
jbreue16 committed Oct 20, 2024
1 parent 280170f commit 9ac4df6
Show file tree
Hide file tree
Showing 3 changed files with 45 additions and 42 deletions.
7 changes: 1 addition & 6 deletions src/libcadet/model/GeneralRateModelDG.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2153,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];
Expand Down
12 changes: 11 additions & 1 deletion src/libcadet/model/GeneralRateModelDG.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -405,13 +405,23 @@ class GeneralRateModelDG : public UnitOperationBase

offsetMetric = VectorXi::Zero(nParType + 1);
for (int parType = 1; parType <= nParType; parType++) {
offsetMetric[parType] += nParCell[parType - 1];
offsetMetric[parType] = offsetMetric[parType - 1] + nParCell[parType - 1];
}

if (firstConfigCall)
{
Dr = new MatrixXd[offsetMetric[nParType]];
Ir = new Vector<active, Dynamic>[offsetMetric[nParType]];
for (int parType = 0; parType < nParType; parType++)
{
for (int cell = 0; cell < nParCell[parType]; cell++)
{
Ir[offsetMetric[parType] + cell].resize(nParNode[parType]);
Ir[offsetMetric[parType] + cell].setZero();
Dr[offsetMetric[parType] + cell].resize(nParNode[parType], nParNode[parType]);
Dr[offsetMetric[parType] + cell].setZero();
}
}
minus_InvMM_ST = new MatrixXd[offsetMetric[nParType]];
parInvMM = new MatrixXd[offsetMetric[nParType]];
secondOrderStiffnessM = new MatrixXd[nParType];
Expand Down
68 changes: 33 additions & 35 deletions test/GeneralRateModelDG.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,42 +248,40 @@ TEST_CASE("GRM_DG with two component linear binding Jacobian", "[GRM],[DG],[Unit
cadet::test::column::testJacobianAD(jpp, std::numeric_limits<float>::epsilon() * 100.0);
}

//// todo fix cant load library with debug info
//TEST_CASE("GRM_DG LWE one vs two identical particle types match", "[GRM],[DG],[Simulation],[ParticleType],[todo]")
//{
// cadet::test::particle::testOneVsTwoIdenticalParticleTypes("GENERAL_RATE_MODEL", "DG", 2e-8, 5e-5);
//}
TEST_CASE("GRM_DG LWE one vs two identical particle types match", "[GRM],[DG],[Simulation],[ParticleType],[CI]")
{
cadet::test::particle::testOneVsTwoIdenticalParticleTypes("GENERAL_RATE_MODEL", "DG", 2e-8, 5e-5);
}

//// todo fix cant load library with debug info in all of the following tests
//TEST_CASE("GRM_DG LWE separate identical particle types match", "[GRM],[DG],[Simulation],[ParticleType],[todo]")
//{
// cadet::test::particle::testSeparateIdenticalParticleTypes("GENERAL_RATE_MODEL", "DG", 1e-15, 1e-15);
//}
//
//TEST_CASE("GRM_DG linear binding single particle matches particle distribution", "[GRM],[DG],[Simulation],[ParticleType],[todo]")
//{
// cadet::test::particle::testLinearMixedParticleTypes("GENERAL_RATE_MODEL", "DG", 5e-8, 5e-5);
//}
//
//TEST_CASE("GRM_DG multiple particle types Jacobian analytic vs AD", "[GRM],[DG],[Jacobian],[AD],[ParticleType],[todo]")
//{
// cadet::test::particle::testJacobianMixedParticleTypes("GENERAL_RATE_MODEL", "DG");
//}
//
//TEST_CASE("GRM_DG multiple particle types time derivative Jacobian vs FD", "[GRM],[DG],[UnitOp],[Residual],[Jacobian],[ParticleType],[todo]")
//{
// cadet::test::particle::testTimeDerivativeJacobianMixedParticleTypesFD("GENERAL_RATE_MODEL", "DG", 1e-6, 0.0, 9e-4);
//}
//
//TEST_CASE("GRM_DG multiple spatially dependent particle types Jacobian analytic vs AD", "[GRM],[DG],[Jacobian],[AD],[ParticleType],[todo]")
//{
// cadet::test::particle::testJacobianSpatiallyMixedParticleTypes("GENERAL_RATE_MODEL", "DG");
//}
//
//TEST_CASE("GRM_DG linear binding single particle matches spatially dependent particle distribution", "[GRM],[DG],[Simulation],[ParticleType],[todo]")
//{
// cadet::test::particle::testLinearSpatiallyMixedParticleTypes("GENERAL_RATE_MODEL", "DG", 5e-8, 5e-5);
//}
TEST_CASE("GRM_DG LWE separate identical particle types match", "[GRM],[DG],[Simulation],[ParticleType],[CI]")
{
cadet::test::particle::testSeparateIdenticalParticleTypes("GENERAL_RATE_MODEL", "DG", 2e-8, 5e-5);
}

TEST_CASE("GRM_DG linear binding single particle matches particle distribution", "[GRM],[DG],[Simulation],[ParticleType],[CI]")
{
cadet::test::particle::testLinearMixedParticleTypes("GENERAL_RATE_MODEL", "DG", 5e-8, 5e-5);
}

TEST_CASE("GRM_DG multiple particle types Jacobian analytic vs AD", "[GRM],[DG],[Jacobian],[AD],[ParticleType],[CI]")
{
cadet::test::particle::testJacobianMixedParticleTypes("GENERAL_RATE_MODEL", "DG");
}

TEST_CASE("GRM_DG multiple particle types time derivative Jacobian vs FD", "[GRM],[DG],[UnitOp],[Residual],[Jacobian],[ParticleType],[CI]")
{
cadet::test::particle::testTimeDerivativeJacobianMixedParticleTypesFD("GENERAL_RATE_MODEL", "DG", 1e-6, 0.0, 9e-4);
}

TEST_CASE("GRM_DG multiple spatially dependent particle types Jacobian analytic vs AD", "[GRM],[DG],[Jacobian],[AD],[ParticleType],[CI]")
{
cadet::test::particle::testJacobianSpatiallyMixedParticleTypes("GENERAL_RATE_MODEL", "DG", 1e-11);
}

TEST_CASE("GRM_DG linear binding single particle matches spatially dependent particle distribution", "[GRM],[DG],[Simulation],[ParticleType],[CI]")
{
cadet::test::particle::testLinearSpatiallyMixedParticleTypes("GENERAL_RATE_MODEL", "DG", 5e-8, 5e-5);
}

TEST_CASE("GRM_DG dynamic reactions Jacobian vs AD bulk", "[GRM],[DG],[Jacobian],[AD],[ReactionModel],[CI]")
{
Expand Down

0 comments on commit 9ac4df6

Please sign in to comment.