Commit fe630600 authored by yy3406's avatar yy3406
Browse files

Merge branch 'recons' into 'master'

Merge recons

See merge request !34

Former-commit-id: 66b4e207
parents ec87706a 4b488032
#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
......
......@@ -14,6 +14,7 @@ class Mesh;
class Config;
class ProblemBase;
class QuadratureBase;
class Reconstructor;
class Solver
{
......@@ -51,6 +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 */
......
......@@ -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
......
......@@ -32,6 +32,8 @@ CFL_NUMBER = 0.7
TIME_FINAL = 0.3
% Maximal Moment degree
MAX_MOMENT_SOLVER = 2
% Reconstruction order
RECONS_ORDER = 2
% ---- Boundary Conditions ----
% Example: BC_DIRICLET = (dummyMarker1, dummyMarker2)
......
......@@ -28,6 +28,8 @@ CONTINUOUS_SLOWING_DOWN = YES
CFL_NUMBER = 0.9
% Final time for simulation
TIME_FINAL = 0.5
% Reconstruction order
RECONS_ORDER = 2
%
CLEAN_FLUX_MATRICES = NO
......
......@@ -358,6 +358,9 @@ void Mesh::ComputePartitioning() {
void Mesh::ComputeSlopes( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const {
for( unsigned k = 0; k < nq; ++k ) {
for( unsigned j = 0; j < _numCells; ++j ) {
psiDerX[j][k] = 0.0;
psiDerY[j][k] = 0.0;
// if( cell->IsBoundaryCell() ) continue; // skip ghost cells
if( _cellBoundaryTypes[j] != 2 ) continue; // skip ghost cells
// compute derivative by summing over cell boundary
......@@ -400,13 +403,15 @@ void Mesh::ReconstructSlopesS( unsigned nq, VectorVector& psiDerX, VectorVector&
void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const {
double phi;
VectorVector dPsiMax = std::vector( _numCells, Vector( nq, 0.0 ) );
VectorVector dPsiMin = std::vector( _numCells, Vector( nq, 0.0 ) );
//VectorVector dPsiMax = std::vector( _numCells, Vector( nq, 0.0 ) );
//VectorVector dPsiMin = std::vector( _numCells, Vector( nq, 0.0 ) );
VectorVector dPsiMax( _numCells, Vector( nq, 0.0 ) );
VectorVector dPsiMin( _numCells, Vector( nq, 0.0 ) );
std::vector<std::vector<Vector>> psiSample( _numCells, std::vector<Vector>( nq, Vector( _numNodesPerCell, 0.0 ) ) );
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 +420,8 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
// skip boundary cells
if( _cellBoundaryTypes[j] != 2 ) continue;
/*
// version 1: original taste
// 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];
......@@ -427,18 +434,19 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
// step 2: choose sample points
// psiSample[j][k][l] = 0.5 * ( psi[j][k] + psi[_cellNeighbors[j][l]][k] ); // interface central points
psiSample[j][k][l] = psi[j][k] + psiDerX[j][k] * ( _nodes[_cells[j][l]][0] - _cellMidPoints[j][0] ) +
psiSample[j][k][l] = psi[j][k] +
psiDerX[j][k] * ( _nodes[_cells[j][l]][0] - _cellMidPoints[j][0] ) +
psiDerY[j][k] * ( _nodes[_cells[j][l]][1] - _cellMidPoints[j][1] ); // vertex points
// 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] = 0.5;
phiSample[j][k][l] = 1.0;
}
}
......@@ -446,12 +454,53 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
phi = min( phiSample[j][k] );
// step 5: limit the slope reconstructed from Gauss theorem
psiDerX[j][k] *= phi;
psiDerY[j][k] *= phi;
*/
double eps = 1e-6;
// version 2: Venkatakrishnan limiter
// 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];
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];
}
for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
psiDerX[j][k] *= phi;
psiDerY[j][k] *= phi;
// step 2: choose sample points
psiSample[j][k][l] = 10.0 * psiDerX[j][k] * (_cellMidPoints[_cellNeighbors[j][l]][0] - _cellMidPoints[j][0]) +
10.0 * psiDerY[j][k] * (_cellMidPoints[_cellNeighbors[j][l]][1] - _cellMidPoints[j][1]);
// step 3: calculate Phi_ij at sample points
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] + eps)/
(dPsiMax[j][k]*dPsiMax[j][k] + dPsiMax[j][k] * psiSample[j][k][l] + 2.0 * psiSample[j][k][l] * psiSample[j][k][l] + eps);
}
else if( psiSample[j][k][l] < 0.0 ) {
phiSample[j][k][l] =
(dPsiMin[j][k]*dPsiMin[j][k] + 2.0 * dPsiMin[j][k] * psiSample[j][k][l] + eps)/
(dPsiMin[j][k]*dPsiMin[j][k] + dPsiMin[j][k] * psiSample[j][k][l] + 2.0 * psiSample[j][k][l] * psiSample[j][k][l] + eps);;
}
else {
phiSample[j][k][l] = 1.0;
}
}
// step 4: find minimum limiter function phi
phi = min( phiSample[j][k] );
//phi = fmin( 0.5, abs(phi) );
// step 5: limit the slope reconstructed from Gauss theorem
psiDerX[j][k] *= phi;
psiDerY[j][k] *= phi;
}
}
}
void Mesh::ComputeBounds() {
......
#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,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,16 +102,59 @@ 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
_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] );
else {
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] );
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] );
}
}
}
}
}
......
......@@ -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] );
}
}
}
}
......
......@@ -29,6 +29,28 @@ Solver::Solver( Config* settings ) : _settings( 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 );
......
#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