Commit 6243cd66 authored by steffen.schotthoefer's avatar steffen.schotthoefer
Browse files

reworked inheritance structure of solver

parent 4ea444ac
Pipeline #97323 passed with stages
in 43 minutes and 1 second
......@@ -10,7 +10,7 @@ class NewtonOptimizer : public OptimizerBase
inline ~NewtonOptimizer() {}
virtual Vector Solve( Vector u ) override;
void Solve( Vector& lambda, Vector& u ) override;
};
#endif // NEWTONOPTIMIZER_H
......@@ -19,7 +19,7 @@ class OptimizerBase
/*! @brief : Computes the optimal Lagrange multilpiers for the dual entropy minimization problem
* @param : Vector u = pointer to vector of given moments. // Maybe use pointer for performance?
* @return : Vector alpha = optimal lagrange multipliers. Has the same length as Vector u. */
virtual Vector Solve( Vector u ) = 0;
virtual void Solve( Vector& lambda, Vector& u ) = 0;
protected:
EntropyBase* _entropy; /*! @brief: Class to handle entropy functional evaluations */
......
#ifndef CSDSNSOLVER_H
#define CSDSNSOLVER_H
#include "solvers/solverbase.h"
#include "solvers/snsolver.h"
class CSDSNSolver : public Solver
class CSDSNSolver : public SNSolver
{
private:
std::vector<double> _dose;
Matrix _scatteringKernel; /*! @brief scattering kernel for the quadrature */
public:
/**
......
......@@ -37,7 +37,10 @@ class MNSolver : public Solver
Vector _scatterMatDiag; /*! @brief: Diagonal of the scattering matrix (its a diagonal matrix by construction) */
QuadratureBase* _quadrature; /*! @brief: Quadrature rule to compute flux Jacobian */ // Shift this to SolverBase!
// quadrature related numbers
VectorVector _quadPoints; /*! @brief quadrature points, dim(_quadPoints) = (_nq,spatialDim) */
Vector _weights; /*! @brief quadrature weights, dim(_weights) = (_nq) */
VectorVector _quadPointsSphere; /*! @brief (my,phi), dim(_quadPoints) = (_nq,2) */
EntropyBase* _entropy; /*! @brief: Class to handle entropy functionals */
VectorVector _alpha; /*! @brief: Lagrange Multipliers for Minimal Entropy problem for each gridCell
......
......@@ -5,9 +5,14 @@
class SNSolver : public Solver
{
private:
protected:
Matrix _scatteringKernel; /*! @brief scattering kernel for the quadrature */
// quadrature related numbers
VectorVector _quadPoints; /*! @brief quadrature points, dim(_quadPoints) = (_nq,spatialDim) */
Vector _weights; /*! @brief quadrature weights, dim(_weights) = (_nq) */
public:
/**
* @brief SNSolver constructor
......
......@@ -3,9 +3,9 @@
#include <mpi.h>
#include "solvers/solverbase.h"
#include "solvers/snsolver.h"
class SNSolverMPI : public Solver
class SNSolverMPI : public SNSolver
{
public:
/**
......
......@@ -33,9 +33,11 @@ class Solver
VectorVector _sigmaT; /*! @brief total cross section for all energies */
// quadrature related numbers
unsigned _nq; /*! @brief number of quadrature points */
VectorVector _quadPoints; /*! @brief quadrature points, dim(_quadPoints) = (_nSystem,spatialDim) */
Vector _weights; /*! @brief quadrature weights, dim(_weights) = (_NCells) */
QuadratureBase* _quadrature; /*! @brief quadrature to create members below */
unsigned _nq; /*! @brief number of quadrature points */
// VectorVector _quadPoints; /*! @brief quadrature points, dim(_quadPoints) = (_nSystem,spatialDim) */
// Vector _weights; /*! @brief quadrature weights, dim(_weights) = (_NCells) */
// Mesh related members
unsigned _nCells; /*! @brief number of spatial cells */
......
......@@ -2,10 +2,83 @@
#include "optimizers/newtonoptimizer.h"
#include "settings/config.h"
Vector NewtonOptimizer::Solve( Vector u ) {
void NewtonOptimizer::Solve( Vector& lambda, Vector& u ) {
// if we have quadratic entropy, then alpha = u;
if( _settings->GetEntropyName() == QUADRATIC ) return u;
if( _settings->GetEntropyName() == QUADRATIC ) return;
return Vector( 1, 0.0 ); // dummy
// Start Newton Algorithm
//
// int maxRefinements = 1000;
// unsigned nTotal = _nTotalForRef[refLevel];
//
// Vector g( _nStates * nTotal );
//
// // check if initial guess is good enough
// Gradient( g, lambda, u, refLevel );
// if( CalcNorm( g, nTotal ) < _settings->GetEpsilon() ) {
// return;
// }
// Matrix H( _nStates * nTotal, _nStates * nTotal );
// Vector dlambdaNew( _nStates * nTotal );
// // calculate initial Hessian and gradient
// Vector dlambda = -g;
// // std::cout << g << std::endl;
// Hessian( H, lambda, refLevel );
//
// // double eps = 0.00001;
// // lambda( 1, 0 ) = lambda( 1, 0 ) + eps;
//
// // Vector gEps( _nStates * nTotal );
//
// // Gradient( gEps, lambda, u, refLevel );
//
// // std::cout << "H_FD = " << ( gEps - g ) / eps << std::endl;
// // std::cout << "H = " << H << std::endl;
// // exit( EXIT_FAILURE );
//
// posv( H, g );
// if( _maxIterations == 1 ) {
// AddMatrixVectorToMatrix( lambda, -_alpha * g, lambda, nTotal );
// return;
// }
// Matrix lambdaNew( _nStates, nTotal );
// AddMatrixVectorToMatrix( lambda, -_alpha * g, lambdaNew, nTotal );
// Gradient( dlambdaNew, lambdaNew, u, refLevel );
// // perform Newton iterations
// for( unsigned l = 0; l < _maxIterations; ++l ) {
// double stepSize = 1.0;
// if( l != 0 ) {
// Gradient( g, lambda, u, refLevel );
// dlambda = -g;
// Hessian( H, lambda, refLevel );
// posv( H, g );
// AddMatrixVectorToMatrix( lambda, -stepSize * _alpha * g, lambdaNew, nTotal );
// Gradient( dlambdaNew, lambdaNew, u, refLevel );
// }
// int refinementCounter = 0;
// // std::cout << CalcNorm( dlambdaNew, nTotal ) << std::endl;
// while( CalcNorm( dlambda, nTotal ) < CalcNorm( dlambdaNew, nTotal ) || !std::isfinite( CalcNorm( dlambdaNew, nTotal ) ) ) {
// stepSize *= 0.5;
// AddMatrixVectorToMatrix( lambda, -stepSize * _alpha * g, lambdaNew, nTotal );
// Gradient( dlambdaNew, lambdaNew, u, refLevel );
// if( CalcNorm( dlambdaNew, nTotal ) < _settings->GetEpsilon() ) {
// lambda = lambdaNew;
// return;
// }
// else if( ++refinementCounter > maxRefinements ) {
// _log->error( "[closure] Newton needed too many refinement steps!" );
// exit( EXIT_FAILURE );
// }
// }
// lambda = lambdaNew;
// if( CalcNorm( dlambdaNew, nTotal ) < _settings->GetEpsilon() ) {
// lambda = lambdaNew;
// return;
// }
// }
// _log->error( "[closure] Newton did not converge!" );
// exit( EXIT_FAILURE );
// return Vector( 1, 0.0 ); // dummy
}
......@@ -5,15 +5,7 @@
#include "solvers/csdsnsolver.h"
CSDSNSolver::CSDSNSolver( Config* settings ) : Solver( settings ) {
_dose = std::vector<double>( _settings->GetNCells(), 0.0 );
QuadratureBase* quad = QuadratureBase::CreateQuadrature( settings->GetQuadName(), settings->GetQuadOrder() );
ScatteringKernel* k = ScatteringKernel::CreateScatteringKernel( settings->GetKernelName(), quad );
_scatteringKernel = k->GetScatteringKernel();
delete quad;
delete k;
}
CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) { _dose = std::vector<double>( _settings->GetNCells(), 0.0 ); }
void CSDSNSolver::Solve() {
auto log = spdlog::get( "event" );
......
......@@ -17,7 +17,13 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ) {
// Is this good (fast) code using a constructor list?
_nMaxMomentsOrder = settings->GetMaxMomentDegree();
_nTotalEntries = GlobalIndex( _nMaxMomentsOrder, int( _nMaxMomentsOrder ) ) + 1;
_quadrature = QuadratureBase::CreateQuadrature( _settings->GetQuadName(), settings->GetQuadOrder() );
// build quadrature object and store quadrature points and weights
_quadPoints = _quadrature->GetPoints();
_weights = _quadrature->GetWeights();
_nq = _quadrature->GetNq();
_quadPointsSphere = _quadrature->GetPointsSphere();
_settings->SetNQuadPoints( _nq );
// transform sigmaT and sigmaS in sigmaA.
_sigmaA = VectorVector( _nEnergies, Vector( _nCells, 0 ) ); // Get rid of this extra vektor!
......@@ -50,7 +56,6 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ) {
}
MNSolver::~MNSolver() {
delete _quadrature;
delete _entropy;
delete _optimizer;
delete _basis;
......@@ -66,8 +71,8 @@ void MNSolver::ComputeMoments() {
double my, phi;
for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
my = _quadrature->GetPointsSphere()[idx_quad][0];
phi = _quadrature->GetPointsSphere()[idx_quad][1];
my = _quadPointsSphere[idx_quad][0];
phi = _quadPointsSphere[idx_quad][1];
_moments[idx_quad] = _basis->ComputeSphericalBasis( my, phi );
}
......@@ -83,13 +88,13 @@ Vector MNSolver::ConstructFlux( unsigned idx_cell ) {
for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
w = _quadrature->GetWeights()[idx_quad];
w = _weights[idx_quad];
entropyFlux = 0.0; // Reset temorary flux
entropyL = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], _moments[idx_quad] ) );
omega = _quadrature->GetPoints()[idx_quad]; // Use Pointer!
omega = _quadPoints[idx_quad]; // Use Pointer!
for( unsigned idx_neigh = 0; idx_neigh < _neighbors[idx_cell].size(); idx_neigh++ ) {
// Store fluxes in psiNew, to save memory
......
......@@ -9,11 +9,12 @@
SNSolver::SNSolver( Config* settings ) : Solver( settings ) {
QuadratureBase* quad = QuadratureBase::CreateQuadrature( settings->GetQuadName(), settings->GetQuadOrder() );
ScatteringKernel* k = ScatteringKernel::CreateScatteringKernel( settings->GetKernelName(), quad );
_scatteringKernel = k->GetScatteringKernel();
_quadPoints = _quadrature->GetPoints();
_weights = _quadrature->GetWeights();
ScatteringKernel* k = ScatteringKernel::CreateScatteringKernel( settings->GetKernelName(), _quadrature );
_scatteringKernel = k->GetScatteringKernel();
delete k;
delete quad;
}
void SNSolver::Solve() {
......
......@@ -4,7 +4,7 @@
#include "solvers/snsolverMPI.h"
#include <mpi.h>
SNSolverMPI::SNSolverMPI( Config* settings ) : Solver( settings ) {}
SNSolverMPI::SNSolverMPI( Config* settings ) : SNSolver( settings ) {}
void SNSolverMPI::Solve() {
/*
......
......@@ -12,7 +12,7 @@
Solver::Solver( Config* settings ) : _settings( settings ) {
// @TODO save parameters from settings class
// build mesh and store all relevant information
// build mesh and store and store frequently used params
_mesh = LoadSU2MeshFromFile( settings );
_areas = _mesh->GetCellAreas();
_neighbors = _mesh->GetNeighbours();
......@@ -20,20 +20,17 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
_nCells = _mesh->GetNumCells();
_settings->SetNCells( _nCells );
// build quadrature object and store quadrature points and weights
QuadratureBase* quad = QuadratureBase::CreateQuadrature( settings->GetQuadName(), settings->GetQuadOrder() );
_quadPoints = quad->GetPoints();
_weights = quad->GetWeights();
_nq = quad->GetNq();
// build quadrature object and store frequently used params
_quadrature = QuadratureBase::CreateQuadrature( settings->GetQuadName(), settings->GetQuadOrder() );
_nq = _quadrature->GetNq();
_settings->SetNQuadPoints( _nq );
delete quad;
// set time step
_dE = ComputeTimeStep( settings->GetCFL() );
_nEnergies = unsigned( settings->GetTEnd() / _dE );
for( unsigned i = 0; i < _nEnergies; ++i ) _energies.push_back( ( i + 1 ) * _dE );
// setup problem
// setup problem and store frequently used params
_problem = ProblemBase::Create( _settings, _mesh );
_sol = _problem->SetupIC();
_s = _problem->GetStoppingPower( _energies );
......@@ -52,6 +49,7 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
}
Solver::~Solver() {
delete _quadrature;
delete _mesh;
delete _problem;
}
......
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