Skip to content

GitLab

Commit 378f938c authored by Tianbai Xiao's avatar Tianbai Xiao
Browse files

Add limiter for PN

parent 1c262c4f
Hide whitespace changes
Inline Side-by-side
code/include/solvers/pnsolver.h
View file @ 378f938c
#ifndef PNSOLVER_H
#define PNSOLVER_H
#include "solverbase.h"
#include "solvers/solverbase.h"
class PNSolver : public Solver
{
......@@ -38,6 +38,9 @@ class PNSolver : public Solver
Vector _scatterMatDiag; /*! @brief: diagonal of the scattering matrix (its a diagonal matrix by construction). Contains eigenvalues of the
scattering kernel. */
VectorVector _solDx;
VectorVector _solDy;
// ---- Member functions ----
// IO
......
code/include/solvers/solverbase.h
View file @ 378f938c
......@@ -52,15 +52,14 @@ class Solver
/*! @brief edge neighbor cell ids, dim(_neighbors) = (_NCells,nEdgesPerCell) */
std::vector<std::vector<unsigned>> _neighbors;
// slope related params
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/src/common/mesh.cpp
View file @ 378f938c
......@@ -406,7 +406,6 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
std::vector<std::vector<Vector>> phiSample( _numCells, std::vector<Vector>( nq, Vector( _numNodesPerCell, 0.0 ) ) );
for( unsigned k = 0; k < nq; ++k ) {
for( unsigned j = 0; j < _numCells; ++j ) {
// reset derivatives
psiDerX[j][k] = 0.0;
......@@ -415,6 +414,7 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
// skip boundary cells
if( _cellBoundaryTypes[j] != 2 ) continue;
/*
// step 1: calculate psi difference around neighbors and theoretical derivatives by Gauss theorem
for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
if( psi[_cellNeighbors[j][l]][k] - psi[j][k] > dPsiMax[j][k] ) dPsiMax[j][k] = psi[_cellNeighbors[j][l]][k] - psi[j][k];
......@@ -432,26 +432,80 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
// step 3: calculate Phi_ij at sample points
if( psiSample[j][k][l] > psi[j][k] ) {
phiSample[j][k][l] = fmin( 0.5, dPsiMax[j][k] / ( psiSample[j][k][l] - psi[j][k] ) );
phiSample[j][k][l] = fmin( 1.0, dPsiMax[j][k] / ( psiSample[j][k][l] - psi[j][k] ) );
}
else if( psiSample[j][l][k] < psi[j][k] ) {
phiSample[j][k][l] = fmin( 0.5, dPsiMin[j][k] / ( psiSample[j][k][l] - psi[j][k] ) );
phiSample[j][k][l] = fmin( 1.0, dPsiMin[j][k] / ( psiSample[j][k][l] - psi[j][k] ) );
}
else {
phiSample[j][k][l] = 1.0;
}
}
// step 4: find minimum limiter function phi
phi = min( phiSample[j][k] );
// step 5: limit the slope reconstructed from Gauss theorem
for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
psiDerX[j][k] *= phi;
psiDerY[j][k] *= phi;
}
*/
// Venkatakrishnan limiter
for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
if( psi[_cellNeighbors[j][l]][k] - psi[j][k] > dPsiMax[j][k] ) dPsiMax[j][k] = psi[_cellNeighbors[j][l]][k] - psi[j][k];
if( psi[_cellNeighbors[j][l]][k] - psi[j][k] < dPsiMin[j][k] ) dPsiMin[j][k] = psi[_cellNeighbors[j][l]][k] - psi[j][k];
psiDerX[j][k] += 0.5 * ( psi[j][k] + psi[_cellNeighbors[j][l]][k] ) * _cellNormals[j][l][0] / _cellAreas[j];
psiDerY[j][k] += 0.5 * ( psi[j][k] + psi[_cellNeighbors[j][l]][k] ) * _cellNormals[j][l][1] / _cellAreas[j];
}
//psiSample[j][k][0] =
//psiDerX[j][k] * ( 0.5 * (_nodes[_cells[j][0]][0] + _nodes[_cells[j][1]][0]) - _cellMidPoints[j][0] ) +
//psiDerY[j][k] * ( 0.5 * (_nodes[_cells[j][0]][1] + _nodes[_cells[j][1]][1]) - _cellMidPoints[j][1] );
//psiSample[j][k][1] =
//psiDerX[j][k] * ( 0.5 * (_nodes[_cells[j][1]][0] + _nodes[_cells[j][2]][0]) - _cellMidPoints[j][0] ) +
//psiDerY[j][k] * ( 0.5 * (_nodes[_cells[j][1]][1] + _nodes[_cells[j][2]][1]) - _cellMidPoints[j][1] );
//psiSample[j][k][2] =
//psiDerX[j][k] * ( 0.5 * (_nodes[_cells[j][2]][0] + _nodes[_cells[j][0]][0]) - _cellMidPoints[j][0] ) +
//psiDerY[j][k] * ( 0.5 * (_nodes[_cells[j][2]][1] + _nodes[_cells[j][0]][1]) - _cellMidPoints[j][1] );
for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
psiSample[j][k][l] = 0.5 * psiDerX[j][k] * (_cellMidPoints[_cellNeighbors[j][l]][0] - _cellMidPoints[j][0]) +
0.5 * psiDerY[j][k] * (_cellMidPoints[_cellNeighbors[j][l]][1] - _cellMidPoints[j][1]);
//psiSample[j][k][l] = 3.0 * psiDerX[j][k] * (_cellMidPoints[_cellNeighbors[j][l]][0] - _cellMidPoints[j][0]) +
// 3.0 * psiDerY[j][k] * (_cellMidPoints[_cellNeighbors[j][l]][1] - _cellMidPoints[j][1]);
if( psiSample[j][k][l] > 0.0 ) {
phiSample[j][k][l] =
(dPsiMax[j][k]*dPsiMax[j][k] + 2.0 * dPsiMax[j][k] * psiSample[j][k][l] + 1e-6)/
(dPsiMax[j][k]*dPsiMax[j][k] + dPsiMax[j][k] * psiSample[j][k][l] + 2.0 * psiSample[j][k][l] * psiSample[j][k][l] + 1e-6);
}
else if( psiSample[j][l][k] < 0.0 ) {
phiSample[j][k][l] =
(dPsiMin[j][k]*dPsiMin[j][k] + 2.0 * dPsiMin[j][k] * psiSample[j][k][l] + 1e-6)/
(dPsiMin[j][k]*dPsiMin[j][k] + dPsiMin[j][k] * psiSample[j][k][l] + 2.0 * psiSample[j][k][l] * psiSample[j][k][l] + 1e-6);;
}
else {
phiSample[j][k][l] = 0.5;
phiSample[j][k][l] = 1.0;
}
}
// step 4: find minimum limiter function phi
phi = min( phiSample[j][k] );
phi = fmin( 1.0, phi );
// step 5: limit the slope reconstructed from Gauss theorem
for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
psiDerX[j][k] *= phi;
psiDerY[j][k] *= phi;
}
}
}
}
void Mesh::ComputeBounds() {
......
code/src/solvers/pnsolver.cpp
View file @ 378f938c
#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,9 @@ PNSolver::PNSolver( Config* settings ) : Solver( settings ) {
// Initialize Scatter Matrix
_scatterMatDiag = Vector( _nTotalEntries, 0 );
_solDx = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
_solDy = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
// Fill System Matrices
ComputeSystemMatrices();
......@@ -71,6 +75,12 @@ void PNSolver::ComputeRadFlux() {
}
void PNSolver::FluxUpdate() {
if( _reconsOrder > 1 ) {
_mesh->ReconstructSlopesU( _nTotalEntries, _solDx, _solDy, _sol ); // unstructured reconstruction
}
Vector solL( _nTotalEntries );
Vector solR( _nTotalEntries );
// Loop over all spatial cells
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
......@@ -89,16 +99,79 @@ 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 {
// 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] );
_solNew[idx_cell] += _g->Flux( _AxPlus,
_AxMinus,
_AyPlus,
_AyMinus,
_AzPlus,
_AzMinus,
_sol[idx_cell],
_sol[_neighbors[idx_cell][idx_neighbor]],
solL,
solR,
_normals[idx_cell][idx_neighbor] );
switch( _reconsOrder ) {
// first order solver
case 1:
_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] );
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 check (if not satisfied, deduce to first order)
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] );
// 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/solverbase.cpp
View file @ 378f938c
......@@ -32,7 +32,7 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
// 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 ) ) );
......
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