From d0f498df90bb90f1ed84fab676a6ab092dff2279 Mon Sep 17 00:00:00 2001 From: Thomas Ceulemans Date: Tue, 3 Dec 2024 16:57:36 +0100 Subject: [PATCH] Implemented inner boundary for imaging --- magritte/setup.py | 33 +++++++++++++++++++++--- src/bindings/pybindings.cpp | 2 +- src/model/geometry/boundary/boundary.cpp | 17 +++++++++++- src/model/geometry/boundary/boundary.hpp | 12 +++++---- src/model/geometry/geometry.hpp | 1 + src/model/geometry/geometry.tpp | 10 +++++++ src/solver/solver.tpp | 24 ++++++++++------- 7 files changed, 80 insertions(+), 19 deletions(-) diff --git a/magritte/setup.py b/magritte/setup.py index a6208352..441ca958 100644 --- a/magritte/setup.py +++ b/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 diff --git a/src/bindings/pybindings.cpp b/src/bindings/pybindings.cpp index 240389b4..ef25e2ce 100644 --- a/src/bindings/pybindings.cpp +++ b/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 diff --git a/src/model/geometry/boundary/boundary.cpp b/src/model/geometry/boundary/boundary.cpp index b8a7a62d..d533dc07 100644 --- a/src/model/geometry/boundary/boundary.cpp +++ b/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]; } diff --git a/src/model/geometry/boundary/boundary.hpp b/src/model/geometry/boundary/boundary.hpp index 83c252ad..29e9abc1 100644 --- a/src/model/geometry/boundary/boundary.hpp +++ b/src/model/geometry/boundary/boundary.hpp @@ -9,17 +9,19 @@ enum BoundaryCondition { Zero, Thermal, CMB }; struct Boundary { std::shared_ptr parameters; - Vector boundary2point; // maps boundary index ∈[0, nboundary-1] to point - // index ∈[0, npoints-1] - Vector point2boundary; // maps point index to boundary index. Maps non-boundary - // points to parameters->npoints() (invalid value) + Vector boundary2point; // maps boundary index ∈[0, nboundary-1] to point + // index ∈[0, npoints-1] + Vector point2boundary; // maps point index to boundary index. Maps non-boundary + // points to parameters->npoints() (invalid value) + Vector is_inner_boundary; // 1 if boundary is an inner boundary, 0 otherwise Vector boundary_condition; Vector boundary_temperature; Boundary(std::shared_ptr 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); diff --git a/src/model/geometry/geometry.hpp b/src/model/geometry/geometry.hpp index 1b12840c..e4305691 100644 --- a/src/model/geometry/geometry.hpp +++ b/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" diff --git a/src/model/geometry/geometry.tpp b/src/model/geometry/geometry.tpp index fdff3578..e711c693 100644 --- a/src/model/geometry/geometry.tpp +++ b/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 diff --git a/src/solver/solver.tpp b/src/solver/solver.tpp index f4a60b0f..d52d72f3 100644 --- a/src/solver/solver.tpp +++ b/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(origin, raydir, nxt, Z, dZ); if (!model.geometry.valid_point(nxt)) {