Skip to content

GitLab

Commit 5388b5cc authored by steffen.schotthoefer's avatar steffen.schotthoefer
Browse files

analytical solution works for gamma=1

parent 379ab8f9
Pipeline #111887 passed with stage
in 15 minutes and 11 seconds
Hide whitespace changes
Inline Side-by-side
code/input/exampleMN.cfg
View file @ 5388b5cc
......@@ -18,7 +18,7 @@ MESH_FILE = linesource.su2
%MESH_FILE = linesource_debug.su2
%
PROBLEM = LINESOURCE
SCATTER_COEFF = 0
SCATTER_COEFF = 1
% ---- Solver specifications ----
%
......@@ -29,7 +29,7 @@ CFL_NUMBER = 0.7
% Final time for simulation
TIME_FINAL = 0.3
% Maximal Moment degree
MAX_MOMENT_SOLVER = 2
MAX_MOMENT_SOLVER = 3
%
%% Entropy settings
ENTROPY_FUNCTIONAL = MAXWELL_BOLZMANN
......@@ -62,3 +62,4 @@ QUAD_ORDER = 8
% ----- Output ----
%
VOLUME_OUTPUT = (ANALYTIC, MINIMAL, MOMENTS, DUAL_MOMENTS)
OUTPUT_FREQUENCY = 1
code/input/examplePN.cfg
View file @ 5388b5cc
......@@ -19,7 +19,7 @@ MESH_FILE = linesource.su2
%
PROBLEM = LINESOURCE
SCATTER_COEFF = 0
SCATTER_COEFF = 1
%
%
% ---- Solver specifications ----
......
code/input/exampleSN.cfg
View file @ 5388b5cc
......@@ -17,7 +17,7 @@ LOG_DIR = ../result/logs
MESH_FILE = linesource.su2
%
PROBLEM = LINESOURCE
SCATTER_COEFF = 0
SCATTER_COEFF = 1
%
% ---- Solver specifications ----
......
code/src/problems/linesource.cpp
View file @ 5388b5cc
......@@ -18,19 +18,24 @@ double LineSource::GetAnalyticalSolution( double x, double y, double t, double s
double solution = 0.0;
double R = sqrt( x * x + y * y );
if( _sigmaS == 0.0 ) {
if( ( t - R ) > 0 ) {
solution = 1 / ( 2 * M_PI * t * sqrt( t * t - R * R ) );
if( t > 0 ) {
if( _sigmaS == 0.0 ) {
if( ( t - R ) > 0 ) {
solution = 1 / ( 2 * M_PI * t * sqrt( t * t - R * R ) );
}
}
}
else {
double gamma = R / t;
else {
double gamma = R / t;
if( ( 1 - gamma ) > 0 ) {
solution = exp( -t ) / ( 2 * M_PI * t * t * sqrt( 1 - gamma * gamma ) ) + 2 * t * HelperIntRho_ptc( R, t );
if( ( 1 - gamma ) > 0 ) {
solution = exp( -t ) / ( 2 * M_PI * t * t * sqrt( 1 - gamma * gamma ) ) + 2 * t * HelperIntRho_ptc( R, t );
}
}
}
return solution;
else {
solution = 0;
}
return ( 4 * M_PI ) * solution;
}
double LineSource::HelperIntRho_ptc( double R, double t ) {
......@@ -42,10 +47,10 @@ double LineSource::HelperIntRho_ptc( double R, double t ) {
// integral is from 0 to sqrt( 1 - gamma * gamma )
double stepsize = sqrt( 1 - gamma * gamma ) / (double)numsteps;
omega = 0.5 * stepsize;
for( int i = 0; i < numsteps; i++ ) {
omega = omega + i * stepsize;
integral += HelperRho_ptc( t * sqrt( gamma * gamma + omega * omega ), t );
omega = i * stepsize + 0.5 * stepsize;
integral += stepsize * HelperRho_ptc( t * sqrt( gamma * gamma + omega * omega ), t );
}
return integral;
}
......@@ -68,6 +73,10 @@ double LineSource::HelperRho_ptc2( double R, double t ) {
// Compute the integralpart with midpoint rule
result = exp( -t ) / ( 32 * M_PI * M_PI * R ) * ( 1 - gamma * gamma ) * HelperIntRho_ptc2( t, gamma );
}
if( __isnan( result ) ) {
double iNan = 0.0;
}
return result;
}
......@@ -86,17 +95,15 @@ double LineSource::HelperIntRho_ptc2( double t, double gamma ) {
double stepsize = M_PI / (double)numsteps;
u = 0.5 * stepsize;
q = ( 1.0 + gamma ) / ( 1.0 - gamma );
for( int i = 0; i < numsteps; i++ ) {
u = u + (double)i * stepsize;
u = i * stepsize + 0.5 * stepsize;
// function evaluation
beta = ( log( q ) + com_one * u ) / ( gamma + com_one * tan( 0.5 * u ) );
complexPart = ( gamma + com_one * tan( 0.5 * u ) ) * beta * beta * beta * exp( 0.5 * t * ( 1 - gamma * gamma ) * beta );
integral += ( 1 / cos( 0.5 * u ) ) * ( 1 / cos( 0.5 * u ) ) * complexPart.real();
integral += stepsize * ( 1 / cos( 0.5 * u ) ) * ( 1 / cos( 0.5 * u ) ) * complexPart.real();
}
return integral;
}
......
code/src/solvers/mnsolver.cpp
View file @ 5388b5cc
......@@ -264,44 +264,39 @@ double MNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
mass += _sol[idx_cell][0] * _areas[idx_cell]; // Should probably go to postprocessing
}
if( _settings->GetOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetOutputFrequency() == 0 ) {
for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
switch( _settings->GetVolumeOutput()[idx_group] ) {
case MINIMAL:
for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
switch( _settings->GetVolumeOutput()[idx_group] ) {
case MINIMAL:
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
_outputFields[idx_group][0][idx_cell] = firstMomentScaleFactor * _sol[idx_cell][0];
}
break;
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][0][idx_cell] = firstMomentScaleFactor * _sol[idx_cell][0];
}
break;
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];
}
}
break;
case DUAL_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] = _alpha[idx_cell][idx_sys];
}
_outputFields[idx_group][idx_sys][idx_cell] = _sol[idx_cell][idx_sys];
}
break;
case ANALYTIC:
// Compute total "mass" of the system ==> to check conservation properties
}
break;
case DUAL_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] = _alpha[idx_cell][idx_sys];
}
}
break;
case ANALYTIC:
// Compute total "mass" of the system ==> to check conservation properties
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
double time = idx_pseudoTime * _dE;
double sigma = 0;
double time = idx_pseudoTime * _dE;
_outputFields[idx_group][0][idx_cell] =
( 4 * M_PI ) * _problem->GetAnalyticalSolution(
_mesh->GetCellMidPoints()[idx_cell][0], _mesh->GetCellMidPoints()[idx_cell][1], time, sigma );
}
break;
_outputFields[idx_group][0][idx_cell] = _problem->GetAnalyticalSolution(
_mesh->GetCellMidPoints()[idx_cell][0], _mesh->GetCellMidPoints()[idx_cell][1], time, _sigmaS[idx_pseudoTime][idx_cell] );
}
break;
default: ErrorMessages::Error( "Volume Output Group not defined for MN Solver!", CURRENT_FUNCTION ); break;
}
default: ErrorMessages::Error( "Volume Output Group not defined for MN Solver!", CURRENT_FUNCTION ); break;
}
}
return mass;
......
code/src/solvers/pnsolver.cpp
View file @ 5388b5cc
......@@ -381,24 +381,21 @@ double PNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
mass += _sol[idx_cell][0] * _areas[idx_cell]; // Should probably go to postprocessing
}
if( _settings->GetOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetOutputFrequency() == 0 ) {
for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
switch( _settings->GetVolumeOutput()[idx_group] ) {
case MINIMAL:
for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) {
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];
}
break;
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][0][idx_cell] = _sol[idx_cell][0];
}
break;
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] = firstMomentScaleFactor * _sol[idx_cell][idx_sys];
}
_outputFields[idx_group][idx_sys][idx_cell] = firstMomentScaleFactor * _sol[idx_cell][idx_sys];
}
break;
default: ErrorMessages::Error( "Volume Output Group not defined for PN Solver!", CURRENT_FUNCTION ); break;
}
}
break;
default: ErrorMessages::Error( "Volume Output Group not defined for PN Solver!", CURRENT_FUNCTION ); break;
}
}
return 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

Imprint