Skip to content

GitLab

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

encapsulated all three base solvers 


Former-commit-id: 64c5cba2
parent ea89c2ee
Hide whitespace changes
Inline Side-by-side
code/include/solvers/mnsolver.h
View file @ 57ffed18
......@@ -79,5 +79,10 @@ class MNSolver : public Solver
/*! @brief Corrects the solution _sol[idx_cell] to be realizable w.r.t. the reconstructed entropy (eta'(alpha*m))
@param idx_cell = cell where the correction happens*/
void ComputeRealizableSolution( unsigned idx_cell );
// Solver
void FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) override;
void FluxUpdate( VectorVector& psiNew ) override;
void IterPreprocessing() override;
};
#endif // MNSOLVER_H
code/include/solvers/pnsolver.h
View file @ 57ffed18
......@@ -52,12 +52,9 @@ class PNSolver : public Solver
double WriteOutputFields( unsigned idx_pseudoTime ) override;
// Solver
/*! @brief Computes the finite volume update for all cells and stores it in psiNew */
void FVMUpdate( VectorVector& psiNew, unsigned idx_energy );
/*! @brief Computes the fluxes for all cells and stores it in psiNew, sets Neumann Boundaries */
void FluxUpdate( VectorVector& psiNew, unsigned idx_energy );
/*! @brief Sets Dirichlet Boundaries and resets temp Variables */
void PrepIteration( VectorVector& psiNew );
void FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) override;
void FluxUpdate( VectorVector& psiNew ) override;
void IterPreprocessing() override;
// Initialization
/*! @brief: parameter functions for setting up system matrix
......
code/include/solvers/snsolver.h
View file @ 57ffed18
......@@ -28,9 +28,11 @@ class SNSolver : public Solver
void PrepareOutputFields() override;
double WriteOutputFields( unsigned idx_pseudoTime ) override;
void PrepIteration( VectorVector& psiNew );
void FluxUpdate( VectorVector& psiNew, unsigned idx_energy );
void FVMUpdate( VectorVector& psiNew, unsigned idx_energy );
// Solver
void FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) override;
void FluxUpdate( VectorVector& psiNew ) override;
void IterPreprocessing() override;
// --- Member variables ---
};
......
code/include/solvers/solverbase.h
View file @ 57ffed18
......@@ -74,6 +74,13 @@ class Solver
///*! @brief Computes the finite volume update for all cells and stores it in psiNew */
// virtual void FVMUpdate( VectorVector& psiNew ) = 0;
/*! @brief Performs preprocessing for the current solver iteration */
virtual void IterPreprocessing() = 0;
/*! @brief Constructs the flux update for the current iteration and stores it in psiNew*/
virtual void FluxUpdate( VectorVector& psiNew ) = 0;
/*! @brief Computes the finite Volume update step for the current iteration */
virtual void FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) = 0;
/**
* @brief ComputeTimeStep calculates the maximal stable time step
* @param cfl is cfl number
......
code/src/solvers/mnsolver.cpp
View file @ 57ffed18
......@@ -132,63 +132,74 @@ void MNSolver::Solve() {
auto log = spdlog::get( "event" );
// angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
VectorVector psiNew = _sol;
double mass = 0;
Save( -1 );
if( rank == 0 ) log->info( "{:10} {:10}", "t", "mass" );
// if( rank == 0 ) log->info( "{:03.8f} {:01.5e} {:01.5e}", -1.0, dFlux, mass1 ); Should be deleted
// Loop over energies (pseudo-time of continuous slowing down approach)
for( unsigned idx_energy = 0; idx_energy < _nEnergies; idx_energy++ ) {
// ------- Reconstruction Step ----------------
// ----- Iteration Preprocessing -----
IterPreprocessing();
_optimizer->SolveMultiCell( _alpha, _sol, _moments );
// ------- Flux Computation Step --------------
FluxUpdate( psiNew );
// Loop over the grid cells
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
// ------ Finite Volume Update Step ------
FVMUpdate( psiNew, idx_energy );
// ------- Relizablity Reconstruction Step ----
// Update Solution
_sol = psiNew;
ComputeRealizableSolution( idx_cell );
// --- VTK and CSV Output ---
mass = WriteOutputFields( idx_energy );
Save( idx_energy );
// WriteNNTrainingData( idx_energy );
// ------- Flux Computation Step --------------
// --- Screen Output ---
if( rank == 0 ) log->info( "{:03.8f} {:01.5e}", _energies[idx_energy], mass );
}
}
// Dirichlet Boundaries are finished now
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
void MNSolver::IterPreprocessing() {
psiNew[idx_cell] = ConstructFlux( idx_cell );
// ------- Reconstruction Step ----------------
// ------ Finite Volume Update Step ------
_optimizer->SolveMultiCell( _alpha, _sol, _moments );
// NEED TO VECTORIZE
for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {
// ------- Relizablity Reconstruction Step ----
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
ComputeRealizableSolution( idx_cell );
}
}
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 */
}
void MNSolver::FluxUpdate( VectorVector& psiNew ) {
// 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 );
}
}
psiNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
}
void MNSolver::FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) {
// 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;
// Update Solution
_sol = psiNew;
for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {
// --- VTK and CSV Output ---
mass = WriteOutputFields( idx_energy );
Save( idx_energy );
// WriteNNTrainingData( idx_energy );
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 */
}
// --- Screen Output ---
if( rank == 0 ) log->info( "{:03.8f} {:01.5e}", _energies[idx_energy], mass );
psiNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
}
}
......
code/src/solvers/pnsolver.cpp
View file @ 57ffed18
......@@ -62,24 +62,16 @@ void PNSolver::Solve() {
double mass = 0;
Save( -1 ); // Save initial condition
unsigned idx_system = 0;
if( rank == 0 ) log->info( "{:10} {:10}", "t", "mass" );
// Remove
mass = WriteOutputFields( 0 );
if( rank == 0 ) log->info( " {:01.5e} {:01.5e}", 0.0, mass );
// Loop over energies (pseudo-time of continuous slowing down approach)
for( unsigned idx_energy = 0; idx_energy < _nEnergies; idx_energy++ ) {
// --- Prepare Boundaries and temp variables
PrepIteration( psiNew );
IterPreprocessing();
// --- Compute Fluxes ---
FluxUpdate( psiNew, idx_energy );
FluxUpdate( psiNew );
// --- Finite Volume Update ---
FVMUpdate( psiNew, idx_energy );
......@@ -96,21 +88,22 @@ void PNSolver::Solve() {
}
}
void PNSolver::PrepIteration( VectorVector& psiNew ) {
void PNSolver::IterPreprocessing() {
// Nothing to preprocess for PNSolver
}
void PNSolver::FluxUpdate( VectorVector& psiNew ) {
// Loop over all spatial cells
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue; // Dirichlet cells stay at IC, farfield assumption
// Dirichlet cells stay at IC, farfield assumption
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
// Reset temporary variable psiNew
for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
psiNew[idx_cell][idx_sys] = 0.0;
}
}
}
void PNSolver::FluxUpdate( VectorVector& psiNew, unsigned idx_energy ) {
// Loop over all spatial cells
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
// 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++ ) {
......@@ -133,10 +126,12 @@ void PNSolver::FluxUpdate( VectorVector& psiNew, unsigned idx_energy ) {
}
void PNSolver::FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) {
// time update angular flux with numerical flux and total scattering cross section
// 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 */
......
code/src/solvers/snsolver.cpp
View file @ 57ffed18
......@@ -25,74 +25,29 @@ SNSolver::SNSolver( Config* settings ) : Solver( settings ) {
}
void SNSolver::Solve() {
auto log = spdlog::get( "event" );
// angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
VectorVector psiNew = _sol;
// reconstruction order
unsigned reconsOrder = _settings->GetReconsOrder();
// 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 ) );
VectorVector psiDy( _nCells, Vector( _nq, 0.0 ) );
// geometric variables for derivatives computation
auto nodes = _mesh->GetNodes();
auto cells = _mesh->GetCells();
auto cellMidPoints = _mesh->GetCellMidPoints();
int rank;
MPI_Comm_rank( MPI_COMM_WORLD, &rank );
// center location of cell interfaces
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] );
}
}
auto log = spdlog::get( "event" );
// distance between cell center to interface center
// VectorVector cellDx( _nCells, Vector( _mesh->GetNumNodesPerCell(), 1e-10 ) );
// VectorVector cellDy( _nCells, Vector( _mesh->GetNumNodesPerCell(), 1e-10 ) );
//
// for( unsigned i = 0; i < _nCells; ++i ) {
// for( unsigned j = 0; j < _mesh->GetNumNodesPerCell(); ++j ) {
// cellDx[i][j] = interfaceMidPoints[i][j][0] - cellMidPoints[i][0];
// cellDy[i][j] = interfaceMidPoints[i][j][1] - cellMidPoints[j][1];
// }
// }
VectorVector psiNew = _sol;
double mass = 0.0;
double dFlux = 1e10;
Vector fluxNew( _nCells, 0.0 );
Vector fluxOld( _nCells, 0.0 );
int rank;
MPI_Comm_rank( MPI_COMM_WORLD, &rank );
if( rank == 0 ) log->info( "{:10} {:10}", "t", "dFlux", "mass" );
// if( rank == 0 ) log->info( "{:01.5e} {:01.5e} {:01.5e}", 0.0, dFlux, mass1 );
// loop over energies (pseudo-time)
for( unsigned idx_energy = 0; idx_energy < _nEnergies; ++idx_energy ) {
// reconstruct slopes for higher order solver
// if( reconsOrder > 1 ) {
// _mesh->ReconstructSlopesU( _nq, psiDx, psiDy, _sol ); // unstructured reconstruction
// //_mesh->ReconstructSlopesS( _nq, psiDx, psiDy, _psi ); // structured reconstruction (not stable currently)
//}
// --- Prepare Boundaries and temp variables
PrepIteration( psiNew );
IterPreprocessing();
// --- Compute Fluxes ---
FluxUpdate( psiNew, idx_energy );
FluxUpdate( psiNew );
// --- Compute Finite Volume Step ---
FVMUpdate( psiNew, idx_energy );
......@@ -115,29 +70,79 @@ void SNSolver::Solve() {
}
}
void SNSolver::PrepIteration( VectorVector& psiNew ) {
// Loop over all spatial cells
for( unsigned idx_cells = 0; idx_cells < _nCells; ++idx_cells ) {
if( _boundaryCells[idx_cells] == BOUNDARY_TYPE::DIRICHLET ) continue; // Dirichlet cells stay at IC, farfield assumption
void SNSolver::IterPreprocessing() {
// Nothing to do for SNSolver
}
// loop over all ordinates
for( unsigned idx_quad = 0; idx_quad < _nq; ++idx_quad ) {
psiNew[idx_cells][idx_quad] = 0.0;
void SNSolver::FluxUpdate( VectorVector& psiNew ) {
// Legacy Code form second order reconstruction:
/*
// reconstruction order
unsigned reconsOrder = _settings->GetReconsOrder();
// 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 ) );
VectorVector psiDy( _nCells, Vector( _nq, 0.0 ) );
// geometric variables for derivatives computation
auto nodes = _mesh->GetNodes();
auto cells = _mesh->GetCells();
auto cellMidPoints = _mesh->GetCellMidPoints();
// center location of cell interfaces
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] );
}
}
}
void SNSolver::FluxUpdate( VectorVector& psiNew, unsigned idx_energy ) {
// loop over all spatial cells
for( unsigned idx_cells = 0; idx_cells < _nCells; ++idx_cells ) {
// loop over all ordinates
distance between cell center to interface center
VectorVector cellDx( _nCells, Vector( _mesh->GetNumNodesPerCell(), 1e-10 ) );
VectorVector cellDy( _nCells, Vector( _mesh->GetNumNodesPerCell(), 1e-10 ) );
for( unsigned i = 0; i < _nCells; ++i ) {
for( unsigned j = 0; j < _mesh->GetNumNodesPerCell(); ++j ) {
cellDx[i][j] = interfaceMidPoints[i][j][0] - cellMidPoints[i][0];
cellDy[i][j] = interfaceMidPoints[i][j][1] - cellMidPoints[j][1];
}
}
*/
/*
// reconstruct slopes for higher order solver
if( reconsOrder > 1 ) {
_mesh->ReconstructSlopesU( _nq, psiDx, psiDy, _sol ); // unstructured reconstruction
//_mesh->ReconstructSlopesS( _nq, psiDx, psiDy, _psi ); // structured reconstruction (not stable currently)
}
*/
// 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;
// Loop over all ordinates
for( unsigned idx_quad = 0; idx_quad < _nq; ++idx_quad ) {
// loop over all neighbor cells (edges) of cell j and compute numerical fluxes
for( unsigned idx_neighbor = 0; idx_neighbor < _neighbors[idx_cells].size(); ++idx_neighbor ) {
// Reset temporary variable
psiNew[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_cells] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cells][idx_neighbor] == _nCells )
psiNew[idx_cells][idx_quad] +=
_g->Flux( _quadPoints[idx_quad], _sol[idx_cells][idx_quad], _sol[idx_cells][idx_quad], _normals[idx_cells][idx_neighbor] );
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_neighbor] == _nCells )
psiNew[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
......@@ -173,10 +178,10 @@ void SNSolver::FluxUpdate( VectorVector& psiNew, unsigned idx_energy ) {
// default: first order solver
default:
*/
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] );
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] );
// }
}
}
......@@ -185,28 +190,30 @@ void SNSolver::FluxUpdate( VectorVector& psiNew, unsigned idx_energy ) {
}
void SNSolver::FVMUpdate( VectorVector& psiNew, unsigned idx_energy ) {
for( unsigned idx_cells = 0; idx_cells < _nCells; ++idx_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;
// 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_cells][idx_quad] = _sol[idx_cells][idx_quad] - ( _dE / _areas[idx_cells] ) * psiNew[idx_cells][idx_quad] -
_dE * _sigmaT[idx_energy][idx_cells] * _sol[idx_cells][idx_quad];
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];
}
// compute scattering effects
psiNew[idx_cells] += _dE * _sigmaS[idx_energy][idx_cells] * _scatteringKernel * _sol[idx_cells]; // multiply scattering matrix with psi
psiNew[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_cells].size() == 1u ) // isotropic source
psiNew[idx_cells] += _dE * _Q[0][idx_cells][0];
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];
else
psiNew[idx_cells] += _dE * _Q[0][idx_cells];
psiNew[idx_cell] += _dE * _Q[0][idx_cell];
}
else {
if( _Q[0][idx_cells].size() == 1u ) // isotropic source
psiNew[idx_cells] += _dE * _Q[idx_energy][idx_cells][0];
if( _Q[0][idx_cell].size() == 1u ) // isotropic source
psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell][0];
else
psiNew[idx_cells] += _dE * _Q[idx_energy][idx_cells];
psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell];
}
}
}
......
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