From e985449cfe569165a70b3b066ca9bf377d80aba8 Mon Sep 17 00:00:00 2001 From: Samuel Roberts Date: Tue, 4 Aug 2026 23:08:04 -0500 Subject: [PATCH 1/7] Add a Scalar template parameter to the heat solve for forward-mode AD Finch's field type was fixed to double. This threads a Scalar template parameter (defaulting to double) through Grid, Boundary, Solver, Layer and SolidificationData, so the temperature field can carry a user-supplied arithmetic type. With double the generated code and temperature fields are unchanged; verified bit-identical temperature output and identical solidification records on the single_line case. The differentiated material and source inputs move into a MaterialProperties struct so a caller can supply values of the field scalar type rather than plain doubles read from the input deck. Quantities that only select a branch (solidus, liquidus) and all mesh geometry stay double. Finch_Scalar.hpp adds exp/fmin/fmax dispatch: arithmetic types forward to Kokkos as before, anything else resolves by ADL into the scalar type's own namespace. Finch therefore carries no dependency on any AD library. The optional finch_sensitivity application demonstrates the result with cpp_oti_lib, obtaining d(QoI)/dp for six material and source parameters from a single solve and validating each against central finite differences. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01PNkQTS2g7P6dNDybhWRkwz --- applications/CMakeLists.txt | 20 ++ applications/Finch_OTI.hpp | 161 ++++++++++++++++ applications/Sensitivity.cpp | 305 +++++++++++++++++++++++++++++++ src/Finch_Boundary.hpp | 11 +- src/Finch_Grid.hpp | 60 +++++- src/Finch_Run.hpp | 13 +- src/Finch_Scalar.hpp | 110 +++++++++++ src/Finch_SolidificationData.hpp | 51 +++--- src/Finch_Solver.hpp | 134 ++++++++++---- 9 files changed, 793 insertions(+), 72 deletions(-) create mode 100644 applications/Finch_OTI.hpp create mode 100644 applications/Sensitivity.cpp create mode 100644 src/Finch_Scalar.hpp diff --git a/applications/CMakeLists.txt b/applications/CMakeLists.txt index 839f979..5610efb 100644 --- a/applications/CMakeLists.txt +++ b/applications/CMakeLists.txt @@ -1,3 +1,23 @@ add_executable(finch SingleLayer.cpp) target_link_libraries(finch Core) install(TARGETS finch DESTINATION ${CMAKE_INSTALL_BINDIR}) + +# Optional forward-mode AD sensitivity application. Built only when the +# header-only cpp_oti_lib is available; Finch itself never depends on it. +find_path(CPP_OTI_LIB_INCLUDE_DIR otinum/otinum.hpp + HINTS + ${CPP_OTI_LIB_DIR} + ${CMAKE_CURRENT_SOURCE_DIR}/../../cpp_oti_lib/include + DOC "Path to the cpp_oti_lib include directory") + +if(CPP_OTI_LIB_INCLUDE_DIR) + message(STATUS "Finch sensitivity: cpp_oti_lib found at ${CPP_OTI_LIB_INCLUDE_DIR}") + add_executable(finch_sensitivity Sensitivity.cpp) + target_link_libraries(finch_sensitivity Core) + target_include_directories(finch_sensitivity PRIVATE + ${CPP_OTI_LIB_INCLUDE_DIR} + ${CMAKE_CURRENT_SOURCE_DIR}) + install(TARGETS finch_sensitivity DESTINATION ${CMAKE_INSTALL_BINDIR}) +else() + message(STATUS "Finch sensitivity: cpp_oti_lib not found, skipping finch_sensitivity") +endif() diff --git a/applications/Finch_OTI.hpp b/applications/Finch_OTI.hpp new file mode 100644 index 0000000..1457e24 --- /dev/null +++ b/applications/Finch_OTI.hpp @@ -0,0 +1,161 @@ +/**************************************************************************** + * OTI (order-truncated imaginary) scalar type glue for Finch. + * + * This header is the ONLY place where Finch and cpp_oti_lib meet. It lives in + * applications/ rather than src/ on purpose: Finch's core carries no dependency + * on any AD library, and the coupling is a single ScalarValue specialization + * plus a parameter-seeding helper. + ****************************************************************************/ + +#ifndef Finch_OTI_H +#define Finch_OTI_H + +#include +#include + +#include "otinum/otinum.hpp" + +#include +#include +#include + +namespace Finch +{ +namespace Math +{ + +// Tell Finch which component of an OTI jet is its numeric value. This is all +// the core library needs to know about the type; every other operation it +// performs (arithmetic, comparison, exp, fmin/fmax) is resolved by ADL into +// namespace oti. +template +struct ScalarValue> +{ + KOKKOS_INLINE_FUNCTION static double + value( const oti::otinum& x ) + { + return static_cast( x.real() ); + } +}; + +} // namespace Math + +namespace Sensitivity +{ + +// The differentiated parameter set for this study. +// +// two_sigma is treated as a single parameter applied to all three source +// directions, which is what the input deck describes for an axisymmetric spot. +// Splitting it into three independent directions is a matter of adding two +// more slots below and seeding them separately. +enum Parameter +{ + Density = 0, + SpecificHeat, + ThermalConductivity, + LatentHeat, + Absorption, + TwoSigma, + NumParameters +}; + +inline const char* name( int p ) +{ + switch ( p ) + { + case Density: + return "density"; + case SpecificHeat: + return "specific_heat"; + case ThermalConductivity: + return "thermal_conductivity"; + case LatentHeat: + return "latent_heat"; + case Absorption: + return "absorption"; + case TwoSigma: + return "two_sigma"; + default: + return "unknown"; + } +} + +inline const char* units( int p ) +{ + switch ( p ) + { + case Density: + return "kg/m^3"; + case SpecificHeat: + return "J/kg/K"; + case ThermalConductivity: + return "W/m/K"; + case LatentHeat: + return "J/kg"; + case Absorption: + return "-"; + case TwoSigma: + return "m"; + default: + return ""; + } +} + +// Nominal value of each parameter, read from the input deck. +inline std::array nominal( const Inputs& db ) +{ + // The single two_sigma slot presumes an axisymmetric source. Fail loudly + // rather than silently differentiate only one direction. + if ( db.source.two_sigma[0] != db.source.two_sigma[1] || + db.source.two_sigma[0] != db.source.two_sigma[2] ) + throw std::runtime_error( + "Sensitivity: two_sigma is seeded as a single parameter, which " + "requires the three components to be equal in the input deck." ); + + std::array p; + p[Density] = db.properties.density; + p[SpecificHeat] = db.properties.specific_heat; + p[ThermalConductivity] = db.properties.thermal_conductivity; + p[LatentHeat] = db.properties.latent_heat; + p[Absorption] = db.source.absorption; + p[TwoSigma] = db.source.two_sigma[0]; + return p; +} + +// Assemble solver properties from an explicit parameter vector. Used for both +// the plain-double finite-difference runs (values perturbed) and the OTI run +// (values seeded as independent variables), so that both paths perturb exactly +// the same quantities. +template +MaterialProperties build( const std::array& p ) +{ + MaterialProperties props; + props.density = p[Density]; + props.specific_heat = p[SpecificHeat]; + props.thermal_conductivity = p[ThermalConductivity]; + props.latent_heat = p[LatentHeat]; + props.absorption = p[Absorption]; + for ( int d = 0; d < 3; ++d ) + props.two_sigma[d] = p[TwoSigma]; + return props; +} + +// Seed every parameter as an independent OTI variable at its nominal value. +template +std::array seed( const Inputs& db ) +{ + static_assert( OTI::nvars == NumParameters, + "OTI algebra must have one variable per differentiated " + "parameter" ); + auto nom = nominal( db ); + std::array p; + for ( int i = 0; i < NumParameters; ++i ) + p[i] = OTI::variable( i, nom[i] ); + return p; +} + +} // namespace Sensitivity +} // namespace Finch + +#endif diff --git a/applications/Sensitivity.cpp b/applications/Sensitivity.cpp new file mode 100644 index 0000000..31309cd --- /dev/null +++ b/applications/Sensitivity.cpp @@ -0,0 +1,305 @@ +/**************************************************************************** + * Parameter sensitivity of the Finch heat solve, by forward-mode AD. + * + * Runs the single-layer heat transport problem once with OTI-valued material + * and source parameters, obtaining d(QoI)/dp for every parameter from that one + * solve, then validates each derivative against a central finite difference + * built from two ordinary double solves per parameter. + ****************************************************************************/ + +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include "Finch_Core.hpp" +#include "Finch_OTI.hpp" + +namespace +{ + +constexpr int NP = Finch::Sensitivity::NumParameters; + +// First-order jets: one variable per differentiated parameter, derivatives +// through order one. Raising the order to 2 here is the only change needed to +// obtain the full parameter Hessian of the same quantities. +using OTI = oti::otinum; + +// Quantities of interest evaluated on the final temperature field. Both are +// smooth functionals of the field, which makes them meaningful finite- +// difference references. +template +struct QoI +{ + // Sum of temperature over all owned nodes. + Scalar sum; + // Temperature at a single fixed node (centre of the owned index space). + Scalar probe; +}; + +// Run the transient solve to completion and evaluate the QoIs. +template +QoI solve( MPI_Comm comm, Finch::Inputs db, + const Finch::MaterialProperties& props ) +{ + std::array bc_types = { "adiabatic", "adiabatic", + "adiabatic", "adiabatic", + "adiabatic", "adiabatic" }; + + Finch::Grid grid( + comm, db.space.cell_size, db.space.global_low_corner, + db.space.global_high_corner, db.space.ranks_per_dim, bc_types, + db.space.initial_temperature ); + + auto fd = Finch::createSolver( db, grid, props ); + + Finch::MovingBeam beam( db.source.scan_path_file ); + + // Time loop. This mirrors Finch::Layer::step without the solidification + // sampling and file output, which are not needed for the QoIs here. + double time = db.time.start_time; + const double dt = db.time.time_step; + + for ( int n = 0; n < db.time.num_steps; ++n ) + { + time += dt; + + beam.move( time ); + double beam_power = beam.power(); + double beam_pos[3]; + for ( std::size_t d = 0; d < 3; ++d ) + beam_pos[d] = beam.position( d ); + + auto T = grid.getTemperature(); + auto T0 = grid.getPreviousTemperature(); + + Kokkos::deep_copy( T0, T ); + + auto owned_space = grid.getIndexSpace(); + fd.solve( ExecSpace(), owned_space, T, T0, beam_power, beam_pos ); + + grid.updateBoundaries(); + grid.gather(); + } + + // Evaluate the QoIs on the host. Copying the field to a host mirror keeps + // this independent of Kokkos reducer support for compound scalar types, + // which is not needed for a once-per-run diagnostic. + auto T = grid.getTemperature(); + auto T_host = Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), T ); + + auto owned = grid.getIndexSpace(); + + QoI q; + q.sum = Scalar( 0 ); + for ( int i = owned.min( 0 ); i < owned.max( 0 ); ++i ) + for ( int j = owned.min( 1 ); j < owned.max( 1 ); ++j ) + for ( int k = owned.min( 2 ); k < owned.max( 2 ); ++k ) + q.sum = q.sum + T_host( i, j, k, 0 ); + + q.probe = T_host( ( owned.min( 0 ) + owned.max( 0 ) ) / 2, + ( owned.min( 1 ) + owned.max( 1 ) ) / 2, + ( owned.min( 2 ) + owned.max( 2 ) ) / 2, 0 ); + + // Reduce the sum across ranks. An OTI jet is a contiguous block of + // coefficients, and summing jets is summing coefficients elementwise, so a + // plain MPI_DOUBLE reduction over ncoeffs entries is exact -- no derived + // datatype or user-defined operator is required. + int comm_size; + MPI_Comm_size( comm, &comm_size ); + if ( comm_size > 1 ) + { + if constexpr ( std::is_same::value ) + { + double local = q.sum; + MPI_Allreduce( &local, &q.sum, 1, MPI_DOUBLE, MPI_SUM, comm ); + } + else + { + Scalar local = q.sum; + MPI_Allreduce( &local[0], &q.sum[0], Scalar::ncoeffs, MPI_DOUBLE, + MPI_SUM, comm ); + } + } + + return q; +} + +void run( MPI_Comm comm, int argc, char* argv[] ) +{ + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = exec_space::memory_space; + + int rank; + MPI_Comm_rank( comm, &rank ); + + Finch::Inputs db( comm, argc, argv ); + + auto nominal = Finch::Sensitivity::nominal( db ); + + // FD step is overridable so the derivative can be checked for step-size + // independence -- the standard way to tell a genuine discrepancy from a + // finite-difference artifact. + double rel_step = 1e-6; + if ( const char* s = std::getenv( "FINCH_FD_STEP" ) ) + rel_step = std::atof( s ); + + // Warm up: the first solve of a process pays for page faults and OpenMP + // thread start-up, which would otherwise be charged to whichever solve + // happens to run first. + { + std::array warm; + for ( int i = 0; i < NP; ++i ) + warm[i] = nominal[i]; + Finch::Inputs warm_db = db; + warm_db.time.num_steps = 2; + solve( + comm, warm_db, Finch::Sensitivity::build( warm ) ); + } + + // ---- One OTI solve gives every derivative ----------------------------- + double t0 = MPI_Wtime(); + auto seeded = Finch::Sensitivity::seed( db ); + auto oti_result = solve( + comm, db, Finch::Sensitivity::build( seeded ) ); + double t_oti = MPI_Wtime() - t0; + + // ---- Reference: one plain double solve, then two per parameter -------- + t0 = MPI_Wtime(); + std::array base; + for ( int i = 0; i < NP; ++i ) + base[i] = nominal[i]; + auto ref = solve( + comm, db, Finch::Sensitivity::build( base ) ); + double t_double = MPI_Wtime() - t0; + + std::array fd_sum, fd_probe; + t0 = MPI_Wtime(); + for ( int i = 0; i < NP; ++i ) + { + double h = rel_step * std::fabs( nominal[i] ); + + auto plus = base; + plus[i] = nominal[i] + h; + auto qp = solve( + comm, db, Finch::Sensitivity::build( plus ) ); + + auto minus = base; + minus[i] = nominal[i] - h; + auto qm = solve( + comm, db, Finch::Sensitivity::build( minus ) ); + + fd_sum[i] = ( qp.sum - qm.sum ) / ( 2.0 * h ); + fd_probe[i] = ( qp.probe - qm.probe ) / ( 2.0 * h ); + } + double t_fd = MPI_Wtime() - t0; + + if ( rank != 0 ) + return; + + // ---- Report ----------------------------------------------------------- + std::printf( "\n" ); + std::printf( "=================================================" + "=================================\n" ); + std::printf( " Finch parameter sensitivity: OTI forward-mode AD vs " + "central finite differences\n" ); + std::printf( "=================================================" + "=================================\n" ); + std::printf( " algebra : oti::otinum<%d,%d> (%d coefficients " + "per node)\n", + OTI::nvars, OTI::order, OTI::ncoeffs ); + std::printf( " time steps : %d\n", db.time.num_steps ); + std::printf( " FD relative step : %g\n", rel_step ); + std::printf( "\n" ); + std::printf( " value check T_sum OTI %.10e double %.10e\n", + oti_result.sum.real(), ref.sum ); + std::printf( " value check T_probe OTI %.10e double %.10e\n", + oti_result.probe.real(), ref.probe ); + std::printf( "\n" ); + + const char* qoi_name[2] = { "T_sum [K]", "T_probe[K]" }; + + for ( int q = 0; q < 2; ++q ) + { + const OTI& jet = ( q == 0 ) ? oti_result.sum : oti_result.probe; + const std::array& fd = ( q == 0 ) ? fd_sum : fd_probe; + + std::printf( " d(%s)/dp\n", qoi_name[q] ); + std::printf( " %-22s %16s %16s %11s %14s\n", "parameter", "OTI", + "central FD", "rel.diff", "p*dQ/dp [K]" ); + std::printf( " %-22s %16s %16s %11s %14s\n", "----------------------", + "----------------", "----------------", "-----------", + "--------------" ); + + // A relative difference is only meaningful if the derivative is itself + // non-negligible. Compare each row against the largest normalized + // sensitivity in the table so that a structurally-zero derivative is + // reported as such rather than as a 100% disagreement. + double ref = 0.0; + for ( int i = 0; i < NP; ++i ) + { + typename OTI::alpha_type alpha{}; + alpha[i] = 1; + ref = + std::max( ref, std::fabs( nominal[i] * jet.partial( alpha ) ) ); + } + + for ( int i = 0; i < NP; ++i ) + { + typename OTI::alpha_type alpha{}; + alpha[i] = 1; + double d_oti = jet.partial( alpha ); + + double scale = std::max( std::fabs( d_oti ), std::fabs( fd[i] ) ); + bool negligible = nominal[i] != 0.0 && + std::fabs( nominal[i] ) * scale < 1e-8 * ref; + + char label[64]; + std::snprintf( label, sizeof( label ), "%s [%s]", + Finch::Sensitivity::name( i ), + Finch::Sensitivity::units( i ) ); + + std::printf( " %-22s %16.6e %16.6e ", label, d_oti, fd[i] ); + if ( negligible ) + std::printf( "%11s ", "negligible" ); + else if ( scale > 0.0 ) + std::printf( "%11.2e ", std::fabs( d_oti - fd[i] ) / scale ); + else + std::printf( "%11s ", "-" ); + std::printf( "%14.4e\n", nominal[i] * d_oti ); + } + std::printf( "\n" ); + } + + std::printf( " cost: 1 OTI solve %.3f s | 1 double solve %.3f s" + " | %d FD solves %.3f s\n", + t_oti, t_double, 2 * NP, t_fd ); + std::printf( " OTI overhead vs one double solve: %.2fx " + "(all %d derivatives)\n", + t_oti / t_double, NP ); + std::printf( " FD cost for the same %d derivatives: %.2fx\n", NP, + t_fd / t_double ); + std::printf( "=================================================" + "=================================\n\n" ); +} + +} // namespace + +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + run( MPI_COMM_WORLD, argc, argv ); + + Kokkos::finalize(); + MPI_Finalize(); + + return 0; +} diff --git a/src/Finch_Boundary.hpp b/src/Finch_Boundary.hpp index 4b6bb93..17b4c3f 100644 --- a/src/Finch_Boundary.hpp +++ b/src/Finch_Boundary.hpp @@ -18,12 +18,16 @@ namespace Finch { +// Boundary values are field values, so they carry the field scalar type. That +// makes a Dirichlet or Neumann value differentiable along with everything else +// when Scalar is an AD type. +template class Boundary { public: // Constructor with BC types and values Boundary( std::array types, - Kokkos::Array values ) + Kokkos::Array values ) : boundary_types( types ) , boundary_values( values ) { @@ -33,7 +37,8 @@ class Boundary // Constructor where no values are needed. Boundary( std::array types ) : boundary_types( types ) - , boundary_values( { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } ) + , boundary_values( { Scalar( 0 ), Scalar( 0 ), Scalar( 0 ), Scalar( 0 ), + Scalar( 0 ), Scalar( 0 ) } ) { for ( int d = 0; d < 6; d++ ) if ( boundary_types[d] == "dirichlet" || @@ -111,7 +116,7 @@ class Boundary //! Boundary types for each plane. std::array boundary_types; //! Boundary values for each plane. - Kokkos::Array boundary_values; + Kokkos::Array boundary_values; //! Boundary types for each plane, converted to int for device. Kokkos::Array boundary_int; //! Boundary indices for each plane. diff --git a/src/Finch_Grid.hpp b/src/Finch_Grid.hpp index b22af68..c3d5148 100644 --- a/src/Finch_Grid.hpp +++ b/src/Finch_Grid.hpp @@ -21,6 +21,9 @@ #include #include +#include + +#include namespace Finch { @@ -30,13 +33,21 @@ namespace Finch if ( comm_rank == 0 ) \ std::cout -template +// Scalar is the type stored in the temperature field. It defaults to double, +// which reproduces the original behavior exactly; supplying another arithmetic +// type propagates it through the field, the halo, and the solver arithmetic. +// Note that the *mesh* stays double: cell size and node coordinates are +// geometry, and are not carried by the field scalar type. +template class Grid { public: // Kokkos memory space using memory_space = MemorySpace; + // Field scalar type + using scalar_type = Scalar; + // Default Kokkos execution space for this memory space using exec_space = typename MemorySpace::execution_space; @@ -46,7 +57,7 @@ class Grid using local_mesh_type = Cabana::Grid::LocalMesh; using array_type = - Cabana::Grid::Array; + Cabana::Grid::Array; using view_type = typename array_type::view_type; int comm_rank, comm_size; @@ -56,8 +67,8 @@ class Grid std::array global_low_corner, std::array global_high_corner, std::array ranks_per_dim, std::array bc_types, - Kokkos::Array bc_values, const double initial_temperature ) - : boundary( Boundary( bc_types, bc_values ) ) + Kokkos::Array bc_values, const double initial_temperature ) + : boundary( Boundary( bc_types, bc_values ) ) { initialize( comm, cell_size, global_low_corner, global_high_corner, ranks_per_dim, initial_temperature ); @@ -76,7 +87,7 @@ class Grid std::array global_high_corner, std::array ranks_per_dim, std::array bc_types, const double initial_temperature ) - : boundary( Boundary( bc_types ) ) + : boundary( Boundary( bc_types ) ) { initialize( comm, cell_size, global_low_corner, global_high_corner, ranks_per_dim, initial_temperature ); @@ -127,13 +138,13 @@ class Grid createArrayLayout( global_grid, halo_width, 1, entity_type() ); std::string name( "temperature" ); - T = Cabana::Grid::createArray( name, layout ); - Cabana::Grid::ArrayOp::assign( *T, initial_temperature, + T = Cabana::Grid::createArray( name, layout ); + Cabana::Grid::ArrayOp::assign( *T, Scalar( initial_temperature ), Cabana::Grid::Ghost() ); // create an array to store previous temperature for explicit update // Note: this is an entirely separate array on purpose (no shallow copy) - T0 = Cabana::Grid::createArray( name, layout ); + T0 = Cabana::Grid::createArray( name, layout ); // create halo halo = createHalo( Cabana::Grid::FaceHaloPattern<3>(), halo_width, *T ); @@ -158,7 +169,36 @@ class Grid void output( const int step, const double time ) { - Cabana::Grid::Experimental::BovWriter::writeTimeStep( step, time, *T ); + if constexpr ( std::is_same::value ) + { + Cabana::Grid::Experimental::BovWriter::writeTimeStep( step, time, + *T ); + } + else + { + // BovWriter requires a value type that has both an MpiTraits and a + // BovFormat specialization, so a compound scalar cannot be written + // directly. Project the field onto a plain double array -- its + // numeric value, as defined by Math::ScalarValue -- and write that. + // The extra components of the scalar (for an AD type, the + // derivatives) are deliberately not part of this file; they are + // reported separately by the application. + auto T_real = Cabana::Grid::createArray( + "temperature", T->layout() ); + auto src = T->view(); + auto dst = T_real->view(); + Kokkos::parallel_for( + "project_field_value", + Kokkos::MDRangePolicy>( + { 0, 0, 0 }, { static_cast( src.extent( 0 ) ), + static_cast( src.extent( 1 ) ), + static_cast( src.extent( 2 ) ) } ), + KOKKOS_LAMBDA( const int i, const int j, const int k ) { + dst( i, j, k, 0 ) = Math::value( src( i, j, k, 0 ) ); + } ); + Cabana::Grid::Experimental::BovWriter::writeTimeStep( step, time, + *T_real ); + } } void updateBoundaries() @@ -189,7 +229,7 @@ class Grid std::shared_ptr T0; //! Boundary conditions. - Boundary boundary; + Boundary boundary; }; } // namespace Finch diff --git a/src/Finch_Run.hpp b/src/Finch_Run.hpp index 86543a7..3acd4d3 100644 --- a/src/Finch_Run.hpp +++ b/src/Finch_Run.hpp @@ -24,15 +24,16 @@ namespace Finch { -template +template class Layer { public: using memory_space = MemorySpace; - using sampling_type = Finch::SolidificationData; + using scalar_type = Scalar; + using sampling_type = Finch::SolidificationData; sampling_type solidification_data_; - Layer( Inputs& inputs, Grid& grid ) + Layer( Inputs& inputs, Grid& grid ) { // Only construct if turned on - will otherwise default and immediately // return from any member functions @@ -43,7 +44,8 @@ class Layer // Run the full timestepped loop template void run( ExecutionSpace exec_space, Inputs& inputs, - Grid& grid, MovingBeam& beam, SolverType& fd ) + Grid& grid, MovingBeam& beam, + SolverType& fd ) { // time stepping double& time = inputs.time.time; @@ -75,7 +77,8 @@ class Layer // Run a single timestep template void step( ExecutionSpace exec_space, double& time, const double dt, - Grid grid, MovingBeam& beam, SolverType& fd ) + Grid grid, MovingBeam& beam, + SolverType& fd ) { time += dt; diff --git a/src/Finch_Scalar.hpp b/src/Finch_Scalar.hpp new file mode 100644 index 0000000..64f7fa3 --- /dev/null +++ b/src/Finch_Scalar.hpp @@ -0,0 +1,110 @@ +/**************************************************************************** + * Copyright (c) 2024 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of Finch. Finch is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +/*! + \file Finch_Scalar.hpp + \brief Scalar-type dispatch so the solver can run on types other than double +*/ + +#ifndef Finch_Scalar_H +#define Finch_Scalar_H + +#include + +#include +#include + +namespace Finch +{ +namespace Math +{ + +// Finch's kernels call these instead of Kokkos:: directly so that the field +// scalar type is not required to be a built-in float/double. +// +// For arithmetic types the call forwards to Kokkos::, which is what the +// solver used before this indirection existed -- same device path, same +// architecture-specific implementation, no change in generated code. +// +// For any other type the call is made unqualified so that argument-dependent +// lookup finds an overload in the scalar type's own namespace. That is the +// extension point: a user-defined arithmetic type (a dual number, an interval, +// a truncated-Taylor AD jet) supplies its own exp/fmin/fmax and needs no edit +// here. Finch therefore carries no dependency on any particular AD library. + +template +KOKKOS_INLINE_FUNCTION auto exp( const T& x ) +{ + if constexpr ( std::is_arithmetic::value ) + { + return Kokkos::exp( x ); + } + else + { + using std::exp; + return exp( x ); + } +} + +template +KOKKOS_INLINE_FUNCTION auto fmin( const T& a, const T& b ) +{ + if constexpr ( std::is_arithmetic::value ) + { + return Kokkos::fmin( a, b ); + } + else + { + using std::fmin; + return fmin( a, b ); + } +} + +template +KOKKOS_INLINE_FUNCTION auto fmax( const T& a, const T& b ) +{ + if constexpr ( std::is_arithmetic::value ) + { + return Kokkos::fmax( a, b ); + } + else + { + using std::fmax; + return fmax( a, b ); + } +} + +// Numeric value of a scalar, used where a plain double is required regardless +// of the field type: file output, MPI reductions over bounds, and the +// solidification event records that are handed to downstream tools. +// +// The primary template covers built-in types. A non-arithmetic scalar type +// specializes this to say which of its components is the "value" -- for a +// forward-mode AD type that is the real/zeroth coefficient. +template +struct ScalarValue +{ + KOKKOS_INLINE_FUNCTION static double value( const T& x ) + { + return static_cast( x ); + } +}; + +template +KOKKOS_INLINE_FUNCTION double value( const T& x ) +{ + return ScalarValue::value( x ); +} + +} // namespace Math +} // namespace Finch + +#endif diff --git a/src/Finch_SolidificationData.hpp b/src/Finch_SolidificationData.hpp index ecc3c33..1abad38 100644 --- a/src/Finch_SolidificationData.hpp +++ b/src/Finch_SolidificationData.hpp @@ -31,11 +31,17 @@ #include #include +#include namespace Finch { -template +// Scalar is the field type of the grid this samples. The recorded events stay +// double: they are the hand-off format to downstream tools (ExaCA), so the +// numeric value of each quantity is projected out via Math::value. Carrying +// derivatives into the event records is a separate, larger change -- it would +// make the event view Scalar-valued and give sensitivities of G and R. +template class SolidificationData { using memory_space = MemorySpace; @@ -68,7 +74,7 @@ class SolidificationData // Default constructor SolidificationData() {} // constructor - SolidificationData( const Inputs& inputs, Grid& grid ) + SolidificationData( const Inputs& inputs, Grid& grid ) : mpi_rank_( grid.comm_rank ) , liquidus_( inputs.properties.liquidus ) , dt_( inputs.time.time_step ) @@ -87,7 +93,7 @@ class SolidificationData capacity, nCmpts ); auto local_grid = grid.getLocalGrid(); - using entity_type = typename Grid::entity_type; + using entity_type = typename Grid::entity_type; auto layout = Cabana::Grid::createArrayLayout( local_grid, 1, entity_type() ); auto tm = @@ -95,20 +101,20 @@ class SolidificationData tm_view = tm->view(); } - void updateEvents( Grid& grid, const double time ) + void updateEvents( Grid& grid, const double time ) { // get local copies from grid auto local_mesh = grid.getLocalMesh(); auto T = grid.getTemperature(); auto T0 = grid.getPreviousTemperature(); - using entity_type = typename Grid::entity_type; + using entity_type = typename Grid::entity_type; Cabana::Grid::grid_parallel_for( "local_grid_for", exec_space(), grid.getIndexSpace(), KOKKOS_CLASS_LAMBDA( const int i, const int j, const int k ) { - double temp = T( i, j, k, 0 ); - double temp0 = T0( i, j, k, 0 ); + Scalar temp = T( i, j, k, 0 ); + Scalar temp0 = T0( i, j, k, 0 ); if ( ( temp <= liquidus_ ) && ( temp0 > liquidus_ ) ) { @@ -129,38 +135,41 @@ class SolidificationData events( current_count, 3 ) = tm_view( i, j, k, 0 ); // event solidification time - double m = ( temp - liquidus_ ) / ( temp - temp0 ); - m = fmin( fmax( m, 0.0 ), 1.0 ); - events( current_count, 4 ) = time - m * dt_; + Scalar m = ( temp - liquidus_ ) / ( temp - temp0 ); + m = Math::fmin( Math::fmax( m, Scalar( 0 ) ), + Scalar( 1 ) ); + events( current_count, 4 ) = + time - Math::value( m ) * dt_; // cooling rate - events( current_count, 5 ) = ( temp0 - temp ) / dt_; + events( current_count, 5 ) = + Math::value( ( temp0 - temp ) / dt_ ); // temperature gradient components - events( current_count, 6 ) = + events( current_count, 6 ) = Math::value( ( T( i + 1, j, k, 0 ) - T( i - 1, j, k, 0 ) ) / - ( 2.0 * cell_size_ ); + ( 2.0 * cell_size_ ) ); - events( current_count, 7 ) = + events( current_count, 7 ) = Math::value( ( T( i, j + 1, k, 0 ) - T( i, j - 1, k, 0 ) ) / - ( 2.0 * cell_size_ ); + ( 2.0 * cell_size_ ) ); - events( current_count, 8 ) = + events( current_count, 8 ) = Math::value( ( T( i, j, k + 1, 0 ) - T( i, j, k - 1, 0 ) ) / - ( 2.0 * cell_size_ ); + ( 2.0 * cell_size_ ) ); } } else if ( ( temp > liquidus_ ) && ( temp0 <= liquidus_ ) ) { - double m = ( temp - liquidus_ ) / ( temp - temp0 ); - m = fmin( fmax( m, 0.0 ), 1.0 ); - tm_view( i, j, k, 0 ) = time - m * dt_; + Scalar m = ( temp - liquidus_ ) / ( temp - temp0 ); + m = Math::fmin( Math::fmax( m, Scalar( 0 ) ), Scalar( 1 ) ); + tm_view( i, j, k, 0 ) = time - Math::value( m ) * dt_; } } ); } // Update the solidification data - void update( Grid& grid, const double time ) + void update( Grid& grid, const double time ) { if ( !enabled_ ) { diff --git a/src/Finch_Solver.hpp b/src/Finch_Solver.hpp index 5653edb..39e36b6 100644 --- a/src/Finch_Solver.hpp +++ b/src/Finch_Solver.hpp @@ -20,6 +20,8 @@ #include #include +#include + namespace Finch { @@ -30,9 +32,49 @@ struct DeviceTag { }; +// The material and source inputs the solver treats as differentiable. Held in +// its own struct, templated on the scalar type, so a caller can hand the solver +// values of a type other than double (see makeProperties below). Quantities the +// solver only ever compares against -- solidus and liquidus -- stay double: +// they select a branch rather than entering the arithmetic. +template +struct MaterialProperties +{ + Scalar density; + Scalar specific_heat; + Scalar thermal_conductivity; + Scalar latent_heat; + Scalar absorption; + Scalar two_sigma[3]; +}; + +// Build the properties for a plain double solve directly from the input deck. +// Scalar defaults to double, so existing callers get exactly the previous +// behavior; callers wanting another scalar type request it explicitly and then +// overwrite the members they care about. +template +MaterialProperties makeProperties( const Inputs& db ) +{ + MaterialProperties props; + props.density = db.properties.density; + props.specific_heat = db.properties.specific_heat; + props.thermal_conductivity = db.properties.thermal_conductivity; + props.latent_heat = db.properties.latent_heat; + props.absorption = db.source.absorption; + for ( std::size_t d = 0; d < 3; ++d ) + props.two_sigma[d] = db.source.two_sigma[d]; + return props; +} + template class Solver { + public: + // The field scalar type is whatever the temperature view holds. No extra + // template parameter is needed: making the Cabana array carry a different + // value type is enough to change the arithmetic throughout the solver. + using scalar_type = typename ViewType::non_const_value_type; + protected: // temperature views are default constructed and updated every step. ViewType T_; @@ -44,28 +86,34 @@ class Solver double dt_; double solidus_; double liquidus_; - double rho_cp_; - double rho_Lf_by_dT_; - double k_by_dx2_; + scalar_type rho_cp_; + scalar_type rho_Lf_by_dT_; + scalar_type k_by_dx2_; // heat source parameters double power_; double position_[3]; - double r_[3]; - double A_inv_[3]; - double I0_; + scalar_type r_[3]; + scalar_type A_inv_[3]; + scalar_type I0_; double w_max_; public: Solver( Inputs db, LocalMeshType local_mesh ) + : Solver( db, local_mesh, makeProperties( db ) ) + { + } + + Solver( Inputs db, LocalMeshType local_mesh, + const MaterialProperties& props ) : local_mesh_( local_mesh ) , power_( 0.0 ) { // solution parameter constants double dx = db.space.cell_size; - double rho = db.properties.density; - double cp = db.properties.specific_heat; - double Lf = db.properties.latent_heat; + scalar_type rho = props.density; + scalar_type cp = props.specific_heat; + scalar_type Lf = props.latent_heat; dt_ = db.time.time_step; @@ -77,7 +125,7 @@ class Solver rho_Lf_by_dT_ = rho * Lf / ( liquidus_ - solidus_ ); - k_by_dx2_ = ( db.properties.thermal_conductivity ) / ( dx * dx ); + k_by_dx2_ = ( props.thermal_conductivity ) / ( dx * dx ); // initialize beam position for ( std::size_t d = 0; d < 3; ++d ) @@ -88,11 +136,11 @@ class Solver // heat source parameter constants for ( std::size_t d = 0; d < 3; ++d ) { - r_[d] = db.source.two_sigma[d] / Kokkos::sqrt( 2.0 ); + r_[d] = props.two_sigma[d] / Kokkos::sqrt( 2.0 ); A_inv_[d] = 1.0 / r_[d] / r_[d]; } - I0_ = ( 2.0 * db.source.absorption ) / + I0_ = ( 2.0 * props.absorption ) / ( M_PI * Kokkos::sqrt( M_PI ) * r_[0] * r_[1] * r_[2] ); // cut off for 3 standard deviations from heat source center @@ -136,13 +184,13 @@ class Solver KOKKOS_INLINE_FUNCTION void operator()( HostTag tag, const int i, const int j, const int k ) const { - double x = T0_( i, j, k, 0 ); + scalar_type x = T0_( i, j, k, 0 ); - double dt_by_rho_cp = ( x >= solidus_ && x <= liquidus_ ) - ? dt_ / ( rho_cp_ + rho_Lf_by_dT_ ) - : dt_ / ( rho_cp_ ); + scalar_type dt_by_rho_cp = ( x >= solidus_ && x <= liquidus_ ) + ? dt_ / ( rho_cp_ + rho_Lf_by_dT_ ) + : dt_ / ( rho_cp_ ); - double rhs = laplacian( i, j, k ) + source( tag, i, j, k ); + scalar_type rhs = laplacian( i, j, k ) + source( tag, i, j, k ); T_( i, j, k, 0 ) = x + rhs * dt_by_rho_cp; } @@ -152,13 +200,13 @@ class Solver void operator()( DeviceTag tag, const int i, const int j, const int k ) const { - double x = T0_( i, j, k, 0 ); + scalar_type x = T0_( i, j, k, 0 ); - double dt_by_rho_cp = + scalar_type dt_by_rho_cp = dt_ / ( rho_cp_ + ( x >= solidus_ ) * ( x <= liquidus_ ) * rho_Lf_by_dT_ ); - double rhs = laplacian( i, j, k ) + source( tag, i, j, k ); + scalar_type rhs = laplacian( i, j, k ) + source( tag, i, j, k ); T_( i, j, k, 0 ) = x + rhs * dt_by_rho_cp; } @@ -176,8 +224,11 @@ class Solver // Normalized weight for the gaussian source term: x in exp(-x) KOKKOS_INLINE_FUNCTION - auto weight( const int i, const int j, const int k ) const + scalar_type weight( const int i, const int j, const int k ) const { + // Mesh coordinates and beam position are geometry, not field values, + // and stay double whatever the field scalar type is. Only the + // accumulation below picks up the scalar type, through A_inv_. double grid_loc[3]; double dist_to_beam[3]; int idx[3] = { i, j, k }; @@ -195,49 +246,66 @@ class Solver // Heating source term, device overload. KOKKOS_INLINE_FUNCTION - auto source( DeviceTag, const int i, const int j, const int k ) const + scalar_type source( DeviceTag, const int i, const int j, const int k ) const { - return I0_ * power_ * Kokkos::exp( -weight( i, j, k ) ); + return I0_ * power_ * Math::exp( -weight( i, j, k ) ); } // Heating source term, host overload. KOKKOS_INLINE_FUNCTION - auto source( HostTag, const int i, const int j, const int k ) const + scalar_type source( HostTag, const int i, const int j, const int k ) const { // performance improvements on host: scoping the exponential if ( power_ ) { - double w = weight( i, j, k ); + scalar_type w = weight( i, j, k ); if ( w < w_max_ ) { - return I0_ * power_ * Kokkos::exp( -w ); + return I0_ * power_ * Math::exp( -w ); } else { - return 0.0; + return scalar_type( 0 ); } } else { - return 0.0; + return scalar_type( 0 ); } } }; // Create a solver based on the grid details and simulation inputs. -template -auto createSolver( Inputs db, Grid grid ) +template +auto createSolver( Inputs db, Grid grid ) { - using entity_type = typename Grid::entity_type; - using view_type = typename Grid::view_type; - using mesh_type = typename Grid::local_mesh_type; + using grid_type = Grid; + using entity_type = typename grid_type::entity_type; + using view_type = typename grid_type::view_type; + using mesh_type = typename grid_type::local_mesh_type; auto local_mesh = grid.getLocalMesh(); return Solver( db, local_mesh ); } +// Create a solver with explicitly supplied material properties. Used when the +// properties carry more than a value -- for example seeded AD variables. +template +auto createSolver( Inputs db, Grid grid, + const MaterialProperties& props ) +{ + using grid_type = Grid; + using entity_type = typename grid_type::entity_type; + using view_type = typename grid_type::view_type; + using mesh_type = typename grid_type::local_mesh_type; + + auto local_mesh = grid.getLocalMesh(); + + return Solver( db, local_mesh, props ); +} + } // namespace Finch #endif From d087a00d338ec7b6297c147f107e3888b08fe91a Mon Sep 17 00:00:00 2001 From: Samm-Py <140015960+Samm-Py@users.noreply.github.com> Date: Fri, 7 Aug 2026 11:22:35 -0500 Subject: [PATCH 2/7] Update applications/Finch_OTI.hpp Co-authored-by: Bruno Turcksin --- applications/Finch_OTI.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/Finch_OTI.hpp b/applications/Finch_OTI.hpp index 1457e24..e948587 100644 --- a/applications/Finch_OTI.hpp +++ b/applications/Finch_OTI.hpp @@ -24,7 +24,7 @@ namespace Finch namespace Math { -// Tell Finch which component of an OTI jet is its numeric value. This is all +// Tell Finch which component of an OTI number is its numeric value. This is all // the core library needs to know about the type; every other operation it // performs (arithmetic, comparison, exp, fmin/fmax) is resolved by ADL into // namespace oti. From f0a907df17d35b1e6beabd6f1067b573f86751e0 Mon Sep 17 00:00:00 2001 From: Samm-Py <140015960+Samm-Py@users.noreply.github.com> Date: Fri, 7 Aug 2026 14:09:54 -0500 Subject: [PATCH 3/7] Update src/Finch_Grid.hpp Co-authored-by: Bruno Turcksin --- src/Finch_Grid.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Finch_Grid.hpp b/src/Finch_Grid.hpp index c3d5148..e606e28 100644 --- a/src/Finch_Grid.hpp +++ b/src/Finch_Grid.hpp @@ -34,7 +34,7 @@ namespace Finch std::cout // Scalar is the type stored in the temperature field. It defaults to double, -// which reproduces the original behavior exactly; supplying another arithmetic +// ; supplying another arithmetic // type propagates it through the field, the halo, and the solver arithmetic. // Note that the *mesh* stays double: cell size and node coordinates are // geometry, and are not carried by the field scalar type. From cad3a85256396113f566db07d7b40f2dee902e55 Mon Sep 17 00:00:00 2001 From: Samm-Py Date: Fri, 7 Aug 2026 14:13:30 -0500 Subject: [PATCH 4/7] Address Sparrow integration review comments --- .github/workflows/CI.yml | 18 ++++ .gitignore | 5 + applications/CMakeLists.txt | 26 ++--- ...tivity.cpp => SingleLayer_Sensitivity.cpp} | 98 +++++++++++++------ .../sparrow/Finch_Sparrow.hpp | 55 +++++++---- src/Finch_Grid.hpp | 2 +- src/Finch_Scalar.hpp | 26 ++--- src/Finch_Solver.hpp | 53 +++++----- 8 files changed, 182 insertions(+), 101 deletions(-) rename applications/{Sensitivity.cpp => SingleLayer_Sensitivity.cpp} (73%) rename applications/Finch_OTI.hpp => integrations/sparrow/Finch_Sparrow.hpp (67%) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 27d18e0..c6318f9 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -79,6 +79,11 @@ jobs: repository: ECP-CoPA/Cabana ref: 0.6.1 path: cabana + - name: Checkout Sparrow + uses: actions/checkout@v3 + with: + repository: ORNL-MDF/Sparrow + path: sparrow - name: Build Cabana working-directory: cabana run: | @@ -96,6 +101,7 @@ jobs: cmake -B build \ -D CMAKE_INSTALL_PREFIX=$HOME/finch \ -D CMAKE_PREFIX_PATH="$HOME/cabana;$HOME/json" \ + -D SPARROW_DIR=$GITHUB_WORKSPACE/sparrow/include \ -D CMAKE_BUILD_TYPE=${{ matrix.cmake_build_type }} \ -D CMAKE_CXX_COMPILER=${{ matrix.cxx }} \ -D CMAKE_CXX_FLAGS="-Wall -Wextra -pedantic -Werror" \ @@ -151,6 +157,11 @@ jobs: repository: ECP-CoPA/Cabana ref: 0.6.1 path: cabana + - name: Checkout Sparrow + uses: actions/checkout@v3 + with: + repository: ORNL-MDF/Sparrow + path: sparrow - name: Build Cabana working-directory: cabana run: | @@ -168,6 +179,7 @@ jobs: cmake -B build \ -D CMAKE_INSTALL_PREFIX=$HOME/finch \ -D CMAKE_PREFIX_PATH="$HOME/cabana;$HOME/json" \ + -D SPARROW_DIR=$GITHUB_WORKSPACE/sparrow/include \ -D CMAKE_BUILD_TYPE=${{ matrix.cmake_build_type }} \ -D CMAKE_CXX_COMPILER=${{ matrix.cxx }} \ -D CMAKE_CXX_FLAGS="-Wall -Wextra -pedantic -Werror -I${MPI_LOCATION}/include" \ @@ -211,6 +223,11 @@ jobs: repository: ECP-CoPA/Cabana ref: 0.6.1 path: cabana + - name: Checkout Sparrow + uses: actions/checkout@v3 + with: + repository: ORNL-MDF/Sparrow + path: sparrow - name: Build Cabana working-directory: cabana run: | @@ -228,6 +245,7 @@ jobs: cmake -B build \ -D CMAKE_INSTALL_PREFIX=$HOME/finch \ -D CMAKE_PREFIX_PATH="$HOME/cabana;$HOME/json" \ + -D SPARROW_DIR=$GITHUB_WORKSPACE/sparrow/include \ -D CMAKE_BUILD_TYPE=${{ matrix.cmake_build_type }} \ -D CMAKE_CXX_FLAGS="-Wall -Wextra -pedantic -Werror" cmake --build build --parallel 2 diff --git a/.gitignore b/.gitignore index fd0b775..65a2461 100644 --- a/.gitignore +++ b/.gitignore @@ -13,3 +13,8 @@ Cabana/ *.dat *.csv path*.txt + +# Spike artifacts +build-baseline/ +spike_runs/ +OTI_INTEGRATION_NOTES.md diff --git a/applications/CMakeLists.txt b/applications/CMakeLists.txt index 5610efb..eb4b5e5 100644 --- a/applications/CMakeLists.txt +++ b/applications/CMakeLists.txt @@ -2,22 +2,24 @@ add_executable(finch SingleLayer.cpp) target_link_libraries(finch Core) install(TARGETS finch DESTINATION ${CMAKE_INSTALL_BINDIR}) -# Optional forward-mode AD sensitivity application. Built only when the -# header-only cpp_oti_lib is available; Finch itself never depends on it. -find_path(CPP_OTI_LIB_INCLUDE_DIR otinum/otinum.hpp +# Optional forward-mode AD sensitivity application. Sparrow is header-only, so +# its include directory is the complete dependency integration; no Sparrow +# library target is built or linked, and Finch itself never depends on it. +find_path(SPARROW_INCLUDE_DIR otinum/otinum.hpp HINTS - ${CPP_OTI_LIB_DIR} - ${CMAKE_CURRENT_SOURCE_DIR}/../../cpp_oti_lib/include - DOC "Path to the cpp_oti_lib include directory") + ${SPARROW_DIR} + ${CMAKE_CURRENT_SOURCE_DIR}/../../Sparrow/include + DOC "Path to the Sparrow include directory") -if(CPP_OTI_LIB_INCLUDE_DIR) - message(STATUS "Finch sensitivity: cpp_oti_lib found at ${CPP_OTI_LIB_INCLUDE_DIR}") - add_executable(finch_sensitivity Sensitivity.cpp) +if(SPARROW_INCLUDE_DIR) + message(STATUS "Finch sensitivity: Sparrow found at ${SPARROW_INCLUDE_DIR}") + add_executable(finch_sensitivity SingleLayer_Sensitivity.cpp) target_link_libraries(finch_sensitivity Core) + target_compile_definitions(finch_sensitivity PRIVATE OTI_ENABLE_KOKKOS) target_include_directories(finch_sensitivity PRIVATE - ${CPP_OTI_LIB_INCLUDE_DIR} - ${CMAKE_CURRENT_SOURCE_DIR}) + ${SPARROW_INCLUDE_DIR} + ${PROJECT_SOURCE_DIR}/integrations/sparrow) install(TARGETS finch_sensitivity DESTINATION ${CMAKE_INSTALL_BINDIR}) else() - message(STATUS "Finch sensitivity: cpp_oti_lib not found, skipping finch_sensitivity") + message(STATUS "Finch sensitivity: Sparrow not found, skipping finch_sensitivity") endif() diff --git a/applications/Sensitivity.cpp b/applications/SingleLayer_Sensitivity.cpp similarity index 73% rename from applications/Sensitivity.cpp rename to applications/SingleLayer_Sensitivity.cpp index 31309cd..8f88174 100644 --- a/applications/Sensitivity.cpp +++ b/applications/SingleLayer_Sensitivity.cpp @@ -1,5 +1,17 @@ /**************************************************************************** - * Parameter sensitivity of the Finch heat solve, by forward-mode AD. + * Copyright (c) 2024 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of Finch. Finch is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +/**************************************************************************** + * Parameter sensitivity of the Finch single-layer heat solve, by forward-mode + * AD. * * Runs the single-layer heat transport problem once with OTI-valued material * and source parameters, obtaining d(QoI)/dp for every parameter from that one @@ -12,6 +24,7 @@ #include #include #include +#include #include #include @@ -19,18 +32,31 @@ #include #include "Finch_Core.hpp" -#include "Finch_OTI.hpp" +#include "Finch_Sparrow.hpp" namespace { constexpr int NP = Finch::Sensitivity::NumParameters; -// First-order jets: one variable per differentiated parameter, derivatives -// through order one. Raising the order to 2 here is the only change needed to -// obtain the full parameter Hessian of the same quantities. +// First-order OTI numbers: one variable per differentiated parameter. This +// application computes and reports gradients only; second derivatives are not +// evaluated. using OTI = oti::otinum; +template +MPI_Datatype mpi_datatype() +{ + static_assert( std::is_same::value || + std::is_same::value, + "Sensitivity MPI reductions support float and double " + "coefficients" ); + if constexpr ( std::is_same::value ) + return MPI_FLOAT; + else + return MPI_DOUBLE; +} + // Quantities of interest evaluated on the final temperature field. Both are // smooth functionals of the field, which makes them meaningful finite- // difference references. @@ -46,7 +72,7 @@ struct QoI // Run the transient solve to completion and evaluate the QoIs. template QoI solve( MPI_Comm comm, Finch::Inputs db, - const Finch::MaterialProperties& props ) + const Finch::SolverParameters& params ) { std::array bc_types = { "adiabatic", "adiabatic", "adiabatic", "adiabatic", @@ -57,7 +83,7 @@ QoI solve( MPI_Comm comm, Finch::Inputs db, db.space.global_high_corner, db.space.ranks_per_dim, bc_types, db.space.initial_temperature ); - auto fd = Finch::createSolver( db, grid, props ); + auto fd = Finch::createSolver( db, grid, params ); Finch::MovingBeam beam( db.source.scan_path_file ); @@ -88,29 +114,37 @@ QoI solve( MPI_Comm comm, Finch::Inputs db, grid.gather(); } - // Evaluate the QoIs on the host. Copying the field to a host mirror keeps - // this independent of Kokkos reducer support for compound scalar types, - // which is not needed for a once-per-run diagnostic. + // Evaluate the sum on the active Kokkos backend. For the OTI solve this + // deliberately reduces a compound Sparrow scalar, exercising its device + // annotations, additive identity, and operator+= integration with Kokkos. auto T = grid.getTemperature(); - auto T_host = Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), T ); - auto owned = grid.getIndexSpace(); QoI q; q.sum = Scalar( 0 ); - for ( int i = owned.min( 0 ); i < owned.max( 0 ); ++i ) - for ( int j = owned.min( 1 ); j < owned.max( 1 ); ++j ) - for ( int k = owned.min( 2 ); k < owned.max( 2 ); ++k ) - q.sum = q.sum + T_host( i, j, k, 0 ); - - q.probe = T_host( ( owned.min( 0 ) + owned.max( 0 ) ) / 2, - ( owned.min( 1 ) + owned.max( 1 ) ) / 2, - ( owned.min( 2 ) + owned.max( 2 ) ) / 2, 0 ); - - // Reduce the sum across ranks. An OTI jet is a contiguous block of - // coefficients, and summing jets is summing coefficients elementwise, so a - // plain MPI_DOUBLE reduction over ncoeffs entries is exact -- no derived - // datatype or user-defined operator is required. + Cabana::Grid::grid_parallel_reduce( + "sensitivity_temperature_sum", ExecSpace(), owned, + KOKKOS_LAMBDA( const int i, const int j, const int k, + Scalar& local_sum ) { + local_sum += T( i, j, k, 0 ); + }, + q.sum ); + + // Copy only the single probe value to the host rather than mirroring the + // complete field. A rank-zero view keeps this path valid for device-only + // memory spaces. + const int probe_i = ( owned.min( 0 ) + owned.max( 0 ) ) / 2; + const int probe_j = ( owned.min( 1 ) + owned.max( 1 ) ) / 2; + const int probe_k = ( owned.min( 2 ) + owned.max( 2 ) ) / 2; + auto probe_device = Kokkos::subview( T, probe_i, probe_j, probe_k, 0 ); + auto probe_host = Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace(), + probe_device ); + q.probe = probe_host(); + + // Reduce the sum across ranks. An OTI number is a contiguous block of + // coefficients, and summing OTI numbers is summing coefficients + // elementwise, so a plain MPI reduction over ncoeffs coefficient values is + // exact -- no derived datatype or user-defined operator is required. int comm_size; MPI_Comm_size( comm, &comm_size ); if ( comm_size > 1 ) @@ -122,9 +156,10 @@ QoI solve( MPI_Comm comm, Finch::Inputs db, } else { + using coeff_type = typename Scalar::coeff_type; Scalar local = q.sum; - MPI_Allreduce( &local[0], &q.sum[0], Scalar::ncoeffs, MPI_DOUBLE, - MPI_SUM, comm ); + MPI_Allreduce( &local[0], &q.sum[0], Scalar::ncoeffs, + mpi_datatype(), MPI_SUM, comm ); } } @@ -227,7 +262,8 @@ void run( MPI_Comm comm, int argc, char* argv[] ) for ( int q = 0; q < 2; ++q ) { - const OTI& jet = ( q == 0 ) ? oti_result.sum : oti_result.probe; + const OTI& oti_number = + ( q == 0 ) ? oti_result.sum : oti_result.probe; const std::array& fd = ( q == 0 ) ? fd_sum : fd_probe; std::printf( " d(%s)/dp\n", qoi_name[q] ); @@ -246,15 +282,15 @@ void run( MPI_Comm comm, int argc, char* argv[] ) { typename OTI::alpha_type alpha{}; alpha[i] = 1; - ref = - std::max( ref, std::fabs( nominal[i] * jet.partial( alpha ) ) ); + ref = std::max( + ref, std::fabs( nominal[i] * oti_number.partial( alpha ) ) ); } for ( int i = 0; i < NP; ++i ) { typename OTI::alpha_type alpha{}; alpha[i] = 1; - double d_oti = jet.partial( alpha ); + double d_oti = oti_number.partial( alpha ); double scale = std::max( std::fabs( d_oti ), std::fabs( fd[i] ) ); bool negligible = nominal[i] != 0.0 && diff --git a/applications/Finch_OTI.hpp b/integrations/sparrow/Finch_Sparrow.hpp similarity index 67% rename from applications/Finch_OTI.hpp rename to integrations/sparrow/Finch_Sparrow.hpp index e948587..7ad9559 100644 --- a/applications/Finch_OTI.hpp +++ b/integrations/sparrow/Finch_Sparrow.hpp @@ -1,14 +1,25 @@ /**************************************************************************** - * OTI (order-truncated imaginary) scalar type glue for Finch. + * Copyright (c) 2024 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of Finch. Finch is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +/**************************************************************************** + * Sparrow OTI (order-truncated imaginary) scalar type glue for Finch. * - * This header is the ONLY place where Finch and cpp_oti_lib meet. It lives in - * applications/ rather than src/ on purpose: Finch's core carries no dependency - * on any AD library, and the coupling is a single ScalarValue specialization - * plus a parameter-seeding helper. + * This header is the ONLY place where Finch and Sparrow meet. It lives in + * integrations/sparrow rather than Finch's core or an application: the core + * carries no dependency on any AD library, and the coupling is a single + * ScalarValue specialization plus a parameter-seeding helper. ****************************************************************************/ -#ifndef Finch_OTI_H -#define Finch_OTI_H +#ifndef Finch_Sparrow_H +#define Finch_Sparrow_H #include #include @@ -31,10 +42,12 @@ namespace Math template struct ScalarValue> { - KOKKOS_INLINE_FUNCTION static double + using value_type = Coeff; + + KOKKOS_INLINE_FUNCTION static value_type value( const oti::otinum& x ) { - return static_cast( x.real() ); + return x.real(); } }; @@ -77,7 +90,8 @@ inline const char* name( int p ) case TwoSigma: return "two_sigma"; default: - return "unknown"; + throw std::out_of_range( + "Sensitivity parameter index is out of range" ); } } @@ -98,7 +112,8 @@ inline const char* units( int p ) case TwoSigma: return "m"; default: - return ""; + throw std::out_of_range( + "Sensitivity parameter index is out of range" ); } } @@ -128,17 +143,17 @@ inline std::array nominal( const Inputs& db ) // (values seeded as independent variables), so that both paths perturb exactly // the same quantities. template -MaterialProperties build( const std::array& p ) +SolverParameters build( const std::array& p ) { - MaterialProperties props; - props.density = p[Density]; - props.specific_heat = p[SpecificHeat]; - props.thermal_conductivity = p[ThermalConductivity]; - props.latent_heat = p[LatentHeat]; - props.absorption = p[Absorption]; + SolverParameters params; + params.density = p[Density]; + params.specific_heat = p[SpecificHeat]; + params.thermal_conductivity = p[ThermalConductivity]; + params.latent_heat = p[LatentHeat]; + params.absorption = p[Absorption]; for ( int d = 0; d < 3; ++d ) - props.two_sigma[d] = p[TwoSigma]; - return props; + params.two_sigma[d] = p[TwoSigma]; + return params; } // Seed every parameter as an independent OTI variable at its nominal value. diff --git a/src/Finch_Grid.hpp b/src/Finch_Grid.hpp index e606e28..c3d5148 100644 --- a/src/Finch_Grid.hpp +++ b/src/Finch_Grid.hpp @@ -34,7 +34,7 @@ namespace Finch std::cout // Scalar is the type stored in the temperature field. It defaults to double, -// ; supplying another arithmetic +// which reproduces the original behavior exactly; supplying another arithmetic // type propagates it through the field, the halo, and the solver arithmetic. // Note that the *mesh* stays double: cell size and node coordinates are // geometry, and are not carried by the field scalar type. diff --git a/src/Finch_Scalar.hpp b/src/Finch_Scalar.hpp index 64f7fa3..b207780 100644 --- a/src/Finch_Scalar.hpp +++ b/src/Finch_Scalar.hpp @@ -36,9 +36,9 @@ namespace Math // // For any other type the call is made unqualified so that argument-dependent // lookup finds an overload in the scalar type's own namespace. That is the -// extension point: a user-defined arithmetic type (a dual number, an interval, -// a truncated-Taylor AD jet) supplies its own exp/fmin/fmax and needs no edit -// here. Finch therefore carries no dependency on any particular AD library. +// extension point: a user-defined scalar type supplies its own exp/fmin/fmax +// and needs no edit here. Finch therefore carries no dependency on any +// particular AD library. template KOKKOS_INLINE_FUNCTION auto exp( const T& x ) @@ -82,9 +82,11 @@ KOKKOS_INLINE_FUNCTION auto fmax( const T& a, const T& b ) } } -// Numeric value of a scalar, used where a plain double is required regardless -// of the field type: file output, MPI reductions over bounds, and the -// solidification event records that are handed to downstream tools. +// Numeric value of a scalar. The trait exposes the underlying arithmetic type +// so an AD scalar backed by float remains float and one backed by double remains +// double. Call sites whose storage format requires double (file output and the +// solidification event records handed to downstream tools) convert at that +// boundary rather than forcing every scalar integration to double here. // // The primary template covers built-in types. A non-arithmetic scalar type // specializes this to say which of its components is the "value" -- for a @@ -92,14 +94,16 @@ KOKKOS_INLINE_FUNCTION auto fmax( const T& a, const T& b ) template struct ScalarValue { - KOKKOS_INLINE_FUNCTION static double value( const T& x ) - { - return static_cast( x ); - } + using value_type = T; + + KOKKOS_INLINE_FUNCTION static value_type value( const T& x ) { return x; } }; template -KOKKOS_INLINE_FUNCTION double value( const T& x ) +using scalar_value_t = typename ScalarValue::value_type; + +template +KOKKOS_INLINE_FUNCTION scalar_value_t value( const T& x ) { return ScalarValue::value( x ); } diff --git a/src/Finch_Solver.hpp b/src/Finch_Solver.hpp index 39e36b6..96cd7fa 100644 --- a/src/Finch_Solver.hpp +++ b/src/Finch_Solver.hpp @@ -34,11 +34,11 @@ struct DeviceTag // The material and source inputs the solver treats as differentiable. Held in // its own struct, templated on the scalar type, so a caller can hand the solver -// values of a type other than double (see makeProperties below). Quantities the -// solver only ever compares against -- solidus and liquidus -- stay double: -// they select a branch rather than entering the arithmetic. +// values of a type other than double (see makeSolverParameters below). +// Quantities the solver only ever compares against -- solidus and liquidus -- +// stay double: they select a branch rather than entering the arithmetic. template -struct MaterialProperties +struct SolverParameters { Scalar density; Scalar specific_heat; @@ -48,22 +48,22 @@ struct MaterialProperties Scalar two_sigma[3]; }; -// Build the properties for a plain double solve directly from the input deck. +// Build the parameters for a plain double solve directly from the input deck. // Scalar defaults to double, so existing callers get exactly the previous // behavior; callers wanting another scalar type request it explicitly and then // overwrite the members they care about. template -MaterialProperties makeProperties( const Inputs& db ) +SolverParameters makeSolverParameters( const Inputs& db ) { - MaterialProperties props; - props.density = db.properties.density; - props.specific_heat = db.properties.specific_heat; - props.thermal_conductivity = db.properties.thermal_conductivity; - props.latent_heat = db.properties.latent_heat; - props.absorption = db.source.absorption; + SolverParameters params; + params.density = db.properties.density; + params.specific_heat = db.properties.specific_heat; + params.thermal_conductivity = db.properties.thermal_conductivity; + params.latent_heat = db.properties.latent_heat; + params.absorption = db.source.absorption; for ( std::size_t d = 0; d < 3; ++d ) - props.two_sigma[d] = db.source.two_sigma[d]; - return props; + params.two_sigma[d] = db.source.two_sigma[d]; + return params; } template @@ -100,20 +100,20 @@ class Solver public: Solver( Inputs db, LocalMeshType local_mesh ) - : Solver( db, local_mesh, makeProperties( db ) ) + : Solver( db, local_mesh, makeSolverParameters( db ) ) { } Solver( Inputs db, LocalMeshType local_mesh, - const MaterialProperties& props ) + const SolverParameters& params ) : local_mesh_( local_mesh ) , power_( 0.0 ) { // solution parameter constants double dx = db.space.cell_size; - scalar_type rho = props.density; - scalar_type cp = props.specific_heat; - scalar_type Lf = props.latent_heat; + scalar_type rho = params.density; + scalar_type cp = params.specific_heat; + scalar_type Lf = params.latent_heat; dt_ = db.time.time_step; @@ -125,7 +125,7 @@ class Solver rho_Lf_by_dT_ = rho * Lf / ( liquidus_ - solidus_ ); - k_by_dx2_ = ( props.thermal_conductivity ) / ( dx * dx ); + k_by_dx2_ = ( params.thermal_conductivity ) / ( dx * dx ); // initialize beam position for ( std::size_t d = 0; d < 3; ++d ) @@ -136,11 +136,11 @@ class Solver // heat source parameter constants for ( std::size_t d = 0; d < 3; ++d ) { - r_[d] = props.two_sigma[d] / Kokkos::sqrt( 2.0 ); + r_[d] = params.two_sigma[d] / Kokkos::sqrt( 2.0 ); A_inv_[d] = 1.0 / r_[d] / r_[d]; } - I0_ = ( 2.0 * props.absorption ) / + I0_ = ( 2.0 * params.absorption ) / ( M_PI * Kokkos::sqrt( M_PI ) * r_[0] * r_[1] * r_[2] ); // cut off for 3 standard deviations from heat source center @@ -290,11 +290,12 @@ auto createSolver( Inputs db, Grid grid ) return Solver( db, local_mesh ); } -// Create a solver with explicitly supplied material properties. Used when the -// properties carry more than a value -- for example seeded AD variables. +// Create a solver with explicitly supplied material and source parameters. +// Used when the parameters carry more than a value -- for example seeded AD +// variables. template auto createSolver( Inputs db, Grid grid, - const MaterialProperties& props ) + const SolverParameters& params ) { using grid_type = Grid; using entity_type = typename grid_type::entity_type; @@ -303,7 +304,7 @@ auto createSolver( Inputs db, Grid grid, auto local_mesh = grid.getLocalMesh(); - return Solver( db, local_mesh, props ); + return Solver( db, local_mesh, params ); } } // namespace Finch From 8ce26e5427653d56ad09565aa471f4e80e29e167 Mon Sep 17 00:00:00 2001 From: Samm-Py Date: Fri, 7 Aug 2026 14:32:30 -0500 Subject: [PATCH 5/7] Clarify scalar template comment --- src/Finch_Grid.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Finch_Grid.hpp b/src/Finch_Grid.hpp index c3d5148..a2425a0 100644 --- a/src/Finch_Grid.hpp +++ b/src/Finch_Grid.hpp @@ -33,9 +33,9 @@ namespace Finch if ( comm_rank == 0 ) \ std::cout -// Scalar is the type stored in the temperature field. It defaults to double, -// which reproduces the original behavior exactly; supplying another arithmetic -// type propagates it through the field, the halo, and the solver arithmetic. +// Scalar is the type stored in the temperature field. It defaults to double; +// supplying another arithmetic type propagates it through the field, the halo, +// and the solver arithmetic. // Note that the *mesh* stays double: cell size and node coordinates are // geometry, and are not carried by the field scalar type. template From 470c38a6f3a89afa14cbe66865992dc6465a2181 Mon Sep 17 00:00:00 2001 From: Samm-Py Date: Fri, 7 Aug 2026 14:38:51 -0500 Subject: [PATCH 6/7] Avoid private Sparrow checkout in Finch CI --- .github/workflows/CI.yml | 18 ------------------ 1 file changed, 18 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index c6318f9..27d18e0 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -79,11 +79,6 @@ jobs: repository: ECP-CoPA/Cabana ref: 0.6.1 path: cabana - - name: Checkout Sparrow - uses: actions/checkout@v3 - with: - repository: ORNL-MDF/Sparrow - path: sparrow - name: Build Cabana working-directory: cabana run: | @@ -101,7 +96,6 @@ jobs: cmake -B build \ -D CMAKE_INSTALL_PREFIX=$HOME/finch \ -D CMAKE_PREFIX_PATH="$HOME/cabana;$HOME/json" \ - -D SPARROW_DIR=$GITHUB_WORKSPACE/sparrow/include \ -D CMAKE_BUILD_TYPE=${{ matrix.cmake_build_type }} \ -D CMAKE_CXX_COMPILER=${{ matrix.cxx }} \ -D CMAKE_CXX_FLAGS="-Wall -Wextra -pedantic -Werror" \ @@ -157,11 +151,6 @@ jobs: repository: ECP-CoPA/Cabana ref: 0.6.1 path: cabana - - name: Checkout Sparrow - uses: actions/checkout@v3 - with: - repository: ORNL-MDF/Sparrow - path: sparrow - name: Build Cabana working-directory: cabana run: | @@ -179,7 +168,6 @@ jobs: cmake -B build \ -D CMAKE_INSTALL_PREFIX=$HOME/finch \ -D CMAKE_PREFIX_PATH="$HOME/cabana;$HOME/json" \ - -D SPARROW_DIR=$GITHUB_WORKSPACE/sparrow/include \ -D CMAKE_BUILD_TYPE=${{ matrix.cmake_build_type }} \ -D CMAKE_CXX_COMPILER=${{ matrix.cxx }} \ -D CMAKE_CXX_FLAGS="-Wall -Wextra -pedantic -Werror -I${MPI_LOCATION}/include" \ @@ -223,11 +211,6 @@ jobs: repository: ECP-CoPA/Cabana ref: 0.6.1 path: cabana - - name: Checkout Sparrow - uses: actions/checkout@v3 - with: - repository: ORNL-MDF/Sparrow - path: sparrow - name: Build Cabana working-directory: cabana run: | @@ -245,7 +228,6 @@ jobs: cmake -B build \ -D CMAKE_INSTALL_PREFIX=$HOME/finch \ -D CMAKE_PREFIX_PATH="$HOME/cabana;$HOME/json" \ - -D SPARROW_DIR=$GITHUB_WORKSPACE/sparrow/include \ -D CMAKE_BUILD_TYPE=${{ matrix.cmake_build_type }} \ -D CMAKE_CXX_FLAGS="-Wall -Wextra -pedantic -Werror" cmake --build build --parallel 2 From aaca526f3ee79fb4289e4440312ce5fb6d1905ae Mon Sep 17 00:00:00 2001 From: Samm-Py Date: Mon, 10 Aug 2026 12:56:44 -0500 Subject: [PATCH 7/7] Add a parameter sensitivity example Replace the ad hoc spike run directories with an example under examples/single_line_sensitivity, laid out like the existing single_line example. run_example.sh performs the reference double solve and then four sensitivity configurations: the default finite difference step, a refined step, latent heat removed, and four MPI ranks. The four rank stage now reuses inputs.json Drop the spike specific .gitignore entries, which belong in a local exclude file rather than in the repository. --- .gitignore | 5 --- examples/single_line_sensitivity/inputs.json | 39 +++++++++++++++++++ .../inputs_nolatent.json | 39 +++++++++++++++++++ .../single_line_sensitivity/run_example.sh | 39 +++++++++++++++++++ .../scan_path_small.txt | 3 ++ 5 files changed, 120 insertions(+), 5 deletions(-) create mode 100644 examples/single_line_sensitivity/inputs.json create mode 100644 examples/single_line_sensitivity/inputs_nolatent.json create mode 100755 examples/single_line_sensitivity/run_example.sh create mode 100644 examples/single_line_sensitivity/scan_path_small.txt diff --git a/.gitignore b/.gitignore index 65a2461..fd0b775 100644 --- a/.gitignore +++ b/.gitignore @@ -13,8 +13,3 @@ Cabana/ *.dat *.csv path*.txt - -# Spike artifacts -build-baseline/ -spike_runs/ -OTI_INTEGRATION_NOTES.md diff --git a/examples/single_line_sensitivity/inputs.json b/examples/single_line_sensitivity/inputs.json new file mode 100644 index 0000000..34b399e --- /dev/null +++ b/examples/single_line_sensitivity/inputs.json @@ -0,0 +1,39 @@ +{ + "time": + { + "Co": 0.125, + "start_time": 0.0, + "end_time": 0.0015, + "total_output_steps": 2, + "total_monitor_steps": 10 + }, + "space": + { + "initial_temperature": 300.0, + "cell_size": 10e-6, + "global_low_corner": [-2e-4, -2e-4, -2e-4], + "global_high_corner": [3e-4, 2e-4, 0.0], + "ranks_per_dim": [1, 1, 1] + }, + "properties": + { + "density": 7500.0, + "specific_heat": 750.0, + "thermal_conductivity": 25.0, + "latent_heat": 2e5, + "solidus": 1410.0, + "liquidus": 1620.0 + }, + "source": + { + "absorption": 0.3, + "two_sigma": [60e-6, 60e-6, 60e-6], + "scan_path_file": "scan_path_small.txt" + }, + "sampling": + { + "type": "solidification_data", + "format": "default", + "directory_name": "solidification" + } +} diff --git a/examples/single_line_sensitivity/inputs_nolatent.json b/examples/single_line_sensitivity/inputs_nolatent.json new file mode 100644 index 0000000..dc1042e --- /dev/null +++ b/examples/single_line_sensitivity/inputs_nolatent.json @@ -0,0 +1,39 @@ +{ + "time": + { + "Co": 0.125, + "start_time": 0.0, + "end_time": 0.0015, + "total_output_steps": 2, + "total_monitor_steps": 10 + }, + "space": + { + "initial_temperature": 300.0, + "cell_size": 10e-6, + "global_low_corner": [-2e-4, -2e-4, -2e-4], + "global_high_corner": [3e-4, 2e-4, 0.0], + "ranks_per_dim": [1, 1, 1] + }, + "properties": + { + "density": 7500.0, + "specific_heat": 750.0, + "thermal_conductivity": 25.0, + "latent_heat": 0.0, + "solidus": 1410.0, + "liquidus": 1620.0 + }, + "source": + { + "absorption": 0.3, + "two_sigma": [60e-6, 60e-6, 60e-6], + "scan_path_file": "scan_path_small.txt" + }, + "sampling": + { + "type": "solidification_data", + "format": "default", + "directory_name": "solidification" + } +} diff --git a/examples/single_line_sensitivity/run_example.sh b/examples/single_line_sensitivity/run_example.sh new file mode 100755 index 0000000..8826360 --- /dev/null +++ b/examples/single_line_sensitivity/run_example.sh @@ -0,0 +1,39 @@ +#!/bin/sh + +# Run from this directory +cd ${0%/*} || exit 1 + +# source executables +FINCH_DIR=`pwd`/../.. +application=$FINCH_DIR/build/install/bin/finch +sensitivity=$FINCH_DIR/build/install/bin/finch_sensitivity + +# reference double solve, writing the temperature field and solidification data +echo "### reference double solve" +$application -i inputs.json + +# one OTI solve gives every parameter derivative, checked against central finite +# differences. Requires Sparrow at configure time; without it the +# finch_sensitivity target is not built and this example cannot run. +echo "### sensitivity, default finite difference step" +$sensitivity -i inputs.json + +# The default step straddles the latent heat branch, so density, specific_heat +# and absorption disagree by ~28% above. The derivative is not wrong: a smaller +# step recovers agreement, which is the opposite of how finite difference error +# normally behaves and is the signature of a discontinuity rather than a bug. +echo "### sensitivity, refined finite difference step" +FINCH_FD_STEP=1e-8 $sensitivity -i inputs.json + +# Removing the latent heat branch entirely makes every parameter agree at the +# default step, confirming the diagnosis above. +echo "### sensitivity, latent heat removed" +$sensitivity -i inputs_nolatent.json + +# Same case and same input file on four ranks: ranks_per_dim is ignored when it +# does not match the communicator size, leaving the decomposition to +# MPI_Dims_create, which gives 2x2x1 here. T_sum and its derivatives match the +# single rank values; T_probe does not, since the probe is the centre of the +# owned index space and so is a different node under domain decomposition. +echo "### sensitivity, four ranks" +mpirun -np 4 $sensitivity -i inputs.json diff --git a/examples/single_line_sensitivity/scan_path_small.txt b/examples/single_line_sensitivity/scan_path_small.txt new file mode 100644 index 0000000..3f2889f --- /dev/null +++ b/examples/single_line_sensitivity/scan_path_small.txt @@ -0,0 +1,3 @@ +Mode X Y Z Power Parameter +1 0.000 0.00 0.0 195 0.0005 +0 0.0002 0.00 0.0 195 0.8