Skip to content

GitLab

Commit 96ad94b5 authored by jannick.wolters's avatar jannick.wolters
Browse files

merged current master 


Former-commit-id: a54b06bb
parents 2360c430 7590242d
Expand all Show whitespace changes
Inline Side-by-side
code/src/solvers/pnsolver.cpp
View file @ 96ad94b5
#include "solvers/pnsolver.h"
#include "common/config.h"
#include "common/io.h"
#include "common/mesh.h"
#include "fluxes/numericalflux.h"
#include "toolboxes/errormessages.h"
#include "toolboxes/textprocessingtoolbox.h"
......@@ -31,6 +32,10 @@ PNSolver::PNSolver( Config* settings ) : Solver( settings ) {
// Initialize Scatter Matrix
_scatterMatDiag = Vector( _nTotalEntries, 0 );
// Initialize temporary storages of solution derivatives
_solDx = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
_solDy = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
// Fill System Matrices
ComputeSystemMatrices();
......@@ -71,6 +76,15 @@ void PNSolver::ComputeRadFlux() {
}
void PNSolver::FluxUpdate() {
if( _reconsOrder > 1 ) {
_mesh->ReconstructSlopesU( _nTotalEntries, _solDx, _solDy, _sol ); // unstructured reconstruction
//_mesh->ComputeSlopes( _nTotalEntries, _solDx, _solDy, _sol ); // unstructured reconstruction
}
//Vector solL( _nTotalEntries );
//Vector solR( _nTotalEntries );
auto solL = _sol[2];
auto solR = _sol[2];
// Loop over all spatial cells
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
......@@ -89,7 +103,10 @@ void PNSolver::FluxUpdate() {
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_neighbor] == _nCells )
_solNew[idx_cell] += _g->Flux(
_AxPlus, _AxMinus, _AyPlus, _AyMinus, _AzPlus, _AzMinus, _sol[idx_cell], _sol[idx_cell], _normals[idx_cell][idx_neighbor] );
else
else {
switch( _reconsOrder ) {
// first order solver
case 1:
_solNew[idx_cell] += _g->Flux( _AxPlus,
_AxMinus,
_AyPlus,
......@@ -99,6 +116,47 @@ void PNSolver::FluxUpdate() {
_sol[idx_cell],
_sol[_neighbors[idx_cell][idx_neighbor]],
_normals[idx_cell][idx_neighbor] );
break;
// second order solver
case 2:
// left status of interface
solL = _sol[idx_cell] +
_solDx[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_neighbor][0] - _cellMidPoints[idx_cell][0] ) +
_solDy[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_neighbor][1] - _cellMidPoints[idx_cell][1] );
// right status of interface
solR = _sol[_neighbors[idx_cell][idx_neighbor]] +
_solDx[_neighbors[idx_cell][idx_neighbor]] * ( _interfaceMidPoints[idx_cell][idx_neighbor][0] - _cellMidPoints[_neighbors[idx_cell][idx_neighbor]][0] ) +
_solDy[_neighbors[idx_cell][idx_neighbor]] * ( _interfaceMidPoints[idx_cell][idx_neighbor][1] - _cellMidPoints[_neighbors[idx_cell][idx_neighbor]][1] );
// positivity checker (if not satisfied, deduce to first order)
// I manually turned it off here since Pn produces negative solutions essentially
//if( min(solL) < 0.0 || min(solR) < 0.0 ) {
// solL = _sol[idx_cell];
// solR = _sol[_neighbors[idx_cell][idx_neighbor]];
//}
// flux evaluation
_solNew[idx_cell] += _g->Flux( _AxPlus,
_AxMinus,
_AyPlus,
_AyMinus,
_AzPlus,
_AzMinus,
solL,
solR,
_normals[idx_cell][idx_neighbor] );
break;
// default: first order solver
default:
_solNew[idx_cell] += _g->Flux( _AxPlus,
_AxMinus,
_AyPlus,
_AyMinus,
_AzPlus,
_AzMinus,
_sol[idx_cell],
_sol[_neighbors[idx_cell][idx_neighbor]],
_normals[idx_cell][idx_neighbor] );
}
}
}
}
}
......
code/src/solvers/snsolver.cpp
View file @ 96ad94b5
......@@ -55,7 +55,6 @@ void SNSolver::FluxUpdate() {
// left and right angular flux of interface, used in numerical flux evaluation
double psiL;
double psiR;
double mass;
// derivatives of angular flux in x and y directions
VectorVector psiDx( _nCells, Vector( _nq, 0.0 ) );
......@@ -98,7 +97,12 @@ void SNSolver::FluxUpdate() {
//_mesh->ReconstructSlopesS( _nq, psiDx, psiDy, _psi ); // structured reconstruction (not stable currently)
}
*/
if( _reconsOrder > 1 ) {
_mesh->ReconstructSlopesU( _nq, _psiDx, _psiDy, _sol ); // unstructured reconstruction
//_mesh->ReconstructSlopesS( _nq, _psiDx, _psiDy, _psi ); // structured reconstruction (not stable currently)
}
double psiL;
double psiR;
// Loop over all spatial cells
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
// Dirichlet cells stay at IC, farfield assumption
......@@ -114,45 +118,40 @@ void SNSolver::FluxUpdate() {
_solNew[idx_cell][idx_quad] +=
_g->Flux( _quadPoints[idx_quad], _sol[idx_cell][idx_quad], _sol[idx_cell][idx_quad], _normals[idx_cell][idx_neighbor] );
else {
/* switch( reconsOrder ) {
switch( _reconsOrder ) {
// first order solver
case 1:
psiNew[idx_cells][idx_quad] += _g->Flux( _quadPoints[idx_quad],
_sol[idx_cells][idx_quad],
_sol[_neighbors[idx_cells][idx_neighbor]][idx_quad],
_normals[idx_cells][idx_neighbor] );
_solNew[idx_cell][idx_quad] += _g->Flux( _quadPoints[idx_quad],
_sol[idx_cell][idx_quad],
_sol[_neighbors[idx_cell][idx_neighbor]][idx_quad],
_normals[idx_cell][idx_neighbor] );
break;
// second order solver
case 2:
// left status of interface
psiL = _sol[idx_cells][idx_quad] +
psiDx[idx_cells][idx_quad] * ( interfaceMidPoints[idx_cells][idx_neighbor][0] -
cellMidPoints[idx_cells][0] ) + psiDy[idx_cells][idx_quad] * ( interfaceMidPoints[idx_cells][idx_neighbor][1] -
cellMidPoints[idx_cells][1] );
psiL = _sol[idx_cell][idx_quad] +
_psiDx[idx_cell][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_neighbor][0] - _cellMidPoints[idx_cell][0] ) +
_psiDy[idx_cell][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_neighbor][1] - _cellMidPoints[idx_cell][1] );
// right status of interface
psiR = _sol[_neighbors[idx_cells][idx_neighbor]][idx_quad] +
psiDx[_neighbors[idx_cells][idx_neighbor]][idx_quad] *
( interfaceMidPoints[idx_cells][idx_neighbor][0] -
cellMidPoints[_neighbors[idx_cells][idx_neighbor]][0] ) + psiDy[_neighbors[idx_cells][idx_neighbor]][idx_quad] * (
interfaceMidPoints[idx_cells][idx_neighbor][1] - cellMidPoints[_neighbors[idx_cells][idx_neighbor]][1] );
psiR = _sol[_neighbors[idx_cell][idx_neighbor]][idx_quad] +
_psiDx[_neighbors[idx_cell][idx_neighbor]][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_neighbor][0] - _cellMidPoints[_neighbors[idx_cell][idx_neighbor]][0] ) +
_psiDy[_neighbors[idx_cell][idx_neighbor]][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_neighbor][1] - _cellMidPoints[_neighbors[idx_cell][idx_neighbor]][1] );
// positivity check (if not satisfied, deduce to first order)
if( psiL < 0.0 || psiR < 0.0 ) {
psiL = _sol[idx_cells][idx_quad];
psiR = _sol[_neighbors[idx_cells][idx_neighbor]][idx_quad];
psiL = _sol[idx_cell][idx_quad];
psiR = _sol[_neighbors[idx_cell][idx_neighbor]][idx_quad];
}
// flux evaluation
psiNew[idx_cells][idx_quad] += _g->Flux( _quadPoints[idx_quad], psiL, psiR,
_normals[idx_cells][idx_neighbor] ); break;
_solNew[idx_cell][idx_quad] += _g->Flux( _quadPoints[idx_quad], psiL, psiR, _normals[idx_cell][idx_neighbor] ); break;
// higher order solver
case 3: std::cout << "higher order is WIP" << std::endl; break;
// default: first order solver
default:
*/
_solNew[idx_cell][idx_quad] += _g->Flux( _quadPoints[idx_quad],
_sol[idx_cell][idx_quad],
_sol[_neighbors[idx_cell][idx_neighbor]][idx_quad],
_normals[idx_cell][idx_neighbor] );
// }
}
}
}
}
......
code/src/solvers/solverbase.cpp
View file @ 96ad94b5
......@@ -25,10 +25,32 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
_settings->SetNCells( _nCells );
// build quadrature object and store frequently used params
_quadrature = QuadratureBase::CreateQuadrature( settings );
_quadrature = QuadratureBase::Create( settings );
_nq = _quadrature->GetNq();
_settings->SetNQuadPoints( _nq );
// build slope related params
_reconstructor = Reconstructor::Create( settings );
_reconsOrder = _reconstructor->GetReconsOrder();
auto nodes = _mesh->GetNodes();
auto cells = _mesh->GetCells();
std::vector<std::vector<Vector>> interfaceMidPoints( _nCells, std::vector<Vector>( _mesh->GetNumNodesPerCell(), Vector( 2, 1e-10 ) ) );
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
for( unsigned k = 0; k < _mesh->GetDim(); ++k ) {
for( unsigned j = 0; j < _neighbors[idx_cell].size() - 1; ++j ) {
interfaceMidPoints[idx_cell][j][k] = 0.5 * ( nodes[cells[idx_cell][j]][k] + nodes[cells[idx_cell][j + 1]][k] );
}
interfaceMidPoints[idx_cell][_neighbors[idx_cell].size() - 1][k] =
0.5 * ( nodes[cells[idx_cell][_neighbors[idx_cell].size() - 1]][k] + nodes[cells[idx_cell][0]][k] );
}
}
_interfaceMidPoints = interfaceMidPoints;
_cellMidPoints = _mesh->GetCellMidPoints();
_psiDx = VectorVector( _nCells, Vector( _nq, 0.0 ) );
_psiDy = VectorVector( _nCells, Vector( _nq, 0.0 ) );
// set time step
_dE = ComputeTimeStep( settings->GetCFL() );
_nEnergies = unsigned( settings->GetTEnd() / _dE );
......
code/src/toolboxes/interpolation.cpp
View file @ 96ad94b5
......@@ -91,19 +91,19 @@ double Interpolation::operator()( double x, double y ) const {
if( _type == cubic ) {
// find closest values to x and y in table (lower bounds)
unsigned xId = IndexOfClosestValue( x, _x );
unsigned yId = IndexOfClosestValue( y, _y );
int xId = IndexOfClosestValue( x, _x );
int yId = IndexOfClosestValue( y, _y );
// store all 16 interpolation points needed
double points[4][4];
for( int i = -1; i < 3; ++i ) {
unsigned idx_y;
idx_y = yId + i < 0 ? 0 : yId + i;
idx_y = yId + i > _data.rows() - 1 ? _data.rows() - 1 : yId + i;
idx_y = yId + i > static_cast<int>( _data.rows() - 1 ) ? _data.rows() - 1 : yId + i;
for( int j = -1; j < 3; ++j ) {
unsigned idx_x;
idx_x = xId + j < 0 ? 0 : xId + j;
idx_x = xId + j > _data.columns() - 1 ? _data.columns() - 1 : xId + j;
idx_x = xId + j > static_cast<int>( _data.columns() - 1 ) ? _data.columns() - 1 : xId + j;
points[i + 1][j + 1] = _data( idx_x, idx_y );
}
......
code/src/toolboxes/reconstructor.cpp
View file @ 96ad94b5
#include "toolboxes/reconstructor.h"
#include "common/config.h"
Reconstructor::Reconstructor( Config* settings ) {}
Reconstructor::Reconstructor( Config* settings ) {
_reconsOrder = settings->GetReconsOrder();
}
Reconstructor* Reconstructor::Create( Config* settings ) {
return new Reconstructor( settings );
}
double FortSign( double a, double b ) {
if( b > 0.0 ) return abs( a );
if( b < 0.0 ) return -abs( a );
if( b > 0.0 ) return std::fabs( a );
if( b < 0.0 ) return -std::fabs( a );
return 0.0;
}
double LMinMod( double sL, double sR ) { return 0.5 * ( FortSign( 1.0, sL ) + FortSign( 1., sR ) ) * fmin( abs( sL ), abs( sR ) ); }
double LMinMod( double sL, double sR ) { return 0.5 * ( FortSign( 1.0, sL ) + FortSign( 1., sR ) ) * fmin( std::fabs( sL ), std::fabs( sR ) ); }
double LVanLeer( double sL, double sR ) {
return ( FortSign( 1.0, sL ) + FortSign( 1.0, sR ) ) * abs( sL ) * abs( sR ) / ( abs( sL ) + abs( sR ) + 0.0000001 );
return ( FortSign( 1.0, sL ) + FortSign( 1.0, sR ) ) * std::fabs( sL ) * std::fabs( sR ) / ( std::fabs( sL ) + std::fabs( sR ) + 0.0000001 );
}
double LSuperBee( double sL, double sR ) {
if( sR >= 0.5 * sL && sR <= 2.0 * sL ) return 0.5 * ( FortSign( 1.0, sL ) + FortSign( 1., sR ) ) * fmax( abs( sL ), abs( sR ) );
if( sR < 0.5 * sL && sR > 2.0 * sL ) return ( FortSign( 1.0, sL ) + FortSign( 1., sR ) ) * fmin( abs( sL ), abs( sR ) );
if( sR >= 0.5 * sL && sR <= 2.0 * sL ) return 0.5 * ( FortSign( 1.0, sL ) + FortSign( 1., sR ) ) * fmax( std::fabs( sL ), std::fabs( sR ) );
if( sR < 0.5 * sL && sR > 2.0 * sL ) return ( FortSign( 1.0, sL ) + FortSign( 1., sR ) ) * fmin( std::fabs( sL ), std::fabs( sR ) );
return 0.0;
}
double LVanAlbaba( double sL, double sR ) { return ( sL * sL * sR + sL * sR * sR ) / ( sL * sL + sR * sR + 0.0000001 ); }
double LWENOJS( double x ) { return 0.0; }
double LWENOJS( double /*x*/ ) { return 0.0; }
double Reconstructor::ReconstructSlopeStruct( double uL, double uC, double uR, double dxL, double dxR, std::string limiter ) const {
double sL = ( uC - uL ) / dxL;
......
code/src/solvers/sphericalharmonics.cpp → code/src/toolboxes/sphericalharmonics.cpp
View file @ 96ad94b5
#include "solvers/sphericalharmonics.h"
#include "toolboxes/sphericalharmonics.h"
#include <math.h>
SphericalHarmonics::SphericalHarmonics( unsigned L_degree ) {
......
code/tests/input/validation_tests/MN_solver/checkerboard_MN.cfg
View file @ 96ad94b5
......@@ -18,7 +18,7 @@ PROBLEM = CHECKERBOARD
% ---- Solver specifications ----
%
CFL_NUMBER = 0.5
TIME_FINAL = 3.2
TIME_FINAL = 0.4
SOLVER = MN_SOLVER
MAX_MOMENT_SOLVER = 2
ENTROPY_FUNCTIONAL = MAXWELL_BOLTZMANN
......
code/tests/input/validation_tests/PN_solver/checkerboard_PN.cfg
View file @ 96ad94b5
......@@ -18,7 +18,7 @@ PROBLEM = CHECKERBOARD
% ---- Solver specifications ----
%
CFL_NUMBER = 0.5
TIME_FINAL = 3.2
TIME_FINAL = 0.4
SOLVER = PN_SOLVER
MAX_MOMENT_SOLVER = 2
%
......
code/tests/input/validation_tests/SN_solver/checkerboard_SN.cfg
View file @ 96ad94b5
......@@ -18,7 +18,7 @@ PROBLEM = CHECKERBOARD
% ---- Solver specifications ----
%
CFL_NUMBER = 0.5
TIME_FINAL = 3.2
TIME_FINAL = 0.4
SOLVER = SN_SOLVER
%
% ---- Boundary Conditions ----
......
code/tests/test_harmonics.cpp
View file @ 96ad94b5
......@@ -2,7 +2,7 @@
#include "catch.hpp"
#include "common/config.h"
#include "quadratures/qgausslegendretensorized.h"
#include "solvers/sphericalharmonics.h"
#include "toolboxes/sphericalharmonics.h"
#include <fstream>
#include <iostream>
......
code/tests/test_kernel.cpp
View file @ 96ad94b5
......@@ -11,7 +11,7 @@ TEST_CASE( "test all scattering kernels", "[kernel]" ) {
// Load Settings from File
Config* config = new Config( filename );
QuadratureBase* quad = QuadratureBase::CreateQuadrature( config ); //@TODO: swap out for different quadrature rule
QuadratureBase* quad = QuadratureBase::Create( config ); //@TODO: swap out for different quadrature rule
SECTION( "isotropic scattering kernel" ) {
......
code/tests/test_optimizer.cpp
View file @ 96ad94b5
......@@ -4,7 +4,7 @@
#include "common/config.h"
#include "optimizers/optimizerbase.h"
#include "quadratures/quadraturebase.h"
#include "solvers/sphericalharmonics.h"
#include "toolboxes/sphericalharmonics.h"
TEST_CASE( "Test the Newton Optimizer", "[optimizers]" ) {
std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/optimizers/unit_optimizerNewton.cfg";
......@@ -16,7 +16,7 @@ TEST_CASE( "Test the Newton Optimizer", "[optimizers]" ) {
SphericalHarmonics basis( config->GetMaxMomentDegree() );
// Get Quadrature
QuadratureBase* quad = QuadratureBase::CreateQuadrature( config );
QuadratureBase* quad = QuadratureBase::Create( config );
// Get Optimizer (Newton)
OptimizerBase* optimizer = OptimizerBase::Create( config );
......
code/tests/test_quadrature.cpp
View file @ 96ad94b5
......@@ -52,7 +52,7 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) {
// Set quadOrder
config->SetQuadOrder( quadratureorder );
QuadratureBase* Q = QuadratureBase::CreateQuadrature( config );
QuadratureBase* Q = QuadratureBase::Create( config );
if( quadraturename == QUAD_GaussLegendre1D ) {
if( !approxequal( Q->SumUpWeights(), 2, lowAccuracyTesting ) ) {
......@@ -79,7 +79,7 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) {
// Special case for Gauss Legendre with half weights
if( quadraturename == QUAD_GaussLegendreTensorized ) {
config->SetSNAllGaussPts( false );
QuadratureBase* Q = QuadratureBase::CreateQuadrature( config );
QuadratureBase* Q = QuadratureBase::Create( config );
if( !approxequal( Q->SumUpWeights(), 4 * M_PI, lowAccuracyTesting ) ) {
printf( "Quadrature %d at order %d . Error : %.15f (low accuracy testing was set to %d). Reduced number of quadrature "
"points used. \n",
......@@ -116,7 +116,7 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) {
bool errorWithinBounds = true;
QuadratureBase* Q = QuadratureBase::CreateQuadrature( config );
QuadratureBase* Q = QuadratureBase::Create( config );
VectorVector points = Q->GetPoints();
for( unsigned i = 0; i < Q->GetNq(); i++ ) {
if( !approxequal( 1.0, norm( points[i] ), lowAccuracyTesting ) ) {
......@@ -134,7 +134,7 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) {
// Special case for Gauss Legendre with half weights
if( quadraturename == QUAD_GaussLegendreTensorized ) {
config->SetSNAllGaussPts( false );
QuadratureBase* Q = QuadratureBase::CreateQuadrature( config );
QuadratureBase* Q = QuadratureBase::Create( config );
VectorVector points = Q->GetPoints();
for( unsigned i = 0; i < Q->GetNq(); i++ ) {
......@@ -169,14 +169,14 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) {
// Set quadOrder
config->SetQuadOrder( quadratureorder );
QuadratureBase* Q = QuadratureBase::CreateQuadrature( config );
QuadratureBase* Q = QuadratureBase::Create( config );
REQUIRE( Q->GetNq() == size( Q->GetWeights() ) );
// Special case for Gauss Legendre with half weights
if( quadraturename == QUAD_GaussLegendreTensorized ) {
config->SetSNAllGaussPts( false );
QuadratureBase* Q = QuadratureBase::CreateQuadrature( config );
QuadratureBase* Q = QuadratureBase::Create( config );
REQUIRE( Q->GetNq() == size( Q->GetWeights() ) );
config->SetSNAllGaussPts( true );
}
......@@ -194,13 +194,13 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) {
// Set quadOrder
config->SetQuadOrder( quadratureorder );
QuadratureBase* Q = QuadratureBase::CreateQuadrature( config );
QuadratureBase* Q = QuadratureBase::Create( config );
REQUIRE( Q->GetNq() == size( Q->GetPoints() ) );
// Special case for Gauss Legendre with half weights
if( quadraturename == QUAD_GaussLegendreTensorized ) {
config->SetSNAllGaussPts( false );
QuadratureBase* Q = QuadratureBase::CreateQuadrature( config );
QuadratureBase* Q = QuadratureBase::Create( config );
REQUIRE( Q->GetNq() == size( Q->GetPoints() ) );
config->SetSNAllGaussPts( true );
}
......@@ -224,7 +224,7 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) {
// Set quadOrder
config->SetQuadOrder( quadratureorder );
QuadratureBase* Q = QuadratureBase::CreateQuadrature( config );
QuadratureBase* Q = QuadratureBase::Create( config );
if( quadraturename == QUAD_GaussLegendre1D ) {
if( !approxequal( Q->Integrate( sin ), 0, lowAccuracyTesting ) ) {
......@@ -252,7 +252,7 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) {
// Special case for Gauss Legendre with half weights
if( quadraturename == QUAD_GaussLegendreTensorized ) {
config->SetSNAllGaussPts( false );
QuadratureBase* Q = QuadratureBase::CreateQuadrature( config );
QuadratureBase* Q = QuadratureBase::Create( config );
if( !approxequal( Q->Integrate( f ), 4.0 * M_PI, lowAccuracyTesting ) ) {
printf(
......
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