MPI started

code/include/snsolver.h

@@ -6,16 +6,17 @@
 
 class SNSolver : public Solver
 {
-public:
+  public:
     /**
      * @brief SNSolver constructor
      * @param settings stores all needed information
      */
-    SNSolver(Settings* settings);
+    SNSolver( Settings* settings );
     /**
      * @brief Solve functions runs main time loop
      */
     virtual void Solve();
+    void SolveMPI();    // can be deleated later
 };
 
-#endif // SNSOLVER_H
+#endif    // SNSOLVER_H

code/include/solver.h

@@ -67,6 +67,11 @@ class Solver
      */
     void ComputeSlopes( VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const;
 
+    /**
+     * @brief SetupIC writes intial condition onto _psi
+     */
+    void SetupIC();
+
   public:
     /**
      * @brief Solver constructor

code/src/snsolver.cpp

 #include "snsolver.h"
+#include <mpi.h>
 
 SNSolver::SNSolver( Settings* settings ) : Solver( settings ) {}
 
 void SNSolver::Solve() {
-    // TODO main SN time loop
+
+    // TODO: set up initial condition
 
     // angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
     VectorVector psiNew = _psi;
@@ -12,18 +14,79 @@ void SNSolver::Solve() {
     for( unsigned n = 0; n < _nTimeSteps; ++n ) {
         // loop over all spatial cells
         for( unsigned j = 0; j < _NCells; ++j ) {
-            Vector scattering = _sigmaSH20[j] * _psi[j];
 
             // loop over all ordinates
             for( unsigned k = 0; k < _nq; ++k ) {
                 psiNew[j][k] = 0.0;
                 // loop over all neighbor cells (edges) of cell j and compute numerical fluxes
                 for( unsigned l = 0; l < _neighbors[j].size(); ++l ) {
+                    // store flux contribution on psiNew to save memory
+                    psiNew[j][k] -= ( _dt / _areas[j] ) * _g->Flux( _quadPoints[k], _psi[j][k], _psi[_neighbors[j][l]][k], _normals[j][l] );
+                }
+                // time update angular flux with numerical flux and total scattering cross section
+                psiNew[j][k] = _psi[j][k] + psiNew[j][k] + _dt * _sigmaTH20[n] * _psi[j][k];
+            }
+            // compute scattering effects
+            psiNew[j] += _sigmaSH20[n] * _psi[j] * _weights;    // multiply scattering matrix with psi
+        }
+        _psi = psiNew;
+        // psiNew.reset();
+    }
+}
+
+void SNSolver::SolveMPI() {
+
+    // setup MPI variables: mype is PE index, npes is number of PEs, ierr is MPI message flag
+    int mype, npes, ierr;
+
+    // determine mype and npes
+    ierr = MPI_Comm_rank( MPI_COMM_WORLD, &mype );
+    ierr = MPI_Comm_size( MPI_COMM_WORLD, &npes );
+
+    // determine size of quadrature array
+    int nqPE = int( ( _nq - 1 ) / npes ) + 1;
+    // nqPEMax needed for allocation
+    int nqPEMax = nqPE;
+    if( mype == npes - 1 ) {
+        nqPE = _nq - mype * nqPE;
+        if( nqPE < 0 ) {
+            nqPE = 0;
+        }
+    }
+
+    int master = 0;    // define PE 0 to be master
+    int tag    = 0;
+    int kStart = mype * ( ( _nq - 1 ) / npes + 1.0 );
+    int kEnd   = kStart + nqPE - 1;
+
+    // TODO: store local psi: dim(psi) = (NCells,nqPE)
+    // TODO: store local scattering matrix: dim(_sigmaSH20) = (nTimeSteps,_nq,nqPE)
+    // TODO: store local weights: dim(_weights) = (nqPE)
+
+    // TODO: set up initial condition
+
+    // angular flux at next time step (maybe store angular flux at all time steps, since time becomes energy?)
+    VectorVector psiNew = _psi;
+
+    // loop over energies (pseudo-time)
+    for( unsigned n = 0; n < _nTimeSteps; ++n ) {
+        // loop over all spatial cells
+        for( unsigned j = 0; j < _NCells; ++j ) {
+
+            // loop over all ordinates
+            for( unsigned k = 0; k < unsigned( nqPE ); ++k ) {
+                psiNew[j][k] = 0.0;
+                // loop over all neighbor cells (edges) of cell j and compute numerical fluxes
+                for( unsigned l = 0; l < _neighbors[j].size(); ++l ) {
+                    // store flux contribution on psiNew to save memory
                     psiNew[j][k] -= ( _dt / _areas[j] ) * _g->Flux( _quadPoints[k], _psi[j][k], _psi[_neighbors[j][l]][k], _normals[j][l] );
                 }
                 // time update angular flux with numerical flux and total scattering cross section
                 psiNew[j][k] = _psi[j][k] + psiNew[j][k] + _dt * _sigmaTH20[n] * _psi[j][k];
             }
+            // compute scattering effects
+            Vector scattering = _sigmaSH20[n] * _psi[j] * _weights;    // multiply scattering matrix with dim(scattering) = _nq -> improve later
+            psiNew[j] += scattering;
         }
         _psi = psiNew;
         // psiNew.reset();

code/src/solver.cpp

@@ -28,6 +28,11 @@ Solver::Solver( Settings* settings ) : _NCells( 1 ) {
     // TODO: load density and stopping power data
     LoadPatientDensity( "test" );
     LoadStoppingPower( "test" );
+    LoadSigmaS( "test" );
+    LoadSigmaT( "test" );
+
+    // TODO: setup initial condtion (distribution at initial energy for CSD)
+    SetupIC();
 }
 
 double Solver::ComputeTimeStep( double cfl ) const {
@@ -74,4 +79,10 @@ void Solver::LoadSigmaT( std::string fileName ) {
     //_sigmaTH20 = ... -> dim(_sigmaTH20) = _nTimeSteps
 }
 
+void Solver::SetupIC() {
+    // @TODO
+    // write initial condition
+    // _psi = ...
+}
+
 Solver* Solver::Create( Settings* settings ) { return new SNSolver( settings ); }