Commit 6b123632 authored by Jonas Kusch's avatar Jonas Kusch
Browse files

fixed segmentation faults in solver

parent dc63188c
Pipeline #83910 passed with stages
......@@ -62,6 +62,11 @@ class Mesh
const std::vector<std::vector<unsigned>>& GetNeighbours() const;
const std::vector<std::vector<Vector>>& GetNormals() const;
/**
* @brief GetBoundaryCellArray returns array with (dirichlet?) boundary cells
*/
const std::vector<bool>& GetBoundaryCellArray() const;
/**
* @brief ComputeSlopes calculates the slope in every cell into x and y direction
* @param nq is number of quadrature points
......
......@@ -25,9 +25,12 @@ class Solver
std::vector<Matrix> _sigmaSH20; // scattering cross section, dim(_sigmaSH20) = (_nTimeSteps,_nq,_nq)
VectorVector _quadPoints; // quadrature points, dim(_quadPoints) = (_nTimeSteps,spatialDim)
Vector _weights; // quadrature weights, dim(_weights) = (_NCells)
std::vector<bool> _boundaryCells; // boundary type for all cells, dim(_boundary) = (_NCells)
// we will have to add a further dimension for quadPoints and weights once we start with multilevel SN
NumericalFlux* _g; // numerical flux function
Mesh* _mesh; // mesh object for writing out information
Settings* _settings;
/**
* @brief LoadPatientDensity loads density of patient from MRT/CT scan and saves it in _density
......
......@@ -307,3 +307,4 @@ const std::vector<double>& Mesh::GetCellAreas() const { return _cellAreas; }
const std::vector<unsigned>& Mesh::GetPartitionIDs() const { return _colors; }
const std::vector<std::vector<unsigned>>& Mesh::GetNeighbours() const { return _cellNeighbors; }
const std::vector<std::vector<Vector>>& Mesh::GetNormals() const { return _cellNormals; }
const std::vector<bool>& Mesh::GetBoundaryCellArray() const { return _isBoundaryCell; }
......@@ -14,7 +14,7 @@ void SNSolver::Solve() {
for( unsigned n = 0; n < _nTimeSteps; ++n ) {
// loop over all spatial cells
for( unsigned j = 0; j < _nCells; ++j ) {
if( _boundaryCells[j] ) continue;
// loop over all ordinates
for( unsigned k = 0; k < _nq; ++k ) {
psiNew[j][k] = 0.0;
......@@ -32,6 +32,10 @@ void SNSolver::Solve() {
_psi = psiNew;
// psiNew.reset();
}
std::vector<std::string> fieldNames;
fieldNames.push_back( "test" );
// ExportVTK( "test", _psi, fieldNames, _settings, _mesh );
}
void SNSolver::SolveMPI() {
......
#include "solver.h"
#include "../../include/quadratures/quadrature.h"
#include "io.h"
#include "mesh.h"
#include "../../include/quadratures/quadrature.h"
#include "snsolver.h"
Solver::Solver( Settings* settings ) {
......@@ -40,6 +40,12 @@ Solver::Solver( Settings* settings ) {
// TODO: setup initial condtion (distribution at initial energy for CSD)
SetupIC();
// setup numerical flux
_g = NumericalFlux::Create( settings );
// boundary type
_boundaryCells = mesh->GetBoundaryCellArray();
}
double Solver::ComputeTimeStep( double cfl ) const {
......@@ -72,11 +78,13 @@ void Solver::LoadStoppingPower( std::string fileName ) {
void Solver::LoadSigmaS( std::string fileName ) {
// @TODO
//_sigmaSH20 = ... -> dim(_sigmaSH20) = (_nTimeSteps,_nq,_nq)
_sigmaSH20 = std::vector( _nCells, Matrix( _nq, _nq, 1.0 ) ); // TODO: double check this!
}
void Solver::LoadSigmaT( std::string fileName ) {
// @TODO
//_sigmaTH20 = ... -> dim(_sigmaTH20) = _nTimeSteps
_sigmaTH20 = std::vector( _nCells, 1.0 );
}
void Solver::SetupIC() {
......
......@@ -3,10 +3,12 @@
UpwindFlux::UpwindFlux( Settings* settings ) : NumericalFlux( settings ) {}
double UpwindFlux::Flux( const Vector& Omega, double psiL, double psiR, const Vector& n ) const {
if( inner( Omega, n ) > 0 ) {
return inner( Omega, n ) * psiL;
// if( inner( Omega, n ) > 0 ) {
double inner = Omega[0] * n[0] + Omega[1] * n[1];
if( inner > 0 ) {
return inner * psiL;
}
else {
return inner( Omega, n ) * psiR;
return inner * psiR;
}
}
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