Commit a4e790e7 authored by steffen.schotthoefer's avatar steffen.schotthoefer
Browse files

added option for scatter coeff for linesource. added compatibility factor for...

added option for scatter coeff for linesource. added compatibility factor for rad flux output and ic for moment solvers
parent 66aa4556
Pipeline #111811 failed with stage
in 14 minutes and 52 seconds
......@@ -58,6 +58,9 @@ class Config
unsigned short _maxMomentDegree; /*!< @brief Maximal Order of Moments for PN and MN Solver */
unsigned short _reconsOrder; /*!< @brief Spatial Order of Accuracy for Solver */
// Linesource
double _sigmaS; /*!< @brief Scattering coeffient for Linesource test case */
/*!< @brief If true, very low entries (10^-10 or smaller) of the flux matrices will be set to zero,
* to improve floating point accuracy */
bool _cleanFluxMat;
......@@ -239,6 +242,8 @@ class Config
bool inline GetSNAllGaussPts() const { return _allGaussPts; }
bool inline GetIsCSD() const { return _csd; }
// Linesource
double inline GetSigmaS() const { return _sigmaS; }
// Optimizer
double inline GetNewtonOptimizerEpsilon() const { return _optimizerEpsilon; }
unsigned inline GetNewtonIter() const { return _newtonIter; }
......
......@@ -8,6 +8,9 @@ class LineSource : public ProblemBase
private:
LineSource() = delete;
protected:
double _sigmaS; /*!< @brief: Scattering coefficient */
public:
LineSource( Config* settings, Mesh* mesh );
......
......@@ -18,7 +18,8 @@ MESH_FILE = linesource.su2
%MESH_FILE = linesource_debug.su2
%
PROBLEM = LINESOURCE
%
SCATTER_COEFF = 0
% ---- Solver specifications ----
%
% Solver type
......@@ -28,7 +29,7 @@ CFL_NUMBER = 0.7
% Final time for simulation
TIME_FINAL = 0.3
% Maximal Moment degree
MAX_MOMENT_SOLVER = 5
MAX_MOMENT_SOLVER = 1
%
%% Entropy settings
ENTROPY_FUNCTIONAL = MAXWELL_BOLZMANN
......@@ -55,9 +56,9 @@ QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
%
% Quadrature Order
QUAD_ORDER = 12
QUAD_ORDER = 8
%
%
% ----- Output ----
%
VOLUME_OUTPUT = ( MINIMAL, MOMENTS, DUAL_MOMENTS)
VOLUME_OUTPUT = (ANALYTIC, MINIMAL, MOMENTS, DUAL_MOMENTS)
......@@ -19,6 +19,8 @@ MESH_FILE = linesource.su2
%
PROBLEM = LINESOURCE
SCATTER_COEFF = 0
%
%
% ---- Solver specifications ----
%
......
% This is a comment
%%%%%%%%%%%%%%%%%%%%%%
% Example config file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Linesource Benchmarking File SN %
% Author <Steffen Schotthöfer> %
% Date 01.10.2020 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% ---- File specifications ----
%
......@@ -14,6 +17,8 @@ LOG_DIR = ../result/logs
MESH_FILE = linesource.su2
%
PROBLEM = LINESOURCE
SCATTER_COEFF = 0
%
% ---- Solver specifications ----
%
......@@ -33,4 +38,4 @@ BC_DIRICHLET = ( void )
% Quadrature Type
QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
% Quadrature Order
QUAD_ORDER = 6
QUAD_ORDER = 8
......@@ -243,6 +243,10 @@ void Config::SetConfigOptions() {
/*! @brief SN_ALL_GAUSS_PTS \n DESCRIPTION: If true, the SN Solver uses all Gauss Quadrature Points for 2d. \n DEFAULT false \ingroup Config */
AddBoolOption( "SN_ALL_GAUSS_PTS", _allGaussPts, false );
// Linesource Testcase Options
/*! @brief SCATTER_COEFF \n DESCRIPTION: Sets the scattering coefficient for the Linesource test case. \n DEFAULT 0.0 \ingroup Config */
AddDoubleOption( "SCATTER_COEFF", _sigmaS, 0.0 );
// 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 );
......
......@@ -7,7 +7,10 @@
// ---- Linesource ----
LineSource::LineSource( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) { _physics = nullptr; }
LineSource::LineSource( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
_physics = nullptr;
_sigmaS = settings->GetSigmaS();
}
LineSource::~LineSource() {}
double LineSource::GetAnalyticalSolution( double x, double y, double t, double sigma_s ) {
......@@ -15,7 +18,7 @@ double LineSource::GetAnalyticalSolution( double x, double y, double t, double s
double solution = 0.0;
double R = sqrt( x * x + y * y );
if( sigma_s == 0.0 ) {
if( _sigmaS == 0.0 ) {
if( ( t - R ) > 0 ) {
solution = 1 / ( 2 * M_PI * t * sqrt( t * t - R * R ) );
}
......@@ -59,21 +62,16 @@ double LineSource::HelperRho_ptc1( double R, double t ) {
}
double LineSource::HelperRho_ptc2( double R, double t ) {
double gamma = R / t;
double result = 0;
double integral = 0;
double gamma = R / t;
double result = 0;
if( 1 - gamma > 0 ) {
// Compute the integralpart with midpoint rule
result = exp( -t ) / ( 32 * M_PI * M_PI * R ) * ( 1 - gamma * gamma ) * HelperIntRho_ptc2( t, gamma );
}
return result;
}
double LineSource::HelperIntRho_ptc2( double t, double gamma ) {
double u = 0;
double integral = 0;
std::complex<double> beta( 0, 0 );
......@@ -109,9 +107,11 @@ LineSource_SN::LineSource_SN( Config* settings, Mesh* mesh ) : LineSource( setti
LineSource_SN::~LineSource_SN() {}
VectorVector LineSource_SN::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) ); }
VectorVector LineSource_SN::GetScatteringXS( const Vector& energies ) {
return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), _sigmaS ) );
}
VectorVector LineSource_SN::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) ); }
VectorVector LineSource_SN::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), _sigmaS ) ); }
std::vector<VectorVector> LineSource_SN::GetExternalSource( const Vector& energies ) {
return std::vector<VectorVector>( 1u, std::vector<Vector>( _mesh->GetNumCells(), Vector( 1u, 0.0 ) ) );
......@@ -168,9 +168,11 @@ LineSource_PN::LineSource_PN( Config* settings, Mesh* mesh ) : LineSource( setti
LineSource_PN::~LineSource_PN() {}
VectorVector LineSource_PN::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) ); }
VectorVector LineSource_PN::GetScatteringXS( const Vector& energies ) {
return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), _sigmaS ) );
}
VectorVector LineSource_PN::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 1.0 ) ); }
VectorVector LineSource_PN::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), _sigmaS ) ); }
std::vector<VectorVector> LineSource_PN::GetExternalSource( const Vector& energies ) {
return std::vector<VectorVector>( 1u, std::vector<Vector>( _mesh->GetNumCells(), Vector( 1u, 0.0 ) ) );
......@@ -185,33 +187,10 @@ VectorVector LineSource_PN::SetupIC() {
// Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments are zero.
double t = 3.2e-4; // pseudo time for gaussian smoothing (Approx to dirac impulse)
// double offset = 0.099;
for( unsigned j = 0; j < cellMids.size(); ++j ) {
double x = cellMids[j][0];
double y = cellMids[j][1]; // (x- 0.5) * (x- 0.5)
psi[j][0] = 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x ) * ( x ) + ( y ) * ( y ) ) / ( 4 * t ) );
//+
// 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x - offset ) * ( x - offset ) + ( y - offset ) * ( y - offset ) ) / ( 4 * t
// ) ) + 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x + offset ) * ( x + offset ) + ( y - offset ) * ( y - offset ) ) / (
// 4 * t ) ) + 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x - offset ) * ( x - offset ) + ( y + offset ) * ( y + offset )
// ) / ( 4 * t ) ) + 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x + offset ) * ( x + offset ) + ( y + offset ) * ( y +
// offset ) ) / ( 4 * t ) ) + 1.0 / ( 4.0 * M_PI * t ) *
// std::exp( -( ( x - 2 * offset ) * ( x - 2 * offset ) + ( y - 2 * offset ) * ( y - 2 * offset ) ) / ( 4 * t ) ) +
// 1.0 / ( 4.0 * M_PI * t ) *
// std::exp( -( ( x + 2 * offset ) * ( x + 2 * offset ) + ( y - 2 * offset ) * ( y - 2 * offset ) ) / ( 4 * t ) ) +
// 1.0 / ( 4.0 * M_PI * t ) *
// std::exp( -( ( x - 2 * offset ) * ( x - 2 * offset ) + ( y + 2 * offset ) * ( y + 2 * offset ) ) / ( 4 * t ) ) +
// 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x + 2 * offset ) * ( x + 2 * offset ) + ( y + 2 * offset ) * ( y + 2 *
// offset ) ) / ( 4 * t ) );
double x = cellMids[j][0];
double y = cellMids[j][1]; // (x- 0.5) * (x- 0.5)
psi[j][0] = sqrt( 4 * M_PI ) * 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) );
}
// Debugging jump test case
// for( unsigned j = 0; j < cellMids.size(); ++j ) {
// if( cellMids[j][0] > -0.2 && cellMids[j][1] > -0.2 && cellMids[j][1] < 0.2 && cellMids[j][0] < 0.2 )
// psi[j][0] = 1.0;
// else
// psi[j][0] = 0.0;
//}
return psi;
}
......@@ -31,14 +31,6 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ) {
_quadPointsSphere = _quadrature->GetPointsSphere();
_settings->SetNQuadPoints( _nq );
// This must be shifted to Problem Class
for( unsigned n = 0; n < _nEnergies; n++ ) {
for( unsigned j = 0; j < _nCells; j++ ) {
_sigmaT[n][j] = 1; //_sigmaT[n][j] - _sigmaS[n][j];
_sigmaS[n][j] = 1;
}
}
// Initialize Scatter Matrix --
_scatterMatDiag = Vector( _nTotalEntries, 0.0 );
ComputeScatterMatrix();
......@@ -263,8 +255,9 @@ void MNSolver::PrepareOutputFields() {
}
double MNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
double mass = 0.0;
unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
double mass = 0.0;
unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
double firstMomentScaleFactor = sqrt( 4 * M_PI );
// Compute total "mass" of the system ==> to check conservation properties
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
......@@ -275,7 +268,7 @@ double MNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
switch( _settings->GetVolumeOutput()[idx_group] ) {
case MINIMAL:
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
_outputFields[idx_group][0][idx_cell] = _sol[idx_cell][0];
_outputFields[idx_group][0][idx_cell] = firstMomentScaleFactor * _sol[idx_cell][0];
}
break;
case MOMENTS:
......@@ -297,10 +290,11 @@ double MNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
double time = idx_pseudoTime * _dE;
double sigma = 1;
double sigma = 0;
_outputFields[idx_group][0][idx_cell] = _problem->GetAnalyticalSolution(
_mesh->GetCellMidPoints()[idx_cell][0], _mesh->GetCellMidPoints()[idx_cell][1], time, sigma );
_outputFields[idx_group][0][idx_cell] =
( 4 * M_PI ) * _problem->GetAnalyticalSolution(
_mesh->GetCellMidPoints()[idx_cell][0], _mesh->GetCellMidPoints()[idx_cell][1], time, sigma );
}
break;
......
......@@ -13,14 +13,6 @@ PNSolver::PNSolver( Config* settings ) : Solver( settings ) {
_LMaxDegree = settings->GetMaxMomentDegree();
_nTotalEntries = GlobalIndex( int( _LMaxDegree ), int( _LMaxDegree ) ) + 1;
// This must be shifted to Problem Class
for( unsigned n = 0; n < _nEnergies; n++ ) {
for( unsigned j = 0; j < _nCells; j++ ) {
_sigmaT[n][j] = 1;
_sigmaS[n][j] = 1;
}
}
// Initialize System Matrices
_Ax = SymMatrix( _nTotalEntries );
_Ay = SymMatrix( _nTotalEntries );
......@@ -42,33 +34,18 @@ PNSolver::PNSolver( Config* settings ) : Solver( settings ) {
// Fill System Matrices
ComputeSystemMatrices();
// std::cout << "System Matrix Set UP!" << std::endl;
// Compute Decomposition in positive and negative (eigenvalue) parts of flux jacobians
ComputeFluxComponents();
// Compute diagonal of the scatter matrix (it's a diagonal matrix)
ComputeScatterMatrix();
// std::cout << "scatterMatrix : " << _scatterMatDiag << "\n";
// AdaptTimeStep();
if( settings->GetCleanFluxMat() ) CleanFluxMatrices();
// std::cout << "--------\n";
// std::cout << "_Ax :\n" << _Ax << "\n "; // _AxP \n" << _AxPlus << "\n _AxM \n" << _AxMinus << "\n";
// std::cout << "_Ay :\n" << _Ay << "\n "; //_AyP \n" << _AyPlus << "\n _AyM \n" << _AyMinus << "\n";
// std::cout << "_Az :\n" << _Az << "\n "; //_AzP \n" << _AzPlus << "\n _AzM \n" << _AzMinus << "\n";
//
// std::cout << "_AxPlus :\n" << _AxPlus << "\n "; // _AxP \n" << _AxPlus << "\n _AxM \n" << _AxMinus << "\n";
// std::cout << "_AyPlus :\n" << _AyPlus << "\n "; //_AyP \n" << _AyPlus << "\n _AyM \n" << _AyMinus << "\n";
// std::cout << "_AzPlus :\n" << _AzPlus << "\n ";
//
// std::cout << "_AxMinus :\n" << _AxMinus << "\n "; // _AxP \n" << _AxPlus << "\n _AxM \n" << _AxMinus << "\n";
// std::cout << "_AyMinus :\n" << _AyMinus << "\n "; //_AyP \n" << _AyPlus << "\n _AyM \n" << _AyMinus << "\n";
// std::cout << "_AzMinus :\n" << _AzMinus << "\n ";
//
// std::cout << "_nCells: " << _nCells << "\n";
// Compute moments of initial condition
// TODO
// Solver output
PrepareOutputFields();
......@@ -91,7 +68,9 @@ void PNSolver::Solve() {
if( rank == 0 ) log->info( "{:10} {:10}", "t", "mass" );
// if( rank == 0 ) log->info( "{:03.8f} {:01.5e} {:01.5e}", -1.0, dFlux, mass1 );
// Remove
mass = WriteOutputFields();
if( rank == 0 ) log->info( " {:01.5e} {:01.5e}", 0.0, mass );
// Loop over energies (pseudo-time of continuous slowing down approach)
for( unsigned idx_energy = 0; idx_energy < _nEnergies; idx_energy++ ) {
......@@ -393,8 +372,9 @@ void PNSolver::PrepareOutputFields() {
}
double PNSolver::WriteOutputFields() {
double mass = 0.0;
unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
double mass = 0.0;
unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
double firstMomentScaleFactor = sqrt( 4 * M_PI );
// Compute total "mass" of the system ==> to check conservation properties
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
......@@ -411,7 +391,7 @@ double PNSolver::WriteOutputFields() {
case MOMENTS:
for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
_outputFields[idx_group][idx_sys][idx_cell] = _sol[idx_cell][idx_sys];
_outputFields[idx_group][idx_sys][idx_cell] = firstMomentScaleFactor * _sol[idx_cell][idx_sys];
}
}
break;
......
......@@ -68,7 +68,18 @@ void SNSolver::Solve() {
Vector fluxOld( _nCells, 0.0 );
int rank;
MPI_Comm_rank( MPI_COMM_WORLD, &rank );
if( rank == 0 ) log->info( "{:10} {:10}", "t", "dFlux" );
if( rank == 0 ) log->info( "{:10} {:10}", "t", "dFlux", "mass" );
// Remove
double mass1 = 0;
for( unsigned i = 0; i < _nCells; ++i ) {
_solverOutput[i] = dot( _sol[i], _weights );
}
// Compute total "mass" of the system ==> to check conservation properties
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
mass1 += _solverOutput[idx_cell] * _areas[idx_cell]; // Should probably go to postprocessing
}
if( rank == 0 ) log->info( "{:01.5e} {:01.5e} {:01.5e}", 0.0, dFlux, mass1 );
// loop over energies (pseudo-time)
for( unsigned n = 0; n < _nEnergies; ++n ) {
......@@ -144,10 +155,18 @@ void SNSolver::Solve() {
fluxNew[i] = dot( _sol[i], _weights );
_solverOutput[i] = fluxNew[i];
}
double mass = 0;
// Compute total "mass" of the system ==> to check conservation properties
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
mass += _solverOutput[idx_cell] * _areas[idx_cell]; // Should probably go to postprocessing
}
// Save( n );
dFlux = blaze::l2Norm( fluxNew - fluxOld );
fluxOld = fluxNew;
if( rank == 0 ) log->info( "{:03.8f} {:01.5e}", _energies[n], dFlux );
if( rank == 0 ) log->info( "{:03.8f} {:01.5e} {:01.5e}", _energies[n], dFlux, mass );
}
}
......
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