Trouble with radiation wavelength grid for simulating CO

[swallstrom]

Hello! I'm trying to simulate a stellar wind with CO and dust emission, so the radiation field runs from 0.1 to 3000 microns. I'm getting errors saying the resolution of the wavelength grid is too small, and even with the nested logarithmic grid I'm having to increase the number of wavelengths to numWavelengthsSubGrid="1000000", at which point SKIRT instead crashes during the setup. How do I deal with narrow CO lines while still having emission from the (cool) star and dust? The errors and my full .ski file are copied below.

Wavelength grid error: 30/04/2024 14:26:50.326 * * Error: Integral of Gaussian line profile over radiation field is inaccurate: 30/04/2024 14:26:50.326 Error: integral equals 0.5316066542 rather than unity 30/04/2024 14:26:50.326 * * Error: over wavelengths from 2600.738724 micron to 2600.776547 micron 30/04/2024 14:26:50.326 Error: --> increase the resolution of the radiation field wavelength grid 30/04/2024 14:26:50.326 * * Error: On line 461 in file /Users/sofiaw/git/SKIRT/git/SKIRT/core/NonLTELineGasMix.cpp 30/04/2024 14:26:50.326 Error: In function updateSpecificState 30/04/2024 14:26:50.326 * * Error: Call stack: 30/04/2024 14:26:50.326 Error: 1 skirt 0x000000010b794f96 _ZN10FatalErrorC2ENSt3112basic_stringIcNS0_11char_traitsIcEENS09allocatorIcEEEEPKciS8 + 1158 30/04/2024 14:26:50.326 * *** Error: 2 skirt 0x000000010b85562b _ZNK16NonLTELineGasMix19updateSpecificStateEP13MaterialStateRKNSt318valarrayIdEE + 7659 30/04/2024 14:26:50.326 * * Error: 3 skirt 0x000000010b83c303 _ZNSt3110function6__funcIZN12MediumSystem23updateDynamicStateMediaEbE3$_4NS_9allocatorIS3EEFvmmEEclEOmS8 + 723 30/04/2024 14:26:50.326 Error: 4 skirt 0x000000010baa45c9 _ZN10ChunkMaker11callForNextERKNSt318functionIFvmmEEE + 73 30/04/2024 14:26:50.326 * * Error: 5 skirt 0x000000010b84b028 _ZN13MultiParallel3runEi + 328 30/04/2024 14:26:50.326 Error: 6 skirt 0x000000010b84ba56 _ZNSt3114__thread_proxyB6v15006INS_5tupleIJNS_10unique_ptrINS_15__thread_structENS_14default_deleteIS3_EEEEM13MultiParallelFviEPS7iEEEEEPvSC + 70 30/04/2024 14:26:50.326 * * Error: 7 libsystem_pthread.dylib 0x00007ff8129f01d3 _pthread_start + 125 30/04/2024 14:26:50.327 Error: 8 libsystem_pthread.dylib 0x00007ff8129ebbd3 thread_start + 15

SKIRT crashing: 30/04/2024 14:37:33.177 Calculating medium properties for 100 cells... 30/04/2024 14:37:33.177 Done calculating cell densities Killed: 9

ski file: <?xml version="1.0" encoding="UTF-8"?>

[petercamps]

This is indeed a tricky problem. Each of the relevant CO lines must be properly resolved, but using that spectral resolution across the full wavelength range is not feasible. Assuming 1 million spatial grid cells with 1 million wavelengths, storing the radiation field would need 8 x 10^6 x 10^6 = 10^15 bytes or 8 TB of memory. Undoubtedly, SKIRT crashes when trying to allocate this array.

The solution is to construct a custom wavelength grid and load it in SKIRT using the FileWavelengthGrid class, placing grid points in the appropriate places. In addition, you need to instruct SKIRT to favor shooting photon packets in the wavelength ranges of interest using a FileWavelengthDistribution.

@Koseimatsu should be able to offer further guidance and provide an example. I'll assign this issue to him as well.

[Koseimatsu]

Hello, thank you for your question. To solve this issue, you have to prepare grids with wide wavelength widths for dust continuum and grids with narrow wavelength widths for CO lines. For example, you can discretize the wavelength grids in the wavelength range from 100 to 3000 um by 100 wavelength widths. Then, you can discretize the wavelength grids in the wavelength range from -30 to 30 km/s around each CO line rest wavelength. Here, you have to be careful to ensure that the wavelength width around each CO line is sufficiently small compared to the microturbulence or thermal velocity. There are two ways to provide grids with different wavelength widths. 1. As Peter suggested, you can construct a custom wavelength grid and load it in SKIRT using the FileWavelengthGrid class. In a benchmark of the SKIRT, you can find the example (van Zadelhoff problem 2; [https://skirt.ugent.be/root/_zadelhoff2002.html]()). 2. You can also use a CompositeWavelengthGrid class. The advantage of this class is that you can define the wavelength grids in the ski file.

                        <radiationFieldWLG type="DisjointWavelengthGrid">
                            <CompositeWavelengthGrid log="false">
                                <wavelengthGrids type="DisjointWavelengthGrid">
                                    <LogWavelengthGrid minWavelength="100 micron" maxWavelength="3000 micron" numWavelengths="100"/>
                                    <LogWavelengthGrid minWavelength="2.602143e+02 micron" maxWavelength="2.602663e+02 micron" numWavelengths="100"/>
                                    <LogWavelengthGrid minWavelength="2.890908e+02 micron" maxWavelength="2.891487e+02 micron" numWavelengths="100"/>
                                    <LogWavelengthGrid minWavelength="3.251927e+02 micron" maxWavelength="3.252578e+02 micron" numWavelengths="100"/>
                                    <LogWavelengthGrid minWavelength="3.716140e+02 micron" maxWavelength="3.716884e+02 micron" numWavelengths="100"/>
                                    <LogWavelengthGrid minWavelength="4.335115e+02 micron" maxWavelength="4.335983e+02 micron" numWavelengths="100"/>
                                    <LogWavelengthGrid minWavelength="5.201807e+02 micron" maxWavelength="5.202848e+02 micron" numWavelengths="100"/>
                                    <LogWavelengthGrid minWavelength="6.501870e+02 micron" maxWavelength="6.503171e+02 micron" numWavelengths="100"/>
                                    <LogWavelengthGrid minWavelength="8.668706e+02 micron" maxWavelength="8.670441e+02 micron" numWavelengths="100"/>
                                    <LogWavelengthGrid minWavelength="1.300279e+03 micron" maxWavelength="1.300539e+03 micron" numWavelengths="100"/>
                                    <LogWavelengthGrid minWavelength="2.600497e+03 micron" maxWavelength="2.601018e+03 micron" numWavelengths="100"/>
                               </wavelengthGrids>
                            </CompositeWavelengthGrid>
                        </radiationFieldWLG>

I also note that it’s better to use these types of wavelength grid classes for the wavelength grids of the instrument system and wavelength bias distribution. I have attached an example ski file for the SKIRT simulation with CO and dust emission, which is used in Section 4 of Matsumoto et al. (2023). Wada2016.ski.txt

[swallstrom]

Thank you very much Kosei for the advice and example ski file! I had to add grids for all the CO lines, but the code runs properly now.