added some basic functions

Former-commit-id: 29ee0750

code/include/solvers/csdpnsolver.h

@@ -34,4 +34,4 @@ class CSDPNSolver : public PNSolver
     void WriteVolumeOutput( unsigned idx_pseudoTime ) override;
 };
 
-#endif // CSDPNSOLVER_H
+#endif    // CSDPNSOLVER_H

code/input/DataGenerator1D.cfg

@@ -7,12 +7,12 @@
 %
 % ---- Datagenerator settings ----
 DATA_GENERATOR_MODE = YES
-TRAINING_SET_SIZE = 100000
+TRAINING_SET_SIZE = 20
 MAX_VALUE_FIRST_MOMENT = 1
 SPATIAL_DIM = 1
 REALIZABLE_SET_EPSILON_U1 = 0.003
 REALIZABLE_SET_EPSILON_U0 = 0.003
-NORMALIZED_SAMPLING = YES
+NORMALIZED_SAMPLING = NO
 MAX_MOMENT_SOLVER = 1
 %
 % ---- File specifications ----

code/src/solvers/csdpnsolver.cpp

-#include "solvers/CSDPNSolver.h"
+#include "solvers/csdpnsolver.h"
 #include "common/config.h"
 #include "common/io.h"
 #include "fluxes/numericalflux.h"
@@ -10,12 +10,110 @@
 #include "spdlog/spdlog.h"
 #include <mpi.h>
 
-CSDPNSolver::CSDPNSolver( Config* settings ) : PNSolver( settings ) {
+CSDPNSolver::CSDPNSolver( Config* settings ) : PNSolver( settings ) {}
+
+void CSDPNSolver::IterPreprocessing( unsigned /*idx_iter*/ ) {
+    // Nothing to preprocess for PNSolver
+}
+
+void CSDPNSolver::IterPostprocessing( unsigned /*idx_iter*/ ) {
+    // --- Update Solution ---
+    _sol = _solNew;
+
+    // --- Compute Flux for solution and Screen Output ---
+    ComputeRadFlux();
+}
+
+void CSDPNSolver::ComputeRadFlux() {
+    double firstMomentScaleFactor = sqrt( 4 * M_PI );
+    for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
+        _fluxNew[idx_cell] = _sol[idx_cell][0] * firstMomentScaleFactor;
+    }
 }
 
-void CSDPNSolver::Solve() {
+void CSDPNSolver::FluxUpdate() {
+    if( _reconsOrder > 1 ) {
+        _mesh->ReconstructSlopesU( _nSystem, _solDx, _solDy, _sol );    // unstructured reconstruction
+        //_mesh->ComputeSlopes( _nTotalEntries, _solDx, _solDy, _sol );    // unstructured reconstruction
+    }
+    // Vector solL( _nTotalEntries );
+    // Vector solR( _nTotalEntries );
+    auto solL = _sol[2];
+    auto solR = _sol[2];
+
+    // Loop over all spatial cells
+    for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
+
+        // Dirichlet cells stay at IC, farfield assumption
+        if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
+
+        // Reset temporary variable psiNew
+        for( unsigned idx_sys = 0; idx_sys < _nSystem; idx_sys++ ) {
+            _solNew[idx_cell][idx_sys] = 0.0;
+        }
+
+        // Loop over all neighbor cells (edges) of cell j and compute numerical fluxes
+        for( unsigned idx_neighbor = 0; idx_neighbor < _neighbors[idx_cell].size(); idx_neighbor++ ) {
+
+            // Compute flux contribution and store in psiNew to save memory
+            if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::NEUMANN && _neighbors[idx_cell][idx_neighbor] == _nCells )
+                _solNew[idx_cell] += _g->Flux(
+                    _AxPlus, _AxMinus, _AyPlus, _AyMinus, _AzPlus, _AzMinus, _sol[idx_cell], _sol[idx_cell], _normals[idx_cell][idx_neighbor] );
+            else {
+                switch( _reconsOrder ) {
+                    // first order solver
+                    case 1:
+                        _solNew[idx_cell] += _g->Flux( _AxPlus,
+                                                       _AxMinus,
+                                                       _AyPlus,
+                                                       _AyMinus,
+                                                       _AzPlus,
+                                                       _AzMinus,
+                                                       _sol[idx_cell],
+                                                       _sol[_neighbors[idx_cell][idx_neighbor]],
+                                                       _normals[idx_cell][idx_neighbor] );
+                        break;
+                    // second order solver
+                    case 2:
+                        // left status of interface
+                        solL = _sol[idx_cell] + _solDx[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_neighbor][0] - _cellMidPoints[idx_cell][0] ) +
+                               _solDy[idx_cell] * ( _interfaceMidPoints[idx_cell][idx_neighbor][1] - _cellMidPoints[idx_cell][1] );
+                        // right status of interface
+                        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 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]];
+                        //}
+
+                        // flux evaluation
+                        _solNew[idx_cell] +=
+                            _g->Flux( _AxPlus, _AxMinus, _AyPlus, _AyMinus, _AzPlus, _AzMinus, solL, solR, _normals[idx_cell][idx_neighbor] );
+                        break;
+                    // default: first order solver
+                    default:
+                        _solNew[idx_cell] += _g->Flux( _AxPlus,
+                                                       _AxMinus,
+                                                       _AyPlus,
+                                                       _AyMinus,
+                                                       _AzPlus,
+                                                       _AzMinus,
+                                                       _sol[idx_cell],
+                                                       _sol[_neighbors[idx_cell][idx_neighbor]],
+                                                       _normals[idx_cell][idx_neighbor] );
+                }
+            }
+        }
+    }
 }
 
+void CSDPNSolver::FVMUpdate( unsigned idx_energy ) {}
+
 void CSDPNSolver::PrepareVolumeOutput() {
     unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
 
@@ -56,7 +154,7 @@ void CSDPNSolver::WriteVolumeOutput( unsigned idx_pseudoTime ) {
     // Compute total "mass" of the system ==> to check conservation properties
     std::vector<double> flux( _nCells, 0.0 );
     for( unsigned idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
-        flux[idx_cell] = dot( _sol[idx_cell], _weights );
+        // flux[idx_cell] = dot( _sol[idx_cell], _weights );
         mass += flux[idx_cell] * _areas[idx_cell];
     }
 

code/src/solvers/pnsolver.cpp

@@ -72,7 +72,7 @@ void PNSolver::ComputeRadFlux() {
     }
 }
 
-void PNSolver::FluxUpdate() {
+void PNSolver::FluxUpdate() override {
     if( _reconsOrder > 1 ) {
         _mesh->ReconstructSlopesU( _nSystem, _solDx, _solDy, _sol );    // unstructured reconstruction
         //_mesh->ComputeSlopes( _nTotalEntries, _solDx, _solDy, _sol );    // unstructured reconstruction
@@ -153,7 +153,7 @@ void PNSolver::FluxUpdate() {
     }
 }
 
-void PNSolver::FVMUpdate( unsigned idx_energy ) {
+void PNSolver::FVMUpdate( unsigned idx_energy ) override {
     // Loop over all spatial cells
     for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
         // Dirichlet cells stay at IC, farfield assumption
@@ -370,7 +370,7 @@ double PNSolver::LegendrePoly( double x, int l ) {    // Legacy. TO BE DELETED
     }
 }
 
-void PNSolver::PrepareVolumeOutput() {
+void PNSolver::PrepareVolumeOutput() override {
     unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
 
     _outputFieldNames.resize( nGroups );
@@ -410,7 +410,7 @@ void PNSolver::PrepareVolumeOutput() {
     }
 }
 
-void PNSolver::WriteVolumeOutput( unsigned idx_pseudoTime ) {
+void PNSolver::WriteVolumeOutput( unsigned idx_pseudoTime ) override {
     unsigned nGroups = (unsigned)_settings->GetNVolumeOutput();
 
     if( ( _settings->GetVolumeOutputFrequency() != 0 && idx_pseudoTime % (unsigned)_settings->GetVolumeOutputFrequency() == 0 ) ||