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

pulled master into nnDataGen 


Former-commit-id: 9d23f5b6
parents ddd36751 7590242d
Show whitespace changes
Inline Side-by-side
code/input/example_csd_1d.cfg
View file @ 66466fe1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checkerboard Benchmarking File SN %
% Author <Jonas Kusch> %
% Date 01.12.2020 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
OUTPUT_DIR = ../result
OUTPUT_FILE = example_csd_1d_10MeV_fine
LOG_DIR = ../result/logs
......
code/src/common/config.cpp
View file @ 66466fe1
......@@ -16,7 +16,6 @@
#include "spdlog/sinks/basic_file_sink.h"
#include "spdlog/sinks/stdout_sinks.h"
#include "spdlog/spdlog.h"
#include <cassert>
#include <filesystem>
#include <fstream>
#include <mpi.h>
......@@ -74,6 +73,12 @@ Config::Config( string case_filename ) {
Config::~Config( void ) {
// Delete all introduced arrays!
// delete _option map values proberly
for( auto const& x : _optionMap ) {
delete x.second;
_optionMap.erase( x.first );
}
}
// ---- Add Options ----
......@@ -500,6 +505,11 @@ void Config::SetPostprocessing() {
CURRENT_FUNCTION );
}
break;
case CSD_SN_NOTRAFO_SOLVER:
case CSD_SN_FOKKERPLANCK_SOLVER:
case CSD_SN_FOKKERPLANCK_TRAFO_SOLVER:
case CSD_SN_FOKKERPLANCK_TRAFO_SOLVER_2D:
case CSD_SN_FOKKERPLANCK_TRAFO_SH_SOLVER_2D:
case CSD_SN_SOLVER:
supportedGroups = { MINIMAL, DOSE };
if( supportedGroups.end() == std::find( supportedGroups.begin(), supportedGroups.end(), _volumeOutput[idx_volOutput] ) ) {
......
code/src/common/mesh.cpp
View file @ 66466fe1
......@@ -358,6 +358,9 @@ void Mesh::ComputePartitioning() {
void Mesh::ComputeSlopes( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const {
for( unsigned k = 0; k < nq; ++k ) {
for( unsigned j = 0; j < _numCells; ++j ) {
psiDerX[j][k] = 0.0;
psiDerY[j][k] = 0.0;
// if( cell->IsBoundaryCell() ) continue; // skip ghost cells
if( _cellBoundaryTypes[j] != 2 ) continue; // skip ghost cells
// compute derivative by summing over cell boundary
......@@ -400,13 +403,15 @@ void Mesh::ReconstructSlopesS( unsigned nq, VectorVector& psiDerX, VectorVector&
void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const {
double phi;
VectorVector dPsiMax = std::vector( _numCells, Vector( nq, 0.0 ) );
VectorVector dPsiMin = std::vector( _numCells, Vector( nq, 0.0 ) );
//VectorVector dPsiMax = std::vector( _numCells, Vector( nq, 0.0 ) );
//VectorVector dPsiMin = std::vector( _numCells, Vector( nq, 0.0 ) );
VectorVector dPsiMax( _numCells, Vector( nq, 0.0 ) );
VectorVector dPsiMin( _numCells, Vector( nq, 0.0 ) );
std::vector<std::vector<Vector>> psiSample( _numCells, std::vector<Vector>( nq, Vector( _numNodesPerCell, 0.0 ) ) );
std::vector<std::vector<Vector>> phiSample( _numCells, std::vector<Vector>( nq, Vector( _numNodesPerCell, 0.0 ) ) );
for( unsigned k = 0; k < nq; ++k ) {
for( unsigned j = 0; j < _numCells; ++j ) {
// reset derivatives
psiDerX[j][k] = 0.0;
......@@ -415,6 +420,8 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
// skip boundary cells
if( _cellBoundaryTypes[j] != 2 ) continue;
/*
// version 1: original taste
// step 1: calculate psi difference around neighbors and theoretical derivatives by Gauss theorem
for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
if( psi[_cellNeighbors[j][l]][k] - psi[j][k] > dPsiMax[j][k] ) dPsiMax[j][k] = psi[_cellNeighbors[j][l]][k] - psi[j][k];
......@@ -427,18 +434,19 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
// step 2: choose sample points
// psiSample[j][k][l] = 0.5 * ( psi[j][k] + psi[_cellNeighbors[j][l]][k] ); // interface central points
psiSample[j][k][l] = psi[j][k] + psiDerX[j][k] * ( _nodes[_cells[j][l]][0] - _cellMidPoints[j][0] ) +
psiSample[j][k][l] = psi[j][k] +
psiDerX[j][k] * ( _nodes[_cells[j][l]][0] - _cellMidPoints[j][0] ) +
psiDerY[j][k] * ( _nodes[_cells[j][l]][1] - _cellMidPoints[j][1] ); // vertex points
// step 3: calculate Phi_ij at sample points
if( psiSample[j][k][l] > psi[j][k] ) {
phiSample[j][k][l] = fmin( 0.5, dPsiMax[j][k] / ( psiSample[j][k][l] - psi[j][k] ) );
phiSample[j][k][l] = fmin( 1.0, dPsiMax[j][k] / ( psiSample[j][k][l] - psi[j][k] ) );
}
else if( psiSample[j][l][k] < psi[j][k] ) {
phiSample[j][k][l] = fmin( 0.5, dPsiMin[j][k] / ( psiSample[j][k][l] - psi[j][k] ) );
phiSample[j][k][l] = fmin( 1.0, dPsiMin[j][k] / ( psiSample[j][k][l] - psi[j][k] ) );
}
else {
phiSample[j][k][l] = 0.5;
phiSample[j][k][l] = 1.0;
}
}
......@@ -446,12 +454,53 @@ void Mesh::ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector&
phi = min( phiSample[j][k] );
// step 5: limit the slope reconstructed from Gauss theorem
for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
psiDerX[j][k] *= phi;
psiDerY[j][k] *= phi;
*/
double eps = 1e-6;
// version 2: Venkatakrishnan limiter
// step 1: calculate psi difference around neighbors and theoretical derivatives by Gauss theorem
for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
if( psi[_cellNeighbors[j][l]][k] - psi[j][k] > dPsiMax[j][k] ) dPsiMax[j][k] = psi[_cellNeighbors[j][l]][k] - psi[j][k];
if( psi[_cellNeighbors[j][l]][k] - psi[j][k] < dPsiMin[j][k] ) dPsiMin[j][k] = psi[_cellNeighbors[j][l]][k] - psi[j][k];
psiDerX[j][k] += 0.5 * ( psi[j][k] + psi[_cellNeighbors[j][l]][k] ) * _cellNormals[j][l][0] / _cellAreas[j];
psiDerY[j][k] += 0.5 * ( psi[j][k] + psi[_cellNeighbors[j][l]][k] ) * _cellNormals[j][l][1] / _cellAreas[j];
}
for( unsigned l = 0; l < _cellNeighbors[j].size(); ++l ) {
// step 2: choose sample points
psiSample[j][k][l] = 10.0 * psiDerX[j][k] * (_cellMidPoints[_cellNeighbors[j][l]][0] - _cellMidPoints[j][0]) +
10.0 * psiDerY[j][k] * (_cellMidPoints[_cellNeighbors[j][l]][1] - _cellMidPoints[j][1]);
// step 3: calculate Phi_ij at sample points
if( psiSample[j][k][l] > 0.0 ) {
phiSample[j][k][l] =
(dPsiMax[j][k]*dPsiMax[j][k] + 2.0 * dPsiMax[j][k] * psiSample[j][k][l] + eps)/
(dPsiMax[j][k]*dPsiMax[j][k] + dPsiMax[j][k] * psiSample[j][k][l] + 2.0 * psiSample[j][k][l] * psiSample[j][k][l] + eps);
}
else if( psiSample[j][k][l] < 0.0 ) {
phiSample[j][k][l] =
(dPsiMin[j][k]*dPsiMin[j][k] + 2.0 * dPsiMin[j][k] * psiSample[j][k][l] + eps)/
(dPsiMin[j][k]*dPsiMin[j][k] + dPsiMin[j][k] * psiSample[j][k][l] + 2.0 * psiSample[j][k][l] * psiSample[j][k][l] + eps);;
}
else {
phiSample[j][k][l] = 1.0;
}
}
// step 4: find minimum limiter function phi
phi = min( phiSample[j][k] );
//phi = fmin( 0.5, abs(phi) );
// step 5: limit the slope reconstructed from Gauss theorem
psiDerX[j][k] *= phi;
psiDerY[j][k] *= phi;
}
}
}
void Mesh::ComputeBounds() {
......
code/src/main.cpp
View file @ 66466fe1
/*! @file: main.cpp
* @brief: Main method to call the KiT-RT solver suite
* @author: J. Kusch, S. Schotthöfer, P. Stammer, J. Wolters, T. Xiao
* @version: 0.1
*/
#include <Python.h>
#include <mpi.h>
#include <string>
#include "common/config.h"
#include "common/io.h"
#include "solvers/solverbase.h"
#include "common/config.h"
#include "solvers/sphericalharmonics.h"
#include <fstream>
#include <iostream>
#include <string>
// ----
#include "optimizers/optimizerbase.h"
#include "quadratures/qgausslegendretensorized.h"
#include "quadratures/qmontecarlo.h"
#include "solvers/sphericalharmonics.h"
#include "toolboxes/datagenerator.h"
double testFunc( double my, double phi ) { return my * my + phi; }
double testFunc2( double x, double y, double z ) { return x + y + z; }
int main( int argc, char** argv ) {
MPI_Init( &argc, &argv );
wchar_t* program = Py_DecodeLocale( argv[0], NULL );
......@@ -32,23 +21,17 @@ int main( int argc, char** argv ) {
// CD Load Settings from File
Config* config = new Config( filename );
// Print input file and run info to file
PrintLogHeader( filename );
if( config->GetDataGeneratorMode() ) {
// Build Data generator
nnDataGenerator* datagen = new nnDataGenerator( config );
// Generate Data and export
datagen->computeTrainingData();
}
else {
// Build solver
Solver* solver = Solver::Create( config );
// Run solver and export
solver->Solve();
solver->PrintVolumeOutput();
}
MPI_Finalize();
return EXIT_SUCCESS;
}
code/src/optimizers/newtonoptimizer.cpp
View file @ 66466fe1
......@@ -6,7 +6,7 @@
#include "toolboxes/errormessages.h"
NewtonOptimizer::NewtonOptimizer( Config* settings ) : OptimizerBase( settings ) {
_quadrature = QuadratureBase::CreateQuadrature( settings );
_quadrature = QuadratureBase::Create( settings );
_nq = _quadrature->GetNq();
_weights = _quadrature->GetWeights();
_quadPointsSphere = _quadrature->GetPointsSphere();
......
code/src/problems/checkerboard.cpp
View file @ 66466fe1
......@@ -27,7 +27,7 @@ VectorVector Checkerboard_SN::GetScatteringXS( const Vector& energies ) { return
VectorVector Checkerboard_SN::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), _totalXS ); }
std::vector<VectorVector> Checkerboard_SN::GetExternalSource( const Vector& energies ) {
std::vector<VectorVector> Checkerboard_SN::GetExternalSource( const Vector& /*energies*/ ) {
VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
auto cellMids = _mesh->GetCellMidPoints();
for( unsigned j = 0; j < cellMids.size(); ++j ) {
......@@ -91,7 +91,7 @@ VectorVector Checkerboard_PN::GetScatteringXS( const Vector& energies ) { return
VectorVector Checkerboard_PN::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), _totalXS ); }
std::vector<VectorVector> Checkerboard_PN::GetExternalSource( const Vector& energies ) {
std::vector<VectorVector> Checkerboard_PN::GetExternalSource( const Vector& /*energies*/ ) {
VectorVector Q( _mesh->GetNumCells(), Vector( 1u, 0.0 ) );
auto cellMids = _mesh->GetCellMidPoints();
for( unsigned j = 0; j < cellMids.size(); ++j ) {
......
code/src/icru.cpp → code/src/problems/icru.cpp
View file @ 66466fe1
#include "icru.h"
#include "problems/icru.h"
void ICRU::angdcs( unsigned ncor, double e, std::vector<double>& dxse, std::vector<double>& dxsi, std::vector<double>& dxs ) {
......
code/src/problems/musclebonelung.cpp
View file @ 66466fe1
......@@ -21,12 +21,12 @@ VectorVector MuscleBoneLung::SetupIC() {
return psi;
}
std::vector<double> MuscleBoneLung::GetDensity( const VectorVector& cellMidPoints ) {
std::vector<double> MuscleBoneLung::GetDensity( const VectorVector& /*cellMidPoints*/ ) {
std::vector<double> densities( 167, 1.04 ); // muscle layer
std::vector<double> bone( 167, 1.85 );
std::vector<double> lung( 665, 0.3 ); // maybe this is just the lung tissue and after that there should be air?
//Concatenate layers
// Concatenate layers
densities.insert( densities.end(), bone.begin(), bone.end() );
densities.insert( densities.end(), lung.begin(), lung.end() );
return densities;
......
code/src/problems/slabgeohg.cpp
View file @ 66466fe1
......@@ -7,21 +7,21 @@
SlabGeoHG::SlabGeoHG( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) { _physics = nullptr; }
SlabGeoHG::~SlabGeoHG(){};
SlabGeoHG::~SlabGeoHG() {}
VectorVector SlabGeoHG::GetScatteringXS( const Vector& energies ) { return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) ); }
VectorVector SlabGeoHG::GetTotalXS( const Vector& energies ) { return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) ); }
std::vector<VectorVector> SlabGeoHG::GetExternalSource( const Vector& energies ) {
std::vector<VectorVector> SlabGeoHG::GetExternalSource( const Vector& /*energies*/ ) {
return std::vector<VectorVector>( 1u, std::vector<Vector>( _mesh->GetNumCells(), Vector( 1u, 0.0 ) ) );
}
VectorVector SlabGeoHG::SetupIC() {
VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) );
auto boundaryCells = _mesh->GetBoundaryTypes();
auto cellMids = _mesh->GetCellMidPoints();
double t = 3.2e-4; // pseudo time for gaussian smoothing
// auto boundaryCells = _mesh->GetBoundaryTypes();
// auto cellMids = _mesh->GetCellMidPoints();
// double t = 3.2e-4; // pseudo time for gaussian smoothing
return psi;
}
code/src/quadratures/productquadrature.cpp → code/src/quadratures/qproduct.cpp
View file @ 66466fe1
#include "quadratures/productquadrature.h"
#include "quadratures/qproduct.h"
#include "toolboxes/errormessages.h"
ProductQuadrature::ProductQuadrature( unsigned order ) : QuadratureBase( order ) {
QProduct::QProduct( unsigned order ) : QuadratureBase( order ) {
SetName();
SetNq();
SetPointsAndWeights();
SetConnectivity();
}
ProductQuadrature::ProductQuadrature( Config* settings ) : QuadratureBase( settings ) {
QProduct::QProduct( Config* settings ) : QuadratureBase( settings ) {
SetName();
SetNq();
SetPointsAndWeights();
SetConnectivity();
}
void ProductQuadrature::SetPointsAndWeights() {
void QProduct::SetPointsAndWeights() {
Vector nodes1D( 2 * _order ), weights1D( 2 * _order );
// construct companion matrix
......@@ -82,13 +82,13 @@ void ProductQuadrature::SetPointsAndWeights() {
}
}
void ProductQuadrature::SetConnectivity() { // TODO
void QProduct::SetConnectivity() { // TODO
// Not initialized for this quadrature.
VectorVectorU connectivity;
_connectivity = connectivity;
}
std::pair<Vector, Matrix> ProductQuadrature::ComputeEigenValTriDiagMatrix( const Matrix& mat ) {
std::pair<Vector, Matrix> QProduct::ComputeEigenValTriDiagMatrix( const Matrix& mat ) {
// copied from 'Numerical Recipes' and updated + modified to work with blaze
unsigned n = mat.rows();
......@@ -155,14 +155,14 @@ std::pair<Vector, Matrix> ProductQuadrature::ComputeEigenValTriDiagMatrix( const
return std::make_pair( d, z );
}
double ProductQuadrature::Pythag( const double a, const double b ) {
double QProduct::Pythag( const double a, const double b ) {
// copied from 'Numerical Recipes'
double absa = std::fabs( a ), absb = std::fabs( b );
return ( absa > absb ? absa * std::sqrt( 1.0 + ( absb / absa ) * ( absb / absa ) )
: ( absb == 0.0 ? 0.0 : absb * std::sqrt( 1.0 + ( absa / absb ) * ( absa / absb ) ) ) );
}
bool ProductQuadrature::CheckOrder() {
bool QProduct::CheckOrder() {
if( _order % 2 == 1 ) { // order needs to be even
ErrorMessages::Error( "ERROR! Order " + std::to_string( _order ) + " for " + GetName() + " not available. \n Order must be an even number. ",
CURRENT_FUNCTION );
......
code/src/quadratures/quadraturebase.cpp
View file @ 66466fe1
#include "quadratures/quadraturebase.h"
#include "common/config.h"
#include "quadratures/productquadrature.h"
#include "quadratures/qgausslegendre1D.h"
#include "quadratures/qgausslegendretensorized.h"
#include "quadratures/qldfesa.h"
#include "quadratures/qlebedev.h"
#include "quadratures/qlevelsymmetric.h"
#include "quadratures/qmontecarlo.h"
#include "quadratures/qproduct.h"
#include "toolboxes/errormessages.h"
QuadratureBase::QuadratureBase( Config* settings ) {
......@@ -16,7 +16,7 @@ QuadratureBase::QuadratureBase( Config* settings ) {
QuadratureBase::QuadratureBase( unsigned order ) : _order( order ) { _settings = nullptr; }
QuadratureBase* QuadratureBase::CreateQuadrature( Config* settings ) {
QuadratureBase* QuadratureBase::Create( Config* settings ) {
QUAD_NAME name = settings->GetQuadName();
switch( name ) {
......@@ -26,12 +26,12 @@ QuadratureBase* QuadratureBase::CreateQuadrature( Config* settings ) {
case QUAD_LevelSymmetric: return new QLevelSymmetric( settings );
case QUAD_LDFESA: return new QLDFESA( settings );
case QUAD_Lebedev: return new QLebedev( settings );
case QUAD_Product: return new ProductQuadrature( settings );
case QUAD_Product: return new QProduct( settings );
default: ErrorMessages::Error( "Creator for the chose quadrature does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
}
}
QuadratureBase* QuadratureBase::CreateQuadrature( QUAD_NAME name, unsigned quadOrder ) {
QuadratureBase* QuadratureBase::Create( QUAD_NAME name, unsigned quadOrder ) {
switch( name ) {
case QUAD_MonteCarlo: return new QMonteCarlo( quadOrder );
......@@ -42,7 +42,7 @@ QuadratureBase* QuadratureBase::CreateQuadrature( QUAD_NAME name, unsigned quadO
case QUAD_LevelSymmetric: return new QLevelSymmetric( quadOrder );
case QUAD_LDFESA: return new QLDFESA( quadOrder );
case QUAD_Lebedev: return new QLebedev( quadOrder );
case QUAD_Product: return new ProductQuadrature( quadOrder );
case QUAD_Product: return new QProduct( quadOrder );
default: ErrorMessages::Error( "Creator for the chose quadrature does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
}
}
......
code/src/solvers/csdsnsolver.cpp
View file @ 66466fe1
......@@ -3,6 +3,7 @@
#include "common/io.h"
#include "fluxes/numericalflux.h"
#include "kernels/scatteringkernelbase.h"
#include "problems/icru.h"
#include "problems/problembase.h"
// externals
......
code/src/solvers/csdsnsolverfp.cpp
View file @ 66466fe1
......@@ -4,6 +4,7 @@
#include "common/mesh.h"
#include "fluxes/numericalflux.h"
#include "kernels/scatteringkernelbase.h"
#include "problems/icru.h"
#include "problems/problembase.h"
#include "quadratures/quadraturebase.h"
......@@ -32,7 +33,7 @@ CSDSNSolverFP::CSDSNSolverFP( Config* settings ) : SNSolver( settings ) {
// create 1D quadrature
unsigned nq = _settings->GetNQuadPoints();
QuadratureBase* quad1D = QuadratureBase::CreateQuadrature( QUAD_GaussLegendre1D, nq );
QuadratureBase* quad1D = QuadratureBase::Create( QUAD_GaussLegendre1D, nq );
Vector w = quad1D->GetWeights();
VectorVector muVec = quad1D->GetPoints();
Vector mu( nq );
......@@ -58,9 +59,6 @@ CSDSNSolverFP::CSDSNSolverFP( Config* settings ) : SNSolver( settings ) {
}
}
// Heney-Greenstein parameter
double g = 0.8;
// determine momente of Heney-Greenstein
_xi = Matrix( 3, _nEnergies );
......
code/src/solvers/csdsnsolvernotrafo.cpp
View file @ 66466fe1
......@@ -3,6 +3,7 @@
#include "common/io.h"
#include "fluxes/numericalflux.h"
#include "kernels/scatteringkernelbase.h"
#include "problems/icru.h"
#include "problems/problembase.h"
// externals
......
code/src/solvers/csdsolvertrafofp.cpp
View file @ 66466fe1
......@@ -3,6 +3,7 @@
#include "common/io.h"
#include "fluxes/numericalflux.h"
#include "kernels/scatteringkernelbase.h"
#include "problems/icru.h"
#include "problems/problembase.h"
#include "quadratures/quadraturebase.h"
......@@ -23,7 +24,7 @@ CSDSolverTrafoFP::CSDSolverTrafoFP( Config* settings ) : SNSolver( settings ) {
// create 1D quadrature
unsigned nq = _settings->GetNQuadPoints();
QuadratureBase* quad1D = QuadratureBase::CreateQuadrature( QUAD_GaussLegendre1D, nq );
QuadratureBase* quad1D = QuadratureBase::Create( QUAD_GaussLegendre1D, nq );
Vector w = quad1D->GetWeights();
VectorVector muVec = quad1D->GetPoints();
Vector mu( nq );
......
code/src/solvers/csdsolvertrafofp2d.cpp
View file @ 66466fe1
......@@ -3,8 +3,9 @@
#include "common/io.h"
#include "fluxes/numericalflux.h"
#include "kernels/scatteringkernelbase.h"
#include "problems/icru.h"
#include "problems/problembase.h"
#include "quadratures/productquadrature.h"
#include "quadratures/qproduct.h"
#include "quadratures/quadraturebase.h"
// externals
......@@ -36,7 +37,7 @@ CSDSolverTrafoFP2D::CSDSolverTrafoFP2D( Config* settings ) : SNSolver( settings
_wa = Vector( 2 * order );
// create quadrature 1D to compute mu grid
QuadratureBase* quad1D = QuadratureBase::CreateQuadrature( QUAD_GaussLegendre1D, 2 * order );
QuadratureBase* quad1D = QuadratureBase::Create( QUAD_GaussLegendre1D, 2 * order );
Vector w = quad1D->GetWeights();
VectorVector muVec = quad1D->GetPoints();
for( unsigned k = 0; k < 2 * order; ++k ) {
......
code/src/solvers/csdsolvertrafofpsh2d.cpp
View file @ 66466fe1
#include "solvers/csdsolvertrafofpsh2d.h"
#include "problems/icru.h"
CSDSolverTrafoFPSH2D::CSDSolverTrafoFPSH2D( Config* settings ) : SNSolver( settings ) {
_dose = std::vector<double>( _settings->GetNCells(), 0.0 );
......@@ -13,7 +14,6 @@ CSDSolverTrafoFPSH2D::CSDSolverTrafoFPSH2D( Config* settings ) : SNSolver( setti
// create quadrature
unsigned order = _quadrature->GetOrder();
unsigned nq = _settings->GetNQuadPoints();
_quadPoints = _quadrature->GetPoints();
_weights = _quadrature->GetWeights();
_quadPointsSphere = _quadrature->GetPointsSphere();
......@@ -25,7 +25,7 @@ CSDSolverTrafoFPSH2D::CSDSolverTrafoFPSH2D( Config* settings ) : SNSolver( setti
_wa = Vector( 2 * order );
// create quadrature 1D to compute mu grid
QuadratureBase* quad1D = QuadratureBase::CreateQuadrature( QUAD_GaussLegendre1D, 2 * order );
QuadratureBase* quad1D = QuadratureBase::Create( QUAD_GaussLegendre1D, 2 * order );
Vector w = quad1D->GetWeights();
VectorVector muVec = quad1D->GetPoints();
for( unsigned k = 0; k < 2 * order; ++k ) {
......
code/src/solvers/mnsolver.cpp
View file @ 66466fe1
......@@ -7,8 +7,8 @@
#include "optimizers/optimizerbase.h"
#include "problems/problembase.h"
#include "quadratures/quadraturebase.h"
#include "solvers/sphericalharmonics.h"
#include "toolboxes/errormessages.h"
#include "toolboxes/sphericalharmonics.h"
#include "toolboxes/textprocessingtoolbox.h"
// externals
......@@ -28,6 +28,10 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ) {
_quadPointsSphere = _quadrature->GetPointsSphere();
//_settings->SetNQuadPoints( _nq );
// Initialize temporary storages of alpha derivatives
_solDx = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
_solDy = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
// Initialize Scatter Matrix --
_scatterMatDiag = Vector( _nTotalEntries, 0.0 );
ComputeScatterMatrix();
......@@ -85,29 +89,49 @@ void MNSolver::ComputeMoments() {
Vector MNSolver::ConstructFlux( unsigned idx_cell ) {
// ---- Integration of Moment of flux ----
//--- Integration of moments of flux ---
double entropyL, entropyR, entropyFlux;
Vector flux( _nTotalEntries, 0.0 );
for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
//--- Temporary storages of reconstructed alpha ---
Vector alphaL( _nTotalEntries, 0.0 );
Vector alphaR( _nTotalEntries, 0.0 );
entropyFlux = 0.0; // Reset temorary flux
entropyL = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], _moments[idx_quad] ) );
for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
entropyFlux = 0.0; // reset temorary flux
for( unsigned idx_neigh = 0; idx_neigh < _neighbors[idx_cell].size(); idx_neigh++ ) {
// Store fluxes in psiNew, to save memory
// Left side reconstruction
if( _reconsOrder > 1 ) {
alphaL = _alpha[idx_cell] +
_solDx[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_neigh][0] - _cellMidPoints[idx_cell][0] ) +
_solDy[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_neigh][1] - _cellMidPoints[idx_cell][1] );
}
else {
alphaL = _alpha[idx_cell];
}
entropyL = _entropy->EntropyPrimeDual( blaze::dot( alphaL, _moments[idx_quad] ) );
// Right side reconstruction
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_neigh] == _nCells )
entropyR = entropyL;
else {
entropyR = _entropy->EntropyPrimeDual( blaze::dot( _alpha[_neighbors[idx_cell][idx_neigh]], _moments[idx_quad] ) );
if( _reconsOrder > 1 ) {
alphaR = _alpha[_neighbors[idx_cell][idx_neigh]] +
_solDx[_neighbors[idx_cell][idx_neigh]] * ( _interfaceMidPoints[idx_cell][idx_neigh][0] - _cellMidPoints[_neighbors[idx_cell][idx_neigh]][0] ) +
_solDy[_neighbors[idx_cell][idx_neigh]] * ( _interfaceMidPoints[idx_cell][idx_neigh][1] - _cellMidPoints[_neighbors[idx_cell][idx_neigh]][1] );
}
else {
alphaR = _alpha[_neighbors[idx_cell][idx_neigh]];
}
entropyR = _entropy->EntropyPrimeDual( blaze::dot( alphaR, _moments[idx_quad] ) );
}
// Entropy flux
entropyFlux += _g->Flux( _quadPoints[idx_quad], entropyL, entropyR, _normals[idx_cell][idx_neigh] );
}
// Solution flux
flux += _moments[idx_quad] * ( _weights[idx_quad] * entropyFlux );
// ------- Relizablity Reconstruction Step ----
}
return flux;
}
......@@ -150,6 +174,10 @@ void MNSolver::ComputeRadFlux() {
}
void MNSolver::FluxUpdate() {
if( _reconsOrder > 1 ) {
_mesh-&