diff --git a/src/Finch_Inputs.hpp b/src/Finch_Inputs.hpp index d740dbd..243966b 100644 --- a/src/Finch_Inputs.hpp +++ b/src/Finch_Inputs.hpp @@ -86,6 +86,15 @@ struct Source std::array two_sigma; std::array r; std::string scan_path_file; + // Heat source shape model: + // "gaussian" (default) — original isotropic 3D Gaussian, backwards + // compatible "dynamic" — super-Gaussian (Coleman et al. 2024), + // requires k and m + std::string shape = "gaussian"; + double k = + 2.0; // radial distribution parameter (k=2: Gaussian, k->inf: uniform) + double m = 2.0; // volumetric shape parameter (m=1: cone, m=2: ellipsoid, + // m->inf: cylinder) }; struct Properties @@ -230,6 +239,12 @@ class Inputs Info << " Y: " << source.two_sigma[1] << std::endl; Info << " Z: " << source.two_sigma[2] << std::endl; Info << " scan path file: " << source.scan_path_file << std::endl; + Info << " shape: " << source.shape << std::endl; + if ( source.shape == "dynamic" ) + { + Info << " k (radial distribution): " << source.k << std::endl; + Info << " m (volumetric shape): " << source.m << std::endl; + } // Print solidification output options Info << "Sampling:" << std::endl; @@ -462,6 +477,26 @@ class Inputs source.two_sigma[2] = fabs( source.two_sigma[2] ); source.scan_path_file = db["source"]["scan_path_file"]; + + // Heat Source Shape Model — optional, defaults to "gaussian" for + // backwards compatibility with existing input files + source.shape = db["source"].value( "shape", "gaussian" ); + + if ( source.shape == "dynamic" ) + { + if ( !db["source"].contains( "k" ) || + !db["source"].contains( "m" ) ) + throw std::runtime_error( + "Error: source.shape = \"dynamic\" requires both " + "\"k\" and \"m\" to be specified in the input file." ); + source.k = db["source"]["k"]; + source.m = db["source"]["m"]; + } + else if ( source.shape != "gaussian" ) + { + throw std::runtime_error( + "Error: source.shape must be \"gaussian\" or \"dynamic\"." ); + } } void readInputSampling( nlohmann::json db ) diff --git a/src/Finch_Solver.hpp b/src/Finch_Solver.hpp index 5653edb..8745705 100644 --- a/src/Finch_Solver.hpp +++ b/src/Finch_Solver.hpp @@ -51,11 +51,23 @@ class Solver // heat source parameters double power_; double position_[3]; + + // Original Gaussian heat source parameters double r_[3]; double A_inv_[3]; + + // Shared double I0_; double w_max_; + // Dynamic super-Gaussian parameters (Coleman et al. 2024) + // Only used when source.shape == "dynamic" + bool dynamic_beam_; + double sigma_[3]; // half-widths for super-Gaussian radial profile + double k_; // radial distribution parameter + double m_; // volumetric shape parameter + double depth_; // heat source depth (= two_sigma[2]) + public: Solver( Inputs db, LocalMeshType local_mesh ) : local_mesh_( local_mesh ) @@ -85,15 +97,52 @@ class Solver position_[d] = 0.0; } - // heat source parameter constants - for ( std::size_t d = 0; d < 3; ++d ) + dynamic_beam_ = ( db.source.shape == "dynamic" ); + + if ( !dynamic_beam_ ) { - r_[d] = db.source.two_sigma[d] / Kokkos::sqrt( 2.0 ); - A_inv_[d] = 1.0 / r_[d] / r_[d]; + // Original isotropic 3D Gaussian + for ( std::size_t d = 0; d < 3; ++d ) + { + r_[d] = db.source.two_sigma[d] / Kokkos::sqrt( 2.0 ); + A_inv_[d] = 1.0 / r_[d] / r_[d]; + } + + I0_ = ( 2.0 * db.source.absorption ) / + ( M_PI * Kokkos::sqrt( M_PI ) * r_[0] * r_[1] * r_[2] ); } + else + { + // Dynamic super-Gaussian (Coleman et al. 2024, Eq. 2-6) + k_ = db.source.k; + m_ = db.source.m; + + // sigma_i = two_sigma_i / (2 * 2^(1/k)) + double two_to_1_over_k = Kokkos::pow( 2.0, 1.0 / k_ ); + for ( std::size_t d = 0; d < 3; ++d ) + { + sigma_[d] = db.source.two_sigma[d] / ( 2.0 * two_to_1_over_k ); + } + + // Heat source depth = two_sigma[2] + depth_ = db.source.two_sigma[2]; + + // Areal normalization A0 (Eq. 3) + double two_to_2_over_k = Kokkos::pow( 2.0, 2.0 / k_ ); + double rx0 = db.source.two_sigma[0] / two_to_2_over_k; + double ry0 = db.source.two_sigma[1] / two_to_2_over_k; + double A0 = M_PI * tgamma( 1.0 + 2.0 / k_ ) * rx0 * ry0; - I0_ = ( 2.0 * db.source.absorption ) / - ( M_PI * Kokkos::sqrt( M_PI ) * r_[0] * r_[1] * r_[2] ); + // Depth integral factor from Eq. 6 + double depth_factor = + ( tgamma( 1.0 + 1.0 / m_ ) * tgamma( 1.0 + 2.0 / m_ ) ) / + tgamma( 1.0 + 3.0 / m_ ); + + // Volume normalization V0 = A0 * depth * depth_factor (Eq. 6) + double V0 = A0 * depth_ * depth_factor; + + I0_ = db.source.absorption / V0; + } // cut off for 3 standard deviations from heat source center w_max_ = Kokkos::log( 3 ) + 2 * Kokkos::log( 10 ); @@ -174,23 +223,48 @@ class Solver k_by_dx2_; } - // Normalized weight for the gaussian source term: x in exp(-x) + // Normalized weight for the source term: returned value w used as exp(-w) + // Branches on dynamic_beam_ set at construction time. KOKKOS_INLINE_FUNCTION auto weight( const int i, const int j, const int k ) const { double grid_loc[3]; - double dist_to_beam[3]; int idx[3] = { i, j, k }; - local_mesh_.coordinates( EntityType(), idx, grid_loc ); - dist_to_beam[0] = grid_loc[0] - position_[0]; - dist_to_beam[1] = grid_loc[1] - position_[1]; - dist_to_beam[2] = grid_loc[2] - position_[2]; + double dx = grid_loc[0] - position_[0]; + double dy = grid_loc[1] - position_[1]; + double dz = grid_loc[2] - position_[2]; - return ( dist_to_beam[0] * dist_to_beam[0] * A_inv_[0] ) + - ( dist_to_beam[1] * dist_to_beam[1] * A_inv_[1] ) + - ( dist_to_beam[2] * dist_to_beam[2] * A_inv_[2] ); + if ( !dynamic_beam_ ) + { + // Original isotropic 3D Gaussian + return ( dx * dx * A_inv_[0] ) + ( dy * dy * A_inv_[1] ) + + ( dz * dz * A_inv_[2] ); + } + else + { + // Dynamic super-Gaussian (Coleman et al. 2024, Eq. 4-5) + // Normalized depth: only apply heat downward from beam center + double dz_norm = Kokkos::fabs( dz ) / depth_; + if ( dz_norm >= 1.0 ) + return w_max_ + 1.0; // below heat source depth: no contribution + + // Radius decays with depth: r(z) = sigma * (1 - |dz/d|^m)^(1/m) + double shape_factor = + Kokkos::pow( 1.0 - Kokkos::pow( dz_norm, m_ ), 1.0 / m_ ); + double rx = sigma_[0] * shape_factor; + double ry = sigma_[1] * shape_factor; + + // Guard against division by zero near the tip + if ( rx < 1e-15 || ry < 1e-15 ) + return w_max_ + 1.0; + + // Super-Gaussian radial exponent: (dx^2/rx^2 + dy^2/ry^2)^(k/2) + double radial = + ( dx * dx ) / ( rx * rx ) + ( dy * dy ) / ( ry * ry ); + return Kokkos::pow( radial, k_ / 2.0 ); + } } // Heating source term, device overload.