Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

core: replaced utils::strcat_alloc with printf #4069

Merged
merged 5 commits into from
Dec 23, 2020
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 9 additions & 14 deletions src/core/electrostatics_magnetostatics/mmm1d.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,6 @@
#include <utils/Vector.hpp>
#include <utils/constants.hpp>
#include <utils/math/sqr.hpp>
#include <utils/strcat_alloc.hpp>

#include <algorithm>
#include <cmath>
Expand All @@ -53,8 +52,6 @@
#include <tuple>
#include <vector>

using Utils::strcat_alloc;

/** How many trial calculations in @ref mmm1d_tune */
#define TEST_INTEGRATIONS 1000

Expand Down Expand Up @@ -339,10 +336,9 @@ double mmm1d_coulomb_pair_energy(double const chpref, Utils::Vector3d const &d,
return chpref * E;
}

int mmm1d_tune(char **log) {
int mmm1d_tune(bool verbose) {
if (MMM1D_sanity_checks())
return ES_ERROR;
char buffer[32 + 2 * ES_DOUBLE_SPACE + ES_INTEGER_SPACE];
double min_time = std::numeric_limits<double>::infinity();
double min_rad = -1;
auto const maxrad = box_geo.length()[2];
Expand Down Expand Up @@ -372,8 +368,9 @@ int mmm1d_tune(char **log) {
if (int_time < 0)
return ES_ERROR;

sprintf(buffer, "r= %f t= %f ms\n", switch_radius, int_time);
*log = strcat_alloc(*log, buffer);
if (verbose) {
std::printf("r= %f t= %f ms\n", switch_radius, int_time);
}

if (int_time < min_time) {
min_time = int_time;
Expand All @@ -385,13 +382,11 @@ int mmm1d_tune(char **log) {
}
switch_radius = min_rad;
mmm1d_params.far_switch_radius_2 = Utils::sqr(switch_radius);
} else {
if (mmm1d_params.far_switch_radius_2 <=
Utils::sqr(bessel_radii[MAXIMAL_B_CUT - 1])) {
// this switching radius is too small for our Bessel series
*log = strcat_alloc(*log, "could not find reasonable bessel cutoff");
return ES_ERROR;
}
} else if (mmm1d_params.far_switch_radius_2 <=
Utils::sqr(bessel_radii[MAXIMAL_B_CUT - 1])) {
// this switching radius is too small for our Bessel series
runtimeErrorMsg() << "could not find reasonable bessel cutoff";
return ES_ERROR;
}

coulomb.method = COULOMB_MMM1D;
Expand Down
4 changes: 2 additions & 2 deletions src/core/electrostatics_magnetostatics/mmm1d.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -76,11 +76,11 @@ double mmm1d_coulomb_pair_energy(double q1q2, Utils::Vector3d const &d,
* @ref MMM1D_struct::bessel_cutoff "Bessel cutoff". Call this only
* on the master node.
*
* @param log contains information about the tuning (tried values and errors)
* @param verbose output information about the tuning (tried values and errors)
* @retval ES_OK
* @retval ES_ERROR
*/
int mmm1d_tune(char **log);
int mmm1d_tune(bool verbose);

#endif
#endif
117 changes: 63 additions & 54 deletions src/core/electrostatics_magnetostatics/p3m-dipolar.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,6 @@
#include <utils/math/int_pow.hpp>
#include <utils/math/sinc.hpp>
#include <utils/math/sqr.hpp>
#include <utils/strcat_alloc.hpp>

#include <boost/mpi/collectives/all_reduce.hpp>
#include <boost/mpi/collectives/reduce.hpp>
Expand All @@ -73,8 +72,6 @@
#include <functional>
#include <vector>

using Utils::strcat_alloc;

/************************************************
* DEFINES
************************************************/
Expand Down Expand Up @@ -831,26 +828,25 @@ static double dp3m_mcr_time(int mesh, int cao, double r_cut_iL,
*
* The @p _r_cut_iL is determined via a simple bisection.
*
* @param[out] log log output
* @param[in] mesh @copybrief P3MParameters::mesh
* @param[in] cao @copybrief P3MParameters::cao
* @param[in] r_cut_iL_min lower bound for @p _r_cut_iL
* @param[in] r_cut_iL_max upper bound for @p _r_cut_iL
* @param[out] _r_cut_iL @copybrief P3MParameters::r_cut_iL
* @param[out] _alpha_L @copybrief P3MParameters::alpha_L
* @param[out] _accuracy @copybrief P3MParameters::accuracy
* @param[in] verbose printf output
*
* @returns The integration time in case of success, otherwise
* -@ref P3M_TUNE_FAIL, -@ref P3M_TUNE_ACCURACY_TOO_LARGE,
* -@ref P3M_TUNE_CAO_TOO_LARGE, -@ref P3M_TUNE_ELCTEST, or
* -@ref P3M_TUNE_CUTOFF_TOO_LARGE
*/
static double dp3m_mc_time(char **log, int mesh, int cao, double r_cut_iL_min,
static double dp3m_mc_time(int mesh, int cao, double r_cut_iL_min,
double r_cut_iL_max, double *_r_cut_iL,
double *_alpha_L, double *_accuracy) {
double *_alpha_L, double *_accuracy, bool verbose) {
double r_cut_iL;
double rs_err, ks_err;
char b[3 * ES_INTEGER_SPACE + 3 * ES_DOUBLE_SPACE + 128];

/* initial checks. */
auto const mesh_size = box_geo.length()[0] / static_cast<double>(mesh);
Expand All @@ -861,8 +857,9 @@ static double dp3m_mc_time(char **log, int mesh, int cao, double r_cut_iL_min,

if (cao >= mesh || k_cut >= std::min(min_box_l, min_local_box_l) - skin) {
/* print result */
sprintf(b, "%-4d %-3d cao too large for this mesh\n", mesh, cao);
*log = strcat_alloc(*log, b);
if (verbose) {
std::printf("%-4d %-3d cao too large for this mesh\n", mesh, cao);
}
return -P3M_TUNE_CAO_TOO_LARGE;
}

Expand All @@ -876,9 +873,11 @@ static double dp3m_mc_time(char **log, int mesh, int cao, double r_cut_iL_min,
return *_accuracy;
if (*_accuracy > dp3m.params.accuracy) {
/* print result */
sprintf(b, "%-4d %-3d %.5e %.5e %.5e %.3e %.3e accuracy not achieved\n",
mesh, cao, r_cut_iL_max, *_alpha_L, *_accuracy, rs_err, ks_err);
*log = strcat_alloc(*log, b);
if (verbose) {
std::printf("%-4d %-3d %.5e %.5e %.5e %.3e %.3e accuracy not achieved\n",
mesh, cao, r_cut_iL_max, *_alpha_L, *_accuracy, rs_err,
ks_err);
}
return -P3M_TUNE_ACCURACY_TOO_LARGE;
}

Expand Down Expand Up @@ -910,7 +909,9 @@ static double dp3m_mc_time(char **log, int mesh, int cao, double r_cut_iL_min,

double const int_time = dp3m_mcr_time(mesh, cao, r_cut_iL, *_alpha_L);
if (int_time == -P3M_TUNE_FAIL) {
*log = strcat_alloc(*log, "tuning failed, test integration not possible\n");
if (verbose) {
std::printf("tuning failed, test integration not possible\n");
}
return int_time;
}

Expand All @@ -921,9 +922,10 @@ static double dp3m_mc_time(char **log, int mesh, int cao, double r_cut_iL_min,
}

/* print result */
sprintf(b, "%-4d %-3d %.5e %.5e %.5e %.3e %.3e %-8.0f\n", mesh, cao, r_cut_iL,
*_alpha_L, *_accuracy, rs_err, ks_err, int_time);
*log = strcat_alloc(*log, b);
if (verbose) {
std::printf("%-4d %-3d %.5e %.5e %.5e %.3e %.3e %-8.0f\n", mesh, cao,
r_cut_iL, *_alpha_L, *_accuracy, rs_err, ks_err, int_time);
}
return int_time;
}

Expand All @@ -932,7 +934,6 @@ static double dp3m_mc_time(char **log, int mesh, int cao, double r_cut_iL_min,
* @p _cao should contain an initial guess, which is then adapted by stepping
* up and down.
*
* @param[out] log log output
* @param[in] mesh @copybrief P3MParameters::mesh
* @param[in] cao_min lower bound for @p _cao
* @param[in] cao_max upper bound for @p _cao
Expand All @@ -943,13 +944,14 @@ static double dp3m_mc_time(char **log, int mesh, int cao, double r_cut_iL_min,
* @param[out] _r_cut_iL @copybrief P3MParameters::r_cut_iL
* @param[out] _alpha_L @copybrief P3MParameters::alpha_L
* @param[out] _accuracy @copybrief P3MParameters::accuracy
* @param[in] verbose printf output
*
* @returns The integration time in case of success, otherwise
* -@ref P3M_TUNE_FAIL or -@ref P3M_TUNE_CAO_TOO_LARGE */
static double dp3m_m_time(char **log, int mesh, int cao_min, int cao_max,
int *_cao, double r_cut_iL_min, double r_cut_iL_max,
static double dp3m_m_time(int mesh, int cao_min, int cao_max, int *_cao,
double r_cut_iL_min, double r_cut_iL_max,
double *_r_cut_iL, double *_alpha_L,
double *_accuracy) {
double *_accuracy, bool verbose) {
double best_time = -1, tmp_r_cut_iL = -1., tmp_alpha_L = 0.0,
tmp_accuracy = 0.0;
/* in which direction improvement is possible. Initially, we don't know it
Expand All @@ -963,8 +965,9 @@ static double dp3m_m_time(char **log, int mesh, int cao_min, int cao_max,
to increase cao to increase the obtainable precision of the far formula. */
double tmp_time;
do {
tmp_time = dp3m_mc_time(log, mesh, cao, r_cut_iL_min, r_cut_iL_max,
&tmp_r_cut_iL, &tmp_alpha_L, &tmp_accuracy);
tmp_time =
dp3m_mc_time(mesh, cao, r_cut_iL_min, r_cut_iL_max, &tmp_r_cut_iL,
&tmp_alpha_L, &tmp_accuracy, verbose);
/* bail out if the force evaluation is not working */
if (tmp_time == -P3M_TUNE_FAIL || tmp_time == -DP3M_RTBISECTION_ERROR)
return tmp_time;
Expand Down Expand Up @@ -1002,8 +1005,8 @@ static double dp3m_m_time(char **log, int mesh, int cao_min, int cao_max,
double dir_times[3];
for (final_dir = -1; final_dir <= 1; final_dir += 2) {
dir_times[final_dir + 1] = tmp_time =
dp3m_mc_time(log, mesh, cao + final_dir, r_cut_iL_min, r_cut_iL_max,
&tmp_r_cut_iL, &tmp_alpha_L, &tmp_accuracy);
dp3m_mc_time(mesh, cao + final_dir, r_cut_iL_min, r_cut_iL_max,
&tmp_r_cut_iL, &tmp_alpha_L, &tmp_accuracy, verbose);
/* bail out on errors, as usual */
if (tmp_time == -P3M_TUNE_FAIL || tmp_time == -DP3M_RTBISECTION_ERROR)
return tmp_time;
Expand Down Expand Up @@ -1053,8 +1056,9 @@ static double dp3m_m_time(char **log, int mesh, int cao_min, int cao_max,

/* move cao into the optimisation direction until we do not gain anymore. */
for (; cao >= cao_min && cao <= cao_max; cao += final_dir) {
tmp_time = dp3m_mc_time(log, mesh, cao, r_cut_iL_min, r_cut_iL_max,
&tmp_r_cut_iL, &tmp_alpha_L, &tmp_accuracy);
tmp_time =
dp3m_mc_time(mesh, cao, r_cut_iL_min, r_cut_iL_max, &tmp_r_cut_iL,
&tmp_alpha_L, &tmp_accuracy, verbose);
/* bail out on errors, as usual */
if (tmp_time == -P3M_TUNE_FAIL || tmp_time == -DP3M_RTBISECTION_ERROR)
return tmp_time;
Expand All @@ -1076,7 +1080,7 @@ static double dp3m_m_time(char **log, int mesh, int cao_min, int cao_max,
return best_time;
}

int dp3m_adaptive_tune(char **logger) {
int dp3m_adaptive_tune(bool verbose) {
/** Tuning of dipolar P3M. The algorithm basically determines the mesh, cao
* and then the real space cutoff, in this nested order.
*
Expand All @@ -1101,26 +1105,24 @@ int dp3m_adaptive_tune(char **logger) {
double alpha_L = -1, tmp_alpha_L = 0.0;
double accuracy = -1, tmp_accuracy = 0.0;
double time_best = 1e20, tmp_time;
char b[3 * ES_INTEGER_SPACE + 3 * ES_DOUBLE_SPACE + 128];

/* preparation */
mpi_call_all(dp3m_count_magnetic_particles);

/* Print Status */
sprintf(b,
"Dipolar P3M tune parameters: Accuracy goal = %.5e prefactor "
"= %.5e\n",
dp3m.params.accuracy, dipole.prefactor);
*logger = strcat_alloc(*logger, b);
sprintf(b, "System: box_l = %.5e # charged part = %d Sum[q_i^2] = %.5e\n",
box_geo.length()[0], dp3m.sum_dip_part, dp3m.sum_mu2);
*logger = strcat_alloc(*logger, b);

if (dp3m.sum_dip_part == 0) {
runtimeErrorMsg() << "no dipolar particles in the system";
return ES_ERROR;
}

/* Print Status */
if (verbose) {
std::printf("Dipolar P3M tune parameters: Accuracy goal = %.5e prefactor "
"= %.5e\n"
"System: box_l = %.5e # charged part = %d Sum[q_i^2] = %.5e\n",
dp3m.params.accuracy, dipole.prefactor, box_geo.length()[0],
dp3m.sum_dip_part, dp3m.sum_mu2);
}

/* parameter ranges */
if (dp3m.params.mesh[0] == 0) {
double expo;
Expand All @@ -1137,8 +1139,9 @@ int dp3m_adaptive_tune(char **logger) {
} else {
tmp_mesh = mesh_max = dp3m.params.mesh[0];

sprintf(b, "fixed mesh %d\n", dp3m.params.mesh[0]);
*logger = strcat_alloc(*logger, b);
if (verbose) {
std::printf("fixed mesh %d\n", dp3m.params.mesh[0]);
}
}

if (dp3m.params.r_cut_iL == 0.0) {
Expand All @@ -1152,8 +1155,9 @@ int dp3m_adaptive_tune(char **logger) {
} else {
r_cut_iL_min = r_cut_iL_max = dp3m.params.r_cut_iL;

sprintf(b, "fixed r_cut_iL %f\n", dp3m.params.r_cut_iL);
*logger = strcat_alloc(*logger, b);
if (verbose) {
std::printf("fixed r_cut_iL %f\n", dp3m.params.r_cut_iL);
}
}

if (dp3m.params.cao == 0) {
Expand All @@ -1163,18 +1167,22 @@ int dp3m_adaptive_tune(char **logger) {
} else {
cao_min = cao_max = cao = dp3m.params.cao;

sprintf(b, "fixed cao %d\n", dp3m.params.cao);
*logger = strcat_alloc(*logger, b);
if (verbose) {
std::printf("fixed cao %d\n", dp3m.params.cao);
}
}

if (verbose) {
std::printf("Dmesh cao Dr_cut_iL Dalpha_L Derr "
" Drs_err Dks_err time [ms]\n");
}
*logger = strcat_alloc(*logger, "Dmesh cao Dr_cut_iL Dalpha_L Derr "
" Drs_err Dks_err time [ms]\n");

/* mesh loop */
for (; tmp_mesh <= mesh_max; tmp_mesh += 2) {
tmp_cao = cao;
tmp_time =
dp3m_m_time(logger, tmp_mesh, cao_min, cao_max, &tmp_cao, r_cut_iL_min,
r_cut_iL_max, &tmp_r_cut_iL, &tmp_alpha_L, &tmp_accuracy);
tmp_time = dp3m_m_time(tmp_mesh, cao_min, cao_max, &tmp_cao, r_cut_iL_min,
r_cut_iL_max, &tmp_r_cut_iL, &tmp_alpha_L,
&tmp_accuracy, verbose);
/* some error occurred during the tuning force evaluation */
if (tmp_time == -P3M_TUNE_FAIL || tmp_time == -DP3M_RTBISECTION_ERROR)
return ES_ERROR;
Expand Down Expand Up @@ -1214,11 +1222,12 @@ int dp3m_adaptive_tune(char **logger) {
/* broadcast tuned p3m parameters */
mpi_bcast_coulomb_params();
/* Tell the user about the outcome */
sprintf(b,
"\nresulting parameters: mesh: %d, cao: %d, r_cut_iL: %.4e,"
"\n alpha_L: %.4e, accuracy: %.4e, time: %.0f\n",
mesh, cao, r_cut_iL, alpha_L, accuracy, time_best);
*logger = strcat_alloc(*logger, b);
if (verbose) {
std::printf(
"\nresulting parameters: mesh: %d, cao: %d, r_cut_iL: %.4e,"
"\n alpha_L: %.4e, accuracy: %.4e, time: %.0f\n",
mesh, cao, r_cut_iL, alpha_L, accuracy, time_best);
}
return ES_OK;
}

Expand Down
4 changes: 2 additions & 2 deletions src/core/electrostatics_magnetostatics/p3m-dipolar.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -157,11 +157,11 @@ void dp3m_deactivate();
* The function is based on routines of the program HE_Q.cpp for charges
* written by M. Deserno.
*
* @param[out] log log output
* @param verbose printf output
* @retval ES_OK
* @retval ES_ERROR
*/
int dp3m_adaptive_tune(char **log);
int dp3m_adaptive_tune(bool verbose);

/** Compute the k-space part of forces and energies for the magnetic
* dipole-dipole interaction
Expand Down
Loading