drvdputt/radagast

Gas module for dust radiative transfer codes

RADAGAST

RADiation and Dust Affecting the GAs STate

Description

RADAGAST is a local (0D) gas model for use in 3D dust radiative transfer simulations. Given a certain radiation field and dust model, the equilibrium state for various quantities of the gas are calculated self-consistently. Currently, the module treats the main species consisting of H (H+, H and H2), and their interactions with the radiation field and the dust grains.

Quantities calculated self-consistently

• Gas equilibrium temperature
• H+, H, and H2 abundance
• Level populations of H and H2
• Grain temperatures and charge distributions

A non-comprehensive list of processes included in this calculation

• Photoionization of H
• Radiative recombination of H+
• Photodissociation of H2
• Formation of H2 on dust grain surfaces
• Heating by the photoelectric effect on dust grains
• Cooling by collisions between gas and dust grains
• Net heating/cooling by collisional transitions of H2 (and H at high temperatures)

From the resulting state, the emissive properties and cross section of the gas can be calculated, so that a radiative transfer code can include these in its model.

• Cross sections for H photoionization, H2 photoionization, and H2 photodissociation
• Continuum emission of H (free-bound, free-free, and two-photon)
• Emission coefficients, opacities, and shapes for the lines of H and H2.

While RADAGAST was mainly developed for use in the 3D Monte Carlo Dust Radiative Transfer code SKIRT, it is designed as a stand-alone module. It is possible to use this module in any C++ code where the necessary details about the radiation field and the dust can be provided, without having to install anything related to SKIRT.

This code was developed during my PhD research, and contains many equations and data from other authors to implement the processes listed above. Consult my PhD thesis (to be published) for details, equations, and citations. Alternatively, one can build the documentation pages or read the comments for each function in the source code.

Third party libraries and data

RADAGAST includes a release of the doctest C++ testing framework, to power a small set of unit tests.

It also uses Eigen3 for storing vectors and matrices, and performing linear algebra, and GSL for solving several non-linear equations. Currently, these libraries are dependencies, but the git history might contain a copy.

Getting started

First read how to build the library.

Then try running the main binary produced at RADAGAST/../cmake_release/src/mains/main. Experiment by giving a few numerical command line arguments to this program. The command main 1e5 1e4 1e3 will run an example model with a gas density of 1e5 cm-3, under the effect of a blackbody source located at a (fixed) distance of 1 pc with a color temperature of 1e4 K and a bolometric luminosity of 1e3 solar luminosities.

When this binary works without problems, copy the library built at RADAGAST/../cmake_release/src/core/libridge.a to your preferred location (or set CMAKE_INSTALL_PREFIX and run make install). Then read the instructions on how to use the API, write a simple C++ script making the calls, and try building and linking your own C++ script to the library.

Integration in SKIRT

To store, use, and update gas properties, we make use of the newly developed dynamic state mechanisms for SKIRT. See pull request. The branch of SKIRT that integrates RADAGAST can be found on my fork, and is called gas-dev. This branch should be compatible with the latest commit of RADAGAST.

The versions of RADAGAST/SKIRT used for my PhD thesis are older. See the thesis release tag for SKIRT, and the thesis release tag for RADAGAST.

radagast's Issues

Update documentation

Needs note about shielding factor, since it's an extra argument at the main interface.

Crash when all wavelengths of grid are > 91.2 nm

The problem lies with the integration of the ionizing radiation. The minimum index for integrating is larger than the maximum index, which leads to a bad access.

From debugger:

skirt`double RADAGAST::TemplatedUtils::integrate<double, std::__1::valarray<double>, std::__1::valarray<double> >(xContainer=0x00007ffeefbfd0c0, yContainer=0x00007ffeefbfcbb8, iMin=100, iMax=99) at TemplatedUtils.tpp:107:86
   104 	                auto iy = std::begin(yContainer) + iMin + 1;
   105 	                // One past the last point
   106 	                auto ixEnd = std::begin(xContainer) + iMax + 1;
-> 107 	                for (; ix != ixEnd; ix++, iy++) answer += 0.5 * (*ix - *(ix - 1)) * (*iy + *(iy - 1));

Full backtrace

* thread #1, queue = 'com.apple.main-thread', stop reason = EXC_BAD_ACCESS (code=1, address=0x107400000)
    frame #0: 0x00000001006772ed skirt`double RADAGAST::TemplatedUtils::integrate<double, std::__1::valarray<double>, std::__1::valarray<double> >(xContainer=0x00007ffeefbfd0c0, yContainer=0x00007ffeefbfcbb8, iMin=100, iMax=99) at TemplatedUtils.tpp:107:86
  * frame #1: 0x00000001006c352b skirt`RADAGAST::Ionization::photoRateCoeff(meanIntensity=0x00007ffeefbfd0a8) at Ionization.cpp:90:25
    frame #2: 0x00000001006bc4d2 skirt`RADAGAST::HModel::HModel(this=0x0000000102517d60, hData=0x000000010251baf0, meanIntensity=0x00007ffeefbfd0a8) at HModel.cpp:16:27
    frame #3: 0x00000001006bc565 skirt`RADAGAST::HModel::HModel(this=0x0000000102517d60, hData=0x000000010251baf0, meanIntensity=0x00007ffeefbfd0a8) at HModel.cpp:15:5
    frame #4: 0x000000010075dbb7 skirt`std::__1::__unique_if<RADAGAST::HModel>::__unique_single std::__1::make_unique<RADAGAST::HModel, RADAGAST::HData*, RADAGAST::Spectrum const*&>(__args=0x00007ffeefbfccf0, __args=0x00007ffeefbfccf8) at memory:3033:32
    frame #5: 0x000000010075db43 skirt`RADAGAST::SpeciesModelManager::makeHModel(this=0x000000010251a2e0, meanIntensity=0x00007ffeefbfd0a8) const at SpeciesModelManager.cpp:22:16
    frame #6: 0x000000010064520a skirt`RADAGAST::GasInterfaceImpl::makeGasSolution(this=0x000000010251a1f0, meanIntensity=0x00007ffeefbfd0a8, gri=0x000000010300be78) const at GasInterfaceImpl.cpp:446:47
    frame #7: 0x0000000100643a54 skirt`RADAGAST::GasInterfaceImpl::solveTemperature(this=0x000000010251a1f0, n=2.4029914185113967, meanIntensity=0x00007ffeefbfd0a8, gri=0x000000010300be78) const at GasInterfaceImpl.cpp:203:25
    frame #8: 0x000000010064394c skirt`RADAGAST::GasInterfaceImpl::updateGasState(this=0x000000010251a1f0, gs=0x000000010302d000, n=2.4029914185113967, meanIntensityv=0x00007ffeefbfd430, grainInfo=0x000000010300be78, gd=0x0000000000000000) const at GasInterfaceImpl.cpp:56:25
    frame #9: 0x000000010063f9e7 skirt`RADAGAST::GasInterface::updateGasState(this=0x0000000102510c10, gs=0x000000010302d000, n=2.4029914185113967, meanIntensityv=0x00007ffeefbfd430, gri=0x000000010300be78, gd=0x0000000000000000) const at GasInterface.cpp:32:17
    frame #10: 0x00000001005d0170 skirt`(anonymous namespace)::updateGasState_impl(m=0, n=2402991.4185113967, meanIntensityv=0x00007ffeefbfd6f0, mixNumberDensv=0x00007ffeefbfd5f8, gasDiagnostics=0x0000000000000000) at Gas.cpp:324:14
    frame #11: 0x00000001005cfd71 skirt`Gas::updateGasState(m=0, n=2402991.4185113967, meanIntensityv=0x00007ffeefbfd6f0, mixNumberDensv=0x00007ffeefbfd5f8) at Gas.cpp:369:5
    frame #12: 0x000000010019a051 skirt`GasRecipe::update(this=0x000000010250f4b0, state=0x00007ffeefbfd6b8, Jv=0x00007ffeefbfd6f0) at GasRecipe.cpp:121:5
    frame #13: 0x00000001001fe333 skirt`MediumSystem::updateDynamicMediumState(this=0x0000000102510cb8, firstIndex=0, numIndices=100)::$_2::operator()(unsigned long, unsigned long) const at MediumSystem.cpp:922:49
    frame #14: 0x00000001001fe1aa skirt`decltype(__f=0x0000000102510cb8, __args=0x00007ffeefbfd8b0, __args=0x00007ffeefbfd8a8)::$_2&>(fp)(std::__1::forward<unsigned long>(fp0), std::__1::forward<unsigned long>(fp0))) std::__1::__invoke<MediumSystem::updateDynamicMediumState()::$_2&, unsigned long, unsigned long>(MediumSystem::updateDynamicMediumState()::$_2&, unsigned long&&, unsigned long&&) at type_traits:3545:1
    frame #15: 0x00000001001fe137 skirt`void std::__1::__invoke_void_return_wrapper<void>::__call<MediumSystem::updateDynamicMediumState(__args=0x0000000102510cb8, __args=0x00007ffeefbfd8b0, __args=0x00007ffeefbfd8a8)::$_2&, unsigned long, unsigned long>(MediumSystem::updateDynamicMediumState()::$_2&, unsigned long&&, unsigned long&&) at __functional_base:348:9
    frame #16: 0x00000001001fe0e7 skirt`std::__1::__function::__alloc_func<MediumSystem::updateDynamicMediumState()::$_2, std::__1::allocator<MediumSystem::updateDynamicMediumState()::$_2>, void (unsigned long, unsigned long)>::operator(this=0x0000000102510cb8, __arg=0x00007ffeefbfd8b0, __arg=0x00007ffeefbfd8a8)(unsigned long&&, unsigned long&&) at functional:1546:16
    frame #17: 0x00000001001fcff8 skirt`std::__1::__function::__func<MediumSystem::updateDynamicMediumState()::$_2, std::__1::allocator<MediumSystem::updateDynamicMediumState()::$_2>, void (unsigned long, unsigned long)>::operator(this=0x0000000102510cb0, __arg=0x00007ffeefbfd8b0, __arg=0x00007ffeefbfd8a8)(unsigned long&&, unsigned long&&) at functional:1720:12
    frame #18: 0x0000000100295caa skirt`std::__1::__function::__value_func<void (unsigned long, unsigned long)>::operator(this=0x00007ffeefbfda80, __args=0x00007ffeefbfd8b0, __args=0x00007ffeefbfd8a8)(unsigned long&&, unsigned long&&) const at functional:1873:16
    frame #19: 0x0000000100295c35 skirt`std::__1::function<void (unsigned long, unsigned long)>::operator(this= Lambda in File MediumSystem.cpp at Line 908, __arg=0, __arg=100)(unsigned long, unsigned long) const at functional:2548:12
    frame #20: 0x0000000100295be0 skirt`SerialParallel::call(this=0x0000000102604150, maxIndex=100, target= Lambda in File MediumSystem.cpp at Line 908)>) at SerialParallel.cpp:17:19
    frame #21: 0x00000001001f27a5 skirt`MediumSystem::updateDynamicMediumState(this=0x000000010250d690) at MediumSystem.cpp:908:36
    frame #22: 0x0000000100203971 skirt`MonteCarloSimulation::runPrimaryEmissionWithDynamicState(this=0x0000000102504a40) at MonteCarloSimulation.cpp:172:41
    frame #23: 0x0000000100203356 skirt`MonteCarloSimulation::runSimulation(this=0x0000000102504a40) at MonteCarloSimulation.cpp:67:13
    frame #24: 0x0000000100296fb3 skirt`Simulation::setupAndRun(this=0x0000000102504a40) at Simulation.cpp:36:5
    frame #25: 0x0000000100007d69 skirt`SkirtCommandLineHandler::doSimulation(this=0x00007ffeefbff068, index=0) at SkirtCommandLineHandler.cpp:368:25
    frame #26: 0x00000001000045a1 skirt`SkirtCommandLineHandler::doBatch(this=0x00007ffeefbff068) at SkirtCommandLineHandler.cpp:185:9
    frame #27: 0x00000001000037ab skirt`SkirtCommandLineHandler::perform(this=0x00007ffeefbff068) at SkirtCommandLineHandler.cpp:58:42
    frame #28: 0x00000001000121b7 skirt`main(argc=4, argv=0x00007ffeefbff238) at SkirtMain.cpp:30:20
    frame #29: 0x00007fff6972acc9 libdyld.dylib`start + 1
    frame #30: 0x00007fff6972acc9 libdyld.dylib`start + 1
(lldb) 

A reasonable solution would be to return 0 in photoRateCoeff if the wavelength grid does not cover the right range.

Instructions to compile when GSL is installed in custom location

Need to find a way to specify where GSL is installed. Ideally through an option in Cmake

C+ cooling

To get realistic temperatures, I need carbon cooling in some shape or form. The C+ abundance is needed for that Different options for recipes:

• Calculate C+ abundance with chemical network
• Fixed C+ to HII, C+ to H, C+ to H2 ratios?
• Some kind of fit based on a grid PDR models, as a function of radiation field and density?

Additionally, O also has an important cooling effect, but with just C+ it will be in the right ball park already.

LineProfile crashes when sigmaGauss is zero

In practice, this occurs when T = 0, and the opacity calculation of the lines (onto a coarse grid) is invoked. While the line profile ()-operator will return a purely lorentzian shape in that case, the problem lies with the recommendedFrequencyGrid function, which makes use of 1/sigmaGauss. The resulting nans eventually cause a segfault because indices get confused.

Expected behavior: While setting T = 0 is never a good idea, LineProfile::addToBinned should not crash. Simply returning 0, and issuing a warning would be better.