Removed redundant copying from calc_self_grav and got rid of calc_self_grav_zero
This commit is contained in:
Yohai Meiron
parent
2c05355eaa
commit
4b8ef0ca61
1 changed files with 36 additions and 147 deletions
183
phigrape.cpp
183
phigrape.cpp
|
|
@ -145,6 +145,8 @@ Last redaction : 2019.04.16 12:55
|
|||
#include <sys/time.h>
|
||||
#include <sys/resource.h>
|
||||
|
||||
#include <algorithm>
|
||||
|
||||
/*
|
||||
double aaa;
|
||||
double aaapars[5];
|
||||
|
|
@ -158,12 +160,12 @@ extern void qwerty_(double *aaa, double *aaapars);
|
|||
#include <mpi.h>
|
||||
|
||||
/* Some "good" functions and constants... */
|
||||
#define SIG(x) ( ((x)<0) ? (-1):(1) )
|
||||
#define ABS(x) ( ((x)<0) ? (-x):(x) )
|
||||
#define MAX(a,b) ( ((a)>(b)) ? (a):(b) )
|
||||
#define MIN(a,b) ( ((a)<(b)) ? (a):(b) )
|
||||
#define SQR(x) ( (x)*(x) )
|
||||
#define POW3(x) ( (x)*SQR(x) )
|
||||
#define SIG(x) (((x)<0) ? (-1):(1) )
|
||||
#define ABS(x) (((x)<0) ? (-x):(x) )
|
||||
#define MAX(a,b) (((a)>(b)) ? (a):(b) )
|
||||
#define MIN(a,b) (((a)<(b)) ? (a):(b) )
|
||||
#define SQR(x) ((x)*(x) )
|
||||
#define POW3(x) ((x)*SQR(x) )
|
||||
|
||||
#define Pi 3.14159265358979323846
|
||||
#define TWOPi 6.283185307179
|
||||
|
|
@ -477,7 +479,7 @@ return - 0 if everything OK
|
|||
|
||||
double my_rand(void)
|
||||
{
|
||||
return( (double)(rand()/(double)RAND_MAX) );
|
||||
return((double)(rand()/(double)RAND_MAX) );
|
||||
}
|
||||
|
||||
double my_rand2(void)
|
||||
|
|
@ -952,101 +954,6 @@ fclose(out);
|
|||
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
|
||||
void calc_self_grav_zero()
|
||||
{
|
||||
|
||||
/* calc the grav for the active particles */
|
||||
|
||||
g6_set_ti(clusterid, time_cur);
|
||||
|
||||
ni = n_act;
|
||||
|
||||
/* define the local phi, a, adot for these active particles */
|
||||
|
||||
for (i=0; i<ni; i+=npipe)
|
||||
{
|
||||
|
||||
nn = npipe;
|
||||
if (ni-i < npipe) nn = ni - i;
|
||||
|
||||
for (ii=0; ii<nn; ii++)
|
||||
{
|
||||
|
||||
index_i[ii] = ind_act[i+ii];
|
||||
h2_i[ii] = eps2;
|
||||
|
||||
for (k=0;k<3;k++)
|
||||
{
|
||||
x_i[ii][k] = x_act[i+ii][k]; v_i[ii][k] = v_act[i+ii][k];
|
||||
} /* k */
|
||||
|
||||
p_i[ii] = -1.0;
|
||||
|
||||
for (k=0;k<3;k++)
|
||||
{
|
||||
a_i[ii][k] = 10.0;
|
||||
jerk_i[ii][k] = 100.0;
|
||||
} /* k */
|
||||
|
||||
} /* ii */
|
||||
|
||||
g6calc_firsthalf(clusterid, n_loc, nn, index_i, x_i, v_i , a_i, jerk_i, p_i, eps2, h2_i);
|
||||
g6calc_lasthalf(clusterid, n_loc, nn, index_i, x_i, v_i, eps2, h2_i, a_i, jerk_i, p_i);
|
||||
|
||||
for (ii=0; ii<nn; ii++)
|
||||
{
|
||||
if (ABS(jerk_i[ii][0]) < 1.0E-05) jerk_i[ii][0] = 1.0E-05;
|
||||
if (ABS(jerk_i[ii][1]) < 1.0E-05) jerk_i[ii][1] = 1.0E-05;
|
||||
if (ABS(jerk_i[ii][2]) < 1.0E-05) jerk_i[ii][2] = 1.0E-05;
|
||||
|
||||
if (ABS(jerk_i[ii][0]) > 100.0) jerk_i[ii][0] = 100.0;
|
||||
if (ABS(jerk_i[ii][1]) > 100.0) jerk_i[ii][1] = 100.0;
|
||||
if (ABS(jerk_i[ii][2]) > 100.0) jerk_i[ii][2] = 100.0;
|
||||
}
|
||||
|
||||
g6calc_firsthalf(clusterid, n_loc, nn, index_i, x_i, v_i , a_i, jerk_i, p_i, eps2, h2_i);
|
||||
g6calc_lasthalf(clusterid, n_loc, nn, index_i, x_i, v_i, eps2, h2_i, a_i, jerk_i, p_i);
|
||||
|
||||
for (ii=0; ii<nn; ii++)
|
||||
{
|
||||
|
||||
pot_act_tmp[i+ii] = p_i[ii];
|
||||
|
||||
for (k=0;k<3;k++)
|
||||
{
|
||||
a_act_tmp[i+ii][k] = a_i[ii][k];
|
||||
adot_act_tmp[i+ii][k] = jerk_i[ii][k];
|
||||
} /* k */
|
||||
|
||||
} /* ii */
|
||||
|
||||
} /* i */
|
||||
|
||||
/* Store the value of the local partial force etc... */
|
||||
|
||||
for (i=0; i<n_act; i++)
|
||||
{
|
||||
|
||||
pot_act_tmp_loc[i] = pot_act_tmp[i];
|
||||
|
||||
for (k=0;k<3;k++)
|
||||
{
|
||||
a_act_tmp_loc[i][k] = a_act_tmp[i][k];
|
||||
adot_act_tmp_loc[i][k] = adot_act_tmp[i][k];
|
||||
} /* k */
|
||||
|
||||
} /* i */
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
void calc_ext_grav_zero()
|
||||
{
|
||||
|
||||
|
|
@ -1194,7 +1101,7 @@ for (i=0; i<ni; i++)
|
|||
dum3 = POW3(dum);
|
||||
dum_g = pow( dum , -g_ext );
|
||||
|
||||
pot_ext[i] -= ( (m_ext/r_ext) / (2.0-g_ext) ) * ( 1.0 - dum2 * dum_g );
|
||||
pot_ext[i] -= ((m_ext/r_ext) / (2.0-g_ext) ) * ( 1.0 - dum2 * dum_g );
|
||||
|
||||
tmp = (m_ext/tmp_r3) * dum3 * dum_g;
|
||||
|
||||
|
|
@ -1205,7 +1112,7 @@ for (i=0; i<ni; i++)
|
|||
rv_ij = ( vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij );
|
||||
|
||||
tmp = ( m_ext/POW3(tmp_r+r_ext) ) * dum_g;
|
||||
tmp_g = ( (r_ext*g_ext + 3.0*tmp_r)/(tmp_r2*(tmp_r+r_ext)) ) * rv_ij;
|
||||
tmp_g = ((r_ext*g_ext + 3.0*tmp_r)/(tmp_r2*(tmp_r+r_ext)) ) * rv_ij;
|
||||
|
||||
adot[i][0] += tmp * ( tmp_g * x_ij - vx_ij );
|
||||
adot[i][1] += tmp * ( tmp_g * y_ij - vy_ij );
|
||||
|
|
@ -1296,17 +1203,15 @@ for (i=0; i<ni; i++)
|
|||
|
||||
}
|
||||
|
||||
void calc_self_grav()
|
||||
void calc_self_grav(double t)
|
||||
{
|
||||
/* calc the new grav for the active particles */
|
||||
|
||||
// TODO possibly redundant copying in this function
|
||||
|
||||
#ifdef TIMING
|
||||
get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
|
||||
#endif
|
||||
|
||||
g6_set_ti(clusterid, min_t);
|
||||
g6_set_ti(clusterid, t);
|
||||
|
||||
ni = n_act;
|
||||
|
||||
|
|
@ -1316,30 +1221,11 @@ void calc_self_grav()
|
|||
nn = npipe;
|
||||
if (ni-i < npipe) nn = ni - i;
|
||||
for (ii=0; ii<nn; ii++) {
|
||||
iii = ind_act[i+ii];
|
||||
index_i[ii] = iii;
|
||||
h2_i[ii] = eps2;
|
||||
for (k=0;k<3;k++) {
|
||||
x_i[ii][k] = x_act_new[i+ii][k];
|
||||
v_i[ii][k] = v_act_new[i+ii][k];
|
||||
} /* k */
|
||||
p_i[ii] = pot_act_tmp_loc[iii];
|
||||
for (k=0;k<3;k++) {
|
||||
a_i[ii][k] = a_act_tmp_loc[iii][k];
|
||||
jerk_i[ii][k] = adot_act_tmp_loc[iii][k];
|
||||
} /* k */
|
||||
h2_i[ii] = eps2; // TODO This should be a global or something
|
||||
} /* ii */
|
||||
g6calc_firsthalf(clusterid, n_loc, nn, index_i, x_i, v_i , a_i, jerk_i, p_i, eps2, h2_i);
|
||||
g6calc_lasthalf(clusterid, n_loc, nn, index_i, x_i, v_i, eps2, h2_i, a_i, jerk_i, p_i);
|
||||
g6calc_firsthalf(clusterid, n_loc, nn, ind_act+i, x_act_new+i, v_act_new+i, a_act_tmp+i, adot_act_tmp+i, pot_act_tmp+i, eps2, h2_i);
|
||||
g6calc_lasthalf( clusterid, n_loc, nn, ind_act+i, x_act_new+i, v_act_new+i, eps2, h2_i, a_act_tmp+i, adot_act_tmp+i, pot_act_tmp+i);
|
||||
g6_calls++;
|
||||
|
||||
for (ii=0; ii<nn; ii++) {
|
||||
pot_act_tmp[i+ii] = p_i[ii];
|
||||
for (k=0;k<3;k++) {
|
||||
a_act_tmp[i+ii][k] = a_i[ii][k];
|
||||
adot_act_tmp[i+ii][k] = jerk_i[ii][k];
|
||||
} /* k */
|
||||
} /* ii */
|
||||
} /* i */
|
||||
|
||||
/* Store the new value of the local partial force etc... */
|
||||
|
|
@ -1348,8 +1234,8 @@ void calc_self_grav()
|
|||
iii = ind_act[i];
|
||||
pot_act_tmp_loc[iii] = pot_act_tmp[i];
|
||||
for (k=0;k<3;k++) {
|
||||
a_act_tmp_loc[iii][k] = a_act_tmp[i][k];
|
||||
adot_act_tmp_loc[iii][k] = adot_act_tmp[i][k];
|
||||
a_act_tmp_loc[iii][k] = a_act_tmp[i][k];
|
||||
adot_act_tmp_loc[iii][k] = adot_act_tmp[i][k];
|
||||
} /* k */
|
||||
} /* i */
|
||||
|
||||
|
|
@ -1510,7 +1396,7 @@ for (i=0; i<ni; i++)
|
|||
dum3 = POW3(dum);
|
||||
dum_g = pow( dum , -g_ext );
|
||||
|
||||
pot_act_ext[i] -= ( (m_ext/r_ext) / (2.0-g_ext) ) * ( 1.0 - dum2 * dum_g );
|
||||
pot_act_ext[i] -= ((m_ext/r_ext) / (2.0-g_ext) ) * ( 1.0 - dum2 * dum_g );
|
||||
|
||||
tmp = (m_ext/tmp_r3) * dum3 * dum_g;
|
||||
|
||||
|
|
@ -1521,7 +1407,7 @@ for (i=0; i<ni; i++)
|
|||
rv_ij = ( vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij );
|
||||
|
||||
tmp = ( m_ext/POW3(tmp_r+r_ext) ) * dum_g;
|
||||
tmp_g = ( (r_ext*g_ext + 3.0*tmp_r)/(tmp_r2*(tmp_r+r_ext)) ) * rv_ij;
|
||||
tmp_g = ((r_ext*g_ext + 3.0*tmp_r)/(tmp_r2*(tmp_r+r_ext)) ) * rv_ij;
|
||||
|
||||
adot_act_new[i][0] += tmp * ( tmp_g * x_ij - vx_ij );
|
||||
adot_act_new[i][1] += tmp * ( tmp_g * y_ij - vy_ij );
|
||||
|
|
@ -1735,8 +1621,8 @@ E_tot = E_pot + E_kin + E_pot_ext;
|
|||
E_tot_corr = E_tot + E_corr;
|
||||
E_tot_corr_sd = E_tot + E_corr + E_sd;
|
||||
|
||||
// if ( (Timesteps == 0.0) && (time_cur == 0.0) )
|
||||
if ( (Timesteps == 0.0) )
|
||||
// if ((Timesteps == 0.0) && (time_cur == 0.0) )
|
||||
if (Timesteps == 0.0)
|
||||
{
|
||||
/* initialize the E_tot_0 + etc... */
|
||||
|
||||
|
|
@ -1997,7 +1883,7 @@ int main(int argc, char *argv[])
|
|||
|
||||
fflush(stdout);
|
||||
|
||||
if ( (diskstep == 0) && (time_cur == 0.0) ) {
|
||||
if ((diskstep == 0) && (time_cur == 0.0)) {
|
||||
out = fopen("contr.dat","w");
|
||||
fclose(out);
|
||||
|
||||
|
|
@ -2249,9 +2135,12 @@ int main(int argc, char *argv[])
|
|||
t_act[i] = t[iii];
|
||||
dt_act[i] = dt[iii];
|
||||
|
||||
} /* i */
|
||||
} /* i */
|
||||
|
||||
calc_self_grav_zero();
|
||||
// NOTE this is where calc_self_grav_zero() used to be. Hopefully the copying from *_act to *_act_new is just a temporary measure.
|
||||
std::copy(&x_act[0][0], &x_act[0][0]+3*N, &x_act_new[0][0]);
|
||||
std::copy(&v_act[0][0], &v_act[0][0]+3*N, &v_act_new[0][0]);
|
||||
calc_self_grav(time_cur);
|
||||
|
||||
/* Wait to all processors to finish his works... */
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
|
|
@ -2604,7 +2493,7 @@ int main(int argc, char *argv[])
|
|||
n_act = 0;
|
||||
|
||||
for (i=0; i<N; i++) {
|
||||
if ( (t[i]+dt[i]) == min_t ) {
|
||||
if (t[i]+dt[i] == min_t ) {
|
||||
ind_act[n_act] = i;
|
||||
n_act++;
|
||||
}
|
||||
|
|
@ -2702,7 +2591,7 @@ int main(int argc, char *argv[])
|
|||
DT_ACT_PRED += (CPU_tmp_user - CPU_tmp_user0);
|
||||
#endif
|
||||
|
||||
calc_self_grav();
|
||||
calc_self_grav(min_t);
|
||||
|
||||
/* Reduce the "global" vectors from "local" on all the nodes */
|
||||
|
||||
|
|
@ -2851,7 +2740,7 @@ int main(int argc, char *argv[])
|
|||
a3dot1abs = sqrt(a3[0]*a3[0]+a3[1]*a3[1]+a3[2]*a3[2]);
|
||||
|
||||
#ifdef ADD_BH2
|
||||
if ( (ind_act[i] == 0) || (ind_act[i] == 1) ) {
|
||||
if ((ind_act[i] == 0) || (ind_act[i] == 1)) {
|
||||
dt_new = sqrt(eta_bh*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
|
||||
} else {
|
||||
dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
|
||||
|
|
@ -2859,7 +2748,7 @@ int main(int argc, char *argv[])
|
|||
#else
|
||||
|
||||
#ifdef ADD_BH1
|
||||
if ( (ind_act[i] == 0) ) {
|
||||
if (ind_act[i] == 0) {
|
||||
dt_new = sqrt(eta_bh*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
|
||||
} else {
|
||||
dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs));
|
||||
|
|
@ -2872,7 +2761,7 @@ int main(int argc, char *argv[])
|
|||
|
||||
if (dt_new < dt_min) dt_tmp = dt_min;
|
||||
|
||||
if ( (dt_new < dt_tmp) && (dt_new > dt_min) ) {
|
||||
if ((dt_new < dt_tmp) && (dt_new > dt_min)) {
|
||||
power = log(dt_new)/log(2.0) - 1;
|
||||
|
||||
dt_tmp = pow(2.0, (double)power);
|
||||
|
|
@ -2965,9 +2854,9 @@ int main(int argc, char *argv[])
|
|||
|
||||
iii = ind_act[i];
|
||||
|
||||
if ( (tmp_r2 < SEMI_a2*R_INF2) ) {
|
||||
if (tmp_r2 < SEMI_a2*R_INF2) {
|
||||
|
||||
if ( (inf_event[iii] == 0) ) {
|
||||
if (inf_event[iii] == 0) {
|
||||
|
||||
fprintf(out,"INF1 %.6E %.16E %07d %07d %.6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E \n",
|
||||
Timesteps, time_cur, i, ind_act[i],
|
||||
|
|
@ -2982,7 +2871,7 @@ int main(int argc, char *argv[])
|
|||
}
|
||||
|
||||
} else {
|
||||
if ( (inf_event[iii] == 1) ) {
|
||||
if (inf_event[iii] == 1) {
|
||||
fprintf(out,"INF2 %.6E %.16E %07d %07d %.6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E %.6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E %.6E \n",
|
||||
Timesteps, time_cur, i, ind_act[i],
|
||||
sqrt(DR2), x_bbhc[0], x_bbhc[1], x_bbhc[2], v_bbhc[0], v_bbhc[1], v_bbhc[2],
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue