Commit 4986f7b5 authored by Tianbai Xiao's avatar Tianbai Xiao
Browse files

Finish PN limiter

parent 378f938c
......@@ -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,8 +403,11 @@ 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 ) ) );
......@@ -415,6 +421,7 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
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,7 +434,8 @@ 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
......@@ -446,13 +454,13 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
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;
}
psiDerX[j][k] *= phi;
psiDerY[j][k] *= phi;
*/
// Venkatakrishnan limiter
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];
......@@ -460,32 +468,22 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
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]);
// 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] + 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);
(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][l][k] < 0.0 ) {
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] + 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);;
(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;
......@@ -494,14 +492,11 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
// step 4: find minimum limiter function phi
phi = min( phiSample[j][k] );
phi = fmin( 1.0, phi );
//phi = fmin( 0.5, abs(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;
}
psiDerX[j][k] *= phi;
psiDerY[j][k] *= phi;
}
}
......
......@@ -77,9 +77,12 @@ 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 );
//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++ ) {
......@@ -100,26 +103,6 @@ void PNSolver::FluxUpdate() {
_solNew[idx_cell] += _g->Flux(
_AxPlus, _AxMinus, _AyPlus, _AyMinus, _AzPlus, _AzMinus, _sol[idx_cell], _sol[idx_cell], _normals[idx_cell][idx_neighbor] );
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,
solL,
solR,
_normals[idx_cell][idx_neighbor] );
switch( _reconsOrder ) {
// first order solver
case 1:
......@@ -139,16 +122,15 @@ void PNSolver::FluxUpdate() {
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]];// +
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]];
}
//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,
......@@ -159,6 +141,7 @@ void PNSolver::FluxUpdate() {
solL,
solR,
_normals[idx_cell][idx_neighbor] );
break;
// default: first order solver
default:
_solNew[idx_cell] += _g->Flux( _AxPlus,
......
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