Skip to content

GitLab

Commit 5df1f107 authored by steffen.schotthoefer's avatar steffen.schotthoefer
Browse files

mn solver started

parent 76518bcf
Hide whitespace changes
Inline Side-by-side
code/include/solvers/mnsolver.h
View file @ 5df1f107
......@@ -2,15 +2,53 @@
#define MNSOLVER_H
#include <pnsolver.h>
#include <sphericalharmonics.h>
class MNSolver : PNSolver
class MNSolver : Solver
{
/*! @brief: compute the legendre function of degree l and order k at point x using the lth degree
* legendre polynomial
* @param: [in] double x : spatial point, -1 <= x <= 1
* [in] int l : degree of associated legendre polynomial, 0 <= l
* [in] int k : order of associated legendre polynomial, -l <= k <= l
* [ou] double : value of legendre polynomial of degree l and order k at point x
public:
/**
* @brief MNSolver constructor
* @param settings stores all needed information
*/
double AssociatedLegendrePoly( double x, int l, int k );
MNSolver( Config* settings );
/**
* @brief Solve functions runs main time loop
*/
void Solve() override;
private:
/*! @brief: Total number of equations in the system */
unsigned _nTotalEntries;
/*! @brief: Absorbtion coefficient for all energies */
VectorVector _sigmaA;
/*! @brief: Class to compute and store current spherical harmonics basis */
SphericalHarmonics _basis;
/*! @brief: Diagonal of the scattering matrix (its a diagonal matrix by construction) */
Vector _scatterMatDiag;
/*! @brief: Diagonal of the system matrix in x direction (its a diagonal matrix by construction) */
Vector _Ax;
/*! @brief: Diagonal of the system matrix in y direction (its a diagonal matrix by construction) */
Vector _Ay;
/*! @brief: Diagonal of the system matrix in z direction (its a diagonal matrix by construction) */
Vector _Az;
/*! @brief : computes the global index of the moment corresponding to basis function (l,k)
* @param : degree l, it must hold: 0 <= l <=_nq
* @param : order k, it must hold: -l <=k <= l
* @returns : global index
*/
int GlobalIndex( int l, int k ) const;
/*! @brief : function for computing and setting up the System Matrices.
The System Matrices are diagonal, so there is no need to diagonalize*/
void ComputeSystemMatrices();
};
#endif // MNSOLVER_H
code/include/solvers/pnsolver.h
View file @ 5df1f107
......@@ -16,13 +16,13 @@ class PNSolver : public Solver
/**
* @brief Solve functions runs main time loop
*/
void Solve() override;
virtual void Solve() override;
/**
* @brief Output solution to VTK file
*/
void Save() const override;
void Save( int currEnergy ) const;
void Save( int currEnergy ) const override;
protected:
// moment orders for P_N
......@@ -87,8 +87,6 @@ class PNSolver : public Solver
void ComputeScatterMatrix();
// Computes Legedre polinomial of oder l at point x
double LegendrePoly( double x, int l );
// Adapt the TimeStep according to the CFL number of the velocities in the AdvectionMatrices
void AdaptTimeStep();
// Sets Entries of FluxMatrices to zero, if they are below double accuracy, to prevent floating point
// inaccuracies lateron
void CleanFluxMatrices();
......
code/include/solvers/snsolver.h
View file @ 5df1f107
......@@ -17,12 +17,12 @@ class SNSolver : public Solver
/**
* @brief Solve functions runs main time loop
*/
virtual void Solve();
void Solve() override;
/**
* @brief Output solution to VTK file
*/
virtual void Save() const;
virtual void Save( int currEnergy ) const;
void Save() const override;
void Save( int currEnergy ) const override;
};
#endif // SNSOLVER_H
code/include/solvers/sphericalharmonics.h
View file @ 5df1f107
......@@ -16,7 +16,7 @@ class SphericalHarmonics
{
public:
/*! @brief : Sets up class for spherical harmonics basis based on legendre
* polynoms and associated legendre polynoms up to degree l.
* polynoms and associated legendre polynoms up to degree L.
* The basis then consists of N = L² +2L basis functions.
* @param : L_degree - maximum degree of spherical harmonics basis, 0 <= L <= 1000 (upper bound
* due to numerical stability)
......@@ -31,6 +31,7 @@ class SphericalHarmonics
std::vector<double> ComputeSphericalBasis( double my, double phi );
private:
/*! @brief: maximal degree of the spherical harmonics basis (this is "L" in the comments)*/
unsigned _LMaxDegree;
/*! @brief: coefficients for the computations of the basis
......
code/src/fluxes/upwindflux.cpp
View file @ 5df1f107
......@@ -13,7 +13,7 @@ double UpwindFlux::Flux( const Vector& Omega, double psiL, double psiR, const Ve
}
/**
* @brief Flux : Computes <VanLeer> upwinding scheme for given flux jacobians of the PN Solver at a given edge and stores it in
* @brief Flux : Computes <Linear> upwinding scheme for given flux jacobians of the PN Solver at a given edge and stores it in
* resultFlux
* @param AxPlus : Positive part of the flux jacobian in x direction
* @param AxMinus : Negative part of the flux jacobian in x direction
......
code/src/solvers/mnsolver.cpp
View file @ 5df1f107
#include <mnsolver.h>
MNSolver::MNSolver( Config* settings )
: Solver( settings ), _basis( _nq ) /* We need max degree +1, since we need the N+1st moment to reconstruct phi */
{
// Is this good (fast) code using a constructor list?
_nTotalEntries = GlobalIndex( _nq, int( _nq ) ) + 1;
// transform sigmaT and sigmaS in sigmaA.
_sigmaA = VectorVector( _nEnergies, Vector( _nCells, 0 ) ); // Get rid of this extra vektor!
for( unsigned n = 0; n < _nEnergies; n++ ) {
for( unsigned j = 0; j < _nCells; j++ ) {
_sigmaA[n][j] = 0; //_sigmaT[n][j] - _sigmaS[n][j];
_sigmaS[n][j] = 1;
}
}
// Initialize Scatter Matrix
_scatterMatDiag = Vector( _nTotalEntries, 1.0 );
_scatterMatDiag[0] = 0.0; // First entry is zero by construction.
// Initialize System Matrices
_Ax = Vector( _nTotalEntries, 0.0 );
_Ay = Vector( _nTotalEntries, 0.0 );
_Az = Vector( _nTotalEntries, 0.0 );
// Fill System Matrices
ComputeSystemMatrices();
}
int MNSolver::GlobalIndex( int l, int k ) const {
int numIndicesPrevLevel = l * l; // number of previous indices untill level l-1
int prevIndicesThisLevel = k + l; // number of previous indices in current level
return numIndicesPrevLevel + prevIndicesThisLevel;
}
void MNSolver::ComputeSystemMatrices() {}
void MNSolver::Solve() {
int rank;
MPI_Comm_rank( MPI_COMM_WORLD, &rank );
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 = _psi;
double dFlux = 1e10;
Vector fluxNew( _nCells, 0.0 );
Vector fluxOld( _nCells, 0.0 );
unsigned idx_system = 0;
double mass1 = 0;
for( unsigned i = 0; i < _nCells; ++i ) {
_solverOutput[i] = _psi[i][0];
mass1 += _psi[i][0];
}
if( rank == 0 ) log->info( "{:10} {:10}", "t", "dFlux" );
if( rank == 0 ) log->info( "{:03.8f} {:01.5e} {:01.5e}", -1.0, dFlux, mass1 );
// Loop over energies (pseudo-time of continuous slowing down approach)
for( unsigned idx_energy = 0; idx_energy < _nEnergies; ++idx_energy ) {
// ------- Reconstruction Step -------
// Todo
// ------- Advection Step ------------
// 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
// Loop over all equations of the system.
for( int idx_lDegree = 0; idx_lDegree <= int( _nq ); idx_lDegree++ ) {
for( int idx_kOrder = -idx_lDegree; idx_kOrder <= idx_lDegree; idx_kOrder++ ) {
idx_system = unsigned( GlobalIndex( idx_lDegree, idx_kOrder ) );
psiNew[idx_cell][idx_system] = 0.0;
// Loop over all neighbor cells (edges) of cell j and compute numerical fluxes
for( unsigned l = 0; l < _neighbors[idx_cell].size(); ++l ) {
// store flux contribution on psiNew_sigmaS to save memory
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][l] == _nCells )
psiNew[idx_cell][idx_system] +=
_g->Flux( _quadPoints[idx_system], _psi[idx_cell][idx_system], _psi[idx_cell][idx_system], _normals[idx_cell][l] );
else
psiNew[idx_cell][idx_system] += _g->Flux( _quadPoints[idx_system],
_psi[idx_cell][idx_system],
_psi[_neighbors[idx_cell][l]][idx_system],
_normals[idx_cell][l] );
}
// time update angular flux with numerical flux and total scattering cross section
psiNew[idx_cell][idx_system] = _psi[idx_cell][idx_system] -
( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system] /* cell averaged flux */
- _dE * _psi[idx_cell][idx_system] *
( _sigmaA[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
}
}
// TODO: figure out a more elegant way
// add external source contribution
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_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];
else
psiNew[idx_cell] += _dE * _Q[idx_energy][idx_cell];
}
}
_psi = psiNew;
for( unsigned i = 0; i < _nCells; ++i ) {
fluxNew[i] = dot( _psi[i], _weights );
_solverOutput[i] = fluxNew[i];
}
// Save( n );
dFlux = blaze::l2Norm( fluxNew - fluxOld );
fluxOld = fluxNew;
if( rank == 0 ) log->info( "{:03.8f} {:01.5e}", _energies[n], dFlux );
}
}
code/src/solvers/pnsolver.cpp
View file @ 5df1f107
#include "solvers/pnsolver.h"
#include "../ext/blaze/blaze/math/blas/cblas/gemm.h"
#include "toolboxes/errormessages.h"
#include "toolboxes/textprocessingtoolbox.h"
#include <mpi.h>
......@@ -107,9 +106,9 @@ void PNSolver::Solve() {
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue; // Dirichlet cells stay at IC, farfield assumption
// Reset temporary variable psiNew
for( int idx_lOrder = 0; idx_lOrder <= int( _nq ); idx_lOrder++ ) {
for( int idx_kOrder = -idx_lOrder; idx_kOrder <= idx_lOrder; idx_kOrder++ ) {
idx_system = unsigned( GlobalIndex( idx_lOrder, idx_kOrder ) );
for( int idx_lDegree = 0; idx_lDegree <= int( _nq ); idx_lDegree++ ) {
for( int idx_kOrder = -idx_lDegree; idx_kOrder <= idx_lDegree; idx_kOrder++ ) {
idx_system = unsigned( GlobalIndex( idx_lDegree, idx_kOrder ) );
psiNew[idx_cell][idx_system] = 0.0;
}
}
......@@ -131,8 +130,6 @@ void PNSolver::Solve() {
_psi[idx_cell],
_psi[_neighbors[idx_cell][idx_neighbor]],
_normals[idx_cell][idx_neighbor] );
// std::cout << "psiNew is " << psiNew[idx_cell] << "\n ";
}
// time update angular flux with numerical flux and total scattering cross section
......@@ -140,15 +137,11 @@ void PNSolver::Solve() {
for( int idx_kOrder = -idx_lOrder; idx_kOrder <= idx_lOrder; idx_kOrder++ ) {
idx_system = unsigned( GlobalIndex( idx_lOrder, idx_kOrder ) );
// std::cout << "psi_old " << _psi[idx_cell][idx_system] << " de " << _dE << " _areas[idx_cell] " << _areas[idx_cell] << " Flux "
// << psiNew[idx_cell][idx_system] << "\n";
psiNew[idx_cell][idx_system] = _psi[idx_cell][idx_system] -
( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_system] /* cell averaged flux */
- _dE * _psi[idx_cell][idx_system] *
( _sigmaA[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
// std::cout << "psi_new " << psiNew[idx_cell][idx_system] << "\n";
}
}
}
......@@ -169,8 +162,6 @@ void PNSolver::Solve() {
}
}
void PNSolver::AdaptTimeStep() { _dE = _dE / _combinedSpectralRadius; }
void PNSolver::ComputeSystemMatrices() {
int idx_col = 0;
unsigned idx_row = 0;
......@@ -345,38 +336,6 @@ void PNSolver::ComputeFluxComponents() {
void PNSolver::ComputeScatterMatrix() {
// right now independent of spatial coordinates
// G[i][i] = 1- g_l
// g_l = 2*pi*\int_-1 ^1 k(x,my)P_l(my)dmy
// use midpoint rule.
// Here is room for optimization!
// Only isotropic scattering: k = 1/(2pi) in 2d !!!
// --- Ansatz for general scattering for later! ---
// unsigned nSteps = 500000;
// double integralSum = 0;
// unsigned idx_global = 0;
//
// for( int idx_lOrder = 0; idx_lOrder <= int( _nq ); idx_lOrder++ ) {
// for( int idx_kOrder = -idx_lOrder; idx_kOrder <= idx_lOrder; idx_kOrder++ ) {
// idx_global = unsigned( GlobalIndex( idx_lOrder, idx_kOrder ) );
// _scatterMatDiag[idx_global] = 0;
//
// // compute value of diagonal
// integralSum = 0.5 * ( Legendre( -1.0, idx_lOrder ) + Legendre( 1.0, idx_lOrder ) ); // Boundary terms
//
// for( unsigned i = 1; i < nSteps - 1; i++ ) {
// integralSum += Legendre( -1.0 + double( i ) * 2 / double( nSteps ), idx_lOrder );
// }
// integralSum *= 2 / double( nSteps );
//
// // isotropic scattering and prefactor of eigenvalue
// integralSum *= 0.5; // 2pi/4pi
// _scatterMatDiag[idx_global] = 1 - integralSum;
// }
//}
// --- Isotropic ---
_scatterMatDiag[0] = 0.0;
for( unsigned idx_diag = 1; idx_diag < _nTotalEntries; idx_diag++ ) {
......
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