update to flux comp and cross sections

code/src/fluxes/upwindflux.cpp

@@ -3,7 +3,7 @@
 UpwindFlux::UpwindFlux() : NumericalFlux() {}
 
 double UpwindFlux::Flux( const Vector& Omega, double psiL, double psiR, const Vector& n ) const {
-    double inner = Omega[0] * n[0] + Omega[1] * n[1]; // Only use x and y axis in 2d case
+    double inner = Omega[0] * n[0] + Omega[1] * n[1];    // Only use x and y axis in 2d case
     if( inner > 0 ) {
         return inner * psiL;
     }

code/src/solvers/csdpnsolver.cpp

@@ -41,6 +41,7 @@ CSDPNSolver::CSDPNSolver( Config* settings ) : PNSolver( settings ) {
     //_sigmaTE  = Vector( _nEnergies, 0.0 );
     //
 
+    // COmpute INitial Conditions
     Vector pos_beam = Vector{ 0.5, 0.5 };
     VectorVector IC( _nCells, Vector( _nSystem ) );
     for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
@@ -52,6 +53,9 @@ CSDPNSolver::CSDPNSolver( Config* settings ) : PNSolver( settings ) {
             IC[idx_cell][idx_sys] = f * StarMAPmoments[idx_sys];    // must be VectorVector
         }
     }
+    // Get Initial Condition, initialize _sol ( Later move to problem)
+    _sol    = IC;
+    _solNew = _sol;
 
     // printf( "%d", sigma_ref.rows() );
 
@@ -59,21 +63,28 @@ CSDPNSolver::CSDPNSolver( Config* settings ) : PNSolver( settings ) {
     // std::cout << sigma_ref.rows() << std::endl;
     // std::cout << _energies.size() << std::endl;
 
+    // Implement Cross Sections
     _sigmaT = VectorVector( _polyDegreeBasis, Vector( _energies.size() ) );
     for( unsigned idx_degree = 0; idx_degree < _polyDegreeBasis; ++idx_degree ) {
         Vector xs_m = blaze::column( sigma_ref, idx_degree );    // Scattering cross section Moments
         Interpolation interp( E_ref, xs_m );
         _sigmaT[idx_degree] = interp( _energies );
     }
+    _sigmaS = _sigmaT;    // Careful! Hardcoded dims: polydegree x nEnergies.
 
-    // std::cout << "here\n";
-    // std::cout << size( sigma_t ) << std::endl;
-    _sigmaT = sigma_t;
-
-    TextProcessingToolbox::PrintMatrix( sigma_t );
-
+    // Implement ??? (Stopping Powers?)
     Interpolation interpS( E_tab, S_tab );
     _s = interpS( _energies );
+
+    // Sanity Check
+    // TextProcessingToolbox::PrintMatrix( _Ax );
+    // TextProcessingToolbox::PrintMatrix( _Ay );
+    // TextProcessingToolbox::PrintMatrix( _Az );
+    //
+    // TextProcessingToolbox::PrintMatrix( _AxPlus );
+    // TextProcessingToolbox::PrintMatrix( _AyPlus );
+    // TextProcessingToolbox::PrintMatrix( _AxMinus );
+    // TextProcessingToolbox::PrintMatrix( _AyMinus );
 }
 
 void CSDPNSolver::SolverPreprocessing() {
@@ -152,11 +163,17 @@ void CSDPNSolver::FVMUpdate( unsigned idx_energy ) {
         // Dirichlet cells stay at IC, farfield assumption
         if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
         // Flux update
-        for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
-            _solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys] - ( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys] /* cell averaged flux */
-                                         - _dE * _sol[idx_cell][idx_sys] *
-                                               ( _sigmaT[idx_energy][idx_cell]                                 /* absorbtion influence */
-                                                 - _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_sys] ); /* scattering influence */
+        // for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
+        for( int idx_l = 0; idx_l <= (int)_polyDegreeBasis; idx_l++ ) {
+            for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
+                int idx_sys = GlobalIndex( idx_l, idx_k );
+
+                _solNew[idx_cell][idx_sys] = _sol[idx_cell][idx_sys] -
+                                             ( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_sys] /* cell averaged flux */
+                                             - _dE * _sol[idx_cell][idx_sys] *
+                                                   ( _sigmaT[idx_energy][idx_l]      /* absorbtion influence */
+                                                     - _sigmaS[idx_energy][idx_l] ); /* scattering influence */
+            }
         }
         // Source Term
         _solNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];

code/src/solvers/pnsolver.cpp

@@ -339,7 +339,7 @@ void PNSolver::ComputeFluxComponents() {
     // std::cout << "Spectral Radius Y " << blaze::max( blaze::abs( eigenValuesY ) ) << "\n";
     // std::cout << "Spectral Radius Z " << blaze::max( blaze::abs( eigenValues ) ) << "\n";
 
-    // _combinedSpectralRadius = blaze::max( blaze::abs( eigenValues + eigenValuesX + eigenValuesY ) );
+    //_combinedSpectralRadius = blaze::max( blaze::abs( eigenValues + eigenValuesX + eigenValuesY ) );
     // std::cout << "Spectral Radius combined " << _combinedSpectralRadius << "\n";
 }