Commit d89104bf authored by pia.stammer's avatar pia.stammer
Browse files

Added testcases for air cavity in water and muscle-bone-lung

parent e757df83
Pipeline #117581 failed with stages
in 26 minutes and 44 seconds
......@@ -73,13 +73,17 @@ enum PROBLEM_NAME {
PROBLEM_ElectronRT,
PROBLEM_WaterPhantom,
PROBLEM_LineSource_Pseudo_1D,
PROBLEM_LineSource_Pseudo_1D_Physics
PROBLEM_LineSource_Pseudo_1D_Physics,
PROBLEM_AirCavity,
PROBLEM_MuscleBoneLung
};
inline std::map<std::string, PROBLEM_NAME> Problem_Map{ { "LINESOURCE", PROBLEM_LineSource },
{ "CHECKERBOARD", PROBLEM_Checkerboard },
{ "ELECTRONRT", PROBLEM_ElectronRT },
{ "WATERPHANTOM", PROBLEM_WaterPhantom },
{ "AIRCAVITY", PROBLEM_AirCavity },
{ "MUSCLEBONELUNG", PROBLEM_MuscleBoneLung },
{ "LINESOURCE_PSEUDO_1D", PROBLEM_LineSource_Pseudo_1D },
{ "LINESOURCE_PSEUDO_1D_PHYSICS", PROBLEM_LineSource_Pseudo_1D_Physics } };
......
#ifndef AIRCAVITY
#define AIRCAVITY
#include "electronrt.h"
class AirCavity1D : public ElectronRT
{
private:
AirCavity1D() = delete;
public:
AirCavity1D( Config* settings, Mesh* mesh );
virtual ~AirCavity1D();
virtual std::vector<VectorVector> GetExternalSource( const Vector& energies );
virtual VectorVector SetupIC();
std::vector<double> GetDensity( const VectorVector& cellMidPoints );
};
#endif // AIRCAVITY
\ No newline at end of file
#ifndef MUSCLEBONELUNG
#define MUSCLEBONELUNG
#include "electronrt.h"
class MuscleBoneLung : public ElectronRT
{
private:
MuscleBoneLung() = delete;
public:
MuscleBoneLung( Config* settings, Mesh* mesh );
virtual ~MuscleBoneLung();
virtual std::vector<VectorVector> GetExternalSource( const Vector& energies );
virtual VectorVector SetupIC();
std::vector<double> GetDensity( const VectorVector& cellMidPoints );
};
#endif // MUSCLEBONELUNG
\ No newline at end of file
OUTPUT_DIR = ../result
OUTPUT_FILE = aircavity_csd_1d_10MeV_fine
LOG_DIR = ../result/logs
MESH_FILE = 1DMesh.su2
PROBLEM = AIRCAVITY
SOLVER = CSD_SN_FOKKERPLANCK_SOLVER
CONTINUOUS_SLOWING_DOWN = YES
HYDROGEN_FILE = ENDL_H.txt
OXYGEN_FILE = ENDL_O.txt
KERNEL = ISOTROPIC_1D
CFL_NUMBER = 0.01
TIME_FINAL = 1.0
CLEAN_FLUX_MATRICES = NO
BC_DIRICHLET = ( dirichlet )
BC_NEUMANN = ( wall_low, wall_up )
QUAD_TYPE = GAUSS_LEGENDRE_1D
QUAD_ORDER = 32
......@@ -37,4 +37,4 @@ BC_DIRICHLET = ( void )
% Quadrature Type
QUAD_TYPE = LEBEDEV
% Quadrature Order
QUAD_ORDER = 4
QUAD_ORDER = 3
OUTPUT_DIR = ../result
OUTPUT_FILE = musclebonelung_csd_1d_10MeV_fine
LOG_DIR = ../result/logs
MESH_FILE = 1DMesh.su2
PROBLEM = MUSCLEBONELUNG
SOLVER = CSD_SN_FOKKERPLANCK_SOLVER
CONTINUOUS_SLOWING_DOWN = YES
HYDROGEN_FILE = ENDL_H.txt
OXYGEN_FILE = ENDL_O.txt
KERNEL = ISOTROPIC_1D
CFL_NUMBER = 0.01
TIME_FINAL = 1.0
CLEAN_FLUX_MATRICES = NO
BC_DIRICHLET = ( dirichlet )
BC_NEUMANN = ( wall_low, wall_up )
QUAD_TYPE = GAUSS_LEGENDRE_1D
QUAD_ORDER = 32
\ No newline at end of file
#include "problems/aircavity1d.h"
#include "common/config.h"
#include "common/mesh.h"
AirCavity1D::AirCavity1D( Config* settings, Mesh* mesh ) : ElectronRT( settings, mesh ) {}
AirCavity1D::~AirCavity1D() { delete _physics; }
std::vector<VectorVector> AirCavity1D::GetExternalSource( const Vector& energies ) {
return std::vector<VectorVector>( energies.size(), std::vector<Vector>( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) ) );
}
VectorVector AirCavity1D::SetupIC() {
VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
auto cellMids = _mesh->GetCellMidPoints();
double s = 0.1;
for( unsigned j = 0; j < cellMids.size(); ++j ) {
double x = cellMids[j][0];
psi[j][_settings->GetNQuadPoints() - 1] = 1.0 / ( s * sqrt( 2 * M_PI ) ) * std::exp( -x * x / ( 2 * s * s ) );
}
return psi;
}
std::vector<double> AirCavity1D::GetDensity( const VectorVector& cellMidPoints ) {
std::vector<double> densities ( 4*cellMidPoints.size()/9, 1.0 );
std::vector<double> air( 2*cellMidPoints.size()/9, 0.00125 );
std::vector<double> water ( 3*cellMidPoints.size()/9, 1.0 );
densities.insert(densities.end(),air.begin(),air.end()-1);
densities.insert(densities.end(),water.begin(),water.end()-1);
return densities; }
#include "problems/musclebonelung.h"
#include "common/config.h"
#include "common/mesh.h"
MuscleBoneLung::MuscleBoneLung( Config* settings, Mesh* mesh ) : ElectronRT( settings, mesh ) {}
MuscleBoneLung::~MuscleBoneLung() { delete _physics; }
std::vector<VectorVector> MuscleBoneLung::GetExternalSource( const Vector& energies ) {
return std::vector<VectorVector>( energies.size(), std::vector<Vector>( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) ) );
}
VectorVector MuscleBoneLung::SetupIC() {
VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
auto cellMids = _mesh->GetCellMidPoints();
double s = 0.1;
for( unsigned j = 0; j < cellMids.size(); ++j ) {
double x = cellMids[j][0];
psi[j][_settings->GetNQuadPoints() - 1] = 1.0 / ( s * sqrt( 2 * M_PI ) ) * std::exp( -x * x / ( 2 * s * s ) );
}
return psi;
}
std::vector<double> MuscleBoneLung::GetDensity( const VectorVector& cellMidPoints ) {
std::vector<double> densities ( 1.5*cellMidPoints.size()/9, 1.05 ); //muscle layer
std::vector<double> bone( 1.5*cellMidPoints.size()/9, 1.92 );
std::vector<double> lung ( 6*cellMidPoints.size()/9, 0.26 );
densities.insert(densities.end(),bone.begin(),bone.end()-1);
densities.insert(densities.end(),lung.begin(),lung.end()-1);
return densities; }
......@@ -5,6 +5,8 @@
#include "problems/linesource.h"
#include "problems/problembase.h"
#include "problems/waterphantom.h"
#include "problems/aircavity1d.h"
#include "problems/musclebonelung.h"
ProblemBase::ProblemBase( Config* settings, Mesh* mesh ) {
_settings = settings;
......@@ -29,6 +31,8 @@ ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) {
return new Checkerboard_SN( settings, mesh ); // default
}
case PROBLEM_ElectronRT: return new ElectronRT( settings, mesh );
case PROBLEM_AirCavity: return new AirCavity1D( settings, mesh );
case PROBLEM_MuscleBoneLung: return new MuscleBoneLung( settings, mesh );
case PROBLEM_WaterPhantom: return new WaterPhantom( settings, mesh );
case PROBLEM_LineSource_Pseudo_1D: return new LineSource_SN_Pseudo1D( settings, mesh );
case PROBLEM_LineSource_Pseudo_1D_Physics: return new LineSource_SN_Pseudo1D_Physics( settings, mesh );
......
......@@ -13,6 +13,7 @@ std::vector<VectorVector> WaterPhantom::GetExternalSource( const Vector& energie
VectorVector WaterPhantom::SetupIC() {
VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
auto cellMids = _mesh->GetCellMidPoints();
std::cout<<"Length of mesh "<< cellMids.size() << std::endl;
double s = 0.1;
for( unsigned j = 0; j < cellMids.size(); ++j ) {
double x = cellMids[j][0];
......
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