Changeset 4571da for src/interface/interface_particles.cpp

Timestamp:

Apr 27, 2012, 11:34:57 PM (14 years ago)

Author:

Julian Iseringhausen <isering@…>

Children:

1a92cf

Parents:

b2154a3

Message:

vmg: Implement fourth-order discretization of the Poisson equation.

git-svn-id: ​https://svn.version.fz-juelich.de/scafacos/trunk@1762 5161e1c8-67bf-11de-9fd5-51895aff932f

File:

• src/interface/interface_particles.cpp (modified) (9 diffs)

src/interface/interface_particles.cpp

@@ -20 +20 @@
 #endif
 
+#include <algorithm>
 #include <cmath>
 #include <cstring>
@@ -85 +86 @@
     for (index.Y()=0; index.Y()<grid.Local().Size().Y(); ++index.Y())
       for (index.Z()=0; index.Z()<grid.Local().Size().Z(); ++index.Z())
-        grid(index + grid.Local().Begin()) = 4.0 * Math::pi * particle_grid.GetVal(index + particle_grid.Local().Begin());
+        grid(index + grid.Local().Begin()) = 4.0 * Math::pi * particle_grid.GetVal(index + particle_grid.Local().Begin());
 
 #ifdef DEBUG_OUTPUT
@@ -99 +100 @@
 void InterfaceParticles::ExportSolution(Grid& grid)
 {
-  Index index;
+  Index i;
   vmg_float length;
-  Vector dist_vec;
 
 #ifdef DEBUG_OUTPUT
@@ -107 +107 @@
   vmg_float e_long = 0.0;
   vmg_float e_self = 0.0;
-  vmg_float e_short = 0.0;
-  vmg_float e_shortcorr = 0.0;
+  vmg_float e_short_peak = 0.0;
+  vmg_float e_short_spline = 0.0;
 #endif
 
@@ -125 +125 @@
   InterpolatePolynomial ip(interpolation_degree);
 
-  const vmg_float r_cut = (near_field_cells+0.5) * grid.Extent().MeshWidth().Max();
-
-  /*
-   * Initialize potential
-   */
-  for (vmg_int i=0; i<num_particles_local; ++i) {
-    p[i] = 0.0;
-    for (vmg_int j=0;j<3;++j)
-      f[3*i+j] = 0.;
-  }
+  const vmg_float r_cut = near_field_cells * grid.Extent().MeshWidth().Max();
+
+  // Initialize field
+  std::fill(f, f+3*num_particles_local, 0.0);
 
   /*
@@ -143 +137 @@
   Grid& particle_grid = comm.GetParticleGrid();
 
-  for (index.X()=0; index.X()<grid.Local().Size().X(); ++index.X())
-    for (index.Y()=0; index.Y()<grid.Local().Size().Y(); ++index.Y())
-      for (index.Z()=0; index.Z()<grid.Local().Size().Z(); ++index.Z())
-        particle_grid(index + particle_grid.Local().Begin()) = grid.GetVal(index + grid.Local().Begin());
+  for (i.X()=0; i.X()<grid.Local().Size().X(); ++i.X())
+    for (i.Y()=0; i.Y()<grid.Local().Size().Y(); ++i.Y())
+      for (i.Z()=0; i.Z()<grid.Local().Size().Z(); ++i.Z())
+        particle_grid(i + particle_grid.Local().Begin()) = grid.GetVal(i + grid.Local().Begin());
 
   comm.CommToGhosts(particle_grid);
@@ -170 +164 @@
 
 #ifdef DEBUG_OUTPUT
-    // TODO The scaling with 1/2 may not be correct. Check that.
-      e_long += 0.5 * (*p1)->Charge() * Helper::InterpolateTrilinear((*p1)->Pos(), particle_grid);
+      e_long += 0.5 * (*p1)->Charge() * ip.Evaluate((*p1)->Pos());
       e_self += 0.5 * (*p1)->Charge() * (*p1)->Charge() * spl.GetAntiDerivativeAtZero();
 #endif
 
-      for (index.X()=-1*near_field_cells-1; index.X()<=near_field_cells+1; ++index.X())
-        for (index.Y()=-1*near_field_cells-1; index.Y()<=near_field_cells+1; ++index.Y())
-          for (index.Z()=-1*near_field_cells-1; index.Z()<=near_field_cells+1; ++index.Z()) {
-
-            std::list<Particle::Particle*>& lc_2 = lc(*iter+index);
+      for (i.X()=-1*near_field_cells; i.X()<=near_field_cells; ++i.X())
+        for (i.Y()=-1*near_field_cells; i.Y()<=near_field_cells; ++i.Y())
+          for (i.Z()=-1*near_field_cells; i.Z()<=near_field_cells; ++i.Z()) {
+
+            std::list<Particle::Particle*>& lc_2 = lc(*iter+i);
 
             for (p2=lc_2.begin(); p2!=lc_2.end(); ++p2)
@@ -186 +179 @@
                 length = ((*p2)->Pos() - (*p1)->Pos()).Length();
 
-                //TODO Rewrite this equation more efficiently
                 if (length < r_cut) {
                   (*p1)->Pot() += (*p2)->Charge() / length * (1.0 + spl.EvaluatePotential(length));
 
 #ifdef DEBUG_OUTPUT
-                  e_short += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length;
-                  e_shortcorr -= 0.5 * (*p1)->Charge() * (*p2)->Charge() / length * spl.EvaluatePotential(length);
+                  e_short_peak += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length;
+                  e_short_spline += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length * spl.EvaluatePotential(length);
 #endif
                 }
@@ -205 +197 @@
   vmg_float* q = factory.GetObjectStorageArray<vmg_float>("PARTICLE_CHARGE_ARRAY");
 
-  for (int i=0; i<num_particles_local; ++i)
-    e += 0.5 * p[i] * q[i];
-
+  e_long = comm.GlobalSumRoot(e_long);
+  e_short_peak = comm.GlobalSumRoot(e_short_peak);
+  e_short_spline = comm.GlobalSumRoot(e_short_spline);
+  e_self = comm.GlobalSumRoot(e_self);
+
+  for (int j=0; j<num_particles_local; ++j)
+    e += 0.5 * p[j] * q[j];
   e = comm.GlobalSumRoot(e);
-  e_long = comm.GlobalSumRoot(e_long);
-  e_self = comm.GlobalSumRoot(e_self);
-  e_short = comm.GlobalSumRoot(e_short);
-  e_shortcorr = comm.GlobalSumRoot(e_shortcorr);
-
-  comm.PrintStringOnce("E_long:  %e", e_long);
-  comm.PrintStringOnce("E_short: %e", e_short);
-  comm.PrintStringOnce("E_shortcorr: %e", e_shortcorr);
-  comm.PrintStringOnce("E_short*: %e", e - e_long + e_self);
-  comm.PrintStringOnce("E_self:  %e", e_self);
-  comm.PrintStringOnce("E_total: %e", e);
+
+  comm.PrintStringOnce("E_long:         %e", e_long);
+  comm.PrintStringOnce("E_short_peak:   %e", e_short_peak);
+  comm.PrintStringOnce("E_short_spline: %e", e_short_spline);
+  comm.PrintStringOnce("E_self:         %e", e_self);
+  comm.PrintStringOnce("E_total:        %e", e);
+  comm.PrintStringOnce("E_total*:       %e", e_long + e_short_peak + e_short_spline - e_self);
 
 #endif /* DEBUG_OUTPUT */