Skip to content

GitLab

Commit c7af2e37 authored by Tianbai Xiao's avatar Tianbai Xiao
Browse files

Adapt 2nd-order scheme onto new solverbase 


Former-commit-id: 8c03cb95
parent ec87706a
Hide whitespace changes
Inline Side-by-side
code/include/solvers/solverbase.h
View file @ c7af2e37
......@@ -14,6 +14,7 @@ class Mesh;
class Config;
class ProblemBase;
class QuadratureBase;
class Reconstructor;
class Solver
{
......@@ -51,6 +52,15 @@ class Solver
/*! @brief edge neighbor cell ids, dim(_neighbors) = (_NCells,nEdgesPerCell) */
std::vector<std::vector<unsigned>> _neighbors;
Reconstructor* _reconstructor; /*! @brief reconstructor class for high order scheme */
unsigned _reconsOrder;
VectorVector _psiDx;
VectorVector _psiDy;
VectorVector _cellMidPoints;
std::vector<std::vector<Vector>> _interfaceMidPoints;
// Solution related members
VectorVector _sol; /*! @brief solution of the PDE, e.g. angular flux or moments */
std::vector<double> _solverOutput; /*! @brief LEGACY: Outputfield for solver ==> Will be replaced by _outputFields in the near future */
......
code/include/toolboxes/reconstructor.h
View file @ c7af2e37
......@@ -15,12 +15,18 @@ class Config;
class Reconstructor
{
protected:
unsigned _reconsOrder;
public:
/**
* @brief Reconstruction
* @param settings
*/
Reconstructor( Config* settings );
static Reconstructor* Create( Config* settings );
unsigned inline GetReconsOrder() { return _reconsOrder; }
/*! Method 1: structured developing
* @brief Slope of angular flux psi inside a given cell
......
code/src/solvers/snsolver.cpp
View file @ c7af2e37
......@@ -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 ) {
// 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] );
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] );
// 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] );
// 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];
}
// flux evaluation
psiNew[idx_cells][idx_quad] += _g->Flux( _quadPoints[idx_quad], psiL, psiR,
_normals[idx_cells][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] );
// }
switch( _reconsOrder ) {
// first order solver
case 1:
_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_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_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_cell][idx_quad];
psiR = _sol[_neighbors[idx_cell][idx_neighbor]][idx_quad];
}
// flux evaluation
_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 @ c7af2e37
......@@ -29,6 +29,30 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
_nq = _quadrature->GetNq();
_settings->SetNQuadPoints( _nq );
auto nodes = _mesh->GetNodes();
auto cells = _mesh->GetCells();
_reconstructor = Reconstructor::Create( settings );
_reconsOrder = _reconstructor->GetReconsOrder();
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();
VectorVector psiDx( _nCells, Vector( _nq, 0.0 ) );
VectorVector psiDy( _nCells, Vector( _nq, 0.0 ) );
_psiDx = psiDx;
_psiDy = psiDy;
// set time step
_dE = ComputeTimeStep( settings->GetCFL() );
_nEnergies = unsigned( settings->GetTEnd() / _dE );
......
code/src/toolboxes/reconstructor.cpp
View file @ c7af2e37
#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 );
......
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