Commit 45f3488d authored by steffen.schotthoefer's avatar steffen.schotthoefer
Browse files

checkerboard IC supports monomial basis now. Renamed to checkerboard_moment...

checkerboard IC supports monomial basis now. Renamed to checkerboard_moment for clarity. mn solver works with 2d quadrature now
parent 21120f99
Pipeline #140154 failed with stage
...@@ -3,6 +3,8 @@ ...@@ -3,6 +3,8 @@
#include "problembase.h" #include "problembase.h"
class SphericalBase;
class Checkerboard_SN : public ProblemBase class Checkerboard_SN : public ProblemBase
{ {
private: private:
...@@ -24,13 +26,14 @@ class Checkerboard_SN : public ProblemBase ...@@ -24,13 +26,14 @@ class Checkerboard_SN : public ProblemBase
virtual VectorVector SetupIC() override; virtual VectorVector SetupIC() override;
}; };
class Checkerboard_PN : public ProblemBase class Checkerboard_Moment : public ProblemBase
{ {
private: private:
Vector _scatteringXS; /*!< @brief Vector of scattering crosssections len: numCells */ Vector _scatteringXS; /*!< @brief Vector of scattering crosssections len: numCells */
Vector _totalXS; /*!< @brief Vector of total crosssections. len: numCells*/ Vector _totalXS; /*!< @brief Vector of total crosssections. len: numCells*/
SphericalBase* _basis; /*!< @brief Class to compute and store current spherical harmonics basis */
Checkerboard_PN() = delete; Checkerboard_Moment() = delete;
bool isAbsorption( const Vector& pos ) const; /*!< @return True if pos is in absorption region, False otherwise */ bool isAbsorption( const Vector& pos ) const; /*!< @return True if pos is in absorption region, False otherwise */
bool isSource( const Vector& pos ) const; /*!< @return True if pos is in source region, False otherwise */ bool isSource( const Vector& pos ) const; /*!< @return True if pos is in source region, False otherwise */
...@@ -46,8 +49,8 @@ class Checkerboard_PN : public ProblemBase ...@@ -46,8 +49,8 @@ class Checkerboard_PN : public ProblemBase
int GlobalIndex( int l, int k ) const; int GlobalIndex( int l, int k ) const;
public: public:
Checkerboard_PN( Config* settings, Mesh* mesh ); Checkerboard_Moment( Config* settings, Mesh* mesh );
virtual ~Checkerboard_PN(); virtual ~Checkerboard_Moment();
virtual VectorVector GetScatteringXS( const Vector& energies ) override; virtual VectorVector GetScatteringXS( const Vector& energies ) override;
virtual VectorVector GetTotalXS( const Vector& energies ) override; virtual VectorVector GetTotalXS( const Vector& energies ) override;
......
...@@ -18,6 +18,7 @@ MESH_FILE = meshes/checkerboard.su2 ...@@ -18,6 +18,7 @@ MESH_FILE = meshes/checkerboard.su2
% ---- Problem specifications ---- % ---- Problem specifications ----
% %
PROBLEM = CHECKERBOARD PROBLEM = CHECKERBOARD
SPATIAL_DIM = 2
% %
% ---- Solver specifications ---- % ---- Solver specifications ----
% %
......
...@@ -3,7 +3,7 @@ ...@@ -3,7 +3,7 @@
UpwindFlux::UpwindFlux() : NumericalFlux() {} UpwindFlux::UpwindFlux() : NumericalFlux() {}
double UpwindFlux::Flux( const Vector& Omega, double psiL, double psiR, const Vector& n ) const { double UpwindFlux::Flux( const Vector& Omega, double psiL, double psiR, const Vector& n ) const {
double inner = Omega[0] * n[0] + Omega[1] * n[1]; double inner = Omega[0] * n[0] + Omega[1] * n[1]; // Only use x and y axis in 2d case
if( inner > 0 ) { if( inner > 0 ) {
return inner * psiL; return inner * psiL;
} }
......
#include "problems/checkerboard.h" #include "problems/checkerboard.h"
#include "common/config.h" #include "common/config.h"
#include "common/mesh.h" #include "common/mesh.h"
#include "toolboxes/sphericalbase.h"
// ---- Checkerboard Sn ---- // ---- Checkerboard Sn ----
// Constructor for Ckeckerboard case with Sn // Constructor for Ckeckerboard case with Sn
...@@ -66,9 +67,12 @@ bool Checkerboard_SN::isSource( const Vector& pos ) const { ...@@ -66,9 +67,12 @@ bool Checkerboard_SN::isSource( const Vector& pos ) const {
//////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////
// ---- Checkerboard Pn ---- // ---- Checkerboard Moments ----
// Constructor for checkerboard case with Pn // Constructor for checkerboard case with Pn
Checkerboard_PN::Checkerboard_PN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) { Checkerboard_Moment::Checkerboard_Moment( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
_basis = SphericalBase::Create(_settings);
_physics = nullptr; _physics = nullptr;
// Initialise crosssections = 1 (scattering) // Initialise crosssections = 1 (scattering)
...@@ -85,13 +89,15 @@ Checkerboard_PN::Checkerboard_PN( Config* settings, Mesh* mesh ) : ProblemBase( ...@@ -85,13 +89,15 @@ Checkerboard_PN::Checkerboard_PN( Config* settings, Mesh* mesh ) : ProblemBase(
} }
} }
Checkerboard_PN::~Checkerboard_PN() {} Checkerboard_Moment::~Checkerboard_Moment() {
delete _basis;
}
VectorVector Checkerboard_PN::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), _scatteringXS ); } VectorVector Checkerboard_Moment::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), _scatteringXS ); }
VectorVector Checkerboard_PN::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), _totalXS ); } VectorVector Checkerboard_Moment::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), _totalXS ); }
std::vector<VectorVector> Checkerboard_PN::GetExternalSource( const Vector& /*energies*/ ) { std::vector<VectorVector> Checkerboard_Moment::GetExternalSource( const Vector& /*energies*/ ) {
VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) ); VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
auto cellMids = _mesh->GetCellMidPoints(); auto cellMids = _mesh->GetCellMidPoints();
for( unsigned j = 0; j < cellMids.size(); ++j ) { for( unsigned j = 0; j < cellMids.size(); ++j ) {
...@@ -100,13 +106,13 @@ std::vector<VectorVector> Checkerboard_PN::GetExternalSource( const Vector& /*en ...@@ -100,13 +106,13 @@ std::vector<VectorVector> Checkerboard_PN::GetExternalSource( const Vector& /*en
return std::vector<VectorVector>( 1u, Q ); return std::vector<VectorVector>( 1u, Q );
} }
VectorVector Checkerboard_PN::SetupIC() { VectorVector Checkerboard_Moment::SetupIC() {
int ntotalEquations = GlobalIndex( _settings->GetMaxMomentDegree(), _settings->GetMaxMomentDegree() ) + 1; int ntotalEquations = _basis->GetBasisSize();
VectorVector psi( _mesh->GetNumCells(), Vector( ntotalEquations, 1e-10 ) ); VectorVector psi( _mesh->GetNumCells(), Vector( ntotalEquations, 1e-10 ) );
return psi; return psi;
} }
bool Checkerboard_PN::isAbsorption( const Vector& pos ) const { bool Checkerboard_Moment::isAbsorption( const Vector& pos ) const {
// Check whether pos is in absorption region // Check whether pos is in absorption region
std::vector<double> lbounds{ 1, 2, 3, 4, 5 }; std::vector<double> lbounds{ 1, 2, 3, 4, 5 };
std::vector<double> ubounds{ 2, 3, 4, 5, 6 }; std::vector<double> ubounds{ 2, 3, 4, 5, 6 };
...@@ -121,16 +127,10 @@ bool Checkerboard_PN::isAbsorption( const Vector& pos ) const { ...@@ -121,16 +127,10 @@ bool Checkerboard_PN::isAbsorption( const Vector& pos ) const {
return false; return false;
} }
bool Checkerboard_PN::isSource( const Vector& pos ) const { bool Checkerboard_Moment::isSource( const Vector& pos ) const {
// Check whether pos is in source region // Check whether pos is in source region
if( pos[0] >= 3 && pos[0] <= 4 && pos[1] >= 3 && pos[1] <= 4 ) if( pos[0] >= 3 && pos[0] <= 4 && pos[1] >= 3 && pos[1] <= 4 )
return true; return true;
else else
return false; 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;
}
...@@ -30,7 +30,7 @@ ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) { ...@@ -30,7 +30,7 @@ ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) {
} }
case PROBLEM_Checkerboard: { case PROBLEM_Checkerboard: {
if( settings->GetSolverName() == PN_SOLVER || settings->GetSolverName() == MN_SOLVER ) if( settings->GetSolverName() == PN_SOLVER || settings->GetSolverName() == MN_SOLVER )
return new Checkerboard_PN( settings, mesh ); return new Checkerboard_Moment( settings, mesh );
else else
return new Checkerboard_SN( settings, mesh ); // default return new Checkerboard_SN( settings, mesh ); // default
} }
......
Markdown is supported
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