Split geometry - Githubissues

bprather commented 5 years ago

The big one.

Add basic tests. These are just flat binary comparison vs previous, canonical version, for 3 problems, 2 of which are images of a particular GRMHD dump (hosted on BH).
Split out most header material from decs.h to individual headers. This allows re-use of C files without disentangling.
Move most useful coordinate/geometry functions out of model.c and tetrads.c in favor of coordinates.c and geometry.c. This could allow sharing with HARM/pyHARM/ebhlight/etc. in future to avoid repeatedly re-implementing coordinates and grids.

TODO:

Not all function declarations are split into headers. More generally, this is the first workable version which passes tests, not a final vision for ipole.
Coordinate, geometry, and model parameters are largely still shared via globals, albeit with marginally improved namespacing. Fixing this in a convenient way is an open problem.
Parameter handling is likewise not devolved to model/coordinate files.
Thin disk problem is still missing, and more generally handling of analytic emissivity distributions is not decided.
Documentation is lacking

gnwong commented 5 years ago

I'm going to reject this request (at commit 6046c95) pending further review. I've run this version of ipole and the master branch version on a dump (bh:/bd4/dumps/tests/dump.h5) with the following set of parameters and got significantly different results. It's unclear where the difference is coming from, but until this is understood, it seems dangerous to merge.

Parameters are

--dump=dump.h5 
--Mbh=6.2e9 
--M_unit=1.e25 
--freqcgs=230.e9 
--thetacam=73 
--outfile=image.h5

Example output from plot_pol.py

Flux [Jy]:    0.73831 0.768636
I,Q,U,V [Jy]: 0.73831 0.0140543 0.00415017 -0.00325547
LP [%]:       1.98483
CP [%]:       -0.440935
EVPA [deg]:   8.22585

Flux [Jy]:    0.874436 0.909812
I,Q,U,V [Jy]: 0.874436 -0.0085912 0.00068545 -0.00465024
LP [%]:       0.985606
CP [%]:       -0.531799
EVPA [deg]:   87.7192

The difference may just be due to my stupidity in incorrectly identifying parameters (if this has changed significantly, documentation is probably important).

Text on standard out is

$ ./ipole_master -par par.par
ipole. githash: 38da6a4
notes:

nthreads = 72
reading data header...
filename: dump.h5
init grid
custom electron model loaded from dump file...
using mixed tp_over_te with trat_small = 1 and trat_large = 40
custom refinement at poles loaded...
start: 0.196852 0 0
stop: 6.90776 1 6.28319
success
MBH: 6.2e+09 [Msun]
L,T,M units: 9.15766e+14 [cm] 30546.6 [s] 1e+25 [g]
rho,u,B units: 1.30211e-20 [g cm^-3] 11.7028 [g cm^-1 s^-2] 12.1269 [G]
reading data...
LOADING DATA
dMact: -22.5114 [code]
Ladv: 5.52615 [code]
Mdot: 7.36952e+21 [g/s]
Mdot: 0.000116845 [MSUN/YR]
Mdot: 8.49589e-07 [Mdotedd]
Mdotedd: 8.67422e+27 [g/s]
Mdotedd: 137.532 [MSUN/YR]
success
Xcam[] = 0.000000e+00 6.907755e+00 2.995930e-01 0.000000e+00
cam_th_cal=0.996947 [deg] th_beg=0.0174 th_len=0 a=0.9375 R0=0 hslope=0.3
L_unit = 9.157655e+14 DX = 4.417204e+01 NX = 160 Dsource = 5.214796e+25 JY = 1.000000e-23
intensity [cgs] to flux per pixel [Jy] conversion: 2.35044
Dsource: 5.2148e+25 [cm]
Dsource: 16900 [kpc]
FOVx, FOVy: 44.172 44.172 [GM/c^2]
FOVx, FOVy: 7.75701e-10 7.75701e-10 [rad]
FOVx, FOVy: 160 160 [muas]
0 2 3 3 4 5 5 9 9 10 10 11 12 12 13 13 14 14 14 14 15 16 16 16 17 18 19 20 22 22 23 24 26 26 26 26 27 28 29 30 31 32 33 34 35 35 36 36 37 38 38 39 39 40 41 42 42 43 44 45 46 47 48 48 49 50 50 52 53 54 55 56 57 57 59 60 61 62 63 64 65 66 66 67 68 69 70 70 71 71 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 93 94 95 96 97 98 99 100 101 102 102 102 103 104 104 105 106 107 107 108 110 111 112 112 113 114 115 116 117 118 119 120 120 121 122 122 123 124 124 125 125 126 127 128 129 129 129 130 131 132 133 134 136 137 138 139 140 140 141 142 143 144 144 145 146 147 149 150 151 152 153 153 154 155 155 156 157 158 159 159
scale = 2.350439e+00
imax=67 jmax=80 Imax=0.0022888 Iavg=1.22701e-05
freq: 2.3e+11 Ftot: 0.73831 (0.768636 unpol) scale=2.35044
nuLnu = 5.80298e+40
Total wallclock time: 27.4122 s

$ ./ipole_split-geometry --dump=dump.h5 --Mbh=6.2e9 --M_unit=1.e25 --freqcgs=230.e9 --thetacam=73 --outfile=out_new.h5
ipole. githash: 6046c95
notes:

nthreads = 72
reading data header...
filename: dump.h5
custom electron model loaded from dump file...
using mixed tp_over_te with trat_small = 1 and trat_large = 30
custom refinement at poles loaded...
start: 0.196852 0 0
stop: 6.90776 1 6.28319
success
MBH: 6.2e+09 [Msun]
L,T,M units: 9.15766e+14 [cm] 30546.6 [s] 1e+25 [g]
rho,u,B units: 1.30211e-20 [g cm^-3] 11.7028 [g cm^-1 s^-2] 12.1269 [G]
reading data...
LOADING DATA
dMact: -22.5114 [code]
Ladv: 5.52615 [code]
Mdot: 7.36952e+21 [g/s]
Mdot: 0.000116845 [MSUN/YR]
Mdot: 8.49589e-07 [Mdotedd]
Mdotedd: 8.67422e+27 [g/s]
Mdotedd: 137.532 [MSUN/YR]
success
Xcam[] = 0.000000e+00 6.907755e+00 2.995932e-01 0.000000e+00
a=0.9375 R0=0 hslope=0.3
No FOV was specified. Using default 160muas!
L_unit = 9.157655e+14 DX = 4.417204e+01 NX = 160 Dsource = 5.214796e+25 JY = 1.000000e-23
intensity [cgs] to flux per pixel [Jy] conversion: 2.35044
Dsource: 5.2148e+25 [cm]
Dsource: 16900 [kpc]
FOVx, FOVy: 44.172 44.172 [GM/c^2]
FOVx, FOVy: 7.75701e-10 7.75701e-10 [rad]
FOVx, FOVy: 160 160 [muas]
0 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 158 159
scale = 2.350439e+00
imax=66 jmax=80 Imax=0.00245458 Iavg=1.45325e-05
freq: 2.3e+11 Ftot: 0.874436 (0.909812 unpol) scale=2.35044
nuLnu = 6.87291e+40
Total wallclock time: 32.5704 s

bprather commented 5 years ago

Default trat_large is now 30 apparently. I thought it was so before, but evidently not. I'll change the default. Note there are two breaking changes we need to make:

Launching geodesics from half-pixels (done)
Re-centering camera to be in ZAMO frame like grtrans (not in this PR)

I've tried to accommodate these in the tests, which spell out problems to avoid breaking on changes in defaults like this

bprather commented 5 years ago

Output from plot_pol.py with the following commandline (both versions, xoff/yoff ignored in master branch). Let me know how to document this -- backward compat section in README?

./ipole --dump=dump.h5 --MBH=6.2e9 --M_unit=1.e25 --freqcgs=230.e9 --thetacam=73 --outfile=image.h5 --xoff=-0.5 --yoff=-0.5

plotting ipole-master/image.h5 Flux [Jy]: 0.73831 0.768636 I,Q,U,V [Jy]: 0.73831 0.0140543 0.0041502 -0.00325547 LP [%]: 1.98484 CP [%]: -0.440935 EVPA [deg]: 8.22588

plotting ipole-split-geom/image.h5 Flux [Jy]: 0.738309 0.768636 I,Q,U,V [Jy]: 0.738309 0.0140568 0.00413986 -0.00325535 LP [%]: 1.98476 CP [%]: -0.44092 EVPA [deg]: 8.20513

AFD-Illinois / ipole

Split geometry #6