diff --git a/examples/TestingGravity/makeIC.py b/examples/TestingGravity/makeIC.py index 58fc51cc96025d13dae14b8c1a5afd811715ad0c..e344842bb268be148111ef126b3fb9a91225d946 100644 --- a/examples/TestingGravity/makeIC.py +++ b/examples/TestingGravity/makeIC.py @@ -34,14 +34,14 @@ PROTON_MASS_IN_CGS = 1.6726231e24 YEAR_IN_CGS = 3.154e+7 # choice of units -const_unit_length_in_cgs = (1000 * PARSEC_IN_CGS ) -const_unit_mass_in_cgs = (SOLAR_MASS_IN_CGS) -const_unit_velocity_in_cgs = (1e5) +const_unit_length_in_cgs = (1000*PARSEC_IN_CGS) +const_unit_mass_in_cgs = (SOLAR_MASS_IN_CGS) +const_unit_velocity_in_cgs = (1e5) # derived units const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs) -const_G = (NEWTON_GRAVITY_CGS * const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/ (const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs)) - +const_G = ((NEWTON_GRAVITY_CGS*const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/(const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs))) +print 'G=', const_G # Parameters diff --git a/src/const.h b/src/const.h index 7e08e6b266c0d3ee0c51052daf1d69298abb4518..7ccd420bdf996500aea9f987d1798569b91fa174 100644 --- a/src/const.h +++ b/src/const.h @@ -80,13 +80,13 @@ /* System of units */ -#define const_unit_length_in_cgs (1000 * PARSEC_IN_CGS ) /* kpc */ -#define const_unit_mass_in_cgs (SOLAR_MASS_IN_CGS) /* solar mass */ -#define const_unit_velocity_in_cgs (1e5) /* km/s */ +#define const_unit_length_in_cgs (1000*PARSEC_IN_CGS) /* kpc */ +#define const_unit_mass_in_cgs (SOLAR_MASS_IN_CGS) /* solar mass */ +#define const_unit_velocity_in_cgs (1e5) /* km/s */ /* Derived constants */ #define const_unit_time_in_cgs (const_unit_length_in_cgs / const_unit_velocity_in_cgs) -#define const_G (NEWTON_GRAVITY_CGS * const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/ (const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs)) +#define const_G ((NEWTON_GRAVITY_CGS*const_unit_mass_in_cgs/(const_unit_velocity_in_cgs*const_unit_velocity_in_cgs*const_unit_length_in_cgs))) /* External Potential Constants */ #define External_Potential_X (50000 * PARSEC_IN_CGS / const_unit_length_in_cgs) diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h index 216ea66ef5baf3ad5279513a413aa1061c170452..89dea853d27c59860710967b1cb783d0ffb70bae 100644 --- a/src/gravity/Default/gravity_part.h +++ b/src/gravity/Default/gravity_part.h @@ -43,6 +43,11 @@ struct gpart { /* Particle time of end of time-step. */ int ti_end; + + /* current time of x, and of v_full */ + float tx=0; + float tv=0; + /* Anonymous union for id/part. */ union { diff --git a/src/gravity/ExternalPotential/gravity_part.h b/src/gravity/ExternalPotential/gravity_part.h index 2e8e5878b1f07d426ef8ce38935d15abcc92d6f9..0505a6a9aaf4540b4e67f76abfb1bc9f910bec8e 100644 --- a/src/gravity/ExternalPotential/gravity_part.h +++ b/src/gravity/ExternalPotential/gravity_part.h @@ -55,6 +55,10 @@ struct gpart { /* Particle time of end of time-step. */ int ti_end; + /* current time of x, and of v_full */ + float tx; + float tv; + /* Anonymous union for id/part. */ union { diff --git a/src/runner.c b/src/runner.c index f99998a90603f703fc47daabbd9fa8a81c04cd45..3d760ea09070f1bb3c62beb6be481dd25ca9839b 100644 --- a/src/runner.c +++ b/src/runner.c @@ -100,6 +100,7 @@ void runner_dograv_external(struct runner *r, struct cell *c) { struct gpart *g, *gparts = c->gparts; float rinv; float L[3], E; + float dx, dy, dz; int i, k, gcount = c->gcount; const int ti_current = r->e->ti_current; @@ -130,21 +131,27 @@ void runner_dograv_external(struct runner *r, struct cell *c) { /* Is this part within the time step? */ if (g->ti_end <= ti_current) { - rinv = 1 / sqrtf((g->x[0]-External_Potential_X)*(g->x[0]-External_Potential_X) + (g->x[1]-External_Potential_Y)*(g->x[1]-External_Potential_Y) + (g->x[2]-External_Potential_Z)*(g->x[2]-External_Potential_Z)); + dx = g->x[0]-External_Potential_X; + dy = g->x[1]-External_Potential_Y; + dz = g->x[2]-External_Potential_Z; + + rinv = 1 / sqrtf(dx*dx + dy*dy + dz*dz); /* check for energy and angular momentum conservation */ E = 0.5 * ((g->v_full[0]*g->v_full[0]) + (g->v_full[1]*g->v_full[1]) + (g->v_full[2]*g->v_full[2])) - const_G * External_Potential_Mass * rinv; - L[0] = (g->x[1] - External_Potential_X)*g->v_full[2] - (g->x[2] - External_Potential_X)*g->v_full[1]; - L[1] = (g->x[2] - External_Potential_Y)*g->v_full[0] - (g->x[0] - External_Potential_Y)*g->v_full[2]; - L[2] = (g->x[0] - External_Potential_Z)*g->v_full[1] - (g->x[1] - External_Potential_Z)*g->v_full[0]; + L[0] = dy * g->v_full[2] - dz * g->v_full[1]; + L[1] = dz * g->v_full[0] - dx * g->v_full[2]; + L[2] = dx * g->v_full[1] - dy * g->v_full[0]; if(g->id == 0) { - message("update %f\t %f\t %f\t %f\t %f\t %f\t %f\t %f\t %f\n", r->e->time, 1./rinv, g->x[0], g->x[1], g->x[2], E, L[0], L[1], L[2]); - message(" ... %f %f %f\n", g->v_full[0], g->v_full[1], g->v_full[2]); - message(" ... %e %e\n", const_G, External_Potential_Mass); + float v2 = g->v_full[0]*g->v_full[0] + g->v_full[1]*g->v_full[1] + g->v_full[2]*g->v_full[2]; + float fg = const_G * External_Potential_Mass * rinv; + //message("grav_external time= %f\t V_c^2= %f GM/r= %f E= %f L[2]= %f x= %f y= %f vx= %f vy= %f\n", r->e->time, v2, fg, E, L[2], g->x[0], g->x[1], g->v_full[0], g->v_full[1]); + message("%f\t %f %f %f %f %f %f %f %f %f %f\n", r->e->time, g->tx, g->tv, v2, fg, E, L[2], g->x[0], g->x[1], g->v_full[0], g->v_full[1]); + // message(" G=%e M=%e\n", const_G, External_Potential_Mass); } - g->a_grav_external[0] = - const_G * External_Potential_Mass * (g->x[0] - External_Potential_X) * rinv * rinv * rinv; - g->a_grav_external[1] = - const_G * External_Potential_Mass * (g->x[1] - External_Potential_Y) * rinv * rinv * rinv; - g->a_grav_external[2] = - const_G * External_Potential_Mass * (g->x[2] - External_Potential_Z) * rinv * rinv * rinv; + g->a_grav_external[0] = - const_G * External_Potential_Mass * dx * rinv * rinv * rinv; + g->a_grav_external[1] = - const_G * External_Potential_Mass * dy * rinv * rinv * rinv; + g->a_grav_external[2] = - const_G * External_Potential_Mass * dz * rinv * rinv * rinv; } } @@ -841,8 +848,13 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) { /* Loop over all gparts in the cell */ for (int k = 0; k < nr_gparts; k++) { - g = &gparts[k]; /* nothing to do */ - g->x[0] *=1; + g = &gparts[k]; + if(g->id == 0) + message(" dt= %f tx= %f tv=%f \n",dt, g->tx, g->tv); + g->x[0] += g->v_full[0] * dt; + g->x[1] += g->v_full[1] * dt; + g->x[2] += g->v_full[2] * dt; + g->tx += dt; } } @@ -1070,11 +1082,15 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) { /* Move particle forward in time */ g->ti_begin = g->ti_end; g->ti_end = g->ti_begin + new_dti; - + + if(g->id == 0) + message(" dt= %f tx= %f tv=%f \n",dt, g->tx, g->tv); + /* Kick particles in momentum space */ g->v_full[0] += g->a_grav_external[0] * dt; g->v_full[1] += g->a_grav_external[1] * dt; g->v_full[2] += g->a_grav_external[2] * dt; + g->tv += dt; /* Minimal time for next end of time-step */ ti_end_min = min(g->ti_end, ti_end_min); diff --git a/src/single_io.c b/src/single_io.c index 52fbb4df0d5f1cfd50dedd5abeb35eeeb6c5ba4c..1d003c850c45fbb5314d29877f7221c912433cc1 100644 --- a/src/single_io.c +++ b/src/single_io.c @@ -411,6 +411,8 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, struct g darkmatter_read_particles(h_grp, *Ndm, *Ndm, 0, *gparts); for(int i=0; i<*Ndm; i++) { (*gparts)[i].id = -abs( (*gparts)[i].id ); + (*gparts)[i].tx = 0; + (*gparts)[i].tv = 0; } }