search

cloudy-astrophysics / bug-tracker-migration-test

Trial run for importing the nublado.org Trac tickets as GitHub issues
0 stars 0 forks source link

line shielding function is poorly behaved for large damping coefficient (trac #272) #274

Closed cloudy-bot closed 5 years ago

cloudy-bot commented 10 years ago

reported by: @CloudyLex

We now use the Federman+79 line self shielding function for all lines in the code. It was derived for small damping coefficients and the damping wings are now treated properly for large damping coefficients.

http://cdsads.u-strasbg.fr/abs/1979ApJ...227..466F

used at rt_continuum_shield_fcn.cpp:63 see google groups thread under "convergence branch" dates 2013 nov 24 starting with Ryan email.

The relevant part of Federman+79 is the bit around A3 and A9 - they do not use standard stellar atmosphere (ie, Chandra) notation but gamma/beta is the ratio of damping to doppler widths, so r is our damping coefficient.

A9 handles the damping wings - looks like the function becomes large when r is large. From Fig 5 of Ferland+13 (Fe UTA) we can get a as large as 1000. we are likely saved by the min at line 111 - otherwise the shield would become quite large at small tau, large a.

Migrated from https://www.nublado.org/ticket/272

{
    "status": "closed",
    "changetime": "2013-12-19T00:21:09Z",
    "_ts": "1387412469000000",
    "description": "We now use the Federman+79 line self shielding function for all lines in the code.  It was derived for small damping coefficients and the damping wings are now treated properly for large damping coefficients.\n\nhttp://cdsads.u-strasbg.fr/abs/1979ApJ...227..466F\n\nused at rt_continuum_shield_fcn.cpp:63\nsee google groups thread under \"convergence branch\" dates 2013 nov 24 starting with Ryan email.\n\nThe relevant part of Federman+79 is the bit around A3 and A9 - they do not use standard stellar atmosphere (ie, Chandra) notation but gamma/beta is the ratio of damping to doppler widths, so r is our damping coefficient.  \n\nA9 handles the damping wings - looks like the function becomes large when r is large.  From Fig 5 of Ferland+13 (Fe UTA) we can get a as large as 1000.  we are likely saved by the min at line 111 - otherwise the shield would become quite large at small tau, large a.",
    "reporter": "gary",
    "cc": "",
    "resolution": "fixed",
    "time": "2013-11-29T18:07:17Z",
    "component": "radiative transfer",
    "summary": "line shielding function is poorly behaved for large damping coefficient",
    "priority": "critical",
    "keywords": "",
    "version": "trunk",
    "milestone": "",
    "owner": "nobody",
    "type": "defect - wrong answer"
}
cloudy-bot commented 10 years ago

@CloudyLex commented:

ran auto at r8626 with attached patch, which caps a at 1e-2 - federman says approximation optimized for 1 = 1e-3

          double DampUsed = MIN2( 1e-2, t.Emis().damp() );
          double t1 = 3.02*pow(DampUsed*1e3,-0.064 );
          double u1 = sqrt(MAX2(tau,0.)*DampUsed )/SDIV(t1);
          wings = DampUsed/SDIV(t1)/sqrt( 0.78540 + POW2(u1) );

following botches result - some big

wolkje auto:cat serious.asr
agn_warm_absorber.out:  ChkMonitor botch>>IRON        9    -1.4061 =    -1.6370  -0.702   0.100  <<BIG BOTCH!!
agn_warm_absorber.out:  ChkMonitor botch>>IRON       10    -0.9514 =    -1.1490  -0.576   0.100  <<BIG BOTCH!!
agn_warm_absorber.out:  ChkMonitor botch>>IRON       11    -0.7946 =    -0.9535  -0.442   0.100  <<BIG BOTCH!!
agn_warm_absorber.out:  ChkMonitor botch>>IRON       12    -0.8528 =    -0.9557  -0.267   0.100
agn_warm_absorber.out:  ChkMonitor botch>>IRON       14    -1.2452 =    -1.1820   0.135   0.100
agn_warm_absorber.out:  ChkMonitor botch>>IRON       15    -1.4551 =    -1.3810   0.157   0.100
agn_warm_absorber.out:  ChkMonitor botch>>IRON       16    -1.5361 =    -1.4340   0.210   0.100
agn_warm_absorber.out:  ChkMonitor botch>>IRON       17    -0.7391 =    -0.5951   0.282   0.100
agn_warm_absorber.out:  ChkMonitor botch>>IRON       18    -0.9801 =    -0.8353   0.284   0.100
agn_warm_absorber.out:  ChkMonitor botch>>IRON       19    -1.5491 =    -1.4040   0.284   0.100
agn_warm_absorber.out:  ChkMonitor botch>>IRON       20    -2.6694 =    -2.5230   0.286   0.100
agn_warm_absorber.out:  ChkMonitor botch>>IRON       21    -4.1280 =    -3.9780   0.292   0.100
blr_level2.out:  ChkMonitor botch>>P  5 1117.98A     6.4717 =     6.4990   0.061   0.050
blr_level2.out:  ChkMonitor botch>>P  5 1128.01A     6.2315 =     6.2600   0.064   0.050
blr_level2.out:  ChkMonitor botch>>PHOS        5    -1.9042 =    -1.8690   0.078   0.050
blr_level2.out:  ChkMonitor botch>>PHOS        6    -2.9959 =    -2.9480   0.104   0.100
blr_n09_p20_Z20.out:  ChkMonitor botch>>totl 1035.00A     7.9767 =     8.0000   0.052   0.050
blr_n13_p18_Z20.out:  ChkMonitor botch>>totl 2326.00A     3.2033 =     3.2270   0.053   0.050
blr_nf84_45deg.out:  ChkMonitor botch>>HELI        3    -2.5043 =    -2.4820   0.050   0.050
blr_rnfa.out:  ChkMonitor botch>>c  3 977.000A     2.4194 =     2.5540   0.053   0.050
blr_rnfa.out:  ChkMonitor botch>>blnd 2798.00A    27.1653 =    25.5300  -0.064   0.050
blr_rnfa.out:  ChkMonitor botch>>Blnd 1888.00A     7.3824 =     7.7900   0.052   0.050
nlr_lex00.out:  ChkMonitor botch>>blnd 2424.00A     0.5035 =     0.5310   0.052   0.050
nlr_liner_grains.out:  ChkMonitor botch>>HELI        2    -0.7413 =    -0.7706  -0.070   0.050
nlr_paris.out:  ChkMonitor botch>>totl 1240.00A     0.2727 =     0.2979   0.085   0.050
nlr_paris.out:  ChkMonitor botch>>totl 1402.00A     0.5242 =     0.5596   0.063   0.050
nlr_paris.out:  ChkMonitor botch>>TOTL 1218.00A     0.2833 =     0.3132   0.095   0.050
nlr_paris.out:  ChkMonitor botch>>TOTL 1035.00A     0.1245 =     0.1452   0.143   0.100
nlr_paris.out:  ChkMonitor botch>>Blnd 1397.00A     0.2380 =     0.2527   0.058   0.050
pdr_leiden_f2.out:  ChkMonitor botch>>Si 2 34.8046m    -4.4268 =    -4.2150   0.386   0.050  <<BIG BOTCH!!
pdr_leiden_f4.out:  ChkMonitor botch>>Si 2 34.8046m    -4.3035 =    -4.2360   0.144   0.050
pdr_leiden_v2.out:  ChkMonitor botch>>Si 2 34.8046m    -3.5455 =    -3.5170   0.063   0.050
pn_fluc.out:  ChkMonitor botch>>blnd 2798.00A     0.9817 =     0.8152  -0.204   0.050  <<BIG BOTCH!!
pn_ots.out:  ChkMonitor botch>>blnd 2798.00A     1.2078 =     0.9973  -0.211   0.050  <<BIG BOTCH!!
pn_paris.out:  ChkMonitor botch>>blnd 2798.00A     1.2327 =     1.0590  -0.164   0.050  <<BIG BOTCH!!
pn_paris_cpre.out:  ChkMonitor botch>>blnd 2798.00A     1.3556 =     1.1280  -0.202   0.050  <<BIG BOTCH!!
cloudy-bot commented 10 years ago

@rjrw commented:

See

http://adsabs.harvard.edu/abs/1974JQSRT..14..319R

for a better approximation.

cloudy-bot commented 10 years ago

r8695 implements fixes to the Federman shielding function to bring it within -15% / +32% of the value obtained from exact integration across a wide range of tau, a (from -80% / +44% previously). As surmised in the original report, the min() does indeed maintain reasonable behaviour at small tau/large a.

With this change, the accuracy isn't substantially worse than the Rodgers & Williams fit used elsewhere in the literature (and now implemented as an option).