core: p3m: Factored out common local mesh calculation from p3m/dp3m

src/core/electrostatics_magnetostatics/p3m-common.cpp

@@ -278,4 +278,88 @@ double p3m_caf(int i, double x, int cao_value) {
  }
  }
 }
+
+void p3m_calc_local_ca_mesh(p3m_local_mesh &local_mesh,
+ const P3MParameters &params,
+ const LocalBox<double> &local_geo, double skin) {
+ int i;
+ int ind[3];
+ /* total skin size */
+ double full_skin[3];
+
+ for (i = 0; i < 3; i++)
+ full_skin[i] = params.cao_cut[i] + skin + params.additional_mesh[i];
+
+ /* inner left down grid point (global index) */
+ for (i = 0; i < 3; i++)
+ local_mesh.in_ld[i] =
+ (int)ceil(local_geo.my_left()[i] * params.ai[i] - params.mesh_off[i]);
+ /* inner up right grid point (global index) */
+ for (i = 0; i < 3; i++)
+ local_mesh.in_ur[i] =
+ (int)floor(local_geo.my_right()[i] * params.ai[i] - params.mesh_off[i]);
+
+ /* correct roundoff errors at boundary */
+ for (i = 0; i < 3; i++) {
+ if ((local_geo.my_right()[i] * params.ai[i] - params.mesh_off[i]) -
+ local_mesh.in_ur[i] <
+ ROUND_ERROR_PREC)
+ local_mesh.in_ur[i]--;
+ if (1.0 + (local_geo.my_left()[i] * params.ai[i] - params.mesh_off[i]) -
+ local_mesh.in_ld[i] <
+ ROUND_ERROR_PREC)
+ local_mesh.in_ld[i]--;
+ }
+ /* inner grid dimensions */
+ for (i = 0; i < 3; i++)
+ local_mesh.inner[i] = local_mesh.in_ur[i] - local_mesh.in_ld[i] + 1;
+ /* index of left down grid point in global mesh */
+ for (i = 0; i < 3; i++)
+ local_mesh.ld_ind[i] =
+ (int)ceil((local_geo.my_left()[i] - full_skin[i]) * params.ai[i] -
+ params.mesh_off[i]);
+ /* left down margin */
+ for (i = 0; i < 3; i++)
+ local_mesh.margin[i * 2] = local_mesh.in_ld[i] - local_mesh.ld_ind[i];
+ /* up right grid point */
+ for (i = 0; i < 3; i++)
+ ind[i] =
+ (int)floor((local_geo.my_right()[i] + full_skin[i]) * params.ai[i] -
+ params.mesh_off[i]);
+ /* correct roundoff errors at up right boundary */
+ for (i = 0; i < 3; i++)
+ if (((local_geo.my_right()[i] + full_skin[i]) * params.ai[i] -
+ params.mesh_off[i]) -
+ ind[i] ==
+ 0)
+ ind[i]--;
+ /* up right margin */
+ for (i = 0; i < 3; i++)
+ local_mesh.margin[(i * 2) + 1] = ind[i] - local_mesh.in_ur[i];
+
+ /* grid dimension */
+ local_mesh.size = 1;
+ for (i = 0; i < 3; i++) {
+ local_mesh.dim[i] = ind[i] - local_mesh.ld_ind[i] + 1;
+ local_mesh.size *= local_mesh.dim[i];
+ }
+ /* reduce inner grid indices from global to local */
+ for (i = 0; i < 3; i++)
+ local_mesh.in_ld[i] = local_mesh.margin[i * 2];
+ for (i = 0; i < 3; i++)
+ local_mesh.in_ur[i] = local_mesh.margin[i * 2] + local_mesh.inner[i];
+
+ local_mesh.q_2_off = local_mesh.dim[2] - params.cao;
+ local_mesh.q_21_off = local_mesh.dim[2] * (local_mesh.dim[1] - params.cao);
+}
+
+void p3m_calc_lm_ld_pos(p3m_local_mesh &local_mesh,
+ const P3MParameters &params) {
+ /* spatial position of left down mesh point */
+ for (int i = 0; i < 3; i++) {
+ local_mesh.ld_pos[i] =
+ (local_mesh.ld_ind[i] + params.mesh_off[i]) * params.a[i];
+ }
+}
+
 #endif /* defined(P3M) || defined(DP3M) */

src/core/electrostatics_magnetostatics/p3m-common.hpp

@@ -36,6 +36,8 @@
 
 #if defined(P3M) || defined(DP3M)
 
+#include "LocalBox.hpp"
+
 /** Error Codes for p3m tuning (version 2) */
 enum P3M_TUNE_ERROR {
  /** force evaluation failed */
@@ -176,6 +178,23 @@ double p3m_analytic_cotangent_sum(int n, double mesh_i, int cao);
  */
 double p3m_caf(int i, double x, int cao_value);
 
+/** Calculate properties of the local FFT mesh for the
+ * charge assignment process.
+ */
+void p3m_calc_local_ca_mesh(p3m_local_mesh &local_mesh,
+ const P3MParameters &params,
+ const LocalBox<double> &local_geo, double skin);
+
+/** Calculate the spatial position of the left down mesh
+ * point of the local mesh, to be stored in
+ * @ref p3m_local_mesh::ld_pos "ld_pos".
+ *
+ * Function called by @ref p3m_calc_local_ca_mesh() once and by
+ * @ref p3m_scaleby_box_l() whenever the box size changes.
+ */
+void p3m_calc_lm_ld_pos(p3m_local_mesh &local_mesh,
+ const P3MParameters &params);
+
 #endif /* P3M || DP3M */
 
 #endif /* _P3M_COMMON_H */

src/core/electrostatics_magnetostatics/p3m-dipolar.cpp

@@ -81,15 +81,6 @@ dp3m_data_struct dp3m;
  */
 static void dp3m_init_a_ai_cao_cut();
 
-/** Calculate for magnetic dipoles the spatial position of the left down mesh
- * point of the local mesh, to be stored in
- * @ref p3m_local_mesh::ld_pos "ld_pos".
- *
- * Function called by @ref dp3m_calc_local_ca_mesh() once and by
- * @ref dp3m_scaleby_box_l() whenever the box size changes.
- */
-static void dp3m_calc_lm_ld_pos();
-
 /** realloc charge assignment fields. */
 static void dp3m_realloc_ca_fields(int newsize);
 
@@ -98,11 +89,6 @@ static void dp3m_realloc_ca_fields(int newsize);
  */
 static bool dp3m_sanity_checks_boxl();
 
-/** Calculate properties of the local FFT mesh for the
- * charge assignment process.
- */
-static void dp3m_calc_local_ca_mesh();
-
 /** Interpolate the P-th order charge assignment function from
  * @cite hockney88a 5-189 (or 8-61). The following charge fractions
  * are also tabulated in @cite deserno98a @cite deserno98b.
@@ -306,7 +292,7 @@ void dp3m_init() {
  dp3m_realloc_ca_fields(CA_INCREMENT);
  }
 
- dp3m_calc_local_ca_mesh();
+ p3m_calc_local_ca_mesh(dp3m.local_mesh, dp3m.params, local_geo, skin);
 
  dp3m.sm.resize(comm_cart, dp3m.local_mesh);
 
@@ -1961,17 +1947,6 @@ double dp3m_rtbisection(double box_size, double prefac, double r_cut_iL,
 
 /************************************************************/
 
-void dp3m_calc_lm_ld_pos() {
- int i;
- for (i = 0; i < 3; i++) {
- dp3m.local_mesh.ld_pos[i] =
- (dp3m.local_mesh.ld_ind[i] + dp3m.params.mesh_off[i]) *
- dp3m.params.a[i];
- }
-}
-
-/************************************************************/
-
 void dp3m_init_a_ai_cao_cut() {
  int i;
  for (i = 0; i < 3; i++) {
@@ -1981,90 +1956,6 @@ void dp3m_init_a_ai_cao_cut() {
  }
 }
 
-/************************************************************/
-
-void dp3m_calc_local_ca_mesh() {
- int i;
- int ind[3];
- /* total skin size */
- double full_skin[3];
-
- for (i = 0; i < 3; i++)
- full_skin[i] =
- dp3m.params.cao_cut[i] + skin + dp3m.params.additional_mesh[i];
-
- /* inner left down grid point (global index) */
- for (i = 0; i < 3; i++)
- dp3m.local_mesh.in_ld[i] = (int)ceil(
- local_geo.my_left()[i] * dp3m.params.ai[i] - dp3m.params.mesh_off[i]);
- /* inner up right grid point (global index) */
- for (i = 0; i < 3; i++)
- dp3m.local_mesh.in_ur[i] = (int)floor(
- local_geo.my_right()[i] * dp3m.params.ai[i] - dp3m.params.mesh_off[i]);
-
- /* correct roundoff errors at boundary */
- for (i = 0; i < 3; i++) {
- if ((local_geo.my_right()[i] * dp3m.params.ai[i] -
- dp3m.params.mesh_off[i]) -
- dp3m.local_mesh.in_ur[i] <
- ROUND_ERROR_PREC)
- dp3m.local_mesh.in_ur[i]--;
- if (1.0 +
- (local_geo.my_left()[i] * dp3m.params.ai[i] -
- dp3m.params.mesh_off[i]) -
- dp3m.local_mesh.in_ld[i] <
- ROUND_ERROR_PREC)
- dp3m.local_mesh.in_ld[i]--;
- }
- /* inner grid dimensions */
- for (i = 0; i < 3; i++)
- dp3m.local_mesh.inner[i] =
- dp3m.local_mesh.in_ur[i] - dp3m.local_mesh.in_ld[i] + 1;
- /* index of left down grid point in global mesh */
- for (i = 0; i < 3; i++)
- dp3m.local_mesh.ld_ind[i] =
- (int)ceil((local_geo.my_left()[i] - full_skin[i]) * dp3m.params.ai[i] -
- dp3m.params.mesh_off[i]);
- /* spacial position of left down mesh point */
- dp3m_calc_lm_ld_pos();
- /* left down margin */
- for (i = 0; i < 3; i++)
- dp3m.local_mesh.margin[i * 2] =
- dp3m.local_mesh.in_ld[i] - dp3m.local_mesh.ld_ind[i];
- /* up right grid point */
- for (i = 0; i < 3; i++)
- ind[i] = (int)floor((local_geo.my_right()[i] + full_skin[i]) *
- dp3m.params.ai[i] -
- dp3m.params.mesh_off[i]);
- /* correct roundoff errors at up right boundary */
- for (i = 0; i < 3; i++)
- if (((local_geo.my_right()[i] + full_skin[i]) * dp3m.params.ai[i] -
- dp3m.params.mesh_off[i]) -
- ind[i] ==
- 0)
- ind[i]--;
- /* up right margin */
- for (i = 0; i < 3; i++)
- dp3m.local_mesh.margin[(i * 2) + 1] = ind[i] - dp3m.local_mesh.in_ur[i];
-
- /* grid dimension */
- dp3m.local_mesh.size = 1;
- for (i = 0; i < 3; i++) {
- dp3m.local_mesh.dim[i] = ind[i] - dp3m.local_mesh.ld_ind[i] + 1;
- dp3m.local_mesh.size *= dp3m.local_mesh.dim[i];
- }
- /* reduce inner grid indices from global to local */
- for (i = 0; i < 3; i++)
- dp3m.local_mesh.in_ld[i] = dp3m.local_mesh.margin[i * 2];
- for (i = 0; i < 3; i++)
- dp3m.local_mesh.in_ur[i] =
- dp3m.local_mesh.margin[i * 2] + dp3m.local_mesh.inner[i];
-
- dp3m.local_mesh.q_2_off = dp3m.local_mesh.dim[2] - dp3m.params.cao;
- dp3m.local_mesh.q_21_off =
- dp3m.local_mesh.dim[2] * (dp3m.local_mesh.dim[1] - dp3m.params.cao);
-}
-
 /*****************************************************************************/
 
 bool dp3m_sanity_checks_boxl() {
@@ -2153,7 +2044,7 @@ void dp3m_scaleby_box_l() {
  dp3m.params.r_cut = dp3m.params.r_cut_iL * box_geo.length()[0];
  dp3m.params.alpha = dp3m.params.alpha_L * (1. / box_geo.length()[0]);
  dp3m_init_a_ai_cao_cut();
- dp3m_calc_lm_ld_pos();
+ p3m_calc_lm_ld_pos(dp3m.local_mesh, dp3m.params);
  dp3m_sanity_checks_boxl();
 
  dp3m_calc_influence_function_force();