Commit 570ea385 authored by Jonas Kusch's avatar Jonas Kusch
Browse files

solver structure added

parent a5a5b102
#ifndef LAXFRIEDRICHSFLUX_H
#define LAXFRIEDRICHSFLUX_H
#include "typedef.h"
#include "numericalflux.h"
class LaxFriedrichsFlux : public NumericalFlux
{
double _dt;
public:
/**
* @brief LaxFriedrichsFlux
* @param settings
*/
LaxFriedrichsFlux(Settings* settings);
/**
* @brief Flux computes flux on edge for fixed ordinate at a given edge
* @param Omega fixed ordinate for flux computation
* @param psiL left solution state
* @param psiR right solution state
* @param n scaled normal vector of given edge
* @return numerical flux value
*/
virtual double Flux(const Vector& Omega, double psiL, double psiR, const Vector& n)const;
};
#endif // LAXFRIEDRICHSFLUX_H
#ifndef NUMERICALFLUX_H
#define NUMERICALFLUX_H
#include "typedef.h"
class NumericalFlux
{
public:
NumericalFlux(Settings* settings);
static NumericalFlux* Create( Settings* settings );
/**
* @brief Flux computes flux on edge for fixed ordinate at a given edge
* @param Omega fixed ordinate for flux computation
* @param psiL left solution state
* @param psiR right solution state
* @param n scaled normal vector of given edge
* @return numerical flux value
*/
virtual double Flux(const Vector& Omega, double psiL, double psiR, const Vector& n)const = 0;
};
#endif // NUMERICALFLUX_H
#ifndef SNSOLVER_H
#define SNSOLVER_H
#include "solver.h"
#include "typedef.h"
class SNSolver : public Solver
{
public:
/**
* @brief SNSolver constructor
* @param settings stores all needed information
*/
SNSolver(Settings* settings);
/**
* @brief Solve functions runs main time loop
*/
virtual void Solve();
};
#endif // SNSOLVER_H
#ifndef SOLVER_H
#define SOLVER_H
// include Matrix, Vector definitions
#include "typedef.h"
#include "numericalflux.h"
#include <string>
class Solver
{
private:
const unsigned _Q; // number of quadrature points
const unsigned _NCells; // number of spatial cells
const unsigned _NTimeSteps; // number of time steps, number of nodal energy values for CSD
Matrix _psi; // angular flux vector, dim(_psi) = (_NCells,_Q)
std::vector<unsigned> _areas; // surface area of all spatial cells, dim(_areas) = _NCells
std::vector<Vector> _normals; // edge normals multiplied by edge length, dim(_normals) = (_NCells,spatialDim)
std::vector<std::vector<unsigned>> _neighbors; // edge normals multiplied by edge length, dim(_neighbors) = (_NCells,nEdgesPerCell)
std::vector<double> _sH20; // stopping power H2O, dim(_sH20) = _nTimeSteps
std::vector<double> _density; // patient density, dim(_density) = _nCells
std::vector<Vector> _quadPoints; // quadrature points, dim(_quadPoints) = (_NCells,spatialDim)
std::vector<double> _weights; // quadrature weights, dim(_weights) = (_NCells)
// we will have to add a further dimension for quadPoints and weights once we start with multilevel SN
NumericalFlux* _g; // numerical flux function
/**
* @brief LoadPatientDensity loads density of patient from MRT/CT scan and saves it in _density
* @param fileName is name of patient file
*/
void LoadPatientDensity(std::string fileName);
/**
* @brief LoadStoppingPower loads stopping power of H20 and saves it in _sH20
* @param fileName is name of stopping power file
*/
void LoadStoppingPower(std::string fileName);
public:
/**
* @brief Solver constructor
* @param settings stores all needed information
*/
Solver( Settings* settings );
/**
* @brief Create constructor
* @param settings stores all needed information
* @return pointer to Solver
*/
static Solver* Create( Settings* settings );
/**
* @brief Solve functions runs main time loop
*/
virtual void Solve( ) = 0;
};
#endif // SOLVER_H
#ifndef TYPEDEFS_H
#define TYPEDEFS_H
#include <vector>
typedef std::vector<std::vector<double>> Matrix;
typedef std::vector<double> Vector;
#endif // TYPEDEFS_H
#ifndef UPWINDFLUX_H
#define UPWINDFLUX_H
#include "typedef.h"
#include "numericalflux.h"
class UpwindFlux : public NumericalFlux
{
public:
/**
* @brief UpwindFlux
* @param settings
*/
UpwindFlux(Settings* settings);
/**
* @brief Flux computes flux on edge for fixed ordinate at a given edge
* @param Omega fixed ordinate for flux computation
* @param psiL left solution state
* @param psiR right solution state
* @param n scaled normal vector of given edge
* @return numerical flux value
*/
virtual double Flux(const Vector& Omega, double psiL, double psiR, const Vector& n)const;
};
#endif // UPWINDFLUX_H
#include "laxfriedrichsflux.h"
LaxFriedrichsFlux::LaxFriedrichsFlux(Settings* settings) : NumericalFlux(settings)
{
}
double LaxFriedrichsFlux::Flux(const Vector& Omega, double psiL, double psiR, const Vector& n)const{
//double normN = norm(n);
//return 0.5 * Omega * ( psiL + psiR ) * n/normN - 0.5 * ( psiR - psiL ) * norm( n ) / _dt;
return -1.0;
}
#include "numericalflux.h"
#include "upwindflux.h"
NumericalFlux::NumericalFlux( Settings* settings ) {}
NumericalFlux* NumericalFlux::Create( Settings* settings ) { return new UpwindFlux( settings ); }
#include "snsolver.h"
SNSolver::SNSolver( Settings* settings ) : Solver(settings)
{
}
void SNSolver::Solve(){
// TODO main SN time loop
}
#include "solver.h"
#include "snsolver.h"
Solver::Solver( Settings* settings ) : _Q( 1 ), _NCells( 1 ), _NTimeSteps( 1 ) {
// @TODO save parameters from settings class
// @TODO create object mesh and store _cells, ...
// @TODO create object quadrature and store _quadPoints, ...
// @TODO create object quadrature and store _quadPoints, ...
// load density and stopping power data
LoadPatientDensity( "test" );
LoadStoppingPower( "test" );
}
void Solver::LoadPatientDensity( std::string fileName ) {
// @TODO
}
void Solver::LoadStoppingPower( std::string fileName ) {
// @TODO
}
Solver* Solver::Create( Settings* settings ) { return new SNSolver( settings ); }
#include "upwindflux.h"
UpwindFlux::UpwindFlux(Settings* settings) : NumericalFlux(settings)
{
}
double UpwindFlux::Flux(const Vector& Omega, double psiL, double psiR, const Vector& n)const{
/*if( Omega * n > 0 ) {
return Omega*psiL * norm(n);
}
else {
return Omega*psiR * norm(n);
}*/
return -1.0;
}
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