Hi,
To try and check that the units are fine you can print out the value of
GravAccel before and after the
addition to see if there is some kind of order of magnitude difference.
But my bet is that there is some particle very close to the origin that
results in a very small R_1,
maybe you can add a check to see if R_1 is smaller than the softening
length, similar to what is
done between lines 1356 ad 1371 of forcetree.c .
Antonio
2015-09-21 4:46 GMT+10:00 Jesus <jesusms_at_astro.ucla.edu>:
> I have written my own potential at the end of gravity_tree(), just as
> instructed in
> http://wwwmpa.mpa-garching.mpg.de/gadget/gadget-list/0232.html
> However, Gadget is giving me the error:
>
> Error: A timestep of size zero was assigned on the integer timeline!
> We better stop.
> Task=5 Part-ID=285634 dt=7.4194e-09 tibase=7.45058e-09 ti_step=0
> ac=1.12005e+12 xyz=(3.24349|-7.37641|0.00861514)
> tree=(-4.50823e+11|1.02527e+12-1.54804e+09)
>
> hydro-frc=(4.55549e+08|2.50445e+08|-2.86563e+08)
> task 5: endrun called with an error level of 818
>
>
> --------------------------------------------------------------------------
> MPI_ABORT was invoked on rank 5 in communicator MPI_COMM_WORLD
> with errorcode 818.
>
> NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
> You may or may not see output from other processes, depending on
> exactly when Open MPI kills them.
> --------------------------------------------------------------------------
> --------------------------------------------------------------------------
> mpirun has exited due to process rank 5 with PID 38354 on….
>
>
> The only problem I can think of is that somehow I did not use the correct
> units that the computes with.
> This is what I have
>
> for(i = 0; i < NumPart; i++)
> if(P[i].Ti_endstep == All.Ti_Current)
> { double R_2 = ( (P[i].Pos[0]*P[i].Pos[0]) + (P[i].Pos[1]*P[i].Pos[1])
> + (P[i].Pos[2]*P[i].Pos[2]) ) ;
> double R_1 = sqrt(R_2);
> double R_3 = R_1 * R_2;
>
>
> P[i].GravAccel[0] += ( (my_potential*conversion) - (All.G*40.0
> *P[i].Pos[0]/R_3) ) ;
> etc
>
> The second term is just the potential of a point mass at the origin with
> mass = 40 mass units (code units).
> my_potential() function is a custom function I added to Gadget to
> calculate another potential term, which result is in km/s^2. The
> "conversion" converts the km/s^2 into code units, (first converting to
> cm/s^2, and then to code units by dividing by 1 acceleration unit, which I
> computed by doing unit_velocity_in_cm_per_s/unit_time_in_seconds).
>
> Am I doing the right unit conversion?
>
> Thank you for your time.
>
> Jesus Salas
>
>
> -----------------------------------------------------------
>
> If you wish to unsubscribe from this mailing, send mail to
> minimalist_at_MPA-Garching.MPG.de with a subject of: unsubscribe gadget-list
> A web-archive of this mailing list is available here:
> http://www.mpa-garching.mpg.de/gadget/gadget-list
>
>
Received on 2015-09-21 01:37:38