GADGET-4
gravity.cc
Go to the documentation of this file.
1/*******************************************************************************
2 * \copyright This file is part of the GADGET4 N-body/SPH code developed
3 * \copyright by Volker Springel. Copyright (C) 2014-2020 by Volker Springel
4 * \copyright (vspringel@mpa-garching.mpg.de) and all contributing authors.
5 *******************************************************************************/
6
12#include "gadgetconfig.h"
13
14#include <math.h>
15#include <mpi.h>
16#include <stdlib.h>
17#include <string.h>
18
19#include "../data/allvars.h"
20#include "../data/dtypes.h"
21#include "../data/intposconvert.h"
22#include "../data/mymalloc.h"
23#include "../domain/domain.h"
24#include "../fmm/fmm.h"
25#include "../gravtree/gravtree.h"
26#include "../logs/logs.h"
27#include "../logs/timer.h"
28#include "../main/simulation.h"
29#include "../mpi_utils/mpi_utils.h"
30#include "../pm/pm.h"
31#include "../system/system.h"
32#include "../time_integration/timestep.h"
33
42{
44
45 mpi_printf("ACCEL: Start tree gravity force computation... (%lld particles)\n", Sp.TimeBinsGravity.GlobalNActiveParticles);
46
48 {
49 GravTree.DoEwald = 0;
50
51#ifdef PMGRID
52 GravTree.DoPM = 0; /* default value */
53#endif
54
55#if defined(PMGRID) && defined(PERIODIC) && \
56 !defined(TREEPM_NOTIMESPLIT) /* classic TreePM in periodic box with time integration split */
57 GravTree.DoEwald = 0;
59
60#ifdef PLACEHIGHRESREGION
61 if(Sp.TimeBinsGravity.GlobalNActiveParticles > All.ActivePartFracForPMinsteadOfEwald * Sp.TotNumPart)
63#endif
64
65#else /* everything else */
66
67#if defined(PMGRID) /* with PM acceleration (we force here TREEPM_NOTIMESPLIT to be set) */
68
69 if(Sp.TimeBinsGravity.GlobalNActiveParticles > All.ActivePartFracForPMinsteadOfEwald * Sp.TotNumPart)
70 {
71 GravTree.DoEwald = 0;
72#ifdef PLACEHIGHRESREGION
74#else
76#endif
77 }
78 else
79 {
80 GravTree.DoPM = 0;
81#if defined(PERIODIC)
82 GravTree.DoEwald = 1;
83#endif
84 }
85
86#else /* here no PM acceleration is used */
87
88#if defined(PERIODIC)
89 GravTree.DoEwald = 1; /* periodic boundaries with Ewald summation */
90#endif
91
92#endif
93
94#endif
95
96#ifdef SECOND_ORDER_LPT_ICS
97 if(All.Ti_Current == 0)
98 second_order_ic_correction();
99#endif
100
102 {
103 /* For the first timestep, we do one gravity calculation up front
104 * with the Barnes & Hut Criterion to allow usage of relative opening
105 * criterion with consistent accuracy.
106 */
107#if defined(PMGRID) && defined(PERIODIC) && !defined(TREEPM_NOTIMESPLIT)
109#endif
110 gravity(timebin);
111
112 gravity_set_oldacc(timebin);
113
114 /* now may switch on relative opening criterion since we have an old acceleration */
117 }
118
119 gravity(timebin); /* computes gravity acceleration */
120
121#ifdef FORCETEST
122#if defined(PMGRID) && defined(PERIODIC) && !defined(TREEPM_NOTIMESPLIT)
123 if(timebin == All.HighestOccupiedGravTimeBin)
124#endif
125 {
127 GravTest.gravity_forcetest(timebin);
129
130 if(FORCETEST >= 2.0) /* this is for a special test where we repeat the calculation */
131 {
132 for(int i = 0; i < (int)(FORCETEST)-1; i++)
133 {
135 Domain.domain_free();
136 Domain.domain_decomposition(STANDARD);
138 NgbTree.treebuild(Sp.NumGas, NULL);
139
140 gravity(timebin);
141
143 GravTest.gravity_forcetest(timebin);
145 }
146 }
147
148 if(FORCETEST > 1.0)
149 endrun();
150 }
151#endif
152 }
153
154 mpi_printf("ACCEL: tree force computation done.\n");
155
156 if(All.TimeLimitCPU == 0)
157 endrun();
158}
159
171void sim::gravity(int timebin)
172{
173 if(timebin == All.HighestOccupiedGravTimeBin)
175 else
177
178 /* let's first initialize the results */
179 for(int i = 0; i < Sp.TimeBinsGravity.NActiveParticles; i++)
180 {
181 int target = Sp.TimeBinsGravity.ActiveParticleList[i];
182#ifdef EVALPOTENTIAL
183 Sp.P[target].Potential = 0;
184#endif
185 for(int j = 0; j < 3; j++)
186 Sp.P[target].GravAccel[j] = 0;
187
189 Sp.P[target].GravCost = 0;
190 }
191
192#ifdef SELFGRAVITY
193 /* set new softening lengths on global steps to take into account possible cosmological time variation */
194 if(timebin == All.HighestOccupiedGravTimeBin)
196
197#if defined(PMGRID) && (defined(TREEPM_NOTIMESPLIT) || defined(PLACEHIGHRESREGION))
199 gravity_pm(timebin);
200
201#if defined(FORCETEST) && defined(PLACEHIGHRESREGION)
202 for(int i = 0; i < Sp.TimeBinsGravity.NActiveParticles; i++)
203 {
204 int target = Sp.TimeBinsGravity.ActiveParticleList[i];
205 Sp.P[target].PotentialHPM = All.G * Sp.P[target].Potential;
206 for(int j = 0; j < 3; j++)
207 Sp.P[target].GravAccelHPM[j] = All.G * Sp.P[target].GravAccel[j];
208 }
209#endif
210
211#endif
212
213#ifdef ALLOW_DIRECT_SUMMATION
215 {
216 GravTree.gravity_direct(&Sp, &Domain, timebin);
217 }
218 else
219#endif
220 {
222
223#ifdef HIERARCHICAL_GRAVITY
225#else
227#endif
228
229#ifdef FMM
230 GravTree.gravity_fmm(timebin);
231#else
232 GravTree.gravity_tree(timebin);
233#endif
234
236 }
237
238 /* now multplify with G and add things for comoving integration */
240
241#endif
242
243#ifdef EXTERNALGRAVITY
244 gravity_external();
245#endif
246}
247
248void sim::gravity_set_oldacc(int timebin)
249{
250#ifdef HIERARCHICAL_GRAVITY
251 if(timebin == All.HighestOccupiedGravTimeBin)
252#endif
253 {
254 mpi_printf("GRAVTREE/FMM: Setting OldAcc!\n");
255
256 particle_data *P = Sp.P;
257
258 double ginv = 1 / All.G;
259
260 for(int idx = 0; idx < Sp.TimeBinsGravity.NActiveParticles; idx++)
261 {
262 int target = Sp.TimeBinsGravity.ActiveParticleList[idx];
263#if defined(PMGRID) && defined(PERIODIC) && !defined(TREEPM_NOTIMESPLIT)
264 double ax = P[target].GravAccel[0] + P[target].GravPM[0];
265 double ay = P[target].GravAccel[1] + P[target].GravPM[1];
266 double az = P[target].GravAccel[2] + P[target].GravPM[2];
267#else
268 double ax = P[target].GravAccel[0];
269 double ay = P[target].GravAccel[1];
270 double az = P[target].GravAccel[2];
271#endif
272 P[target].OldAcc = sqrt(ax * ax + ay * ay + az * az) * ginv;
273 }
274 }
275}
276
278{
279 particle_data *P = Sp.P;
280
281#ifndef PERIODIC
283 {
284 /* here we carry out an integration in comoving coordinates but in a non-periodic space (i.e. the 'big sphere setup') */
285 double fac = 0.5 * All.Hubble * All.Hubble * All.Omega0 / All.G;
286
287 for(int i = 0; i < Sp.TimeBinsGravity.NActiveParticles; i++)
288 {
289 int target = Sp.TimeBinsGravity.ActiveParticleList[i];
290
291 double pos[3];
292 Sp.intpos_to_pos(P[target].IntPos, pos); /* converts the integer distance to floating point */
293
294 for(int j = 0; j < 3; j++)
295 P[target].GravAccel[j] += fac * pos[j];
296
297#ifdef EVALPOTENTIAL
298 double r2 = 0;
299 for(int k = 0; k < 3; k++)
300 r2 += pos[k] * pos[k];
301 P[target].Potential -= fac * r2;
302#endif
303 }
304 }
305#endif
306
307 /* muliply by G */
308 for(int idx = 0; idx < Sp.TimeBinsGravity.NActiveParticles; idx++)
309 {
310 int target = Sp.TimeBinsGravity.ActiveParticleList[idx];
311
312 for(int j = 0; j < 3; j++)
313 P[target].GravAccel[j] *= All.G;
314
315#if defined(EVALPOTENTIAL) && defined(PMGRID) && defined(PERIODIC)
316 /* To get correct zero point in potential for TreePM calculation, need to add this term,
317 * because we cannot include the pi/(V*alpha^2) term in the correction potential in real space
318 * since we only touch a restricted set of particles in the tree calculation
319 */
320#ifndef GRAVITY_TALLBOX
321 /* note: for the tallbox, the constant potential term is comuted as part of the long-range force and not the short-range force,
322 * unlike in the ordinary periodic case */
323 if(GravTree.DoPM)
324 {
325 double alpha = 0.5 / Sp.Asmth[0];
326 P[target].Potential +=
327 All.TotalMass * M_PI / (alpha * alpha * All.BoxSize * All.BoxSize * All.BoxSize) * (LONG_X * LONG_Y * LONG_Z);
328 }
329#endif
330#endif
331
332#if defined(EVALPOTENTIAL)
333#ifndef FMM
334 /* remove self-interaction */
335#if NSOFTCLASSES > 1
336 P[target].Potential += P[target].getMass() / (All.ForceSoftening[P[target].getSofteningClass()] / 2.8);
337#else
338 P[target].Potential += P[target].getMass() / (All.ForceSoftening[0] / 2.8);
339#endif
340
341#endif
342
343#if defined(FMM) && defined(PERIODIC) && !defined(PMGRID)
344 /* in FMM case, add in interaction with other images in periodic grid */
345 P[target].Potential += P[target].getMass() * Ewald.ewald_gridlookup_origin_D0();
346#endif
347
348 P[target].Potential *= All.G;
349
350#if defined(PMGRID) && !defined(FORCETEST_TESTFORCELAW)
351 P[target].Potential += P[target].PM_Potential; /* add in long-range potential */
352#endif
353#endif
354
356 {
357#ifdef PERIODIC
358 Terminate(
359 "You specified a periodic simulation in physical coordinates but with a non-zero cosmological constant - this can't be "
360 "run");
361#endif
362 /* Finally, the following factor allows a computation of a cosmological simulation
363 with vacuum energy in physical coordinates */
364
365 double pos[3];
366 Sp.intpos_to_pos(P[target].IntPos, pos); /* converts the integer distance to floating point */
367
368 double fac = All.OmegaLambda * All.Hubble * All.Hubble;
369
370 for(int j = 0; j < 3; j++)
371 Sp.P[target].GravAccel[j] += fac * pos[j];
372
373#ifdef EVALPOTENTIAL
374 double r2 = 0;
375 for(int k = 0; k < 3; k++)
376 r2 += pos[k] * pos[k];
377 P[target].Potential -= 0.5 * fac * r2;
378#endif
379 }
380 }
381}
382
383#if defined(PMGRID) && (defined(TREEPM_NOTIMESPLIT) || defined(PLACEHIGHRESREGION))
384void sim::gravity_pm(int timebin)
385{
386 double tstart = Logs.second();
387 TIMER_START(CPU_PM_GRAVITY);
388 mpi_printf("TREEPM: Starting PM part of force calculation. (timebin=%d)\n", timebin);
389
390#if !defined(PERIODIC) || defined(PLACEHIGHRESREGION)
391 PM.pm_init_regionsize();
392#endif
393
394 if((GravTree.DoPM & TREE_DO_BASE_PM))
395 {
396#ifdef PERIODIC
397 PM.pmforce_periodic(LOW_MESH, NULL);
398#else
399 /* non periodic PM mesh */
400 PM.pmforce_nonperiodic(LOW_MESH);
401#endif
402 }
403
404#ifdef PLACEHIGHRESREGION
405 if((GravTree.DoPM & TREE_DO_HIGHRES_PM))
406 PM.pmforce_nonperiodic(HIGH_MESH);
407#endif
408
409 mpi_printf("TREEPM: Finished PM part of force calculation.\n");
410 TIMER_STOP(CPU_PM_GRAVITY);
411 double tend = Logs.second();
412 All.CPUForLastPMExecution = Logs.timediff(tstart, tend);
413}
414#endif
415
419#if defined(PMGRID) && defined(PERIODIC) && !defined(TREEPM_NOTIMESPLIT)
421{
422#ifndef SELFGRAVITY
423 return;
424#endif
425
426 double tstart = Logs.second();
427 TIMER_START(CPU_PM_GRAVITY);
428
429 for(int i = 0; i < Sp.NumPart; i++)
430 {
431 Sp.P[i].GravPM[0] = Sp.P[i].GravPM[1] = Sp.P[i].GravPM[2] = 0;
432#ifdef EVALPOTENTIAL
433 Sp.P[i].PM_Potential = 0;
434#endif
435 }
436
437 PM.pmforce_periodic(0, NULL);
438
439 /* multiply with the gravitational constant */
440 for(int i = 0; i < Sp.NumPart; i++)
441 {
442 for(int j = 0; j < 3; j++)
443 Sp.P[i].GravPM[j] *= All.G;
444#ifdef EVALPOTENTIAL
445 Sp.P[i].PM_Potential *= All.G;
446#endif
447 }
448
449 TIMER_STOP(CPU_PM_GRAVITY);
450 double tend = Logs.second();
451 All.CPUForLastPMExecution = Logs.timediff(tstart, tend);
452
453 Sp.find_long_range_step_constraint();
454}
455#endif
global_data_all_processes All
Definition: main.cc:40
char DoEwald
Definition: gravtree.h:361
void set_softenings(void)
This function sets the (comoving) softening length of all particle types in the table All....
char MeasureCostFlag
Definition: gravtree.h:291
void gravity_tree(int timebin)
This function computes the gravitational forces for all active particles.
Definition: gwalk.cc:600
void intpos_to_pos(MyIntPosType *intpos, T *posdiff)
double timediff(double t0, double t1)
Definition: logs.cc:488
double second(void)
Definition: logs.cc:471
void mpi_printf(const char *fmt,...)
Definition: setcomm.h:55
MPI_Comm Communicator
Definition: setcomm.h:31
void compute_grav_accelerations(int timebin)
This routine computes the gravitational accelerations for all active particles.
Definition: gravity.cc:41
domain< simparticles > Domain
Definition: simulation.h:58
void gravity_comoving_factors(int timebin)
Definition: gravity.cc:277
gwalk GravTree
Definition: simulation.h:65
void gravity(int timebin)
main driver routine of gravity tree/fmm force calculation
Definition: gravity.cc:171
void gravity_pm(int timebin)
void gravity_long_range_force(void)
sph NgbTree
Definition: simulation.h:68
void gravity_set_oldacc(int timebin)
Definition: gravity.cc:248
void endrun(void)
This function aborts the simulations.
Definition: begrun.cc:395
simparticles Sp
Definition: simulation.h:56
long long TotNumPart
Definition: simparticles.h:46
TimeBinData TimeBinsGravity
Definition: simparticles.h:205
particle_data * P
Definition: simparticles.h:54
void treeallocate(int max_partindex, partset *Pptr, domain< partset > *Dptr)
Definition: tree.cc:776
void treefree(void)
Definition: tree.cc:1232
int treebuild(int ninsert, int *indexlist)
Definition: tree.cc:52
#define LOW_MESH
Definition: constants.h:33
#define HIGH_MESH
Definition: constants.h:34
#define DIRECT_SUMMATION_THRESHOLD
Definition: constants.h:39
#define M_PI
Definition: constants.h:56
@ STANDARD
Definition: domain.h:23
#define LONG_Y
Definition: dtypes.h:369
#define LONG_X
Definition: dtypes.h:361
#define LONG_Z
Definition: dtypes.h:377
ewald Ewald
Definition: main.cc:42
#define TREE_ACTIVE_CUTTOFF_HIGHRES_PM
Definition: gravtree.h:24
#define TREE_DO_BASE_PM
Definition: gravtree.h:21
#define TREE_DO_HIGHRES_PM
Definition: gravtree.h:22
#define TREE_ACTIVE_CUTTOFF_BASE_PM
Definition: gravtree.h:23
logs Logs
Definition: main.cc:43
#define TIMER_START(counter)
Starts the timer counter.
Definition: logs.h:197
#define TIMER_STOP(counter)
Stops the timer counter.
Definition: logs.h:220
#define Terminate(...)
Definition: macros.h:15
void sumup_large_ints(int n, int *src, long long *res, MPI_Comm comm)
expr sqrt(half arg)
Definition: half.hpp:2777
int NActiveParticles
Definition: timestep.h:18
int * ActiveParticleList
Definition: timestep.h:20
long long GlobalNActiveParticles
Definition: timestep.h:19
integertime Ti_Current
Definition: allvars.h:188
double ForceSoftening[NSOFTCLASSES+NSOFTCLASSES_HYDRO+2]
Definition: allvars.h:258
MyDouble getMass(void)
unsigned char getSofteningClass(void)
vector< MyFloat > GravAccel
Definition: particle_data.h:55