-
Notifications
You must be signed in to change notification settings - Fork 7
/
poissondrivers.c
157 lines (130 loc) · 4.38 KB
/
poissondrivers.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <string.h>
#include <fftw3.h>
#include <mpi.h>
#include <hdf5.h>
#include <gsl/gsl_math.h>
#include <string.h>
#include "raytrace.h"
void fullsky_partdist_poissondriver(void)
{
double time;
long doPoissonSolve;
long totNlensPlaneParts;
//read parts
if(!rayTraceData.UseHEALPixLensPlaneMaps)
{
time = -MPI_Wtime();
logProfileTag(PROFILETAG_PARTIO);
read_lcparts_at_planenum_fullsky_partdist(rayTraceData.CurrentPlaneNum);
get_smoothing_lengths();
logProfileTag(PROFILETAG_PARTIO);
time += MPI_Wtime();
if(ThisTask == 0)
fprintf(stderr,"read %ld parts in %g seconds.\n",NlensPlaneParts,time);
}
/*
if we read particles, then solve poisson eqn. Otherwise, even if
backdens is zero or non-zero, we do not need to solve the poisson eqn
*/
MPI_Allreduce(&NlensPlaneParts,&totNlensPlaneParts,1,MPI_LONG,MPI_SUM,MPI_COMM_WORLD);
if(totNlensPlaneParts > 0)
doPoissonSolve = 1;
else
doPoissonSolve = 0;
#ifdef DEBUG
#if DEBUG_LEVEL > 1
if(NlensPlaneParts > 0)
fprintf(stderr,"%d: pos,mass = %f|%f|%f|%e, total # of parts = %ld, # of parts on this task = %ld\n",
ThisTask,lensPlaneParts[0].pos[0],lensPlaneParts[0].pos[1],lensPlaneParts[0].pos[2],lensPlaneParts[0].mass,
totNlensPlaneParts,NlensPlaneParts);
else
fprintf(stderr,"%d: # of parts on this task = %ld\n",ThisTask,NlensPlaneParts);
#endif
#endif
//run poisson solver
if(doPoissonSolve)
{
do_healpix_sht_poisson_solve(rayTraceData.densfact,rayTraceData.backdens);
//do MG step if needed
#ifndef SHTONLY
logProfileTag(PROFILETAG_PARTIO);
read_lcparts_at_planenum(rayTraceData.CurrentPlaneNum);
get_smoothing_lengths();
logProfileTag(PROFILETAG_PARTIO);
#ifdef DEBUG_IO_DD
write_bundlecells2ascii("preMGPS");
#endif
mgpoissonsolve(rayTraceData.densfact,rayTraceData.backdens);
#endif
}
//free parts
destroy_parts();
}
void cutsky_partdist_poissondriver(void)
{
double time;
long doPoissonSolve;
//read parts
#ifdef SHTONLY
if(!rayTraceData.UseHEALPixLensPlaneMaps)
{
#endif
time = -MPI_Wtime();
logProfileTag(PROFILETAG_PARTIO);
read_lcparts_at_planenum(rayTraceData.CurrentPlaneNum);
get_smoothing_lengths();
logProfileTag(PROFILETAG_PARTIO);
time += MPI_Wtime();
if(ThisTask == 0)
fprintf(stderr,"read %ld parts in %g seconds.\n",NlensPlaneParts,time);
#ifdef SHTONLY
}
#endif
/* if there are particles read into memory, then we need to solve the poisson equation
otherwise, there are two cases
1) if backdens == 0, then rho == 0 everywhere and phi = 0 as well
2) if backdens != 0, then rho is negative in the region being ray traced and zero elsewhere, so still need to solve the poisson equation
*/
#ifdef NOBACKDENS
long totNlensPlaneParts;
MPI_Allreduce(&NlensPlaneParts,&totNlensPlaneParts,1,MPI_LONG,MPI_SUM,MPI_COMM_WORLD);
if(totNlensPlaneParts > 0)
doPoissonSolve = 1;
else
doPoissonSolve = 0;
#else
doPoissonSolve = 1;
#endif
#ifdef DEBUG
#if DEBUG_LEVEL > 1
#ifndef NOBACKDENS
long totNlensPlaneParts;
MPI_Allreduce(&NlensPlaneParts,&totNlensPlaneParts,1,MPI_LONG,MPI_SUM,MPI_COMM_WORLD);
#endif
if(NlensPlaneParts > 0)
fprintf(stderr,"%d: pos,mass = %f|%f|%f|%e, total # of parts = %ld, # of parts on this task = %ld\n",
ThisTask,lensPlaneParts[0].pos[0],lensPlaneParts[0].pos[1],lensPlaneParts[0].pos[2],lensPlaneParts[0].mass,
totNlensPlaneParts,NlensPlaneParts);
else
fprintf(stderr,"%d: # of parts on this task = %ld\n",ThisTask,NlensPlaneParts);
#endif
#endif
//run poisson solver
if(doPoissonSolve)
{
do_healpix_sht_poisson_solve(rayTraceData.densfact,rayTraceData.backdens);
//do MG step if needed
#ifndef SHTONLY
#ifdef DEBUG_IO_DD
write_bundlecells2ascii("preMGPS");
#endif
mgpoissonsolve(rayTraceData.densfact,rayTraceData.backdens);
#endif
}
//free parts since we do not need them anymore
destroy_parts();
}