The use of CellMedium

[kpolsen]

Description This is a question about the new CellMedium class: https://skirt.ugent.be/skirt9/class_cell_medium.html

Is there a way to keep the cell structure of the original simulation, preventing SKIRT from "re-gridding" the input medium?

I'm trying to run SKIRT on a pretty big ENZO AMR simulation of a galaxy with ~29 million gas cells and ~10,000 stars. The initial run finished in around 2 hours with 1e7 photon packages. However, the total number of cells in the skirt output got reduced to 1,887,642, and I wonder if this is due to my skirt input file or not controlled by the user.

Data Below is my entire ski file - in particular, I wasn't sure what to choose under grid type="SpatialGrid":

<?xml version="1.0" encoding="UTF-8"?>
<!-- A SKIRT parameter file © Astronomical Observatory, Ghent University -->
<skirt-simulation-hierarchy type="MonteCarloSimulation" format="9" producer="SKIRT v9.0 (git b1ccc21-dirty built on 04/08/2020 at 11:57:27)" time="2020-08-04T12:12:37.681">
    <MonteCarloSimulation userLevel="Regular" simulationMode="DustEmissionWithSelfAbsorption" numPackets="1e7">
        <random type="Random">
            <Random seed="0"/>
        </random>
        <units type="Units">
            <ExtragalacticUnits fluxOutputStyle="Wavelength"/>
        </units>
        <cosmology type="Cosmology">
            <LocalUniverseCosmology/>
        </cosmology>
        <sourceSystem type="SourceSystem">
            <SourceSystem minWavelength="0.001 micron" maxWavelength="160 micron" wavelengths="0.55 micron" sourceBias="0.5">
                <sources type="Source">
                    <ParticleSource filename="enzo_15kpc.stars" importVelocity="false" importVelocityDispersion="false" useColumns="" sourceWeight="1" wavelengthBias="0.5">
                        <smoothingKernel type="SmoothingKernel">
                            <CubicSplineSmoothingKernel/>
                        </smoothingKernel>
                        <sedFamily type="SEDFamily">
                            <BruzualCharlotSEDFamily imf="Chabrier" resolution="Low"/>
                        </sedFamily>
                        <wavelengthBiasDistribution type="WavelengthDistribution">
                            <LogWavelengthDistribution minWavelength="0.0001 micron" maxWavelength="1e6 micron"/>
                        </wavelengthBiasDistribution>
                    </ParticleSource>
                </sources>
            </SourceSystem>
        </sourceSystem>
        <mediumSystem type="MediumSystem">
            <MediumSystem numDensitySamples="100">
                <photonPacketOptions type="PhotonPacketOptions">
                    <PhotonPacketOptions minWeightReduction="1e4" minScattEvents="0" pathLengthBias="0.5"/>
                </photonPacketOptions>
                <dustEmissionOptions type="DustEmissionOptions">
                    <DustEmissionOptions dustEmissionType="Stochastic" includeHeatingByCMB="false" storeEmissionRadiationField="false" secondaryPacketsMultiplier="1" spatialBias="0.5" wavelengthBias="0.5">
                        <cellLibrary type="SpatialCellLibrary">
                            <AllCellsLibrary/>
                        </cellLibrary>
                        <radiationFieldWLG type="DisjointWavelengthGrid">
                            <NestedLogWavelengthGrid minWavelengthBaseGrid="0.001 micron" maxWavelengthBaseGrid="160 micron" numWavelengthsBaseGrid="60" minWavelengthSubGrid="0.0912 micron" maxWavelengthSubGrid="0.207 micron" numWavelengthsSubGrid="25"/>
                        </radiationFieldWLG>
                        <dustEmissionWLG type="DisjointWavelengthGrid">
                            <NestedLogWavelengthGrid minWavelengthBaseGrid="1 micron" maxWavelengthBaseGrid="1000 micron" numWavelengthsBaseGrid="300" minWavelengthSubGrid="3 micron" maxWavelengthSubGrid="25 micron" numWavelengthsSubGrid="100"/>
                        </dustEmissionWLG>
                        <wavelengthBiasDistribution type="WavelengthDistribution">
                            <LogWavelengthDistribution minWavelength="0.0001 micron" maxWavelength="1e6 micron"/>
                        </wavelengthBiasDistribution>
                    </DustEmissionOptions>
                </dustEmissionOptions>
                <dustSelfAbsorptionOptions type="DustSelfAbsorptionOptions">
                    <DustSelfAbsorptionOptions minIterations="1" maxIterations="10" maxFractionOfPrimary="0.01" maxFractionOfPrevious="0.03" iterationPacketsMultiplier="1"/>
                </dustSelfAbsorptionOptions>
                <media type="Medium">
                    <CellMedium filename="enzo_15kpc.gas" massType="Mass" massFraction="0.2" importMetallicity="true" importTemperature="false" maxTemperature="0 K" importVelocity="false" importMagneticField="false" importVariableMixParams="false" useColumns="">
                        <materialMix type="MaterialMix">
                            <ThemisDustMix numSilicateSizes="15" numHydrocarbonSizes="15"/>
                        </materialMix>
                    </CellMedium>
                </media>
                <grid type="SpatialGrid">
                    <PolicyTreeSpatialGrid minX="-14963 pc" maxX="14963 pc" minY="-14963 pc" maxY="14963 pc" minZ="-4961 pc" maxZ="4961 pc" treeType="OctTree">
                        <policy type="TreePolicy">
                            <SiteListTreePolicy minLevel="3" maxLevel="7" numExtraLevels="0"/>
                        </policy>
                    </PolicyTreeSpatialGrid>
                </grid>
            </MediumSystem>
        </mediumSystem>
        <instrumentSystem type="InstrumentSystem">
            <InstrumentSystem>
                <defaultWavelengthGrid type="WavelengthGrid">
                    <LogWavelengthGrid minWavelength="0.001 micron" maxWavelength="160 micron" numWavelengths="750"/>
                </defaultWavelengthGrid>
                <instruments type="Instrument">
                    <SEDInstrument instrumentName="xy" distance="10 Mpc" inclination="0 deg" azimuth="0 deg" roll="90 deg" recordComponents="true" numScatteringLevels="0" recordPolarization="false" recordStatistics="false"/>
                </instruments>
            </InstrumentSystem>
        </instrumentSystem>
        <probeSystem type="ProbeSystem">
            <ProbeSystem>
                <probes type="Probe">
                    <RadiationFieldPerCellProbe probeName="rfpc" writeWavelengthGrid="true"/>
                    <DefaultRadiationFieldCutsProbe probeName="rfcut" writeWavelengthGrid="true"/>
                </probes>
            </ProbeSystem>
        </probeSystem>
    </MonteCarloSimulation>
</skirt-simulation-hierarchy>

System SKIRT v9, Linux, Ubuntu

Context It would be easier in post-processing if everything followed the original grid so I can compare with other cell properties in the simulation.

[petercamps]

With CellMedium, the structure of the imported mesh (i.e. the relationships between the neighbouring cells) is lost, so you cannot use the original grid for the radiative transfer. A new grid must be constructed in SKIRT, and this new grid will never be identical to the original grid.

If you want the SKIRT grid to be identical to the input grid, you need to use AdapativeMeshMedium. Because the cells are then supplied in Morton order, the mesh structure can be reconstructed and used for radiative transfer.

If you're OK with regridding, you can use CellMedium. The PolicyTreeSpatialGrid with OctTree tree type is a good choice. There are two options for the regridding policy. You could use the SiteListTreePolicy, as in your ski file. This builds an oct-tree grid based on the cell centers in the input mesh (the grid will approximate but not be identical to the original grid). Your ski file has a maximum of only 7 subdivision levels, which is why the SKIRT grid has fewer cells. You will need to increase this to 11 or 12 levels, I guess.

The second option is the DensityTreePolicy, which builds the oct-tree grid based on the dust density in the input distribution, regardless of the original gridding. This is a good option, but it does require fine-tuning of the parameters. See the tutorial on this subject.

[kpolsen]

I've decided to stick with CellMedium for now, however I wanted to know how to increase the number of levels used and in particular decrease the smallest grid level volume? It seems that the first 5 levels are not used, but changing minLevel to "5" doesn't seem to have an effect:

26/08/2020 15:52:58.902   Constructing intermediate 61x61x61 grid for cells...
26/08/2020 15:53:03.315     Smallest number of cells per grid cell: 3
26/08/2020 15:53:03.315     Largest  number of cells per grid cell: 4437
26/08/2020 15:53:03.315     Average  number of cells per grid cell: 250.0
26/08/2020 15:53:03.329   Constructing the spatial tree grid...
26/08/2020 15:53:03.329   Subdividing level 0: 1 nodes
26/08/2020 15:53:03.329   Subdividing level 1: 8 nodes
26/08/2020 15:53:03.329   Subdividing level 2: 64 nodes
26/08/2020 15:53:03.330   Subdividing tree to insert 28938614 sites
26/08/2020 15:53:08.331   Inserting site 22169999: 76.6%
26/08/2020 15:53:11.055   Finished construction of the spatial tree grid
26/08/2020 15:53:11.055   Number of cells at each level in the tree hierarchy:
26/08/2020 15:53:11.055     Level  0:        0 (  0.0%)  |
26/08/2020 15:53:11.055     Level  1:        0 (  0.0%)  |
26/08/2020 15:53:11.055     Level  2:        0 (  0.0%)  |
26/08/2020 15:53:11.055     Level  3:        0 (  0.0%)  |
26/08/2020 15:53:11.055     Level  4:        0 (  0.0%)  |
26/08/2020 15:53:11.055     Level  5:      781 (  0.0%)  |
26/08/2020 15:53:11.055     Level  6:    22901 (  1.2%)  |
26/08/2020 15:53:11.055     Level  7:  1863960 ( 98.7%)  |********************
26/08/2020 15:53:11.055     TOTAL   :  1887642 (100.0%)

Sorry to be returning to this!

[petercamps]

As I wrote in the previous response, you need to increase the maximum allowed subdivision level from 7 to a higher number. I recommend that you work through the tutorial indicated in the previous response. Even though the tutorial focuses on the density policy, it will give you some intuition that you can also use with the site list policy.