Implemented inner boundary for imaging

magritte/setup.py

@@ -419,7 +419,7 @@ def set_boundary_condition_zero (model):
         Updated Magritte object.
     """
     for b in range(model.parameters.nboundary()):
-        model.geometry.boundary.set_boundary_condition (b, BoundaryCondition.Zero)
+        model.geometry.boundary.set_boundary_condition (b, BoundaryCondition.Zero, False)
     # Done
     return model
 
@@ -439,12 +439,39 @@ def set_boundary_condition_CMB (model):
         Updated Magritte object.
     """
     for b in range(model.parameters.nboundary()):
-        model.geometry.boundary.set_boundary_condition (b, BoundaryCondition.CMB)
+        model.geometry.boundary.set_boundary_condition (b, BoundaryCondition.CMB, False)
     model.geometry.boundary.boundary_temperature.set([T_CMB for _ in range(model.parameters.nboundary())])
     # Done
     return model
 
 
+def set_inner_boundary_condition_thermal (model, T_in, bdy_end, bdy_start = 0):
+    """
+    Setter for incoming black body radiation boundary condition at the inner boundary point.
+
+    Parameters
+    ----------
+    model : Magritte model object
+        Magritte model object to set.
+    T_in : float
+        Boundary temperature at the inner boundary.
+    bdy_end : int
+        Boundary index of the last boundary point where the boundary condition should be set + 1.
+    bdy_start : int = 0
+        Boundary index of the first boundary point where the boundary condition should be set.
+    Returns
+    -------
+    out : Magritte model object
+        Updated Magritte object.
+    """
+    Tmp = np.array([T_CMB for _ in range(model.parameters.nboundary())])
+    for b in range(bdy_start, bdy_end):
+        Tmp[b] = T_in
+        model.geometry.boundary.set_boundary_condition (b, BoundaryCondition.Thermal, True)
+    model.geometry.boundary.boundary_temperature.set(Tmp)
+    return model
+
+
 def set_boundary_condition_1D (model, T_in=T_CMB, T_out=T_CMB):
     """
     Setter for incoming black body radiation boundary condition at each boundary point.
@@ -470,7 +497,7 @@ def set_boundary_condition_1D (model, T_in=T_CMB, T_out=T_CMB):
     else:
         # Set all boundary conditions to Thermal
         for b in range(model.parameters.nboundary()):
-            model.geometry.boundary.set_boundary_condition (b, BoundaryCondition.Thermal)
+            model.geometry.boundary.set_boundary_condition (b, BoundaryCondition.Thermal, False)
         # Set inner and outer temperature
         model.geometry.boundary.boundary_temperature.set([T_in, T_out])
     # Done

src/bindings/pybindings.cpp

@@ -453,7 +453,7 @@ PYBIND11_MODULE(core, module) {
             "(only relevant for thermal boundary conditions).")
         // functions
         .def("set_boundary_condition", &Boundary::set_boundary_condition,
-            "Setter for the boundary condition.")
+            "Setter for the boundary condition.", py::arg("b"), py::arg("cd"), py::arg("is_inner"))
         .def("get_boundary_condition", &Boundary::get_boundary_condition,
             "Getter for the boundary condition.")
         // io

src/model/geometry/boundary/boundary.cpp

@@ -28,6 +28,7 @@ void Boundary ::read(const Io& io) {
     // Set boundary conditions
     boundary_condition.resize(parameters->nboundary());
     boundary_temperature.resize(parameters->nboundary());
+    is_inner_boundary.resize(parameters->nboundary());
 
     Size1 boundary_condition_int(parameters->nboundary());
 
@@ -40,6 +41,16 @@ void Boundary ::read(const Io& io) {
         if (boundary_condition_int[b] == 2) boundary_condition[b] = CMB;
     }
 
+    // Older versions of magritte do not specify whether a boundary is an inner boundary. Therefore,
+    // we set all boundaries to be outer boundaries by default.
+    int err = io.read_list(prefix + "is_inner_boundary", is_inner_boundary);
+
+    if (err != 0) {
+        for (Size b = 0; b < parameters->nboundary(); b++) {
+            is_inner_boundary[b] = 0;
+        }
+    }
+
     boundary2point.copy_vec_to_ptr();
     point2boundary.copy_vec_to_ptr();
 
@@ -62,11 +73,15 @@ void Boundary ::write(const Io& io) const {
 
     io.write_list(prefix + "boundary_temperature", boundary_temperature);
     io.write_list(prefix + "boundary_condition", boundary_condition_int);
+    io.write_list(prefix + "is_inner_boundary", is_inner_boundary);
 }
 
-BoundaryCondition Boundary ::set_boundary_condition(const Size b, const BoundaryCondition cd) {
+BoundaryCondition Boundary ::set_boundary_condition(
+    const Size b, const BoundaryCondition cd, const bool is_inner) {
     boundary_condition.resize(parameters->nboundary());
     boundary_condition[b] = cd;
+    is_inner_boundary.resize(parameters->nboundary());
+    is_inner_boundary[b] = (int)is_inner;
 
     return boundary_condition[b];
 }

src/model/geometry/boundary/boundary.hpp

@@ -9,17 +9,19 @@ enum BoundaryCondition { Zero, Thermal, CMB };
 struct Boundary {
     std::shared_ptr<Parameters> parameters;
 
-    Vector<Size> boundary2point; // maps boundary index ∈[0, nboundary-1] to point
-                                 // index ∈[0, npoints-1]
-    Vector<Size> point2boundary; // maps point index to boundary index. Maps non-boundary
-                                 // points to parameters->npoints() (invalid value)
+    Vector<Size> boundary2point;    // maps boundary index ∈[0, nboundary-1] to point
+                                    // index ∈[0, npoints-1]
+    Vector<Size> point2boundary;    // maps point index to boundary index. Maps non-boundary
+                                    // points to parameters->npoints() (invalid value)
+    Vector<Size> is_inner_boundary; // 1 if boundary is an inner boundary, 0 otherwise
 
     Vector<BoundaryCondition> boundary_condition;
     Vector<Real> boundary_temperature;
 
     Boundary(std::shared_ptr<Parameters> params) : parameters(params){};
 
-    BoundaryCondition set_boundary_condition(const Size b, const BoundaryCondition cd);
+    BoundaryCondition set_boundary_condition(
+        const Size b, const BoundaryCondition cd, const bool is_inner);
     BoundaryCondition get_boundary_condition(const Size b) const;
 
     void read(const Io& io);

src/model/geometry/geometry.hpp

@@ -82,6 +82,7 @@ struct Geometry {
 
     inline bool valid_point(const Size p) const;
     inline bool not_on_boundary(const Size p) const;
+    inline bool not_inner_boundary(const Size p) const;
 };
 
 #include "geometry.tpp"

src/model/geometry/geometry.tpp

@@ -234,6 +234,16 @@ inline bool Geometry ::not_on_boundary(const Size p) const {
     return (boundary.point2boundary[p] == parameters->npoints());
 }
 
+///  Check whether a point is not an inner boundary point
+///    @param[in] p : point index
+///    @returns true if p is not on the boundary
+inline bool Geometry ::not_inner_boundary(const Size p) const {
+    if (not_on_boundary(p)) {
+        return true;
+    }
+    return (boundary.is_inner_boundary[boundary.point2boundary[p]] == 0);
+}
+
 ///  Getter for the number of the next cell on ray and its distance along ray in
 ///  the general case without any further assumptions
 ///    @param[in]      origin : position from which the ray originates

src/solver/solver.tpp

@@ -865,9 +865,11 @@ accel inline Size Solver ::trace_ray_imaging(const Model& model, const Vector3D&
 
         set_data(model, crt, nxt, shift_crt, shift_nxt, dZ, increment, id1, id2);
 
-        while (model.geometry.not_on_boundary(nxt)
-               || (model.geometry.points.position[nxt] - model.geometry.points.center).dot(raydir)
-                      < 0) {
+        while (
+            model.geometry.not_on_boundary(nxt)
+            || (model.geometry.not_inner_boundary(nxt)
+                && (model.geometry.points.position[nxt] - model.geometry.points.center).dot(raydir)
+                       < 0)) {
             crt       = nxt;
             shift_crt = shift_nxt;
 
@@ -925,9 +927,11 @@ accel inline Size Solver ::trace_ray_imaging_for_line(const Model& model, const
 
         set_data_for_line(model, l, crt, nxt, shift_crt, shift_nxt, dZ, increment, id1, id2);
 
-        while (model.geometry.not_on_boundary(nxt)
-               || (model.geometry.points.position[nxt] - model.geometry.points.center).dot(raydir)
-                      < 0) {
+        while (
+            model.geometry.not_on_boundary(nxt)
+            || (model.geometry.not_inner_boundary(nxt)
+                && (model.geometry.points.position[nxt] - model.geometry.points.center).dot(raydir)
+                       < 0)) {
             crt       = nxt;
             shift_crt = shift_nxt;
 
@@ -973,9 +977,11 @@ accel inline Size Solver ::get_ray_length_new_imager(
     if (model.geometry.valid_point(nxt)) {
         l += interp_helper.get_n_interp(model, crt, nxt);
 
-        while (model.geometry.not_on_boundary(nxt)
-               || (model.geometry.points.position[nxt] - model.geometry.points.center).dot(raydir)
-                      < 0) {
+        while (
+            model.geometry.not_on_boundary(nxt)
+            || (model.geometry.not_inner_boundary(nxt)
+                && (model.geometry.points.position[nxt] - model.geometry.points.center).dot(raydir)
+                       < 0)) {
             crt = nxt;
             nxt = model.geometry.get_next<Imagetracer>(origin, raydir, nxt, Z, dZ);
             if (!model.geometry.valid_point(nxt)) {