Commit fe81a279 authored by Steffen Schotthöfer's avatar Steffen Schotthöfer
Browse files

clean up solvers... replaced temp variable psiNew with _solNEw

parent 21109255
......@@ -52,8 +52,8 @@ class MNSolver : public Solver
void WriteOutputFields( unsigned idx_pseudoTime ) override;
// Solver
void FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) override;
void FluxUpdate( VectorVector& psiNew ) override;
void FVMUpdate( unsigned idx_energy ) override;
void FluxUpdate() override;
void IterPreprocessing() override;
void IterPostprocessing();
/*! @brief : Construct flux by computing the Moment of the sum of FVM discretization at the interface of cell
......
......@@ -48,8 +48,8 @@ class PNSolver : public Solver
void WriteOutputFields( unsigned idx_pseudoTime ) override;
// Solver
void FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) override;
void FluxUpdate( VectorVector& psiNew ) override;
void FVMUpdate( unsigned idx_energy ) override;
void FluxUpdate() override;
void IterPreprocessing() override;
void IterPostprocessing();
......
......@@ -25,8 +25,8 @@ class SNSolver : public Solver
void WriteOutputFields( unsigned idx_pseudoTime ) override;
// Solver
void FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) override;
void FluxUpdate( VectorVector& psiNew ) override;
void FVMUpdate( unsigned idx_energy ) override;
void FluxUpdate() override;
void IterPreprocessing() override;
void IterPostprocessing() override;
......
......@@ -89,9 +89,9 @@ class Solver
/*! @brief Performs postprocessing for the current solver iteration */
virtual void IterPostprocessing() = 0;
/*! @brief Constructs the flux update for the current iteration and stores it in psiNew*/
virtual void FluxUpdate( VectorVector& psiNew ) = 0;
virtual void FluxUpdate() = 0;
/*! @brief Computes the finite Volume update step for the current iteration */
virtual void FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) = 0;
virtual void FVMUpdate( unsigned idx_energy ) = 0;
// Helper
/**
......
......@@ -152,16 +152,16 @@ void MNSolver::ComputeRadFlux() {
}
}
void MNSolver::FluxUpdate( VectorVector& psiNew ) {
void MNSolver::FluxUpdate() {
// Loop over the grid cells
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
// Dirichlet Boundaries stay
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
psiNew[idx_cell] = ConstructFlux( idx_cell );
_solNew[idx_cell] = ConstructFlux( idx_cell );
}
}
void MNSolver::FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) {
void MNSolver::FVMUpdate( unsigned idx_energy ) {
// Loop over the grid cells
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
// Dirichlet Boundaries stay
......@@ -169,14 +169,14 @@ void MNSolver::FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) {
for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {
psiNew[idx_cell][idx_system] = _sol[idx_cell][idx_system] -
( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system] /* cell averaged flux */
- _dE * _sol[idx_cell][idx_system] *
( _sigmaT[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
_solNew[idx_cell][idx_system] = _sol[idx_cell][idx_system] -
( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_system] /* cell averaged flux */
- _dE * _sol[idx_cell][idx_system] *
( _sigmaT[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
}
psiNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
_solNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
}
}
......
......@@ -70,7 +70,7 @@ void PNSolver::ComputeRadFlux() {
}
}
void PNSolver::FluxUpdate( VectorVector& psiNew ) {
void PNSolver::FluxUpdate() {
// Loop over all spatial cells
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
......@@ -79,7 +79,7 @@ void PNSolver::FluxUpdate( VectorVector& psiNew ) {
// Reset temporary variable psiNew
for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
psiNew[idx_cell][idx_sys] = 0.0;
_solNew[idx_cell][idx_sys] = 0.0;
}
// Loop over all neighbor cells (edges) of cell j and compute numerical fluxes
......@@ -87,35 +87,36 @@ void PNSolver::FluxUpdate( VectorVector& psiNew ) {
// Compute flux contribution and store in psiNew to save memory
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_neighbor] == _nCells )
psiNew[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] );
else
psiNew[idx_cell] += _g->Flux( _AxPlus,
_AxMinus,
_AyPlus,
_AyMinus,
_AzPlus,
_AzMinus,
_sol[idx_cell],
_sol[_neighbors[idx_cell][idx_neighbor]],
_normals[idx_cell][idx_neighbor] );
_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] );
}
}
}
void PNSolver::FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) {
void PNSolver::FVMUpdate( unsigned idx_energy ) {
// Loop over all spatial cells
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
// Dirichlet cells stay at IC, farfield assumption
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
// Flux update
for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
psiNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys] - ( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_sys] /* cell averaged flux */
- _dE * _sol[idx_cell][idx_sys] *
( _sigmaT[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_sys] ); /* scattering influence */
_solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys] - ( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys] /* cell averaged flux */
- _dE * _sol[idx_cell][idx_sys] *
( _sigmaT[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_sys] ); /* scattering influence */
}
psiNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
// Source Term
_solNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
}
}
......
......@@ -42,7 +42,7 @@ void SNSolver::ComputeRadFlux() {
}
}
void SNSolver::FluxUpdate( VectorVector& psiNew ) {
void SNSolver::FluxUpdate() {
// Legacy Code form second order reconstruction:
......@@ -104,12 +104,12 @@ void SNSolver::FluxUpdate( VectorVector& psiNew ) {
// Loop over all ordinates
for( unsigned idx_quad = 0; idx_quad < _nq; ++idx_quad ) {
// Reset temporary variable
psiNew[idx_cell][idx_quad] = 0.0;
_solNew[idx_cell][idx_quad] = 0.0;
// Loop over all neighbor cells (edges) of cell j and compute numerical fluxes
for( unsigned idx_neighbor = 0; idx_neighbor < _neighbors[idx_cell].size(); ++idx_neighbor ) {
// store flux contribution on psiNew_sigmaS to save memory
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_neighbor] == _nCells )
psiNew[idx_cell][idx_quad] +=
_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 ) {
......@@ -146,10 +146,10 @@ void SNSolver::FluxUpdate( VectorVector& psiNew ) {
// default: first order solver
default:
*/
psiNew[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] );
_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] );
// }
}
}
......@@ -157,31 +157,31 @@ void SNSolver::FluxUpdate( VectorVector& psiNew ) {
}
}
void SNSolver::FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) {
void SNSolver::FVMUpdate( unsigned idx_energy ) {
for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
// 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 ) {
// time update angular flux with numerical flux and total scattering cross section
psiNew[idx_cell][idx_quad] = _sol[idx_cell][idx_quad] - ( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_quad] -
_dE * _sigmaT[idx_energy][idx_cell] * _sol[idx_cell][idx_quad];
_solNew[idx_cell][idx_quad] = _sol[idx_cell][idx_quad] - ( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_quad] -
_dE * _sigmaT[idx_energy][idx_cell] * _sol[idx_cell][idx_quad];
}
// compute scattering effects
psiNew[idx_cell] += _dE * _sigmaS[idx_energy][idx_cell] * _scatteringKernel * _sol[idx_cell]; // multiply scattering matrix with psi
_solNew[idx_cell] += _dE * _sigmaS[idx_energy][idx_cell] * _scatteringKernel * _sol[idx_cell]; // multiply scattering matrix with psi
// Source Term
if( _Q.size() == 1u ) { // constant source for all energies
if( _Q[0][idx_cell].size() == 1u ) // isotropic source
psiNew[idx_cell] += _dE * _Q[0][idx_cell][0];
_solNew[idx_cell] += _dE * _Q[0][idx_cell][0];
else
psiNew[idx_cell] += _dE * _Q[0][idx_cell];
_solNew[idx_cell] += _dE * _Q[0][idx_cell];
}
else {
if( _Q[0][idx_cell].size() == 1u ) // isotropic source
psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell][0];
_solNew[idx_cell] += _dE * _Q[idx_energy][idx_cell][0];
else
psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell];
_solNew[idx_cell] += _dE * _Q[idx_energy][idx_cell];
}
}
}
......
......@@ -164,10 +164,10 @@ void Solver::Solve() {
IterPreprocessing();
// --- Compute Fluxes ---
FluxUpdate( _solNew );
FluxUpdate();
// --- Finite Volume Update ---
FVMUpdate( _solNew, idx_energy );
FVMUpdate( idx_energy );
// --- Postprocessing ---
IterPostprocessing();
......
Supports Markdown
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