tackle realizability issue

code/include/solvers/mnsolver.h

@@ -61,5 +61,9 @@ class MNSolver : public Solver
 
     /*! @brief : Pre-Compute Moments at all quadrature points. */
     void ComputeMoments();
+
+    /*! @brief Corrects the solution _sol[idx_cell] to be realizable w.r.t. the reconstructed entropy (eta'(alpha*m))
+        @param idx_cell = cell where the correction happens*/
+    void ComputeRealizableSolution( unsigned idx_cell );
 };
 #endif    // MNSOLVER_H

code/src/optimizers/newtonoptimizer.cpp

@@ -68,10 +68,10 @@ void NewtonOptimizer::Solve( Vector& lambda, Vector& sol, VectorVector& moments
      */
 
     // if we  have quadratic entropy, then alpha = u;
-    if( _settings->GetEntropyName() == QUADRATIC && _settings->GetNewtonFastMode() ) {
-        lambda = sol;
-        return;
-    }
+    // if( _settings->GetEntropyName() == QUADRATIC && _settings->GetNewtonFastMode() ) {
+    //    lambda = sol;
+    //    return;
+    //}
 
     // Start Newton Algorithm
 

code/src/solvers/mnsolver.cpp

@@ -105,6 +105,14 @@ Vector MNSolver::ConstructFlux( unsigned idx_cell ) {
     return flux;
 }
 
+void MNSolver::ComputeRealizableSolution( unsigned idx_cell ) {
+    double entropyReconstruction = 0.0;
+
+    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
+        entropyReconstruction = _entropy->EntropyPrimeDual( blaze::dot( _alpha[idx_cell], _moments[idx_quad] ) );
+    }
+}
+
 void MNSolver::Solve() {
 
     int rank;
@@ -143,19 +151,22 @@ void MNSolver::Solve() {
         // Loop over the grid cells
         for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
 
-            // ------- Reconstruction Step -------
+            // ------- Reconstruction Step ----------------
 
-            // _alpha[idx_cell] = _sol[idx_cell];
             _optimizer->Solve( _alpha[idx_cell], _sol[idx_cell], _moments );
 
-            // ------- Flux Computation Step ---------
+            // ------- Relizablity Reconstruction Step ----
+
+            ComputeRealizableSolution( idx_cell );
+
+            // ------- Flux Computation Step --------------
 
             // Dirichlet Boundaries are finished now
             if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
 
             psiNew[idx_cell] = ConstructFlux( idx_cell );
 
-            // ------ FVM Step ------
+            // ------ Finite Volume Update Step ------
 
             // NEED TO VECTORIZE
             for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {