Merge branch 'recons' into 'master'

Merge recons into master

See merge request !36

Former-commit-id: 2e59500d

code/include/solvers/mnsolver.h

@@ -42,6 +42,9 @@ class MNSolver : public Solver
                                            Layout: _nCells x _nTotalEntries*/
     OptimizerBase* _optimizer; /*! @brief: Class to solve minimal entropy problem */
 
+    VectorVector _solDx; /*! @brief:  temporary storage of x-derivatives of alpha */
+    VectorVector _solDy; /*! @brief:  temporary storage of y-derivatives of alpha */
+
     // ---- Private Member functions ---
 
     // IO

code/include/solvers/pnsolver.h

@@ -38,8 +38,8 @@ class PNSolver : public Solver
     Vector _scatterMatDiag; /*! @brief: diagonal of the scattering matrix (its a diagonal matrix by construction). Contains eigenvalues of the
                                scattering kernel.  */
 
-    VectorVector _solDx;
-    VectorVector _solDy;
+    VectorVector _solDx; /*! @brief:  temporary storage of x-derivatives of solution */
+    VectorVector _solDy; /*! @brief:  temporary storage of y-derivatives of solution */
 
     // ---- Member functions ----
 

code/include/solvers/solverbase.h

@@ -53,12 +53,12 @@ class Solver
     std::vector<std::vector<unsigned>> _neighbors;
 
     // slope related params
-    Reconstructor* _reconstructor; /*! @brief reconstructor class for high order scheme */
-    unsigned _reconsOrder;
-    VectorVector _psiDx;
-    VectorVector _psiDy;
-    VectorVector _cellMidPoints;
-    std::vector<std::vector<Vector>> _interfaceMidPoints;
+    Reconstructor* _reconstructor; /*! @brief reconstructor object for high-order scheme */
+    unsigned _reconsOrder; /*! @brief reconstruction order (current: 1 & 2) */
+    VectorVector _psiDx; /*! @brief slope of solutions in X direction */
+    VectorVector _psiDy; /*! @brief slope of solutions in Y direction */
+    VectorVector _cellMidPoints; /*! @brief middle point locations of elements */
+    std::vector<std::vector<Vector>> _interfaceMidPoints; /*! @brief middle point locations of edges */
 
     // Solution related members
     VectorVector _sol;                 /*! @brief solution of the PDE, e.g. angular flux or moments */

code/input/exampleMN.cfg

@@ -36,6 +36,8 @@ MAX_MOMENT_SOLVER = 0
 %% Entropy settings
 ENTROPY_FUNCTIONAL = MAXWELL_BOLTZMANN
 ENTROPY_OPTIMIZER = NEWTON
+% Reconstruction order
+RECONS_ORDER = 2
 %
 % ----- Newton Solver Specifications ----
 %

code/src/solvers/mnsolver.cpp

@@ -31,6 +31,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();
@@ -88,29 +92,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++ ) {
-
-        entropyFlux = 0.0;    // Reset temorary flux
+    //--- Temporary storages of reconstructed alpha ---
+    Vector alphaL( _nTotalEntries, 0.0 );
+    Vector alphaR( _nTotalEntries, 0.0 );
 
-        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;
 }
@@ -153,6 +177,10 @@ void MNSolver::ComputeRadFlux() {
 }
 
 void MNSolver::FluxUpdate() {
+    if( _reconsOrder > 1 ) {
+        _mesh->ReconstructSlopesU( _nTotalEntries, _solDx, _solDy, _alpha );    // unstructured reconstruction
+    }
+
     // Loop over the grid cells
     for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
         // Dirichlet Boundaries stay

code/src/solvers/pnsolver.cpp

@@ -32,6 +32,7 @@ PNSolver::PNSolver( Config* settings ) : Solver( settings ) {
     // Initialize Scatter Matrix
     _scatterMatDiag = Vector( _nTotalEntries, 0 );
 
+    // Initialize temporary storages of solution derivatives
     _solDx = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
     _solDy = VectorVector( _nCells, Vector( _nTotalEntries, 0.0 ) );
 
@@ -126,7 +127,8 @@ void PNSolver::FluxUpdate() {
                         solR = _sol[_neighbors[idx_cell][idx_neighbor]] +
                                _solDx[_neighbors[idx_cell][idx_neighbor]] * ( _interfaceMidPoints[idx_cell][idx_neighbor][0] - _cellMidPoints[_neighbors[idx_cell][idx_neighbor]][0] ) + 
                                _solDy[_neighbors[idx_cell][idx_neighbor]] * ( _interfaceMidPoints[idx_cell][idx_neighbor][1] - _cellMidPoints[_neighbors[idx_cell][idx_neighbor]][1] );
-                        // positivity check (if not satisfied, deduce to first order)
+                        // positivity checker (if not satisfied, deduce to first order)
+                        // I manually turned it off here since Pn produces negative solutions essentially
                         //if( min(solL) < 0.0 || min(solR) < 0.0 ) {
                         //    solL = _sol[idx_cell];
                         //    solR = _sol[_neighbors[idx_cell][idx_neighbor]];