search

SKIRT / SKIRT9

SKIRT version 9 -- advanced radiative transfer in dusty systems
http://www.skirt.ugent.be
GNU Affero General Public License v3.0
37 stars 31 forks source link

Unable to reach convergence between input and gridded mass #209

Closed GervitKT closed 7 months ago

GervitKT commented 7 months ago

Description I am getting 99% discrepancy between the input and the gridded medium mass. I have tried doing this with both particle media and cell media but the problem persists. The distribution of cells as indicated by the log files seems pretty reasonable. Is there any way to judge why the mass is getting lost, and can SKIRT also produce the residual maps for medium density directly, so as to see the regions in which the mass is lost.

Data I upload here log file for cell media as that is my preferred choice and also the convergence plots for both particle and cell media. gal503_standard_cell_media_log.txt gal503_standard_dns_dust_cell_media.pdf gal503_standard_dns_dust_particle_media.pdf

and this is the ski file:

<?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 e790b9f built on 01/03/2024 at 13:50:55)" time="2024-03-15T16:36:53.938">
    <MonteCarloSimulation userLevel="Regular" simulationMode="DustEmission" iteratePrimaryEmission="false" iterateSecondaryEmission="false" numPackets="1e7">
        <random type="Random">
            <Random seed="0"/>
        </random>
        <units type="Units">
            <ExtragalacticUnits wavelengthOutputStyle="Wavelength" fluxOutputStyle="Wavelength"/>
        </units>
        <cosmology type="Cosmology">
            <LocalUniverseCosmology/>
        </cosmology>
        <sourceSystem type="SourceSystem">
            <SourceSystem minWavelength="0.01 micron" maxWavelength="1000 micron" wavelengths="0.55 micron" sourceBias="0.5">
                <sources type="Source">
                    <ParticleSource filename="/vol/aibn1089/data1/gtrehan/Skirt input files/input_files/gal503/skirt_src_infile_v9.txt" importVelocity="true" importVelocityDispersion="false" importCurrentMass="false" useColumns="" sourceWeight="1" wavelengthBias="0.5">
                        <smoothingKernel type="SmoothingKernel">
                            <CubicSplineSmoothingKernel/>
                        </smoothingKernel>
                        <sedFamily type="SEDFamily">
                            <FSPSSEDFamily imf="Salpeter"/>
                        </sedFamily>
                        <wavelengthBiasDistribution type="WavelengthDistribution">
                            <DefaultWavelengthDistribution/>
                        </wavelengthBiasDistribution>
                    </ParticleSource>
                </sources>
            </SourceSystem>
        </sourceSystem>
        <mediumSystem type="MediumSystem">
            <MediumSystem>
                <photonPacketOptions type="PhotonPacketOptions">
                    <PhotonPacketOptions explicitAbsorption="false" forceScattering="true" minWeightReduction="1e4" minScattEvents="0" pathLengthBias="0.5"/>
                </photonPacketOptions>
                <radiationFieldOptions type="RadiationFieldOptions">
                    <RadiationFieldOptions storeRadiationField="true">
                        <radiationFieldWLG type="DisjointWavelengthGrid">
                            <LogWavelengthGrid minWavelength="0.01 micron" maxWavelength="1000 micron" numWavelengths="500"/>
                        </radiationFieldWLG>
                    </RadiationFieldOptions>
                </radiationFieldOptions>
                <secondaryEmissionOptions type="SecondaryEmissionOptions">
                    <SecondaryEmissionOptions storeEmissionRadiationField="false" secondaryPacketsMultiplier="1" spatialBias="0.5" sourceBias="0.5"/>
                </secondaryEmissionOptions>
                <dustEmissionOptions type="DustEmissionOptions">
                    <DustEmissionOptions dustEmissionType="Equilibrium" includeHeatingByCMB="false" maxFractionOfPrimary="0.01" maxFractionOfPrevious="0.03" sourceWeight="1" wavelengthBias="0.5">
                        <cellLibrary type="SpatialCellLibrary">
                            <AllCellsLibrary/>
                        </cellLibrary>
                        <dustEmissionWLG type="DisjointWavelengthGrid">
                            <LogWavelengthGrid minWavelength="0.01 micron" maxWavelength="1000 micron" numWavelengths="500"/>
                        </dustEmissionWLG>
                        <wavelengthBiasDistribution type="WavelengthDistribution">
                            <DefaultWavelengthDistribution/>
                        </wavelengthBiasDistribution>
                    </DustEmissionOptions>
                </dustEmissionOptions>
                <media type="Medium">
                    <CellMedium filename="/vol/aibn1089/data1/gtrehan/Skirt input files/input_files/gal503/skirt_dust_infile_v9.txt" massType="MassDensity" massFraction="1" importMetallicity="false" importTemperature="false" maxTemperature="0 K" importVelocity="true" importMagneticField="false" importVariableMixParams="false" useColumns="">
                        <materialMix type="MaterialMix">
                            <ZubkoDustMix numSilicateSizes="5" numGraphiteSizes="5" numPAHSizes="5"/>
                        </materialMix>
                    </CellMedium>
                </media>
                <samplingOptions type="SamplingOptions">
                    <SamplingOptions numDensitySamples="100" numPropertySamples="1" aggregateVelocity="Average"/>
                </samplingOptions>
                <grid type="SpatialGrid">
                    <PolicyTreeSpatialGrid minX="-5108.841 pc" maxX="5108.841 pc" minY="-5108.841 pc" maxY="5108.841 pc" minZ="-5108.841 pc" maxZ="5108.841 pc" treeType="OctTree">
                        <policy type="TreePolicy">
                            <DensityTreePolicy minLevel="0" maxLevel="20" maxDustFraction="1e-6" maxDustOpticalDepth="0" wavelength="0.55 micron" maxDustDensityDispersion="0" maxElectronFraction="1e-6" maxGasFraction="1e-6"/>
                        </policy>
                    </PolicyTreeSpatialGrid>
                </grid>
            </MediumSystem>
        </mediumSystem>
        <instrumentSystem type="InstrumentSystem">
            <InstrumentSystem>
                <defaultWavelengthGrid type="WavelengthGrid">
                    <LogWavelengthGrid minWavelength="1e-2 micron" maxWavelength="1e3 micron" numWavelengths="250"/>
                </defaultWavelengthGrid>
                <instruments type="Instrument">
                    <FullInstrument instrumentName="i00" distance="10 Mpc" inclination="0 deg" azimuth="0 deg" roll="90 deg" fieldOfViewX="10217.6829 pc" numPixelsX="800" centerX="0 pc" fieldOfViewY="10217.6829 pc" numPixelsY="800" centerY="0 pc" recordComponents="false" numScatteringLevels="0" recordPolarization="false" recordStatistics="false">
                        <wavelengthGrid type="WavelengthGrid">
                            <LogWavelengthGrid minWavelength="1e-2 micron" maxWavelength="1e3 micron" numWavelengths="250"/>
                        </wavelengthGrid>
                    </FullInstrument>
                    <FullInstrument instrumentName="i90" distance="10 Mpc" inclination="90 deg" azimuth="0 deg" roll="90 deg" fieldOfViewX="10217.6829 pc" numPixelsX="800" centerX="0 pc" fieldOfViewY="10217.6829 pc" numPixelsY="800" centerY="0 pc" recordComponents="false" numScatteringLevels="0" recordPolarization="false" recordStatistics="false">
                        <wavelengthGrid type="WavelengthGrid">
                            <LogWavelengthGrid minWavelength="1e-2 micron" maxWavelength="1e3 micron" numWavelengths="250"/>
                        </wavelengthGrid>
                    </FullInstrument>
                    <FullInstrument instrumentName="i45" distance="10 Mpc" inclination="45 deg" azimuth="0 deg" roll="90 deg" fieldOfViewX="10217.6829 pc" numPixelsX="800" centerX="0 pc" fieldOfViewY="10217.6829 pc" numPixelsY="800" centerY="0 pc" recordComponents="false" numScatteringLevels="0" recordPolarization="false" recordStatistics="false">
                        <wavelengthGrid type="WavelengthGrid">
                            <LogWavelengthGrid minWavelength="1e-2 micron" maxWavelength="1e3 micron" numWavelengths="250"/>
                        </wavelengthGrid>
                    </FullInstrument>
                    <FullInstrument instrumentName="i30" distance="10 Mpc" inclination="30 deg" azimuth="0 deg" roll="90 deg" fieldOfViewX="10217.6829 pc" numPixelsX="800" centerX="0 pc" fieldOfViewY="10217.6829 pc" numPixelsY="800" centerY="0 pc" recordComponents="false" numScatteringLevels="0" recordPolarization="false" recordStatistics="false">
                        <wavelengthGrid type="WavelengthGrid">
                            <LogWavelengthGrid minWavelength="1e-2 micron" maxWavelength="1e3 micron" numWavelengths="250"/>
                        </wavelengthGrid>
                    </FullInstrument>
                    <FullInstrument instrumentName="i70" distance="10 Mpc" inclination="70 deg" azimuth="0 deg" roll="90 deg" fieldOfViewX="31089.765191355 pc" numPixelsX="800" centerX="0 pc" fieldOfViewY="31089.765191355 pc" numPixelsY="800" centerY="0 pc" recordComponents="false" numScatteringLevels="0" recordPolarization="false" recordStatistics="false">
                        <wavelengthGrid type="WavelengthGrid">
                            <LogWavelengthGrid minWavelength="1e-2 micron" maxWavelength="1e3 micron" numWavelengths="250"/>
                        </wavelengthGrid>
                    </FullInstrument>
                </instruments>
            </InstrumentSystem>
        </instrumentSystem>
        <probeSystem type="ProbeSystem">
            <ProbeSystem>
                <probes type="Probe">
                    <ConvergenceInfoProbe probeName="cnv" wavelength="0.0001 micron" probeAfter="Setup"/>
                    <ConvergenceCutsProbe probeName="dns" probeAfter="Setup"/>
                    <SpatialCellPropertiesProbe probeName="prp" wavelength="0.0001 micron" probeAfter="Setup"/>
                    <DensityProbe probeName="dns" aggregation="Type" probeAfter="Setup">
                        <form type="Form">
                        <PerCellForm/>
                        </form>
                    </DensityProbe>
                    <TemperatureProbe probeName="tmp" aggregation="Type" probeAfter="Run">
                        <form type="Form">
                        <DefaultCutsForm/>
                        </form>
                    </TemperatureProbe>
                    <TemperatureProbe probeName="tmp" aggregation="Type" probeAfter="Run">
                        <form type="Form">
                        <PerCellForm/>
                        </form>
                    </TemperatureProbe>
                    <LuminosityProbe probeName="src_lum">
                        <wavelengthGrid type="WavelengthGrid">
                            <LogWavelengthGrid minWavelength="0.01 micron" maxWavelength="1e3 micron" numWavelengths="500"/>
                        </wavelengthGrid>
                    </LuminosityProbe>
                    <DensityProbe probeName="dns_med" aggregation="Type" probeAfter="Setup">
                        <form type="Form">
                            <DefaultCutsForm/>
                        </form>
                    </DensityProbe>
                    <SpatialGridPlotProbe probeName="grd"/>
                    <DustGrainPopulationsProbe probeName="dust_prop"/>
                    <DustGrainSizeDistributionProbe probeName="grain_size" numSamples="250"/>
                </probes>
            </ProbeSystem>
        </probeSystem>
    </MonteCarloSimulation>
</skirt-simulation-hierarchy>
petercamps commented 7 months ago

The input mass accumulates all particles/cells in the input file, regardless of their position. Your input file may thus contain particles/cells outside the grid boundary that are counted towards the input mass but are not gridded. Another possibility is that the input file contains extremely tiny (or even zero sized) particles/cells that are not resolved by the grid and thus "missed" by the random sampling mechanism.

SKIRT does not offer a way to visualize the residual, but you can easily write a python script to plot a histogram of the mass distribution in your input file as a function of radius (allocating all mass of each particle/cell to its center). This plot should tell you whether any large masses appear in unexpected areas, and whether SKIRT has captured the relevant mass or not.

If such an analysis does not resolve your issue, let us know., and we'll look into it further.

GervitKT commented 7 months ago

Thanks for the quick response. As you said, my input file does have cells outside the specified domain for the grid. I did not expect it and probably is due to some error in reading the cell positions from the simulation snapshot. Thanks for the help.