Add Nested Density Tree Policy

SKIRT/core/DensityTreePolicy.cpp

@@ -20,6 +20,8 @@
 
 void DensityTreePolicy::setupSelfBefore()
 {
+    TreePolicy::setupSelfBefore();
+
     // get the random generator and the number of samples (not a big effort, so we do this even if we don't need it)
     _random = find<Random>();
     _numSamples = find<Configuration>()->numDensitySamples();
@@ -85,6 +87,10 @@ void DensityTreePolicy::setupSelfBefore()
 
 bool DensityTreePolicy::needsSubdivide(TreeNode* node)
 {
+    // handle minimum and maximum level
+    if (node->level() < minLevel()) return true;
+    if (node->level() >= maxLevel()) return false;
+
     // results for the sampled mass or number densities, if applicable
     double rho = 0.;          // dust mass density
     double rhomin = DBL_MAX;  // smallest sample for dust mass density
@@ -179,28 +185,13 @@ vector<TreeNode*> DensityTreePolicy::constructTree(TreeNode* root)
     // initialize the tree node list with the root node as the first item
     vector<TreeNode*> nodev{root};
 
-    // recursively subdivide the root node until the minimum level has been reached
+    // initialize iteration variables to level 0
     int level = 0;    // current level
     size_t lbeg = 0;  // node index range for the current level;
     size_t lend = 1;  // at level 0, the node list contains just the root node
-    while (level != minLevel())
-    {
-        log->info("Subdividing level " + std::to_string(level) + ": " + std::to_string(lend - lbeg) + " nodes");
-        log->infoSetElapsed(lend - lbeg);
-        for (size_t l = lbeg; l != lend; ++l)
-        {
-            nodev[l]->subdivide(nodev);
-            if ((l + 1) % logDivideChunkSize == 0)
-                log->infoIfElapsed("Subdividing level " + std::to_string(level) + ": ", logDivideChunkSize);
-        }
-        // update iteration variables to the next level
-        level++;
-        lbeg = lend;
-        lend = nodev.size();
-    }
 
-    // recursively subdivide the nodes beyond the minimum level until all nodes satisfy the configured criteria
-    while (level != maxLevel() && lend != lbeg)
+    // recursively subdivide nodes until all nodes satisfy the configured criteria
+    while (lend != lbeg)
     {
         size_t numEvalNodes = lend - lbeg;
         log->info("Subdividing level " + std::to_string(level) + ": " + std::to_string(numEvalNodes) + " nodes");
@@ -238,6 +229,7 @@ vector<TreeNode*> DensityTreePolicy::constructTree(TreeNode* root)
                     log->infoIfElapsed("Subdivision for level " + std::to_string(level) + ": ", logDivideChunkSize);
             }
         }
+
         // update iteration variables to the next level
         level++;
         lbeg = lend;

SKIRT/core/DensityTreePolicy.hpp

@@ -123,26 +123,23 @@ class DensityTreePolicy : public TreePolicy, public MaterialWavelengthRangeInter
         evaluate the configured criteria. */
     void setupSelfBefore() override;
 
-private:
+public:
     /** This function returns true if the given node needs to be subdivided according to the
-        criteria configured for this policy, and false otherwise. The minimum and maximum level are
-        not checked, because this function is never called for nodes that don't conform to the
-        level criteria. */
-    bool needsSubdivide(TreeNode* node);
+        criteria configured for this policy, including minimum and maximum level, and false
+        otherwise. */
+    virtual bool needsSubdivide(TreeNode* node);
 
-public:
     /** This function constructs the hierarchical tree and all (interconnected) nodes forming the
         tree as described for the corresponding pure virtual function in the base class. The
-        implementation for this class loops over the tree subdivision levels (up to the maximum
-        level configured in the TreePolicy base class). For each level, the function alternates
-        between evaluating all of the nodes (i.e. determining which nodes need subdivision) and
-        actually subdividing the nodes that need it.
+        implementation for this class loops over the tree subdivision levels. For each level, the
+        function alternates between evaluating all of the nodes (i.e. determining which nodes need
+        subdivision) and actually subdividing the nodes that need it.
 
         These operations are split over two phases because the first one can be parallelized (the
         only output is a Boolean flag), while the second one cannot (the tree structure is updated
         in various ways). Parallelizing the first operation is often meaningful, because
-        determining whether a node needs subdivision can be resource-intensive (for example, it may
-        require sampling densities in the source distribution). */
+        determining whether a node needs subdivision can be resource-intensive. For example, it may
+        require sampling densities in the source distribution. */
     vector<TreeNode*> constructTree(TreeNode* root) override;
 
     //======================== Other Functions =======================
@@ -177,7 +174,7 @@ class DensityTreePolicy : public TreePolicy, public MaterialWavelengthRangeInter
     bool _hasElectronFraction{false};
     bool _hasGasFraction{false};
 
-    // cashed values for each material type (valid if corresponding flag is enabled)
+    // cached values for each material type (valid if corresponding flag is enabled)
     double _dustMass{0.};
     double _dustKappa{0.};
     double _electronNumber{0.};

SKIRT/core/NestedDensityTreePolicy.cpp

@@ -0,0 +1,42 @@
+/*//////////////////////////////////////////////////////////////////
+////     The SKIRT project -- advanced radiative transfer       ////
+////       © Astronomical Observatory, Ghent University         ////
+///////////////////////////////////////////////////////////////// */
+
+#include "NestedDensityTreePolicy.hpp"
+#include "FatalError.hpp"
+#include "Log.hpp"
+#include "SpatialGrid.hpp"
+#include "TreeNode.hpp"
+
+////////////////////////////////////////////////////////////////////
+
+void NestedDensityTreePolicy::setupSelfBefore()
+{
+    DensityTreePolicy::setupSelfBefore();
+
+    // verify that the inner box is not empty
+    if (_innerMaxX <= _innerMinX) throw FATALERROR("The extent of the inner box should be positive in the X direction");
+    if (_innerMaxY <= _innerMinY) throw FATALERROR("The extent of the inner box should be positive in the Y direction");
+    if (_innerMaxZ <= _innerMinZ) throw FATALERROR("The extent of the inner box should be positive in the Z direction");
+
+    // copy the coordinates to a Box instance for ease of use
+    _inner = Box(_innerMinX, _innerMinY, _innerMinZ, _innerMaxX, _innerMaxY, _innerMaxZ);
+
+    // verify that the inner box is inside the spatial grid
+    auto grid = find<SpatialGrid>(false);
+    if (!grid->boundingBox().contains(_inner))
+        find<Log>()->warning("The nested density tree policy region is not fully inside the spatial grid domain");
+}
+
+////////////////////////////////////////////////////////////////////
+
+bool NestedDensityTreePolicy::needsSubdivide(TreeNode* node)
+{
+    if (_inner.intersects(node->extent()))
+        return _innerPolicy->needsSubdivide(node);
+    else
+        return DensityTreePolicy::needsSubdivide(node);
+}
+
+////////////////////////////////////////////////////////////////////

SKIRT/core/NestedDensityTreePolicy.hpp

@@ -0,0 +1,88 @@
+/*//////////////////////////////////////////////////////////////////
+////     The SKIRT project -- advanced radiative transfer       ////
+////       © Astronomical Observatory, Ghent University         ////
+///////////////////////////////////////////////////////////////// */
+
+#ifndef NESTEDDENSITYTREEPOLICY_HPP
+#define NESTEDDENSITYTREEPOLICY_HPP
+
+#include "Box.hpp"
+#include "DensityTreePolicy.hpp"
+
+/** NestedDensityTreePolicy is a DensityTreePolicy policy that allows defining separate subdivision
+    criteria in a given subregion of the spatial grid. This can be used, for example, to specify a
+    higher resolution in a given region of interest.
+
+    The class inherits from DensityTreePolicy and uses all the inherited properties to classify the
+    outer region. Additional \em innerMin/innerMax coordinate properties define the bounding box of
+    the inner region. Finally, the additional property \em innerPolicy, which is also expected to
+    be a DensityTreePolicy instance, specifies the subdivision criteria for the inner region.
+
+    If a TreeNode intersects with the inner region the node will be subdivided based on the
+    criteria defined by the \em innerPolicy. This means that even TreeNodes that are almost fully
+    outside the inner region can be subdivided by those criteria if they have a non-zero
+    intersection with the inner region.
+
+    It is not meaningful to specify an inner region that extends outside of the spatial grid
+    bounding box, because none of the tree nodes will intersect those outside areas. Therefore, a
+    warning is issued if the inner region is not fully inside the spatial grid domain.
+
+    <b>Recursive nesting</b>
+
+    By default, the inner policy is a regular DensityTreePolicy instance. However, because
+    NestedDensityTreePolicy inherits DensityTreePolicy, it is also possible to again select a
+    NestedDensityTreePolicy instance as the inner policy, leading to recursive nesting. While this
+    is not a recommended use case, it would allow specifying recursively increasing resolution in
+    nested, successively smaller regions of the spatial domain. */
+class NestedDensityTreePolicy : public DensityTreePolicy
+{
+    ITEM_CONCRETE(NestedDensityTreePolicy, DensityTreePolicy,
+                  "a tree grid construction policy using a nested density tree policy")
+        ATTRIBUTE_TYPE_DISPLAYED_IF(NestedDensityTreePolicy, "Level2")
+
+        PROPERTY_DOUBLE(innerMinX, "the start point of the inner box in the X direction")
+        ATTRIBUTE_QUANTITY(innerMinX, "length")
+
+        PROPERTY_DOUBLE(innerMaxX, "the end point of the inner box in the X direction")
+        ATTRIBUTE_QUANTITY(innerMaxX, "length")
+
+        PROPERTY_DOUBLE(innerMinY, "the start point of the inner box in the Y direction")
+        ATTRIBUTE_QUANTITY(innerMinY, "length")
+
+        PROPERTY_DOUBLE(innerMaxY, "the end point of the inner box in the Y direction")
+        ATTRIBUTE_QUANTITY(innerMaxY, "length")
+
+        PROPERTY_DOUBLE(innerMinZ, "the start point of the inner box in the Z direction")
+        ATTRIBUTE_QUANTITY(innerMinZ, "length")
+
+        PROPERTY_DOUBLE(innerMaxZ, "the end point of the inner box in the Z direction")
+        ATTRIBUTE_QUANTITY(innerMaxZ, "length")
+
+        PROPERTY_ITEM(innerPolicy, DensityTreePolicy, "the density tree policy for the inner region")
+        ATTRIBUTE_DEFAULT_VALUE(innerPolicy, "DensityTreePolicy")
+
+    ITEM_END()
+
+    //============= Construction - Setup - Destruction =============
+
+protected:
+    /** This function checks whether the inner region is inside the outer region. */
+    void setupSelfBefore() override;
+
+    //======================== Other Functions =======================
+
+public:
+    /** This function determines whether the given node needs to be subdivided. Depending on
+        whether the node intersects with the inner region or not, the request is passed to the
+        needsSubdivide() function of the nested policy or to the needsSubdivide() function of our
+        base class. */
+    bool needsSubdivide(TreeNode* node) override;
+
+    //======================== Data Members ========================
+
+private:
+    // the inner region box
+    Box _inner;
+};
+
+#endif

SKIRT/core/SimulationItemRegistry.cpp

@@ -172,6 +172,7 @@
 #include "MollweideProjection.hpp"
 #include "MonteCarloSimulation.hpp"
 #include "MultiGaussianExpansionGeometry.hpp"
+#include "NestedDensityTreePolicy.hpp"
 #include "NestedLogWavelengthGrid.hpp"
 #include "NetzerAngularDistribution.hpp"
 #include "NoPolarizationProfile.hpp"
@@ -484,6 +485,7 @@ SimulationItemRegistry::SimulationItemRegistry(string version, string format)
     // spatial grid policies
     ItemRegistry::add<TreePolicy>();
     ItemRegistry::add<DensityTreePolicy>();
+    ItemRegistry::add<NestedDensityTreePolicy>();
     ItemRegistry::add<SiteListTreePolicy>();
 
     // one-dimensional meshes for spatial grids

SKIRT/utils/Box.hpp

@@ -108,6 +108,13 @@ class Box
         return x >= _xmin && x <= _xmax && y >= _ymin && y <= _ymax && z >= _zmin && z <= _zmax;
     }
 
+    /** This function returns true if the given box is inside this box, false otherwise. */
+    inline bool contains(const Box& box) const
+    {
+        return xmin() <= box.xmin() && xmax() >= box.xmax() && ymin() <= box.ymin() && ymax() >= box.ymax()
+               && zmin() <= box.zmin() && zmax() >= box.zmax();
+    }
+
     /** This function returns the volume \f$(x_\text{max}-x_\text{min}) \times
         (y_\text{max}-y_\text{min}) \times (z_\text{max}-z_\text{min})\f$ of the box. */
     inline double volume() const { return (_xmax - _xmin) * (_ymax - _ymin) * (_zmax - _zmin); }
@@ -159,6 +166,16 @@ class Box
         k = std::max(0, std::min(nz - 1, static_cast<int>(nz * (r.z() - _zmin) / (_zmax - _zmin))));
     }
 
+    /** This function returns true if the given box and this box have a non-zero intersection,
+        false otherwise. */
+    inline bool intersects(const Box& box) const
+    {
+        if (xmax() < box.xmin() || box.xmax() < xmin()) return false;
+        if (ymax() < box.ymin() || box.ymax() < ymin()) return false;
+        if (zmax() < box.zmin() || box.zmax() < zmin()) return false;
+        return true;
+    }
+
     /** This function intersects the receiving axis-aligned bounding box with a ray (half-line)
         defined by the specified starting position \f$\bf{r}\f$ and direction \f$\bf{k}\f$. If the
         ray intersects the box, the function returns true after storing the near and far