PN Runs, needs testing

code/include/mesh.h

@@ -66,6 +66,7 @@ class Mesh
 
     const std::vector<BOUNDARY_TYPE>& GetBoundaryTypes() const;
 
+    double GetDistanceToOrigin( unsigned idx_cell ) const;
     /**
      * @brief ComputeSlopes calculates the slope in every cell into x and y direction
      * @param nq is number of quadrature points

code/include/numericalflux.h

@@ -44,6 +44,17 @@ class NumericalFlux
                        const Vector&,
                        const Vector& n,
                        Vector& resultFlux ) const = 0;
+
+    virtual void FluxVanLeer( const Matrix& Ax,
+                              const Matrix& AxAbs,
+                              const Matrix& Ay,
+                              const Matrix& AyAbs,
+                              const Matrix& Az,
+                              const Matrix& AzAbs,
+                              const Vector& psiL,
+                              const Vector& psiR,
+                              const Vector& n,
+                              Vector& resultFlux ) const = 0;
 };
 
 #endif    // NUMERICALFLUX_H

code/include/pnsolver.h

@@ -25,12 +25,13 @@ class PNSolver : public Solver
      */
     void Save() const override;
 
+    void Save( int currEnergy ) const;
+
   protected:
     // moment orders for P_N
     unsigned _nTotalEntries;    // total number of equations in the system
 
     VectorVector _sigmaA;    // absorbtion coefficient for all energies
-
     // System Matrix for x, y and z flux
     //    ==> not needed after computation of A+ and A- ==> maybe safe only temporarly and remove as member?
     SymMatrix _Ax;
@@ -40,10 +41,13 @@ class PNSolver : public Solver
     // Upwinding Matrices
     Matrix _AxPlus;
     Matrix _AxMinus;
+    Matrix _AxAbs;
     Matrix _AyPlus;
     Matrix _AyMinus;
+    Matrix _AyAbs;
     Matrix _AzPlus;
     Matrix _AzMinus;
+    Matrix _AzAbs;
 
     Vector _scatterMatDiag;    // diagonal of the scattering matrix (its a diagonal matrix by construction)
 
@@ -64,8 +68,19 @@ class PNSolver : public Solver
     int Sgn( int k ) const;
     int kPlus( int k ) const;
     int kMinus( int k ) const;
+
+    /*! @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 : Checks, if index invariant for global index holds
+     *  @returns : True, if invariant holds, else false
+     */
+    bool CheckIndex( int l, int k ) const;
+
     // function for computing and setting up system matrices
     void ComputeSystemMatrices();
     // function for computing and setting up flux matrices for upwinding

code/include/snsolver.h

@@ -25,6 +25,7 @@ class SNSolver : public Solver
      * @brief Output solution to VTK file
      */
     virtual void Save() const;
+    virtual void Save( int currEnergy ) const;
 };
 
 #endif    // SNSOLVER_H

code/include/solver.h

@@ -33,6 +33,8 @@ class Solver
     VectorVector _quadPoints;                         // quadrature points, dim(_quadPoints) = (_nTimeSteps,spatialDim)
     Vector _weights;                                  // quadrature weights, dim(_weights) = (_NCells)
     std::vector<BOUNDARY_TYPE> _boundaryCells;        // boundary type for all cells, dim(_boundary) = (_NCells)
+    std::vector<double> _solverOutput;                // PROTOTYPE: Outputfield for solver
+
     // we will have to add a further dimension for quadPoints and weights once we start with multilevel SN
 
     NumericalFlux* _g;    // numerical flux function
@@ -69,6 +71,8 @@ class Solver
      * @brief Output solution to VTK file
      */
     virtual void Save() const = 0;
+
+    virtual void Save( int currEnergy ) const = 0;
 };
 
 #endif    // SOLVER_H

code/include/upwindflux.h

@@ -48,6 +48,32 @@ class UpwindFlux : public NumericalFlux
                const Vector& psiR,
                const Vector& n,
                Vector& resultFlux ) const override;
+
+    /**
+     * @brief Flux      : Computes <VanLeer> 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
+     * @param AyPlus    : Positive part of the flux jacobian in y direction
+     * @param AyMinus   : Negative part of the flux jacobian in y direction
+     * @param AzPlus    : Positive part of the flux jacobian in z direction
+     * @param AzMinus   : Negative part of the flux jacobian in z direction
+     * @param psiL      : Solution state of left hand side control volume
+     * @param psiR      : Solution state of right hand side control volume
+     * @param n         : Normal vector at the edge between left and right control volume
+     * @param resultFlux: Vector with resulting flux.
+     * @return          : void
+     */
+    void FluxVanLeer( const Matrix& Ax,
+                      const Matrix& AxAbs,
+                      const Matrix& Ay,
+                      const Matrix& AyAbs,
+                      const Matrix& Az,
+                      const Matrix& AzAbs,
+                      const Vector& psiL,
+                      const Vector& psiR,
+                      const Vector& n,
+                      Vector& resultFlux ) const override;
 };
 
 #endif    // UPWINDFLUX_H

code/input/example.cfg

@@ -20,7 +20,7 @@ PROBLEM = LINESOURCE
 % Solver type
 SOLVER = PN_SOLVER
 % CFL number
-CFL_NUMBER = 0.4
+CFL_NUMBER = 0.3
 % Final time for simulation
 TIME_FINAL = 0.3
 
@@ -33,4 +33,4 @@ BC_DIRICHLET = ( void )
 % Quadrature Type
 QUAD_TYPE = MONTE_CARLO
 % Quadrature Order
-QUAD_ORDER = 2 
+QUAD_ORDER = 1

code/src/mesh.cpp

@@ -33,7 +33,7 @@ void Mesh::ComputeConnectivity() {
         sortedBoundaries.push_back( _boundaries[i].second );
         std::sort( sortedBoundaries[i].begin(), sortedBoundaries[i].end() );
     }
-    //#pragma omp parallel for
+#pragma omp parallel for
     for( unsigned i = mpiCellStart; i < mpiCellEnd; ++i ) {
         std::vector<unsigned>* cellsI = &sortedCells[i];
         for( unsigned j = 0; j < _numCells; ++j ) {
@@ -333,3 +333,11 @@ 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<BOUNDARY_TYPE>& Mesh::GetBoundaryTypes() const { return _cellBoundaryTypes; }
+
+double Mesh::GetDistanceToOrigin( unsigned idx_cell ) const {
+    double distance = 0.0;
+    for( unsigned idx_dim = 0; idx_dim < _dim; idx_dim++ ) {
+        distance += _cellMidPoints[idx_cell][idx_dim] * _cellMidPoints[idx_cell][idx_dim];
+    }
+    return sqrt( distance );
+}

code/src/problems/linesource.cpp

@@ -30,7 +30,7 @@ VectorVector LineSource_SN::SetupIC() {
     for( unsigned j = 0; j < cellMids.size(); ++j ) {
         double x = cellMids[j][0];
         double y = cellMids[j][1];
-        psi[j]   = 1 / ( 4 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) ) / ( 4 * M_PI );
+        psi[j]   = 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) ) / ( 4 * M_PI );
     }
     return psi;
 }
@@ -76,7 +76,11 @@ VectorVector LineSource_PN::SetupIC() {
     for( unsigned j = 0; j < cellMids.size(); ++j ) {
         double x  = cellMids[j][0];
         double y  = cellMids[j][1];
-        psi[j][0] = 1 / ( 4 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) ) / ( 4 * M_PI );
+        psi[j][0] = 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) ) / ( 4 * M_PI );
     }
+
+    // for( unsigned j = 0; j < cellMids.size(); ++j ) {
+    //    psi[j][0] = 1;    // / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 4 * t ) ) / ( 4 * M_PI );
+    //}
     return psi;
 }

code/src/settings/config.cpp

@@ -323,7 +323,7 @@ void Config::SetPostprocessing() {
 #endif
 
     // Check if there are contradictive or false options set.
-
+    _boundaries.resize( _nMarkerDirichlet + _nMarkerNeumann );
     // Regroup Boundary Conditions to  std::vector<std::pair<std::string, BOUNDARY_TYPE>> _boundaries;
     for( int i = 0; i < _nMarkerDirichlet; i++ ) {
         _boundaries.push_back( std::pair<std::string, BOUNDARY_TYPE>( _MarkerDirichlet[i], DIRICHLET ) );

code/src/snsolver.cpp

@@ -53,8 +53,10 @@ void SNSolver::Solve() {
         }
         _psi = psiNew;
         for( unsigned i = 0; i < _nCells; ++i ) {
-            fluxNew[i] = dot( _psi[i], _weights );
+            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 );
@@ -71,3 +73,10 @@ void SNSolver::Save() const {
     std::vector<std::vector<std::vector<double>>> results{ scalarField };
     ExportVTK( _settings->GetOutputFile(), results, fieldNames, _mesh );
 }
+
+void SNSolver::Save( int currEnergy ) const {
+    std::vector<std::string> fieldNames{ "flux" };
+    std::vector<std::vector<double>> scalarField( 1, _solverOutput );
+    std::vector<std::vector<std::vector<double>>> results{ scalarField };
+    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( currEnergy ), results, fieldNames, _mesh );
+}

code/src/solver.cpp

@@ -44,6 +44,9 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
 
     // boundary type
     _boundaryCells = _mesh->GetBoundaryTypes();
+
+    // Solver Output
+    _solverOutput.resize( _nCells );    // Currently only Flux
 }
 
 double Solver::ComputeTimeStep( double cfl ) const {
@@ -64,3 +67,16 @@ Solver* Solver::Create( Config* settings ) {
         default: return new SNSolver( settings );
     }
 }
+
+void Solver::Save() const {
+    std::vector<std::string> fieldNames{ "flux" };
+    std::vector<double> flux;
+    flux.resize( _nCells );
+
+    for( unsigned i = 0; i < _nCells; ++i ) {
+        flux[i] = _psi[i][0];
+    }
+    std::vector<std::vector<double>> scalarField( 1, flux );
+    std::vector<std::vector<std::vector<double>>> results{ scalarField };
+    ExportVTK( _settings->GetOutputFile() + "_" + std::to_string( _nEnergies ), results, fieldNames, _mesh );
+}

code/src/upwindflux.cpp

@@ -23,11 +23,40 @@ void UpwindFlux::Flux( const Matrix& AxPlus,
                        const Vector& n,
                        Vector& resultFlux ) const {
     // 2d only atm!!!
-    std::cout << "AxPlus" << AxPlus << std::endl;
-    std::cout << "psiL" << psiL << std::endl;
-    Vector temp = AxPlus * psiL;
-    std::cout << "AxPlus * psiL" << temp << std::endl;
-    std::cout << "n_x *AxPlus * psiL" << n[0] * temp << std::endl;
+    // std::cout << "AxPlus" << AxPlus << std::endl;
+    // std::cout << "psiL" << psiL << std::endl;
+    // Vector temp = AxPlus * psiL;
+    // std::cout << "AxPlus * psiL" << temp << std::endl;
+    // std::cout << "n_x *AxPlus * psiL" << n[0] * temp << std::endl;
 
-    resultFlux = ( n[0] * AxPlus + n[1] * AyPlus ) * psiL + ( n[0] * AxMinus + n[1] * AyMinus ) * psiR;
+    resultFlux += ( n[0] * AyPlus + n[1] * AzPlus ) * psiL + ( n[0] * AyMinus + n[1] * AzMinus ) * psiR;
+}
+
+/**
+ * @brief Flux      : Computes <VanLeer> 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
+ * @param AyPlus    : Positive part of the flux jacobian in y direction
+ * @param AyMinus   : Negative part of the flux jacobian in y direction
+ * @param AzPlus    : Positive part of the flux jacobian in z direction
+ * @param AzMinus   : Negative part of the flux jacobian in z direction
+ * @param psiL      : Solution state of left hand side control volume
+ * @param psiR      : Solution state of right hand side control volume
+ * @param n         : Normal vector at the edge between left and right control volume
+ * @param resultFlux: Vector with resulting flux.
+ * @return          : void
+ */
+void UpwindFlux::FluxVanLeer( const Matrix& Ax,
+                              const Matrix& AxAbs,
+                              const Matrix& Ay,
+                              const Matrix& AyAbs,
+                              const Matrix& Az,
+                              const Matrix& AzAbs,
+                              const Vector& psiL,
+                              const Vector& psiR,
+                              const Vector& n,
+                              Vector& resultFlux ) const {
+    resultFlux += n[0] * 0.5 * Ax * ( psiL + psiR ) - abs( n[0] ) * AxAbs * ( psiR - psiL ) + n[1] * 0.5 * Az * ( psiL + psiR ) -
+                  abs( n[1] ) * AzAbs * ( psiR - psiL );
 }