fixed bug in MN Solver, added multi-cell newton optimization

code/include/optimizers/newtonoptimizer.h

@@ -12,8 +12,8 @@ class NewtonOptimizer : public OptimizerBase
 
     inline ~NewtonOptimizer() {}
 
-    void Solve( Vector& lambda, Vector& u, VectorVector& moments, unsigned idx_cell = 0 ) override;
-    void SolveMultiCell( VectorVector& lambda, VectorVector& u, VectorVector& moments ) override {}
+    void Solve( Vector& lambda, Vector& sol, VectorVector& moments, unsigned idx_cell = 0 ) override;
+    void SolveMultiCell( VectorVector& lambda, VectorVector& sol, VectorVector& moments ) override;
 
   private:
     /*! @brief: Computes gradient of objective function and stores it in grad

code/src/optimizers/newtonoptimizer.cpp

@@ -59,6 +59,23 @@ void NewtonOptimizer::ComputeHessian( Vector& alpha, VectorVector& moments, Matr
     }
 }
 
+void NewtonOptimizer::SolveMultiCell( VectorVector& lambda, VectorVector& sol, VectorVector& moments ) {
+
+    unsigned nCells = lambda.size();
+
+    // if we  have quadratic entropy, then alpha = u;
+    if( _settings->GetEntropyName() == QUADRATIC && _settings->GetNewtonFastMode() ) {
+        for( unsigned idx_cell = 0; idx_cell < nCells; idx_cell++ ) {
+            lambda[idx_cell] = sol[idx_cell];
+        }
+        return;
+    }
+
+    for( unsigned idx_cell = 0; idx_cell < nCells; idx_cell++ ) {
+        Solve( lambda[idx_cell], sol[idx_cell], moments );
+    }
+}
+
 void NewtonOptimizer::Solve( Vector& lambda, Vector& sol, VectorVector& moments, unsigned idx_cell ) {
 
     /* solve the problem argmin ( <eta(alpha*m)>-alpha*u))

code/src/problems/linesource.cpp

@@ -89,11 +89,23 @@ VectorVector LineSource_PN::SetupIC() {
     VectorVector cellMids = _mesh->GetCellMidPoints();
 
     // 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 t      = 3.2e-4;    // pseudo time for gaussian smoothing (Approx to dirac impulse)
+    double offset = 0.13;
     for( unsigned j = 0; j < cellMids.size(); ++j ) {
-        double x  = cellMids[j][0];
-        double y  = cellMids[j][1];
-        psi[j][0] = 1.0 / ( 4.0 * M_PI * t ) * std::exp( -( x * x + y * y ) / ( 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

@@ -156,18 +156,17 @@ void MNSolver::Solve() {
     // auto end   = chrono::steady_clock::now();
 
     // Loop over energies (pseudo-time of continuous slowing down approach)
-
     for( unsigned idx_energy = 0; idx_energy < _nEnergies; idx_energy++ ) {
 
         // Loop over the grid cells
-        for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
+        //        for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
 
-            // ------- Reconstruction Step ----------------
+        // solTimesArea[idx_cell] = _sol[idx_cell] * _areas[idx_cell];    // reconstrucor need moments, not control volume averaged moments!
+        // }
 
-            solTimesArea[idx_cell] = _sol[idx_cell] * _areas[idx_cell];    // reconstrucor need moments, not control volume averaged moments!
-        }
+        // ------- Reconstruction Step ----------------
 
-        _optimizer->SolveMultiCell( _alpha, solTimesArea, _moments );
+        _optimizer->SolveMultiCell( _alpha, _sol, _moments );
 
         // Loop over the grid cells
         for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
@@ -249,11 +248,11 @@ void MNSolver::WriteNNTrainingData( unsigned idx_pseudoTime ) {
         for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
             myfile << "," << _sol[idx_cell][idx_sys];
         }
-        myfile << "\n" << 1 << ", " << _nTotalEntries << "," << idx_pseudoTime;
+        myfile << " \n" << 1 << ", " << _nTotalEntries << "," << idx_pseudoTime;
         for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
             myfile << "," << _alpha[idx_cell][idx_sys];
         }
-        if( idx_cell < _nTotalEntries - 1 ) myfile << "\n";
+        myfile << "\n";
     }
     myfile.close();
 }