From 345acf900b1d036f7133933cba1db3c5b50f0eea Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 13 Dec 2024 11:05:00 +0700 Subject: [PATCH 01/14] add moving surface to surface.py --- openmc/surface.py | 233 +++++++++++++++++++++++++++++++++++----------- 1 file changed, 181 insertions(+), 52 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index c95025f949b..b1335799547 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -148,7 +148,13 @@ class Surface(IDManagerMixin, ABC): Name of the surface type : str Type of the surface - + moving : bool + Whether or not the surface is moving + moving_velocities : iterable of iterable of float + 3D velocities in which the surface moves + moving_durations : iterable of float + Durations during which the surface moves with the associated + velocities. """ next_id = 1 @@ -167,6 +173,10 @@ def __init__(self, surface_id=None, boundary_type='transmission', # Value - coefficient value self._coefficients = {} + self._moving = False + self._moving_velocities = np.array([[]]) + self._moving_durations = np.array([]) + def __neg__(self): return Halfspace(self, '-') @@ -192,6 +202,24 @@ def __repr__(self): string += coefficients + # TODO: print moving durations and velocities + if self._moving: + ''' + string += '{0: <20}{1}{2}\n'.format( + 'moving velocities', '=\t', np.array_repr(self._moving_velocities).replace('\n', ' ')) + string += '{0: <20}{1}{2}\n'.format('moving durations', + '=\t', self._moving_durations) + ''' + string += "\tMove parameters\n" + N_move = len(self._moving_durations) + string += "T Vx Vy Vz\n" + for i in range(N_move): + string += '{0: <8} {1: <8} {2: <8} {3: <8}\n'.format( + self._moving_durations[i], + self._moving_velocities[i,0], + self._moving_velocities[i,1], + self._moving_velocities[i,2]) + return string @property @@ -234,6 +262,18 @@ def albedo(self, albedo): def coefficients(self): return self._coefficients + @property + def moving(self): + return self._moving + + @property + def moving_velicities(self): + return self._moving_velocities + + @property + def moving_durations(self): + return self._moving_durations + def bounding_box(self, side): """Determine an axis-aligned bounding box. @@ -327,6 +367,69 @@ def is_equal(self, other): return np.allclose(coeffs1, coeffs2, rtol=0., atol=self._atol) + def move(self, velocities, durations): + """Set the surface continuous movement based on the given velocities + and durations. + + Parameters + ========== + velocities : iterable of iterable of float + 3D velocities in which the surface moves + durations : iterable of float + Durations during which the surface moves with the associated + velocities. + """ + self._moving = True + self._moving_velocities = np.array(velocities) + self._moving_durations = np.array(durations) + + # Add the final static condition + self._moving_velocities = np.append(self._moving_velocities, np.zeros((1,3)), axis=0) + self._moving_durations = np.append(self._moving_durations, np.inf) + + # Set moving time grid and translations for evaluation convenience + N_move = len(self._moving_durations) + self._moving_time_grid = np.zeros(N_move + 1) + self._moving_translations = np.zeros((N_move + 1, 3)) + for n in range(N_move - 1): + duration = self._moving_durations[n] + velocity = self._moving_velocities[n] + + t_start = self._moving_time_grid[n] + self._moving_time_grid[n + 1] = t_start + duration + + trans_start = self._moving_translations[n] + self._moving_translations[n + 1] = trans_start + velocity * duration + # Manual assignment to avoid Numpy's warning + self._moving_time_grid[N_move] = np.inf + self._moving_translations[N_move] = self._moving_translations[N_move - 1] + + + def _adjust_point(self, point, time): + """Adjust point based on the surface movement""" + if not self._moving: + return point + + # Get moving index + N_move = len(self._moving_durations) + time_grid = self._moving_time_grid + idx = np.searchsorted(time_grid, time, side='right') - 1 + + # Get moving translation, velocity, and starting time + trans_0 = self._moving_translations[idx] + time_0 = self._moving_time_grid[idx] + V = self._moving_velocities[idx] + + # Translate the particle + t_local = time - time_0 + x = point[0] - (trans_0[0] + V[0] * t_local) + y = point[1] - (trans_0[1] + V[1] * t_local) + z = point[2] - (trans_0[2] + V[2] * t_local) + + print(time_grid) + print(time, idx) + return (x, y, z) + @abstractmethod def _get_base_coeffs(self): """Return polynomial coefficients representing the implicit surface @@ -335,19 +438,22 @@ def _get_base_coeffs(self): """ @abstractmethod - def evaluate(self, point): - """Evaluate the surface equation at a given point. + def evaluate(self, point, time=0.0): + """Evaluate the surface equation at a given point and time. Parameters ---------- point : 3-tuple of float The Cartesian coordinates, :math:`(x',y',z')`, at which the surface equation should be evaluated. + time : float, optional + The time :math:`t'` at which the surface equation should be evaluated Returns ------- float Evaluation of the surface polynomial at point :math:`(x',y',z')` + and time :math:`t'` """ @@ -432,6 +538,12 @@ def to_xml_element(self): element.set("coeffs", ' '.join([str(self._coefficients.setdefault(key, 0.0)) for key in self._coeff_keys])) + if self._moving: + element.set("moving_velocities", ' '.join( + [str(value) for value in self._moving_velocities[:-1]])) + element.set("moving_durations", ' '.join( + [str(value) for value in self._moving_durations[:-1]])) + return element @staticmethod @@ -574,14 +686,16 @@ def bounding_box(self, side): return BoundingBox(ll, ur) - def evaluate(self, point): - """Evaluate the surface equation at a given point. + def evaluate(self, point, time=0.0): + """Evaluate the surface equation at a given point and time. Parameters ---------- point : 3-tuple of float The Cartesian coordinates, :math:`(x',y',z')`, at which the surface equation should be evaluated. + time : float, optional + The time :math:`t'` at which the surface equation should be evaluated Returns ------- @@ -590,7 +704,7 @@ def evaluate(self, point): """ - x, y, z = point + x, y, z = self._adjust_point(point, time) a, b, c, d = self._get_base_coeffs() return a*x + b*y + c*z - d @@ -872,8 +986,9 @@ def __init__(self, x0=0., *args, **kwargs): c = SurfaceCoefficient(0.) d = x0 - def evaluate(self, point): - return point[0] - self.x0 + def evaluate(self, point, time=0.0): + x, y, z = self._adjust_point(point, time) + return x - self.x0 class YPlane(PlaneMixin, Surface): @@ -937,8 +1052,9 @@ def __init__(self, y0=0., *args, **kwargs): c = SurfaceCoefficient(0.) d = y0 - def evaluate(self, point): - return point[1] - self.y0 + def evaluate(self, point, time=0.0): + x, y, z = self._adjust_point(point, time) + return y - self.y0 class ZPlane(PlaneMixin, Surface): @@ -1002,8 +1118,9 @@ def __init__(self, z0=0., *args, **kwargs): c = SurfaceCoefficient(1.) d = z0 - def evaluate(self, point): - return point[2] - self.z0 + def evaluate(self, point, time=0.0): + x, y, z = self._adjust_point(point, time) + return z - self.z0 class QuadricMixin: @@ -1060,14 +1177,16 @@ def eigh(self, coeffs=None): """ return np.linalg.eigh(self.get_Abc(coeffs=coeffs)[0]) - def evaluate(self, point): - """Evaluate the surface equation at a given point. + def evaluate(self, point, time=0.0): + """Evaluate the surface equation at a given point and time. Parameters ---------- point : 3-tuple of float The Cartesian coordinates, :math:`(x',y',z')`, in [cm] at which the surface equation should be evaluated. + time : float, optional + The time :math:`t'` at which the surface equation should be evaluated Returns ------- @@ -1076,7 +1195,7 @@ def evaluate(self, point): Jz' + K = 0` """ - x = np.asarray(point) + x = np.asarray(self._adjust_point(point, time)) A, b, c = self.get_Abc() return x.T @ A @ x + b.T @ x + c @@ -1452,9 +1571,10 @@ def bounding_box(self, side): elif side == '+': return BoundingBox.infinite() - def evaluate(self, point): - y = point[1] - self.y0 - z = point[2] - self.z0 + def evaluate(self, point, time=0.0): + x, y, z = self._adjust_point(point, time) + y -= self.y0 + z -= self.z0 return y*y + z*z - self.r**2 @@ -1550,9 +1670,10 @@ def bounding_box(self, side): elif side == '+': return BoundingBox.infinite() - def evaluate(self, point): - x = point[0] - self.x0 - z = point[2] - self.z0 + def evaluate(self, point, time=0.0): + x, y, z = self._adjust_point(point, time) + x -= self.x0 + z -= self.z0 return x*x + z*z - self.r**2 @@ -1648,9 +1769,10 @@ def bounding_box(self, side): elif side == '+': return BoundingBox.infinite() - def evaluate(self, point): - x = point[0] - self.x0 - y = point[1] - self.y0 + def evaluate(self, point, time=0.0): + x, y, z = self._adjust_point(point, time) + x -= self.x0 + y -= self.y0 return x*x + y*y - self.r**2 @@ -1745,10 +1867,11 @@ def bounding_box(self, side): elif side == '+': return BoundingBox.infinite() - def evaluate(self, point): - x = point[0] - self.x0 - y = point[1] - self.y0 - z = point[2] - self.z0 + def evaluate(self, point, time=0.0): + x, y, z = self._adjust_point(point, time) + x -= self.x0 + y -= self.y0 + z -= self.z0 return x*x + y*y + z*z - self.r**2 @@ -2004,10 +2127,11 @@ def _get_base_coeffs(self): return (a, b, c, d, e, f, g, h, j, k) - def evaluate(self, point): - x = point[0] - self.x0 - y = point[1] - self.y0 - z = point[2] - self.z0 + def evaluate(self, point, time=0.0): + x, y, z = self._adjust_point(point, time) + x -= self.x0 + y -= self.y0 + z -= self.z0 return y*y + z*z - self.r2*x*x @@ -2105,10 +2229,11 @@ def _get_base_coeffs(self): return (a, b, c, d, e, f, g, h, j, k) - def evaluate(self, point): - x = point[0] - self.x0 - y = point[1] - self.y0 - z = point[2] - self.z0 + def evaluate(self, point, time=0.0): + x, y, z = self._adjust_point(point, time) + x -= self.x0 + y -= self.y0 + z -= self.z0 return x*x + z*z - self.r2*y*y @@ -2206,10 +2331,11 @@ def _get_base_coeffs(self): return (a, b, c, d, e, f, g, h, j, k) - def evaluate(self, point): - x = point[0] - self.x0 - y = point[1] - self.y0 - z = point[2] - self.z0 + def evaluate(self, point, time=0.0): + x, y, z = self._adjust_point(point, time) + x -= self.x0 + y -= self.y0 + z -= self.z0 return x*x + y*y - self.r2*z*z @@ -2408,10 +2534,11 @@ class XTorus(TorusMixin, Surface): """ _type = 'x-torus' - def evaluate(self, point): - x = point[0] - self.x0 - y = point[1] - self.y0 - z = point[2] - self.z0 + def evaluate(self, point, time=0.0): + x, y, z = self._adjust_point(point, time) + x -= self.x0 + y -= self.y0 + z -= self.z0 a = self.a b = self.b c = self.c @@ -2483,10 +2610,11 @@ class YTorus(TorusMixin, Surface): """ _type = 'y-torus' - def evaluate(self, point): - x = point[0] - self.x0 - y = point[1] - self.y0 - z = point[2] - self.z0 + def evaluate(self, point, time=0.0): + x, y, z = self._adjust_point(point, time) + x -= self.x0 + y -= self.y0 + z -= self.z0 a = self.a b = self.b c = self.c @@ -2558,10 +2686,11 @@ class ZTorus(TorusMixin, Surface): _type = 'z-torus' - def evaluate(self, point): - x = point[0] - self.x0 - y = point[1] - self.y0 - z = point[2] - self.z0 + def evaluate(self, point, time=0.0): + x, y, z = self._adjust_point(point, time) + x -= self.x0 + y -= self.y0 + z -= self.z0 a = self.a b = self.b c = self.c From 14030f289499cfca6b2d4915d77f542ed79f7117 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Sat, 14 Dec 2024 16:01:22 -0800 Subject: [PATCH 02/14] add moving surface input parser in surface.cpp --- include/openmc/surface.h | 6 +++++ src/surface.cpp | 56 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 62 insertions(+) diff --git a/include/openmc/surface.h b/include/openmc/surface.h index af235301c14..93f23f46f05 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -41,6 +41,12 @@ class Surface { GeometryType geom_type_; //!< Geometry type indicator (CSG or DAGMC) bool surf_source_ {false}; //!< Activate source banking for the surface? + //!< Moving surface parameters + bool moving_ {false}; + vector moving_time_grid_; // Time grid points [0.0, ..., INFTY] + vector moving_translations_; // Translations at the time grid points + vector moving_velocities_; // Velocities within time bins + explicit Surface(pugi::xml_node surf_node); Surface(); diff --git a/src/surface.cpp b/src/surface.cpp index 50ef2a12830..4b01072b21f 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -191,9 +191,65 @@ CSGSurface::CSGSurface() : Surface {} { geom_type_ = GeometryType::CSG; }; + CSGSurface::CSGSurface(pugi::xml_node surf_node) : Surface {surf_node} { geom_type_ = GeometryType::CSG; + + // Not moving? + if (!(check_for_node(surf_node, "moving_velocities") || + check_for_node(surf_node, "moving_durations"))){ + moving_ = false; + return; + } + + // Now, set the surface moving parameters + moving_ = true; + + // Moving durations + auto durations = get_node_array(surf_node, "moving_durations"); + const int N_move = durations.size() + 1; + + // Moving time grids + moving_time_grid_.resize(N_move + 1); + moving_time_grid_[0] = 0.0; + for (int n = 0; n < N_move - 1; n++) { + moving_time_grid_[n + 1] = moving_time_grid_[n] + durations[n]; + } + moving_time_grid_[N_move] = INFTY; + + // Moving velocities + moving_velocities_.resize(N_move); + std::string velocities_spec = get_node_value(surf_node, "moving_velocities"); + // Parse + std::vector numbers; + for (int i = 0; i < velocities_spec.size();) { + if (velocities_spec[i] == '-' || std::isdigit(velocities_spec[i])) { + int j = i + 1; + while (j < velocities_spec.size() && std::isdigit(velocities_spec[j])) { + j++; + } + numbers.push_back(std::stod(velocities_spec.substr(i, j - i))); + i = j; + } + i++; + } + // Assign to velocities + for (int n = 0; n < N_move - 1; n++) { + int idx = 3 * n; + moving_velocities_[n][0] = numbers[idx]; + moving_velocities_[n][1] = numbers[idx + 1]; + moving_velocities_[n][2] = numbers[idx + 2]; + } + moving_velocities_[N_move - 1] *= 0.0; + + // Moving translations + moving_translations_.resize(N_move + 1); + moving_translations_[0] *= 0.0; + for (int n = 0; n < N_move - 1; n++) { + moving_translations_[n + 1] = moving_translations_[n] + moving_velocities_[n] * durations[n]; + } + moving_translations_[N_move] = moving_translations_[N_move - 1]; }; //============================================================================== From c0519737f7072b1f78e48ac9209355eca87b17c5 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Sat, 14 Dec 2024 22:34:59 -0800 Subject: [PATCH 03/14] only CSG surfaces can move --- include/openmc/surface.h | 12 ++++++------ src/surface.cpp | 11 ++++++----- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 93f23f46f05..b02d9cf4cdf 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -41,12 +41,6 @@ class Surface { GeometryType geom_type_; //!< Geometry type indicator (CSG or DAGMC) bool surf_source_ {false}; //!< Activate source banking for the surface? - //!< Moving surface parameters - bool moving_ {false}; - vector moving_time_grid_; // Time grid points [0.0, ..., INFTY] - vector moving_translations_; // Translations at the time grid points - vector moving_velocities_; // Velocities within time bins - explicit Surface(pugi::xml_node surf_node); Surface(); @@ -103,6 +97,12 @@ class Surface { class CSGSurface : public Surface { public: + //!< Moving surface parameters + bool moving_ {false}; + vector moving_time_grid_; // Time grid points [0.0, ..., INFTY] + vector moving_translations_; // Translations at the time grid points + vector moving_velocities_; // Velocities within time bins + explicit CSGSurface(pugi::xml_node surf_node); CSGSurface(); }; diff --git a/src/surface.cpp b/src/surface.cpp index 4b01072b21f..2f4b5d8cfba 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -197,12 +197,12 @@ CSGSurface::CSGSurface(pugi::xml_node surf_node) : Surface {surf_node} geom_type_ = GeometryType::CSG; // Not moving? - if (!(check_for_node(surf_node, "moving_velocities") || - check_for_node(surf_node, "moving_durations"))){ + if (!(check_for_node(surf_node, "moving_velocities") || + check_for_node(surf_node, "moving_durations"))) { moving_ = false; return; } - + // Now, set the surface moving parameters moving_ = true; @@ -217,7 +217,7 @@ CSGSurface::CSGSurface(pugi::xml_node surf_node) : Surface {surf_node} moving_time_grid_[n + 1] = moving_time_grid_[n] + durations[n]; } moving_time_grid_[N_move] = INFTY; - + // Moving velocities moving_velocities_.resize(N_move); std::string velocities_spec = get_node_value(surf_node, "moving_velocities"); @@ -247,7 +247,8 @@ CSGSurface::CSGSurface(pugi::xml_node surf_node) : Surface {surf_node} moving_translations_.resize(N_move + 1); moving_translations_[0] *= 0.0; for (int n = 0; n < N_move - 1; n++) { - moving_translations_[n + 1] = moving_translations_[n] + moving_velocities_[n] * durations[n]; + moving_translations_[n + 1] = + moving_translations_[n] + moving_velocities_[n] * durations[n]; } moving_translations_[N_move] = moving_translations_[N_move - 1]; }; From aef0321f333e86b33b77fc66ab1fa9a2290d9636 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Sun, 15 Dec 2024 01:56:14 -0800 Subject: [PATCH 04/14] add moving surface evaluation. add time to GeometryState --- include/openmc/cell.h | 14 ++++--- include/openmc/dagmc.h | 2 +- include/openmc/particle_data.h | 18 ++++---- include/openmc/surface.h | 76 +++++++++++++++++++++++++--------- include/openmc/universe.h | 2 +- openmc/surface.py | 1 - src/boundary_condition.cpp | 8 ++-- src/cell.cpp | 14 +++---- src/dagmc.cpp | 2 +- src/geometry.cpp | 7 ++-- src/surface.cpp | 59 ++++++++++++++++++-------- src/universe.cpp | 9 ++-- 12 files changed, 139 insertions(+), 73 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 70843140bad..30a0bb5d89c 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -77,10 +77,11 @@ class Region { //! \param r The 3D Cartesian coordinate to check. //! \param u A direction used to "break ties" the coordinates are very //! close to a surface. + //! \param t The time coordinate to check. //! \param on_surface The signed index of a surface that the coordinate is //! known to be on. This index takes precedence over surface sense //! calculations. - bool contains(Position r, Direction u, int32_t on_surface) const; + bool contains(Position r, Direction u, double t, int32_t on_surface) const; //! Find the oncoming boundary of this cell. std::pair distance( @@ -110,14 +111,14 @@ class Region { //! Determine if a particle is inside the cell for a simple cell (only //! intersection operators) - bool contains_simple(Position r, Direction u, int32_t on_surface) const; + bool contains_simple(Position r, Direction u, double t, int32_t on_surface) const; //! Determine if a particle is inside the cell for a complex cell. //! //! Uses the comobination of half-spaces and binary operators to determine //! if short circuiting can be used. Short cicuiting uses the relative and //! absolute depth of parenthases in the expression. - bool contains_complex(Position r, Direction u, int32_t on_surface) const; + bool contains_complex(Position r, Direction u, double t, int32_t on_surface) const; //! BoundingBox if the paritcle is in a simple cell. BoundingBox bounding_box_simple() const; @@ -177,10 +178,11 @@ class Cell { //! \param r The 3D Cartesian coordinate to check. //! \param u A direction used to "break ties" the coordinates are very //! close to a surface. + //! \param t The time coordinate to check. //! \param on_surface The signed index of a surface that the coordinate is //! known to be on. This index takes precedence over surface sense //! calculations. - virtual bool contains(Position r, Direction u, int32_t on_surface) const = 0; + virtual bool contains(Position r, Direction u, double t, int32_t on_surface) const = 0; //! Find the oncoming boundary of this cell. virtual std::pair distance( @@ -376,9 +378,9 @@ class CSGCell : public Cell { return region_.distance(r, u, on_surface); } - bool contains(Position r, Direction u, int32_t on_surface) const override + bool contains(Position r, Direction u, double t, int32_t on_surface) const override { - return region_.contains(r, u, on_surface); + return region_.contains(r, u, t, on_surface); } BoundingBox bounding_box() const override diff --git a/include/openmc/dagmc.h b/include/openmc/dagmc.h index 2facf4fc05e..14a62e1b4f3 100644 --- a/include/openmc/dagmc.h +++ b/include/openmc/dagmc.h @@ -40,7 +40,7 @@ class DAGSurface : public Surface { moab::EntityHandle mesh_handle() const; - double evaluate(Position r) const override; + double evaluate(Position r, double t) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; Direction reflect(Position r, Direction u, GeometryState* p) const override; diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 164148cce10..53b8268a94a 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -296,6 +296,12 @@ class GeometryState { Direction& u_local() { return coord_[n_coord_ - 1].u; } const Direction& u_local() const { return coord_[n_coord_ - 1].u; } + // Accessors for time (units are seconds). + double& time() { return time_; } + const double& time() const { return time_; } + double& time_last() { return time_last_; } + const double& time_last() const { return time_last_; } + // Surface that the particle is on int& surface() { return surface_; } const int& surface() const { return surface_; } @@ -337,6 +343,9 @@ class GeometryState { Position r_last_; //!< previous coordinates Direction u_last_; //!< previous direction coordinates + double time_; //!< time + double time_last_; //!< previous time + int surface_ {0}; //!< index for surface particle is on BoundaryInfo boundary_; //!< Info about the next intersection @@ -407,8 +416,6 @@ class ParticleData : public GeometryState { double wgt_ {1.0}; double mu_; - double time_ {0.0}; - double time_last_ {0.0}; double wgt_last_ {1.0}; bool fission_ {false}; @@ -515,13 +522,6 @@ class ParticleData : public GeometryState { double& mu() { return mu_; } const double& mu() const { return mu_; } - // Tracks the time of a particle as it traverses the problem. - // Units are seconds. - double& time() { return time_; } - const double& time() const { return time_; } - double& time_last() { return time_last_; } - const double& time_last() const { return time_last_; } - // What event took place, described in greater detail below TallyEvent& event() { return event_; } const TallyEvent& event() const { return event_; } diff --git a/include/openmc/surface.h b/include/openmc/surface.h index b02d9cf4cdf..dfb3cab4aa8 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -50,9 +50,10 @@ class Surface { //! \param r The 3D Cartesian coordinate of a point. //! \param u A direction used to "break ties" and pick a sense when the //! point is very close to the surface. + //! \param t The time for the evaluation //! \return true if the point is on the "positive" side of the surface and //! false otherwise. - bool sense(Position r, Direction u) const; + bool sense(Position r, Direction u, double t) const; //! Determine the direction of a ray reflected from the surface. //! \param[in] r The point at which the ray is incident. @@ -67,10 +68,12 @@ class Surface { //! Evaluate the equation describing the surface. //! - //! Surfaces can be described by some function f(x, y, z) = 0. This member - //! function evaluates that mathematical function. + //! Surfaces can be described by some function f(x, y, z) = 0. This member + //! function evaluates that mathematical function. The time variable t is + //! needed for evaluation of moving surfaces. //! \param r A 3D Cartesian coordinate. - virtual double evaluate(Position r) const = 0; + //! \param t The time for the evaluation. + virtual double evaluate(Position r, double t) const = 0; //! Compute the distance between a point and the surface along a ray. //! \param r A 3D Cartesian coordinate. @@ -105,6 +108,11 @@ class CSGSurface : public Surface { explicit CSGSurface(pugi::xml_node surf_node); CSGSurface(); + double evaluate(Position r, double t) const override; + +protected: + //! Static CSG surface evaluation. + virtual double _evaluate(Position r) const = 0; }; //============================================================================== @@ -116,13 +124,15 @@ class CSGSurface : public Surface { class SurfaceXPlane : public CSGSurface { public: explicit SurfaceXPlane(pugi::xml_node surf_node); - double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; BoundingBox bounding_box(bool pos_side) const override; double x0_; + +protected: + double _evaluate(Position r) const override; }; //============================================================================== @@ -134,13 +144,15 @@ class SurfaceXPlane : public CSGSurface { class SurfaceYPlane : public CSGSurface { public: explicit SurfaceYPlane(pugi::xml_node surf_node); - double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; BoundingBox bounding_box(bool pos_side) const override; double y0_; + +protected: + double _evaluate(Position r) const override; }; //============================================================================== @@ -152,13 +164,15 @@ class SurfaceYPlane : public CSGSurface { class SurfaceZPlane : public CSGSurface { public: explicit SurfaceZPlane(pugi::xml_node surf_node); - double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; BoundingBox bounding_box(bool pos_side) const override; double z0_; + +protected: + double _evaluate(Position r) const override; }; //============================================================================== @@ -170,12 +184,14 @@ class SurfaceZPlane : public CSGSurface { class SurfacePlane : public CSGSurface { public: explicit SurfacePlane(pugi::xml_node surf_node); - double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; double A_, B_, C_, D_; + +protected: + double _evaluate(Position r) const override; }; //============================================================================== @@ -188,13 +204,15 @@ class SurfacePlane : public CSGSurface { class SurfaceXCylinder : public CSGSurface { public: explicit SurfaceXCylinder(pugi::xml_node surf_node); - double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; BoundingBox bounding_box(bool pos_side) const override; double y0_, z0_, radius_; + +protected: + double _evaluate(Position r) const override; }; //============================================================================== @@ -207,13 +225,15 @@ class SurfaceXCylinder : public CSGSurface { class SurfaceYCylinder : public CSGSurface { public: explicit SurfaceYCylinder(pugi::xml_node surf_node); - double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; BoundingBox bounding_box(bool pos_side) const override; double x0_, z0_, radius_; + +protected: + double _evaluate(Position r) const override; }; //============================================================================== @@ -226,13 +246,15 @@ class SurfaceYCylinder : public CSGSurface { class SurfaceZCylinder : public CSGSurface { public: explicit SurfaceZCylinder(pugi::xml_node surf_node); - double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; BoundingBox bounding_box(bool pos_side) const override; double x0_, y0_, radius_; + +protected: + double _evaluate(Position r) const override; }; //============================================================================== @@ -245,13 +267,15 @@ class SurfaceZCylinder : public CSGSurface { class SurfaceSphere : public CSGSurface { public: explicit SurfaceSphere(pugi::xml_node surf_node); - double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; BoundingBox bounding_box(bool pos_side) const override; double x0_, y0_, z0_, radius_; + +protected: + double _evaluate(Position r) const override; }; //============================================================================== @@ -264,12 +288,14 @@ class SurfaceSphere : public CSGSurface { class SurfaceXCone : public CSGSurface { public: explicit SurfaceXCone(pugi::xml_node surf_node); - double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; double x0_, y0_, z0_, radius_sq_; + +protected: + double _evaluate(Position r) const override; }; //============================================================================== @@ -282,12 +308,14 @@ class SurfaceXCone : public CSGSurface { class SurfaceYCone : public CSGSurface { public: explicit SurfaceYCone(pugi::xml_node surf_node); - double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; double x0_, y0_, z0_, radius_sq_; + +protected: + double _evaluate(Position r) const override; }; //============================================================================== @@ -300,12 +328,14 @@ class SurfaceYCone : public CSGSurface { class SurfaceZCone : public CSGSurface { public: explicit SurfaceZCone(pugi::xml_node surf_node); - double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; double x0_, y0_, z0_, radius_sq_; + +protected: + double _evaluate(Position r) const override; }; //============================================================================== @@ -318,13 +348,15 @@ class SurfaceZCone : public CSGSurface { class SurfaceQuadric : public CSGSurface { public: explicit SurfaceQuadric(pugi::xml_node surf_node); - double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; // Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0 double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_; + +protected: + double _evaluate(Position r) const override; }; //============================================================================== @@ -336,12 +368,14 @@ class SurfaceQuadric : public CSGSurface { class SurfaceXTorus : public CSGSurface { public: explicit SurfaceXTorus(pugi::xml_node surf_node); - double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; double x0_, y0_, z0_, A_, B_, C_; + +protected: + double _evaluate(Position r) const override; }; //============================================================================== @@ -353,12 +387,14 @@ class SurfaceXTorus : public CSGSurface { class SurfaceYTorus : public CSGSurface { public: explicit SurfaceYTorus(pugi::xml_node surf_node); - double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; double x0_, y0_, z0_, A_, B_, C_; + +protected: + double _evaluate(Position r) const override; }; //============================================================================== @@ -370,12 +406,14 @@ class SurfaceYTorus : public CSGSurface { class SurfaceZTorus : public CSGSurface { public: explicit SurfaceZTorus(pugi::xml_node surf_node); - double evaluate(Position r) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; double x0_, y0_, z0_, A_, B_, C_; + +protected: + double _evaluate(Position r) const override; }; //============================================================================== diff --git a/include/openmc/universe.h b/include/openmc/universe.h index 9fea06bccba..218b6039d2c 100644 --- a/include/openmc/universe.h +++ b/include/openmc/universe.h @@ -60,7 +60,7 @@ class UniversePartitioner { explicit UniversePartitioner(const Universe& univ); //! Return the list of cells that could contain the given coordinates. - const vector& get_cells(Position r, Direction u) const; + const vector& get_cells(Position r, Direction u, double t) const; private: //! A sorted vector of indices to surfaces that partition the universe diff --git a/openmc/surface.py b/openmc/surface.py index b1335799547..61f87b29692 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -411,7 +411,6 @@ def _adjust_point(self, point, time): return point # Get moving index - N_move = len(self._moving_durations) time_grid = self._moving_time_grid idx = np.searchsorted(time_grid, time, side='right') - 1 diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index b58054dce8c..74d0b7604b3 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -101,7 +101,7 @@ TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf) Position origin {0, 0, 0}; Direction u = surf1.normal(origin); double d1; - double e1 = surf1.evaluate(origin); + double e1 = surf1.evaluate(origin, 0.0); // Time = 0.0 if (e1 > FP_COINCIDENT) { d1 = -surf1.distance(origin, -u, false); } else if (e1 < -FP_COINCIDENT) { @@ -112,7 +112,7 @@ TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf) // Compute the distance from the second surface to the origin. double d2; - double e2 = surf2.evaluate(origin); + double e2 = surf2.evaluate(origin, 0.0); // Time = 0.0 if (e2 > FP_COINCIDENT) { d2 = -surf2.distance(origin, -u, false); } else if (e2 < -FP_COINCIDENT) { @@ -217,14 +217,14 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf) } // Make sure both surfaces intersect the origin - if (std::abs(surf1.evaluate({0, 0, 0})) > FP_COINCIDENT) { + if (std::abs(surf1.evaluate({0, 0, 0}, 0.0)) > FP_COINCIDENT) { // Time = 0.0 throw std::invalid_argument(fmt::format( "Rotational periodic BCs are only " "supported for rotations about the origin, but surface {} does not " "intersect the origin.", surf1.id_)); } - if (std::abs(surf2.evaluate({0, 0, 0})) > FP_COINCIDENT) { + if (std::abs(surf2.evaluate({0, 0, 0}, 0.0)) > FP_COINCIDENT) { // Time = 0.0 throw std::invalid_argument(fmt::format( "Rotational periodic BCs are only " "supported for rotations about the origin, but surface {} does not " diff --git a/src/cell.cpp b/src/cell.cpp index 88876678706..9b56a5b6123 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -812,18 +812,18 @@ std::pair Region::distance( //============================================================================== -bool Region::contains(Position r, Direction u, int32_t on_surface) const +bool Region::contains(Position r, Direction u, double t, int32_t on_surface) const { if (simple_) { - return contains_simple(r, u, on_surface); + return contains_simple(r, u, t, on_surface); } else { - return contains_complex(r, u, on_surface); + return contains_complex(r, u, t, on_surface); } } //============================================================================== -bool Region::contains_simple(Position r, Direction u, int32_t on_surface) const +bool Region::contains_simple(Position r, Direction u, double t, int32_t on_surface) const { for (int32_t token : expression_) { // Assume that no tokens are operators. Evaluate the sense of particle with @@ -835,7 +835,7 @@ bool Region::contains_simple(Position r, Direction u, int32_t on_surface) const return false; } else { // Note the off-by-one indexing - bool sense = model::surfaces[abs(token) - 1]->sense(r, u); + bool sense = model::surfaces[abs(token) - 1]->sense(r, u, t); if (sense != (token > 0)) { return false; } @@ -846,7 +846,7 @@ bool Region::contains_simple(Position r, Direction u, int32_t on_surface) const //============================================================================== -bool Region::contains_complex(Position r, Direction u, int32_t on_surface) const +bool Region::contains_complex(Position r, Direction u, double t, int32_t on_surface) const { bool in_cell = true; int total_depth = 0; @@ -865,7 +865,7 @@ bool Region::contains_complex(Position r, Direction u, int32_t on_surface) const in_cell = false; } else { // Note the off-by-one indexing - bool sense = model::surfaces[abs(token) - 1]->sense(r, u); + bool sense = model::surfaces[abs(token) - 1]->sense(r, u, t); in_cell = (sense == (token > 0)); } } else if ((token == OP_UNION && in_cell == true) || diff --git a/src/dagmc.cpp b/src/dagmc.cpp index b79676c3626..fa86631f758 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -727,7 +727,7 @@ moab::EntityHandle DAGSurface::mesh_handle() const return dagmc_ptr()->entity_by_index(2, dag_index()); } -double DAGSurface::evaluate(Position r) const +double DAGSurface::evaluate(Position r, double t) const { return 0.0; } diff --git a/src/geometry.cpp b/src/geometry.cpp index 445b19faac1..0cf954a6557 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -43,7 +43,7 @@ bool check_cell_overlap(GeometryState& p, bool error) // Loop through each cell on this level for (auto index_cell : univ.cells_) { Cell& c = *model::cells[index_cell]; - if (c.contains(p.coord(j).r, p.coord(j).u, p.surface())) { + if (c.contains(p.coord(j).r, p.coord(j).u, p.time(), p.surface())) { if (index_cell != p.coord(j).cell) { if (error) { fatal_error( @@ -110,7 +110,7 @@ bool find_cell_inner( for (auto it = neighbor_list->cbegin(); it != neighbor_list->cend(); ++it) { i_cell = *it; - // Make sure the search cell is in the same universe. + // Make sure the search cell is in the same `universe. int i_universe = p.lowest_coord().universe; if (model::cells[i_cell]->universe_ != i_universe) continue; @@ -118,8 +118,9 @@ bool find_cell_inner( // Check if this cell contains the particle. Position r {p.r_local()}; Direction u {p.u_local()}; + double t {p.time()}; auto surf = p.surface(); - if (model::cells[i_cell]->contains(r, u, surf)) { + if (model::cells[i_cell]->contains(r, u, t, surf)) { p.lowest_coord().cell = i_cell; found = true; break; diff --git a/src/surface.cpp b/src/surface.cpp index 2f4b5d8cfba..c64bbfeb2f5 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -16,6 +16,7 @@ #include "openmc/hdf5_interface.h" #include "openmc/math_functions.h" #include "openmc/random_lcg.h" +#include "openmc/search.h" #include "openmc/settings.h" #include "openmc/string_utils.h" #include "openmc/xml_interface.h" @@ -114,11 +115,11 @@ Surface::Surface(pugi::xml_node surf_node) } } -bool Surface::sense(Position r, Direction u) const +bool Surface::sense(Position r, Direction u, double t) const { // Evaluate the surface equation at the particle's coordinates to determine // which side the particle is on. - const double f = evaluate(r); + const double f = evaluate(r, t); // Check which side of surface the point is on. if (std::abs(f) < FP_COINCIDENT) { @@ -253,6 +254,30 @@ CSGSurface::CSGSurface(pugi::xml_node surf_node) : Surface {surf_node} moving_translations_[N_move] = moving_translations_[N_move - 1]; }; +double CSGSurface::evaluate(Position r, double t) const +{ + if (!moving_) { + return _evaluate(r); + } + // The surface moves + + // Get moving index + int idx = lower_bound_index( + moving_time_grid_.begin(), moving_time_grid_.end(), t); + + // Get moving translation, velocity, and starting time + Position translation = moving_translations_[idx]; + double time_0 = moving_time_grid_[idx]; + Position velocity = moving_velocities_[idx]; + + // Move the position relative to the surface movement + double t_local = t - time_0; + Position r_moved = r - (translation + velocity * t_local); + + // Evaluate the moved position + return _evaluate(r_moved); +} + //============================================================================== // Generic functions for x-, y-, and z-, planes. //============================================================================== @@ -280,7 +305,7 @@ SurfaceXPlane::SurfaceXPlane(pugi::xml_node surf_node) : CSGSurface(surf_node) read_coeffs(surf_node, id_, {&x0_}); } -double SurfaceXPlane::evaluate(Position r) const +double SurfaceXPlane::_evaluate(Position r) const { return r.x - x0_; } @@ -320,7 +345,7 @@ SurfaceYPlane::SurfaceYPlane(pugi::xml_node surf_node) : CSGSurface(surf_node) read_coeffs(surf_node, id_, {&y0_}); } -double SurfaceYPlane::evaluate(Position r) const +double SurfaceYPlane::_evaluate(Position r) const { return r.y - y0_; } @@ -360,7 +385,7 @@ SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node) : CSGSurface(surf_node) read_coeffs(surf_node, id_, {&z0_}); } -double SurfaceZPlane::evaluate(Position r) const +double SurfaceZPlane::_evaluate(Position r) const { return r.z - z0_; } @@ -400,7 +425,7 @@ SurfacePlane::SurfacePlane(pugi::xml_node surf_node) : CSGSurface(surf_node) read_coeffs(surf_node, id_, {&A_, &B_, &C_, &D_}); } -double SurfacePlane::evaluate(Position r) const +double SurfacePlane::_evaluate(Position r) const { return A_ * r.x + B_ * r.y + C_ * r.z - D_; } @@ -519,7 +544,7 @@ SurfaceXCylinder::SurfaceXCylinder(pugi::xml_node surf_node) read_coeffs(surf_node, id_, {&y0_, &z0_, &radius_}); } -double SurfaceXCylinder::evaluate(Position r) const +double SurfaceXCylinder::_evaluate(Position r) const { return axis_aligned_cylinder_evaluate<1, 2>(r, y0_, z0_, radius_); } @@ -562,7 +587,7 @@ SurfaceYCylinder::SurfaceYCylinder(pugi::xml_node surf_node) read_coeffs(surf_node, id_, {&x0_, &z0_, &radius_}); } -double SurfaceYCylinder::evaluate(Position r) const +double SurfaceYCylinder::_evaluate(Position r) const { return axis_aligned_cylinder_evaluate<0, 2>(r, x0_, z0_, radius_); } @@ -606,7 +631,7 @@ SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node) read_coeffs(surf_node, id_, {&x0_, &y0_, &radius_}); } -double SurfaceZCylinder::evaluate(Position r) const +double SurfaceZCylinder::_evaluate(Position r) const { return axis_aligned_cylinder_evaluate<0, 1>(r, x0_, y0_, radius_); } @@ -649,7 +674,7 @@ SurfaceSphere::SurfaceSphere(pugi::xml_node surf_node) : CSGSurface(surf_node) read_coeffs(surf_node, id_, {&x0_, &y0_, &z0_, &radius_}); } -double SurfaceSphere::evaluate(Position r) const +double SurfaceSphere::_evaluate(Position r) const { const double x = r.x - x0_; const double y = r.y - y0_; @@ -815,7 +840,7 @@ SurfaceXCone::SurfaceXCone(pugi::xml_node surf_node) : CSGSurface(surf_node) read_coeffs(surf_node, id_, {&x0_, &y0_, &z0_, &radius_sq_}); } -double SurfaceXCone::evaluate(Position r) const +double SurfaceXCone::_evaluate(Position r) const { return axis_aligned_cone_evaluate<0, 1, 2>(r, x0_, y0_, z0_, radius_sq_); } @@ -847,7 +872,7 @@ SurfaceYCone::SurfaceYCone(pugi::xml_node surf_node) : CSGSurface(surf_node) read_coeffs(surf_node, id_, {&x0_, &y0_, &z0_, &radius_sq_}); } -double SurfaceYCone::evaluate(Position r) const +double SurfaceYCone::_evaluate(Position r) const { return axis_aligned_cone_evaluate<1, 0, 2>(r, y0_, x0_, z0_, radius_sq_); } @@ -879,7 +904,7 @@ SurfaceZCone::SurfaceZCone(pugi::xml_node surf_node) : CSGSurface(surf_node) read_coeffs(surf_node, id_, {&x0_, &y0_, &z0_, &radius_sq_}); } -double SurfaceZCone::evaluate(Position r) const +double SurfaceZCone::_evaluate(Position r) const { return axis_aligned_cone_evaluate<2, 0, 1>(r, z0_, x0_, y0_, radius_sq_); } @@ -912,7 +937,7 @@ SurfaceQuadric::SurfaceQuadric(pugi::xml_node surf_node) : CSGSurface(surf_node) surf_node, id_, {&A_, &B_, &C_, &D_, &E_, &F_, &G_, &H_, &J_, &K_}); } -double SurfaceQuadric::evaluate(Position r) const +double SurfaceQuadric::_evaluate(Position r) const { const double x = r.x; const double y = r.y; @@ -1078,7 +1103,7 @@ void SurfaceXTorus::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -double SurfaceXTorus::evaluate(Position r) const +double SurfaceXTorus::_evaluate(Position r) const { double x = r.x - x0_; double y = r.y - y0_; @@ -1131,7 +1156,7 @@ void SurfaceYTorus::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -double SurfaceYTorus::evaluate(Position r) const +double SurfaceYTorus::_evaluate(Position r) const { double x = r.x - x0_; double y = r.y - y0_; @@ -1184,7 +1209,7 @@ void SurfaceZTorus::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -double SurfaceZTorus::evaluate(Position r) const +double SurfaceZTorus::_evaluate(Position r) const { double x = r.x - x0_; double y = r.y - y0_; diff --git a/src/universe.cpp b/src/universe.cpp index b4ef6b4b265..d549dee95be 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -40,10 +40,11 @@ void Universe::to_hdf5(hid_t universes_group) const bool Universe::find_cell(GeometryState& p) const { const auto& cells { - !partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())}; + !partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local(), p.time())}; Position r {p.r_local()}; Position u {p.u_local()}; + double t {p.time()}; auto surf = p.surface(); int32_t i_univ = p.lowest_coord().universe; @@ -51,7 +52,7 @@ bool Universe::find_cell(GeometryState& p) const if (model::cells[i_cell]->universe_ != i_univ) continue; // Check if this cell contains the particle - if (model::cells[i_cell]->contains(r, u, surf)) { + if (model::cells[i_cell]->contains(r, u, t, surf)) { p.lowest_coord().cell = i_cell; return true; } @@ -178,7 +179,7 @@ UniversePartitioner::UniversePartitioner(const Universe& univ) } const vector& UniversePartitioner::get_cells( - Position r, Direction u) const + Position r, Direction u, double t) const { // Perform a binary search for the partition containing the given coordinates. int left = 0; @@ -187,7 +188,7 @@ const vector& UniversePartitioner::get_cells( while (true) { // Check the sense of the coordinates for the current surface. const auto& surf = *model::surfaces[surfs_[middle]]; - if (surf.sense(r, u)) { + if (surf.sense(r, u, t)) { // The coordinates lie in the positive halfspace. Recurse if there are // more surfaces to check. Otherwise, return the cells on the positive // side of this surface. From f94a40097ec6252bfaac8a8e3fe8a26d0baecd40 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Sun, 15 Dec 2024 01:56:40 -0800 Subject: [PATCH 05/14] clang-format --- include/openmc/cell.h | 12 ++++++++---- include/openmc/particle_data.h | 6 +++--- src/cell.cpp | 9 ++++++--- src/surface.cpp | 6 +++--- src/universe.cpp | 5 +++-- 5 files changed, 23 insertions(+), 15 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 30a0bb5d89c..1af863cbc6b 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -111,14 +111,16 @@ class Region { //! Determine if a particle is inside the cell for a simple cell (only //! intersection operators) - bool contains_simple(Position r, Direction u, double t, int32_t on_surface) const; + bool contains_simple( + Position r, Direction u, double t, int32_t on_surface) const; //! Determine if a particle is inside the cell for a complex cell. //! //! Uses the comobination of half-spaces and binary operators to determine //! if short circuiting can be used. Short cicuiting uses the relative and //! absolute depth of parenthases in the expression. - bool contains_complex(Position r, Direction u, double t, int32_t on_surface) const; + bool contains_complex( + Position r, Direction u, double t, int32_t on_surface) const; //! BoundingBox if the paritcle is in a simple cell. BoundingBox bounding_box_simple() const; @@ -182,7 +184,8 @@ class Cell { //! \param on_surface The signed index of a surface that the coordinate is //! known to be on. This index takes precedence over surface sense //! calculations. - virtual bool contains(Position r, Direction u, double t, int32_t on_surface) const = 0; + virtual bool contains( + Position r, Direction u, double t, int32_t on_surface) const = 0; //! Find the oncoming boundary of this cell. virtual std::pair distance( @@ -378,7 +381,8 @@ class CSGCell : public Cell { return region_.distance(r, u, on_surface); } - bool contains(Position r, Direction u, double t, int32_t on_surface) const override + bool contains( + Position r, Direction u, double t, int32_t on_surface) const override { return region_.contains(r, u, t, on_surface); } diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 53b8268a94a..196ee8469a3 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -301,7 +301,7 @@ class GeometryState { const double& time() const { return time_; } double& time_last() { return time_last_; } const double& time_last() const { return time_last_; } - + // Surface that the particle is on int& surface() { return surface_; } const int& surface() const { return surface_; } @@ -343,8 +343,8 @@ class GeometryState { Position r_last_; //!< previous coordinates Direction u_last_; //!< previous direction coordinates - double time_; //!< time - double time_last_; //!< previous time + double time_; //!< time + double time_last_; //!< previous time int surface_ {0}; //!< index for surface particle is on diff --git a/src/cell.cpp b/src/cell.cpp index 9b56a5b6123..8d51dcffde8 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -812,7 +812,8 @@ std::pair Region::distance( //============================================================================== -bool Region::contains(Position r, Direction u, double t, int32_t on_surface) const +bool Region::contains( + Position r, Direction u, double t, int32_t on_surface) const { if (simple_) { return contains_simple(r, u, t, on_surface); @@ -823,7 +824,8 @@ bool Region::contains(Position r, Direction u, double t, int32_t on_surface) con //============================================================================== -bool Region::contains_simple(Position r, Direction u, double t, int32_t on_surface) const +bool Region::contains_simple( + Position r, Direction u, double t, int32_t on_surface) const { for (int32_t token : expression_) { // Assume that no tokens are operators. Evaluate the sense of particle with @@ -846,7 +848,8 @@ bool Region::contains_simple(Position r, Direction u, double t, int32_t on_surfa //============================================================================== -bool Region::contains_complex(Position r, Direction u, double t, int32_t on_surface) const +bool Region::contains_complex( + Position r, Direction u, double t, int32_t on_surface) const { bool in_cell = true; int total_depth = 0; diff --git a/src/surface.cpp b/src/surface.cpp index c64bbfeb2f5..745732cf8f5 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -262,8 +262,8 @@ double CSGSurface::evaluate(Position r, double t) const // The surface moves // Get moving index - int idx = lower_bound_index( - moving_time_grid_.begin(), moving_time_grid_.end(), t); + int idx = + lower_bound_index(moving_time_grid_.begin(), moving_time_grid_.end(), t); // Get moving translation, velocity, and starting time Position translation = moving_translations_[idx]; @@ -273,7 +273,7 @@ double CSGSurface::evaluate(Position r, double t) const // Move the position relative to the surface movement double t_local = t - time_0; Position r_moved = r - (translation + velocity * t_local); - + // Evaluate the moved position return _evaluate(r_moved); } diff --git a/src/universe.cpp b/src/universe.cpp index d549dee95be..3c266f1f839 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -39,8 +39,9 @@ void Universe::to_hdf5(hid_t universes_group) const bool Universe::find_cell(GeometryState& p) const { - const auto& cells { - !partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local(), p.time())}; + const auto& cells {!partitioner_ ? cells_ + : partitioner_->get_cells( + p.r_local(), p.u_local(), p.time())}; Position r {p.r_local()}; Position u {p.u_local()}; From 3cb037e32b2345ab5f0701e3cabd34e03765532e Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Sun, 15 Dec 2024 16:48:29 -0800 Subject: [PATCH 06/14] add dot_normal to consider moving speed. add speed to GeometryState. --- include/openmc/cell.h | 15 +++++++++------ include/openmc/dagmc.h | 1 + include/openmc/particle_data.h | 6 ++++++ include/openmc/surface.h | 13 ++++++++++++- include/openmc/universe.h | 2 +- src/cell.cpp | 14 +++++++------- src/dagmc.cpp | 5 +++++ src/geometry.cpp | 5 +++-- src/surface.cpp | 32 ++++++++++++++++++++++++++++++-- src/universe.cpp | 9 +++++---- 10 files changed, 79 insertions(+), 23 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 1af863cbc6b..9d76de6acf7 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -78,10 +78,11 @@ class Region { //! \param u A direction used to "break ties" the coordinates are very //! close to a surface. //! \param t The time coordinate to check. + //! \param speed Particle speed to "break ties" for moving surface. //! \param on_surface The signed index of a surface that the coordinate is //! known to be on. This index takes precedence over surface sense //! calculations. - bool contains(Position r, Direction u, double t, int32_t on_surface) const; + bool contains(Position r, Direction u, double t, double speed, int32_t on_surface) const; //! Find the oncoming boundary of this cell. std::pair distance( @@ -112,7 +113,8 @@ class Region { //! Determine if a particle is inside the cell for a simple cell (only //! intersection operators) bool contains_simple( - Position r, Direction u, double t, int32_t on_surface) const; + + Position r, Direction u, double t, double speed, int32_t on_surface) const; //! Determine if a particle is inside the cell for a complex cell. //! @@ -120,7 +122,7 @@ class Region { //! if short circuiting can be used. Short cicuiting uses the relative and //! absolute depth of parenthases in the expression. bool contains_complex( - Position r, Direction u, double t, int32_t on_surface) const; + Position r, Direction u, double t, double speed, int32_t on_surface) const; //! BoundingBox if the paritcle is in a simple cell. BoundingBox bounding_box_simple() const; @@ -181,11 +183,12 @@ class Cell { //! \param u A direction used to "break ties" the coordinates are very //! close to a surface. //! \param t The time coordinate to check. + //! \param speed Particle speed to "break ties" for moving surface. //! \param on_surface The signed index of a surface that the coordinate is //! known to be on. This index takes precedence over surface sense //! calculations. virtual bool contains( - Position r, Direction u, double t, int32_t on_surface) const = 0; + Position r, Direction u, double t, double speed, int32_t on_surface) const = 0; //! Find the oncoming boundary of this cell. virtual std::pair distance( @@ -382,9 +385,9 @@ class CSGCell : public Cell { } bool contains( - Position r, Direction u, double t, int32_t on_surface) const override + Position r, Direction u, double t, double speed, int32_t on_surface) const override { - return region_.contains(r, u, t, on_surface); + return region_.contains(r, u, t, speed, on_surface); } BoundingBox bounding_box() const override diff --git a/include/openmc/dagmc.h b/include/openmc/dagmc.h index 14a62e1b4f3..5b4374f44a9 100644 --- a/include/openmc/dagmc.h +++ b/include/openmc/dagmc.h @@ -43,6 +43,7 @@ class DAGSurface : public Surface { double evaluate(Position r, double t) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; + double dot_normal(Position r, Direction u, double t, double speed) const override; Direction reflect(Position r, Direction u, GeometryState* p) const override; inline void to_hdf5_inner(hid_t group_id) const override {}; diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 196ee8469a3..8b281aeb894 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -302,6 +302,10 @@ class GeometryState { double& time_last() { return time_last_; } const double& time_last() const { return time_last_; } + // Accessors for speed (cm/s). + double& speed() { return speed_; } + const double& speed() const { return speed_; } + // Surface that the particle is on int& surface() { return surface_; } const int& surface() const { return surface_; } @@ -345,6 +349,8 @@ class GeometryState { double time_; //!< time double time_last_; //!< previous time + + double speed_; int surface_ {0}; //!< index for surface particle is on diff --git a/include/openmc/surface.h b/include/openmc/surface.h index dfb3cab4aa8..6cd6a43b2b8 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -51,9 +51,10 @@ class Surface { //! \param u A direction used to "break ties" and pick a sense when the //! point is very close to the surface. //! \param t The time for the evaluation + //! \param speed The point speed to "break ties" with moving surface. //! \return true if the point is on the "positive" side of the surface and //! false otherwise. - bool sense(Position r, Direction u, double t) const; + bool sense(Position r, Direction u, double t, double speed) const; //! Determine the direction of a ray reflected from the surface. //! \param[in] r The point at which the ray is incident. @@ -86,6 +87,15 @@ class Surface { //! \param r A 3D Cartesian coordinate. //! \return Normal direction virtual Direction normal(Position r) const = 0; + + //! Compute the dot product of the local outward normal direction of the + //! surface to a geometry coordinate. + //! \param r A 3D Cartesian coordinate. + //! \param u The direction of the ray. + //! \param t The time for the evaluation. + //! \param speed The speed of the particle. + //! \return The dot product + virtual double dot_normal(Position r, Direction u, double t, double speed) const = 0; //! Write all information needed to reconstruct the surface to an HDF5 group. //! \param group_id An HDF5 group id. @@ -109,6 +119,7 @@ class CSGSurface : public Surface { explicit CSGSurface(pugi::xml_node surf_node); CSGSurface(); double evaluate(Position r, double t) const override; + double dot_normal(Position r, Direction u, double t, double speed) const override; protected: //! Static CSG surface evaluation. diff --git a/include/openmc/universe.h b/include/openmc/universe.h index 218b6039d2c..06863e38691 100644 --- a/include/openmc/universe.h +++ b/include/openmc/universe.h @@ -60,7 +60,7 @@ class UniversePartitioner { explicit UniversePartitioner(const Universe& univ); //! Return the list of cells that could contain the given coordinates. - const vector& get_cells(Position r, Direction u, double t) const; + const vector& get_cells(Position r, Direction u, double t, double speed) const; private: //! A sorted vector of indices to surfaces that partition the universe diff --git a/src/cell.cpp b/src/cell.cpp index 8d51dcffde8..b52671f3130 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -813,19 +813,19 @@ std::pair Region::distance( //============================================================================== bool Region::contains( - Position r, Direction u, double t, int32_t on_surface) const + Position r, Direction u, double t, double speed, int32_t on_surface) const { if (simple_) { - return contains_simple(r, u, t, on_surface); + return contains_simple(r, u, t, speed, on_surface); } else { - return contains_complex(r, u, t, on_surface); + return contains_complex(r, u, t, speed, on_surface); } } //============================================================================== bool Region::contains_simple( - Position r, Direction u, double t, int32_t on_surface) const + Position r, Direction u, double t, double speed, int32_t on_surface) const { for (int32_t token : expression_) { // Assume that no tokens are operators. Evaluate the sense of particle with @@ -837,7 +837,7 @@ bool Region::contains_simple( return false; } else { // Note the off-by-one indexing - bool sense = model::surfaces[abs(token) - 1]->sense(r, u, t); + bool sense = model::surfaces[abs(token) - 1]->sense(r, u, t, speed); if (sense != (token > 0)) { return false; } @@ -849,7 +849,7 @@ bool Region::contains_simple( //============================================================================== bool Region::contains_complex( - Position r, Direction u, double t, int32_t on_surface) const + Position r, Direction u, double t, double speed, int32_t on_surface) const { bool in_cell = true; int total_depth = 0; @@ -868,7 +868,7 @@ bool Region::contains_complex( in_cell = false; } else { // Note the off-by-one indexing - bool sense = model::surfaces[abs(token) - 1]->sense(r, u, t); + bool sense = model::surfaces[abs(token) - 1]->sense(r, u, t, speed); in_cell = (sense == (token > 0)); } } else if ((token == OP_UNION && in_cell == true) || diff --git a/src/dagmc.cpp b/src/dagmc.cpp index fa86631f758..907aca93be7 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -758,6 +758,11 @@ Direction DAGSurface::normal(Position r) const return dir; } +double DAGSurface::dot_normal(Position r, Direction u, double t, double speed) const +{ + return u.dot(normal(r)); +} + Direction DAGSurface::reflect(Position r, Direction u, GeometryState* p) const { Expects(p); diff --git a/src/geometry.cpp b/src/geometry.cpp index 0cf954a6557..c33618aed3e 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -43,7 +43,7 @@ bool check_cell_overlap(GeometryState& p, bool error) // Loop through each cell on this level for (auto index_cell : univ.cells_) { Cell& c = *model::cells[index_cell]; - if (c.contains(p.coord(j).r, p.coord(j).u, p.time(), p.surface())) { + if (c.contains(p.coord(j).r, p.coord(j).u, p.time(), p.surface(), p.speed())) { if (index_cell != p.coord(j).cell) { if (error) { fatal_error( @@ -119,8 +119,9 @@ bool find_cell_inner( Position r {p.r_local()}; Direction u {p.u_local()}; double t {p.time()}; + double speed {p.speed()}; auto surf = p.surface(); - if (model::cells[i_cell]->contains(r, u, t, surf)) { + if (model::cells[i_cell]->contains(r, u, t, speed, surf)) { p.lowest_coord().cell = i_cell; found = true; break; diff --git a/src/surface.cpp b/src/surface.cpp index 745732cf8f5..acdbda11468 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -115,7 +115,7 @@ Surface::Surface(pugi::xml_node surf_node) } } -bool Surface::sense(Position r, Direction u, double t) const +bool Surface::sense(Position r, Direction u, double t, double speed) const { // Evaluate the surface equation at the particle's coordinates to determine // which side the particle is on. @@ -126,7 +126,7 @@ bool Surface::sense(Position r, Direction u, double t) const // Particle may be coincident with this surface. To determine the sense, we // look at the direction of the particle relative to the surface normal (by // default in the positive direction) via their dot product. - return u.dot(normal(r)) > 0.0; + return dot_normal(r, u, t, speed) > 0.0; } return f > 0.0; } @@ -278,6 +278,34 @@ double CSGSurface::evaluate(Position r, double t) const return _evaluate(r_moved); } +double CSGSurface::dot_normal(Position r, Direction u, double t, double speed) const +{ + if (!moving_) { + return u.dot(normal(r)); + } + // The surface moves + + // Get moving index + int idx = + lower_bound_index(moving_time_grid_.begin(), moving_time_grid_.end(), t); + + // Get moving translation, velocity, and starting time + Position translation = moving_translations_[idx]; + double time_0 = moving_time_grid_[idx]; + Position velocity = moving_velocities_[idx]; + + // Move the position relative to the surface movement + double t_local = t - time_0; + Position r_moved = r - (translation + velocity * t_local); + + // Get the relative direction + Direction u_relative = u - velocity / speed; + + // Get the dot product + return u_relative.dot(normal(r_moved)); +} + + //============================================================================== // Generic functions for x-, y-, and z-, planes. //============================================================================== diff --git a/src/universe.cpp b/src/universe.cpp index 3c266f1f839..4a0f847910e 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -41,11 +41,12 @@ bool Universe::find_cell(GeometryState& p) const { const auto& cells {!partitioner_ ? cells_ : partitioner_->get_cells( - p.r_local(), p.u_local(), p.time())}; + p.r_local(), p.u_local(), p.time(), p.speed())}; Position r {p.r_local()}; Position u {p.u_local()}; double t {p.time()}; + double speed {p.speed()}; auto surf = p.surface(); int32_t i_univ = p.lowest_coord().universe; @@ -53,7 +54,7 @@ bool Universe::find_cell(GeometryState& p) const if (model::cells[i_cell]->universe_ != i_univ) continue; // Check if this cell contains the particle - if (model::cells[i_cell]->contains(r, u, t, surf)) { + if (model::cells[i_cell]->contains(r, u, t, speed, surf)) { p.lowest_coord().cell = i_cell; return true; } @@ -180,7 +181,7 @@ UniversePartitioner::UniversePartitioner(const Universe& univ) } const vector& UniversePartitioner::get_cells( - Position r, Direction u, double t) const + Position r, Direction u, double t, double speed) const { // Perform a binary search for the partition containing the given coordinates. int left = 0; @@ -189,7 +190,7 @@ const vector& UniversePartitioner::get_cells( while (true) { // Check the sense of the coordinates for the current surface. const auto& surf = *model::surfaces[surfs_[middle]]; - if (surf.sense(r, u, t)) { + if (surf.sense(r, u, t, speed)) { // The coordinates lie in the positive halfspace. Recurse if there are // more surfaces to check. Otherwise, return the cells on the positive // side of this surface. From 9edd2a2b92d32dd4d04f13ed1abb12cd92e84292 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Sun, 15 Dec 2024 16:48:44 -0800 Subject: [PATCH 07/14] clang-format --- include/openmc/cell.h | 11 ++++++----- include/openmc/dagmc.h | 3 ++- include/openmc/particle_data.h | 2 +- include/openmc/surface.h | 10 ++++++---- include/openmc/universe.h | 3 ++- src/dagmc.cpp | 3 ++- src/geometry.cpp | 3 ++- src/surface.cpp | 6 +++--- src/universe.cpp | 4 ++-- 9 files changed, 26 insertions(+), 19 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 9d76de6acf7..c0b034a0ca6 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -82,7 +82,8 @@ class Region { //! \param on_surface The signed index of a surface that the coordinate is //! known to be on. This index takes precedence over surface sense //! calculations. - bool contains(Position r, Direction u, double t, double speed, int32_t on_surface) const; + bool contains( + Position r, Direction u, double t, double speed, int32_t on_surface) const; //! Find the oncoming boundary of this cell. std::pair distance( @@ -187,8 +188,8 @@ class Cell { //! \param on_surface The signed index of a surface that the coordinate is //! known to be on. This index takes precedence over surface sense //! calculations. - virtual bool contains( - Position r, Direction u, double t, double speed, int32_t on_surface) const = 0; + virtual bool contains(Position r, Direction u, double t, double speed, + int32_t on_surface) const = 0; //! Find the oncoming boundary of this cell. virtual std::pair distance( @@ -384,8 +385,8 @@ class CSGCell : public Cell { return region_.distance(r, u, on_surface); } - bool contains( - Position r, Direction u, double t, double speed, int32_t on_surface) const override + bool contains(Position r, Direction u, double t, double speed, + int32_t on_surface) const override { return region_.contains(r, u, t, speed, on_surface); } diff --git a/include/openmc/dagmc.h b/include/openmc/dagmc.h index 5b4374f44a9..263996d694f 100644 --- a/include/openmc/dagmc.h +++ b/include/openmc/dagmc.h @@ -43,7 +43,8 @@ class DAGSurface : public Surface { double evaluate(Position r, double t) const override; double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; - double dot_normal(Position r, Direction u, double t, double speed) const override; + double dot_normal( + Position r, Direction u, double t, double speed) const override; Direction reflect(Position r, Direction u, GeometryState* p) const override; inline void to_hdf5_inner(hid_t group_id) const override {}; diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 8b281aeb894..9e524eb8dc0 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -349,7 +349,7 @@ class GeometryState { double time_; //!< time double time_last_; //!< previous time - + double speed_; int surface_ {0}; //!< index for surface particle is on diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 6cd6a43b2b8..82fa8187ddd 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -87,15 +87,16 @@ class Surface { //! \param r A 3D Cartesian coordinate. //! \return Normal direction virtual Direction normal(Position r) const = 0; - - //! Compute the dot product of the local outward normal direction of the + + //! Compute the dot product of the local outward normal direction of the //! surface to a geometry coordinate. //! \param r A 3D Cartesian coordinate. //! \param u The direction of the ray. //! \param t The time for the evaluation. //! \param speed The speed of the particle. //! \return The dot product - virtual double dot_normal(Position r, Direction u, double t, double speed) const = 0; + virtual double dot_normal( + Position r, Direction u, double t, double speed) const = 0; //! Write all information needed to reconstruct the surface to an HDF5 group. //! \param group_id An HDF5 group id. @@ -119,7 +120,8 @@ class CSGSurface : public Surface { explicit CSGSurface(pugi::xml_node surf_node); CSGSurface(); double evaluate(Position r, double t) const override; - double dot_normal(Position r, Direction u, double t, double speed) const override; + double dot_normal( + Position r, Direction u, double t, double speed) const override; protected: //! Static CSG surface evaluation. diff --git a/include/openmc/universe.h b/include/openmc/universe.h index 06863e38691..f93e4800094 100644 --- a/include/openmc/universe.h +++ b/include/openmc/universe.h @@ -60,7 +60,8 @@ class UniversePartitioner { explicit UniversePartitioner(const Universe& univ); //! Return the list of cells that could contain the given coordinates. - const vector& get_cells(Position r, Direction u, double t, double speed) const; + const vector& get_cells( + Position r, Direction u, double t, double speed) const; private: //! A sorted vector of indices to surfaces that partition the universe diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 907aca93be7..88cef3ef73c 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -758,7 +758,8 @@ Direction DAGSurface::normal(Position r) const return dir; } -double DAGSurface::dot_normal(Position r, Direction u, double t, double speed) const +double DAGSurface::dot_normal( + Position r, Direction u, double t, double speed) const { return u.dot(normal(r)); } diff --git a/src/geometry.cpp b/src/geometry.cpp index c33618aed3e..cc73611d29c 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -43,7 +43,8 @@ bool check_cell_overlap(GeometryState& p, bool error) // Loop through each cell on this level for (auto index_cell : univ.cells_) { Cell& c = *model::cells[index_cell]; - if (c.contains(p.coord(j).r, p.coord(j).u, p.time(), p.surface(), p.speed())) { + if (c.contains( + p.coord(j).r, p.coord(j).u, p.time(), p.surface(), p.speed())) { if (index_cell != p.coord(j).cell) { if (error) { fatal_error( diff --git a/src/surface.cpp b/src/surface.cpp index acdbda11468..5d0234e679a 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -278,13 +278,14 @@ double CSGSurface::evaluate(Position r, double t) const return _evaluate(r_moved); } -double CSGSurface::dot_normal(Position r, Direction u, double t, double speed) const +double CSGSurface::dot_normal( + Position r, Direction u, double t, double speed) const { if (!moving_) { return u.dot(normal(r)); } // The surface moves - + // Get moving index int idx = lower_bound_index(moving_time_grid_.begin(), moving_time_grid_.end(), t); @@ -305,7 +306,6 @@ double CSGSurface::dot_normal(Position r, Direction u, double t, double speed) c return u_relative.dot(normal(r_moved)); } - //============================================================================== // Generic functions for x-, y-, and z-, planes. //============================================================================== diff --git a/src/universe.cpp b/src/universe.cpp index 4a0f847910e..89cbe09dd93 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -40,8 +40,8 @@ void Universe::to_hdf5(hid_t universes_group) const bool Universe::find_cell(GeometryState& p) const { const auto& cells {!partitioner_ ? cells_ - : partitioner_->get_cells( - p.r_local(), p.u_local(), p.time(), p.speed())}; + : partitioner_->get_cells(p.r_local(), + p.u_local(), p.time(), p.speed())}; Position r {p.r_local()}; Position u {p.u_local()}; From c4f70f01f137e22ae3be6044c7bce329aaf4562c Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Sun, 15 Dec 2024 16:59:14 -0800 Subject: [PATCH 08/14] rename speed() method as get_speed() --- include/openmc/particle.h | 2 +- src/particle.cpp | 4 +++- src/simulation.cpp | 2 ++ 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 6a2e67049fd..f08a02ae318 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -37,7 +37,7 @@ class Particle : public ParticleData { //========================================================================== // Methods - double speed() const; + double get_speed() const; //! moves the particle by the distance length to its next location //! \param length the distance the particle is moved diff --git a/src/particle.cpp b/src/particle.cpp index 64c50c9438f..34771a40c65 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -42,7 +42,7 @@ namespace openmc { // Particle implementation //============================================================================== -double Particle::speed() const +double Particle::get_speed() const { // Determine mass in eV/c^2 double mass; @@ -150,6 +150,8 @@ void Particle::event_calculate_xs() r_last() = r(); time_last() = time(); + speed() = get_speed(); + // Reset event variables event() = TallyEvent::KILL; event_nuclide() = NUCLIDE_NONE; diff --git a/src/simulation.cpp b/src/simulation.cpp index 09b3764732b..9e5a48beb2a 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -572,6 +572,8 @@ void initialize_history(Particle& p, int64_t index_source) } p.current_work() = index_source; + p.speed() = p.get_speed(); + // set identifier for particle p.id() = simulation::work_index[mpi::rank] + index_source; From 51531934f526ac9cf333a529c2ccc4a8e779d41f Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Wed, 23 Apr 2025 11:17:39 -0700 Subject: [PATCH 09/14] reapply moving surface from ilhamv repo --- include/openmc/cell.h | 8 ++-- include/openmc/surface.h | 38 +++++++++------- src/boundary_condition.cpp | 8 ++-- src/cell.cpp | 4 +- src/geometry.cpp | 4 +- src/plot.cpp | 2 +- src/surface.cpp | 92 +++++++++++++++++++++++++++++++------- 7 files changed, 112 insertions(+), 44 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index c0b034a0ca6..8dedb68600e 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -87,7 +87,7 @@ class Region { //! Find the oncoming boundary of this cell. std::pair distance( - Position r, Direction u, int32_t on_surface) const; + Position r, Direction u, double t, double speed, int32_t on_surface) const; //! Get the BoundingBox for this cell. BoundingBox bounding_box(int32_t cell_id) const; @@ -193,7 +193,7 @@ class Cell { //! Find the oncoming boundary of this cell. virtual std::pair distance( - Position r, Direction u, int32_t on_surface, GeometryState* p) const = 0; + Position r, Direction u, double t, double speed, int32_t on_surface, GeometryState* p) const = 0; //! Write all information needed to reconstruct the cell to an HDF5 group. //! \param group_id An HDF5 group id. @@ -380,9 +380,9 @@ class CSGCell : public Cell { vector surfaces() const override { return region_.surfaces(); } std::pair distance(Position r, Direction u, - int32_t on_surface, GeometryState* p) const override + double t, double speed, int32_t on_surface, GeometryState* p) const override { - return region_.distance(r, u, on_surface); + return region_.distance(r, u, t, speed, on_surface); } bool contains(Position r, Direction u, double t, double speed, diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 82fa8187ddd..9309010e884 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -79,9 +79,11 @@ class Surface { //! Compute the distance between a point and the surface along a ray. //! \param r A 3D Cartesian coordinate. //! \param u The direction of the ray. + //! \param t The time of the moving point (for moving surface). + //! \param speed The speed of moving point (for moving surface). //! \param coincident A hint to the code that the given point should lie //! exactly on the surface. - virtual double distance(Position r, Direction u, bool coincident) const = 0; + virtual double distance(Position r, Direction u, double time, double speed, bool coincident) const = 0; //! Compute the local outward normal direction of the surface. //! \param r A 3D Cartesian coordinate. @@ -122,10 +124,12 @@ class CSGSurface : public Surface { double evaluate(Position r, double t) const override; double dot_normal( Position r, Direction u, double t, double speed) const override; + double distance(Position r, Direction u, double time, double speed, bool coincident) const override; protected: - //! Static CSG surface evaluation. + //! Static CSG surface functions virtual double _evaluate(Position r) const = 0; + virtual double _distance(Position r, Direction u, bool coincident) const = 0; }; //============================================================================== @@ -137,7 +141,6 @@ class CSGSurface : public Surface { class SurfaceXPlane : public CSGSurface { public: explicit SurfaceXPlane(pugi::xml_node surf_node); - double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; BoundingBox bounding_box(bool pos_side) const override; @@ -146,6 +149,7 @@ class SurfaceXPlane : public CSGSurface { protected: double _evaluate(Position r) const override; + double _distance(Position r, Direction u, bool coincident) const override; }; //============================================================================== @@ -157,7 +161,6 @@ class SurfaceXPlane : public CSGSurface { class SurfaceYPlane : public CSGSurface { public: explicit SurfaceYPlane(pugi::xml_node surf_node); - double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; BoundingBox bounding_box(bool pos_side) const override; @@ -166,6 +169,7 @@ class SurfaceYPlane : public CSGSurface { protected: double _evaluate(Position r) const override; + double _distance(Position r, Direction u, bool coincident) const override; }; //============================================================================== @@ -177,7 +181,6 @@ class SurfaceYPlane : public CSGSurface { class SurfaceZPlane : public CSGSurface { public: explicit SurfaceZPlane(pugi::xml_node surf_node); - double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; BoundingBox bounding_box(bool pos_side) const override; @@ -186,6 +189,7 @@ class SurfaceZPlane : public CSGSurface { protected: double _evaluate(Position r) const override; + double _distance(Position r, Direction u, bool coincident) const override; }; //============================================================================== @@ -197,7 +201,6 @@ class SurfaceZPlane : public CSGSurface { class SurfacePlane : public CSGSurface { public: explicit SurfacePlane(pugi::xml_node surf_node); - double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; @@ -205,6 +208,7 @@ class SurfacePlane : public CSGSurface { protected: double _evaluate(Position r) const override; + double _distance(Position r, Direction u, bool coincident) const override; }; //============================================================================== @@ -217,7 +221,6 @@ class SurfacePlane : public CSGSurface { class SurfaceXCylinder : public CSGSurface { public: explicit SurfaceXCylinder(pugi::xml_node surf_node); - double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; BoundingBox bounding_box(bool pos_side) const override; @@ -226,6 +229,7 @@ class SurfaceXCylinder : public CSGSurface { protected: double _evaluate(Position r) const override; + double _distance(Position r, Direction u, bool coincident) const override; }; //============================================================================== @@ -238,7 +242,6 @@ class SurfaceXCylinder : public CSGSurface { class SurfaceYCylinder : public CSGSurface { public: explicit SurfaceYCylinder(pugi::xml_node surf_node); - double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; BoundingBox bounding_box(bool pos_side) const override; @@ -247,6 +250,7 @@ class SurfaceYCylinder : public CSGSurface { protected: double _evaluate(Position r) const override; + double _distance(Position r, Direction u, bool coincident) const override; }; //============================================================================== @@ -259,7 +263,6 @@ class SurfaceYCylinder : public CSGSurface { class SurfaceZCylinder : public CSGSurface { public: explicit SurfaceZCylinder(pugi::xml_node surf_node); - double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; BoundingBox bounding_box(bool pos_side) const override; @@ -268,6 +271,7 @@ class SurfaceZCylinder : public CSGSurface { protected: double _evaluate(Position r) const override; + double _distance(Position r, Direction u, bool coincident) const override; }; //============================================================================== @@ -280,7 +284,6 @@ class SurfaceZCylinder : public CSGSurface { class SurfaceSphere : public CSGSurface { public: explicit SurfaceSphere(pugi::xml_node surf_node); - double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; BoundingBox bounding_box(bool pos_side) const override; @@ -289,6 +292,7 @@ class SurfaceSphere : public CSGSurface { protected: double _evaluate(Position r) const override; + double _distance(Position r, Direction u, bool coincident) const override; }; //============================================================================== @@ -301,7 +305,6 @@ class SurfaceSphere : public CSGSurface { class SurfaceXCone : public CSGSurface { public: explicit SurfaceXCone(pugi::xml_node surf_node); - double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; @@ -309,6 +312,7 @@ class SurfaceXCone : public CSGSurface { protected: double _evaluate(Position r) const override; + double _distance(Position r, Direction u, bool coincident) const override; }; //============================================================================== @@ -321,7 +325,6 @@ class SurfaceXCone : public CSGSurface { class SurfaceYCone : public CSGSurface { public: explicit SurfaceYCone(pugi::xml_node surf_node); - double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; @@ -329,6 +332,7 @@ class SurfaceYCone : public CSGSurface { protected: double _evaluate(Position r) const override; + double _distance(Position r, Direction u, bool coincident) const override; }; //============================================================================== @@ -341,7 +345,6 @@ class SurfaceYCone : public CSGSurface { class SurfaceZCone : public CSGSurface { public: explicit SurfaceZCone(pugi::xml_node surf_node); - double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; @@ -349,6 +352,7 @@ class SurfaceZCone : public CSGSurface { protected: double _evaluate(Position r) const override; + double _distance(Position r, Direction u, bool coincident) const override; }; //============================================================================== @@ -361,7 +365,6 @@ class SurfaceZCone : public CSGSurface { class SurfaceQuadric : public CSGSurface { public: explicit SurfaceQuadric(pugi::xml_node surf_node); - double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; @@ -370,6 +373,7 @@ class SurfaceQuadric : public CSGSurface { protected: double _evaluate(Position r) const override; + double _distance(Position r, Direction u, bool coincident) const override; }; //============================================================================== @@ -381,7 +385,6 @@ class SurfaceQuadric : public CSGSurface { class SurfaceXTorus : public CSGSurface { public: explicit SurfaceXTorus(pugi::xml_node surf_node); - double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; @@ -389,6 +392,7 @@ class SurfaceXTorus : public CSGSurface { protected: double _evaluate(Position r) const override; + double _distance(Position r, Direction u, bool coincident) const override; }; //============================================================================== @@ -400,7 +404,6 @@ class SurfaceXTorus : public CSGSurface { class SurfaceYTorus : public CSGSurface { public: explicit SurfaceYTorus(pugi::xml_node surf_node); - double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; @@ -408,6 +411,7 @@ class SurfaceYTorus : public CSGSurface { protected: double _evaluate(Position r) const override; + double _distance(Position r, Direction u, bool coincident) const override; }; //============================================================================== @@ -419,7 +423,6 @@ class SurfaceYTorus : public CSGSurface { class SurfaceZTorus : public CSGSurface { public: explicit SurfaceZTorus(pugi::xml_node surf_node); - double distance(Position r, Direction u, bool coincident) const override; Direction normal(Position r) const override; void to_hdf5_inner(hid_t group_id) const override; @@ -427,6 +430,7 @@ class SurfaceZTorus : public CSGSurface { protected: double _evaluate(Position r) const override; + double _distance(Position r, Direction u, bool coincident) const override; }; //============================================================================== diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 74d0b7604b3..f726140ef86 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -103,9 +103,9 @@ TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf) double d1; double e1 = surf1.evaluate(origin, 0.0); // Time = 0.0 if (e1 > FP_COINCIDENT) { - d1 = -surf1.distance(origin, -u, false); + d1 = -surf1.distance(origin, -u, 0.0, 0.0, false); // Time = speed = 0.0 } else if (e1 < -FP_COINCIDENT) { - d1 = surf1.distance(origin, u, false); + d1 = surf1.distance(origin, u, 0.0, 0.0, false); // Time = speed = 0.0 } else { d1 = 0.0; } @@ -114,9 +114,9 @@ TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf) double d2; double e2 = surf2.evaluate(origin, 0.0); // Time = 0.0 if (e2 > FP_COINCIDENT) { - d2 = -surf2.distance(origin, -u, false); + d2 = -surf2.distance(origin, -u, 0.0, 0.0, false); // Time = speed = 0.0 } else if (e2 < -FP_COINCIDENT) { - d2 = surf2.distance(origin, u, false); + d2 = surf2.distance(origin, u, 0.0, 0.0, false); // Time = speed = 0.0 } else { d2 = 0.0; } diff --git a/src/cell.cpp b/src/cell.cpp index b52671f3130..a274de37c73 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -783,7 +783,7 @@ std::string Region::str() const //============================================================================== std::pair Region::distance( - Position r, Direction u, int32_t on_surface) const + Position r, Direction u, double t, double speed, int32_t on_surface) const { double min_dist {INFTY}; int32_t i_surf {std::numeric_limits::max()}; @@ -796,7 +796,7 @@ std::pair Region::distance( // Calculate the distance to this surface. // Note the off-by-one indexing bool coincident {std::abs(token) == std::abs(on_surface)}; - double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)}; + double d {model::surfaces[abs(token) - 1]->distance(r, u, t, speed, coincident)}; // Check if this distance is the new minimum. if (d < min_dist) { diff --git a/src/geometry.cpp b/src/geometry.cpp index cc73611d29c..2b275e48c27 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -372,10 +372,12 @@ BoundaryInfo distance_to_boundary(GeometryState& p) const auto& coord {p.coord(i)}; const Position& r {coord.r}; const Direction& u {coord.u}; + const double t {p.time()}; + const double speed {p.speed()}; Cell& c {*model::cells[coord.cell]}; // Find the oncoming surface in this cell and the distance to it. - auto surface_distance = c.distance(r, u, p.surface(), &p); + auto surface_distance = c.distance(r, u, t, speed, p.surface(), &p); d_surf = surface_distance.first; level_surf_cross = surface_distance.second; diff --git a/src/plot.cpp b/src/plot.cpp index 348138570c1..4a793788b60 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -1098,7 +1098,7 @@ int ProjectionPlot::advance_to_boundary_from_void(GeometryState& p) Universe* uni = model::universes[model::root_universe].get(); int intersected_surface = -1; for (auto c_i : uni->cells_) { - auto dist = model::cells.at(c_i)->distance(coord.r, coord.u, 0, &p); + auto dist = model::cells.at(c_i)->distance(coord.r, coord.u, p.time(), p.speed(), 0, &p); if (dist.first < min_dist) { min_dist = dist.first; intersected_surface = dist.second; diff --git a/src/surface.cpp b/src/surface.cpp index 5d0234e679a..7177ec84189 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -306,6 +306,68 @@ double CSGSurface::dot_normal( return u_relative.dot(normal(r_moved)); } +double CSGSurface::distance(Position r, Direction u, double t, double speed, bool coincident) const +{ + if (!moving_) { + return _distance(r, u, coincident); + } + // The surface moves + + // Store the origin coordinate + Position r_origin {r}; + Position u_origin {u}; + double t_origin {t}; + + // Get moving interval index + int idx = + lower_bound_index(moving_time_grid_.begin(), moving_time_grid_.end(), t); + + // Distance accumulator + double distance_total = 0.0; + + // Evaluate the current and the subsequent intervals until intersecting + while (idx < moving_velocities_.size()) { + // Get moving translation, velocity, and starting time + Position translation = moving_translations_[idx]; + double time_0 = moving_time_grid_[idx]; + Position velocity = moving_velocities_[idx]; + + // Adjust the position based on the surface movement + double t_local = t - time_0; + r -= (translation + velocity * t_local); + + // Adjust to relative direction + u -= velocity / speed; + + // Get distance using the static function based on the adjusted position + // and direction + double distance = _distance(r, u, coincident); + + // Intersection within the interval? + double t_distance = distance / speed; + double t_interval = moving_time_grid_[idx + 1] - t; + if (t_distance < t_interval) { + // Return the total distance + return distance_total + distance; + } + // Not intersecting. + // Prepare to check the next interval. + + // Accumulate distance + distance_total += t_interval * speed; + + // March forward the coordinate + r = r_origin + distance_total * u_origin; + u = u_origin; + t = moving_time_grid_[idx + 1]; + + idx++; + } + + // No intersection + return INFTY; +} + //============================================================================== // Generic functions for x-, y-, and z-, planes. //============================================================================== @@ -338,7 +400,7 @@ double SurfaceXPlane::_evaluate(Position r) const return r.x - x0_; } -double SurfaceXPlane::distance(Position r, Direction u, bool coincident) const +double SurfaceXPlane::_distance(Position r, Direction u, bool coincident) const { return axis_aligned_plane_distance<0>(r, u, coincident, x0_); } @@ -378,7 +440,7 @@ double SurfaceYPlane::_evaluate(Position r) const return r.y - y0_; } -double SurfaceYPlane::distance(Position r, Direction u, bool coincident) const +double SurfaceYPlane::_distance(Position r, Direction u, bool coincident) const { return axis_aligned_plane_distance<1>(r, u, coincident, y0_); } @@ -418,7 +480,7 @@ double SurfaceZPlane::_evaluate(Position r) const return r.z - z0_; } -double SurfaceZPlane::distance(Position r, Direction u, bool coincident) const +double SurfaceZPlane::_distance(Position r, Direction u, bool coincident) const { return axis_aligned_plane_distance<2>(r, u, coincident, z0_); } @@ -458,7 +520,7 @@ double SurfacePlane::_evaluate(Position r) const return A_ * r.x + B_ * r.y + C_ * r.z - D_; } -double SurfacePlane::distance(Position r, Direction u, bool coincident) const +double SurfacePlane::_distance(Position r, Direction u, bool coincident) const { const double f = A_ * r.x + B_ * r.y + C_ * r.z - D_; const double projection = A_ * u.x + B_ * u.y + C_ * u.z; @@ -577,7 +639,7 @@ double SurfaceXCylinder::_evaluate(Position r) const return axis_aligned_cylinder_evaluate<1, 2>(r, y0_, z0_, radius_); } -double SurfaceXCylinder::distance( +double SurfaceXCylinder::_distance( Position r, Direction u, bool coincident) const { return axis_aligned_cylinder_distance<0, 1, 2>( @@ -620,7 +682,7 @@ double SurfaceYCylinder::_evaluate(Position r) const return axis_aligned_cylinder_evaluate<0, 2>(r, x0_, z0_, radius_); } -double SurfaceYCylinder::distance( +double SurfaceYCylinder::_distance( Position r, Direction u, bool coincident) const { return axis_aligned_cylinder_distance<1, 0, 2>( @@ -664,7 +726,7 @@ double SurfaceZCylinder::_evaluate(Position r) const return axis_aligned_cylinder_evaluate<0, 1>(r, x0_, y0_, radius_); } -double SurfaceZCylinder::distance( +double SurfaceZCylinder::_distance( Position r, Direction u, bool coincident) const { return axis_aligned_cylinder_distance<2, 0, 1>( @@ -710,7 +772,7 @@ double SurfaceSphere::_evaluate(Position r) const return x * x + y * y + z * z - radius_ * radius_; } -double SurfaceSphere::distance(Position r, Direction u, bool coincident) const +double SurfaceSphere::_distance(Position r, Direction u, bool coincident) const { const double x = r.x - x0_; const double y = r.y - y0_; @@ -873,7 +935,7 @@ double SurfaceXCone::_evaluate(Position r) const return axis_aligned_cone_evaluate<0, 1, 2>(r, x0_, y0_, z0_, radius_sq_); } -double SurfaceXCone::distance(Position r, Direction u, bool coincident) const +double SurfaceXCone::_distance(Position r, Direction u, bool coincident) const { return axis_aligned_cone_distance<0, 1, 2>( r, u, coincident, x0_, y0_, z0_, radius_sq_); @@ -905,7 +967,7 @@ double SurfaceYCone::_evaluate(Position r) const return axis_aligned_cone_evaluate<1, 0, 2>(r, y0_, x0_, z0_, radius_sq_); } -double SurfaceYCone::distance(Position r, Direction u, bool coincident) const +double SurfaceYCone::_distance(Position r, Direction u, bool coincident) const { return axis_aligned_cone_distance<1, 0, 2>( r, u, coincident, y0_, x0_, z0_, radius_sq_); @@ -937,7 +999,7 @@ double SurfaceZCone::_evaluate(Position r) const return axis_aligned_cone_evaluate<2, 0, 1>(r, z0_, x0_, y0_, radius_sq_); } -double SurfaceZCone::distance(Position r, Direction u, bool coincident) const +double SurfaceZCone::_distance(Position r, Direction u, bool coincident) const { return axis_aligned_cone_distance<2, 0, 1>( r, u, coincident, z0_, x0_, y0_, radius_sq_); @@ -974,7 +1036,7 @@ double SurfaceQuadric::_evaluate(Position r) const z * (C_ * z + F_ * x + J_) + K_; } -double SurfaceQuadric::distance( +double SurfaceQuadric::_distance( Position r, Direction ang, bool coincident) const { const double& x = r.x; @@ -1140,7 +1202,7 @@ double SurfaceXTorus::_evaluate(Position r) const std::pow(std::sqrt(y * y + z * z) - A_, 2) / (C_ * C_) - 1.; } -double SurfaceXTorus::distance(Position r, Direction u, bool coincident) const +double SurfaceXTorus::_distance(Position r, Direction u, bool coincident) const { double x = r.x - x0_; double y = r.y - y0_; @@ -1193,7 +1255,7 @@ double SurfaceYTorus::_evaluate(Position r) const std::pow(std::sqrt(x * x + z * z) - A_, 2) / (C_ * C_) - 1.; } -double SurfaceYTorus::distance(Position r, Direction u, bool coincident) const +double SurfaceYTorus::_distance(Position r, Direction u, bool coincident) const { double x = r.x - x0_; double y = r.y - y0_; @@ -1246,7 +1308,7 @@ double SurfaceZTorus::_evaluate(Position r) const std::pow(std::sqrt(x * x + y * y) - A_, 2) / (C_ * C_) - 1.; } -double SurfaceZTorus::distance(Position r, Direction u, bool coincident) const +double SurfaceZTorus::_distance(Position r, Direction u, bool coincident) const { double x = r.x - x0_; double y = r.y - y0_; From 1ad03d83e56c8e82771542cbe15e0f2c72ad296a Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Wed, 4 Jun 2025 17:26:15 +0700 Subject: [PATCH 10/14] corrected coincidence treatment --- src/cell.cpp | 5 +---- src/surface.cpp | 17 +++++++++++++++-- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index 6b1b53054f4..5c6c6cec35c 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -827,7 +827,6 @@ bool Region::contains_simple( // is moving, we need to factor in the surface-particle relative velocity, // so the surface-sense calculation needs to be done. - /* if (model::surfaces[abs(token) - 1]->moving_) { // Note the off-by-one indexing bool sense = model::surfaces[abs(token) - 1]->sense(r, u, t, speed); @@ -837,7 +836,6 @@ bool Region::contains_simple( continue; } // Surface is not moving - */ if (token == on_surface) { } else if (-token == on_surface) { @@ -869,8 +867,7 @@ bool Region::contains_complex( // If the token is a union or intersection check to // short circuit if (token < OP_UNION) { - //if (model::surfaces[abs(token) - 1]->moving_) { - if (false) { + if (model::surfaces[abs(token) - 1]->moving_) { // Note the off-by-one indexing bool sense = model::surfaces[abs(token) - 1]->sense(r, u, t, speed); in_cell = (sense == (token > 0)); diff --git a/src/surface.cpp b/src/surface.cpp index 2a0b551d74c..0b4b3ebe949 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -247,6 +247,8 @@ double Surface::distance(Position r, Direction u, double t, double speed, bool c return _distance(r, u, coincident); } // The surface moves + + coincident = false; // Store the origin coordinate Position r_origin {r}; @@ -273,10 +275,17 @@ double Surface::distance(Position r, Direction u, double t, double speed, bool c // Adjust to relative direction u -= velocity / speed; + + // Normalize scaled relative direction + double norm = u.norm(); + u /= norm; // Get distance using the static function based on the adjusted position // and direction - double distance = _distance(r, u, coincident); + double distance = _distance(r, u, coincident) * norm; + + // Rescale relative direction + u *= norm; // Intersection within the interval? double t_distance = distance / speed; @@ -326,7 +335,11 @@ double Surface::dot_normal( // Get the relative direction Direction u_relative = u - velocity / speed; - + + // Normalize scaled relative direction + double norm = u_relative.norm(); + u_relative /= norm; + // Get the dot product return u_relative.dot(normal(r_moved)); } From 54efcd9aea3adbec3f9d9079ab75457ea4acc2c8 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Mon, 9 Jun 2025 13:16:12 +0700 Subject: [PATCH 11/14] add time to plotter --- include/openmc/plot.h | 2 ++ openmc/model/model.py | 3 ++- openmc/plots.py | 18 ++++++++++++++++++ src/plot.cpp | 2 ++ 4 files changed, 24 insertions(+), 1 deletion(-) diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 69e801c5fc3..8073a51e83e 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -171,6 +171,7 @@ class SlicePlotBase { Position origin_; //!< Plot origin in geometry Position width_; //!< Plot width in geometry PlotBasis basis_; //!< Plot basis (XY/XZ/YZ) + double time_; //!< Plot time array pixels_; //!< Plot size in pixels bool slice_color_overlaps_; //!< Show overlapping cells? int slice_level_ {-1}; //!< Plot universe level @@ -223,6 +224,7 @@ T SlicePlotBase::get_map() const GeometryState p; p.r() = xyz; p.u() = dir; + p.time() = time_; p.coord(0).universe = model::root_universe; int level = slice_level_; int j {}; diff --git a/openmc/model/model.py b/openmc/model/model.py index c19c8ac9f15..24b3f29cb8f 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -909,6 +909,7 @@ def plot( width: Sequence[float] | None = None, pixels: int | Sequence[int] = 40000, basis: str = 'xy', + time: float = 0.0, color_by: str = 'cell', colors: dict | None = None, seed: int | None = None, @@ -996,6 +997,7 @@ def plot( plot.width = width plot.pixels = pixels plot.basis = basis + plot.time = time plot.color_by = color_by plot.show_overlaps = show_overlaps if overlap_color is not None: @@ -1091,7 +1093,6 @@ def plot( if outline != 'only': axes.imshow(img, extent=(x_min, x_max, y_min, y_max), **kwargs) - if n_samples: # Sample external source particles particles = self.sample_external_source(n_samples) diff --git a/openmc/plots.py b/openmc/plots.py index 40fcf0c78da..405776c8206 100644 --- a/openmc/plots.py +++ b/openmc/plots.py @@ -185,6 +185,8 @@ the image aspect ratio. basis : {'xy', 'xz', 'yz'} The basis directions for the plot + time : float + Simulation time at which geometry is evaluated color_by : {'cell', 'material'} Indicate whether the plot should be colored by cell or by material colors : dict @@ -673,6 +675,8 @@ class Plot(PlotBase): The type of the plot basis : {'xy', 'xz', 'yz'} The basis directions for the plot + time : float + Time at which geometry is evaluated meshlines : dict Dictionary defining type, id, linewidth and color of a mesh to be plotted on top of a plot @@ -685,6 +689,7 @@ def __init__(self, plot_id=None, name=''): self._origin = [0., 0., 0.] self._type = 'slice' self._basis = 'xy' + self._time = 0.0 self._meshlines = None @property @@ -725,6 +730,15 @@ def basis(self, basis): cv.check_value('plot basis', basis, _BASES) self._basis = basis + @property + def time(self): + return self._time + + @time.setter + def time(self, time): + cv.check_type('plot time', time, Real) + self._time = time + @property def meshlines(self): return self._meshlines @@ -766,6 +780,7 @@ def __repr__(self): string += '{: <16}=\t{}\n'.format('\tBasis', self._basis) string += '{: <16}=\t{}\n'.format('\tWidth', self._width) string += '{: <16}=\t{}\n'.format('\tOrigin', self._origin) + string += '{: <16}=\t{}\n'.format('\tTime', self._time) string += '{: <16}=\t{}\n'.format('\tPixels', self._pixels) string += '{: <16}=\t{}\n'.format('\tColor by', self._color_by) string += '{: <16}=\t{}\n'.format('\tBackground', self._background) @@ -899,6 +914,9 @@ def to_xml_element(self): subelement = ET.SubElement(element, "width") subelement.text = ' '.join(map(str, self._width)) + subelement = ET.SubElement(element, "time") + subelement.text = str(self._time) + if self._colors: self._colors_to_xml(element) diff --git a/src/plot.cpp b/src/plot.cpp index dbc25b21e4d..a8a18a6fdbb 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -145,6 +145,7 @@ void Plot::print_info() const // Plot parameters fmt::print("Origin: {} {} {}\n", origin_[0], origin_[1], origin_[2]); + fmt::print("Time: {}\n", time_); if (PlotType::slice == type_) { fmt::print("Width: {:4} {:4}\n", width_[0], width_[1]); @@ -721,6 +722,7 @@ Plot::Plot(pugi::xml_node plot_node, PlotType type) { set_output_path(plot_node); set_basis(plot_node); + time_ = get_node_array(plot_node, "time")[0]; set_origin(plot_node); set_width(plot_node); set_meshlines(plot_node); From 4a7b114e2639b73bd79bf8c785c9ff30a1ea3181 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Mon, 9 Jun 2025 13:16:49 +0700 Subject: [PATCH 12/14] handle speed exception for MG --- src/particle.cpp | 14 +++++++++++++- src/simulation.cpp | 10 +++++++++- 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index 0f62e97b1fe..a0b489cbf4e 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -154,7 +154,14 @@ void Particle::event_calculate_xs() r_last() = r(); time_last() = time(); - speed() = get_speed(); + // Calculate speed. + // In MG mode, we can't determine particle speed if cell, and thus material, + // is not known. Just set to zero for now. + if (settings::run_CE || !(lowest_coord().cell == C_NONE)) { + this->speed() = this->get_speed(); + } else { + this->speed() = 0.0; + } // Reset event variables event() = TallyEvent::KILL; @@ -180,6 +187,11 @@ void Particle::event_calculate_xs() cell_last(j) = coord(j).cell; } n_coord_last() = n_coord(); + + // Calculate speed for MG mode + if (!settings::run_CE) { + this->speed() = 0.0; + } } // Write particle track. diff --git a/src/simulation.cpp b/src/simulation.cpp index a5de9f522ea..13f51c51a8f 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -577,7 +577,15 @@ void initialize_history(Particle& p, int64_t index_source) } p.current_work() = index_source; - p.speed() = p.get_speed(); + // Speed will be needed if the particle happens to be coincident on a moving + // surface. In MG mode, we get a chicken-egg situation as we need material to + // get speed, but we need speed to determine cell and thus material. + // Put zero for now. + if (settings::run_CE) { + p.speed() = p.get_speed(); + } else { + p.speed() = 0.0; + } // set identifier for particle p.id() = simulation::work_index[mpi::rank] + index_source; From e834f45f33d4b0c0599fa24599195a423466d719 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Wed, 11 Jun 2025 09:24:13 +0700 Subject: [PATCH 13/14] fix surface velocities xml parser --- src/surface.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/surface.cpp b/src/surface.cpp index 0b4b3ebe949..443778e1478 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -143,7 +143,7 @@ Surface::Surface(pugi::xml_node surf_node) for (int i = 0; i < velocities_spec.size();) { if (velocities_spec[i] == '-' || std::isdigit(velocities_spec[i])) { int j = i + 1; - while (j < velocities_spec.size() && std::isdigit(velocities_spec[j])) { + while (j < velocities_spec.size() && (std::isdigit(velocities_spec[j]) || velocities_spec[j] == '.')) { j++; } numbers.push_back(std::stod(velocities_spec.substr(i, j - i))); From 1b7823fa16f37d6baecd6a1d8b9558f93dfd2127 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Wed, 11 Jun 2025 13:02:06 +0700 Subject: [PATCH 14/14] correctly get speed for new particle --- src/particle.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/particle.cpp b/src/particle.cpp index a0b489cbf4e..5b0b4adbce7 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -190,7 +190,7 @@ void Particle::event_calculate_xs() // Calculate speed for MG mode if (!settings::run_CE) { - this->speed() = 0.0; + this->speed() = this->get_speed(); } }