aircavity tested

code/src/problems/aircavity1d.cpp

@@ -21,10 +21,16 @@ VectorVector AirCavity1D::SetupIC() {
     return psi;
 }
 
-std::vector<double> AirCavity1D::GetDensity( const VectorVector& cellMidPoints ) { 
-    std::vector<double> densities ( 4*cellMidPoints.size()/9, 1.0 );
-    std::vector<double> air( 2*cellMidPoints.size()/9, 0.00125  );
-    std::vector<double> water ( 3*cellMidPoints.size()/9, 1.0 );
-    densities.insert(densities.end(),air.begin(),air.end());
-    densities.insert(densities.end(),water.begin(),water.end());
-    return densities; }
+std::vector<double> AirCavity1D::GetDensity( const VectorVector& cellMidPoints ) {
+    std::vector<double> densities( cellMidPoints.size(), 1.0 );
+    /*std::vector<double> densities( 4 * cellMidPoints.size() / 9, 1.0 );
+    std::vector<double> air( 2 * cellMidPoints.size() / 9, 0.00125 );
+    std::vector<double> water( 3 * cellMidPoints.size() / 9, 1.0 );
+    densities.insert( densities.end(), air.begin(), air.end() );
+    densities.insert( densities.end(), water.begin(), water.end() );
+    return densities;*/
+    for( unsigned j = 0; j < cellMidPoints.size(); ++j ) {
+        if( cellMidPoints[j][0] > 1.5 - 2.5 && cellMidPoints[j][0] < 2.0 - 2.5 ) densities[j] = 0.01;
+    }
+    return densities;
+}

code/src/solvers/csdsolvertrafofp.cpp

@@ -102,8 +102,10 @@ void CSDSolverTrafoFP::Solve() {
     auto log = spdlog::get( "event" );
 
     auto energiesOrig = _energies;
+    //_density          = std::vector<double>( _density.size(), 1.0 );
 
     // setup incoming BC on left
+    _sol = VectorVector( _density.size(), Vector( _settings->GetNQuadPoints(), 0.0 ) );
     for( unsigned k = 0; k < _nq; ++k ) {
         if( _quadPoints[k][0] > 0 && !_RT ) _sol[0][k] = 1e5 * exp( -10.0 * pow( 1.0 - _quadPoints[k][0], 2 ) );
     }
@@ -133,8 +135,6 @@ void CSDSolverTrafoFP::Solve() {
         }
     }
 
-    VectorVector psi1 = _sol;
-
     // store transformed energies ETilde instead of E in _energies vector (cf. Dissertation Kerstion Kuepper, Eq. 1.25)
     double tmp   = 0.0;
     _energies[0] = 0.0;
@@ -188,11 +188,8 @@ void CSDSolverTrafoFP::Solve() {
         if( _RT ) {
             for( unsigned k = 0; k < _nq; ++k ) {
                 if( _quadPoints[k][0] > 0 ) {
-                    //_sol[0][k] = 1e5 * exp( -200.0 * pow( 1.0 - _quadPoints[k][0], 2 ) ) * exp(-50*pow(_energies[0]-_energies[n],2))* _density[0] *
-                    //_s[_nEnergies - n- 1];
                     _sol[0][k] = 1e5 * exp( -200.0 * pow( 1.0 - _quadPoints[k][0], 2 ) ) *
                                  exp( -50.0 * pow( _energyMax - energiesOrig[_nEnergies - n - 1], 2 ) ) * _density[0] * _s[_nEnergies - n - 1];
-                    // psiNew[0][k] = _sol[0][k];
                 }
             }
         }
@@ -223,7 +220,7 @@ void CSDSolverTrafoFP::Solve() {
                                         _areas[j];
                 }
                 // time update angular flux with numerical flux and total scattering cross section
-                psiNew[j][i] = _sol[j][i] - _dE * psiNew[j][i];    // + _dE * psi1[j][i];
+                psiNew[j][i] = _sol[j][i] - _dE * psiNew[j][i];
             }
         }
 
@@ -281,14 +278,14 @@ void CSDSolverTrafoFP::GenerateEnergyGrid( bool refinement ) {
         }
     }
     else {
-        double energySwitch    = 0.11;
+        double energySwitch    = 0.58;
         double energySwitchMin = 0.03;
         // write equidistant energy grid
 
         _dE                 = ComputeTimeStep( _settings->GetCFL() );
         unsigned nEnergies1 = unsigned( ( _energyMax - energySwitch ) / _dE );
-        unsigned nEnergies2 = unsigned( ( energySwitch - energySwitchMin ) / ( _dE ) );
-        unsigned nEnergies3 = unsigned( ( energySwitchMin - _energyMin ) / ( _dE / 1 ) );
+        unsigned nEnergies2 = unsigned( ( energySwitch - energySwitchMin ) / ( _dE / 2 ) );
+        unsigned nEnergies3 = unsigned( ( energySwitchMin - _energyMin ) / ( _dE / 2 ) );
         _nEnergies          = nEnergies1 + nEnergies2 + nEnergies3 - 2;
         std::cout << "nEnergies1 = " << nEnergies1 << std::endl;
         std::cout << "nEnergies2 = " << nEnergies2 << std::endl;

code/src/solvers/solverbase.cpp

@@ -56,7 +56,6 @@ Solver::Solver( Config* settings ) : _settings( settings ) {
     // store density
     VectorVector cellMids = _mesh->GetCellMidPoints();
     _density              = _problem->GetDensity( cellMids );
-    std::cout << "Density set with size " << _density.size() << std::endl;
 }
 
 Solver::~Solver() {