fix couplong constants

src/PartBunch/PartBunch.cpp

@@ -2,6 +2,62 @@
 #include <boost/numeric/ublas/io.hpp>
 #include "Utilities/Util.h"
 
+template<>
+void PartBunch<double,3>::spaceChargeEFieldCheck() {
+
+ Inform msg("EParticleStats");
+
+ auto pE_view = this->pcontainer_m->E.getView();
+ 
+ double avgE          = 0.0;
+ double minEComponent = std::numeric_limits<double>::max();
+ double maxEComponent = std::numeric_limits<double>::min();
+ double minE          = std::numeric_limits<double>::max();
+ double maxE          = std::numeric_limits<double>::min();
+ double cc            = getCouplingConstant();
+ 
+ int myRank = ippl::Comm->rank();
+ Kokkos::parallel_reduce(
+                         "check e-field", this->getLocalNum(),
+                         KOKKOS_LAMBDA(const int i, double& loc_avgE, double& loc_minEComponent,
+                                       double& loc_maxEComponent, double& loc_minE, double& loc_maxE) {
+                             double EX    = pE_view[i][0]*cc;
+                             double EY    = pE_view[i][1]*cc;
+                             double EZ    = pE_view[i][2]*cc;
+
+                             double ENorm = Kokkos::sqrt(EX*EX + EY*EY + EZ*EZ);
+                             
+                             loc_avgE += ENorm;
+
+                             loc_minEComponent = EX < loc_minEComponent ? EX : loc_minEComponent;
+                             loc_minEComponent = EY < loc_minEComponent ? EY : loc_minEComponent;
+                             loc_minEComponent = EZ < loc_minEComponent ? EZ : loc_minEComponent;
+                             
+                             loc_maxEComponent = EX > loc_maxEComponent ? EX : loc_maxEComponent;
+                             loc_maxEComponent = EY > loc_maxEComponent ? EY : loc_maxEComponent;
+                             loc_maxEComponent = EZ > loc_maxEComponent ? EZ : loc_maxEComponent;
+
+                             loc_minE = ENorm < loc_minE ? ENorm : loc_minE;
+                             loc_maxE = ENorm > loc_maxE ? ENorm : loc_maxE;
+                         },
+                         Kokkos::Sum<double>(avgE), Kokkos::Min<double>(minEComponent),
+                         Kokkos::Max<double>(maxEComponent), Kokkos::Min<double>(minE),
+                         Kokkos::Max<double>(maxE));
+
+ MPI_Reduce(myRank == 0 ? MPI_IN_PLACE : &avgE, &avgE, 1, MPI_DOUBLE, MPI_SUM, 0, ippl::Comm->getCommunicator());
+ MPI_Reduce(myRank == 0 ? MPI_IN_PLACE : &minEComponent, &minEComponent, 1, MPI_DOUBLE, MPI_MIN, 0, ippl::Comm->getCommunicator());
+ MPI_Reduce(myRank == 0 ? MPI_IN_PLACE : &maxEComponent, &maxEComponent, 1, MPI_DOUBLE, MPI_MAX, 0, ippl::Comm->getCommunicator());
+ MPI_Reduce(myRank == 0 ? MPI_IN_PLACE : &minE, &minE, 1, MPI_DOUBLE, MPI_MIN, 0, ippl::Comm->getCommunicator());
+ MPI_Reduce(myRank == 0 ? MPI_IN_PLACE : &maxE, &maxE, 1, MPI_DOUBLE, MPI_MAX, 0, ippl::Comm->getCommunicator());
+ 
+ avgE /= this->getTotalNum();
+ msg << "avgENorm = " << avgE << endl;
+ msg << "minEComponent = " << minEComponent << endl;
+ msg << "maxEComponent = " << maxEComponent << endl;
+ msg << "minE = " << minE << endl;
+ msg << "maxE = " << maxE << endl;
+}
+
 template<>
 void PartBunch<double,3>::calcBeamParameters() {
         // Usefull constants
@@ -97,10 +153,7 @@ void PartBunch<double,3>::calcBeamParameters() {
             rmax_m(i) = rmax[i];
             rmin_m(i) = rmin[i];
         }
-
         ippl::Comm->barrier();
-
-
     }
 
 template<>
@@ -285,9 +338,7 @@ void PartBunch<double,3>::computeSelfFields() {
     //divide charge by a 'grid-cube' volume to get [C/m^3]
     rho = rho *tmp2;
         
-    double Npoints = nr_m[0] * nr_m[1] * nr_m[2];
 
-    // auto rhoView = rho.getView();
     localDensity2 = 0.;
     ippl::parallel_reduce(
                           "Density stats", ippl::getRangePolicy(rhoView),
@@ -297,18 +348,22 @@ void PartBunch<double,3>::computeSelfFields() {
                           },
                           Kokkos::Sum<double>(localDensity2) );
     ippl::Comm->reduce(localDensity2, Density2, 1, std::plus<double>());
-    
+
+    double Npoints = nr_m[0] * nr_m[1] * nr_m[2];
     rmsDensity_m = std::sqrt( (1.0 /Npoints) * Density2 / Physics::q_e / Physics::q_e );
 
     this->calcDebyeLength();
 
     m << "gammaz = "  << gammaz << endl;
     m << "hr_scaled = " << hr_scaled << endl;
-
+    m << "coupling constant= " << getCouplingConstant() << endl;
+    
     this->fsolver_m->runSolver();    
 
     gather(this->pcontainer_m->E, this->fcontainer_m->getE(), this->pcontainer_m->R);
-    
+
+    spaceChargeEFieldCheck();
+
     IpplTimings::stopTimer(SolveTimer);
 }
 

src/PartBunch/PartBunch.hpp

@@ -62,6 +62,14 @@ diagnostics(): calculate statistics and maybe write tp h5 and stat files
 
 #include "Algorithms/PartData.h"
 
+
+template <typename T>
+KOKKOS_INLINE_FUNCTION typename T::value_type L2Norm(T& x) {
+    return sqrt(dot(x, x).apply());
+}
+
+
+
 using view_type = typename ippl::detail::ViewType<ippl::Vector<double, 3>, 1>::view_type;
 
 
@@ -383,10 +391,11 @@ public:
     */
 
     double getCouplingConstant() const {
-        return 1.0;
+        return couplingConstant_m;
     }
-
+    
     void setCouplingConstant(double c) {
+        couplingConstant_m = c;
     }
 
     void calcBeamParameters();
@@ -850,6 +859,10 @@ public:
     double calculateAngle(double x, double y) {
         return 0.0;
     }
+
+    // Sanity check functions
+    void spaceChargeEFieldCheck();
+
 };
 
 /* \todo