Dehnen's code bugs

[greenyellowcat]

Dear Peter Teuben, Sorry to bother you, but I encountered some errors after upgrading to the new version of nemo. Specifically, mkgalaxy utilite stopped working. Error log here:

file tst.err

mkgalaxy error output

disc parameters: Md = 1 (disc mass) Nd = 100 (number of disc bodies) Rd = 1 (disc scale radius) Rdmax = 4.5 (disc truncation radius) Zd = 0.05 (disc scale height) Rsig = 2 (if !=0 : scale radius for sigma_R) Q = 1.2 (Toomre's Q: constant if Rsig=0, otherwise Q(Rsig)) Nbpo = 50 (number of disc bodies sampled per orbit) ni = 4 (number of iterations in disc sampling) epsd = 0.01 (gravitational softening length for disc bodies) halo parameters: Mh = 8 (halo mass) Nh = 400000 (number of halo bodies) innerh = 1 (halo inner logarithmic density slope) outerh = 3 (halo outer logarithmic density slope) etah = 1 (halo transition exponent) Rcoreh = 0 (halo core radius) Rh = 6 (halo scale length) Rth = 100 (halo truncation radius) betah = 0 (halo anisotropy parameter) r_ah = 0 (halo anisotropy radius; 0 maps to infinity) epsh = 0.02 (gravitational softening length for halo bodies) bulge parameters: Mb = 0.2 (bulge mass) Nb = 20 (number of bulge bodies) innerb = 1 (bulge inner density exponent) outerb = 4 (bulge outer density exponent) etab = 1 (bulge transition exponent) Rcoreb = 0 (bulge core radius) Rb = 0.2 (bulge scale radius) Rtb = 0 (bulge truncation radius) betab = 0 (bulge anisotropy parameter) r_ab = 0 (bulge anisotropy radius; 0 maps to infinity) epsb = 0.01 (gravitational softening length for bulge bodies) parameters controlling code: kmax = 3 (maximum timestep = 2^-kmax) kmin = 7 (minimum timestep = 2^-kmin) fac = 0.01 (time step control: tau < fac/|acc|) fph = 0.04 (time step control: tau < fph/|phi|) tgrow = 40 (disc growth time) seed = 1 (seed for RNGs) nmax = 12 (maximum n in potential expansion) lmax = 8 (maximum l in potential expansion) debug = 2 (debug level used to run all falcON programs)

======================================== mkgalaxy: issuing command: rm -f tst.h tst.b tst.s tst.prm >& /dev/null

======================================== mkgalaxy: issuing commands: echo accname=DiscPot > tst.prm echo accpars=0,37681.3,1,-0.025,0,0 >> tst.prm

======================================== mkgalaxy: issuing command: mkhalo out=tst.h! nbody=200000 M=8 inner=1 outer=3 eta=1 r_s=6 r_t=100 r_c=0 b=0 r_a=0
seed=1 eps=0.02 giveF=f accname=Halo+Monopole
accpars=0,0.2,0.2,1,4,1,0,0;1,10
accfile=;tst.prm debug=2 >>& tst.err

nemo Debug Info: Keyword file will be mkhalo.def

nemo Debug Info: savehist: progname=mkhalo help_level=0

nemo Debug Info: Found -1 output keywords

nemo Debug Info: set_xrandom(NUMREC portable) seed=1

nemo Debug Info: hasvalue: checking indexing on eps

nemo Debug Info: hasvalue: checking indexing on accname

nemo Debug Info: hasvalue: checking indexing on accpars

nemo Debug Info: hasvalue: checking indexing on accfile

nemo Debug Info: get_acceleration("Halo+Monopole","0,0.2,0.2,1,4,1,0,0;1,10",";tst.prm")

nemo Debug Info: single_acceleration("Halo", "0,0.2,0.2,1,4,1,0,0", "(null)")

nemo Debug Info: MySymbols: NULL code in loadobjDL

nemo Debug Info: loadobj: /home/ssd/opt/nemo_newest/nemo/usr/dehnen/falcON/acc//Halo.so

nemo Debug Info: findfn: looking up iniacceleration

falcON Debug Info [inc/public/halo.h:445]: ModifiedDoublePowerLawHalo: rh0=3.978874

falcON Debug Info [src/public/lib/halo.cc:434]: HaloPotential: minimum & maximum tabulated radii = 0.000390625, 102.4

falcON Debug Info [src/public/acc/Halo.cc:87]: external potential "Halo" recognizes 8 parameters:

omega pattern speed (ignored) [0] r_s scale radius [1] m_t total mass [1] g_i inner power-law slope [7/9] g_o outer power-law slope [31/9] eta transition steepness [4/9] r_t truncation radius (0->oo) [0] r_c core radius [0] The halo density proportional to

  Model(x) * sech(r/r_t)

with x=sqrt(r^2+r_c^2)/r_s, and -gi eta [gi-go]/eta Model(x) = x (x + 1).

If file is given, it's interpreted as follows: Plummer: gi=0, go=5, eta=2 Jaffe: gi=2, go=4, eta=1 Hernquist: gi=1, go=4, eta=1 Dehnen: go=4, eta=1 NFW: gi=1, go=4, eta=1 Moore: gi=3/2, go=3, eta=3/2 DM: gi=7/9, go=31/9, eta=4/9 and differing values are ignored.

falcON Debug Info [src/public/acc/Halo.cc:117]: external potential "Halo" initialized with:

r_s = 0.2 m_t = 0.2 g_i = 1 g_o = 4 eta = 1 r_t = 0 r_c = 0

nemo Debug Info: single_acceleration("Monopole", "1,10", "tst.prm")

nemo Debug Info: loadobj: /home/ssd/opt/nemo_newest/nemo/usr/dehnen/falcON/acc//Monopole.so

nemo Debug Info: findfn: looking up iniacceleration

nemo Debug info: Monopole:

provides the acceleration field due to the potential

Phi_0(x) + A(t) * [Phi(x) - Phi_0(x)],

where Phi(x) is a given conservative potential, Phi_0(x) its monopole and A(t) an adiabatic growth factor with A=0 at t<t0 and A=1 at t>t0+tau. A data file is required and must be of the form accname=NAME [accpars=PARS] [accfile=FILE] characters after (and including) '#' are ignored. Parameters: 5 with the meanings: par[0] = t0: start time for growth [0] par[1] = tau: time scale for growth [1] par[2] = innermost radius in table for monopole [0.001] par[3] = outermost radius in table for monopole [1000] par[4] = number of points in table for monopole [1001]

nemo Debug Info: get_acceleration("DiscPot","0,37681.3,1,-0.025,0,0","(null)")

nemo Debug Info: single_acceleration("DiscPot", "0,37681.3,1,-0.025,0,0", "(null)")

nemo Debug Info: loadobj: /home/ssd/opt/nemo_newest/nemo/usr/dehnen/falcON/acc//DiscPot.so

nemo Debug Info: findfn: looking up iniacceleration

nemo debug info:

external potential "DiscPot" recognizes 5 parameters: omega pattern speed (ignored) Sig_0 central surface density [1] R_d disc scale length [1] z_d disc scale height [0.1] R_0 radius of central hole [0] eps cosine modulation term [0] The disc density is given by rho(R,z)=Sigma(R)*h(z) with

Sigma(R) = Sig_0 exp(-R_0/R-R/R_d+eps*cos[R/R_d])

      exp(-z/z_d) / (2 z_d)          if z_d > 0

h(z) = delta(z) if z_d = 0 sech^2(z/2|z_d|) / (4|z_d|) if z_d < 0

nemo debug info:

external potential "DiscPot" initialized with: Sig_0 = 37681.3 R_d = 1 z_d = -0.025 R_0 = 0 eps = 0

nemo Debug Info: hasvalue: checking indexing on inner

nemo Debug Info: hasvalue: checking indexing on outer

nemo Debug Info: hasvalue: checking indexing on eta

nemo Debug Info: hasvalue: checking indexing on r_2

falcON Debug Info [inc/public/halo.h:445]: ModifiedDoublePowerLawHalo: rh0=0.001375

nemo Debug Info: hasvalue: checking indexing on max_r

falcON Debug Info [src/public/lib/halo.cc:434]: HaloPotential: minimum & maximum tabulated radii = 0.0234375, 12288

falcON Error: [src/public/lib/halo.cc:500]: HaloPotential: external Rh(98.3572)=-9.323822386764764e-13

terminate called after throwing an instance of 'WDutils::exception' what(): [src/public/lib/halo.cc:500]: HaloPotential: external Rh(98.3572)=-9.323822386764764e-13

PS1: I varied the parameters of the halo. The results are all the same. The bug appears. PS2: I am a regular user of your software and am very grateful to you for your cumbersome work. It would be great if these issues could be resolved.

[teuben]

I probably should know what OS you are using. And secondly, how did you upgrade. You can either install a new one, or use "git pull" and then "mknemo gyrfalcON".

If you start from scratch, then the doc/install_nemo script is the easiest.

Coming back to "mkgalaxy", I also have a core dump on my system, mkhalo is the one that core dumps on my system (Ubuntu 18.04). The executable is fine, the help= parameter works, so it's not a shared library issue.

[greenyellowcat]

OS version: openSUSE Leap 15.0 Upgrade: I started from scratch and used the instructions written in your README file here on GitHub. Tried to recompile gyrfalcON with ``mknemo gyrfalcON'' command. Checked mkgalaxy after, still got the same error.

[teuben]

I started a fresh install, and indeed, mkgalaxy dies in mkhalo (core dump) for me. The reason why wasn't immediately obvious to me. This command reproduces the problem (which I've added to the testsuite for 'dehen')

mkhalo out=mkhalo1.dat nbody=6000 M=24 inner=7/9 outer=31/9 eta=4/9 r_s=6 r_t=60 r_c=0 b=0 r_a=0 seed=1 eps=0.02 giveF=f accname=Halo accpars=0,0.2,0.2,1,4,1,0,0 debug=2

[greenyellowcat]

I'm sorry for the late reply. The test is certainly good. But is there any hope of solving the problem itself?

[teuben]

I've taken a very unconventional approach to solving this: the code complained numbers were very close to 0, and overshooting. So I patched them at 0. The resulting model then "looks" ok, but this clearly needs to be better investigated.

For now, there will be a big fat warning, something along the lines

falcON Warning [src/public/lib/halo.cc:512]: HaloModel: PJT patched at 1386/1746

it's quite possible the model that is produced is bogus. So i'm not closing this issue yet and let you play with this.