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

analytical solution works for gamma=1


Former-commit-id: 5388b5cc
parent df903b84
...@@ -18,7 +18,7 @@ MESH_FILE = linesource.su2 ...@@ -18,7 +18,7 @@ MESH_FILE = linesource.su2
%MESH_FILE = linesource_debug.su2 %MESH_FILE = linesource_debug.su2
% %
PROBLEM = LINESOURCE PROBLEM = LINESOURCE
SCATTER_COEFF = 0 SCATTER_COEFF = 1
% ---- Solver specifications ---- % ---- Solver specifications ----
% %
...@@ -29,7 +29,7 @@ CFL_NUMBER = 0.7 ...@@ -29,7 +29,7 @@ CFL_NUMBER = 0.7
% Final time for simulation % Final time for simulation
TIME_FINAL = 0.3 TIME_FINAL = 0.3
% Maximal Moment degree % Maximal Moment degree
MAX_MOMENT_SOLVER = 2 MAX_MOMENT_SOLVER = 3
% %
%% Entropy settings %% Entropy settings
ENTROPY_FUNCTIONAL = MAXWELL_BOLZMANN ENTROPY_FUNCTIONAL = MAXWELL_BOLZMANN
...@@ -62,3 +62,4 @@ QUAD_ORDER = 8 ...@@ -62,3 +62,4 @@ QUAD_ORDER = 8
% ----- Output ---- % ----- Output ----
% %
VOLUME_OUTPUT = (ANALYTIC, MINIMAL, MOMENTS, DUAL_MOMENTS) VOLUME_OUTPUT = (ANALYTIC, MINIMAL, MOMENTS, DUAL_MOMENTS)
OUTPUT_FREQUENCY = 1
...@@ -19,7 +19,7 @@ MESH_FILE = linesource.su2 ...@@ -19,7 +19,7 @@ MESH_FILE = linesource.su2
% %
PROBLEM = LINESOURCE PROBLEM = LINESOURCE
SCATTER_COEFF = 0 SCATTER_COEFF = 1
% %
% %
% ---- Solver specifications ---- % ---- Solver specifications ----
......
...@@ -17,7 +17,7 @@ LOG_DIR = ../result/logs ...@@ -17,7 +17,7 @@ LOG_DIR = ../result/logs
MESH_FILE = linesource.su2 MESH_FILE = linesource.su2
% %
PROBLEM = LINESOURCE PROBLEM = LINESOURCE
SCATTER_COEFF = 0 SCATTER_COEFF = 1
% %
% ---- Solver specifications ---- % ---- Solver specifications ----
......
...@@ -18,19 +18,24 @@ double LineSource::GetAnalyticalSolution( double x, double y, double t, double s ...@@ -18,19 +18,24 @@ double LineSource::GetAnalyticalSolution( double x, double y, double t, double s
double solution = 0.0; double solution = 0.0;
double R = sqrt( x * x + y * y ); double R = sqrt( x * x + y * y );
if( _sigmaS == 0.0 ) { if( t > 0 ) {
if( ( t - R ) > 0 ) { if( _sigmaS == 0.0 ) {
solution = 1 / ( 2 * M_PI * t * sqrt( t * t - R * R ) ); if( ( t - R ) > 0 ) {
solution = 1 / ( 2 * M_PI * t * sqrt( t * t - R * R ) );
}
} }
} else {
else { double gamma = R / t;
double gamma = R / t;
if( ( 1 - gamma ) > 0 ) { if( ( 1 - gamma ) > 0 ) {
solution = exp( -t ) / ( 2 * M_PI * t * t * sqrt( 1 - gamma * gamma ) ) + 2 * t * HelperIntRho_ptc( R, t ); 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 ) { double LineSource::HelperIntRho_ptc( double R, double t ) {
...@@ -42,10 +47,10 @@ 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 ) // integral is from 0 to sqrt( 1 - gamma * gamma )
double stepsize = sqrt( 1 - gamma * gamma ) / (double)numsteps; double stepsize = sqrt( 1 - gamma * gamma ) / (double)numsteps;
omega = 0.5 * stepsize;
for( int i = 0; i < numsteps; i++ ) { for( int i = 0; i < numsteps; i++ ) {
omega = omega + i * stepsize; omega = i * stepsize + 0.5 * stepsize;
integral += HelperRho_ptc( t * sqrt( gamma * gamma + omega * omega ), t );
integral += stepsize * HelperRho_ptc( t * sqrt( gamma * gamma + omega * omega ), t );
} }
return integral; return integral;
} }
...@@ -68,6 +73,10 @@ double LineSource::HelperRho_ptc2( double R, double t ) { ...@@ -68,6 +73,10 @@ double LineSource::HelperRho_ptc2( double R, double t ) {
// Compute the integralpart with midpoint rule // Compute the integralpart with midpoint rule
result = exp( -t ) / ( 32 * M_PI * M_PI * R ) * ( 1 - gamma * gamma ) * HelperIntRho_ptc2( t, gamma ); 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; return result;
} }
...@@ -86,17 +95,15 @@ double LineSource::HelperIntRho_ptc2( double t, double gamma ) { ...@@ -86,17 +95,15 @@ double LineSource::HelperIntRho_ptc2( double t, double gamma ) {
double stepsize = M_PI / (double)numsteps; double stepsize = M_PI / (double)numsteps;
u = 0.5 * stepsize;
q = ( 1.0 + gamma ) / ( 1.0 - gamma ); q = ( 1.0 + gamma ) / ( 1.0 - gamma );
for( int i = 0; i < numsteps; i++ ) { for( int i = 0; i < numsteps; i++ ) {
u = u + (double)i * stepsize; u = i * stepsize + 0.5 * stepsize;
// function evaluation // function evaluation
beta = ( log( q ) + com_one * u ) / ( gamma + com_one * tan( 0.5 * u ) ); 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 ); 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; return integral;
} }
......
...@@ -264,44 +264,39 @@ double MNSolver::WriteOutputFields( unsigned idx_pseudoTime ) { ...@@ -264,44 +264,39 @@ double MNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
mass += _sol[idx_cell][0] * _areas[idx_cell]; // Should probably go to postprocessing 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] ) {
for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) { case MINIMAL:
switch( _settings->GetVolumeOutput()[idx_group] ) { for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
case MINIMAL: _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 ) { for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
_outputFields[idx_group][0][idx_cell] = firstMomentScaleFactor * _sol[idx_cell][0]; _outputFields[idx_group][idx_sys][idx_cell] = _sol[idx_cell][idx_sys];
}
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];
}
} }
break; }
case ANALYTIC: break;
// Compute total "mass" of the system ==> to check conservation properties case DUAL_MOMENTS:
for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) { 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 time = idx_pseudoTime * _dE;
double sigma = 0;
_outputFields[idx_group][0][idx_cell] = _outputFields[idx_group][0][idx_cell] = _problem->GetAnalyticalSolution(
( 4 * M_PI ) * _problem->GetAnalyticalSolution( _mesh->GetCellMidPoints()[idx_cell][0], _mesh->GetCellMidPoints()[idx_cell][1], time, _sigmaS[idx_pseudoTime][idx_cell] );
_mesh->GetCellMidPoints()[idx_cell][0], _mesh->GetCellMidPoints()[idx_cell][1], time, sigma ); }
} break;
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; return mass;
......
...@@ -381,24 +381,21 @@ double PNSolver::WriteOutputFields( unsigned idx_pseudoTime ) { ...@@ -381,24 +381,21 @@ double PNSolver::WriteOutputFields( unsigned idx_pseudoTime ) {
mass += _sol[idx_cell][0] * _areas[idx_cell]; // Should probably go to postprocessing 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] ) {
for( unsigned idx_group = 0; idx_group < nGroups; idx_group++ ) { case MINIMAL:
switch( _settings->GetVolumeOutput()[idx_group] ) { for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
case MINIMAL: _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 ) { for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
_outputFields[idx_group][0][idx_cell] = _sol[idx_cell][0]; _outputFields[idx_group][idx_sys][idx_cell] = firstMomentScaleFactor * _sol[idx_cell][idx_sys];
}
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];
}
} }
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; return mass;
......
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