diff --git a/src/core/electrostatics_magnetostatics/p3m-common.cpp b/src/core/electrostatics_magnetostatics/p3m-common.cpp index c49ed644c37..359cdd241ed 100644 --- a/src/core/electrostatics_magnetostatics/p3m-common.cpp +++ b/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 ¶ms, + const LocalBox &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 ¶ms) { + /* 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) */ diff --git a/src/core/electrostatics_magnetostatics/p3m-common.hpp b/src/core/electrostatics_magnetostatics/p3m-common.hpp index 0bf0e358060..f29fa76b726 100644 --- a/src/core/electrostatics_magnetostatics/p3m-common.hpp +++ b/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 ¶ms, + const LocalBox &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 ¶ms); + #endif /* P3M || DP3M */ #endif /* _P3M_COMMON_H */ diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp index 888d4269b86..12c57f9e910 100644 --- a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp +++ b/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(); diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp index afec9d01820..f0b90f45443 100644 --- a/src/core/electrostatics_magnetostatics/p3m.cpp +++ b/src/core/electrostatics_magnetostatics/p3m.cpp @@ -90,14 +90,6 @@ p3m_data_struct p3m; */ static void p3m_init_a_ai_cao_cut(); -/** 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 length changes. - */ -static void p3m_calc_lm_ld_pos(); - /** Calculates the dipole term */ static double p3m_calc_dipole_term(bool force_flag, bool energy_flag, const ParticleRange &particles); @@ -114,11 +106,6 @@ static bool p3m_sanity_checks_system(const Utils::Vector3i &grid); */ static bool p3m_sanity_checks_boxl(); -/** Calculates properties of the local FFT mesh for the charge assignment - * process. - */ -static void p3m_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. @@ -250,7 +237,7 @@ void p3m_init() { p3m_realloc_ca_fields(CA_INCREMENT); #endif - p3m_calc_local_ca_mesh(); + p3m_calc_local_ca_mesh(p3m.local_mesh, p3m.params, local_geo, skin); p3m.sm.resize(comm_cart, p3m.local_mesh); @@ -1815,93 +1802,6 @@ void p3m_tune_aliasing_sums(int nx, int ny, int nz, const int mesh[3], } } -void p3m_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] = p3m.params.cao_cut[i] + skin + p3m.params.additional_mesh[i]; - - /* inner left down grid point (global index) */ - for (i = 0; i < 3; i++) - p3m.local_mesh.in_ld[i] = (int)ceil( - local_geo.my_left()[i] * p3m.params.ai[i] - p3m.params.mesh_off[i]); - /* inner up right grid point (global index) */ - for (i = 0; i < 3; i++) - p3m.local_mesh.in_ur[i] = (int)floor( - local_geo.my_right()[i] * p3m.params.ai[i] - p3m.params.mesh_off[i]); - - /* correct roundoff errors at boundary */ - for (i = 0; i < 3; i++) { - if ((local_geo.my_right()[i] * p3m.params.ai[i] - p3m.params.mesh_off[i]) - - p3m.local_mesh.in_ur[i] < - ROUND_ERROR_PREC) - p3m.local_mesh.in_ur[i]--; - if (1.0 + - (local_geo.my_left()[i] * p3m.params.ai[i] - - p3m.params.mesh_off[i]) - - p3m.local_mesh.in_ld[i] < - ROUND_ERROR_PREC) - p3m.local_mesh.in_ld[i]--; - } - /* inner grid dimensions */ - for (i = 0; i < 3; i++) - p3m.local_mesh.inner[i] = - p3m.local_mesh.in_ur[i] - p3m.local_mesh.in_ld[i] + 1; - /* index of left down grid point in global mesh */ - for (i = 0; i < 3; i++) - p3m.local_mesh.ld_ind[i] = - (int)ceil((local_geo.my_left()[i] - full_skin[i]) * p3m.params.ai[i] - - p3m.params.mesh_off[i]); - /* left down margin */ - for (i = 0; i < 3; i++) - p3m.local_mesh.margin[i * 2] = - p3m.local_mesh.in_ld[i] - p3m.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]) * p3m.params.ai[i] - - p3m.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]) * p3m.params.ai[i] - - p3m.params.mesh_off[i]) - - ind[i] == - 0) - ind[i]--; - /* up right margin */ - for (i = 0; i < 3; i++) - p3m.local_mesh.margin[(i * 2) + 1] = ind[i] - p3m.local_mesh.in_ur[i]; - - /* grid dimension */ - p3m.local_mesh.size = 1; - for (i = 0; i < 3; i++) { - p3m.local_mesh.dim[i] = ind[i] - p3m.local_mesh.ld_ind[i] + 1; - p3m.local_mesh.size *= p3m.local_mesh.dim[i]; - } - /* reduce inner grid indices from global to local */ - for (i = 0; i < 3; i++) - p3m.local_mesh.in_ld[i] = p3m.local_mesh.margin[i * 2]; - for (i = 0; i < 3; i++) - p3m.local_mesh.in_ur[i] = - p3m.local_mesh.margin[i * 2] + p3m.local_mesh.inner[i]; - - p3m.local_mesh.q_2_off = p3m.local_mesh.dim[2] - p3m.params.cao; - p3m.local_mesh.q_21_off = - p3m.local_mesh.dim[2] * (p3m.local_mesh.dim[1] - p3m.params.cao); -} - -void p3m_calc_lm_ld_pos() { - int i; - /* spatial position of left down mesh point */ - for (i = 0; i < 3; i++) { - p3m.local_mesh.ld_pos[i] = - (p3m.local_mesh.ld_ind[i] + p3m.params.mesh_off[i]) * p3m.params.a[i]; - } -} - void p3m_init_a_ai_cao_cut() { int i; for (i = 0; i < 3; i++) { @@ -2002,7 +1902,7 @@ void p3m_scaleby_box_l() { p3m.params.r_cut = p3m.params.r_cut_iL * box_geo.length()[0]; p3m.params.alpha = p3m.params.alpha_L * (1. / box_geo.length()[0]); p3m_init_a_ai_cao_cut(); - p3m_calc_lm_ld_pos(); + p3m_calc_lm_ld_pos(p3m.local_mesh, p3m.params); p3m_sanity_checks_boxl(); p3m_calc_influence_function_force(); p3m_calc_influence_function_energy();