Skip to content

GitLab

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

Finish PN limiter 


Former-commit-id: 4986f7b5
parent 9168817a
Hide whitespace changes
Inline Side-by-side
code/src/common/mesh.cpp
View file @ a073d4d9
...@@ -358,6 +358,9 @@ void Mesh::ComputePartitioning() { ...@@ -358,6 +358,9 @@ void Mesh::ComputePartitioning() {
void Mesh::ComputeSlopes( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const { void Mesh::ComputeSlopes( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const {
for( unsigned k = 0; k < nq; ++k ) { for( unsigned k = 0; k < nq; ++k ) {
for( unsigned j = 0; j < _numCells; ++j ) { 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( cell->IsBoundaryCell() ) continue; // skip ghost cells
if( _cellBoundaryTypes[j] != 2 ) continue; // skip ghost cells if( _cellBoundaryTypes[j] != 2 ) continue; // skip ghost cells
// compute derivative by summing over cell boundary // compute derivative by summing over cell boundary
...@@ -400,8 +403,11 @@ void Mesh::ReconstructSlopesS( unsigned nq, VectorVector& psiDerX, VectorVector& ...@@ -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 { void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const {
double phi; double phi;
VectorVector dPsiMax = 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 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>> 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 ) ) ); 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& ...@@ -415,6 +421,7 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
if( _cellBoundaryTypes[j] != 2 ) continue; if( _cellBoundaryTypes[j] != 2 ) continue;
/* /*
// version 1: original taste
// step 1: calculate psi difference around neighbors and theoretical derivatives by Gauss theorem // step 1: calculate psi difference around neighbors and theoretical derivatives by Gauss theorem
for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) { 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] > 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& ...@@ -427,7 +434,8 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) { for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
// step 2: choose sample points // 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] = 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 psiDerY[j][k] * ( _nodes[_cells[j][l]][1] - _cellMidPoints[j][1] ); // vertex points
// step 3: calculate Phi_ij at sample points // step 3: calculate Phi_ij at sample points
...@@ -446,13 +454,13 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector& ...@@ -446,13 +454,13 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
phi = min( phiSample[j][k] ); phi = min( phiSample[j][k] );
// step 5: limit the slope reconstructed from Gauss theorem // step 5: limit the slope reconstructed from Gauss theorem
for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) { psiDerX[j][k] *= phi;
psiDerX[j][k] *= phi; psiDerY[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 ) { 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] > 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]; 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& ...@@ -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]; 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]; 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 ) { 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]) + // step 2: choose sample points
0.5 * psiDerY[j][k] * (_cellMidPoints[_cellNeighbors[j][l]][1] - _cellMidPoints[j][1]); psiSample[j][k][l] = 10.0 * psiDerX[j][k] * (_cellMidPoints[_cellNeighbors[j][l]][0] - _cellMidPoints[j][0]) +
//psiSample[j][k][l] = 3.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]);
// 3.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 ) { if( psiSample[j][k][l] > 0.0 ) {
phiSample[j][k][l] = 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] + 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] + 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] + eps);
} }
else if( psiSample[j][l][k] < 0.0 ) { else if( psiSample[j][k][l] < 0.0 ) {
phiSample[j][k][l] = 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] + 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] + 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] + eps);;
} }
else { else {
phiSample[j][k][l] = 1.0; phiSample[j][k][l] = 1.0;
...@@ -494,14 +492,11 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector& ...@@ -494,14 +492,11 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
// step 4: find minimum limiter function phi // step 4: find minimum limiter function phi
phi = min( phiSample[j][k] ); 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 // step 5: limit the slope reconstructed from Gauss theorem
for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) { psiDerX[j][k] *= phi;
psiDerX[j][k] *= phi; psiDerY[j][k] *= phi;
psiDerY[j][k] *= phi;
}
} }
} }
......
code/src/solvers/pnsolver.cpp
View file @ a073d4d9
...@@ -77,9 +77,12 @@ void PNSolver::ComputeRadFlux() { ...@@ -77,9 +77,12 @@ void PNSolver::ComputeRadFlux() {
void PNSolver::FluxUpdate() { void PNSolver::FluxUpdate() {
if( _reconsOrder > 1 ) { if( _reconsOrder > 1 ) {
_mesh->ReconstructSlopesU( _nTotalEntries, _solDx, _solDy, _sol ); // unstructured reconstruction _mesh->ReconstructSlopesU( _nTotalEntries, _solDx, _solDy, _sol ); // unstructured reconstruction
//_mesh->ComputeSlopes( _nTotalEntries, _solDx, _solDy, _sol ); // unstructured reconstruction
} }
Vector solL( _nTotalEntries ); //Vector solL( _nTotalEntries );
Vector solR( _nTotalEntries ); //Vector solR( _nTotalEntries );
auto solL = _sol[2];
auto solR = _sol[2];
// Loop over all spatial cells // Loop over all spatial cells
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) { for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
...@@ -100,26 +103,6 @@ void PNSolver::FluxUpdate() { ...@@ -100,26 +103,6 @@ void PNSolver::FluxUpdate() {
_solNew[idx_cell] += _g->Flux( _solNew[idx_cell] += _g->Flux(
_AxPlus, _AxMinus, _AyPlus, _AyMinus, _AzPlus, _AzMinus, _sol[idx_cell], _sol[idx_cell], _normals[idx_cell][idx_neighbor] ); _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,
solL,
solR,
_normals[idx_cell][idx_neighbor] );
switch( _reconsOrder ) { switch( _reconsOrder ) {
// first order solver // first order solver
case 1: case 1:
...@@ -139,16 +122,15 @@ void PNSolver::FluxUpdate() { ...@@ -139,16 +122,15 @@ void PNSolver::FluxUpdate() {
solL = _sol[idx_cell] + solL = _sol[idx_cell] +
_solDx[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_neighbor][0] - _cellMidPoints[idx_cell][0] ) + _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] ); _solDy[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_neighbor][1] - _cellMidPoints[idx_cell][1] );
// right status of interface // 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] ) + _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] ); _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) // positivity check (if not satisfied, deduce to first order)
if( min(solL) < 0.0 || min(solR) < 0.0 ) { //if( min(solL) < 0.0 || min(solR) < 0.0 ) {
solL = _sol[idx_cell]; // solL = _sol[idx_cell];
solR = _sol[_neighbors[idx_cell][idx_neighbor]]; // solR = _sol[_neighbors[idx_cell][idx_neighbor]];
} //}
// flux evaluation // flux evaluation
_solNew[idx_cell] += _g->Flux( _AxPlus, _solNew[idx_cell] += _g->Flux( _AxPlus,
_AxMinus, _AxMinus,
...@@ -159,6 +141,7 @@ void PNSolver::FluxUpdate() { ...@@ -159,6 +141,7 @@ void PNSolver::FluxUpdate() {
solL, solL,
solR, solR,
_normals[idx_cell][idx_neighbor] ); _normals[idx_cell][idx_neighbor] );
break;
// default: first order solver // default: first order solver
default: default:
_solNew[idx_cell] += _g->Flux( _AxPlus, _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

Imprint