diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 954ccf20..2cdc1229 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -6,6 +6,7 @@ set(Adamantine_HEADERS ${CMAKE_CURRENT_SOURCE_DIR}/DataAssimilator.hh ${CMAKE_CURRENT_SOURCE_DIR}/ElectronBeamHeatSource.hh ${CMAKE_CURRENT_SOURCE_DIR}/ExperimentalData.hh + ${CMAKE_CURRENT_SOURCE_DIR}/GaussianHeatSource.hh ${CMAKE_CURRENT_SOURCE_DIR}/Geometry.hh ${CMAKE_CURRENT_SOURCE_DIR}/GoldakHeatSource.hh ${CMAKE_CURRENT_SOURCE_DIR}/HeatSource.hh @@ -42,6 +43,7 @@ set(Adamantine_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/CubeHeatSource.cc ${CMAKE_CURRENT_SOURCE_DIR}/DataAssimilator.cc ${CMAKE_CURRENT_SOURCE_DIR}/ElectronBeamHeatSource.cc + ${CMAKE_CURRENT_SOURCE_DIR}/GaussianHeatSource.cc ${CMAKE_CURRENT_SOURCE_DIR}/Geometry.cc ${CMAKE_CURRENT_SOURCE_DIR}/GoldakHeatSource.cc ${CMAKE_CURRENT_SOURCE_DIR}/MaterialPropertyInstDev.cc diff --git a/source/GaussianHeatSource.cc b/source/GaussianHeatSource.cc new file mode 100644 index 00000000..f504f878 --- /dev/null +++ b/source/GaussianHeatSource.cc @@ -0,0 +1,211 @@ +/* SPDX-FileCopyrightText: Copyright (c) 2026, the adamantine authors. + * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception + */ + +#include +#include +#include + +#include + +namespace adamantine +{ + +template +GaussianHeatSource::GaussianHeatSource( + boost::property_tree::ptree const &beam_database, + boost::optional const + &units_optional_database) + : HeatSource(beam_database, units_optional_database), + _five_axis(this->_scan_path.is_five_axis()) +{ + _A = beam_database.get("A"); + _B = beam_database.get("B"); +} + +template +void GaussianHeatSource::update_time(double time) +{ + double const segment_power_modifier = + this->_scan_path.get_power_modifier(time); + this->_source_on = (segment_power_modifier > 0.0); + + dealii::Point<3> const &path = this->_scan_path.value(time); + _quaternion = this->_scan_path.get_current_quaternion(); + + for (unsigned int d = 0; d < 3; ++d) + { + _beam_center(d) = path(d); + } + + _depth = this->_beam.depth; + _radius_squared = this->_beam.radius_squared; + + double const aspect_ratio = std::max(1.0, _depth[0] / this->_beam.radius); + double const n = + std::min(std::max(0.0, _A * std::log2(aspect_ratio) + _B), 9.0); + _k = std::pow(2.0, n); + + double const V0 = + 0.5 * dealii::numbers::PI * _radius_squared[0] * _depth[0] * + std::tgamma(1.0 / _k) / (_k * std::pow(3.0, 1.0 / _k)); + + double const effective_power = + this->_beam.absorption_efficiency * this->_beam.max_power * + segment_power_modifier; + + _alpha = effective_power / V0; +} + +template +double GaussianHeatSource::value(dealii::Point const &point) const +{ + dealii::Point rotated_point; + if constexpr (dim == 2) + { + rotated_point = point; + } + else + { + rotated_point = _five_axis ? _quaternion.rotate(point) : point; + } + + double const z = rotated_point[axis::z] - _beam_center[axis::z][0]; + + if (z > 0. || z < -_depth[0]) + { + return 0.; + } + + double xpy_squared = dealii::Utilities::fixed_power<2>( + rotated_point[axis::x] - _beam_center[axis::x][0]); + if constexpr (dim == 3) + { + xpy_squared += dealii::Utilities::fixed_power<2>( + rotated_point[axis::y] - _beam_center[axis::y][0]); + } + + if (xpy_squared > 5. * _radius_squared[0]) + { + return 0.; + } + + double const radial_component = + std::exp(-2.0 * xpy_squared / _radius_squared[0]); + + double const depth_component = + std::exp(-3.0 * std::pow(std::abs(z / _depth[0]), _k)); + + return _alpha[0] * radial_component * depth_component; +} + +template +dealii::VectorizedArray GaussianHeatSource::value( + dealii::Point> const &points) const +{ + dealii::Point> rotated_points; + if constexpr (dim == 2) + { + rotated_points = points; + } + else + { + rotated_points = _five_axis ? _quaternion.rotate(points) : points; + } + + auto const z = rotated_points[axis::z] - _beam_center[axis::z]; + + dealii::VectorizedArray depth_mask; + for (unsigned int i = 0; i < depth_mask.size(); ++i) + { + depth_mask[i] = (z[i] > 0. || z[i] < -_depth[i]) ? 0. : 1.; + } + if (depth_mask.sum() == 0) + { + return depth_mask; + } + + auto xpy_squared = dealii::Utilities::fixed_power<2>( + rotated_points[axis::x] - _beam_center[axis::x]); + if constexpr (dim == 3) + { + auto const y_squared = dealii::Utilities::fixed_power<2>( + rotated_points[axis::y] - _beam_center[axis::y]); + xpy_squared += y_squared; + } + + auto const xpy_area = 5. * _radius_squared - xpy_squared; + dealii::VectorizedArray xpy_mask; + for (unsigned int i = 0; i < xpy_mask.size(); ++i) + { + xpy_mask[i] = xpy_area[i] < 0. ? 0. : 1.; + } + if (xpy_mask.sum() == 0) + { + return xpy_mask; + } + + dealii::VectorizedArray shape = _k; + + auto const radial_component = + std::exp(-2.0 * xpy_squared / _radius_squared); + + auto const depth_component = + std::exp(-3.0 * std::pow(std::abs(z / _depth), shape)); + + return depth_mask * _alpha * radial_component * depth_component; +} + +template +dealii::BoundingBox +GaussianHeatSource::get_bounding_box(double const time, + double const scaling_factor) const +{ + dealii::Point<3> const &beam_center = this->_scan_path.value(time, false); + if constexpr (dim == 2) + { + return {{{beam_center[axis::x] - scaling_factor * this->_beam.radius, + beam_center[axis::z] - scaling_factor * this->_beam.depth}, + {beam_center[axis::x] + scaling_factor * this->_beam.radius, + beam_center[axis::z]}}}; + } + else + { + dealii::Point<3> max_corner( + beam_center[axis::x] + scaling_factor * this->_beam.radius, + beam_center[axis::y] + scaling_factor * this->_beam.radius, + beam_center[axis::z]); + dealii::Point<3> min_corner( + beam_center[axis::x] - scaling_factor * this->_beam.radius, + beam_center[axis::y] - scaling_factor * this->_beam.radius, + beam_center[axis::z] - scaling_factor * this->_beam.depth); + + if (this->_scan_path.is_five_axis()) + { + dealii::Point<3> rotated_max_corner = + this->_scan_path.rotate(time, max_corner); + dealii::Point<3> rotated_min_corner = + this->_scan_path.rotate(time, min_corner); + + dealii::Point<3> new_max_corner; + dealii::Point<3> new_min_corner; + for (int d = 0; d < 3; ++d) + { + new_max_corner[d] = + std::max(rotated_max_corner[d], rotated_min_corner[d]); + new_min_corner[d] = + std::min(rotated_max_corner[d], rotated_min_corner[d]); + } + + return {{new_min_corner, new_max_corner}}; + } + else + { + return {{min_corner, max_corner}}; + } + } +} + +} // namespace adamantine + +INSTANTIATE_DIM(GaussianHeatSource) diff --git a/source/GaussianHeatSource.hh b/source/GaussianHeatSource.hh new file mode 100644 index 00000000..6f944e1e --- /dev/null +++ b/source/GaussianHeatSource.hh @@ -0,0 +1,87 @@ +/* SPDX-FileCopyrightText: Copyright (c) 2026, the adamantine authors. + * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception + */ + +#ifndef GAUSSIAN_HEAT_SOURCE_HH +#define GAUSSIAN_HEAT_SOURCE_HH + +#include +#include + +#include + +namespace adamantine +{ +/** + * A derived class from HeatSource for the Gaussian model of a laser heat source. + */ +template +class GaussianHeatSource final : public HeatSource +{ +public: + /** + * Constructor. + * \param[in] beam_database requires the following entries: + * - absorption_efficiency: double in \f$[0,1]\f$ + * - depth: double in \f$[0,\infty)\f$ + * - diameter: double in \f$[0,\infty)\f$ + * - max_power: double in \f$[0, \infty)\f$ + * - scan_path_file: name of the file that contains the scan path + * segments + * - scan_path_file_format: format of the scan path file + * - A: double, empirical coefficient for the model + * - B: double, empirical coefficient for the model + * \param[in] units_optional_database may contain the following entries: + * - heat_source.dimension + * - heat_source.power + */ + GaussianHeatSource( + boost::property_tree::ptree const &beam_database, + boost::optional const + &units_optional_database); + + /** + * Set the time variable. This updates the beam position and recalculates + * the shape factor (_k) and normalization constant (_alpha). + */ + void update_time(double time) final; + + /** + * Returns the value of a Gaussian heat source at a specified point. + */ + double value(dealii::Point const &point) const final; + + /** + * Same function as above but it uses vectorized data. + */ + dealii::VectorizedArray + value(dealii::Point> const &points) + const final; + + /** + * Returns the bounding box of the heat source at a given time. + */ + dealii::BoundingBox + get_bounding_box(double time, double const scaling_factor) const final; + +private: + double _A; + double _B; + + bool const _five_axis; + Quaternion _quaternion; + + dealii::Point<3, dealii::VectorizedArray> _beam_center; + dealii::VectorizedArray _depth = + std::numeric_limits::signaling_NaN(); + dealii::VectorizedArray _radius_squared = + std::numeric_limits::signaling_NaN(); + + dealii::VectorizedArray _alpha = + std::numeric_limits::signaling_NaN(); + + double _k = std::numeric_limits::signaling_NaN(); +}; +} // namespace adamantine + +#endif diff --git a/source/ThermalPhysics.templates.hh b/source/ThermalPhysics.templates.hh index ca9693e7..ead50412 100644 --- a/source/ThermalPhysics.templates.hh +++ b/source/ThermalPhysics.templates.hh @@ -5,6 +5,7 @@ #ifndef THERMAL_PHYSICS_TEMPLATES_HH #define THERMAL_PHYSICS_TEMPLATES_HH +#include #include #include #include @@ -402,6 +403,11 @@ ThermalPhysics>( beam_database, units_optional_database); } + else if (type == "gaussian") + { + _heat_sources[i] = std::make_shared>( + beam_database, units_optional_database); + } else { ASSERT_THROW(false, "Beam type '" + diff --git a/source/ensemble_management.cc b/source/ensemble_management.cc index d1fe6cdd..c0a27f6d 100644 --- a/source/ensemble_management.cc +++ b/source/ensemble_management.cc @@ -2,6 +2,7 @@ * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception */ +#include #include #include #include @@ -158,6 +159,16 @@ std::vector>> get_bounding_heat_sources( bounding_heat_sources[i] = std::make_shared>( beam_database, units_optional_database); } + else if (type == "gaussian") + { + boost::property_tree::ptree modified_beam_database = beam_database; + modified_beam_database.put("diameter", diameter_max[bounding_source]); + modified_beam_database.put("depth", depth_max[bounding_source]); + ++bounding_source; + + bounding_heat_sources[i] = std::make_shared>( + modified_beam_database, units_optional_database); + } } return bounding_heat_sources; diff --git a/source/validate_input_database.cc b/source/validate_input_database.cc index 2467fd33..ffe244a9 100644 --- a/source/validate_input_database.cc +++ b/source/validate_input_database.cc @@ -369,10 +369,11 @@ void validate_input_database(boost::property_tree::ptree &database) "sources.beam_" + std::to_string(beam_index) + ".type"); ASSERT_THROW(boost::iequals(beam_type, "goldak") || boost::iequals(beam_type, "electron_beam") || - boost::iequals(beam_type, "cube"), + boost::iequals(beam_type, "cube") || + boost::iequals(beam_type, "gaussian"), "Beam type, '" + beam_type + "', is not recognized. Valid options are: 'goldak', " - "'electron_beam', and 'cube'."); + "'electron_beam', 'cube', and 'gaussian'."); ASSERT_THROW(database.get_child("sources") .get_child("beam_" + std::to_string(beam_index)) .count("scan_path_file") != 0,