Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

accelerated sn solver by 25% through memory access optimization in Fl… #45

Merged
merged 1 commit into from
Mar 1, 2024
Merged

accelerated sn solver by 25% through memory access optimization in Fl… #45

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 12 additions & 0 deletions include/fluxes/numericalflux.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,18 @@ class NumericalFluxBase
*/
virtual double Flux( const Vector& Omega, double psiL, double psiR, const Vector& n ) const = 0;

/**
* @brief Flux computes flux on edge for fixed ordinate at a given edge using x and y axis
* @param quadPts ordinates for flux computation
* @param psiL left solution state
* @param psiR right solution state
* @param n scaled normal vector of given edge
* @param n_sys number of quadpts
* @param numerical flux value
*/
virtual void Flux( const VectorVector& quadPts, const Vector& psiL, const Vector& psiR, Vector& flux, const Vector& n, unsigned n_sys ) = 0;


/**
* @brief Flux computes flux on edge for fixed ordinate at a given edge (in XZ plane)
* @param Omega fixed ordinate for flux computation
Expand Down
12 changes: 12 additions & 0 deletions include/fluxes/upwindflux.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,18 @@ class UpwindFlux : public NumericalFluxBase
* @return numerical flux value
*/
double Flux( const Vector& Omega, double psiL, double psiR, const Vector& n ) const override;

/**
* @brief Flux computes flux on edge for fixed ordinate at a given edge using x and y axis
* @param quadPts ordinates for flux computation
* @param psiL left solution state
* @param psiR right solution state
* @param n scaled normal vector of given edge
* @param n_sys number of quadpts
* @param flux numerical flux value
*/
void Flux( const VectorVector& quadPts, const Vector& psiL, const Vector& psiR, Vector& flux, const Vector& n, unsigned n_sys ) override;

/**
* @brief FluxXZ computes flux on edge for fixed ordinate at a given edge using x and z axis
* @param Omega fixed ordinate for flux computation
Expand Down
14 changes: 14 additions & 0 deletions src/fluxes/upwindflux.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,20 @@ double UpwindFlux::Flux( const Vector& Omega, double psiL, double psiR, const Ve
}
}

void UpwindFlux::Flux( const VectorVector& quadPts, const Vector& psiL, const Vector& psiR, Vector& flux, const Vector& n, unsigned n_sys ) {
double inner; // Only use x and y axis in 2d case, minus because y axis is flipped for some reason

for (unsigned idx_q = 0; idx_q < n_sys; idx_q++){
inner = quadPts[idx_q][0] * n[0] + quadPts[idx_q][1] * n[1];
if( inner > 0 ) {
flux[idx_q]+= inner * psiL[idx_q];
}
else {
flux[idx_q] += inner * psiR[idx_q];
}
}
}

double UpwindFlux::FluxXZ( const Vector& Omega, double psiL, double psiR, const Vector& n ) const {
double inner = Omega[0] * n[0] + Omega[2] * n[1]; // Only use x and z axis in 2d case
if( inner > 0 ) {
Expand Down
94 changes: 44 additions & 50 deletions src/solvers/snsolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,68 +116,62 @@ void SNSolver::FluxUpdatePseudo2D() {
double solR;
// Dirichlet cells stay at IC, farfield assumption
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
// Loop over all ordinates
for( unsigned idx_quad = 0; idx_quad < _nq; ++idx_quad ) {
// Reset temporary variable
_solNew[idx_cell][idx_quad] = 0.0;
_solNew[idx_cell]*= 0.0; //blaze op

// Loop over all neighbor cells (edges) of cell j and compute numerical
// fluxes
for( unsigned idx_nbr = 0; idx_nbr < _neighbors[idx_cell].size(); ++idx_nbr ) {
// store flux contribution on psiNew_sigmaS to save memory
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_nbr] == _nCells ) {
_solNew[idx_cell][idx_quad] += _g->Flux( _quadPoints[idx_quad],
_sol[idx_cell][idx_quad],
_problem->GetGhostCellValue( idx_cell, _sol[idx_cell] )[idx_quad],
_normals[idx_cell][idx_nbr] );
_g->Flux( _quadPoints, _sol[idx_cell], _problem->GetGhostCellValue( idx_cell, _sol[idx_cell] ),
_solNew[idx_cell],
_normals[idx_cell][idx_nbr], _nq );
}
else {
unsigned int nbr_glob = _neighbors[idx_cell][idx_nbr]; // global idx of neighbor cell
// first order solver
_solNew[idx_cell][idx_quad] +=
_g->Flux( _quadPoints[idx_quad], _sol[idx_cell][idx_quad], _sol[nbr_glob][idx_quad], _normals[idx_cell][idx_nbr] );
_g->Flux( _quadPoints, _sol[idx_cell], _sol[ _neighbors[idx_cell][idx_nbr]], _solNew[idx_cell],_normals[idx_cell][idx_nbr], _nq );
}
}
}

}
}
else if( _reconsOrder == 2 ) {
// Loop over all spatial cells
#pragma omp parallel for
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
double solL;
double solR;
Vector solL;
Vector solR;
// Dirichlet cells stay at IC, farfield assumption
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
// Loop over all ordinates
for( unsigned idx_quad = 0; idx_quad < _nq; ++idx_quad ) {
// Reset temporary variable
_solNew[idx_cell][idx_quad] = 0.0;
_solNew[idx_cell] *= 0.0; // blaze operation
// Loop over all neighbor cells (edges) of cell j and compute numerical
// fluxes
for( unsigned idx_nbr = 0; idx_nbr < _neighbors[idx_cell].size(); ++idx_nbr ) {
// store flux contribution on psiNew_sigmaS to save memory
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_nbr] == _nCells ) {
_solNew[idx_cell][idx_quad] += _g->Flux( _quadPoints[idx_quad],
_sol[idx_cell][idx_quad],
_problem->GetGhostCellValue( idx_cell, _sol[idx_cell] )[idx_quad],
_normals[idx_cell][idx_nbr] );
}
else {
unsigned int nbr_glob = _neighbors[idx_cell][idx_nbr]; // global idx of neighbor cell
for( unsigned idx_nbr = 0; idx_nbr < _neighbors[idx_cell].size(); ++idx_nbr ) {
// store flux contribution on psiNew_sigmaS to save memory
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_nbr] == _nCells ) {
_g->Flux( _quadPoints, _sol[idx_cell], _problem->GetGhostCellValue( idx_cell, _sol[idx_cell] ),
_solNew[idx_cell],
_normals[idx_cell][idx_nbr], _nq );
}
else {
unsigned int nbr_glob = _neighbors[idx_cell][idx_nbr]; // global idx of neighbor cell

// left status of interface
solL = _sol[idx_cell][idx_quad] +
_limiter[idx_cell][idx_quad] *
( _solDx[idx_cell][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[idx_cell][0] ) +
_solDy[idx_cell][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[idx_cell][1] ) );
solR = _sol[nbr_glob][idx_quad] +
_limiter[nbr_glob][idx_quad] *
( _solDx[nbr_glob][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[nbr_glob][0] ) +
_solDy[nbr_glob][idx_quad] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[nbr_glob][1] ) );
// left status of interface
solL = _sol[idx_cell] +
_limiter[idx_cell] *
(_solDx[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[idx_cell][0] ) +
_solDy[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[idx_cell][1] ) );


// flux evaluation
_solNew[idx_cell][idx_quad] += _g->Flux( _quadPoints[idx_quad], solL, solR, _normals[idx_cell][idx_nbr] );
}
solR = _sol[nbr_glob] +
_limiter[nbr_glob] *
( _solDx[nbr_glob]* ( _interfaceMidPoints[idx_cell][idx_nbr][0] - _cellMidPoints[nbr_glob][0] ) +
_solDy[nbr_glob] * ( _interfaceMidPoints[idx_cell][idx_nbr][1] - _cellMidPoints[nbr_glob][1] ) );

// flux evaluation
_g->Flux( _quadPoints, solL, solR, _solNew[idx_cell], _normals[idx_cell][idx_nbr], _nq );
}
}
}
Expand All @@ -191,19 +185,19 @@ void SNSolver::FVMUpdate( unsigned idx_iter ) {
// Dirichlet cells stay at IC, farfield assumption
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
// loop over all ordinates
for( unsigned idx_quad = 0; idx_quad < _nq; ++idx_quad ) {
//for( unsigned idx_quad = 0; idx_quad < _nq; ++idx_quad ) {
// time update angular flux with numerical flux and total scattering cross
// section
_solNew[idx_cell][idx_quad] = _sol[idx_cell][idx_quad] - ( _dT / _areas[idx_cell] ) * _solNew[idx_cell][idx_quad] -
_dT * _sigmaT[0][idx_cell] * _sol[idx_cell][idx_quad];
}
_solNew[idx_cell] = _sol[idx_cell] - ( _dT / _areas[idx_cell] ) * _solNew[idx_cell] -
_dT * _sigmaT[0][idx_cell] * _sol[idx_cell];
//}
// compute scattering effects
_solNew[idx_cell] += _dT * _sigmaS[0][idx_cell] * _scatteringKernel * _sol[idx_cell]; // multiply scattering matrix with psi
// Source Term
if( _Q[0][idx_cell].size() == 1u ) // isotropic source
_solNew[idx_cell] += _dT * _Q[0][idx_cell][0];
else
_solNew[idx_cell] += _dT * _Q[0][idx_cell];
//if( _Q[0][idx_cell].size() == 1u ) // isotropic source
_solNew[idx_cell] += _dT * _Q[0][idx_cell][0];
//else
// _solNew[idx_cell] += _dT * _Q[0][idx_cell];
}
}
else {
Expand All @@ -212,12 +206,12 @@ void SNSolver::FVMUpdate( unsigned idx_iter ) {
// Dirichlet cells stay at IC, farfield assumption
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
// loop over all ordinates
for( unsigned idx_quad = 0; idx_quad < _nq; ++idx_quad ) {
//for( unsigned idx_quad = 0; idx_quad < _nq; ++idx_quad ) {
// time update angular flux with numerical flux and total scattering cross
// section
_solNew[idx_cell][idx_quad] = _sol[idx_cell][idx_quad] - ( _dT / _areas[idx_cell] ) * _solNew[idx_cell][idx_quad] -
_dT * _sigmaT[idx_iter][idx_cell] * _sol[idx_cell][idx_quad];
}
_solNew[idx_cell] = _sol[idx_cell] - ( _dT / _areas[idx_cell] ) * _solNew[idx_cell] -
_dT * _sigmaT[idx_iter][idx_cell] * _sol[idx_cell];
// }
// compute scattering effects
_solNew[idx_cell] += _dT * _sigmaS[idx_iter][idx_cell] * _scatteringKernel * _sol[idx_cell]; // multiply scattering matrix with psi

Expand Down