Commit 54c8a262 authored by jannick.wolters's avatar jannick.wolters
Browse files

added checkerboard test case to mn and pn

parent f449aac2
Pipeline #112310 failed with stage
in 12 minutes and 52 seconds
......@@ -65,10 +65,10 @@ inline std::map<std::string, SOLVER_NAME> Solver_Map{
{ "SN_SOLVER", SN_SOLVER }, { "CSD_SN_SOLVER", CSD_SN_SOLVER }, { "PN_SOLVER", PN_SOLVER }, { "MN_SOLVER", MN_SOLVER } };
// Entropy functional
enum ENTROPY_NAME { QUADRATIC, MAXWELL_BOLZMANN, BOSE_EINSTEIN, FERMI_DIRAC };
enum ENTROPY_NAME { QUADRATIC, MAXWELL_BOLTZMANN, BOSE_EINSTEIN, FERMI_DIRAC };
inline std::map<std::string, ENTROPY_NAME> Entropy_Map{
{ "QUADRATIC", QUADRATIC }, { "MAXWELL_BOLZMANN", MAXWELL_BOLZMANN }, { "BOSE_EINSTEIN", BOSE_EINSTEIN }, { "FERMI_DIRAC", FERMI_DIRAC } };
{ "QUADRATIC", QUADRATIC }, { "MAXWELL_BOLTZMANN", MAXWELL_BOLTZMANN }, { "BOSE_EINSTEIN", BOSE_EINSTEIN }, { "FERMI_DIRAC", FERMI_DIRAC } };
// Optimizer
enum OPTIMIZER_NAME { NEWTON, ML };
......
......@@ -3,20 +3,42 @@
#include "problembase.h"
class Checkerboard : public ProblemBase
class Checkerboard_SN : public ProblemBase
{
private:
Vector _scatteringXS;
Vector _totalXS;
Checkerboard() = delete;
Checkerboard_SN() = delete;
bool isAbsorption( const Vector& pos ) const;
bool isSource( const Vector& pos ) const;
public:
Checkerboard( Config* settings, Mesh* mesh );
virtual ~Checkerboard();
Checkerboard_SN( Config* settings, Mesh* mesh );
virtual ~Checkerboard_SN();
virtual VectorVector GetScatteringXS( const Vector& energies );
virtual VectorVector GetTotalXS( const Vector& energies );
virtual std::vector<VectorVector> GetExternalSource( const Vector& energies );
virtual VectorVector SetupIC();
};
class Checkerboard_PN : public ProblemBase
{
private:
Vector _scatteringXS;
Vector _totalXS;
Checkerboard_PN() = delete;
bool isAbsorption( const Vector& pos ) const;
bool isSource( const Vector& pos ) const;
int GlobalIndex( int l, int k ) const;
public:
Checkerboard_PN( Config* settings, Mesh* mesh );
virtual ~Checkerboard_PN();
virtual VectorVector GetScatteringXS( const Vector& energies );
virtual VectorVector GetTotalXS( const Vector& energies );
......
......@@ -250,6 +250,7 @@ void Config::SetConfigOptions() {
// Entropy related options
/*! @brief Entropy Functional \n DESCRIPTION: Entropy functional used for the MN_Solver \n DEFAULT QUADRTATIC @ingroup Config. */
AddEnumOption( "ENTROPY_FUNCTIONAL", _entropyName, Entropy_Map, QUADRATIC );
AddEnumOption( "ENTROPY_FUNCTIONAL", _entropyName, Entropy_Map, MAXWELL_BOLTZMANN );
/*! @brief Optimizer Name \n DESCRIPTION: Optimizer used to determine the minimal Entropy reconstruction \n DEFAULT NEWTON \ingroup Config */
AddEnumOption( "ENTROPY_OPTIMIZER", _entropyOptimizerName, Optimizer_Map, NEWTON );
......
......@@ -7,7 +7,7 @@ EntropyBase* EntropyBase::Create( Config* settings ) {
switch( settings->GetEntropyName() ) {
case QUADRATIC: return new QuadraticEntropy();
case MAXWELL_BOLZMANN: return new MaxwellBoltzmannEntropy();
case MAXWELL_BOLTZMANN: return new MaxwellBoltzmannEntropy();
default: return new QuadraticEntropy();
}
}
......@@ -2,7 +2,7 @@
#include "common/config.h"
#include "common/mesh.h"
Checkerboard::Checkerboard( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
Checkerboard_SN::Checkerboard_SN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
_physics = nullptr;
_scatteringXS = Vector( _mesh->GetNumCells(), 1.0 );
_totalXS = Vector( _mesh->GetNumCells(), 1.0 );
......@@ -15,13 +15,13 @@ Checkerboard::Checkerboard( Config* settings, Mesh* mesh ) : ProblemBase( settin
}
}
Checkerboard::~Checkerboard() {}
Checkerboard_SN::~Checkerboard_SN() {}
VectorVector Checkerboard::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), _scatteringXS ); }
VectorVector Checkerboard_SN::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), _scatteringXS ); }
VectorVector Checkerboard::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), _totalXS ); }
VectorVector Checkerboard_SN::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), _totalXS ); }
std::vector<VectorVector> Checkerboard::GetExternalSource( const Vector& energies ) {
std::vector<VectorVector> Checkerboard_SN::GetExternalSource( const Vector& energies ) {
VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
auto cellMids = _mesh->GetCellMidPoints();
for( unsigned j = 0; j < cellMids.size(); ++j ) {
......@@ -30,12 +30,12 @@ std::vector<VectorVector> Checkerboard::GetExternalSource( const Vector& energie
return std::vector<VectorVector>( 1u, Q );
}
VectorVector Checkerboard::SetupIC() {
VectorVector Checkerboard_SN::SetupIC() {
VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
return psi;
}
bool Checkerboard::isAbsorption( const Vector& pos ) const {
bool Checkerboard_SN::isAbsorption( const Vector& pos ) const {
std::vector<double> lbounds{ 1, 2, 3, 4, 5 };
std::vector<double> ubounds{ 2, 3, 4, 5, 6 };
for( unsigned k = 0; k < lbounds.size(); ++k ) {
......@@ -49,9 +49,72 @@ bool Checkerboard::isAbsorption( const Vector& pos ) const {
return false;
}
bool Checkerboard::isSource( const Vector& pos ) const {
bool Checkerboard_SN::isSource( const Vector& pos ) const {
if( pos[0] >= 3 && pos[0] <= 4 && pos[1] >= 3 && pos[1] <= 4 )
return true;
else
return false;
}
////////////////////////////////////////////////////////////////////////////////////////////////////
Checkerboard_PN::Checkerboard_PN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
_physics = nullptr;
_scatteringXS = Vector( _mesh->GetNumCells(), 1.0 );
_totalXS = Vector( _mesh->GetNumCells(), 1.0 );
auto cellMids = _mesh->GetCellMidPoints();
for( unsigned j = 0; j < cellMids.size(); ++j ) {
if( isAbsorption( cellMids[j] ) ) {
_scatteringXS[j] = 0.0;
_totalXS[j] = 10.0;
}
}
}
Checkerboard_PN::~Checkerboard_PN() {}
VectorVector Checkerboard_PN::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), _scatteringXS ); }
VectorVector Checkerboard_PN::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), _totalXS ); }
std::vector<VectorVector> Checkerboard_PN::GetExternalSource( const Vector& energies ) {
VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
auto cellMids = _mesh->GetCellMidPoints();
for( unsigned j = 0; j < cellMids.size(); ++j ) {
if( isSource( cellMids[j] ) ) Q[j] = 1.0 / std::sqrt( 4 * M_PI ); // isotropic source
}
return std::vector<VectorVector>( 1u, Q );
}
VectorVector Checkerboard_PN::SetupIC() {
int ntotalEquations = GlobalIndex( _settings->GetMaxMomentDegree(), _settings->GetMaxMomentDegree() ) + 1;
VectorVector psi( _mesh->GetNumCells(), Vector( ntotalEquations, 1e-10 ) );
return psi;
}
bool Checkerboard_PN::isAbsorption( const Vector& pos ) const {
std::vector<double> lbounds{ 1, 2, 3, 4, 5 };
std::vector<double> ubounds{ 2, 3, 4, 5, 6 };
for( unsigned k = 0; k < lbounds.size(); ++k ) {
for( unsigned l = 0; l < lbounds.size(); ++l ) {
if( ( l + k ) % 2 == 1 || ( k == 2 && l == 2 ) || ( k == 2 && l == 4 ) ) continue;
if( pos[0] >= lbounds[k] && pos[0] <= ubounds[k] && pos[1] >= lbounds[l] && pos[1] <= ubounds[l] ) {
return true;
}
}
}
return false;
}
bool Checkerboard_PN::isSource( const Vector& pos ) const {
if( pos[0] >= 3 && pos[0] <= 4 && pos[1] >= 3 && pos[1] <= 4 )
return true;
else
return false;
}
int Checkerboard_PN::GlobalIndex( int l, int k ) const {
int numIndicesPrevLevel = l * l; // number of previous indices untill level l-1
int prevIndicesThisLevel = k + l; // number of previous indices in current level
return numIndicesPrevLevel + prevIndicesThisLevel;
}
......@@ -22,7 +22,12 @@ ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) {
else
return new LineSource_SN( settings, mesh ); // default
}
case PROBLEM_Checkerboard: return new Checkerboard( settings, mesh );
case PROBLEM_Checkerboard: {
if( settings->GetSolverName() == PN_SOLVER || settings->GetSolverName() == MN_SOLVER )
return new Checkerboard_PN( settings, mesh );
else
return new Checkerboard_SN( settings, mesh ); // default
}
case PROBLEM_ElectronRT: return new ElectronRT( settings, mesh );
case PROBLEM_WaterPhantom: return new WaterPhantom( settings, mesh );
case PROBLEM_LineSource_Pseudo_1D: return new LineSource_SN_Pseudo1D( settings, mesh );
......
......@@ -175,6 +175,8 @@ void MNSolver::Solve() {
( _sigmaT[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
}
psiNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
}
// Update Solution
......
......@@ -117,6 +117,7 @@ void PNSolver::Solve() {
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
}
}
psiNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
}
// Update Solution
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment