From 9b764c38b4d4cfb5d2f3dd1ba4b15a0766f8d802 Mon Sep 17 00:00:00 2001 From: MattRolchigo <65y@ornl.gov> Date: Thu, 19 Dec 2024 15:18:44 -0500 Subject: [PATCH 1/2] Option for temperature-dependent thermal conductivity --- src/Finch_Inputs.hpp | 37 ++++++++++---- src/Finch_Solver.hpp | 117 +++++++++++++++++++++++++++++++++++-------- 2 files changed, 123 insertions(+), 31 deletions(-) diff --git a/src/Finch_Inputs.hpp b/src/Finch_Inputs.hpp index e30895a..a50f760 100644 --- a/src/Finch_Inputs.hpp +++ b/src/Finch_Inputs.hpp @@ -91,7 +91,8 @@ struct Properties { double density; double specific_heat; - double thermal_conductivity; + double thermal_conductivity_0; + double thermal_conductivity_T; double thermal_diffusivity; double latent_heat; double solidus; @@ -211,7 +212,8 @@ class Inputs Info << "Properties:" << std::endl; Info << " Density: " << properties.density << std::endl; Info << " Specific Heat: " << properties.specific_heat << std::endl; - Info << " Thermal Conductivity: " << properties.thermal_conductivity + Info << " Thermal Conductivity: " << properties.thermal_conductivity_0 + << " + " << properties.thermal_conductivity_T << " * T" << std::endl; Info << " Latent Heat: " << properties.latent_heat << std::endl; Info << " Solidus: " << properties.solidus << std::endl; @@ -264,19 +266,22 @@ class Inputs return filename_s; } - void parseInputFile( MPI_Comm comm, const std::string filename ) + void parseInputFile( MPI_Comm comm, const std::string filename, + const double est_max_temperature = 5000.0 ) { readInput( filename ); write(); - // create auxiliary properties - properties.thermal_diffusivity = - ( properties.thermal_conductivity ) / + // create auxiliary properties - estimate time step based on a maximum + // temperature and the parsed thermal conducitivity values + const double est_max_thermal_diffusivity = + ( properties.thermal_conductivity_0 + + properties.thermal_conductivity_T * est_max_temperature ) / ( properties.density * properties.specific_heat ); time.time_step = ( time.Co * space.cell_size * space.cell_size ) / - ( properties.thermal_diffusivity ); + ( est_max_thermal_diffusivity ); Info << "Calculated time step: " << time.time_step << std::endl; @@ -335,8 +340,22 @@ class Inputs // Read properties components properties.density = db["properties"]["density"]; properties.specific_heat = db["properties"]["specific_heat"]; - properties.thermal_conductivity = - db["properties"]["thermal_conductivity"]; + // Thermal conductvity as a function of temperature in the form a + b * + // Temperature. If only one value is given, a fixed thermal conductivity + // is used + if ( db["properties"]["thermal_conductivity"].size() > 1 ) + { + properties.thermal_conductivity_0 = + db["properties"]["thermal_conductivity"][0]; + properties.thermal_conductivity_T = + db["properties"]["thermal_conductivity"][1]; + } + else + { + properties.thermal_conductivity_0 = + db["properties"]["thermal_conductivity"]; + properties.thermal_conductivity_T = 0.0; + } properties.latent_heat = db["properties"]["latent_heat"]; properties.solidus = db["properties"]["solidus"]; properties.liquidus = db["properties"]["liquidus"]; diff --git a/src/Finch_Solver.hpp b/src/Finch_Solver.hpp index 5653edb..e8d90c0 100644 --- a/src/Finch_Solver.hpp +++ b/src/Finch_Solver.hpp @@ -46,7 +46,9 @@ class Solver double liquidus_; double rho_cp_; double rho_Lf_by_dT_; - double k_by_dx2_; + double inv_dx2_; + double k_0_; + double k_T_; // heat source parameters double power_; @@ -77,7 +79,10 @@ class Solver rho_Lf_by_dT_ = rho * Lf / ( liquidus_ - solidus_ ); - k_by_dx2_ = ( db.properties.thermal_conductivity ) / ( dx * dx ); + inv_dx2_ = 1.0 / ( dx * dx ); + + k_0_ = db.properties.thermal_conductivity_0; + k_T_ = db.properties.thermal_conductivity_T; // initialize beam position for ( std::size_t d = 0; d < 3; ++d ) @@ -136,15 +141,41 @@ 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 ); - - double dt_by_rho_cp = ( x >= solidus_ && x <= liquidus_ ) - ? dt_ / ( rho_cp_ + rho_Lf_by_dT_ ) - : dt_ / ( rho_cp_ ); + // First nearest neighbor stencil for cell at i,j,k: temperature and + // thermal conductivity + const double temp_local = T0_( i, j, k, 0 ); + const double temp_px = T0_( i + 1, j, k, 0 ); + const double temp_nx = T0_( i - 1, j, k, 0 ); + const double temp_py = T0_( i, j + 1, k, 0 ); + const double temp_ny = T0_( i, j - 1, k, 0 ); + const double temp_pz = T0_( i, j, k + 1, 0 ); + const double temp_nz = T0_( i, j, k - 1, 0 ); + + const double kappa_local = kappa_of_temperature( temp_local ); + const double kappa_px = kappa_of_temperature( temp_px ); + const double kappa_nx = kappa_of_temperature( temp_nx ); + const double kappa_py = kappa_of_temperature( temp_py ); + const double kappa_ny = kappa_of_temperature( temp_ny ); + const double kappa_pz = kappa_of_temperature( temp_pz ); + const double kappa_nz = kappa_of_temperature( temp_nz ); - double rhs = laplacian( i, j, k ) + source( tag, i, j, k ); - - T_( i, j, k, 0 ) = x + rhs * dt_by_rho_cp; + double dt_by_rho_cp = + ( temp_local >= solidus_ && temp_local <= liquidus_ ) + ? dt_ / ( rho_cp_ + rho_Lf_by_dT_ ) + : dt_ / ( rho_cp_ ); + + const double laplacian_x = laplacian_k( + temp_local, temp_px, temp_nx, kappa_local, kappa_px, kappa_nx ); + const double laplacian_y = laplacian_k( + temp_local, temp_py, temp_ny, kappa_local, kappa_py, kappa_ny ); + const double laplacian_z = laplacian_k( + temp_local, temp_pz, temp_nz, kappa_local, kappa_pz, kappa_nz ); + + const double rhs = + ( laplacian_x + laplacian_y + laplacian_z ) * inv_dx2_ + + source( tag, i, j, k ); + + T_( i, j, k, 0 ) = temp_local + rhs * dt_by_rho_cp; } // Device tagged version of the temperature solver @@ -152,26 +183,68 @@ class Solver void operator()( DeviceTag tag, const int i, const int j, const int k ) const { - double x = T0_( i, j, k, 0 ); + // First nearest neighbor stencil for cell at i,j,k: temperature and + // thermal conductivity + const double temp_local = T0_( i, j, k, 0 ); + const double temp_px = T0_( i + 1, j, k, 0 ); + const double temp_nx = T0_( i - 1, j, k, 0 ); + const double temp_py = T0_( i, j + 1, k, 0 ); + const double temp_ny = T0_( i, j - 1, k, 0 ); + const double temp_pz = T0_( i, j, k + 1, 0 ); + const double temp_nz = T0_( i, j, k - 1, 0 ); + + const double kappa_local = kappa_of_temperature( temp_local ); + const double kappa_px = kappa_of_temperature( temp_px ); + const double kappa_nx = kappa_of_temperature( temp_nx ); + const double kappa_py = kappa_of_temperature( temp_py ); + const double kappa_ny = kappa_of_temperature( temp_ny ); + const double kappa_pz = kappa_of_temperature( temp_pz ); + const double kappa_nz = kappa_of_temperature( temp_nz ); double dt_by_rho_cp = - dt_ / ( rho_cp_ + - ( x >= solidus_ ) * ( x <= liquidus_ ) * rho_Lf_by_dT_ ); + dt_ / ( rho_cp_ + ( temp_local >= solidus_ ) * + ( temp_local <= liquidus_ ) * rho_Lf_by_dT_ ); - double rhs = laplacian( i, j, k ) + source( tag, i, j, k ); + const double laplacian_x = laplacian_k( + temp_local, temp_px, temp_nx, kappa_local, kappa_px, kappa_nx ); + const double laplacian_y = laplacian_k( + temp_local, temp_py, temp_ny, kappa_local, kappa_py, kappa_ny ); + const double laplacian_z = laplacian_k( + temp_local, temp_pz, temp_nz, kappa_local, kappa_pz, kappa_nz ); + + const double rhs = + ( laplacian_x + laplacian_y + laplacian_z ) * inv_dx2_ + + source( tag, i, j, k ); + + T_( i, j, k, 0 ) = temp_local + rhs * dt_by_rho_cp; + } - T_( i, j, k, 0 ) = x + rhs * dt_by_rho_cp; + // Get temperature-dependent thermal conductivity + KOKKOS_INLINE_FUNCTION + auto kappa_of_temperature( const double temp_local ) const + { + return k_0_ + k_T_ * temp_local; + } + + // Get harmonic average of two values + KOKKOS_INLINE_FUNCTION + auto harmonic_mean( const double x1, const double x2 ) const + { + return 2.0 * x1 * x2 / ( x1 + x2 ); } - // First-order centered space laplacian stencil + // First-order centered space laplacian stencil for one direction - + // temperature-dependent thermal conductivity KOKKOS_INLINE_FUNCTION - auto laplacian( const int i, const int j, const int k ) const + auto laplacian_k( const double temp_local, const double temp_positive, + const double temp_negative, const double kappa_local, + const double kappa_positive, + const double kappa_negative ) const { - return ( T0_( i - 1, j, k, 0 ) + T0_( i + 1, j, k, 0 ) + - T0_( i, j - 1, k, 0 ) + T0_( i, j + 1, k, 0 ) + - T0_( i, j, k - 1, 0 ) + T0_( i, j, k + 1, 0 ) - - 6.0 * T0_( i, j, k, 0 ) ) * - k_by_dx2_; + return ( harmonic_mean( kappa_local, kappa_positive ) * + ( temp_positive - temp_local ) - + harmonic_mean( kappa_local, kappa_negative ) * + ( temp_local - temp_negative ) ); } // Normalized weight for the gaussian source term: x in exp(-x) From b7f3dc66fa44399466344197476f1f9dca6f8f90 Mon Sep 17 00:00:00 2001 From: MattRolchigo <65y@ornl.gov> Date: Thu, 19 Dec 2024 17:55:25 -0500 Subject: [PATCH 2/2] Separate conductivity for liquid and solid phases --- src/Finch_Inputs.hpp | 58 ++++++++++++++++++++---------- src/Finch_Solver.hpp | 86 +++++++++++++++++++++++++++++++------------- 2 files changed, 101 insertions(+), 43 deletions(-) diff --git a/src/Finch_Inputs.hpp b/src/Finch_Inputs.hpp index a50f760..92971a8 100644 --- a/src/Finch_Inputs.hpp +++ b/src/Finch_Inputs.hpp @@ -91,9 +91,11 @@ struct Properties { double density; double specific_heat; - double thermal_conductivity_0; - double thermal_conductivity_T; - double thermal_diffusivity; + double thermal_conductivity_solid_0; + double thermal_conductivity_solid_T; + double thermal_conductivity_liquid_0; + double thermal_conductivity_liquid_T; + double vaporization_temperature; double latent_heat; double solidus; double liquidus; @@ -212,9 +214,12 @@ class Inputs Info << "Properties:" << std::endl; Info << " Density: " << properties.density << std::endl; Info << " Specific Heat: " << properties.specific_heat << std::endl; - Info << " Thermal Conductivity: " << properties.thermal_conductivity_0 - << " + " << properties.thermal_conductivity_T << " * T" - << std::endl; + Info << " Thermal Conductivity, Solid: " + << properties.thermal_conductivity_solid_0 << " + " + << properties.thermal_conductivity_solid_T << " * T" << std::endl; + Info << " Thermal Conductivity, Liquid: " + << properties.thermal_conductivity_liquid_0 << " + " + << properties.thermal_conductivity_liquid_T << " * T" << std::endl; Info << " Latent Heat: " << properties.latent_heat << std::endl; Info << " Solidus: " << properties.solidus << std::endl; Info << " Liquidus: " << properties.liquidus << std::endl; @@ -266,8 +271,7 @@ class Inputs return filename_s; } - void parseInputFile( MPI_Comm comm, const std::string filename, - const double est_max_temperature = 5000.0 ) + void parseInputFile( MPI_Comm comm, const std::string filename ) { readInput( filename ); @@ -276,8 +280,9 @@ class Inputs // create auxiliary properties - estimate time step based on a maximum // temperature and the parsed thermal conducitivity values const double est_max_thermal_diffusivity = - ( properties.thermal_conductivity_0 + - properties.thermal_conductivity_T * est_max_temperature ) / + ( properties.thermal_conductivity_liquid_0 + + properties.thermal_conductivity_liquid_T * + properties.vaporization_temperature ) / ( properties.density * properties.specific_heat ); time.time_step = ( time.Co * space.cell_size * space.cell_size ) / @@ -341,20 +346,35 @@ class Inputs properties.density = db["properties"]["density"]; properties.specific_heat = db["properties"]["specific_heat"]; // Thermal conductvity as a function of temperature in the form a + b * - // Temperature. If only one value is given, a fixed thermal conductivity - // is used - if ( db["properties"]["thermal_conductivity"].size() > 1 ) + // Temperature for solid and liquid. If only one value is given, fixed + // thermal conductivity for both phases are used + if ( db["properties"].contains( "thermal_conductivity_solid" ) && + db["properties"].contains( "thermal_conductivity_liquid" ) ) { - properties.thermal_conductivity_0 = - db["properties"]["thermal_conductivity"][0]; - properties.thermal_conductivity_T = - db["properties"]["thermal_conductivity"][1]; + properties.thermal_conductivity_solid_0 = + db["properties"]["thermal_conductivity_solid"][0]; + properties.thermal_conductivity_solid_T = + db["properties"]["thermal_conductivity_solid"][1]; + properties.thermal_conductivity_liquid_0 = + db["properties"]["thermal_conductivity_liquid"][0]; + properties.thermal_conductivity_liquid_T = + db["properties"]["thermal_conductivity_liquid"][1]; + // Also expects a max temperature (vaporization temperature) for + // thermal conductivity calculations + properties.vaporization_temperature = + db["properties"]["vaporization_temperature"]; } else { - properties.thermal_conductivity_0 = + properties.thermal_conductivity_solid_0 = + db["properties"]["thermal_conductivity"]; + properties.thermal_conductivity_liquid_0 = db["properties"]["thermal_conductivity"]; - properties.thermal_conductivity_T = 0.0; + properties.thermal_conductivity_solid_T = 0.0; + properties.thermal_conductivity_liquid_T = 0.0; + // Default value for a max temperature (unused since + // thermal_conductivity_T = 0) + properties.vaporization_temperature = 5000.0; } properties.latent_heat = db["properties"]["latent_heat"]; properties.solidus = db["properties"]["solidus"]; diff --git a/src/Finch_Solver.hpp b/src/Finch_Solver.hpp index e8d90c0..d0457fe 100644 --- a/src/Finch_Solver.hpp +++ b/src/Finch_Solver.hpp @@ -44,11 +44,15 @@ class Solver double dt_; double solidus_; double liquidus_; + double inv_freezing_range_; double rho_cp_; double rho_Lf_by_dT_; double inv_dx2_; - double k_0_; - double k_T_; + double k_solid_0_; + double k_solid_T_; + double k_liquid_0_; + double k_liquid_T_; + double temp_max_; // heat source parameters double power_; @@ -75,14 +79,20 @@ class Solver liquidus_ = db.properties.liquidus; + inv_freezing_range_ = 1.0 / ( liquidus_ - solidus_ ); + + temp_max_ = db.properties.vaporization_temperature; + rho_cp_ = rho * cp; rho_Lf_by_dT_ = rho * Lf / ( liquidus_ - solidus_ ); inv_dx2_ = 1.0 / ( dx * dx ); - k_0_ = db.properties.thermal_conductivity_0; - k_T_ = db.properties.thermal_conductivity_T; + k_solid_0_ = db.properties.thermal_conductivity_solid_0; + k_solid_T_ = db.properties.thermal_conductivity_solid_T; + k_liquid_0_ = db.properties.thermal_conductivity_liquid_0; + k_liquid_T_ = db.properties.thermal_conductivity_liquid_T; // initialize beam position for ( std::size_t d = 0; d < 3; ++d ) @@ -151,16 +161,27 @@ class Solver const double temp_pz = T0_( i, j, k + 1, 0 ); const double temp_nz = T0_( i, j, k - 1, 0 ); - const double kappa_local = kappa_of_temperature( temp_local ); - const double kappa_px = kappa_of_temperature( temp_px ); - const double kappa_nx = kappa_of_temperature( temp_nx ); - const double kappa_py = kappa_of_temperature( temp_py ); - const double kappa_ny = kappa_of_temperature( temp_ny ); - const double kappa_pz = kappa_of_temperature( temp_pz ); - const double kappa_nz = kappa_of_temperature( temp_nz ); + // Liquid fraction + const double liquid_fraction = + Kokkos::fmax( 1.0, Kokkos::fmin( 0.0, ( temp_local - solidus_ ) * + inv_freezing_range_ ) ); + const double kappa_local = + kappa_of_temperature( temp_local, liquid_fraction ); + const double kappa_px = + kappa_of_temperature( temp_px, liquid_fraction ); + const double kappa_nx = + kappa_of_temperature( temp_nx, liquid_fraction ); + const double kappa_py = + kappa_of_temperature( temp_py, liquid_fraction ); + const double kappa_ny = + kappa_of_temperature( temp_ny, liquid_fraction ); + const double kappa_pz = + kappa_of_temperature( temp_pz, liquid_fraction ); + const double kappa_nz = + kappa_of_temperature( temp_nz, liquid_fraction ); double dt_by_rho_cp = - ( temp_local >= solidus_ && temp_local <= liquidus_ ) + ( liquid_fraction >= 0.0 && liquid_fraction <= 1.0 ) ? dt_ / ( rho_cp_ + rho_Lf_by_dT_ ) : dt_ / ( rho_cp_ ); @@ -193,17 +214,28 @@ class Solver const double temp_pz = T0_( i, j, k + 1, 0 ); const double temp_nz = T0_( i, j, k - 1, 0 ); - const double kappa_local = kappa_of_temperature( temp_local ); - const double kappa_px = kappa_of_temperature( temp_px ); - const double kappa_nx = kappa_of_temperature( temp_nx ); - const double kappa_py = kappa_of_temperature( temp_py ); - const double kappa_ny = kappa_of_temperature( temp_ny ); - const double kappa_pz = kappa_of_temperature( temp_pz ); - const double kappa_nz = kappa_of_temperature( temp_nz ); + // Liquid fraction + const double liquid_fraction = + Kokkos::fmax( 1.0, Kokkos::fmin( 0.0, ( temp_local - solidus_ ) * + inv_freezing_range_ ) ); + const double kappa_local = + kappa_of_temperature( temp_local, liquid_fraction ); + const double kappa_px = + kappa_of_temperature( temp_px, liquid_fraction ); + const double kappa_nx = + kappa_of_temperature( temp_nx, liquid_fraction ); + const double kappa_py = + kappa_of_temperature( temp_py, liquid_fraction ); + const double kappa_ny = + kappa_of_temperature( temp_ny, liquid_fraction ); + const double kappa_pz = + kappa_of_temperature( temp_pz, liquid_fraction ); + const double kappa_nz = + kappa_of_temperature( temp_nz, liquid_fraction ); double dt_by_rho_cp = - dt_ / ( rho_cp_ + ( temp_local >= solidus_ ) * - ( temp_local <= liquidus_ ) * rho_Lf_by_dT_ ); + dt_ / ( rho_cp_ + ( liquid_fraction >= 0.0 ) * + ( liquid_fraction <= 1.0 ) * rho_Lf_by_dT_ ); const double laplacian_x = laplacian_k( temp_local, temp_px, temp_nx, kappa_local, kappa_px, kappa_nx ); @@ -219,11 +251,17 @@ class Solver T_( i, j, k, 0 ) = temp_local + rhs * dt_by_rho_cp; } - // Get temperature-dependent thermal conductivity + // Get temperature-dependent thermal conductivity, using liquid and solid + // values. Conductivity is capped at the value where temp_local = temp_max_ KOKKOS_INLINE_FUNCTION - auto kappa_of_temperature( const double temp_local ) const + auto kappa_of_temperature( const double temp_local, + const double liquid_fraction ) const { - return k_0_ + k_T_ * temp_local; + return ( 1.0 - liquid_fraction ) * + ( k_solid_0_ + k_solid_T_ * temp_local ) + + liquid_fraction * + ( k_liquid_0_ + + k_liquid_T_ * Kokkos::fmin( temp_local, temp_max_ ) ); } // Get harmonic average of two values