Skip to content
GitLab
Commit 432d0701 authored by Steffen Schotthöfer's avatar Steffen Schotthöfer
Browse files

some more pushes

parent 8728d2cf
Pipeline #142820 failed with stage
in 8 minutes and 21 seconds
Hide whitespace changes
Inline Side-by-side
code/include/solvers/csdpnsolver.h
View file @ 432d0701
......@@ -31,6 +31,8 @@ class CSDPNSolver : public PNSolver
private:
void IterPreprocessing( unsigned /*idx_iter*/ ) override;
void IterPostprocessing( unsigned /*idx_iter*/ ) override;
void ComputeRadFlux() override;
void FluxUpdate() override;
void FVMUpdate( unsigned idx_energy ) override;
......
code/include/solvers/pnsolver.h
View file @ 432d0701
......@@ -14,7 +14,7 @@ class PNSolver : public SolverBase
/*! @brief PNSolver destructor */
virtual ~PNSolver() {}
private:
protected:
unsigned _nSystem; /*!< @brief total number of equations in the system */
unsigned _polyDegreeBasis; /*!< @brief maximal degree of the spherical harmonics basis*/
......
code/src/solvers/csdpnsolver.cpp
View file @ 432d0701
#include "solvers/csdpnsolver.h"
#include "common/config.h"
#include "common/io.h"
#include "common/mesh.h"
#include "fluxes/numericalflux.h"
#include "kernels/scatteringkernelbase.h"
#include "problems/icru.h"
......@@ -10,7 +11,15 @@
#include "spdlog/spdlog.h"
#include <mpi.h>
CSDPNSolver::CSDPNSolver( Config* settings ) : PNSolver( settings ) {}
CSDPNSolver::CSDPNSolver( Config* settings ) : PNSolver( settings ) {
_dose = std::vector<double>( _settings->GetNCells(), 0.0 );
// only dummies for compilation
_energies = Vector( _nEnergies, 0.0 ); // equidistant
_angle = _energies;
_sigmaSE = { Matrix( _nEnergies, 0.0 ) };
_sigmaTE = Vector( _nEnergies, 0.0 );
}
void CSDPNSolver::IterPreprocessing( unsigned /*idx_iter*/ ) {
// Nothing to preprocess for PNSolver
......@@ -32,84 +41,85 @@ void CSDPNSolver::ComputeRadFlux() {
}
void CSDPNSolver::FluxUpdate() {
if( _reconsOrder > 1 ) {
_mesh->ReconstructSlopesU( _nSystem, _solDx, _solDy, _sol ); // unstructured reconstruction
//_mesh->ComputeSlopes( _nTotalEntries, _solDx, _solDy, _sol ); // unstructured reconstruction
}
// 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++ ) {
// 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 < _nSystem; idx_sys++ ) {
_solNew[idx_cell][idx_sys] = 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++ ) {
// Compute flux contribution and store in psiNew to save memory
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_neighbor] == _nCells )
_solNew[idx_cell] += _g->Flux(
_AxPlus, _AxMinus, _AyPlus, _AyMinus, _AzPlus, _AzMinus, _sol[idx_cell], _sol[idx_cell], _normals[idx_cell][idx_neighbor] );
else {
switch( _reconsOrder ) {
// first order solver
case 1:
_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] );
break;
// second order solver
case 2:
// 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] );
// positivity checker (if not satisfied, deduce to first order)
// I manually turned it off here since Pn produces negative solutions essentially
// 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, _AyPlus, _AyMinus, _AzPlus, _AzMinus, solL, solR, _normals[idx_cell][idx_neighbor] );
break;
// default: first order solver
default:
_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] );
}
}
}
}
// if( _reconsOrder > 1 ) {
// _mesh->ReconstructSlopesU( _nSystem, _solDx, _solDy, _sol ); // unstructured reconstruction
// //_mesh->ComputeSlopes( _nTotalEntries, _solDx, _solDy, _sol ); // unstructured reconstruction
//}
//// 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++ ) {
//
// // 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 < _nSystem; idx_sys++ ) {
// _solNew[idx_cell][idx_sys] = 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++ ) {
//
// // Compute flux contribution and store in psiNew to save memory
// if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_neighbor] == _nCells )
// _solNew[idx_cell] += _g->Flux(
// _AxPlus, _AxMinus, _AyPlus, _AyMinus, _AzPlus, _AzMinus, _sol[idx_cell], _sol[idx_cell], _normals[idx_cell][idx_neighbor] );
// else {
// switch( _reconsOrder ) {
// // first order solver
// case 1:
// _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] );
// break;
// // second order solver
// case 2:
// // 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] );
// // positivity checker (if not satisfied, deduce to first order)
// // I manually turned it off here since Pn produces negative solutions essentially
// // 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, _AyPlus, _AyMinus, _AzPlus, _AzMinus, solL, solR, _normals[idx_cell][idx_neighbor] );
// break;
// // default: first order solver
// default:
// _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 CSDPNSolver::FVMUpdate( unsigned idx_energy ) {}
......
code/src/solvers/pnsolver.cpp
View file @ 432d0701
......@@ -72,7 +72,7 @@ void PNSolver::ComputeRadFlux() {
}
}
void PNSolver::FluxUpdate() override {
void PNSolver::FluxUpdate() {
if( _reconsOrder > 1 ) {
_mesh->ReconstructSlopesU( _nSystem, _solDx, _solDy, _sol ); // unstructured reconstruction
//_mesh->ComputeSlopes( _nTotalEntries, _solDx, _solDy, _sol ); // unstructured reconstruction
......@@ -153,7 +153,7 @@ void PNSolver::FluxUpdate() override {
}
}
void PNSolver::FVMUpdate( unsigned idx_energy ) override {
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
......@@ -370,7 +370,7 @@ double PNSolver::LegendrePoly( double x, int l ) { // Legacy. TO BE DELETED
}
}
void PNSolver::PrepareVolumeOutput() override {
void PNSolver::PrepareVolumeOutput() {
unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
_outputFieldNames.resize( nGroups );
......@@ -410,7 +410,7 @@ void PNSolver::PrepareVolumeOutput() override {
}
}
void PNSolver::WriteVolumeOutput( unsigned idx_pseudoTime ) override {
void PNSolver::WriteVolumeOutput( unsigned idx_pseudoTime ) {
unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||
......
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

Imprint