GADGET-4
lightcone.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#ifdef LIGHTCONE
15
16#include <gsl/gsl_rng.h>
17#include <hdf5.h>
18#include <math.h>
19#include <mpi.h>
20#include <stdio.h>
21#include <stdlib.h>
22#include <string.h>
23#include <sys/stat.h>
24#include <sys/types.h>
25#include <unistd.h>
26
27#include "../data/allvars.h"
28#include "../data/dtypes.h"
29#include "../data/intposconvert.h"
30#include "../data/mymalloc.h"
31#include "../io/hdf5_util.h"
32#include "../lightcone/lightcone.h"
33#include "../lightcone/lightcone_massmap_io.h"
34#include "../main/main.h"
35#include "../main/simulation.h"
36#include "../sort/cxxsort.h"
37#include "../system/system.h"
38
39void lightcone::lightcone_init_intposconverter(double linklength)
40{
41 double maxlength = 0;
42
43#ifdef LIGHTCONE_PARTICLES
44 if(ConeGlobComDistStart > maxlength)
45 maxlength = ConeGlobComDistStart;
46#endif
47
48#ifdef LIGHTCONE_MASSMAPS
49 for(int i = 0; i < NumMassMapBoundaries; i++)
50 if(MassMapBoundariesComDist[i] > maxlength)
51 maxlength = MassMapBoundariesComDist[i];
52#endif
53
54 maxlength += linklength;
55
56#ifdef LIGHTCONE_PARTICLES
57 Lp->RegionLen = 2.0 * maxlength;
58 Lp->FacCoordToInt = pow(2.0, BITS_FOR_POSITIONS) / Lp->RegionLen;
59 Lp->FacIntToCoord = Lp->RegionLen / pow(2.0, BITS_FOR_POSITIONS);
60#endif
61}
62
63#ifdef LIGHTCONE_MASSMAPS
64int lightcone::lightcone_add_position_massmaps(particle_data *P, double *pos, double ascale)
65{
66 int buffer_full_flag = 0;
67
68 long ipnest;
69 vec2pix_ring(All.LightConeMassMapsNside, pos, &ipnest);
70
71 int task;
72 subdivide_evenly_get_bin(Mp->Npix, NTask, ipnest, &task);
73
74 if(Mp->NumPart >= Mp->MaxPart)
75 {
76 int new_maxpart = std::max<int>(Mp->NumPart + 1, (1 + ALLOC_TOLERANCE) * (1 + ALLOC_TOLERANCE) * Mp->MaxPart);
77 if(Mp->NumPart >= new_maxpart)
78 Terminate("Mp->NumPart >= new_maxpart");
79
80 Mp->reallocate_memory_maxpart(new_maxpart);
81
82 buffer_full_flag = 1;
83 }
84
85 int p = Mp->NumPart++;
86
87 Mp->P[p].Ascale = ascale;
88 Mp->P[p].setMass(P->getMass());
89
90 Mp->P[p].PixIndex = ipnest;
91 Mp->P[p].Task = task;
92
93 if(task < 0 || task >= NTask || ipnest < 0 || ipnest >= Mp->Npix)
94 Terminate("strange assignment: task=%d NTask=%d pos=(%g|%g|%g) ipnest=%d\n", task, NTask, pos[0], pos[1], pos[2], (int)ipnest);
95
96 return buffer_full_flag;
97}
98#endif
99
100#ifdef LIGHTCONE_PARTICLES
101void lightcone::lightcone_add_position_particles(particle_data *P, double *pos, double ascale)
102{
103 if(Lp->NumPart >= Lp->MaxPart)
104 {
105 int new_maxpart = (1 + ALLOC_TOLERANCE) * Lp->MaxPart;
106 if(Lp->NumPart >= new_maxpart)
107 Terminate("Lp->NumPart >= new_maxpart. Lp->NumPart=%d Lp->MaxPart=%d new_maxpart=%d", Lp->NumPart, Lp->MaxPart, new_maxpart);
108
109 Lp->reallocate_memory_maxpart(new_maxpart);
110 }
111
112 int q = Lp->NumPart++;
113
114 MyIntPosType intpos[3];
115
116 Lp->pos_to_signedintpos(pos, (MySignedIntPosType *)intpos);
117
118 for(int i = 0; i < 3; i++)
119 Lp->P[q].IntPos[i] = intpos[i];
120
121 Lp->P[q].Ascale = ascale;
122
123 /* we change to physical velocity here */
124 for(int i = 0; i < 3; i++)
125 Lp->P[q].Vel[i] = P->Vel[i] / ascale;
126
127#ifdef LIGHTCONE_OUTPUT_ACCELERATIONS
128 for(int i = 0; i < 3; i++)
129 Lp->P[q].GravAccel[i] = P->GravAccel[i];
130#endif
131
132 Lp->P[q].setType(P->getType());
133 Lp->P[q].setMass(P->getMass());
134 Lp->P[q].ID.set(P->ID.get());
135#if defined(LIGHTCONE_PARTICLES_GROUPS) && defined(FOF)
136 Lp->P[q].setSofteningClass(P->getSofteningClass());
137#endif
138
139 if(P->ID.is_previously_most_bound())
140 Lp->P[q].ID.mark_as_formerly_most_bound();
141
142#if defined(LIGHTCONE_PARTICLES_GROUPS) && defined(FOF)
143 Lp->P[q].setFlagSaveDistance();
144#endif
145}
146#endif
147
148int lightcone::lightcone_test_for_particle_addition(particle_data *P, integertime time0, integertime time1, double dt_drift)
149{
150 int buffer_full_flag = 0;
151
152 int flag = 0;
153
154#ifdef LIGHTCONE_PARTICLES
155 if(time0 < ConeGlobTime_end && time1 > ConeGlobTime_start) /* we have a potential overlap */
156 flag = 1;
157#endif
158
159#ifdef LIGHTCONE_MASSMAPS
160 if(time0 < MassMapBoundariesTime[NumMassMapBoundaries - 1] && time1 > MassMapBoundariesTime[0])
161 flag = 1;
162#endif
163
164 if(flag == 0)
165 return buffer_full_flag;
166
167 static integertime last_time0 = -1, last_time1 = -1;
168 static double last_R0 = -1, last_R1 = -1;
169
170 if(time0 != last_time0 || time1 != last_time1)
171 {
172 last_R0 = Driftfac.get_comoving_distance(time0);
173 last_R1 = Driftfac.get_comoving_distance(time1);
174
175 last_time0 = time0;
176 last_time1 = time1;
177 }
178
179 double R0 = last_R0;
180 double R1 = last_R1;
181
182 double R0_squared = R0 * R0;
183 double R1_squared = R1 * R1;
184
185 double R1prime = R1 - (R0 - R1);
186
187 double pos[3];
188 Sp->intpos_to_pos(P->IntPos, pos);
189
190#ifdef LIGHTCONE_PARTICLES
191 bool previously = P->ID.is_previously_most_bound();
192#endif
193
194 NumLastCheck = 0;
195
196 for(int n = 0; n < NumBoxes; n++)
197 {
198 if(R0 < BoxList[n].Rmin)
199 continue;
200
201 if(R1prime > BoxList[n].Rmax)
202 break;
203
204 NumLastCheck++;
205
206 int i = BoxList[n].i;
207 int j = BoxList[n].j;
208 int k = BoxList[n].k;
209
210 double PosA[3];
211
212 PosA[0] = pos[0] + i * All.BoxSize;
213 PosA[1] = pos[1] + j * All.BoxSize;
214 PosA[2] = pos[2] + k * All.BoxSize;
215
216 double rA2 = PosA[0] * PosA[0] + PosA[1] * PosA[1] + PosA[2] * PosA[2];
217
218 if(rA2 < R0_squared)
219 {
220 double PosB[3];
221 double diffBminusA[3];
222
223 for(int q = 0; q < 3; q++)
224 {
225 diffBminusA[q] = P->Vel[q] * dt_drift * Sp->FacIntToCoord;
226 PosB[q] = PosA[q] + diffBminusA[q];
227 }
228
229 double rB2 = PosB[0] * PosB[0] + PosB[1] * PosB[1] + PosB[2] * PosB[2];
230
231 if(rB2 > R1_squared)
232 {
233 /* ok, particle crossed the lightcone. Interpolate the coordinate of the crossing */
234
235 double dr2 = diffBminusA[0] * diffBminusA[0] + diffBminusA[1] * diffBminusA[1] + diffBminusA[2] * diffBminusA[2];
236
237 double a = pow(R1 - R0, 2) - dr2;
238 double b = 2 * R0 * (R1 - R0) - 2 * (PosA[0] * diffBminusA[0] + PosA[1] * diffBminusA[1] + PosA[2] * diffBminusA[2]);
239 double c = R0 * R0 - rA2;
240
241 double det = b * b - 4 * a * c;
242
243 if(det < 0)
244 Terminate(
245 "det=%g R0=%g R1=%g rA=%g rB=%g dr=%g dt_drift=%g dx=(%g|%g|%g) vel=(%g|%g|%g) "
246 "posA=(%g|%g|%g) posB=(%g|%g|%g)\n",
247 det, R0, R1, sqrt(rA2), sqrt(rB2), sqrt(dr2), dt_drift, P->Vel[0] * dt_drift * Sp->FacIntToCoord,
248 P->Vel[1] * dt_drift * Sp->FacIntToCoord, P->Vel[2] * dt_drift * Sp->FacIntToCoord, P->Vel[0], P->Vel[1],
249 P->Vel[2], PosA[0], PosA[1], PosA[2], PosB[0], PosB[1], PosB[2]);
250
251 double fac = (-b - sqrt(det)) / (2 * a);
252
253 vector<double> Pos;
254
255 for(int q = 0; q < 3; q++)
256 Pos[q] = PosA[q] + fac * diffBminusA[q];
257
258 double ascale = All.TimeBegin * exp((time0 + (time1 - time0) * fac) * All.Timebase_interval);
259
260 /* now we can add particle at position Pos[] to the lightcone, provided it fits into the angular mask */
261
262 if(fac < 0 || fac > 1)
263 {
264 warn(
265 "ascale=%g fac=%g fac-1%g R0=%g R1=%g rA=%g rB=%g dr=%g dt_drift=%g dx=(%g|%g|%g) vel=(%g|%g|%g) "
266 "posA=(%g|%g|%g) posB=(%g|%g|%g)\n",
267 ascale, fac, fac - 1, R0, R1, sqrt(rA2), sqrt(rB2), sqrt(dr2), dt_drift,
268 P->Vel[0] * dt_drift * Sp->FacIntToCoord, P->Vel[1] * dt_drift * Sp->FacIntToCoord,
269 P->Vel[2] * dt_drift * Sp->FacIntToCoord, P->Vel[0], P->Vel[1], P->Vel[2], PosA[0], PosA[1], PosA[2], PosB[0],
270 PosB[1], PosB[2]);
271 }
272 else
273 {
274#ifdef LIGHTCONE_PARTICLES
275 if(ascale >= ConeGlobAstart && ascale < ConeGlobAend)
276 for(int cone = 0; cone < Nlightcones; cone++)
277 if(lightcone_is_cone_member_basic(ascale, Pos, previously, cone))
278 {
279 /* we only add the particle once if it is at least contained in one of the cones */
280 lightcone_add_position_particles(P, Pos.da, ascale);
281 break;
282 }
283#endif
284
285#ifdef LIGHTCONE_MASSMAPS
286 if(ascale >= MassMapBoundariesAscale[0] && ascale < MassMapBoundariesAscale[NumMassMapBoundaries - 1])
287 buffer_full_flag |= lightcone_add_position_massmaps(P, Pos.da, ascale);
288#endif
289 }
290 }
291 }
292 }
293
294 return buffer_full_flag;
295}
296
297#ifdef LIGHTCONE_PARTICLES
298
299bool lightcone::lightcone_is_cone_member(int i, int cone)
300{
301#if defined(LIGHTCONE_PARTICLES_GROUPS) && defined(FOF)
302 if(!Lp->P[i].getFlagSaveDistance())
303 return false;
304#endif
305
306 vector<double> pos;
307
308 if(i >= Lp->NumPart)
309 Terminate("i=%d Lp->NumPart=%d\n", i, Lp->NumPart);
310
311 Lp->signedintpos_to_pos((MySignedIntPosType *)Lp->P[i].IntPos, pos.da);
312
313 return lightcone_is_cone_member_basic(Lp->P[i].Ascale, pos, Lp->P[i].ID.is_previously_most_bound(), cone);
314}
315
316bool lightcone::lightcone_is_cone_member_basic(double ascale, vector<double> &pos, bool previously, int cone)
317{
318 if(ascale < Cones[cone].Astart || ascale > Cones[cone].Aend)
319 return false;
320
321 if(Cones[cone].OnlyMostBoundFlag == 1 && previously == false)
322 return false;
323
324 if(Cones[cone].LcType == LC_TYPE_FULLSKY)
325 {
326 return true;
327 }
328 else if(Cones[cone].LcType == LC_TYPE_OCTANT)
329 {
330 double rad = pos.norm();
331 double phi = atan2(pos[1], pos[0]);
332 double theta = acos(pos[2] / rad);
333
334 if(phi < 0)
335 phi += 2 * M_PI;
336
337 int octnr = phi / (0.5 * M_PI);
338
339 if(octnr >= 4)
340 octnr = 3;
341
342 if(octnr < 0)
343 octnr = 0;
344
345 if(theta > 0.5 * M_PI)
346 octnr += 4;
347
348 if(octnr == Cones[cone].OctantNr)
349 return true;
350 else
351 return false;
352 }
353 else if(Cones[cone].LcType == LC_TYPE_PENCIL)
354 {
355 double rad = pos.norm();
356
357 // note: PencilDir is already normalized
358
359 double angle = acos((pos * Cones[cone].PencilDir) / rad);
360
361 if(angle < Cones[cone].PencilAngleRad)
362 return true;
363 else
364 return false;
365 }
366 else if(Cones[cone].LcType == LC_TYPE_DISK)
367 {
368 double dist = pos * Cones[cone].DiskNormal;
369
370 if(fabs(dist) < 0.5 * Cones[cone].DiskThickness)
371 return true;
372 else
373 return false;
374 }
375 else if(Cones[cone].LcType == LC_TYPE_SQUAREMAP)
376 {
377 double x = pos * Cones[cone].SquareMapXdir;
378 double y = pos * Cones[cone].SquareMapYdir;
379 double z = pos * Cones[cone].SquareMapZdir;
380
381 /* not angle has been converted from degrees to rad */
382 if(z > 0 && fabs(x) < Cones[cone].SquareMapAngleRad * z && fabs(y) < Cones[cone].SquareMapAngleRad * z)
383 return true;
384 else
385 return false;
386 }
387 return false;
388}
389
390void lightcone::lightcone_init_geometry(char *fname)
391{
392 if(!All.ComovingIntegrationOn)
393 Terminate("LIGHTCONE_PARTICLES: makes only sense for cosmological simulations with ComovingIntegrationOn enabled\n");
394
395 if(ThisTask == 0)
396 {
397 for(int iter = 0; iter < 2; iter++)
398 {
399 Nlightcones = 0;
400 FILE *fd;
401
402 if(!(fd = fopen(fname, "r")))
403 Terminate("LIGHTCONE_PARTICLES: cannot read lightcone geometries in file `%s'\n", fname);
404
405 if(iter == 0)
406 {
407 while(1)
408 {
409 int lc_type;
410 if(fscanf(fd, "%d", &lc_type) != 1)
411 break;
412
413 double dummy;
414
415 switch(lc_type)
416 {
417 case LC_TYPE_FULLSKY:
418 if(fscanf(fd, "%lg %lg %lg", &dummy, &dummy, &dummy) != 3)
419 Terminate("LIGHTCONE_PARTICLES: incorrect data for lightcone type %d in file '%s'", lc_type, fname);
420 break;
421
422 case LC_TYPE_OCTANT:
423 if(fscanf(fd, "%lg %lg %lg %lg", &dummy, &dummy, &dummy, &dummy) != 4)
424 Terminate("LIGHTCONE_PARTICLES: incorrect data for lightcone type %d in file '%s'", lc_type, fname);
425 break;
426
427 case LC_TYPE_PENCIL:
428 if(fscanf(fd, "%lg %lg %lg %lg %lg %lg %lg", &dummy, &dummy, &dummy, &dummy, &dummy, &dummy, &dummy) != 7)
429 Terminate("LIGHTCONE_PARTICLES: incorrect data for lightcone type %d in file '%s'", lc_type, fname);
430 break;
431
432 case LC_TYPE_DISK:
433 if(fscanf(fd, "%lg %lg %lg %lg %lg %lg %lg", &dummy, &dummy, &dummy, &dummy, &dummy, &dummy, &dummy) != 7)
434 Terminate("LIGHTCONE_PARTICLES: incorrect data for lightcone type %d in file '%s'", lc_type, fname);
435 break;
436
437 case LC_TYPE_SQUAREMAP:
438 if(fscanf(fd, "%lg %lg %lg %lg %lg %lg %lg %lg %lg %lg", &dummy, &dummy, &dummy, &dummy, &dummy, &dummy,
439 &dummy, &dummy, &dummy, &dummy) != 10)
440 Terminate("LIGHTCONE_PARTICLES: incorrect data for squaremap type %d in file '%s'", lc_type, fname);
441 break;
442
443 default:
444 Terminate("LIGHTCONE_PARTICLES: unknown lightcone type %d in file '%s'", lc_type, fname);
445 break;
446 }
447
448 Nlightcones++;
449 }
450
451 Cones = (cone_data *)Mem.mymalloc("Cones", Nlightcones * sizeof(cone_data));
452
453 mpi_printf("LIGHTCONE_PARTICLES: read definitions with %d entries in file `%s'.\n", Nlightcones, fname);
454 }
455 else if(iter == 1)
456 {
457 while(1)
458 {
459 int lc_type;
460 if(fscanf(fd, "%d", &lc_type) != 1)
461 break;
462
463 Cones[Nlightcones].LcType = lc_type;
464
465 switch(lc_type)
466 {
467 case LC_TYPE_FULLSKY:
468 fscanf(fd, "%d %lg %lg", &Cones[Nlightcones].OnlyMostBoundFlag, &Cones[Nlightcones].Astart,
469 &Cones[Nlightcones].Aend);
470 sprintf(Cones[Nlightcones].Tag, "Full-sky");
471 break;
472
473 case LC_TYPE_OCTANT:
474 fscanf(fd, "%d %lg %lg", &Cones[Nlightcones].OnlyMostBoundFlag, &Cones[Nlightcones].Astart,
475 &Cones[Nlightcones].Aend);
476 fscanf(fd, "%d", &Cones[Nlightcones].OctantNr);
477 sprintf(Cones[Nlightcones].Tag, "Octant");
478 break;
479
480 case LC_TYPE_PENCIL:
481 fscanf(fd, "%d %lg %lg", &Cones[Nlightcones].OnlyMostBoundFlag, &Cones[Nlightcones].Astart,
482 &Cones[Nlightcones].Aend);
483 fscanf(fd, "%lg %lg %lg", &Cones[Nlightcones].PencilDir[0], &Cones[Nlightcones].PencilDir[1],
484 &Cones[Nlightcones].PencilDir[2]);
485 fscanf(fd, "%lg", &Cones[Nlightcones].PencilAngle);
486
487 /* normalize the normal vector in case it is not normalized yet */
488 Cones[Nlightcones].PencilDir *= 1.0 / Cones[Nlightcones].PencilDir.norm();
489
490 /* convert to rad */
491 Cones[Nlightcones].PencilAngleRad = Cones[Nlightcones].PencilAngle * M_PI / 180.0;
492
493 sprintf(Cones[Nlightcones].Tag, "Pencil-Beam");
494 break;
495
496 case LC_TYPE_DISK:
497 fscanf(fd, "%d %lg %lg", &Cones[Nlightcones].OnlyMostBoundFlag, &Cones[Nlightcones].Astart,
498 &Cones[Nlightcones].Aend);
499 fscanf(fd, "%lg %lg %lg", &Cones[Nlightcones].DiskNormal[0], &Cones[Nlightcones].DiskNormal[1],
500 &Cones[Nlightcones].DiskNormal[2]);
501 fscanf(fd, "%lg", &Cones[Nlightcones].DiskThickness);
502
503 /* normalize the normal vector in case it is not normalized yet */
504 Cones[Nlightcones].DiskNormal *= 1.0 / Cones[Nlightcones].DiskNormal.norm();
505 sprintf(Cones[Nlightcones].Tag, "Disk (for image)");
506 break;
507
508 case LC_TYPE_SQUAREMAP:
509 fscanf(fd, "%d %lg %lg", &Cones[Nlightcones].OnlyMostBoundFlag, &Cones[Nlightcones].Astart,
510 &Cones[Nlightcones].Aend);
511 fscanf(fd, "%lg %lg %lg", &Cones[Nlightcones].SquareMapZdir[0], &Cones[Nlightcones].SquareMapZdir[1],
512 &Cones[Nlightcones].SquareMapZdir[2]);
513 fscanf(fd, "%lg %lg %lg", &Cones[Nlightcones].SquareMapXdir[0], &Cones[Nlightcones].SquareMapXdir[1],
514 &Cones[Nlightcones].SquareMapXdir[2]);
515 fscanf(fd, "%lg", &Cones[Nlightcones].SquareMapAngle);
516
517 /* establish coordinate system */
518
519 Cones[Nlightcones].SquareMapZdir *= 1.0 / Cones[Nlightcones].SquareMapZdir.norm();
520
521 // cross product
522 Cones[Nlightcones].SquareMapYdir = Cones[Nlightcones].SquareMapZdir ^ Cones[Nlightcones].SquareMapXdir;
523
524 Cones[Nlightcones].SquareMapYdir *= 1.0 / Cones[Nlightcones].SquareMapYdir.norm();
525
526 // cross product
527 Cones[Nlightcones].SquareMapXdir = Cones[Nlightcones].SquareMapYdir ^ Cones[Nlightcones].SquareMapZdir;
528
529 Cones[Nlightcones].SquareMapAngleRad = Cones[Nlightcones].SquareMapAngle * M_PI / 180.0;
530
531 mpi_printf("LIGHTCONE_SQUAREMAP: cone=%2d x-axis = %15g %15g %15g\n", Nlightcones,
532 Cones[Nlightcones].SquareMapXdir[0], Cones[Nlightcones].SquareMapXdir[1],
533 Cones[Nlightcones].SquareMapXdir[2]);
534 mpi_printf("LIGHTCONE_SQUAREMAP: cone=%2d y-axis = %15g %15g %15g\n", Nlightcones,
535 Cones[Nlightcones].SquareMapYdir[0], Cones[Nlightcones].SquareMapYdir[1],
536 Cones[Nlightcones].SquareMapYdir[2]);
537 mpi_printf("LIGHTCONE_SQUAREMAP: cone=%2d z-axis = %15g %15g %15g\n", Nlightcones,
538 Cones[Nlightcones].SquareMapZdir[0], Cones[Nlightcones].SquareMapZdir[1],
539 Cones[Nlightcones].SquareMapZdir[2]);
540 sprintf(Cones[Nlightcones].Tag, "Square-map");
541 break;
542 default:
543 Terminate("odd");
544 }
545
546 Nlightcones++;
547 }
548 }
549 fclose(fd);
550 }
551
552 for(int i = 0; i < Nlightcones; i++)
553 mpi_printf("LIGHTCONE_PARTICLES: lightcone #%2d: %18s %20s Astart=%10g Aend=%10g\n", i, Cones[i].Tag,
554 Cones[i].OnlyMostBoundFlag ? "(only most bound)" : "(all particles)", Cones[i].Astart, Cones[i].Aend);
555
556 mpi_printf("\n");
557 }
558
559 /* now tell also all other ranks about the lightcones */
560
561 MPI_Bcast(&Nlightcones, 1, MPI_INT, 0, Communicator);
562
563 if(ThisTask != 0)
564 Cones = (cone_data *)Mem.mymalloc("Cones", Nlightcones * sizeof(cone_data));
565
566 MPI_Bcast(Cones, Nlightcones * sizeof(cone_data), MPI_BYTE, 0, Communicator);
567}
568
569int lightcone::lightcone_init_times(void)
570{
571 for(int i = 0; i < Nlightcones; i++)
572 {
573 Cones[i].Time_start = log(Cones[i].Astart / All.TimeBegin) / All.Timebase_interval;
574 Cones[i].Time_end = log(Cones[i].Aend / All.TimeBegin) / All.Timebase_interval;
575
576 Cones[i].ComDistStart = Driftfac.get_comoving_distance(Cones[i].Time_start);
577 Cones[i].ComDistEnd = Driftfac.get_comoving_distance(Cones[i].Time_end);
578 }
579
580 double astart = 1.0e10;
581 double aend = 0;
582
583 for(int i = 0; i < Nlightcones; i++)
584 {
585 if(Cones[i].Astart < astart)
586 astart = Cones[i].Astart;
587
588 if(Cones[i].Aend > aend)
589 aend = Cones[i].Aend;
590 }
591
592 ConeGlobAstart = astart;
593 ConeGlobAend = aend;
594
595 ConeGlobTime_start = log(ConeGlobAstart / All.TimeBegin) / All.Timebase_interval;
596 ConeGlobTime_end = log(ConeGlobAend / All.TimeBegin) / All.Timebase_interval;
597
598 ConeGlobComDistStart = Driftfac.get_comoving_distance(ConeGlobTime_start);
599 ConeGlobComDistEnd = Driftfac.get_comoving_distance(ConeGlobTime_end);
600
601 int n = ceil(ConeGlobComDistStart / All.BoxSize + 1);
602
603 for(int rep = 0; rep < 2; rep++)
604 {
605 if(rep == 1)
606 // BoxList = (boxlist *) Mem.mymalloc_movable(&BoxList, "BoxList", NumBoxes * sizeof(boxlist));
607 BoxList = Mem.alloc_movable<boxlist> MMM(BoxList, NumBoxes);
608
609 NumBoxes = 0;
610
611 for(int i = -n; i <= n; i++)
612 for(int j = -n; j <= n; j++)
613 for(int k = -n; k <= n; k++)
614 {
615 double corner[3];
616
617 corner[0] = i * All.BoxSize;
618 corner[1] = j * All.BoxSize;
619 corner[2] = k * All.BoxSize;
620
621 double Rmin, Rmax;
622
623 if(lightcone_box_at_corner_overlaps_at_least_with_one_cone(corner, Rmin, Rmax))
624 {
625 if(rep == 1)
626 {
627 BoxList[NumBoxes].i = i;
628 BoxList[NumBoxes].j = j;
629 BoxList[NumBoxes].k = k;
630 BoxList[NumBoxes].Rmin = Rmin;
631 BoxList[NumBoxes].Rmax = Rmax;
632 }
633
634 NumBoxes++;
635 }
636 }
637 }
638
639 mycxxsort(BoxList, BoxList + NumBoxes, lightcone_compare_BoxList_Rmax);
640
641 lightcone_clear_boxlist(All.Time);
642
643 NumLastCheck = 0;
644
645 double fac = (4 * M_PI / 3.0) * pow(ConeGlobComDistStart, 3) / pow(All.BoxSize, 3);
646
647 mpi_printf(
648 "LIGHTCONE_PARTICLES: scale_factor: %10g to %10g comoving distance: %10g to %10g covered volume in units of box "
649 "volume=%g\n",
650 ConeGlobAstart, ConeGlobAend, ConeGlobComDistStart, ConeGlobComDistEnd, fac);
651
652 mpi_printf("LIGHTCONE_PARTICLES: number of box replicas to check for this lightcone geometry settings = %d\n", NumBoxes);
653
654 if(NumBoxes > LIGHTCONE_MAX_BOXREPLICAS)
655 {
656 mpi_printf(
657 "\nLIGHTCONE_PARTICLES: Your lightcone extends to such high redshift that the box needs to be replicated a huge number of "
658 "times to cover it,\n"
659 "more than the prescribed limit of LIGHTCONE_MAX_BOXREPLICAS=%d. We better don't do such an inefficient run, unless you "
660 "override this constant.\n",
661 LIGHTCONE_MAX_BOXREPLICAS);
662 return 1;
663 }
664
665 return 0;
666}
667
668void lightcone::lightcone_clear_boxlist(double ascale)
669{
670 double time_start = log(ascale / All.TimeBegin) / All.Timebase_interval;
671
672 double dist = Driftfac.get_comoving_distance(time_start);
673
674 int count = 0;
675
676 for(int i = 0; i < NumBoxes; i++)
677 {
678 if(dist < BoxList[i].Rmin)
679 {
680 BoxList[i] = BoxList[--NumBoxes];
681 i--;
682 count++;
683 }
684 }
685
686 if(count)
687 {
688 mpi_printf("LIGHTCONE: Eliminated %d entries from BoxList\n", count);
689 mycxxsort(BoxList, BoxList + NumBoxes, lightcone_compare_BoxList_Rmax);
690 }
691}
692
693bool lightcone::lightcone_box_at_corner_overlaps_at_least_with_one_cone(double *corner, double &Rmin, double &Rmax)
694{
695 Rmin = MAX_DOUBLE_NUMBER;
696 Rmax = 0;
697
698 for(int ii = 0; ii <= 1; ii++)
699 for(int jj = 0; jj <= 1; jj++)
700 for(int kk = 0; kk <= 1; kk++)
701 {
702 double crn[3];
703 crn[0] = corner[0] + ii * All.BoxSize;
704 crn[1] = corner[1] + jj * All.BoxSize;
705 crn[2] = corner[2] + kk * All.BoxSize;
706
707 double r = sqrt(crn[0] * crn[0] + crn[1] * crn[1] + crn[2] * crn[2]);
708 if(Rmin > r)
709 Rmin = r;
710 if(Rmax < r)
711 Rmax = r;
712 }
713
714 for(int cone = 0; cone < Nlightcones; cone++)
715 {
716 if(Rmin < Cones[cone].ComDistStart && Rmax > Cones[cone].ComDistEnd)
717 {
718 if(Cones[cone].LcType == LC_TYPE_FULLSKY)
719 {
720 return true;
721 }
722 else if(Cones[cone].LcType == LC_TYPE_OCTANT)
723 {
724 return true; /* still need to make this more selective */
725 }
726 else if(Cones[cone].LcType == LC_TYPE_PENCIL)
727 {
728 return true; /* still need to make this more selective */
729 }
730 else if(Cones[cone].LcType == LC_TYPE_DISK)
731 {
732 double mindist = MAX_DOUBLE_NUMBER;
733 double first_dist = 0;
734
735 for(int ii = 0; ii <= 1; ii++)
736 for(int jj = 0; jj <= 1; jj++)
737 for(int kk = 0; kk <= 1; kk++)
738 {
739 double crn[3];
740 crn[0] = corner[0] + ii * All.BoxSize;
741 crn[1] = corner[1] + jj * All.BoxSize;
742 crn[2] = corner[2] + kk * All.BoxSize;
743
744 double dist =
745 crn[0] * Cones[cone].DiskNormal[0] + crn[1] * Cones[cone].DiskNormal[1] + crn[2] * Cones[cone].DiskNormal[2];
746
747 if(ii == 0 && jj == 0 && kk == 0) // upon first iteration
748 first_dist = dist;
749 else
750 {
751 if(first_dist * dist < 0) // points on opposite side imply overlap
752 return true;
753 }
754
755 if(fabs(dist) < mindist)
756 mindist = fabs(dist);
757
758 if(mindist < 0.5 * Cones[cone].DiskThickness)
759 return true;
760 }
761 }
762 }
763 }
764
765 return false;
766}
767
768#endif
769
770#ifdef LIGHTCONE_MASSMAPS
771
772void lightcone::lightcone_init_massmaps(void)
773{
774 if(!All.ComovingIntegrationOn)
775 Terminate("LIGHTCONE_MASSMAPS: Makes only sense for cosmological simulations with ComovingIntegrationOn enabled\n");
776
777 Mp->Npix = nside2npix(All.LightConeMassMapsNside);
778
779 subdivide_evenly(Mp->Npix, NTask, ThisTask, &Mp->FirstPix, &Mp->NpixLoc);
780
781 if(ThisTask == 0)
782 {
783 for(int iter = 0; iter < 2; iter++)
784 {
785 NumMassMapBoundaries = 0;
786
787 if(iter == 0)
788 {
789 double zend = 0;
790 double com_dist_end = 0;
791 NumMassMapBoundaries++;
792
793 while(zend <= All.LightConeMassMapMaxRedshift)
794 {
795 com_dist_end += All.LightConeMassMapThickness;
796
797 double aend = Driftfac.get_scalefactor_for_comoving_distance(com_dist_end);
798
799 zend = 1 / aend - 1;
800
801 NumMassMapBoundaries++;
802 }
803
804 MassMapBoundariesAscale = (double *)Mem.mymalloc_movable(&MassMapBoundariesAscale, "MassMapBoundariesAscale",
805 NumMassMapBoundaries * sizeof(double));
806
807 mpi_printf("LIGHTCONE_MASSMAPS: %d entries\n", NumMassMapBoundaries);
808
809 if(NumMassMapBoundaries < 2)
810 Terminate("Less than two boundaries detected");
811 }
812 else if(iter == 1)
813 {
814 double zend = 0;
815 double com_dist_end = 0;
816 MassMapBoundariesAscale[NumMassMapBoundaries] = 1.0;
817 NumMassMapBoundaries++;
818
819 while(zend <= All.LightConeMassMapMaxRedshift)
820 {
821 com_dist_end += All.LightConeMassMapThickness;
822
823 double aend = Driftfac.get_scalefactor_for_comoving_distance(com_dist_end);
824
825 MassMapBoundariesAscale[NumMassMapBoundaries] = aend;
826
827 zend = 1 / aend - 1;
828
829 NumMassMapBoundaries++;
830 }
831 }
832 }
833
834 mycxxsort(MassMapBoundariesAscale, MassMapBoundariesAscale + NumMassMapBoundaries, compare_doubles);
835 }
836
837 /* now tell also all other ranks about the lightcones */
838
839 MPI_Bcast(&NumMassMapBoundaries, 1, MPI_INT, 0, Communicator);
840
841 if(ThisTask != 0)
842 MassMapBoundariesAscale =
843 (double *)Mem.mymalloc_movable(&MassMapBoundariesAscale, "MassMapBoundariesAscale", NumMassMapBoundaries * sizeof(double));
844
845 MPI_Bcast(MassMapBoundariesAscale, NumMassMapBoundaries * sizeof(double), MPI_BYTE, 0, Communicator);
846
847 MassMapBoundariesTime =
848 (integertime *)Mem.mymalloc_movable(&MassMapBoundariesTime, "MassMapBoundariesTime", NumMassMapBoundaries * sizeof(integertime));
849 MassMapBoundariesComDist =
850 (double *)Mem.mymalloc_movable(&MassMapBoundariesComDist, "MassMapBoundariesComDist", NumMassMapBoundaries * sizeof(double));
851}
852
853int lightcone::lightcone_massmap_report_boundaries(void)
854{
855 double fac_max = 0;
856 const double allowed_fac_max = LIGHTCONE_MAX_BOXREPLICAS;
857
858 for(int i = 0; i < NumMassMapBoundaries; i++)
859 {
860 MassMapBoundariesTime[i] = log(MassMapBoundariesAscale[i] / All.TimeBegin) / All.Timebase_interval;
861 MassMapBoundariesComDist[i] = Driftfac.get_comoving_distance(MassMapBoundariesTime[i]);
862 }
863
864 for(int i = 0; i < NumMassMapBoundaries - 1; i++)
865 {
866 MassMapBoundariesTime[i] = log(MassMapBoundariesAscale[i] / All.TimeBegin) / All.Timebase_interval;
867 MassMapBoundariesComDist[i] = Driftfac.get_comoving_distance(MassMapBoundariesTime[i]);
868
869 double fac = (4 * M_PI / 3.0) * pow(MassMapBoundariesComDist[i], 3) / pow(All.BoxSize, 3);
870 double shell = fac - (4 * M_PI / 3.0) * pow(MassMapBoundariesComDist[i + 1], 3) / pow(All.BoxSize, 3);
871
872 mpi_printf(
873 "LIGHTCONE_MASSMAPS: entry=%3d scale_factor=%10.6f redshift=%10.6f comoving distance=%12.3f shell-volume in "
874 "units of box "
875 "volume=%g\n",
876 i, MassMapBoundariesAscale[i], 1 / MassMapBoundariesAscale[i] - 1, MassMapBoundariesComDist[i], shell);
877
878 if(fac > fac_max)
879 fac_max = fac;
880 }
881
882 if(fac_max > allowed_fac_max)
883 {
884 mpi_printf(
885 "\nLIGHTCONE_MASSMAPS: Your lightcone mass maps extend to such high redshift that the box needs to be replicated a huge "
886 "number of times to cover it (volume ratio %g). We better don't do such an inefficient run.\n"
887 "Setting LIGHTCONE_MAX_BOXREPLICAS=%g to a higher value (at least %g) can override this, however.\n",
888 fac_max, (double)LIGHTCONE_MAX_BOXREPLICAS, fac_max);
889 return 1;
890 }
891
892 return 0;
893}
894
895void lightcone::lightcone_massmap_flush(int dump_allowed_flag)
896{
897 lightcone_massmap_binning();
898
899 if(dump_allowed_flag)
900 {
901 while(All.CurrentMassMapBoundary < NumMassMapBoundaries - 1 &&
902 (All.Time >= MassMapBoundariesAscale[All.CurrentMassMapBoundary + 1] || All.Ti_Current >= TIMEBASE))
903 {
904 lightcone_massmap_io Lcone(Mp, this, Communicator, All.SnapFormat); /* get an I/O object */
905 Lcone.lightcone_massmap_save(All.CurrentMassMapBoundary++);
906
907 lightcone_massmap_binning();
908 }
909 }
910}
911
912void lightcone::lightcone_massmap_binning(void)
913{
914 double t0 = Logs.second();
915
916 int *Send_count = (int *)Mem.mymalloc_movable(&Send_count, "Send_count", sizeof(int) * NTask);
917 int *Send_offset = (int *)Mem.mymalloc_movable(&Send_offset, "Send_offset", sizeof(int) * NTask);
918 int *Recv_count = (int *)Mem.mymalloc_movable(&Recv_count, "Recv_count", sizeof(int) * NTask);
919 int *Recv_offset = (int *)Mem.mymalloc_movable(&Recv_offset, "Recv_offset", sizeof(int) * NTask);
920
921 /* count how many we have of each task */
922 for(int i = 0; i < NTask; i++)
923 Send_count[i] = 0;
924
925 for(int i = 0; i < Mp->NumPart; i++)
926 {
927 int target = Mp->P[i].Task;
928
929 if(target < 0 || target >= NTask)
930 Terminate("i=%d: target=%d target < 0 || target >= NTask", i, target);
931
932 if(target != ThisTask)
933 Send_count[target]++;
934 }
935
936 myMPI_Alltoall(Send_count, 1, MPI_INT, Recv_count, 1, MPI_INT, Communicator);
937
938 Recv_offset[0] = Send_offset[0] = 0;
939 int nexport = 0, nimport = 0;
940
941 for(int i = 0; i < NTask; i++)
942 {
943 nimport += Recv_count[i];
944 nexport += Send_count[i];
945
946 if(i > 0)
947 {
948 Send_offset[i] = Send_offset[i - 1] + Send_count[i - 1];
949 Recv_offset[i] = Recv_offset[i - 1] + Recv_count[i - 1];
950 }
951 }
952
953 lightcone_massmap_data *send_P =
954 (lightcone_massmap_data *)Mem.mymalloc_movable(&send_P, "send_P", nexport * sizeof(lightcone_massmap_data));
955
956 for(int i = 0; i < NTask; i++)
957 Send_count[i] = 0;
958
959 for(int i = 0; i < Mp->NumPart; i++)
960 {
961 int target = Mp->P[i].Task;
962
963 if(target != ThisTask)
964 {
965 send_P[Send_offset[target] + Send_count[target]] = Mp->P[i];
966 Send_count[target]++;
967
968 Mp->P[i] = Mp->P[Mp->NumPart - 1];
969 Mp->NumPart--;
970 i--;
971 }
972 }
973
974 if(Mp->NumPart + nimport > Mp->MaxPart)
975 Mp->reallocate_memory_maxpart(Mp->NumPart + nimport);
976
977 for(int ngrp = 1; ngrp < (1 << PTask); ngrp++)
978 {
979 int recvTask = ThisTask ^ ngrp;
980
981 if(recvTask < NTask)
982 if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
983 myMPI_Sendrecv(&send_P[Send_offset[recvTask]], Send_count[recvTask] * sizeof(lightcone_massmap_data), MPI_BYTE, recvTask,
984 TAG_DENS_A, &Mp->P[Mp->NumPart + Recv_offset[recvTask]],
985 Recv_count[recvTask] * sizeof(lightcone_massmap_data), MPI_BYTE, recvTask, TAG_DENS_A, Communicator,
986 MPI_STATUS_IGNORE);
987 }
988
989 Mp->NumPart += nimport;
990
991 Mem.myfree_movable(send_P);
992
993 Mem.myfree(Recv_offset);
994 Mem.myfree(Recv_count);
995 Mem.myfree(Send_offset);
996 Mem.myfree(Send_count);
997
998 mpi_printf("LIGHTCONE_MASSMAPS: before binning %9d (maxpart=%9d, memory for buffer %g MB)\n", Mp->NumPart, Mp->MaxPart,
999 (double)Mp->MaxPart * sizeof(lightcone_massmap_data) * TO_MBYTE_FAC);
1000
1001 int expunge = 0;
1002
1003 /* now bin on the map */
1004 for(int i = 0; i < Mp->NumPart; i++)
1005 {
1006 if(Mp->P[i].Task != ThisTask)
1007 Terminate("can't be");
1008
1009 if(Mp->P[i].Ascale >= MassMapBoundariesAscale[All.CurrentMassMapBoundary] &&
1010 Mp->P[i].Ascale < MassMapBoundariesAscale[All.CurrentMassMapBoundary + 1])
1011 {
1012 int pix = Mp->P[i].PixIndex - Mp->FirstPix;
1013 if(pix < 0 || pix >= Mp->NpixLoc)
1014 Terminate("wrong pixel");
1015
1016 MassMap[pix] += Mp->P[i].getMass();
1017
1018 Mp->P[i] = Mp->P[Mp->NumPart - 1];
1019 Mp->NumPart--;
1020 i--;
1021 }
1022 else if(Mp->P[i].Ascale < MassMapBoundariesAscale[All.CurrentMassMapBoundary])
1023 {
1024 Mp->P[i] = Mp->P[Mp->NumPart - 1];
1025 Mp->NumPart--;
1026 i--;
1027
1028 expunge++;
1029 }
1030 }
1031
1032 if(Mp->MaxPart > LIGHTCONE_MASSMAP_ALLOC_FAC * (Sp->TotNumPart / NTask) &&
1033 Mp->NumPart < LIGHTCONE_MAX_FILLFACTOR * LIGHTCONE_MASSMAP_ALLOC_FAC * (Sp->TotNumPart / NTask))
1034 Mp->reallocate_memory_maxpart(LIGHTCONE_MASSMAP_ALLOC_FAC * (Sp->TotNumPart / NTask));
1035
1036 double t1 = Logs.second();
1037
1038 mpi_printf("LIGHTCONE_MASSMAPS: after binning %9d (maxpart=%9d, memory for buffer %g MB) took=%g sec expunge=%d\n", Mp->NumPart,
1039 Mp->MaxPart, (double)Mp->MaxPart * sizeof(lightcone_massmap_data) * TO_MBYTE_FAC, Logs.timediff(t0, t1), expunge);
1040}
1041
1042#endif
1043
1044#endif
All
global_data_all_processes All
Definition: main.cc:40
MyIDStorage::get
MyIDType get(void) const
Definition: idstorage.h:87
MyIDStorage::is_previously_most_bound
bool is_previously_most_bound(void)
Definition: idstorage.h:117
driftfac::get_scalefactor_for_comoving_distance
double get_scalefactor_for_comoving_distance(double dist)
Definition: driftfac.cc:209
driftfac::get_comoving_distance
double get_comoving_distance(integertime time0)
Definition: driftfac.cc:191
logs::timediff
double timediff(double t0, double t1)
Definition: logs.cc:488
logs::second
double second(void)
Definition: logs.cc:471
memory::alloc_movable
T * alloc_movable(T *&ptr, const char *varname, size_t n, const char *func, const char *file, int linenr)
Definition: mymalloc.h:80
vector::norm
double norm(void)
Definition: symtensors.h:189
vector::da
T da[3]
Definition: symtensors.h:106
ALLOC_TOLERANCE
#define ALLOC_TOLERANCE
Definition: constants.h:74
LIGHTCONE_MAX_FILLFACTOR
#define LIGHTCONE_MAX_FILLFACTOR
Definition: constants.h:70
TO_MBYTE_FAC
#define TO_MBYTE_FAC
Definition: constants.h:59
LIGHTCONE_MASSMAP_ALLOC_FAC
#define LIGHTCONE_MASSMAP_ALLOC_FAC
Definition: constants.h:66
integertime
int integertime
Definition: constants.h:331
MAX_DOUBLE_NUMBER
#define MAX_DOUBLE_NUMBER
Definition: constants.h:81
TIMEBASE
#define TIMEBASE
Definition: constants.h:333
M_PI
#define M_PI
Definition: constants.h:56
mycxxsort
double mycxxsort(T *begin, T *end, Tcomp comp)
Definition: cxxsort.h:39
MyIntPosType
uint32_t MyIntPosType
Definition: dtypes.h:35
MySignedIntPosType
int32_t MySignedIntPosType
Definition: dtypes.h:36
BITS_FOR_POSITIONS
#define BITS_FOR_POSITIONS
Definition: dtypes.h:37
Logs
logs Logs
Definition: main.cc:43
Terminate
#define Terminate(...)
Definition: macros.h:15
warn
#define warn(...)
Definition: macros.h:30
Driftfac
driftfac Driftfac
Definition: main.cc:41
TAG_DENS_A
#define TAG_DENS_A
Definition: mpi_utils.h:50
myMPI_Alltoall
int myMPI_Alltoall(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, int recvcount, MPI_Datatype recvtype, MPI_Comm comm)
Definition: myalltoall.cc:301
myMPI_Sendrecv
int myMPI_Sendrecv(void *sendbuf, size_t sendcount, MPI_Datatype sendtype, int dest, int sendtag, void *recvbuf, size_t recvcount, MPI_Datatype recvtype, int source, int recvtag, MPI_Comm comm, MPI_Status *status)
Definition: sizelimited_sendrecv.cc:25
Mem
memory Mem
Definition: main.cc:44
MMM
#define MMM(x, y)
Definition: mymalloc.h:23
half_float::detail::acos
expr acos(half arg)
Definition: half.hpp:2844
half_float::detail::exp
expr exp(half arg)
Definition: half.hpp:2724
half_float::detail::fabs
half fabs(half arg)
Definition: half.hpp:2627
half_float::detail::atan2
expr atan2(half x, half y)
Definition: half.hpp:2859
half_float::detail::ceil
half ceil(half arg)
Definition: half.hpp:2950
half_float::detail::log
expr log(half arg)
Definition: half.hpp:2745
half_float::detail::pow
expr pow(half base, half exp)
Definition: half.hpp:2803
half_float::detail::sqrt
expr sqrt(half arg)
Definition: half.hpp:2777
global_data_all_processes::BoxSize
double BoxSize
Definition: allvars.h:144
global_data_all_processes::SnapFormat
int SnapFormat
Definition: allvars.h:46
global_data_all_processes::Ti_Current
integertime Ti_Current
Definition: allvars.h:188
global_data_all_processes::ComovingIntegrationOn
int ComovingIntegrationOn
Definition: allvars.h:148
global_data_all_processes::TimeBegin
double TimeBegin
Definition: allvars.h:180
global_data_all_processes::Time
double Time
Definition: allvars.h:179
global_data_all_processes::Timebase_interval
double Timebase_interval
Definition: allvars.h:187
particle_data
Definition: particle_data.h:35
particle_data::getMass
MyDouble getMass(void)
Definition: particle_data.h:176
particle_data::getSofteningClass
unsigned char getSofteningClass(void)
Definition: particle_data.h:196
particle_data::Vel
MyFloat Vel[3]
Definition: particle_data.h:54
particle_data::ID
MyIDStorage ID
Definition: particle_data.h:70
particle_data::getType
unsigned char getType(void)
Definition: particle_data.h:140
particle_data::GravAccel
vector< MyFloat > GravAccel
Definition: particle_data.h:55
particle_data::IntPos
MyIntPosType IntPos[3]
Definition: particle_data.h:53
subdivide_evenly_get_bin
void subdivide_evenly_get_bin(int N, int pieces, int index, int *bin)
Definition: system.cc:70
subdivide_evenly
void subdivide_evenly(long long N, int pieces, int index_bin, long long *first, int *count)
Definition: system.cc:51