From 649fad9904fab7266e39735feb3959b93f51d768 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Wed, 22 Jan 2020 15:09:11 +0100 Subject: [PATCH 01/15] core: p3m: Reorder forward decls --- .../electrostatics_magnetostatics/p3m.cpp | 171 ++++++++---------- 1 file changed, 75 insertions(+), 96 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp index 4f028eae22d..59f050502b7 100644 --- a/src/core/electrostatics_magnetostatics/p3m.cpp +++ b/src/core/electrostatics_magnetostatics/p3m.cpp @@ -137,19 +137,6 @@ static void p3m_calc_influence_function_force(); */ static void p3m_calc_influence_function_energy(); -/** Calculate the aliasing sums for the optimal influence function. - * - * Calculate the aliasing sums in the nominator and denominator of - * the expression for the optimal influence function (see - * @cite hockney88a : 8-22, p. 275). - * - * \param n n-vector for which the aliasing sum is to be performed. - * \param nominator aliasing sums in the nominator. - * \retval denominator aliasing sum in the denominator - */ -double p3m_perform_aliasing_sums_force(int n[3], double nominator[3]); -double p3m_perform_aliasing_sums_energy(int n[3]); - /*@}*/ /** @name P3M tuning helper functions */ @@ -188,17 +175,6 @@ static void p3m_tune_aliasing_sums(int nx, int ny, int nz, const int mesh[3], double alpha_L_i, double *alias1, double *alias2); -/** Template parameterized calculation of the charge assignment to be called by - * wrapper. - * \tparam cao charge assignment order. - */ -template -static void p3m_do_charge_assign(const ParticleRange &particles); - -template -void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, - int cp_cnt); - p3m_data_struct::p3m_data_struct() { /* local_mesh is uninitialized */ /* sm is uninitialized */ @@ -382,78 +358,6 @@ void p3m_interpolate_charge_assignment_function() { } } -/* Template wrapper for p3m_do_charge_assign() */ -void p3m_charge_assign(const ParticleRange &particles) { - switch (p3m.params.cao) { - case 1: - p3m_do_charge_assign<1>(particles); - break; - case 2: - p3m_do_charge_assign<2>(particles); - break; - case 3: - p3m_do_charge_assign<3>(particles); - break; - case 4: - p3m_do_charge_assign<4>(particles); - break; - case 5: - p3m_do_charge_assign<5>(particles); - break; - case 6: - p3m_do_charge_assign<6>(particles); - break; - case 7: - p3m_do_charge_assign<7>(particles); - break; - } -} - -/** Assign the charges */ -template void p3m_do_charge_assign(const ParticleRange &particles) { - /* charged particle counter, charge fraction counter */ - int cp_cnt = 0; - /* prepare local FFT mesh */ - for (int i = 0; i < p3m.local_mesh.size; i++) - p3m.rs_mesh[i] = 0.0; - - for (auto &p : particles) { - if (p.p.q != 0.0) { - p3m_do_assign_charge(p.p.q, p.r.p, cp_cnt); - cp_cnt++; - } - } - - p3m_shrink_wrap_charge_grid(cp_cnt); -} - -/* Template wrapper for p3m_do_assign_charge() */ -void p3m_assign_charge(double q, const Utils::Vector3d &real_pos, int cp_cnt) { - switch (p3m.params.cao) { - case 1: - p3m_do_assign_charge<1>(q, real_pos, cp_cnt); - break; - case 2: - p3m_do_assign_charge<2>(q, real_pos, cp_cnt); - break; - case 3: - p3m_do_assign_charge<3>(q, real_pos, cp_cnt); - break; - case 4: - p3m_do_assign_charge<4>(q, real_pos, cp_cnt); - break; - case 5: - p3m_do_assign_charge<5>(q, real_pos, cp_cnt); - break; - case 6: - p3m_do_assign_charge<6>(q, real_pos, cp_cnt); - break; - case 7: - p3m_do_assign_charge<7>(q, real_pos, cp_cnt); - break; - } -} - template void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, int cp_cnt) { @@ -543,6 +447,81 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, } } +/** Template parameterized calculation of the charge assignment to be called by + * wrapper. + * \tparam cao charge assignment order. + */ +template void p3m_do_charge_assign(const ParticleRange &particles) { + /* charged particle counter, charge fraction counter */ + int cp_cnt = 0; + /* prepare local FFT mesh */ + for (int i = 0; i < p3m.local_mesh.size; i++) + p3m.rs_mesh[i] = 0.0; + + for (auto &p : particles) { + if (p.p.q != 0.0) { + p3m_do_assign_charge(p.p.q, p.r.p, cp_cnt); + cp_cnt++; + } + } + + p3m_shrink_wrap_charge_grid(cp_cnt); +} + +/* Template wrapper for p3m_do_charge_assign() */ +void p3m_charge_assign(const ParticleRange &particles) { + switch (p3m.params.cao) { + case 1: + p3m_do_charge_assign<1>(particles); + break; + case 2: + p3m_do_charge_assign<2>(particles); + break; + case 3: + p3m_do_charge_assign<3>(particles); + break; + case 4: + p3m_do_charge_assign<4>(particles); + break; + case 5: + p3m_do_charge_assign<5>(particles); + break; + case 6: + p3m_do_charge_assign<6>(particles); + break; + case 7: + p3m_do_charge_assign<7>(particles); + break; + } +} + +/* Template wrapper for p3m_do_assign_charge() */ +void p3m_assign_charge(double q, const Utils::Vector3d &real_pos, int cp_cnt) { + switch (p3m.params.cao) { + case 1: + p3m_do_assign_charge<1>(q, real_pos, cp_cnt); + break; + case 2: + p3m_do_assign_charge<2>(q, real_pos, cp_cnt); + break; + case 3: + p3m_do_assign_charge<3>(q, real_pos, cp_cnt); + break; + case 4: + p3m_do_assign_charge<4>(q, real_pos, cp_cnt); + break; + case 5: + p3m_do_assign_charge<5>(q, real_pos, cp_cnt); + break; + case 6: + p3m_do_assign_charge<6>(q, real_pos, cp_cnt); + break; + case 7: + p3m_do_assign_charge<7>(q, real_pos, cp_cnt); + break; + } +} + void p3m_shrink_wrap_charge_grid(int n_charges) { /* we do not really want to export these */ if (n_charges < p3m.ca_num) From 06c6ee81873f8b875e144369bdedf7b2d60adf46 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Wed, 22 Jan 2020 15:54:38 +0100 Subject: [PATCH 02/15] maintainer: Fixed syntax of p3m benchmark --- maintainer/benchmarks/p3m.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/maintainer/benchmarks/p3m.py b/maintainer/benchmarks/p3m.py index 43fe781596d..2c9bd7b63ab 100644 --- a/maintainer/benchmarks/p3m.py +++ b/maintainer/benchmarks/p3m.py @@ -56,7 +56,7 @@ import espressomd -from espressomd import electrostatics +from espressomd import electrostatics, minimize_energy if args.visualizer: from espressomd import visualization @@ -140,10 +140,10 @@ energy = system.analysis.energy() print("Before Minimization: E_total = {}".format(energy["total"])) -system.minimize_energy.init(f_max=1000, gamma=30.0, +minimize_energy.steepest_descent(system, f_max=1000, gamma=30.0, + max_steps=1000, max_displacement=0.05) +minimize_energy.steepest_descent(system, f_max=1000, gamma=30.0, max_steps=1000, max_displacement=0.05) -system.minimize_energy.minimize() -system.minimize_energy.minimize() energy = system.analysis.energy() print("After Minimization: E_total = {}".format(energy["total"])) From 324e84143abe6b927b18f859f2623c582641b309 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Wed, 22 Jan 2020 16:57:37 +0100 Subject: [PATCH 03/15] core: p3m: Use less memory bandwidth --- .../electrostatics_magnetostatics/p3m.cpp | 35 ++++++++++++------- src/utils/include/utils/Vector.hpp | 2 +- 2 files changed, 23 insertions(+), 14 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp index 59f050502b7..1d706af8414 100644 --- a/src/core/electrostatics_magnetostatics/p3m.cpp +++ b/src/core/electrostatics_magnetostatics/p3m.cpp @@ -58,6 +58,7 @@ using Utils::strcat_alloc; #include #include +#include #include #include #include @@ -371,9 +372,6 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, // make sure we have enough space if (cp_cnt >= p3m.ca_num) p3m_realloc_ca_fields(cp_cnt + 1); - // do it here, since p3m_realloc_ca_fields may change the address of - // p3m.ca_frac - double *cur_ca_frac = p3m.ca_frac.data() + cao * cao * cao * cp_cnt; for (int d = 0; d < 3; d++) { /* particle position in mesh coordinates */ @@ -436,15 +434,18 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, for (int i2 = 0; i2 < cao; i2++) { auto const cur_ca_frac_val = q * tmp1 * w_z[i2]; p3m.rs_mesh[q_ind] += cur_ca_frac_val; - /* store current ca frac */ - if (cp_cnt >= 0) - *(cur_ca_frac++) = cur_ca_frac_val; q_ind++; } q_ind += p3m.local_mesh.q_2_off; } q_ind += p3m.local_mesh.q_21_off; } + + if (cp_cnt >= 0) { + boost::copy(w_x, p3m.ca_frac.begin() + (3 * cp_cnt + 0) * cao); + boost::copy(w_y, p3m.ca_frac.begin() + (3 * cp_cnt + 1) * cao); + boost::copy(w_z, p3m.ca_frac.begin() + (3 * cp_cnt + 2) * cao); + } } /** Template parameterized calculation of the charge assignment to be called by @@ -532,25 +533,33 @@ void p3m_shrink_wrap_charge_grid(int n_charges) { template static void P3M_assign_forces(double force_prefac, const ParticleRange &particles) { + using Utils::make_const_span; + using Utils::Span; + using Utils::Vector; + /* charged particle counter, charge fraction counter */ int cp_cnt = 0; - int cf_cnt = 0; - /* index, index jumps for rs_mesh array */ - int q_ind = 0; for (auto &p : particles) { auto const q = p.p.q; if (q != 0.0) { - q_ind = p3m.ca_fmp[cp_cnt]; + Span w_x(p3m.ca_frac.data() + (3 * cp_cnt + 0) * cao, cao); + const Vector w_y( + make_const_span(p3m.ca_frac.data() + (3 * cp_cnt + 1) * cao, cao)); + const Vector w_z( + make_const_span(p3m.ca_frac.data() + (3 * cp_cnt + 2) * cao, cao)); + + /* index, index jumps for rs_mesh array */ + auto q_ind = p3m.ca_fmp[cp_cnt]; for (int i0 = 0; i0 < cao; i0++) { + auto const fx = q * w_x[i0]; for (int i1 = 0; i1 < cao; i1++) { + auto const fxy = w_y[i1] * fx; for (int i2 = 0; i2 < cao; i2++) { - p.f.f -= force_prefac * p3m.ca_frac[cf_cnt] * + p.f.f -= force_prefac * fxy * w_z[i2] * Utils::Vector3d{p3m.E_mesh[0][q_ind], p3m.E_mesh[1][q_ind], p3m.E_mesh[2][q_ind]}; - q_ind++; - cf_cnt++; } q_ind += p3m.local_mesh.q_2_off; } diff --git a/src/utils/include/utils/Vector.hpp b/src/utils/include/utils/Vector.hpp index bdebb05c061..fa4b0158ee7 100644 --- a/src/utils/include/utils/Vector.hpp +++ b/src/utils/include/utils/Vector.hpp @@ -83,7 +83,7 @@ template class Vector : public Array { } template - Vector(InputIterator first, InputIterator last) { + Vector(InputIterator first, InputIterator last) : Base() { if (std::distance(first, last) == N) { std::copy_n(first, N, begin()); } else { From c44b3fa9d67770e5a9bc08912e0f8573341169d2 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 23 Jan 2020 15:48:13 +0100 Subject: [PATCH 04/15] core: p3m: Removed interpolation function interpolation --- maintainer/benchmarks/p3m.py | 4 +- .../p3m-common.hpp | 8 +- .../p3m-dipolar.cpp | 141 ++++-------------- .../p3m-dipolar.hpp | 3 - .../electrostatics_magnetostatics/p3m.cpp | 77 ++-------- .../electrostatics_magnetostatics/p3m.hpp | 12 +- src/python/espressomd/electrostatics.pxd | 29 +--- src/python/espressomd/electrostatics.pyx | 21 +-- src/python/espressomd/magnetostatics.pxd | 1 - src/python/espressomd/magnetostatics.pyx | 1 - src/python/espressomd/p3m_common.pxd | 2 - testsuite/python/analyze_energy.py | 2 +- testsuite/python/coulomb_cloud_wall.py | 21 +-- testsuite/python/electrostaticInteractions.py | 6 +- 14 files changed, 57 insertions(+), 271 deletions(-) diff --git a/maintainer/benchmarks/p3m.py b/maintainer/benchmarks/p3m.py index 2c9bd7b63ab..7c75363f921 100644 --- a/maintainer/benchmarks/p3m.py +++ b/maintainer/benchmarks/p3m.py @@ -141,9 +141,9 @@ energy = system.analysis.energy() print("Before Minimization: E_total = {}".format(energy["total"])) minimize_energy.steepest_descent(system, f_max=1000, gamma=30.0, - max_steps=1000, max_displacement=0.05) + max_steps=1000, max_displacement=0.05) minimize_energy.steepest_descent(system, f_max=1000, gamma=30.0, - max_steps=1000, max_displacement=0.05) + max_steps=1000, max_displacement=0.05) energy = system.analysis.energy() print("After Minimization: E_total = {}".format(energy["total"])) diff --git a/src/core/electrostatics_magnetostatics/p3m-common.hpp b/src/core/electrostatics_magnetostatics/p3m-common.hpp index 3eddddf3400..361f4952928 100644 --- a/src/core/electrostatics_magnetostatics/p3m-common.hpp +++ b/src/core/electrostatics_magnetostatics/p3m-common.hpp @@ -112,8 +112,6 @@ typedef struct { double mesh_off[3] = {P3M_MESHOFF, P3M_MESHOFF, P3M_MESHOFF}; /** charge assignment order ([0,7]). */ int cao = 0; - /** number of interpolation points for charge assignment function */ - int inter = P3M_N_INTERPOL; /** accuracy of the actual parameter set. */ double accuracy = 0.0; @@ -131,8 +129,6 @@ typedef struct { /** unscaled @ref P3MParameters::r_cut_iL "r_cut_iL" for use with fast * inline functions only */ double r_cut = -1.; - /** full size of the interpolated assignment function */ - int inter2 = 0; /** number of points unto which a single charge is interpolated, i.e. * p3m.cao^3 */ int cao3 = 0; @@ -142,8 +138,8 @@ typedef struct { template void serialize(Archive &ar, long int) { ar &tuning &alpha_L &r_cut_iL &mesh; - ar &mesh_off &cao &inter &accuracy &epsilon &cao_cut; - ar &a &ai &alpha &r_cut &inter2 &cao3 &additional_mesh; + ar &mesh_off &cao &accuracy &epsilon &cao_cut; + ar &a &ai &alpha &r_cut &cao3 &additional_mesh; } } P3MParameters; diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp index 4d80c4d3dbb..92861564327 100644 --- a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp +++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp @@ -90,12 +90,6 @@ static void dp3m_realloc_ca_fields(int newsize); */ static bool dp3m_sanity_checks_boxl(); -/** 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. - */ -static void dp3m_interpolate_dipole_assignment_function(); - /** Shift the mesh points by mesh/2 */ static void dp3m_calc_meshift(); @@ -298,9 +292,6 @@ void dp3m_init() { /* fix box length dependent constants */ dp3m_scaleby_box_l(); - if (dp3m.params.inter > 0) - dp3m_interpolate_dipole_assignment_function(); - dp3m.pos_shift = std::floor((dp3m.params.cao - 1) / 2.0) - (dp3m.params.cao % 2) / 2.0; @@ -419,9 +410,6 @@ void dp3m_set_tune_params(double r_cut, int mesh, int cao, double alpha, if (accuracy >= 0) dp3m.params.accuracy = accuracy; - - if (n_interpol != -1) - dp3m.params.inter = n_interpol; } /*****************************************************************************/ @@ -486,40 +474,6 @@ int dp3m_set_eps(double eps) { return ES_OK; } -int dp3m_set_ninterpol(int n) { - if (n < 0) - return ES_ERROR; - - dp3m.params.inter = n; - - mpi_bcast_coulomb_params(); - - return ES_OK; -} - -/*****************************************************************************/ - -void dp3m_interpolate_dipole_assignment_function() { - double dInterpol = 0.5 / (double)dp3m.params.inter; - int i; - long j; - - if (dp3m.params.inter == 0) - return; - - dp3m.params.inter2 = 2 * dp3m.params.inter + 1; - - for (i = 0; i < dp3m.params.cao; i++) { - /* allocate memory for interpolation array */ - dp3m.int_caf[i].resize(2 * dp3m.params.inter + 1); - - /* loop over all interpolation points */ - for (j = -dp3m.params.inter; j <= dp3m.params.inter; j++) - dp3m.int_caf[i][j + dp3m.params.inter] = - Utils::bspline(i, j * dInterpol, dp3m.params.cao); - } -} - void dp3m_dipole_assign(const ParticleRange &particles) { /* magnetic particle counter, dipole fraction counter */ int cp_cnt = 0; @@ -551,8 +505,6 @@ void dp3m_assign_dipole(double const real_pos[3], double mu, int nmp; /* distance to nearest mesh point */ double dist[3]; - /* index for caf interpolation grid */ - int arg[3]; /* index, index jumps for dp3m.rs_mesh array */ int q_ind = 0; double cur_ca_frac_val, *cur_ca_frac; @@ -564,74 +516,39 @@ void dp3m_assign_dipole(double const real_pos[3], double mu, // dp3m.ca_frac cur_ca_frac = dp3m.ca_frac.data() + dp3m.params.cao3 * cp_cnt; - if (dp3m.params.inter == 0) { - for (d = 0; d < 3; d++) { - /* particle position in mesh coordinates */ - pos = ((real_pos[d] - dp3m.local_mesh.ld_pos[d]) * dp3m.params.ai[d]) - - dp3m.pos_shift; - /* nearest mesh point */ - nmp = (int)pos; - /* distance to nearest mesh point */ - dist[d] = (pos - nmp) - 0.5; - /* 3d-array index of nearest mesh point */ - q_ind = (d == 0) ? nmp : nmp + dp3m.local_mesh.dim[d] * q_ind; - } - - if (cp_cnt >= 0) - dp3m.ca_fmp[cp_cnt] = q_ind; - - for (i0 = 0; i0 < dp3m.params.cao; i0++) { - tmp0 = Utils::bspline(i0, dist[0], dp3m.params.cao); - for (i1 = 0; i1 < dp3m.params.cao; i1++) { - tmp1 = tmp0 * Utils::bspline(i1, dist[1], dp3m.params.cao); - for (i2 = 0; i2 < dp3m.params.cao; i2++) { - cur_ca_frac_val = tmp1 * Utils::bspline(i2, dist[2], dp3m.params.cao); - if (cp_cnt >= 0) - *(cur_ca_frac++) = cur_ca_frac_val; - if (mu != 0.0) { - dp3m.rs_mesh_dip[0][q_ind] += dip[0] * cur_ca_frac_val; - dp3m.rs_mesh_dip[1][q_ind] += dip[1] * cur_ca_frac_val; - dp3m.rs_mesh_dip[2][q_ind] += dip[2] * cur_ca_frac_val; - } - q_ind++; - } - q_ind += dp3m.local_mesh.q_2_off; - } - q_ind += dp3m.local_mesh.q_21_off; - } - } else { + for (d = 0; d < 3; d++) { /* particle position in mesh coordinates */ - for (d = 0; d < 3; d++) { - pos = ((real_pos[d] - dp3m.local_mesh.ld_pos[d]) * dp3m.params.ai[d]) - - dp3m.pos_shift; - nmp = (int)pos; - arg[d] = (int)((pos - nmp) * dp3m.params.inter2); - /* for the first dimension, q_ind is always zero, so this shifts correctly - */ - q_ind = nmp + dp3m.local_mesh.dim[d] * q_ind; - } - if (cp_cnt >= 0) - dp3m.ca_fmp[cp_cnt] = q_ind; - - for (i0 = 0; i0 < dp3m.params.cao; i0++) { - tmp0 = dp3m.int_caf[i0][arg[0]]; - for (i1 = 0; i1 < dp3m.params.cao; i1++) { - tmp1 = tmp0 * dp3m.int_caf[i1][arg[1]]; - for (i2 = 0; i2 < dp3m.params.cao; i2++) { - cur_ca_frac_val = tmp1 * dp3m.int_caf[i2][arg[2]]; - if (cp_cnt >= 0) - *(cur_ca_frac++) = cur_ca_frac_val; - if (mu != 0.0) { - dp3m.rs_mesh_dip[0][q_ind] += dip[0] * cur_ca_frac_val; - dp3m.rs_mesh_dip[1][q_ind] += dip[1] * cur_ca_frac_val; - dp3m.rs_mesh_dip[2][q_ind] += dip[2] * cur_ca_frac_val; - } - q_ind++; + pos = ((real_pos[d] - dp3m.local_mesh.ld_pos[d]) * dp3m.params.ai[d]) - + dp3m.pos_shift; + /* nearest mesh point */ + nmp = (int)pos; + /* distance to nearest mesh point */ + dist[d] = (pos - nmp) - 0.5; + /* 3d-array index of nearest mesh point */ + q_ind = (d == 0) ? nmp : nmp + dp3m.local_mesh.dim[d] * q_ind; + } + + if (cp_cnt >= 0) + dp3m.ca_fmp[cp_cnt] = q_ind; + + for (i0 = 0; i0 < dp3m.params.cao; i0++) { + tmp0 = Utils::bspline(i0, dist[0], dp3m.params.cao); + for (i1 = 0; i1 < dp3m.params.cao; i1++) { + tmp1 = tmp0 * Utils::bspline(i1, dist[1], dp3m.params.cao); + for (i2 = 0; i2 < dp3m.params.cao; i2++) { + cur_ca_frac_val = tmp1 * Utils::bspline(i2, dist[2], dp3m.params.cao); + if (cp_cnt >= 0) + *(cur_ca_frac++) = cur_ca_frac_val; + if (mu != 0.0) { + dp3m.rs_mesh_dip[0][q_ind] += dip[0] * cur_ca_frac_val; + dp3m.rs_mesh_dip[1][q_ind] += dip[1] * cur_ca_frac_val; + dp3m.rs_mesh_dip[2][q_ind] += dip[2] * cur_ca_frac_val; } - q_ind += dp3m.local_mesh.q_2_off; + q_ind++; } - q_ind += dp3m.local_mesh.q_21_off; + q_ind += dp3m.local_mesh.q_2_off; } + q_ind += dp3m.local_mesh.q_21_off; } } diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp index 4cffb17a23a..9165e7285fe 100644 --- a/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp +++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp @@ -115,9 +115,6 @@ void dp3m_set_tune_params(double r_cut, int mesh, int cao, double alpha, int dp3m_set_params(double r_cut, int mesh, int cao, double alpha, double accuracy); -/** @copydoc p3m_set_ninterpol */ -int dp3m_set_ninterpol(int n); - /** @copydoc p3m_set_mesh_offset */ int dp3m_set_mesh_offset(double x, double y, double z); diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp index 1d706af8414..a7a9868a829 100644 --- a/src/core/electrostatics_magnetostatics/p3m.cpp +++ b/src/core/electrostatics_magnetostatics/p3m.cpp @@ -106,12 +106,6 @@ static bool p3m_sanity_checks_system(const Utils::Vector3i &grid); */ static bool p3m_sanity_checks_boxl(); -/** 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. - */ -static void p3m_interpolate_charge_assignment_function(); - /** Shift the mesh points by mesh/2 */ static void p3m_calc_meshift(); @@ -227,9 +221,6 @@ void p3m_init() { /* fix box length dependent constants */ p3m_scaleby_box_l(); - if (p3m.params.inter > 0) - p3m_interpolate_charge_assignment_function(); - /* position offset for calc. of first meshpoint */ p3m.pos_shift = std::floor((p3m.params.cao - 1) / 2.0) - (p3m.params.cao % 2) / 2.0; @@ -239,7 +230,7 @@ void p3m_init() { } void p3m_set_tune_params(double r_cut, const int mesh[3], int cao, double alpha, - double accuracy, int n_interpol) { + double accuracy) { if (r_cut >= 0) { p3m.params.r_cut = r_cut; p3m.params.r_cut_iL = r_cut * (1. / box_geo.length()[0]); @@ -261,9 +252,6 @@ void p3m_set_tune_params(double r_cut, const int mesh[3], int cao, double alpha, if (accuracy >= 0) p3m.params.accuracy = accuracy; - - if (n_interpol != -1) - p3m.params.inter = n_interpol; } /*@}*/ @@ -327,46 +315,11 @@ int p3m_set_eps(double eps) { return ES_OK; } -int p3m_set_ninterpol(int n) { - if (n < 0) - return ES_ERROR; - - p3m.params.inter = n; - - mpi_bcast_coulomb_params(); - - return ES_OK; -} - -void p3m_interpolate_charge_assignment_function() { - double dInterpol = 0.5 / (double)p3m.params.inter; - int i; - long j; - - if (p3m.params.inter == 0) - return; - - p3m.params.inter2 = 2 * p3m.params.inter + 1; - - for (i = 0; i < p3m.params.cao; i++) { - /* allocate memory for interpolation array */ - p3m.int_caf[i].resize(2 * p3m.params.inter + 1); - - /* loop over all interpolation points */ - for (j = -p3m.params.inter; j <= p3m.params.inter; j++) - p3m.int_caf[i][j + p3m.params.inter] = - Utils::bspline(i, j * dInterpol, p3m.params.cao); - } -} - template void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, int cp_cnt) { - auto const inter = not(p3m.params.inter == 0); /* distance to nearest mesh point */ double dist[3]; - /* index for caf interpolation grid */ - int arg[3]; /* index, index jumps for rs_mesh array */ int q_ind = 0; // make sure we have enough space @@ -383,12 +336,8 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, /* 3d-array index of nearest mesh point */ q_ind = (d == 0) ? nmp : nmp + p3m.local_mesh.dim[d] * q_ind; - if (!inter) - /* distance to nearest mesh point */ - dist[d] = (pos - nmp) - 0.5; - else - /* distance to nearest mesh point for interpolation */ - arg[d] = (int)((pos - nmp) * p3m.params.inter2); + /* distance to nearest mesh point */ + dist[d] = (pos - nmp) - 0.5; #ifdef ADDITIONAL_CHECKS if (pos < -skin * p3m.params.ai[d]) { @@ -411,20 +360,12 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, Utils::Array w_x, w_y, w_z; - if (!inter) { - for (int i = 0; i < cao; i++) { - using Utils::bspline; + for (int i = 0; i < cao; i++) { + using Utils::bspline; - w_x[i] = bspline(i, dist[0]); - w_y[i] = bspline(i, dist[1]); - w_z[i] = bspline(i, dist[2]); - } - } else { - for (int i = 0; i < cao; i++) { - w_x[i] = p3m.int_caf[i][arg[0]]; - w_y[i] = p3m.int_caf[i][arg[1]]; - w_z[i] = p3m.int_caf[i][arg[2]]; - } + w_x[i] = bspline(i, dist[0]); + w_y[i] = bspline(i, dist[1]); + w_z[i] = bspline(i, dist[2]); } for (int i0 = 0; i0 < cao; i0++) { @@ -815,7 +756,7 @@ void p3m_realloc_ca_fields(int newsize) { newsize = CA_INCREMENT; p3m.ca_num = newsize; - p3m.ca_frac.resize(p3m.params.cao3 * p3m.ca_num); + p3m.ca_frac.resize(3 * p3m.params.cao * p3m.ca_num); p3m.ca_fmp.resize(p3m.ca_num); } diff --git a/src/core/electrostatics_magnetostatics/p3m.hpp b/src/core/electrostatics_magnetostatics/p3m.hpp index fd733a52071..05dd3bf898e 100644 --- a/src/core/electrostatics_magnetostatics/p3m.hpp +++ b/src/core/electrostatics_magnetostatics/p3m.hpp @@ -68,9 +68,6 @@ struct p3m_data_struct { /** square of sum of charges (only on master node). */ double square_sum_q; - /** interpolation of the charge assignment function. */ - std::array, 7> int_caf; - /** position shift for calc. of first assignment mesh point. */ double pos_shift; /** help variable for calculation of aliasing sums */ @@ -226,10 +223,9 @@ inline void p3m_add_pair_force(double q1q2, Utils::Vector3d const &d, * @param[in] cao @copybrief P3MParameters::cao * @param[in] alpha @copybrief P3MParameters::alpha * @param[in] accuracy @copybrief P3MParameters::accuracy - * @param[in] n_interpol @copybrief P3MParameters::inter */ void p3m_set_tune_params(double r_cut, const int mesh[3], int cao, double alpha, - double accuracy, int n_interpol); + double accuracy); /** Set custom parameters * @@ -256,12 +252,6 @@ int p3m_set_mesh_offset(double x, double y, double z); */ int p3m_set_eps(double eps); -/** Set @ref P3MParameters::inter "inter" parameter - * - * @param[in] n @copybrief P3MParameters::inter - */ -int p3m_set_ninterpol(int n); - /** Calculate real space contribution of Coulomb pair energy. */ inline double p3m_pair_energy(double chgfac, double dist) { if (dist < p3m.params.r_cut && dist != 0) { diff --git a/src/python/espressomd/electrostatics.pxd b/src/python/espressomd/electrostatics.pxd index ba0c0ba02d8..328a70201e0 100644 --- a/src/python/espressomd/electrostatics.pxd +++ b/src/python/espressomd/electrostatics.pxd @@ -65,31 +65,13 @@ IF ELECTROSTATICS: void deactivate_method() IF P3M: - cdef extern from "electrostatics_magnetostatics/p3m-common.hpp": - ctypedef struct P3MParameters: - double alpha_L - double r_cut_iL - int mesh[3] - double mesh_off[3] - int cao - int inter - double accuracy - double epsilon - double cao_cut[3] - double a[3] - double ai[3] - double alpha - double r_cut - int inter2 - int cao3 - double additional_mesh[3] + from p3m_common cimport P3MParameters cdef extern from "electrostatics_magnetostatics/p3m.hpp": int p3m_set_params(double r_cut, int * mesh, int cao, double alpha, double accuracy) - void p3m_set_tune_params(double r_cut, int mesh[3], int cao, double alpha, double accuracy, int n_interpol) + void p3m_set_tune_params(double r_cut, int mesh[3], int cao, double alpha, double accuracy) int p3m_set_mesh_offset(double x, double y, double z) int p3m_set_eps(double eps) - int p3m_set_ninterpol(int n) int p3m_adaptive_tune(char ** log) ctypedef struct p3m_data_struct: @@ -153,19 +135,18 @@ IF ELECTROSTATICS: return p3m_set_params(r_cut, mesh, cao, alpha, accuracy) - cdef inline python_p3m_set_tune_params(p_r_cut, p_mesh, p_cao, p_alpha, p_accuracy, p_n_interpol): + cdef inline python_p3m_set_tune_params(p_r_cut, p_mesh, p_cao, p_alpha, p_accuracy): # cdef inline python_p3m_set_tune_params(): cdef int mesh[3] cdef double r_cut cdef int cao cdef double alpha cdef double accuracy - cdef int n_interpol r_cut = p_r_cut cao = p_cao alpha = p_alpha accuracy = p_accuracy - n_interpol = p_n_interpol + if is_valid_type(p_mesh, int): mesh[0] = p_mesh mesh[1] = p_mesh @@ -173,7 +154,7 @@ IF ELECTROSTATICS: else: mesh = p_mesh - p3m_set_tune_params(r_cut, mesh, cao, alpha, accuracy, n_interpol) + p3m_set_tune_params(r_cut, mesh, cao, alpha, accuracy) cdef extern from "electrostatics_magnetostatics/debye_hueckel.hpp": ctypedef struct Debye_hueckel_params: diff --git a/src/python/espressomd/electrostatics.pyx b/src/python/espressomd/electrostatics.pyx index ceb0bf59717..884d931740f 100644 --- a/src/python/espressomd/electrostatics.pyx +++ b/src/python/espressomd/electrostatics.pyx @@ -267,10 +267,6 @@ IF P3M == 1: or self._params["epsilon"] == "metallic"): raise ValueError("epsilon should be a double or 'metallic'") - if not (self._params["inter"] == default_params["inter"] - or self._params["inter"] >= 0): - raise ValueError("inter should be a positive integer") - if not (self._params["mesh_off"] == default_params["mesh_off"] or len(self._params) != 3): raise ValueError( @@ -283,14 +279,13 @@ IF P3M == 1: def valid_keys(self): return ["mesh", "cao", "accuracy", "epsilon", "alpha", "r_cut", - "prefactor", "tune", "check_neutrality", "inter"] + "prefactor", "tune", "check_neutrality"] def required_keys(self): return ["prefactor", "accuracy"] def default_params(self): return {"cao": 0, - "inter": -1, "r_cut": -1, "alpha": 0, "accuracy": 0, @@ -319,8 +314,6 @@ IF P3M == 1: # which resets r_cut if lb is zero. OK. # Sets eps, bcast p3m_set_eps(self._params["epsilon"]) - # Sets ninterpol, bcast - p3m_set_ninterpol(self._params["inter"]) python_p3m_set_mesh_offset(self._params["mesh_off"]) def _tune(self): @@ -329,8 +322,7 @@ IF P3M == 1: self._params["mesh"], self._params["cao"], -1.0, - self._params["accuracy"], - self._params["inter"]) + self._params["accuracy"]) resp = python_p3m_adaptive_tune() if resp: raise Exception( @@ -419,10 +411,6 @@ IF P3M == 1: raise ValueError( "epsilon should be a double or 'metallic'") - if not (self._params["inter"] == default_params["inter"] - or self._params["inter"] > 0): - raise ValueError("inter should be a positive integer") - if not (self._params["mesh_off"] == default_params["mesh_off"] or len(self._params) != 3): raise ValueError( @@ -437,7 +425,6 @@ IF P3M == 1: def default_params(self): return {"cao": 0, - "inter": -1, "r_cut": -1, "alpha": 0, "accuracy": 0, @@ -460,8 +447,7 @@ IF P3M == 1: self._params["mesh"], self._params["cao"], -1.0, - self._params["accuracy"], - self._params["inter"]) + self._params["accuracy"]) resp = python_p3m_adaptive_tune() if resp: raise Exception( @@ -485,7 +471,6 @@ IF P3M == 1: self._params["alpha"], self._params["accuracy"]) p3m_set_eps(self._params["epsilon"]) - p3m_set_ninterpol(self._params["inter"]) python_p3m_set_mesh_offset(self._params["mesh_off"]) handle_errors("p3m gpu init") diff --git a/src/python/espressomd/magnetostatics.pxd b/src/python/espressomd/magnetostatics.pxd index 4373c0a509e..14ef94f12a0 100644 --- a/src/python/espressomd/magnetostatics.pxd +++ b/src/python/espressomd/magnetostatics.pxd @@ -63,7 +63,6 @@ IF DP3M == 1: void dp3m_set_tune_params(double r_cut, int mesh, int cao, double alpha, double accuracy, int n_interpol) int dp3m_set_mesh_offset(double x, double y, double z) int dp3m_set_eps(double eps) - int dp3m_set_ninterpol(int n) int dp3m_adaptive_tune(char ** log) int dp3m_deactivate() diff --git a/src/python/espressomd/magnetostatics.pyx b/src/python/espressomd/magnetostatics.pyx index a5c18babd14..9c3651ce25a 100644 --- a/src/python/espressomd/magnetostatics.pyx +++ b/src/python/espressomd/magnetostatics.pyx @@ -165,7 +165,6 @@ IF DP3M == 1: def _set_params_in_es_core(self): self.set_magnetostatics_prefactor() dp3m_set_eps(self._params["epsilon"]) - dp3m_set_ninterpol(self._params["inter"]) self.python_dp3m_set_mesh_offset(self._params["mesh_off"]) self.python_dp3m_set_params( self._params["r_cut"], self._params["mesh"], diff --git a/src/python/espressomd/p3m_common.pxd b/src/python/espressomd/p3m_common.pxd index c6c99b4be71..8b9bb28f12f 100644 --- a/src/python/espressomd/p3m_common.pxd +++ b/src/python/espressomd/p3m_common.pxd @@ -24,7 +24,6 @@ IF P3M == 1 or DP3M == 1: int mesh[3] double mesh_off[3] int cao - int inter double accuracy double epsilon double cao_cut[3] @@ -32,6 +31,5 @@ IF P3M == 1 or DP3M == 1: double ai[3] double alpha double r_cut - int inter2 int cao3 double additional_mesh[3] diff --git a/testsuite/python/analyze_energy.py b/testsuite/python/analyze_energy.py index 03000fb49ed..a9768b078fe 100644 --- a/testsuite/python/analyze_energy.py +++ b/testsuite/python/analyze_energy.py @@ -180,7 +180,7 @@ def test_electrostatics(self): self.assertAlmostEqual(energy["kinetic"], 0., delta=1e-7) self.assertAlmostEqual(energy["bonded"], 0., delta=1e-7) self.assertAlmostEqual(energy["non_bonded"], 0, delta=1e-7) - self.assertAlmostEqual(energy["coulomb"], u_p3m, delta=1e-7) + self.assertAlmostEqual(energy["coulomb"], u_p3m, delta=1e-5) self.system.part[0].q = 0 self.system.part[1].q = 0 self.system.part[0].pos = [1, 2, 2] diff --git a/testsuite/python/coulomb_cloud_wall.py b/testsuite/python/coulomb_cloud_wall.py index e2fa89e9053..ab72e11b1b2 100644 --- a/testsuite/python/coulomb_cloud_wall.py +++ b/testsuite/python/coulomb_cloud_wall.py @@ -88,26 +88,9 @@ def compare(self, method_name, energy=True, prefactor=None): # Tests for individual methods @utx.skipIfMissingFeatures(["P3M"]) - def test_p3m_direct_caf(self): + def test_p3m_direct(self): """ - This checks P3M with using the charge assignment - function (window function) directly by setting the - `inter` parameter to zero. - - """ - - self.S.actors.add( - espressomd.electrostatics.P3M( - prefactor=3, r_cut=1.001, accuracy=1e-3, - mesh=64, cao=7, alpha=2.70746, tune=False, inter=0)) - self.S.integrator.run(0) - self.compare("p3m", energy=True, prefactor=3) - - @utx.skipIfMissingFeatures(["P3M"]) - def test_p3m_interpolated_caf(self): - """ - This checks P3M with using an interpolated charge assignment - function (window function), which is the default. + This checks P3M. """ diff --git a/testsuite/python/electrostaticInteractions.py b/testsuite/python/electrostaticInteractions.py index 13eb6c5047c..7dcbf6301dc 100644 --- a/testsuite/python/electrostaticInteractions.py +++ b/testsuite/python/electrostaticInteractions.py @@ -93,13 +93,13 @@ def test_p3m(self): tune=False) self.system.actors.add(p3m) self.assertAlmostEqual(self.system.analysis.energy()['coulomb'], - p3m_energy) + p3m_energy, places=5) # need to update forces self.system.integrator.run(0) np.testing.assert_allclose(np.copy(self.system.part[0].f), - [p3m_force, 0, 0], atol=1E-5) + [p3m_force, 0, 0], atol=1E-4) np.testing.assert_allclose(np.copy(self.system.part[1].f), - [-p3m_force, 0, 0], atol=1E-10) + [-p3m_force, 0, 0], atol=1E-5) self.system.actors.remove(p3m) def test_dh(self): From 815b106a7933adbc03159e9ff44932c38996f6b4 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 23 Jan 2020 15:58:00 +0100 Subject: [PATCH 05/15] core: p3m: Switched to Utils::get_linear_index for index calc --- .../electrostatics_magnetostatics/fft.cpp | 2 +- .../electrostatics_magnetostatics/fft.hpp | 2 +- .../p3m-common.hpp | 2 +- .../electrostatics_magnetostatics/p3m.cpp | 32 ++++++------------- .../p3m_send_mesh.cpp | 12 +++---- .../p3m_send_mesh.hpp | 10 +++--- 6 files changed, 25 insertions(+), 35 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/fft.cpp b/src/core/electrostatics_magnetostatics/fft.cpp index 45a5620f958..7499281a655 100644 --- a/src/core/electrostatics_magnetostatics/fft.cpp +++ b/src/core/electrostatics_magnetostatics/fft.cpp @@ -489,7 +489,7 @@ void calc_2d_grid(int n, int grid[3]) { } } // namespace -int fft_init(int const *ca_mesh_dim, int const *ca_mesh_margin, +int fft_init(const Utils::Vector3i &ca_mesh_dim, int const *ca_mesh_margin, int *global_mesh_dim, double *global_mesh_off, int *ks_pnum, fft_data_struct &fft, const Utils::Vector3i &grid, const boost::mpi::communicator &comm) { diff --git a/src/core/electrostatics_magnetostatics/fft.hpp b/src/core/electrostatics_magnetostatics/fft.hpp index 8694974487d..4bf414768c8 100644 --- a/src/core/electrostatics_magnetostatics/fft.hpp +++ b/src/core/electrostatics_magnetostatics/fft.hpp @@ -186,7 +186,7 @@ struct fft_data_struct { * \param comm MPI communicator. * \return Maximal size of local fft mesh (needed for allocation of ca_mesh). */ -int fft_init(int const *ca_mesh_dim, int const *ca_mesh_margin, +int fft_init(const Utils::Vector3i &ca_mesh_dim, int const *ca_mesh_margin, int *global_mesh_dim, double *global_mesh_off, int *ks_pnum, fft_data_struct &fft, const Utils::Vector3i &grid, const boost::mpi::communicator &comm); diff --git a/src/core/electrostatics_magnetostatics/p3m-common.hpp b/src/core/electrostatics_magnetostatics/p3m-common.hpp index 361f4952928..6b823ca0a95 100644 --- a/src/core/electrostatics_magnetostatics/p3m-common.hpp +++ b/src/core/electrostatics_magnetostatics/p3m-common.hpp @@ -71,7 +71,7 @@ enum P3M_TUNE_ERROR { typedef struct { /* local mesh characterization. */ /** dimension (size) of local mesh. */ - int dim[3]; + Utils::Vector3i dim; /** number of local mesh points. */ int size; /** index of lower left corner of the diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp index a7a9868a829..bcb1ea9c84a 100644 --- a/src/core/electrostatics_magnetostatics/p3m.cpp +++ b/src/core/electrostatics_magnetostatics/p3m.cpp @@ -63,6 +63,7 @@ using Utils::strcat_alloc; #include #include #include +#include /************************************************ * variables @@ -320,40 +321,27 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, int cp_cnt) { /* distance to nearest mesh point */ double dist[3]; - /* index, index jumps for rs_mesh array */ - int q_ind = 0; // make sure we have enough space if (cp_cnt >= p3m.ca_num) p3m_realloc_ca_fields(cp_cnt + 1); + /* nearest mesh point */ + Utils::Vector3i nmp; + for (int d = 0; d < 3; d++) { /* particle position in mesh coordinates */ auto const pos = ((real_pos[d] - p3m.local_mesh.ld_pos[d]) * p3m.params.ai[d]) - p3m.pos_shift; - /* nearest mesh point */ - auto const nmp = (int)pos; - /* 3d-array index of nearest mesh point */ - q_ind = (d == 0) ? nmp : nmp + p3m.local_mesh.dim[d] * q_ind; + + nmp[d] = (int)pos; /* distance to nearest mesh point */ - dist[d] = (pos - nmp) - 0.5; - -#ifdef ADDITIONAL_CHECKS - if (pos < -skin * p3m.params.ai[d]) { - fprintf(stderr, "%d: rs_mesh underflow! (pos %f)\n", this_node, - real_pos[d]); - fprintf(stderr, "%d: allowed coordinates: %f - %f\n", this_node, - local_geo.my_left()[d] - skin, local_geo.my_right()[d] + skin); - } - if ((nmp + cao) > p3m.local_mesh.dim[d]) { - fprintf(stderr, "%d: rs_mesh overflow! (pos %f, nmp=%d)\n", this_node, - real_pos[d], nmp); - fprintf(stderr, "%d: allowed coordinates: %f - %f\n", this_node, - local_geo.my_left()[d] - skin, local_geo.my_right()[d] + skin); - } -#endif + dist[d] = (pos - nmp[d]) - 0.5; } + /* 3d-array index of nearest mesh point */ + auto q_ind = Utils::get_linear_index(nmp, p3m.local_mesh.dim, + Utils::MemoryOrder::ROW_MAJOR); if (cp_cnt >= 0) p3m.ca_fmp[cp_cnt] = q_ind; diff --git a/src/core/electrostatics_magnetostatics/p3m_send_mesh.cpp b/src/core/electrostatics_magnetostatics/p3m_send_mesh.cpp index 674a83c5a2a..8984a53989e 100644 --- a/src/core/electrostatics_magnetostatics/p3m_send_mesh.cpp +++ b/src/core/electrostatics_magnetostatics/p3m_send_mesh.cpp @@ -103,7 +103,7 @@ void p3m_send_mesh::resize(const boost::mpi::communicator &comm, void p3m_send_mesh::gather_grid(Utils::Span meshes, const boost::mpi::communicator &comm, - const int *dim) { + const Utils::Vector3i &dim) { auto const node_neighbors = Utils::Mpi::cart_neighbors<3>(comm); send_grid.resize(max * meshes.size()); recv_grid.resize(max * meshes.size()); @@ -116,7 +116,7 @@ void p3m_send_mesh::gather_grid(Utils::Span meshes, if (s_size[s_dir] > 0) for (size_t i = 0; i < meshes.size(); i++) { fft_pack_block(meshes[i], send_grid.data() + i * s_size[s_dir], - s_ld[s_dir], s_dim[s_dir], dim, 1); + s_ld[s_dir], s_dim[s_dir], dim.data(), 1); } /* communication */ @@ -133,7 +133,7 @@ void p3m_send_mesh::gather_grid(Utils::Span meshes, if (r_size[r_dir] > 0) { for (size_t i = 0; i < meshes.size(); i++) { p3m_add_block(recv_grid.data() + i * r_size[r_dir], meshes[i], - r_ld[r_dir], r_dim[r_dir], dim); + r_ld[r_dir], r_dim[r_dir], dim.data()); } } } @@ -141,7 +141,7 @@ void p3m_send_mesh::gather_grid(Utils::Span meshes, void p3m_send_mesh::spread_grid(Utils::Span meshes, const boost::mpi::communicator &comm, - const int *dim) { + const Utils::Vector3i &dim) { auto const node_neighbors = Utils::Mpi::cart_neighbors<3>(comm); send_grid.resize(max * meshes.size()); recv_grid.resize(max * meshes.size()); @@ -154,7 +154,7 @@ void p3m_send_mesh::spread_grid(Utils::Span meshes, if (r_size[r_dir] > 0) for (size_t i = 0; i < meshes.size(); i++) { fft_pack_block(meshes[i], send_grid.data() + i * r_size[r_dir], - r_ld[r_dir], r_dim[r_dir], dim, 1); + r_ld[r_dir], r_dim[r_dir], dim.data(), 1); } /* communication */ if (node_neighbors[r_dir] != comm.rank()) { @@ -170,7 +170,7 @@ void p3m_send_mesh::spread_grid(Utils::Span meshes, if (s_size[s_dir] > 0) { for (size_t i = 0; i < meshes.size(); i++) { fft_unpack_block(recv_grid.data() + i * s_size[s_dir], meshes[i], - s_ld[s_dir], s_dim[s_dir], dim, 1); + s_ld[s_dir], s_dim[s_dir], dim.data(), 1); } } } diff --git a/src/core/electrostatics_magnetostatics/p3m_send_mesh.hpp b/src/core/electrostatics_magnetostatics/p3m_send_mesh.hpp index 6197d90cd8c..d4e1106d1a0 100644 --- a/src/core/electrostatics_magnetostatics/p3m_send_mesh.hpp +++ b/src/core/electrostatics_magnetostatics/p3m_send_mesh.hpp @@ -66,15 +66,17 @@ class p3m_send_mesh { void resize(const boost::mpi::communicator &comm, const p3m_local_mesh &local_mesh); void gather_grid(Utils::Span meshes, - const boost::mpi::communicator &comm, const int *dim); + const boost::mpi::communicator &comm, + const Utils::Vector3i &dim); void gather_grid(double *mesh, const boost::mpi::communicator &comm, - const int *dim) { + const Utils::Vector3i &dim) { gather_grid(Utils::make_span(&mesh, 1), comm, dim); } void spread_grid(Utils::Span meshes, - const boost::mpi::communicator &comm, const int *dim); + const boost::mpi::communicator &comm, + const Utils::Vector3i &dim); void spread_grid(double *mesh, const boost::mpi::communicator &comm, - const int *dim) { + const Utils::Vector3i &dim) { spread_grid(Utils::make_span(&mesh, 1), comm, dim); } }; From 0bd823daaed30de28c25e406a63786853a6ac7e0 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Thu, 23 Jan 2020 16:05:56 +0100 Subject: [PATCH 06/15] core: p3m: Calculate pos_shift otf --- src/core/electrostatics_magnetostatics/p3m.cpp | 17 +++++++---------- src/core/electrostatics_magnetostatics/p3m.hpp | 2 -- 2 files changed, 7 insertions(+), 12 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp index bcb1ea9c84a..5174cd5d249 100644 --- a/src/core/electrostatics_magnetostatics/p3m.cpp +++ b/src/core/electrostatics_magnetostatics/p3m.cpp @@ -179,8 +179,6 @@ p3m_data_struct::p3m_data_struct() { sum_q2 = 0.0; square_sum_q = 0.0; - pos_shift = 0.0; - ca_num = 0; ks_pnum = 0; } @@ -222,10 +220,6 @@ void p3m_init() { /* fix box length dependent constants */ p3m_scaleby_box_l(); - /* position offset for calc. of first meshpoint */ - p3m.pos_shift = - std::floor((p3m.params.cao - 1) / 2.0) - (p3m.params.cao % 2) / 2.0; - p3m_count_charged_particles(); } } @@ -319,11 +313,11 @@ int p3m_set_eps(double eps) { template void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, int cp_cnt) { + /** position shift for calc. of first assignment mesh point. */ + static auto const pos_shift = std::floor((cao - 1) / 2.0) - (cao % 2) / 2.0; + /* distance to nearest mesh point */ double dist[3]; - // make sure we have enough space - if (cp_cnt >= p3m.ca_num) - p3m_realloc_ca_fields(cp_cnt + 1); /* nearest mesh point */ Utils::Vector3i nmp; @@ -332,7 +326,7 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, /* particle position in mesh coordinates */ auto const pos = ((real_pos[d] - p3m.local_mesh.ld_pos[d]) * p3m.params.ai[d]) - - p3m.pos_shift; + pos_shift; nmp[d] = (int)pos; @@ -343,6 +337,9 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, auto q_ind = Utils::get_linear_index(nmp, p3m.local_mesh.dim, Utils::MemoryOrder::ROW_MAJOR); + // make sure we have enough space + if (cp_cnt >= p3m.ca_num) + p3m_realloc_ca_fields(cp_cnt + 1); if (cp_cnt >= 0) p3m.ca_fmp[cp_cnt] = q_ind; diff --git a/src/core/electrostatics_magnetostatics/p3m.hpp b/src/core/electrostatics_magnetostatics/p3m.hpp index 05dd3bf898e..4d5893c926e 100644 --- a/src/core/electrostatics_magnetostatics/p3m.hpp +++ b/src/core/electrostatics_magnetostatics/p3m.hpp @@ -68,8 +68,6 @@ struct p3m_data_struct { /** square of sum of charges (only on master node). */ double square_sum_q; - /** position shift for calc. of first assignment mesh point. */ - double pos_shift; /** help variable for calculation of aliasing sums */ std::vector meshift_x; std::vector meshift_y; From 1cb9096c447f7051bbd7bc4d1711ebbbbc2b4e21 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Fri, 24 Jan 2020 09:47:01 +0100 Subject: [PATCH 07/15] core: p3m: Facotred out interpolation weight storage --- .../electrostatics_magnetostatics/elc.cpp | 13 +-- .../electrostatics_magnetostatics/p3m.cpp | 108 ++++++------------ .../electrostatics_magnetostatics/p3m.hpp | 21 ++-- .../p3m_interpolation.hpp | 72 ++++++++++++ 4 files changed, 120 insertions(+), 94 deletions(-) create mode 100644 src/core/electrostatics_magnetostatics/p3m_interpolation.hpp diff --git a/src/core/electrostatics_magnetostatics/elc.cpp b/src/core/electrostatics_magnetostatics/elc.cpp index a98126bbad6..e5ee7cac199 100644 --- a/src/core/electrostatics_magnetostatics/elc.cpp +++ b/src/core/electrostatics_magnetostatics/elc.cpp @@ -1301,7 +1301,7 @@ void assign_image_charge(const Particle &p) { auto const q_eff = elc_params.delta_mid_bot * p.p.q; auto const pos = Utils::Vector3d{p.r.p[0], p.r.p[1], -p.r.p[2]}; - p3m_assign_charge(q_eff, pos, -1); + p3m_assign_charge(q_eff, pos, nullptr); } if (p.r.p[2] > (elc_params.h - elc_params.space_layer)) { @@ -1309,27 +1309,24 @@ void assign_image_charge(const Particle &p) { auto const pos = Utils::Vector3d{p.r.p[0], p.r.p[1], 2 * elc_params.h - p.r.p[2]}; - p3m_assign_charge(q_eff, pos, -1); + p3m_assign_charge(q_eff, pos, nullptr); } } } // namespace void ELC_p3m_charge_assign_both(const ParticleRange &particles) { - /* charged particle counter, charge fraction counter */ - int cp_cnt = 0; + p3m.inter_weights.reset(p3m.params.cao); + /* prepare local FFT mesh */ for (int i = 0; i < p3m.local_mesh.size; i++) p3m.rs_mesh[i] = 0.0; for (auto &p : particles) { if (p.p.q != 0.0) { - p3m_assign_charge(p.p.q, p.r.p, cp_cnt); + p3m_assign_charge(p.p.q, p.r.p, &p3m.inter_weights); assign_image_charge(p); - - cp_cnt++; } } - p3m_shrink_wrap_charge_grid(cp_cnt); } void ELC_p3m_charge_assign_image(const ParticleRange &particles) { diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp index 5174cd5d249..b2dc4361c89 100644 --- a/src/core/electrostatics_magnetostatics/p3m.cpp +++ b/src/core/electrostatics_magnetostatics/p3m.cpp @@ -58,12 +58,10 @@ using Utils::strcat_alloc; #include #include -#include #include #include #include #include -#include /************************************************ * variables @@ -97,9 +95,6 @@ p3m_data_struct p3m; */ static void p3m_init_a_ai_cao_cut(); -/** Realloc charge assignment fields. */ -static void p3m_realloc_ca_fields(int newsize); - static bool p3m_sanity_checks_system(const Utils::Vector3i &grid); /** Checks for correctness for charges in P3M of the cao_cut, @@ -179,7 +174,6 @@ p3m_data_struct::p3m_data_struct() { sum_q2 = 0.0; square_sum_q = 0.0; - ca_num = 0; ks_pnum = 0; } @@ -198,10 +192,6 @@ void p3m_init() { * the cutoff for charge assignment p3m.params.cao_cut */ p3m_init_a_ai_cao_cut(); - /* initialize ca fields to size CA_INCREMENT: p3m.ca_frac and p3m.ca_fmp */ - p3m.ca_num = 0; - p3m_realloc_ca_fields(CA_INCREMENT); - p3m_calc_local_ca_mesh(p3m.local_mesh, p3m.params, local_geo, skin); p3m.sm.resize(comm_cart, p3m.local_mesh); @@ -312,7 +302,7 @@ int p3m_set_eps(double eps) { template void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, - int cp_cnt) { + p3m_interpolation_weights *inter_weights = nullptr) { /** position shift for calc. of first assignment mesh point. */ static auto const pos_shift = std::floor((cao - 1) / 2.0) - (cao % 2) / 2.0; @@ -337,12 +327,6 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, auto q_ind = Utils::get_linear_index(nmp, p3m.local_mesh.dim, Utils::MemoryOrder::ROW_MAJOR); - // make sure we have enough space - if (cp_cnt >= p3m.ca_num) - p3m_realloc_ca_fields(cp_cnt + 1); - if (cp_cnt >= 0) - p3m.ca_fmp[cp_cnt] = q_ind; - Utils::Array w_x, w_y, w_z; for (int i = 0; i < cao; i++) { @@ -353,6 +337,10 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, w_z[i] = bspline(i, dist[2]); } + if (inter_weights) { + inter_weights->store(q_ind, w_x, w_y, w_z); + } + for (int i0 = 0; i0 < cao; i0++) { auto const tmp0 = w_x[i0]; for (int i1 = 0; i1 < cao; i1++) { @@ -366,12 +354,6 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, } q_ind += p3m.local_mesh.q_21_off; } - - if (cp_cnt >= 0) { - boost::copy(w_x, p3m.ca_frac.begin() + (3 * cp_cnt + 0) * cao); - boost::copy(w_y, p3m.ca_frac.begin() + (3 * cp_cnt + 1) * cao); - boost::copy(w_z, p3m.ca_frac.begin() + (3 * cp_cnt + 2) * cao); - } } /** Template parameterized calculation of the charge assignment to be called by @@ -379,24 +361,21 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, * \tparam cao charge assignment order. */ template void p3m_do_charge_assign(const ParticleRange &particles) { - /* charged particle counter, charge fraction counter */ - int cp_cnt = 0; - /* prepare local FFT mesh */ - for (int i = 0; i < p3m.local_mesh.size; i++) - p3m.rs_mesh[i] = 0.0; - for (auto &p : particles) { if (p.p.q != 0.0) { - p3m_do_assign_charge(p.p.q, p.r.p, cp_cnt); - cp_cnt++; + p3m_do_assign_charge(p.p.q, p.r.p, &p3m.inter_weights); } } - - p3m_shrink_wrap_charge_grid(cp_cnt); } /* Template wrapper for p3m_do_charge_assign() */ void p3m_charge_assign(const ParticleRange &particles) { + p3m.inter_weights.reset(p3m.params.cao); + + /* prepare local FFT mesh */ + for (int i = 0; i < p3m.local_mesh.size; i++) + p3m.rs_mesh[i] = 0.0; + switch (p3m.params.cao) { case 1: p3m_do_charge_assign<1>(particles); @@ -423,38 +402,33 @@ void p3m_charge_assign(const ParticleRange &particles) { } /* Template wrapper for p3m_do_assign_charge() */ -void p3m_assign_charge(double q, const Utils::Vector3d &real_pos, int cp_cnt) { +void p3m_assign_charge(double q, const Utils::Vector3d &real_pos, + p3m_interpolation_weights *inter_weights) { switch (p3m.params.cao) { case 1: - p3m_do_assign_charge<1>(q, real_pos, cp_cnt); + p3m_do_assign_charge<1>(q, real_pos, inter_weights); break; case 2: - p3m_do_assign_charge<2>(q, real_pos, cp_cnt); + p3m_do_assign_charge<2>(q, real_pos, inter_weights); break; case 3: - p3m_do_assign_charge<3>(q, real_pos, cp_cnt); + p3m_do_assign_charge<3>(q, real_pos, inter_weights); break; case 4: - p3m_do_assign_charge<4>(q, real_pos, cp_cnt); + p3m_do_assign_charge<4>(q, real_pos, inter_weights); break; case 5: - p3m_do_assign_charge<5>(q, real_pos, cp_cnt); + p3m_do_assign_charge<5>(q, real_pos, inter_weights); break; case 6: - p3m_do_assign_charge<6>(q, real_pos, cp_cnt); + p3m_do_assign_charge<6>(q, real_pos, inter_weights); break; case 7: - p3m_do_assign_charge<7>(q, real_pos, cp_cnt); + p3m_do_assign_charge<7>(q, real_pos, inter_weights); break; } } -void p3m_shrink_wrap_charge_grid(int n_charges) { - /* we do not really want to export these */ - if (n_charges < p3m.ca_num) - p3m_realloc_ca_fields(n_charges); -} - /* Assign the forces obtained from k-space */ template static void P3M_assign_forces(double force_prefac, @@ -463,26 +437,30 @@ static void P3M_assign_forces(double force_prefac, using Utils::Span; using Utils::Vector; + assert(cao == p3m.inter_weights.cao()); + /* charged particle counter, charge fraction counter */ int cp_cnt = 0; for (auto &p : particles) { auto const q = p.p.q; if (q != 0.0) { - Span w_x(p3m.ca_frac.data() + (3 * cp_cnt + 0) * cao, cao); - const Vector w_y( - make_const_span(p3m.ca_frac.data() + (3 * cp_cnt + 1) * cao, cao)); - const Vector w_z( - make_const_span(p3m.ca_frac.data() + (3 * cp_cnt + 2) * cao, cao)); - - /* index, index jumps for rs_mesh array */ - auto q_ind = p3m.ca_fmp[cp_cnt]; + auto const pref = q * force_prefac; + + Span w_x, cw_y, cw_z; + int q_ind; + std::tie(q_ind, w_x, cw_y, cw_z) = p3m.inter_weights.load(cp_cnt); + + /* We only pin the directions that are used multiple times. */ + const Vector w_y(cw_y); + const Vector w_z(cw_z); + for (int i0 = 0; i0 < cao; i0++) { - auto const fx = q * w_x[i0]; + auto const fx = w_x[i0]; for (int i1 = 0; i1 < cao; i1++) { auto const fxy = w_y[i1] * fx; for (int i2 = 0; i2 < cao; i2++) { - p.f.f -= force_prefac * fxy * w_z[i2] * + p.f.f -= pref * fxy * w_z[i2] * Utils::Vector3d{p3m.E_mesh[0][q_ind], p3m.E_mesh[1][q_ind], p3m.E_mesh[2][q_ind]}; q_ind++; @@ -613,7 +591,7 @@ double p3m_calc_kspace_forces(bool force_flag, bool energy_flag, // Note: after these calls, the grids are in the order yzx and not xyz // anymore!!! - /* The dipol moment is only needed if we don't have metallic boundaries. */ + /* The dipole moment is only needed if we don't have metallic boundaries. */ auto const box_dipole = (p3m.params.epsilon != P3M_EPSILON_METALLIC) ? boost::make_optional(calc_dipole_moment( comm_cart, particles, box_geo)) @@ -731,20 +709,6 @@ double p3m_calc_kspace_forces(bool force_flag, bool energy_flag, return 0.0; } -/************************************************************/ - -void p3m_realloc_ca_fields(int newsize) { - newsize = ((newsize + CA_INCREMENT - 1) / CA_INCREMENT) * CA_INCREMENT; - if (newsize == p3m.ca_num) - return; - if (newsize < CA_INCREMENT) - newsize = CA_INCREMENT; - - p3m.ca_num = newsize; - p3m.ca_frac.resize(3 * p3m.params.cao * p3m.ca_num); - p3m.ca_fmp.resize(p3m.ca_num); -} - void p3m_calc_meshift() { p3m.meshift_x.resize(p3m.params.mesh[0]); p3m.meshift_y.resize(p3m.params.mesh[1]); diff --git a/src/core/electrostatics_magnetostatics/p3m.hpp b/src/core/electrostatics_magnetostatics/p3m.hpp index 4d5893c926e..e06a8de9c66 100644 --- a/src/core/electrostatics_magnetostatics/p3m.hpp +++ b/src/core/electrostatics_magnetostatics/p3m.hpp @@ -39,6 +39,7 @@ #include "fft.hpp" #include "p3m-common.hpp" +#include "p3m_interpolation.hpp" #include "p3m_send_mesh.hpp" #include @@ -81,12 +82,7 @@ struct p3m_data_struct { /** Energy optimised influence function (k-space) */ std::vector g_energy; - /** number of charged particles on the node. */ - int ca_num; - /** Charge fractions for mesh assignment. */ - std::vector ca_frac; - /** index of first mesh point for charge assignment. */ - std::vector ca_fmp; + p3m_interpolation_weights inter_weights; /** number of permutations in k_space */ int ks_pnum; @@ -178,15 +174,12 @@ void p3m_charge_assign(const ParticleRange &particles); * * @param[in] q %Particle charge * @param[in] real_pos %Particle position in real space - * @param[in] cp_cnt The running index, which may be smaller than 0, in - * which case the charge is assumed to be virtual and - * is not stored in the @ref p3m_data_struct::ca_frac - * "ca_frac" arrays + * @param[in] inter_weights The running index, which may be smaller than 0, + * in which case the charge is assumed to be virtual and is not stored in the + * @ref p3m_data_struct::ca_frac "ca_frac" arrays */ -void p3m_assign_charge(double q, const Utils::Vector3d &real_pos, int cp_cnt); - -/** Shrink wrap the charge grid */ -void p3m_shrink_wrap_charge_grid(int n_charges); +void p3m_assign_charge(double q, const Utils::Vector3d &real_pos, + p3m_interpolation_weights *inter_weights); /** Calculate real space contribution of Coulomb pair forces. */ inline void p3m_add_pair_force(double q1q2, Utils::Vector3d const &d, diff --git a/src/core/electrostatics_magnetostatics/p3m_interpolation.hpp b/src/core/electrostatics_magnetostatics/p3m_interpolation.hpp new file mode 100644 index 00000000000..36bc160e624 --- /dev/null +++ b/src/core/electrostatics_magnetostatics/p3m_interpolation.hpp @@ -0,0 +1,72 @@ +#ifndef ESPRESSO_P3M_INTERPOLATION_HPP +#define ESPRESSO_P3M_INTERPOLATION_HPP + +#include +#include + +#include + +#include +#include +#include + +struct p3m_interpolation_weights { + size_t m_cao = 0; + /** Charge fractions for mesh assignment. */ + std::vector ca_frac; + /** index of first mesh point for charge assignment. */ + std::vector ca_fmp; + + /** + * @brief Number of points in the cache. + * @return Number of points currently in the cache. + */ + auto size() const { return ca_fmp.size(); } + + /** + * @brief Charge assignment order the weights are for. + * @return The charge assigment order. + */ + auto cao() const { return m_cao; } + + /** + * @brief Push back weights for one point. + * + * @param q_ind Mesh index. + * @param w_x Weights in first direction. + * @param w_y Weights in second direction. + * @param w_z Weights in third direction. + */ + void store(int q_ind, Utils::Span w_x, + Utils::Span w_y, Utils::Span w_z) { + ca_fmp.push_back(q_ind); + auto it = std::back_inserter(ca_frac); + boost::copy(w_x, it); + boost::copy(w_y, it); + boost::copy(w_z, it); + } + + auto load(size_t i) const { + using Utils::make_const_span; + + assert(i < size()); + + auto const ind = ca_frac.data() + 3 * i * m_cao; + return std::make_tuple(ca_fmp[i], make_const_span(ind + 0 * m_cao, m_cao), + make_const_span(ind + 1 * m_cao, m_cao), + make_const_span(ind + 2 * m_cao, m_cao)); + } + + /** + * @brief Reset the cache. + * + * @param cao Interpolation order. + */ + void reset(int cao) { + m_cao = cao; + ca_frac.clear(); + ca_fmp.clear(); + } +}; + +#endif // ESPRESSO_P3M_INTERPOLATION_HPP From bd4061a2e2ce49a7bde1ca505029aef18d0c1ddb Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Fri, 24 Jan 2020 09:53:07 +0100 Subject: [PATCH 08/15] core: dp3m: Removed unused code --- .../p3m-common.hpp | 4 +- .../p3m-dipolar.cpp | 66 ------------------- .../electrostatics_magnetostatics/p3m.cpp | 36 ++++++---- src/python/espressomd/p3m_common.pxd | 2 - 4 files changed, 25 insertions(+), 83 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/p3m-common.hpp b/src/core/electrostatics_magnetostatics/p3m-common.hpp index 6b823ca0a95..87da9c0665f 100644 --- a/src/core/electrostatics_magnetostatics/p3m-common.hpp +++ b/src/core/electrostatics_magnetostatics/p3m-common.hpp @@ -38,6 +38,8 @@ #include "LocalBox.hpp" +#include + /** Error Codes for p3m tuning (version 2) */ enum P3M_TUNE_ERROR { /** force evaluation failed */ @@ -122,7 +124,7 @@ typedef struct { /** mesh constant. */ double a[3] = {}; /** inverse mesh constant. */ - double ai[3] = {}; + Utils::Vector3d ai = {}; /** unscaled @ref P3MParameters::alpha_L "alpha_L" for use with fast * inline functions only */ double alpha = 0.0; diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp index 92861564327..7161244cd4b 100644 --- a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp +++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp @@ -166,72 +166,6 @@ static double dp3m_average_dipolar_self_energy(double box_l, int mesh); static double dp3m_perform_aliasing_sums_dipolar_self_energy(const int n[3]); /************************************************************/ -/* functions related to the correction of the dipolar p3m-energy */ - -/* - -// Do the sum over k<>0 where k={kx,ky,kz} with kx,ky,kz INTEGERS, of -// exp(-Utils::pi()**2*k**2/alpha**2/L**2) -static double dp3m_sumi1(double alpha_L){ - int k2,kx,ky,kz,kx2,ky2,limit=60; - double suma,alpha_L2; - - alpha_L2= alpha_L* alpha_L; - - //fprintf(stderr,"alpha_L=%le\n",alpha_L); - //fprintf(stderr,"Utils::pi()=%le\n",Utils::pi()); - - - suma=0.0; - for(kx=-limit;kx<=limit;kx++){ - kx2=kx*kx; - for(ky=-limit;ky<=limit;ky++){ - ky2=ky*ky; - for(kz=-limit;kz<=limit;kz++){ - k2=kx2+ky2+kz*kz; - suma+=exp(-Utils::pi()*Utils::pi()*k2/(alpha_L*alpha_L)); - }}} - suma-=1; //It's easier to subtract the term k=0 later than put an if -inside the loops - - - //fprintf(stderr,"suma=%le\n",suma); - - - return suma; -} - -*/ - -/************************************************************/ - -/* -// Do the sum over n<>0 where n={nx*L,ny*L,nz*L} with nx,ny,nz INTEGERS, of -// exp(-alpha_iL**2*n**2) -static double dp3m_sumi2(double alpha_L){ - int n2,nx,ny,nz,nx2,ny2,limit=60; - double suma; - - - - suma=0.0; - for(nx=-limit;nx<=limit;nx++){ - nx2=nx*nx; - for(ny=-limit;ny<=limit;ny++){ - ny2=ny*ny; - for(nz=-limit;nz<=limit;nz++){ - n2=nx2+ny2+nz*nz; - suma+=exp(-alpha_L*alpha_L*n2); - }}} - suma-=1; //It's easier to subtract the term n=0 later than put an if -inside the loops - - - - return suma; -} - -*/ dp3m_data_struct::dp3m_data_struct() { params.epsilon = P3M_EPSILON_MAGNETIC; diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp index b2dc4361c89..ef66060619f 100644 --- a/src/core/electrostatics_magnetostatics/p3m.cpp +++ b/src/core/electrostatics_magnetostatics/p3m.cpp @@ -302,6 +302,8 @@ int p3m_set_eps(double eps) { template void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, + const Utils::Vector3d &ai, + p3m_local_mesh const &local_mesh, p3m_interpolation_weights *inter_weights = nullptr) { /** position shift for calc. of first assignment mesh point. */ static auto const pos_shift = std::floor((cao - 1) / 2.0) - (cao % 2) / 2.0; @@ -314,9 +316,7 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, for (int d = 0; d < 3; d++) { /* particle position in mesh coordinates */ - auto const pos = - ((real_pos[d] - p3m.local_mesh.ld_pos[d]) * p3m.params.ai[d]) - - pos_shift; + auto const pos = ((real_pos[d] - local_mesh.ld_pos[d]) * ai[d]) - pos_shift; nmp[d] = (int)pos; @@ -324,7 +324,7 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, dist[d] = (pos - nmp[d]) - 0.5; } /* 3d-array index of nearest mesh point */ - auto q_ind = Utils::get_linear_index(nmp, p3m.local_mesh.dim, + auto q_ind = Utils::get_linear_index(nmp, local_mesh.dim, Utils::MemoryOrder::ROW_MAJOR); Utils::Array w_x, w_y, w_z; @@ -350,9 +350,9 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, p3m.rs_mesh[q_ind] += cur_ca_frac_val; q_ind++; } - q_ind += p3m.local_mesh.q_2_off; + q_ind += local_mesh.q_2_off; } - q_ind += p3m.local_mesh.q_21_off; + q_ind += local_mesh.q_21_off; } } @@ -363,7 +363,8 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, template void p3m_do_charge_assign(const ParticleRange &particles) { for (auto &p : particles) { if (p.p.q != 0.0) { - p3m_do_assign_charge(p.p.q, p.r.p, &p3m.inter_weights); + p3m_do_assign_charge(p.p.q, p.r.p, p3m.params.ai, p3m.local_mesh, + &p3m.inter_weights); } } } @@ -406,25 +407,32 @@ void p3m_assign_charge(double q, const Utils::Vector3d &real_pos, p3m_interpolation_weights *inter_weights) { switch (p3m.params.cao) { case 1: - p3m_do_assign_charge<1>(q, real_pos, inter_weights); + p3m_do_assign_charge<1>(q, real_pos, p3m.params.ai, p3m.local_mesh, + inter_weights); break; case 2: - p3m_do_assign_charge<2>(q, real_pos, inter_weights); + p3m_do_assign_charge<2>(q, real_pos, p3m.params.ai, p3m.local_mesh, + inter_weights); break; case 3: - p3m_do_assign_charge<3>(q, real_pos, inter_weights); + p3m_do_assign_charge<3>(q, real_pos, p3m.params.ai, p3m.local_mesh, + inter_weights); break; case 4: - p3m_do_assign_charge<4>(q, real_pos, inter_weights); + p3m_do_assign_charge<4>(q, real_pos, p3m.params.ai, p3m.local_mesh, + inter_weights); break; case 5: - p3m_do_assign_charge<5>(q, real_pos, inter_weights); + p3m_do_assign_charge<5>(q, real_pos, p3m.params.ai, p3m.local_mesh, + inter_weights); break; case 6: - p3m_do_assign_charge<6>(q, real_pos, inter_weights); + p3m_do_assign_charge<6>(q, real_pos, p3m.params.ai, p3m.local_mesh, + inter_weights); break; case 7: - p3m_do_assign_charge<7>(q, real_pos, inter_weights); + p3m_do_assign_charge<7>(q, real_pos, p3m.params.ai, p3m.local_mesh, + inter_weights); break; } } diff --git a/src/python/espressomd/p3m_common.pxd b/src/python/espressomd/p3m_common.pxd index 8b9bb28f12f..3f7b685038c 100644 --- a/src/python/espressomd/p3m_common.pxd +++ b/src/python/espressomd/p3m_common.pxd @@ -28,8 +28,6 @@ IF P3M == 1 or DP3M == 1: double epsilon double cao_cut[3] double a[3] - double ai[3] double alpha double r_cut - int cao3 double additional_mesh[3] From e2688129884d37be42f1a1da8350c35e7ce9c97f Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Fri, 31 Jan 2020 15:52:07 +0100 Subject: [PATCH 09/15] core: p3m: Factored out charge interpolation --- .../electrostatics_magnetostatics/p3m.cpp | 64 ++---------- .../p3m_interpolation.hpp | 98 ++++++++++++++++--- 2 files changed, 96 insertions(+), 66 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp index ef66060619f..a00d9c82069 100644 --- a/src/core/electrostatics_magnetostatics/p3m.cpp +++ b/src/core/electrostatics_magnetostatics/p3m.cpp @@ -305,55 +305,15 @@ void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, const Utils::Vector3d &ai, p3m_local_mesh const &local_mesh, p3m_interpolation_weights *inter_weights = nullptr) { - /** position shift for calc. of first assignment mesh point. */ - static auto const pos_shift = std::floor((cao - 1) / 2.0) - (cao % 2) / 2.0; - - /* distance to nearest mesh point */ - double dist[3]; - - /* nearest mesh point */ - Utils::Vector3i nmp; - - for (int d = 0; d < 3; d++) { - /* particle position in mesh coordinates */ - auto const pos = ((real_pos[d] - local_mesh.ld_pos[d]) * ai[d]) - pos_shift; - - nmp[d] = (int)pos; - - /* distance to nearest mesh point */ - dist[d] = (pos - nmp[d]) - 0.5; - } - /* 3d-array index of nearest mesh point */ - auto q_ind = Utils::get_linear_index(nmp, local_mesh.dim, - Utils::MemoryOrder::ROW_MAJOR); - - Utils::Array w_x, w_y, w_z; - - for (int i = 0; i < cao; i++) { - using Utils::bspline; - - w_x[i] = bspline(i, dist[0]); - w_y[i] = bspline(i, dist[1]); - w_z[i] = bspline(i, dist[2]); - } + auto const w = + p3m_calculate_interpolation_weights(real_pos, ai, local_mesh); if (inter_weights) { - inter_weights->store(q_ind, w_x, w_y, w_z); + inter_weights->store(w); } - for (int i0 = 0; i0 < cao; i0++) { - auto const tmp0 = w_x[i0]; - for (int i1 = 0; i1 < cao; i1++) { - auto const tmp1 = tmp0 * w_y[i1]; - for (int i2 = 0; i2 < cao; i2++) { - auto const cur_ca_frac_val = q * tmp1 * w_z[i2]; - p3m.rs_mesh[q_ind] += cur_ca_frac_val; - q_ind++; - } - q_ind += local_mesh.q_2_off; - } - q_ind += local_mesh.q_21_off; - } + p3m_interpolate(local_mesh, w, + [q](int ind, double w) { p3m.rs_mesh[ind] += w * q; }); } /** Template parameterized calculation of the charge assignment to be called by @@ -455,20 +415,16 @@ static void P3M_assign_forces(double force_prefac, if (q != 0.0) { auto const pref = q * force_prefac; - Span w_x, cw_y, cw_z; - int q_ind; - std::tie(q_ind, w_x, cw_y, cw_z) = p3m.inter_weights.load(cp_cnt); + auto const w = p3m.inter_weights.load(cp_cnt); - /* We only pin the directions that are used multiple times. */ - const Vector w_y(cw_y); - const Vector w_z(cw_z); + auto q_ind = w.ind; for (int i0 = 0; i0 < cao; i0++) { - auto const fx = w_x[i0]; + auto const fx = w.w_x[i0]; for (int i1 = 0; i1 < cao; i1++) { - auto const fxy = w_y[i1] * fx; + auto const fxy = w.w_y[i1] * fx; for (int i2 = 0; i2 < cao; i2++) { - p.f.f -= pref * fxy * w_z[i2] * + p.f.f -= pref * fxy * w.w_z[i2] * Utils::Vector3d{p3m.E_mesh[0][q_ind], p3m.E_mesh[1][q_ind], p3m.E_mesh[2][q_ind]}; q_ind++; diff --git a/src/core/electrostatics_magnetostatics/p3m_interpolation.hpp b/src/core/electrostatics_magnetostatics/p3m_interpolation.hpp index 36bc160e624..5378d8ae868 100644 --- a/src/core/electrostatics_magnetostatics/p3m_interpolation.hpp +++ b/src/core/electrostatics_magnetostatics/p3m_interpolation.hpp @@ -3,6 +3,7 @@ #include #include +#include #include @@ -10,6 +11,13 @@ #include #include +template struct InterpolationWeights { + /** Linear index of the corner of the interpolation cube. */ + int ind; + /** Weights for the directions */ + Utils::Array w_x, w_y, w_z; +}; + struct p3m_interpolation_weights { size_t m_cao = 0; /** Charge fractions for mesh assignment. */ @@ -37,24 +45,31 @@ struct p3m_interpolation_weights { * @param w_y Weights in second direction. * @param w_z Weights in third direction. */ - void store(int q_ind, Utils::Span w_x, - Utils::Span w_y, Utils::Span w_z) { - ca_fmp.push_back(q_ind); + template void store(const InterpolationWeights &w) { + assert(cao == m_cao); + + ca_fmp.push_back(w.ind); auto it = std::back_inserter(ca_frac); - boost::copy(w_x, it); - boost::copy(w_y, it); - boost::copy(w_z, it); + boost::copy(w.w_x, it); + boost::copy(w.w_y, it); + boost::copy(w.w_z, it); } - auto load(size_t i) const { - using Utils::make_const_span; + template InterpolationWeights load(size_t i) const { + assert(cao == m_cao); + using Utils::make_const_span; assert(i < size()); - auto const ind = ca_frac.data() + 3 * i * m_cao; - return std::make_tuple(ca_fmp[i], make_const_span(ind + 0 * m_cao, m_cao), - make_const_span(ind + 1 * m_cao, m_cao), - make_const_span(ind + 2 * m_cao, m_cao)); + InterpolationWeights ret; + ret.ind = ca_fmp[i]; + + auto const offset = ca_frac.data() + 3 * i * m_cao; + boost::copy(make_const_span(offset + 0 * m_cao, m_cao), ret.w_x.begin()); + boost::copy(make_const_span(offset + 1 * m_cao, m_cao), ret.w_y.begin()); + boost::copy(make_const_span(offset + 2 * m_cao, m_cao), ret.w_z.begin()); + + return ret; } /** @@ -69,4 +84,63 @@ struct p3m_interpolation_weights { } }; +template +InterpolationWeights +p3m_calculate_interpolation_weights(const Utils::Vector3d &real_pos, + const Utils::Vector3d &ai, + p3m_local_mesh const &local_mesh) { + /** position shift for calc. of first assignment mesh point. */ + static auto const pos_shift = std::floor((cao - 1) / 2.0) - (cao % 2) / 2.0; + + /* distance to nearest mesh point */ + double dist[3]; + + /* nearest mesh point */ + Utils::Vector3i nmp; + + for (int d = 0; d < 3; d++) { + /* particle position in mesh coordinates */ + auto const pos = ((real_pos[d] - local_mesh.ld_pos[d]) * ai[d]) - pos_shift; + + nmp[d] = (int)pos; + + /* distance to nearest mesh point */ + dist[d] = (pos - nmp[d]) - 0.5; + } + + InterpolationWeights ret; + + /* 3d-array index of nearest mesh point */ + ret.ind = Utils::get_linear_index(nmp, local_mesh.dim, + Utils::MemoryOrder::ROW_MAJOR); + for (int i = 0; i < cao; i++) { + using Utils::bspline; + + ret.w_x[i] = bspline(i, dist[0]); + ret.w_y[i] = bspline(i, dist[1]); + ret.w_z[i] = bspline(i, dist[2]); + } + + return ret; +} + +template +void p3m_interpolate(p3m_local_mesh const &local_mesh, + InterpolationWeights const &w, Kernel kernel) { + auto q_ind = w.ind; + for (int i0 = 0; i0 < cao; i0++) { + auto const tmp0 = w.w_x[i0]; + for (int i1 = 0; i1 < cao; i1++) { + auto const tmp1 = tmp0 * w.w_y[i1]; + for (int i2 = 0; i2 < cao; i2++) { + kernel(q_ind, tmp1 * w.w_z[i2]); + + q_ind++; + } + q_ind += local_mesh.q_2_off; + } + q_ind += local_mesh.q_21_off; + } +} + #endif // ESPRESSO_P3M_INTERPOLATION_HPP From 1b49e39a2f5bbff7edd704d554eda24142444400 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Fri, 31 Jan 2020 16:21:17 +0100 Subject: [PATCH 10/15] core: p3m: Factored out force interpolation --- .../electrostatics_magnetostatics/p3m.cpp | 25 ++++++------------- 1 file changed, 7 insertions(+), 18 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp index a00d9c82069..9d243bbb4b6 100644 --- a/src/core/electrostatics_magnetostatics/p3m.cpp +++ b/src/core/electrostatics_magnetostatics/p3m.cpp @@ -414,26 +414,15 @@ static void P3M_assign_forces(double force_prefac, auto const q = p.p.q; if (q != 0.0) { auto const pref = q * force_prefac; + auto const w = p3m.inter_weights.load(cp_cnt++); - auto const w = p3m.inter_weights.load(cp_cnt); - - auto q_ind = w.ind; + Utils::Vector3d E{}; + p3m_interpolate(p3m.local_mesh, w, [&E](int ind, double w) { + E += w * Utils::Vector3d{p3m.E_mesh[0][ind], p3m.E_mesh[1][ind], + p3m.E_mesh[2][ind]}; + }); - for (int i0 = 0; i0 < cao; i0++) { - auto const fx = w.w_x[i0]; - for (int i1 = 0; i1 < cao; i1++) { - auto const fxy = w.w_y[i1] * fx; - for (int i2 = 0; i2 < cao; i2++) { - p.f.f -= pref * fxy * w.w_z[i2] * - Utils::Vector3d{p3m.E_mesh[0][q_ind], p3m.E_mesh[1][q_ind], - p3m.E_mesh[2][q_ind]}; - q_ind++; - } - q_ind += p3m.local_mesh.q_2_off; - } - q_ind += p3m.local_mesh.q_21_off; - } - cp_cnt++; + p.f.f -= pref * E; } } } From 0e7cab4ccc14197429bc873571d7b1d09314e2e5 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Fri, 31 Jan 2020 20:22:30 +0100 Subject: [PATCH 11/15] core: dp3m: Modernize --- .../p3m-dipolar.cpp | 58 +++++++++++-------- .../p3m-dipolar.hpp | 16 ----- 2 files changed, 34 insertions(+), 40 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp index 7161244cd4b..fb5465f269f 100644 --- a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp +++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp @@ -408,27 +408,18 @@ int dp3m_set_eps(double eps) { return ES_OK; } -void dp3m_dipole_assign(const ParticleRange &particles) { - /* magnetic particle counter, dipole fraction counter */ - int cp_cnt = 0; - - /* prepare local FFT mesh */ - for (auto &i : dp3m.rs_mesh_dip) - for (int j = 0; j < dp3m.local_mesh.size; j++) - i[j] = 0.0; - - for (auto const &p : particles) { - if (p.p.dipm != 0.0) { - dp3m_assign_dipole(p.r.p.data(), p.p.dipm, p.calc_dip().data(), cp_cnt); - cp_cnt++; - } - } - - dp3m_shrink_wrap_dipole_grid(cp_cnt); -} - -void dp3m_assign_dipole(double const real_pos[3], double mu, - double const dip[3], int cp_cnt) { +/** Assign a single dipole into the current dipole grid. + * + * @param[in] real_pos %Particle position in real space + * @param[in] mu %Particle magnetic dipole magnitude + * @param[in] dip %Particle magnetic dipole vector + * @param[in] cp_cnt The running index, which may be smaller than 0, in + * which case the dipole is assumed to be virtual and + * is not stored in the @ref dp3m_data_struct::ca_frac + * "ca_frac" arrays + */ +static void dp3m_assign_dipole(Utils::Vector3d const &real_pos, double mu, + Utils::Vector3d const &dip, int cp_cnt) { /* we do not really want to export these, but this function should be inlined */ int d, i0, i1, i2; @@ -465,11 +456,11 @@ void dp3m_assign_dipole(double const real_pos[3], double mu, if (cp_cnt >= 0) dp3m.ca_fmp[cp_cnt] = q_ind; - for (i0 = 0; i0 < dp3m.params.cao; i0++) { + for (int i0 = 0; i0 < dp3m.params.cao; i0++) { tmp0 = Utils::bspline(i0, dist[0], dp3m.params.cao); - for (i1 = 0; i1 < dp3m.params.cao; i1++) { + for (int i1 = 0; i1 < dp3m.params.cao; i1++) { tmp1 = tmp0 * Utils::bspline(i1, dist[1], dp3m.params.cao); - for (i2 = 0; i2 < dp3m.params.cao; i2++) { + for (int i2 = 0; i2 < dp3m.params.cao; i2++) { cur_ca_frac_val = tmp1 * Utils::bspline(i2, dist[2], dp3m.params.cao); if (cp_cnt >= 0) *(cur_ca_frac++) = cur_ca_frac_val; @@ -486,6 +477,25 @@ void dp3m_assign_dipole(double const real_pos[3], double mu, } } +void dp3m_dipole_assign(const ParticleRange &particles) { + /* magnetic particle counter, dipole fraction counter */ + int cp_cnt = 0; + + /* prepare local FFT mesh */ + for (auto &i : dp3m.rs_mesh_dip) + for (int j = 0; j < dp3m.local_mesh.size; j++) + i[j] = 0.0; + + for (auto const &p : particles) { + if (p.p.dipm != 0.0) { + dp3m_assign_dipole(p.r.p, p.p.dipm, p.calc_dip(), cp_cnt); + cp_cnt++; + } + } + + dp3m_shrink_wrap_dipole_grid(cp_cnt); +} + void dp3m_shrink_wrap_dipole_grid(int n_dipoles) { if (n_dipoles < dp3m.ca_num) dp3m_realloc_ca_fields(n_dipoles); diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp index 9165e7285fe..bd4b20ec79d 100644 --- a/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp +++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp @@ -65,9 +65,6 @@ struct dp3m_data_struct { /** Sum of square of magnetic dipoles (only on master node). */ double sum_mu2; - /** interpolation of the charge assignment function. */ - std::vector> int_caf{7}; - /** position shift for calc. of first assignment mesh point. */ double pos_shift; /** help variable for calculation of aliasing sums */ @@ -192,19 +189,6 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag, */ void dp3m_count_magnetic_particles(); -/** Assign a single dipole into the current dipole grid. - * - * @param[in] real_pos %Particle position in real space - * @param[in] mu %Particle magnetic dipole magnitude - * @param[in] dip %Particle magnetic dipole vector - * @param[in] cp_cnt The running index, which may be smaller than 0, in - * which case the dipole is assumed to be virtual and - * is not stored in the @ref dp3m_data_struct::ca_frac - * "ca_frac" arrays - */ -void dp3m_assign_dipole(double const real_pos[3], double mu, - double const dip[3], int cp_cnt); - /** Shrink wrap the dipoles grid */ void dp3m_shrink_wrap_dipole_grid(int n_dipoles); From a41eb3224497eb00506e1fd2c9e66457802ef202 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Mon, 3 Feb 2020 10:30:24 +0100 Subject: [PATCH 12/15] core: p3m: Docstrings --- .../p3m-dipolar.cpp | 211 +++++------------- .../p3m-dipolar.hpp | 11 +- .../electrostatics_magnetostatics/p3m.cpp | 4 +- .../electrostatics_magnetostatics/p3m.hpp | 12 +- .../p3m_interpolation.hpp | 59 ++++- 5 files changed, 106 insertions(+), 191 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp index fb5465f269f..cfa10667e88 100644 --- a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp +++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp @@ -82,9 +82,6 @@ dp3m_data_struct dp3m; */ static void dp3m_init_a_ai_cao_cut(); -/** realloc charge assignment fields. */ -static void dp3m_realloc_ca_fields(int newsize); - /** Checks for correctness for magnetic dipoles in P3M of the cao_cut, * necessary when the box length changes */ @@ -176,8 +173,6 @@ dp3m_data_struct::dp3m_data_struct() { sum_mu2 = 0.0; pos_shift = 0.0; - - ca_num = 0; ks_pnum = 0; energy_correction = 0.0; @@ -195,8 +190,6 @@ void dp3m_deactivate() { } void dp3m_init() { - int n; - if (dipole.prefactor <= 0.0) { // dipolar prefactor is zero: magnetostatics switched off dp3m.params.r_cut = 0.0; @@ -211,14 +204,6 @@ void dp3m_init() { * and the cutoff for charge assignment dp3m.params.cao_cut */ dp3m_init_a_ai_cao_cut(); - /* initialize ca fields to size CA_INCREMENT: dp3m.ca_frac and dp3m.ca_fmp - */ - dp3m.ca_num = 0; - if (dp3m.ca_num < CA_INCREMENT) { - dp3m.ca_num = 0; - dp3m_realloc_ca_fields(CA_INCREMENT); - } - p3m_calc_local_ca_mesh(dp3m.local_mesh, dp3m.params, local_geo, skin); dp3m.sm.resize(comm_cart, dp3m.local_mesh); @@ -226,9 +211,6 @@ void dp3m_init() { /* fix box length dependent constants */ dp3m_scaleby_box_l(); - dp3m.pos_shift = - std::floor((dp3m.params.cao - 1) / 2.0) - (dp3m.params.cao % 2) / 2.0; - int ca_mesh_size = fft_init(dp3m.local_mesh.dim, dp3m.local_mesh.margin, dp3m.params.mesh, dp3m.params.mesh_off, &dp3m.ks_pnum, dp3m.fft, node_grid, comm_cart); @@ -411,75 +393,24 @@ int dp3m_set_eps(double eps) { /** Assign a single dipole into the current dipole grid. * * @param[in] real_pos %Particle position in real space - * @param[in] mu %Particle magnetic dipole magnitude * @param[in] dip %Particle magnetic dipole vector - * @param[in] cp_cnt The running index, which may be smaller than 0, in - * which case the dipole is assumed to be virtual and - * is not stored in the @ref dp3m_data_struct::ca_frac - * "ca_frac" arrays */ -static void dp3m_assign_dipole(Utils::Vector3d const &real_pos, double mu, - Utils::Vector3d const &dip, int cp_cnt) { - /* we do not really want to export these, but this function should be inlined - */ - int d, i0, i1, i2; - double tmp0, tmp1; - /* position of a particle in local mesh units */ - double pos; - /* 1d-index of nearest mesh point */ - int nmp; - /* distance to nearest mesh point */ - double dist[3]; - /* index, index jumps for dp3m.rs_mesh array */ - int q_ind = 0; - double cur_ca_frac_val, *cur_ca_frac; - - // make sure we have enough space - if (cp_cnt >= dp3m.ca_num) - dp3m_realloc_ca_fields(cp_cnt + 1); - // do it here, since p3m_realloc_ca_fields may change the address of - // dp3m.ca_frac - cur_ca_frac = dp3m.ca_frac.data() + dp3m.params.cao3 * cp_cnt; - - for (d = 0; d < 3; d++) { - /* particle position in mesh coordinates */ - pos = ((real_pos[d] - dp3m.local_mesh.ld_pos[d]) * dp3m.params.ai[d]) - - dp3m.pos_shift; - /* nearest mesh point */ - nmp = (int)pos; - /* distance to nearest mesh point */ - dist[d] = (pos - nmp) - 0.5; - /* 3d-array index of nearest mesh point */ - q_ind = (d == 0) ? nmp : nmp + dp3m.local_mesh.dim[d] * q_ind; - } - - if (cp_cnt >= 0) - dp3m.ca_fmp[cp_cnt] = q_ind; - - for (int i0 = 0; i0 < dp3m.params.cao; i0++) { - tmp0 = Utils::bspline(i0, dist[0], dp3m.params.cao); - for (int i1 = 0; i1 < dp3m.params.cao; i1++) { - tmp1 = tmp0 * Utils::bspline(i1, dist[1], dp3m.params.cao); - for (int i2 = 0; i2 < dp3m.params.cao; i2++) { - cur_ca_frac_val = tmp1 * Utils::bspline(i2, dist[2], dp3m.params.cao); - if (cp_cnt >= 0) - *(cur_ca_frac++) = cur_ca_frac_val; - if (mu != 0.0) { - dp3m.rs_mesh_dip[0][q_ind] += dip[0] * cur_ca_frac_val; - dp3m.rs_mesh_dip[1][q_ind] += dip[1] * cur_ca_frac_val; - dp3m.rs_mesh_dip[2][q_ind] += dip[2] * cur_ca_frac_val; - } - q_ind++; - } - q_ind += dp3m.local_mesh.q_2_off; - } - q_ind += dp3m.local_mesh.q_21_off; - } +template +static void dp3m_assign_dipole(Utils::Vector3d const &real_pos, + Utils::Vector3d const &dip) { + auto const weights = p3m_calculate_interpolation_weights( + real_pos, dp3m.params.ai, dp3m.local_mesh); + p3m_interpolate(dp3m.local_mesh, weights, [&dip](int ind, double w) { + dp3m.rs_mesh_dip[0][ind] += w * dip[0]; + dp3m.rs_mesh_dip[1][ind] += w * dip[1]; + dp3m.rs_mesh_dip[2][ind] += w * dip[2]; + }); + + dp3m.inter_weights.store(weights); } void dp3m_dipole_assign(const ParticleRange &particles) { - /* magnetic particle counter, dipole fraction counter */ - int cp_cnt = 0; + dp3m.inter_weights.reset(dp3m.params.cao); /* prepare local FFT mesh */ for (auto &i : dp3m.rs_mesh_dip) @@ -488,74 +419,52 @@ void dp3m_dipole_assign(const ParticleRange &particles) { for (auto const &p : particles) { if (p.p.dipm != 0.0) { - dp3m_assign_dipole(p.r.p, p.p.dipm, p.calc_dip(), cp_cnt); - cp_cnt++; + switch (dp3m.params.cao) { + case 1: + dp3m_assign_dipole<1>(p.r.p, p.calc_dip()); + break; + case 2: + dp3m_assign_dipole<2>(p.r.p, p.calc_dip()); + break; + case 3: + dp3m_assign_dipole<3>(p.r.p, p.calc_dip()); + break; + case 4: + dp3m_assign_dipole<4>(p.r.p, p.calc_dip()); + break; + case 5: + dp3m_assign_dipole<5>(p.r.p, p.calc_dip()); + break; + case 6: + dp3m_assign_dipole<6>(p.r.p, p.calc_dip()); + break; + case 7: + dp3m_assign_dipole<7>(p.r.p, p.calc_dip()); + break; + default: + throw std::runtime_error("Unknown order."); + } } } - - dp3m_shrink_wrap_dipole_grid(cp_cnt); -} - -void dp3m_shrink_wrap_dipole_grid(int n_dipoles) { - if (n_dipoles < dp3m.ca_num) - dp3m_realloc_ca_fields(n_dipoles); } #ifdef ROTATION /** Assign the torques obtained from k-space */ +template static void P3M_assign_torques(double prefac, int d_rs, const ParticleRange &particles) { - /* particle counter, charge fraction counter */ - int cp_cnt = 0, cf_cnt = 0; - /* index, index jumps for dp3m.rs_mesh array */ - int q_ind; - auto const q_m_off = (dp3m.local_mesh.dim[2] - dp3m.params.cao); - auto const q_s_off = - dp3m.local_mesh.dim[2] * (dp3m.local_mesh.dim[1] - dp3m.params.cao); - + /* particle counter */ + int cp_cnt = 0; for (auto &p : particles) { - if ((p.p.dipm) != 0.0) { - const Utils::Vector3d dip = p.calc_dip(); - q_ind = dp3m.ca_fmp[cp_cnt]; - for (int i0 = 0; i0 < dp3m.params.cao; i0++) { - for (int i1 = 0; i1 < dp3m.params.cao; i1++) { - for (int i2 = 0; i2 < dp3m.params.cao; i2++) { - /* - The following line would fill the torque with the k-space electric - field (without the self-field term) [notice the minus sign!]: - p.f.torque[d_rs] -= - prefac*dp3m.ca_frac[cf_cnt]*dp3m.rs_mesh[q_ind]; - Since the torque is the dipole moment cross-product with E, we - have: - */ - switch (d_rs) { - case 0: // E_x - p.f.torque[1] -= - dip[2] * prefac * dp3m.ca_frac[cf_cnt] * dp3m.rs_mesh[q_ind]; - p.f.torque[2] += - dip[1] * prefac * dp3m.ca_frac[cf_cnt] * dp3m.rs_mesh[q_ind]; - break; - case 1: // E_y - p.f.torque[0] += - dip[2] * prefac * dp3m.ca_frac[cf_cnt] * dp3m.rs_mesh[q_ind]; - p.f.torque[2] -= - dip[0] * prefac * dp3m.ca_frac[cf_cnt] * dp3m.rs_mesh[q_ind]; - break; - case 2: // E_z - p.f.torque[0] -= - dip[1] * prefac * dp3m.ca_frac[cf_cnt] * dp3m.rs_mesh[q_ind]; - p.f.torque[1] += - dip[0] * prefac * dp3m.ca_frac[cf_cnt] * dp3m.rs_mesh[q_ind]; - } - q_ind++; - cf_cnt++; - } - q_ind += q_m_off; - } - q_ind += q_s_off; - } - cp_cnt++; - } + auto const w = dp3m.inter_weights.load(cp_cnt++); + + Utils::Vector3d E{}; + + p3m_interpolate(dp3m.local_mesh, w, [&E, d_rs](int ind, double w) { + E[d_rs] += w * dp3m.rs_mesh[ind]; + }); + + p.f.torque += vector_product(p.calc_dip(), prefac * E); } } #endif @@ -935,20 +844,6 @@ double calc_surface_term(bool force_flag, bool energy_flag, /*****************************************************************************/ -void dp3m_realloc_ca_fields(int newsize) { - newsize = ((newsize + CA_INCREMENT - 1) / CA_INCREMENT) * CA_INCREMENT; - if (newsize == dp3m.ca_num) - return; - if (newsize < CA_INCREMENT) - newsize = CA_INCREMENT; - - dp3m.ca_num = newsize; - dp3m.ca_frac.resize(dp3m.params.cao3 * dp3m.ca_num); - dp3m.ca_fmp.resize(dp3m.ca_num); -} - -/*****************************************************************************/ - void dp3m_calc_meshift() { int i; double dmesh; @@ -1848,11 +1743,7 @@ bool dp3m_sanity_checks(const Utils::Vector3i &grid) { runtimeErrorMsg() << "dipolar P3M requires periodicity 1 1 1"; ret = true; } - /* - if (n_nodes != 1) { - runtimeErrorMsg() <<"dipolar P3M does not run in parallel"; - ret = true; - } */ + if (cell_structure.type != CELL_STRUCTURE_DOMDEC) { runtimeErrorMsg() << "dipolar P3M at present requires the domain " "decomposition cell system"; diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp index bd4b20ec79d..6e7b565eb74 100644 --- a/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp +++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp @@ -40,6 +40,7 @@ #include "electrostatics_magnetostatics/dipole.hpp" #include "fft.hpp" #include "p3m-common.hpp" +#include "p3m_interpolation.hpp" #include "p3m_send_mesh.hpp" #include @@ -78,13 +79,8 @@ struct dp3m_data_struct { /** Energy optimised influence function (k-space) */ std::vector g_energy; - /** number of charged particles on the node. */ - int ca_num; + p3m_interpolation_cache inter_weights; - /** Charge fractions for mesh assignment. */ - std::vector ca_frac; - /** index of first mesh point for charge assignment. */ - std::vector ca_fmp; /** number of permutations in k_space */ int ks_pnum; @@ -189,9 +185,6 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag, */ void dp3m_count_magnetic_particles(); -/** Shrink wrap the dipoles grid */ -void dp3m_shrink_wrap_dipole_grid(int n_dipoles); - /** Calculate real space contribution of p3m dipolar pair forces and torques. * If NPT is compiled in, it returns the energy, which is needed for NPT. */ diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp index 9d243bbb4b6..de66cec3c41 100644 --- a/src/core/electrostatics_magnetostatics/p3m.cpp +++ b/src/core/electrostatics_magnetostatics/p3m.cpp @@ -304,7 +304,7 @@ template void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos, const Utils::Vector3d &ai, p3m_local_mesh const &local_mesh, - p3m_interpolation_weights *inter_weights = nullptr) { + p3m_interpolation_cache *inter_weights = nullptr) { auto const w = p3m_calculate_interpolation_weights(real_pos, ai, local_mesh); @@ -364,7 +364,7 @@ void p3m_charge_assign(const ParticleRange &particles) { /* Template wrapper for p3m_do_assign_charge() */ void p3m_assign_charge(double q, const Utils::Vector3d &real_pos, - p3m_interpolation_weights *inter_weights) { + p3m_interpolation_cache *inter_weights) { switch (p3m.params.cao) { case 1: p3m_do_assign_charge<1>(q, real_pos, p3m.params.ai, p3m.local_mesh, diff --git a/src/core/electrostatics_magnetostatics/p3m.hpp b/src/core/electrostatics_magnetostatics/p3m.hpp index e06a8de9c66..c375555b379 100644 --- a/src/core/electrostatics_magnetostatics/p3m.hpp +++ b/src/core/electrostatics_magnetostatics/p3m.hpp @@ -82,7 +82,7 @@ struct p3m_data_struct { /** Energy optimised influence function (k-space) */ std::vector g_energy; - p3m_interpolation_weights inter_weights; + p3m_interpolation_cache inter_weights; /** number of permutations in k_space */ int ks_pnum; @@ -164,9 +164,6 @@ bool p3m_sanity_checks(); void p3m_count_charged_particles(); /** Assign the physical charges using the tabulated charge assignment function. - * The charge fractions are buffered - * in @ref p3m_data_struct::ca_fmp "ca_fmp" and @ref p3m_data_struct::ca_frac - * "ca_frac". */ void p3m_charge_assign(const ParticleRange &particles); @@ -174,12 +171,11 @@ void p3m_charge_assign(const ParticleRange &particles); * * @param[in] q %Particle charge * @param[in] real_pos %Particle position in real space - * @param[in] inter_weights The running index, which may be smaller than 0, - * in which case the charge is assumed to be virtual and is not stored in the - * @ref p3m_data_struct::ca_frac "ca_frac" arrays + * @param[in] inter_weights Cached interpolation weights to be used, can be + * null. */ void p3m_assign_charge(double q, const Utils::Vector3d &real_pos, - p3m_interpolation_weights *inter_weights); + p3m_interpolation_cache *inter_weights); /** Calculate real space contribution of Coulomb pair forces. */ inline void p3m_add_pair_force(double q1q2, Utils::Vector3d const &d, diff --git a/src/core/electrostatics_magnetostatics/p3m_interpolation.hpp b/src/core/electrostatics_magnetostatics/p3m_interpolation.hpp index 5378d8ae868..00b255e4cb8 100644 --- a/src/core/electrostatics_magnetostatics/p3m_interpolation.hpp +++ b/src/core/electrostatics_magnetostatics/p3m_interpolation.hpp @@ -11,6 +11,13 @@ #include #include +/** + * @brief Interpolation weights for one point. + * + * Interpolation weights and grid offset for one point. + * + * @tparam cao Interpolation order. + */ template struct InterpolationWeights { /** Linear index of the corner of the interpolation cube. */ int ind; @@ -18,13 +25,20 @@ template struct InterpolationWeights { Utils::Array w_x, w_y, w_z; }; -struct p3m_interpolation_weights { +/** + * @brief Cache for interpolation weights. + * + * This is a storage container for interpolation weights of + * type InterpolationWeights. + */ +class p3m_interpolation_cache { size_t m_cao = 0; /** Charge fractions for mesh assignment. */ std::vector ca_frac; /** index of first mesh point for charge assignment. */ std::vector ca_fmp; +public: /** * @brief Number of points in the cache. * @return Number of points currently in the cache. @@ -40,10 +54,9 @@ struct p3m_interpolation_weights { /** * @brief Push back weights for one point. * - * @param q_ind Mesh index. - * @param w_x Weights in first direction. - * @param w_y Weights in second direction. - * @param w_z Weights in third direction. + * @tparam cao Interpolation order has to match the order + * set at last call to @ref p3m_interpolation_cache::reset. + * @param w Interpolation weights to store. */ template void store(const InterpolationWeights &w) { assert(cao == m_cao); @@ -55,6 +68,17 @@ struct p3m_interpolation_weights { boost::copy(w.w_z, it); } + /** + * @brief Load entry from the cache. + * + * This loads an entry at an index from the cache, + * the entries are indexed by the order they were stored. + * + * @tparam cao Interpolation order has to match the order + * set at last call to @ref p3m_interpolation_cache::reset. + * @param i Index of the entry to load. + * @return i-it interpolation weights. + */ template InterpolationWeights load(size_t i) const { assert(cao == m_cao); @@ -86,7 +110,7 @@ struct p3m_interpolation_weights { template InterpolationWeights -p3m_calculate_interpolation_weights(const Utils::Vector3d &real_pos, +p3m_calculate_interpolation_weights(const Utils::Vector3d &position, const Utils::Vector3d &ai, p3m_local_mesh const &local_mesh) { /** position shift for calc. of first assignment mesh point. */ @@ -100,7 +124,7 @@ p3m_calculate_interpolation_weights(const Utils::Vector3d &real_pos, for (int d = 0; d < 3; d++) { /* particle position in mesh coordinates */ - auto const pos = ((real_pos[d] - local_mesh.ld_pos[d]) * ai[d]) - pos_shift; + auto const pos = ((position[d] - local_mesh.ld_pos[d]) * ai[d]) - pos_shift; nmp[d] = (int)pos; @@ -124,16 +148,27 @@ p3m_calculate_interpolation_weights(const Utils::Vector3d &real_pos, return ret; } +/** + * @brief P3M grid interpolation. + * + * This runs an kernel for every interpolation point + * in a set of interpolation weights with the linear + * grid index and the weight of the point as arguments. + * + * @param local_mesh Mesh info. + * @param weights Set of weights + * @param kernel The kernel to run. + */ template void p3m_interpolate(p3m_local_mesh const &local_mesh, - InterpolationWeights const &w, Kernel kernel) { - auto q_ind = w.ind; + InterpolationWeights const &weights, Kernel kernel) { + auto q_ind = weights.ind; for (int i0 = 0; i0 < cao; i0++) { - auto const tmp0 = w.w_x[i0]; + auto const tmp0 = weights.w_x[i0]; for (int i1 = 0; i1 < cao; i1++) { - auto const tmp1 = tmp0 * w.w_y[i1]; + auto const tmp1 = tmp0 * weights.w_y[i1]; for (int i2 = 0; i2 < cao; i2++) { - kernel(q_ind, tmp1 * w.w_z[i2]); + kernel(q_ind, tmp1 * weights.w_z[i2]); q_ind++; } From 054f1fbd6e49a6b067766f2113f5e0fe0dec9360 Mon Sep 17 00:00:00 2001 From: Florian Weik Date: Mon, 20 Apr 2020 14:17:30 +0200 Subject: [PATCH 13/15] utils: integral_parameter + test --- .../include/utils/integral_parameter.hpp | 63 +++++++++++++++++++ src/utils/tests/CMakeLists.txt | 1 + src/utils/tests/integral_parameter_test.cpp | 46 ++++++++++++++ 3 files changed, 110 insertions(+) create mode 100644 src/utils/include/utils/integral_parameter.hpp create mode 100644 src/utils/tests/integral_parameter_test.cpp diff --git a/src/utils/include/utils/integral_parameter.hpp b/src/utils/include/utils/integral_parameter.hpp new file mode 100644 index 00000000000..1df9595b0c2 --- /dev/null +++ b/src/utils/include/utils/integral_parameter.hpp @@ -0,0 +1,63 @@ +/* + * Copyright (C) 2020 The ESPResSo project + * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + * Max-Planck-Institute for Polymer Research, Theory Group + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +#ifndef ESPRESSO_INTEGRAL_PARAMETER_HPP +#define ESPRESSO_INTEGRAL_PARAMETER_HPP + +#include + +namespace Utils { +namespace detail { +template