solve #60 (realizability)

code/include/optimizers/newtonoptimizer.h

@@ -16,6 +16,10 @@ class NewtonOptimizer : public OptimizerBase
     void SolveMultiCell( VectorVector& lambda, VectorVector& sol, VectorVector& moments ) override;
 
   private:
+    /*! @brief: Computes the objective function
+                grad = <eta(alpha*m)> - alpha*sol */
+    double ComputeObjFunc( Vector& alpha, Vector& sol, VectorVector& moments );
+
     /*! @brief: Computes gradient of objective function and stores it in grad
                 grad = <m*eta*'(alpha*m)> - sol */
     void ComputeGradient( Vector& alpha, Vector& sol, VectorVector& moments, Vector& grad );
@@ -24,8 +28,6 @@ class NewtonOptimizer : public OptimizerBase
                 grad = <mXm*eta*'(alpha*m)> */
     void ComputeHessian( Vector& alpha, VectorVector& moments, Matrix& hessian );
 
-    double ComputeObjFunc( Vector& alpha, Vector& sol, VectorVector& moments );
-
     QuadratureBase* _quadrature; /*! @brief: used quadrature */    // THis is memory doubling! Try to use a pointer.
     unsigned _nq;                                                  /*! @brief: number of quadrature points */
     Vector _weights;                                               /*!  @brief quadrature weights, dim(_weights) = (_nq) */

code/src/problems/linesource.cpp

@@ -90,22 +90,21 @@ VectorVector LineSource_PN::SetupIC() {
 
     // Initial condition is dirac impulse at (x,y) = (0,0) ==> constant in angle ==> all moments are zero.
     double t      = 3.2e-4;    // pseudo time for gaussian smoothing (Approx to dirac impulse)
-    double offset = 0.13;
+    double offset = 0;
     for( unsigned j = 0; j < cellMids.size(); ++j ) {
-        double x = cellMids[j][0];
-        double y = cellMids[j][1];    // (x- 0.5) * (x- 0.5)
-        psi[j][0] =
-            1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x - offset ) * ( x - offset ) + ( y - offset ) * ( y - offset ) ) / ( 4 * t ) ) +
-            1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x + offset ) * ( x + offset ) + ( y - offset ) * ( y - offset ) ) / ( 4 * t ) ) +
-            1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x - offset ) * ( x - offset ) + ( y + offset ) * ( y + offset ) ) / ( 4 * t ) ) +
-            1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x + offset ) * ( x + offset ) + ( y + offset ) * ( y + offset ) ) / ( 4 * t ) ) +
-            1.0 / ( 4.0 * M_PI * t ) *
-                std::exp( -( ( x - 2 * offset ) * ( x - 2 * offset ) + ( y - 2 * offset ) * ( y - 2 * offset ) ) / ( 4 * t ) ) +
-            1.0 / ( 4.0 * M_PI * t ) *
-                std::exp( -( ( x + 2 * offset ) * ( x + 2 * offset ) + ( y - 2 * offset ) * ( y - 2 * offset ) ) / ( 4 * t ) ) +
-            1.0 / ( 4.0 * M_PI * t ) *
-                std::exp( -( ( x - 2 * offset ) * ( x - 2 * offset ) + ( y + 2 * offset ) * ( y + 2 * offset ) ) / ( 4 * t ) ) +
-            1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x + 2 * offset ) * ( x + 2 * offset ) + ( y + 2 * offset ) * ( y + 2 * offset ) ) / ( 4 * t ) );
+        double x  = cellMids[j][0];
+        double y  = cellMids[j][1];    // (x- 0.5) * (x- 0.5)
+        psi[j][0] = 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x - offset ) * ( x - offset ) + ( y - offset ) * ( y - offset ) ) / ( 4 * t ) );    //+
+        // 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x + offset ) * ( x + offset ) + ( y - offset ) * ( y - offset ) ) / ( 4 * t ) ) +
+        // 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x - offset ) * ( x - offset ) + ( y + offset ) * ( y + offset ) ) / ( 4 * t ) ) +
+        // 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x + offset ) * ( x + offset ) + ( y + offset ) * ( y + offset ) ) / ( 4 * t ) ) +
+        // 1.0 / ( 4.0 * M_PI * t ) *
+        //    std::exp( -( ( x - 2 * offset ) * ( x - 2 * offset ) + ( y - 2 * offset ) * ( y - 2 * offset ) ) / ( 4 * t ) ) +
+        // 1.0 / ( 4.0 * M_PI * t ) *
+        //    std::exp( -( ( x + 2 * offset ) * ( x + 2 * offset ) + ( y - 2 * offset ) * ( y - 2 * offset ) ) / ( 4 * t ) ) +
+        // 1.0 / ( 4.0 * M_PI * t ) *
+        //    std::exp( -( ( x - 2 * offset ) * ( x - 2 * offset ) + ( y + 2 * offset ) * ( y + 2 * offset ) ) / ( 4 * t ) ) +
+        // 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( ( x + 2 * offset ) * ( x + 2 * offset ) + ( y + 2 * offset ) * ( y + 2 * offset ) ) / ( 4 * t ) );
     }
 
     // Debugging jump test case

code/src/solvers/mnsolver.cpp

@@ -59,6 +59,31 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ) {
 
     // Solver output
     _outputFields = std::vector( _nTotalEntries, std::vector( _nCells, 0.0 ) );
+
+    // Delete that
+    /*
+        std::string filename = "obFunc.csv";
+        std::ofstream myfile;
+        myfile.open( filename );
+
+        NewtonOptimizer* optimizer = new NewtonOptimizer( settings );
+
+        Vector alpha( 4, 0.0 );
+        Vector u( 4, 0.0 );
+        u[0] = -1.11;
+        u[1] = 0;
+        u[2] = 0;
+        u[4] = 0;
+
+        for( double a1 = -1000; a1 <= 100; a1 += 10.0 / 100.0 ) {
+            for( double a2 = -60; a2 <= 100; a2 += 60.0 / 40.0 ) {
+                alpha[0] = a1;
+                alpha[1] = a2;
+                myfile << a1 << ", " << a2 << ", " << optimizer->ComputeObjFunc( alpha, u, _moments ) << "\n";
+            }
+        }
+        myfile.close();
+     */
 }
 
 MNSolver::~MNSolver() {
@@ -173,7 +198,7 @@ void MNSolver::Solve() {
 
             // ------- Relizablity Reconstruction Step ----
 
-            // ComputeRealizableSolution( idx_cell );
+            ComputeRealizableSolution( idx_cell );
 
             // ------- Flux Computation Step --------------