Changeset 716da7 for src/interface/interface_particles.cpp

Timestamp:

Apr 24, 2012, 2:26:14 PM (14 years ago)

Author:

Julian Iseringhausen <isering@…>

Children:

b51c3b

Parents:

e3dbbf

Message:

Fix energy calculation.

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

File:

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

src/interface/interface_particles.cpp

@@ -25 +25 @@
 #include "base/helper.hpp"
 #include "base/index.hpp"
+#include "base/interpolate_polynomial.hpp"
 #include "base/linked_cell_list.hpp"
 #include "base/math.hpp"
@@ -101 +102 @@
   vmg_float length;
   Vector dist_vec;
+  InterpolatePolynomial ip(8);
 
 #ifdef DEBUG_OUTPUT
@@ -119 +121 @@
   const int& near_field_cells = factory.GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS");
   vmg_float* p = factory.GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");
-  vmg_float* f = factory.GetObjectStorageArray<vmg_float>("PARTICLE_FORCE_ARRAY");
+  vmg_float* f = factory.GetObjectStorageArray<vmg_float>("PARTICLE_FIELD_ARRAY");
 
   const vmg_float r_cut = (near_field_cells+0.5) * grid.Extent().MeshWidth().Max();
@@ -152 +154 @@
   for (p_iter=particles.begin(); p_iter!=particles.end(); ++p_iter) {
 
+    const Index global_index = (p_iter->Pos() - particle_grid.Extent().Begin()) / particle_grid.Extent().MeshWidth();
+    const Index local_index = global_index - particle_grid.Global().LocalBegin() + particle_grid.Local().Begin();
+
+    ip.ComputeCoefficients(particle_grid, local_index);
+
     // Interpolate long range part of potential minus self-interaction
-    p_iter->Pot() = Helper::InterpolateTrilinear(p_iter->Pos(), particle_grid)
-      - p_iter->Charge() * spl.GetAntiDerivativeAtZero();
+    p_iter->Pot() = ip.Evaluate(p_iter->Pos()) - p_iter->Charge() * spl.GetAntiDerivativeAtZero();
+    p_iter->Field() = ip.EvaluateGradient(p_iter->Pos());
 
 #ifdef DEBUG_OUTPUT
@@ -171 +178 @@
   Grid::iterator iter;
 
-  comm.CommLCListToGhosts(grid, lc);
+  comm.CommLCListToGhosts(lc);
 
   for (iter=lc.Iterators().Local().Begin(); iter!=lc.Iterators().Local().End(); ++iter) {
@@ -179 +186 @@
         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 (p2=lc_2.begin(); p2!=lc_2.end(); ++p2)
               if (*p1 != *p2) {
@@ -187 +196 @@
                 //TODO Rewrite this equation more efficiently
                 if (length < r_cut) {
-                  (*p1)->Pot() += (*p2)->Charge() / length * (1.0 - spl.EvaluatePotential(length));
+                  (*p1)->Pot() += (*p2)->Charge() / length * (1.0 + spl.EvaluatePotential(length));
 
 #ifdef DEBUG_OUTPUT
@@ -193 +202 @@
                   e_shortcorr -= 0.5 * (*p1)->Charge() * (*p2)->Charge() / length * spl.EvaluatePotential(length);
 #endif
-
                 }
               }