From a6a97d6c8f5c262936f49640febc1d4f4efc645c Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 18 Jun 2024 13:12:01 -0400 Subject: [PATCH 1/9] Add heat transfer --- .../inputs/thermal_deformation.json | 24 ++-- src/CabanaPD_ForceModels.hpp | 40 ++++++ src/CabanaPD_HeatTransfer.hpp | 127 +++++++++++++++++ src/CabanaPD_Particles.hpp | 23 ++- src/CabanaPD_Solver.hpp | 18 ++- src/CabanaPD_Types.hpp | 10 +- src/force/CabanaPD_ForceModels_PMB.hpp | 132 ++++++++++++++++++ 7 files changed, 359 insertions(+), 15 deletions(-) create mode 100644 src/CabanaPD_HeatTransfer.hpp diff --git a/examples/thermomechanics/inputs/thermal_deformation.json b/examples/thermomechanics/inputs/thermal_deformation.json index 4f0f6908..48af0aa8 100644 --- a/examples/thermomechanics/inputs/thermal_deformation.json +++ b/examples/thermomechanics/inputs/thermal_deformation.json @@ -1,13 +1,15 @@ { - "num_cells" : {"value": [100, 30, 3]}, - "system_size" : {"value": [1.0, 0.3, 0.03], "unit": "m"}, - "density" : {"value": 3980, "unit": "kg/m^3"}, - "elastic_modulus" : {"value": 370e+9, "unit": "Pa"}, - "thermal_coefficient" : {"value": 7.5E-6, "unit": "oC^{-1}"}, - "reference_temperature" : {"value": 20.0, "unit": "oC"}, - "horizon" : {"value": 0.03, "unit": "m"}, - "final_time" : {"value": 0.01, "unit": "s"}, - "timestep" : {"value": 7.5E-7, "unit": "s"}, - "output_frequency" : {"value": 100}, - "output_reference" : {"value": true} + "num_cells" : {"value": [100, 30, 3]}, + "system_size" : {"value": [1.0, 0.3, 0.03], "unit": "m"}, + "density" : {"value": 3980, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 370e+9, "unit": "Pa"}, + "thermal_coefficient" : {"value": 7.5E-6, "unit": "oC^{-1}"}, + "thermal_conductivity" : {"value": 31, "unit": "W/(m.K)"}, + "specific_heat_capacity" : {"value": 880, "unit": "J/(kg.K)"}, + "reference_temperature" : {"value": 20.0, "unit": "oC"}, + "horizon" : {"value": 0.03, "unit": "m"}, + "final_time" : {"value": 0.01, "unit": "s"}, + "timestep" : {"value": 7.5E-7, "unit": "s"}, + "output_frequency" : {"value": 100}, + "output_reference" : {"value": true} } diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index 6685de1c..ea63c8fe 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -12,6 +12,7 @@ #ifndef FORCE_MODELS_H #define FORCE_MODELS_H +#include #include namespace CabanaPD @@ -71,6 +72,45 @@ struct BaseTemperatureModel } }; +// This class stores temperature parameters needed for heat transfer, but not +// the temperature itself (stored instead in the static temperature class +// above). +struct BaseDynamicTemperatureModel +{ + double delta; + + double thermal_coeff; + double kappa; + double cp; + bool constant_microconductivity; + + BaseDynamicTemperatureModel( const double _delta, const double _kappa, + const double _cp, + const bool _constant_microconductivity = true ) + { + set_param( _delta, _kappa, _cp, _constant_microconductivity ); + } + + void set_param( const double _delta, const double _kappa, const double _cp, + const bool _constant_microconductivity ) + { + delta = _delta; + kappa = _kappa; + cp = _cp; + const double d3 = _delta * _delta * _delta; + thermal_coeff = 9.0 / 2.0 * _kappa / pi / d3; + constant_microconductivity = _constant_microconductivity; + } + + KOKKOS_INLINE_FUNCTION double microconductivity_function( double r ) const + { + if ( constant_microconductivity ) + return thermal_coeff; + else + return 4.0 * thermal_coeff * ( 1.0 - r / delta ); + } +}; + template struct ForceModel; diff --git a/src/CabanaPD_HeatTransfer.hpp b/src/CabanaPD_HeatTransfer.hpp new file mode 100644 index 00000000..75a7531c --- /dev/null +++ b/src/CabanaPD_HeatTransfer.hpp @@ -0,0 +1,127 @@ +/**************************************************************************** + * Copyright (c) 2022-2023 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of CabanaPD. CabanaPD 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 * + ****************************************************************************/ + +#ifndef HEATTRANSFER_H +#define HEATTRANSFER_H + +#include + +namespace CabanaPD +{ + +// Peridynamic heat transfer with forward-Euler time integration. +// Inherits only because this is a similar neighbor-based kernel. +template +class HeatTransfer : public Force +{ + protected: + using base_type = Force; + using base_type::_half_neigh; + using base_type::_timer; + using model_type = ModelType; + static_assert( + std::is_same_v ); + + Timer _euler_timer = base_type::_energy_timer; + ModelType _model; + + public: + HeatTransfer( const bool half_neigh, const ModelType model ) + : base_type( half_neigh ) + , _model( model ) + { + } + + template + void + computeHeatTransferFull( TemperatureType& conduction, const PosType& x, + const PosType& u, const ParticleType& particles, + const NeighListType& neigh_list, const int n_local, + ParallelType& neigh_op_tag ) + { + _timer.start(); + + auto model = _model; + const auto vol = particles.sliceVolume(); + const auto temp = particles.sliceTemperature(); + + auto temp_func = KOKKOS_LAMBDA( const int i, const int j ) + { + double xi, r, s; + getDistance( x, u, i, j, xi, r, s ); + + const double coeff = model.microconductivity_function( xi ); + conduction( i ) += + coeff * ( temp( j ) - temp( i ) ) / xi / xi * vol( j ); + }; + + Kokkos::RangePolicy policy( 0, n_local ); + Cabana::neighbor_parallel_for( + policy, temp_func, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, "CabanaPD::HeatTransfer::computeFull" ); + + _timer.stop(); + } + + template + void forwardEuler( const ParticleType& particles, const double dt ) + { + _euler_timer.start(); + auto model = _model; + const auto rho = particles.sliceDensity(); + const auto conduction = particles.sliceTemperatureConduction(); + auto temp = particles.sliceTemperature(); + auto n_local = particles.n_local; + auto euler_func = KOKKOS_LAMBDA( const int i ) + { + temp( i ) += dt / rho( i ) / model.cp * conduction( i ); + }; + Kokkos::RangePolicy policy( 0, n_local ); + Kokkos::parallel_for( "CabanaPD::HeatTransfer::forwardEuler", policy, + euler_func ); + _euler_timer.stop(); + } +}; + +// Heat transfer free function. +template +void computeHeatTransfer( HeatTransferType& heat_transfer, + ParticleType& particles, + const NeighListType& neigh_list, + const ParallelType& neigh_op_tag, const double dt ) +{ + auto n_local = particles.n_local; + auto x = particles.sliceReferencePosition(); + auto u = particles.sliceDisplacement(); + auto conduction = particles.sliceTemperatureConduction(); + auto conduction_a = particles.sliceTemperatureConductionAtomic(); + + // Reset temperature conduction. + Cabana::deep_copy( conduction, 0.0 ); + + // Temperature only needs to be atomic if using team threading. + if ( std::is_same::value ) + heat_transfer.computeHeatTransferFull( + conduction_a, x, u, particles, neigh_list, n_local, neigh_op_tag ); + else + heat_transfer.computeHeatTransferFull( + conduction, x, u, particles, neigh_list, n_local, neigh_op_tag ); + Kokkos::fence(); + + heat_transfer.forwardEuler( particles, dt ); + Kokkos::fence(); +} + +} // namespace CabanaPD + +#endif diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index c1bdb472..82e8f693 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -88,6 +88,7 @@ class Particles public: using self_type = Particles; + using thermal_type = TemperatureIndependent; using memory_space = MemorySpace; using execution_space = typename memory_space::execution_space; static constexpr int dim = Dimension; @@ -502,6 +503,7 @@ class Particles Particles; using base_type = Particles; + using thermal_type = TemperatureIndependent; using memory_space = typename base_type::memory_space; using base_type::dim; @@ -662,6 +664,7 @@ class Particles Particles; using base_type = Particles; + using thermal_type = TemperatureDependent; using memory_space = typename base_type::memory_space; using base_type::dim; @@ -673,8 +676,8 @@ class Particles // These are split since weighted volume only needs to be communicated once // and dilatation only needs to be communicated for LPS. - using scalar_type = typename base_type::scalar_type; - using aosoa_temp_type = Cabana::AoSoA; + using temp_types = Cabana::MemberTypes; + using aosoa_temp_type = Cabana::AoSoA; // Per type. using base_type::n_types; @@ -729,6 +732,22 @@ class Particles { return Cabana::slice<0>( _aosoa_temp, "temperature" ); } + auto sliceTemperatureConduction() + { + return Cabana::slice<1>( _aosoa_temp, "temperature_conduction" ); + } + auto sliceTemperatureConduction() const + { + return Cabana::slice<1>( _aosoa_temp, "temperature_conduction" ); + } + auto sliceTemperatureConductionAtomic() + { + auto temp = sliceTemperature(); + using slice_type = decltype( temp ); + using atomic_type = typename slice_type::atomic_access_slice; + atomic_type temp_a = temp; + return temp_a; + } void resize( int new_local, int new_ghost ) { diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 04329965..55073ee2 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -74,6 +74,7 @@ #include #include #include +#include #include #include #include @@ -103,13 +104,16 @@ class SolverElastic using force_model_type = ForceModel; using force_type = Force; using comm_type = Comm; + typename particle_type::thermal_type>; using neighbor_type = Cabana::VerletList; using neigh_iter_tag = Cabana::SerialOpTag; using input_type = InputType; + // Optional modules. + using heat_transfer_type = HeatTransfer; + SolverElastic( input_type _inputs, std::shared_ptr _particles, force_model_type force_model ) @@ -133,6 +137,12 @@ class SolverElastic TemperatureDependent>::value ) force_model.update( particles->sliceTemperature() ); + // Create heat transfer if needed. + if constexpr ( std::is_same::value ) + heat_transfer = std::make_shared( + inputs["half_neigh"], force_model ); + force = std::make_shared( inputs["half_neigh"], force_model ); @@ -265,6 +275,10 @@ class SolverElastic if constexpr ( std::is_same::value ) comm->gatherTemperature(); + if constexpr ( std::is_same::value ) + computeHeatTransfer( *heat_transfer, *particles, *neighbors, + neigh_iter_tag{}, dt ); // Update ghost particles. comm->gatherDisplacement(); @@ -437,6 +451,8 @@ class SolverElastic std::shared_ptr integrator; std::shared_ptr force; std::shared_ptr neighbors; + // Optional. + std::shared_ptr heat_transfer; std::string output_file; std::string error_file; diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index 62076735..b1fa4399 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -14,19 +14,27 @@ namespace CabanaPD { +// Mechanics types. struct Elastic { }; struct Fracture { }; + +// Thermal types. +struct TemperatureIndependent +{ +}; struct TemperatureDependent { }; -struct TemperatureIndependent +struct DynamicTemperature : public TemperatureDependent { + using base_type = TemperatureDependent; }; +// Model types. struct PMB { }; diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp index c80acd63..f0200ec3 100644 --- a/src/force/CabanaPD_ForceModels_PMB.hpp +++ b/src/force/CabanaPD_ForceModels_PMB.hpp @@ -252,6 +252,138 @@ auto createForceModel( PMB, Fracture, ParticleType particles, delta, K, G0, temp, alpha, temp0 ); } +template +struct ForceModel + : public ForceModel, + BaseDynamicTemperatureModel +{ + using base_type = + ForceModel; + using base_temperature_type = BaseDynamicTemperatureModel; + using base_model = PMB; + using fracture_type = Elastic; + using thermal_type = DynamicTemperature; + + using base_type::c; + using base_type::delta; + using base_type::K; + + // Thermal parameters + using base_temperature_type::cp; + using base_temperature_type::kappa; + using base_temperature_type::thermal_coeff; + using base_type::alpha; + using base_type::temp0; + using base_type::temperature; + + // Explicitly use the temperature-dependent stretch. + using base_type::thermalStretch; + + // ForceModel(){}; + ForceModel( const double _delta, const double _K, + const TemperatureType _temp, const double _kappa, + const double _cp, const double _alpha, + const double _temp0 = 0.0, + const bool _constant_microconductivity = true ) + : base_type( _delta, _K, _temp, _alpha, _temp0 ) + , base_temperature_type( _delta, _kappa, _cp, + _constant_microconductivity ) + { + set_param( _delta, _K, _kappa, _cp, _alpha, _temp0, + _constant_microconductivity ); + } + + void set_param( const double _delta, const double _K, const double _kappa, + const double _cp, const double _alpha, const double _temp0, + const bool _constant_microconductivity ) + { + base_type::set_param( _delta, _K, _alpha, _temp0 ); + base_temperature_type::set_param( _delta, _kappa, _cp, + _constant_microconductivity ); + } +}; + +template +auto createForceModel( PMB, Elastic, ParticleType particles, const double delta, + const double K, const double kappa, const double cp, + const double alpha, const double temp0, + const bool constant_microconductivity = true ) +{ + auto temp = particles.sliceTemperature(); + using temp_type = decltype( temp ); + return ForceModel( + delta, K, temp, kappa, cp, alpha, temp0, constant_microconductivity ); +} + +template +struct ForceModel + : public ForceModel, + BaseDynamicTemperatureModel +{ + using base_type = + ForceModel; + using base_temperature_type = BaseDynamicTemperatureModel; + using base_model = typename base_type::base_model; + using fracture_type = typename base_type::fracture_type; + using thermal_type = DynamicTemperature; + + using base_type::c; + using base_type::delta; + using base_type::K; + + // Does not use the base bond_break_coeff. + using base_type::G0; + using base_type::s0; + + // Thermal parameters + using base_temperature_type::cp; + using base_temperature_type::kappa; + using base_temperature_type::thermal_coeff; + using base_type::alpha; + using base_type::temp0; + using base_type::temperature; + + // Explicitly use the temperature-dependent stretch. + using base_type::thermalStretch; + + // ForceModel(){}; + ForceModel( const double _delta, const double _K, const double _G0, + const TemperatureType _temp, const double _kappa, + const double _cp, const double _alpha, + const double _temp0 = 0.0, + const bool _constant_microconductivity = true ) + : base_type( _delta, _K, _G0, _temp, _alpha, _temp0 ) + { + set_param( _delta, _K, _G0, _kappa, _cp, _alpha, _temp0 ); + base_temperature_type::set_param( _delta, _kappa, _cp, + _constant_microconductivity ); + } + + void set_param( const double _delta, const double _K, const double _G0, + const double _kappa, const double _cp, const double _alpha, + const double _temp0, + const bool _constant_microconductivity ) + { + base_type::set_param( _delta, _K, _G0, _alpha, _temp0 ); + base_temperature_type::set_param( _delta, _kappa, _cp, + _constant_microconductivity ); + } +}; + +template +auto createForceModel( PMB, Fracture, ParticleType particles, + const double delta, const double K, const double G0, + const double kappa, const double cp, const double alpha, + const double temp0, + const bool constant_microconductivity = true ) +{ + auto temp = particles.sliceTemperature(); + using temp_type = decltype( temp ); + return ForceModel( + delta, K, G0, temp, kappa, cp, alpha, temp0, + constant_microconductivity ); +} + } // namespace CabanaPD #endif From cd00c096f92f41b63bedc548b8e9b147cc5c8630 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 3 Jul 2024 13:02:19 -0400 Subject: [PATCH 2/9] fixup: solver comments --- src/CabanaPD_Solver.hpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 55073ee2..a9bd4774 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -99,6 +99,7 @@ class SolverElastic using memory_space = MemorySpace; using exec_space = typename memory_space::execution_space; + // Core module types - required for all problems. using particle_type = ParticleType; using integrator_type = Integrator; using force_model_type = ForceModel; @@ -111,7 +112,7 @@ class SolverElastic using neigh_iter_tag = Cabana::SerialOpTag; using input_type = InputType; - // Optional modules. + // Optional module types. using heat_transfer_type = HeatTransfer; SolverElastic( input_type _inputs, @@ -445,19 +446,21 @@ class SolverElastic double dt; protected: + // Core modules. input_type inputs; std::shared_ptr particles; std::shared_ptr comm; std::shared_ptr integrator; std::shared_ptr force; std::shared_ptr neighbors; - // Optional. + // Optional modules. std::shared_ptr heat_transfer; + // Output files. std::string output_file; std::string error_file; - // Combined from many class timers. + // Note: init_time is combined from many class timers. double _init_time; Timer _init_timer; Timer _neighbor_timer; From 152bb2757da5b19ee47e75c6125ca426b528c183 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 31 Jul 2024 13:40:53 -0400 Subject: [PATCH 3/9] Add new heat transfer example Includes various b.c. options Modified b.c. on thermal deformation heat transfer example Changed example output profile dimension --- examples/thermomechanics/CMakeLists.txt | 5 +- .../inputs/thermal_deformation.json | 26 +- .../inputs/thermal_deformation_cube.json | 15 ++ .../inputs/thermal_deformation_cube_temp.json | 15 ++ .../thermomechanics/thermal_deformation.cpp | 2 +- .../thermal_deformation_heat_transfer.cpp | 244 ++++++++++++++++++ 6 files changed, 292 insertions(+), 15 deletions(-) create mode 100644 examples/thermomechanics/inputs/thermal_deformation_cube.json create mode 100644 examples/thermomechanics/inputs/thermal_deformation_cube_temp.json create mode 100644 examples/thermomechanics/thermal_deformation_heat_transfer.cpp diff --git a/examples/thermomechanics/CMakeLists.txt b/examples/thermomechanics/CMakeLists.txt index 73a932b3..f08ec6e1 100644 --- a/examples/thermomechanics/CMakeLists.txt +++ b/examples/thermomechanics/CMakeLists.txt @@ -4,4 +4,7 @@ target_link_libraries(ThermalDeformation LINK_PUBLIC CabanaPD) add_executable(ThermalCrack thermal_crack.cpp) target_link_libraries(ThermalCrack LINK_PUBLIC CabanaPD) -install(TARGETS ThermalDeformation ThermalCrack DESTINATION ${CMAKE_INSTALL_BINDIR}) +add_executable(ThermalDeformationHeatTransfer thermal_deformation_heat_transfer.cpp) +target_link_libraries(ThermalDeformationHeatTransfer LINK_PUBLIC CabanaPD) + +install(TARGETS ThermalDeformation ThermalCrack ThermalDeformationHeatTransfer DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/examples/thermomechanics/inputs/thermal_deformation.json b/examples/thermomechanics/inputs/thermal_deformation.json index 48af0aa8..83c8874e 100644 --- a/examples/thermomechanics/inputs/thermal_deformation.json +++ b/examples/thermomechanics/inputs/thermal_deformation.json @@ -1,15 +1,15 @@ { - "num_cells" : {"value": [100, 30, 3]}, - "system_size" : {"value": [1.0, 0.3, 0.03], "unit": "m"}, - "density" : {"value": 3980, "unit": "kg/m^3"}, - "elastic_modulus" : {"value": 370e+9, "unit": "Pa"}, - "thermal_coefficient" : {"value": 7.5E-6, "unit": "oC^{-1}"}, - "thermal_conductivity" : {"value": 31, "unit": "W/(m.K)"}, - "specific_heat_capacity" : {"value": 880, "unit": "J/(kg.K)"}, - "reference_temperature" : {"value": 20.0, "unit": "oC"}, - "horizon" : {"value": 0.03, "unit": "m"}, - "final_time" : {"value": 0.01, "unit": "s"}, - "timestep" : {"value": 7.5E-7, "unit": "s"}, - "output_frequency" : {"value": 100}, - "output_reference" : {"value": true} + "num_cells" : {"value": [101, 31, 3]}, + "system_size" : {"value": [1.0, 0.3, 0.03], "unit": "m"}, + "density" : {"value": 3980, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 370e+9, "unit": "Pa"}, + "thermal_expansion_coeff" : {"value": 7.5E-6, "unit": "oC^{-1}"}, + "thermal_conductivity" : {"value": 31, "unit": "W/(m.K)"}, + "specific_heat_capacity" : {"value": 880, "unit": "J/(kg.K)"}, + "reference_temperature" : {"value": 20.0, "unit": "oC"}, + "horizon" : {"value": 0.03, "unit": "m"}, + "final_time" : {"value": 0.01, "unit": "s"}, + "timestep" : {"value": 7.5E-7, "unit": "s"}, + "output_frequency" : {"value": 100}, + "output_reference" : {"value": true} } diff --git a/examples/thermomechanics/inputs/thermal_deformation_cube.json b/examples/thermomechanics/inputs/thermal_deformation_cube.json new file mode 100644 index 00000000..b7e054ec --- /dev/null +++ b/examples/thermomechanics/inputs/thermal_deformation_cube.json @@ -0,0 +1,15 @@ +{ + "num_cells" : {"value": [40, 40, 40]}, + "system_size" : {"value": [0.1, 0.1, 0.1], "unit": "m"}, + "density" : {"value": 8915, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 115e+9, "unit": "Pa"}, + "thermal_expansion_coeff" : {"value": 17e-6, "unit": "oC^{-1}"}, + "thermal_conductivity" : {"value": 387, "unit": "W/(m.K)"}, + "specific_heat_capacity" : {"value": 385, "unit": "J/(kg.K)"}, + "reference_temperature" : {"value": 100.0, "unit": "oC"}, + "horizon" : {"value": 0.0075, "unit": "m"}, + "final_time" : {"value": 8, "unit": "s"}, + "timestep" : {"value": 5.1e-8, "unit": "s"}, + "output_frequency" : {"value": 1}, + "output_reference" : {"value": true} +} diff --git a/examples/thermomechanics/inputs/thermal_deformation_cube_temp.json b/examples/thermomechanics/inputs/thermal_deformation_cube_temp.json new file mode 100644 index 00000000..163195b0 --- /dev/null +++ b/examples/thermomechanics/inputs/thermal_deformation_cube_temp.json @@ -0,0 +1,15 @@ +{ + "num_cells" : {"value": [51, 51, 51]}, + "system_size" : {"value": [0.5, 0.5, 0.5], "unit": "m"}, + "density" : {"value": 3980, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 370e+9, "unit": "Pa"}, + "thermal_expansion_coeff" : {"value": 7.5E-6, "unit": "oC^{-1}"}, + "thermal_conductivity" : {"value": 31, "unit": "W/(m.K)"}, + "specific_heat_capacity" : {"value": 880, "unit": "J/(kg.K)"}, + "reference_temperature" : {"value": 100.0, "unit": "oC"}, + "horizon" : {"value": 0.03, "unit": "m"}, + "final_time" : {"value": 0.01, "unit": "s"}, + "timestep" : {"value": 1.0E-7, "unit": "s"}, + "output_frequency" : {"value": 100}, + "output_reference" : {"value": true} +} diff --git a/examples/thermomechanics/thermal_deformation.cpp b/examples/thermomechanics/thermal_deformation.cpp index 643cd21d..74491ed3 100644 --- a/examples/thermomechanics/thermal_deformation.cpp +++ b/examples/thermomechanics/thermal_deformation.cpp @@ -41,7 +41,7 @@ void thermalDeformationExample( const std::string filename ) double nu = 0.25; double K = E / ( 3 * ( 1 - 2 * nu ) ); double delta = inputs["horizon"]; - double alpha = inputs["thermal_coefficient"]; + double alpha = inputs["thermal_expansion_coeff"]; // Problem parameters double temp0 = inputs["reference_temperature"]; diff --git a/examples/thermomechanics/thermal_deformation_heat_transfer.cpp b/examples/thermomechanics/thermal_deformation_heat_transfer.cpp new file mode 100644 index 00000000..1259beba --- /dev/null +++ b/examples/thermomechanics/thermal_deformation_heat_transfer.cpp @@ -0,0 +1,244 @@ +/**************************************************************************** + * Copyright (c) 2022 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of CabanaPD. CabanaPD 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 * + ****************************************************************************/ + +#include +#include + +#include "mpi.h" + +#include + +#include + +// Simulate thermally-induced deformation in a rectangular plate. +void thermalDeformationHeatTransferExample( const std::string filename ) +{ + // ==================================================== + // Use default Kokkos spaces + // ==================================================== + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = typename exec_space::memory_space; + + // ==================================================== + // Read inputs + // ==================================================== + CabanaPD::Inputs inputs( filename ); + + // ==================================================== + // Material and problem parameters + // ==================================================== + // Material parameters + double rho0 = inputs["density"]; + double E = inputs["elastic_modulus"]; + double nu = 0.25; + double K = E / ( 3 * ( 1 - 2 * nu ) ); + double delta = inputs["horizon"]; + double alpha = inputs["thermal_expansion_coeff"]; + double kappa = inputs["thermal_conductivity"]; + double cp = inputs["specific_heat_capacity"]; + + // Problem parameters + double temp0 = inputs["reference_temperature"]; + + // ==================================================== + // Discretization + // ==================================================== + // FIXME: set halo width based on delta + std::array low_corner = inputs["low_corner"]; + std::array high_corner = inputs["high_corner"]; + std::array num_cells = inputs["num_cells"]; + int m = std::floor( delta / + ( ( high_corner[0] - low_corner[0] ) / num_cells[0] ) ); + int halo_width = m + 1; // Just to be safe. + + // ==================================================== + // Force model type + // ==================================================== + using model_type = CabanaPD::PMB; + using thermal_type = CabanaPD::DynamicTemperature; + + // ==================================================== + // Particle generation + // ==================================================== + // Does not set displacements, velocities, etc. + auto particles = + std::make_shared>( + exec_space(), low_corner, high_corner, num_cells, halo_width ); + + // ==================================================== + // Custom particle initialization + // ==================================================== + auto rho = particles->sliceDensity(); + auto temp = particles->sliceTemperature(); + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + // Density + rho( pid ) = rho0; + // Temperature + temp( pid ) = temp0; + }; + particles->updateParticles( exec_space{}, init_functor ); + + // ==================================================== + // Force model + // ==================================================== + auto force_model = CabanaPD::createForceModel( + model_type{}, CabanaPD::Elastic{}, *particles, delta, K, kappa, cp, + alpha, temp0 ); + + // ==================================================== + // Create solver + // ==================================================== + auto cabana_pd = CabanaPD::createSolverElastic( + inputs, particles, force_model ); + + // ==================================================== + // Boundary condition + // ==================================================== + // EXAMPLE 1: Temperature profile imposed over entire domain + using plane_type = CabanaPD::RegionBoundary; + /* + plane_type plane( low_corner[0], high_corner[0], + low_corner[1], high_corner[1], + low_corner[2], high_corner[2] ); + + std::vector planes = { plane }; + */ + + // EXAMPLE 2: Temperature profile imposed on top, bottom, left, and right + // surfaces + double dx = particles->dx[0]; + double dy = particles->dx[1]; + double dz = particles->dx[2]; + + // Left surface: x-direction + plane_type plane1( low_corner[0] - dx, low_corner[0] + dx, low_corner[1], + high_corner[1], low_corner[2], high_corner[2] ); + + // Right surface: x-direction + plane_type plane2( high_corner[0] - dx, high_corner[0] + dx, low_corner[1], + high_corner[1], low_corner[2], high_corner[2] ); + + // Top surface: y-direction + plane_type plane3( low_corner[0], high_corner[0], high_corner[1] - dy, + high_corner[1] + dy, low_corner[2], high_corner[2] ); + + // Bottom surface: y-direction + plane_type plane4( low_corner[0], high_corner[0], low_corner[1] - dy, + low_corner[1] + dy, low_corner[2], high_corner[2] ); + + // Front surface: z-direction + plane_type plane5( low_corner[0], high_corner[0], low_corner[1], + high_corner[1], low_corner[2] - dz, low_corner[2] + dz ); + + // Back surface: z-direction + plane_type plane6( low_corner[0], high_corner[0], low_corner[1], + high_corner[1], high_corner[2] - dz, + high_corner[2] + dz ); + + // std::vector planes = { plane1, plane2, plane3, plane4 }; + // std::vector planes = { plane1, plane2, plane3, plane4, + // plane5, plane6 }; std::vector planes = { plane1, plane2 }; + std::vector planes = { plane3, plane4 }; + /* + // EXAMPLE 3: Temperature profile imposed on top, bottom, left, and + right + // nonlocal boundaries (width delta) Top surface + plane_type plane1( + low_corner[0], high_corner[0], high_corner[1] - delta, + high_corner[1] + delta, low_corner[2], high_corner[2] ); + + // Bottom surface + plane_type plane2( + low_corner[0], high_corner[0], low_corner[1] - delta, + low_corner[1] + delta, low_corner[2], high_corner[2] ); + + // Left surface + plane_type plane3( + low_corner[0] - delta, low_corner[0] + delta, low_corner[1], + high_corner[1], low_corner[2], high_corner[2] ); + + // Right surface + plane_type plane4( + high_corner[0] - delta, high_corner[0] + delta, low_corner[1], + high_corner[1], low_corner[2], high_corner[2] ); + + std::vector planes = { plane1, plane2, plane3, + plane4 }; + */ + + auto x = particles->sliceReferencePosition(); + const double low_corner_y = low_corner[1]; + // This is purposely delayed until after solver init so that ghosted + // particles are correctly taken into account for lambda capture here. + auto temp_func = KOKKOS_LAMBDA( const int pid, const double t ) + { + // temp( pid ) = temp0 + 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; + temp( pid ) = 0.0; + }; + auto bc = CabanaPD::createBoundaryCondition( + temp_func, exec_space{}, *particles, planes, false, 1.0 ); + + // ==================================================== + // Simulation run + // ==================================================== + cabana_pd->init( bc ); + cabana_pd->run( bc ); + + // ==================================================== + // Outputs + // ==================================================== + + // Output temperature along the y-axis + auto temp_output = particles->sliceTemperature(); + auto value = KOKKOS_LAMBDA( const int pid ) { return temp_output( pid ); }; + + int profile_dim = 1; + std::string file_name = "temperature_yaxis_profile.txt"; + createOutputProfile( MPI_COMM_WORLD, num_cells[1], profile_dim, file_name, + *particles, value ); + + /* + // Output y-displacement along the x-axis + createDisplacementProfile( MPI_COMM_WORLD, + "ydisplacement_xaxis_profile.txt", *particles, + num_cells[0], 0, 1 ); + + // Output y-displacement along the y-axis + createDisplacementProfile( MPI_COMM_WORLD, + "ydisplacement_yaxis_profile.txt", *particles, + num_cells[1], 1, 1 ); + + // Output displacement magnitude along the x-axis + createDisplacementMagnitudeProfile( + MPI_COMM_WORLD, "displacement_magnitude_xaxis_profile.txt", *particles, + num_cells[0], 0 ); + + // Output displacement magnitude along the y-axis + createDisplacementMagnitudeProfile( + MPI_COMM_WORLD, "displacement_magnitude_yaxis_profile.txt", *particles, + num_cells[1], 1 ); + + */ +} + +// Initialize MPI+Kokkos. +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + thermalDeformationHeatTransferExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} From c1db9f7dd9b8dd08d07216ddbbed8e22954c7d30 Mon Sep 17 00:00:00 2001 From: pabloseleson Date: Tue, 6 Aug 2024 15:44:35 -0400 Subject: [PATCH 4/9] Rearrange order of operations in solver for heat transfer --- src/CabanaPD_Solver.hpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index a9bd4774..884958a9 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -269,20 +269,20 @@ class SolverElastic // Integrate - velocity Verlet first half. integrator->initialHalfStep( *particles ); - // Add non-force boundary condition. - if ( !boundary_condition.forceUpdate() ) - boundary_condition.apply( exec_space(), *particles, step * dt ); + // Update ghost particles. + comm->gatherDisplacement(); - if constexpr ( std::is_same::value ) - comm->gatherTemperature(); if constexpr ( std::is_same::value ) computeHeatTransfer( *heat_transfer, *particles, *neighbors, neigh_iter_tag{}, dt ); + if constexpr ( std::is_same::value ) + comm->gatherTemperature(); - // Update ghost particles. - comm->gatherDisplacement(); + // Add non-force boundary condition. + if ( !boundary_condition.forceUpdate() ) + boundary_condition.apply( exec_space(), *particles, step * dt ); // Compute internal forces. updateForce(); From 7c3fd65d44c6971d78338277f08f9b8fc0da0a3c Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 20 Aug 2024 11:37:48 -0400 Subject: [PATCH 5/9] readme: update features --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index bf9a79e1..83d31d97 100644 --- a/README.md +++ b/README.md @@ -115,6 +115,7 @@ CabanaPD currently includes the following: - Elastic (no failure) - Brittle fracture - Thermomechanics (bond-based only) + - Optional heat transfer (elastic only) - Time integration - Velocity Verlet - Pre-crack creation From 341b7d45a28156a6892846d443292988156752e7 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 9 Sep 2024 13:19:21 -0400 Subject: [PATCH 6/9] Subcycle thermal step --- .../inputs/thermal_deformation_cube.json | 1 + src/CabanaPD_Solver.hpp | 12 +++++++++--- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/examples/thermomechanics/inputs/thermal_deformation_cube.json b/examples/thermomechanics/inputs/thermal_deformation_cube.json index b7e054ec..9d78a057 100644 --- a/examples/thermomechanics/inputs/thermal_deformation_cube.json +++ b/examples/thermomechanics/inputs/thermal_deformation_cube.json @@ -10,6 +10,7 @@ "horizon" : {"value": 0.0075, "unit": "m"}, "final_time" : {"value": 8, "unit": "s"}, "timestep" : {"value": 5.1e-8, "unit": "s"}, + "thermal_subcycle_steps": {"value": 100}, "output_frequency" : {"value": 1}, "output_reference" : {"value": true} } diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 884958a9..a12e173f 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -141,9 +141,11 @@ class SolverElastic // Create heat transfer if needed. if constexpr ( std::is_same::value ) + { + thermal_subcycle_steps = inputs["thermal_subcycle_steps"]; heat_transfer = std::make_shared( inputs["half_neigh"], force_model ); - + } force = std::make_shared( inputs["half_neigh"], force_model ); @@ -274,8 +276,11 @@ class SolverElastic if constexpr ( std::is_same::value ) - computeHeatTransfer( *heat_transfer, *particles, *neighbors, - neigh_iter_tag{}, dt ); + { + if ( step % thermal_subcycle_steps == 0 ) + computeHeatTransfer( *heat_transfer, *particles, *neighbors, + neigh_iter_tag{}, dt ); + } if constexpr ( std::is_same::value ) comm->gatherTemperature(); @@ -444,6 +449,7 @@ class SolverElastic int output_frequency; bool output_reference; double dt; + int thermal_subcycle_steps; protected: // Core modules. From 7e92e5baa7892940825d9749296237e08f64078b Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 24 Sep 2024 16:03:24 -0400 Subject: [PATCH 7/9] Update heat transfer examples --- .../thermomechanics/inputs/heat_transfer.json | 16 +++ .../thermomechanics/inputs/thermal_crack.json | 3 +- .../inputs/thermal_deformation.json | 1 + .../inputs/thermal_deformation_cube_temp.json | 15 --- ...=> thermal_deformation_heat_transfer.json} | 10 +- examples/thermomechanics/thermal_crack.cpp | 2 +- .../thermal_deformation_heat_transfer.cpp | 112 ++---------------- 7 files changed, 38 insertions(+), 121 deletions(-) create mode 100644 examples/thermomechanics/inputs/heat_transfer.json delete mode 100644 examples/thermomechanics/inputs/thermal_deformation_cube_temp.json rename examples/thermomechanics/inputs/{thermal_deformation_cube.json => thermal_deformation_heat_transfer.json} (67%) diff --git a/examples/thermomechanics/inputs/heat_transfer.json b/examples/thermomechanics/inputs/heat_transfer.json new file mode 100644 index 00000000..dcda7d0b --- /dev/null +++ b/examples/thermomechanics/inputs/heat_transfer.json @@ -0,0 +1,16 @@ +{ + "num_cells" : {"value": [49, 49, 49]}, + "system_size" : {"value": [0.1, 0.1, 0.1], "unit": "m"}, + "density" : {"value": 8915, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 115e+9, "unit": "Pa"}, + "thermal_expansion_coeff" : {"value": 17e-6, "unit": "oC^{-1}"}, + "thermal_conductivity" : {"value": 387, "unit": "W/(m.K)"}, + "specific_heat_capacity" : {"value": 385, "unit": "J/(kg.K)"}, + "reference_temperature" : {"value": 100.0, "unit": "oC"}, + "horizon" : {"value": 0.00615, "unit": "m"}, + "final_time" : {"value": 8, "unit": "s"}, + "timestep" : {"value": 0.02, "unit": "s", "note": "This timestep is too large for mechanics to be stable. Use thermal_deformation_heat_transfer.json instead if coupled thermomechanics is desired."}, + "thermal_subcycle_steps" : {"value": 1}, + "output_frequency" : {"value": 10}, + "output_reference" : {"value": true} +} diff --git a/examples/thermomechanics/inputs/thermal_crack.json b/examples/thermomechanics/inputs/thermal_crack.json index 4427f6c6..bb15a543 100644 --- a/examples/thermomechanics/inputs/thermal_crack.json +++ b/examples/thermomechanics/inputs/thermal_crack.json @@ -4,13 +4,14 @@ "density" : {"value": 3980, "unit": "kg/m^3"}, "elastic_modulus" : {"value": 370e+9, "unit": "Pa"}, "fracture_energy" : {"value": 24.32, "unit": "J/m^2"}, - "thermal_coefficient" : {"value": 7.5E-6, "unit": "oC^{-1}"}, + "thermal_expansion_coeff" : {"value": 7.5E-6, "unit": "oC^{-1}"}, "reference_temperature" : {"value": 300.0, "unit": "oC"}, "background_temperature" : {"value": 20.0, "unit": "oC"}, "surface_temperature_ramp_time" : {"value": 0.001, "unit": "s"}, "horizon" : {"value": 3.0e-4, "unit": "m"}, "final_time" : {"value": 7.5e-4, "unit": "s"}, "timestep" : {"value": 7.5E-9, "unit": "s"}, + "thermal_subcycle_steps" : {"value": 100}, "output_frequency" : {"value": 200}, "output_reference" : {"value": true} } diff --git a/examples/thermomechanics/inputs/thermal_deformation.json b/examples/thermomechanics/inputs/thermal_deformation.json index 83c8874e..7249c936 100644 --- a/examples/thermomechanics/inputs/thermal_deformation.json +++ b/examples/thermomechanics/inputs/thermal_deformation.json @@ -10,6 +10,7 @@ "horizon" : {"value": 0.03, "unit": "m"}, "final_time" : {"value": 0.01, "unit": "s"}, "timestep" : {"value": 7.5E-7, "unit": "s"}, + "thermal_subcycle_steps": {"value": 100}, "output_frequency" : {"value": 100}, "output_reference" : {"value": true} } diff --git a/examples/thermomechanics/inputs/thermal_deformation_cube_temp.json b/examples/thermomechanics/inputs/thermal_deformation_cube_temp.json deleted file mode 100644 index 163195b0..00000000 --- a/examples/thermomechanics/inputs/thermal_deformation_cube_temp.json +++ /dev/null @@ -1,15 +0,0 @@ -{ - "num_cells" : {"value": [51, 51, 51]}, - "system_size" : {"value": [0.5, 0.5, 0.5], "unit": "m"}, - "density" : {"value": 3980, "unit": "kg/m^3"}, - "elastic_modulus" : {"value": 370e+9, "unit": "Pa"}, - "thermal_expansion_coeff" : {"value": 7.5E-6, "unit": "oC^{-1}"}, - "thermal_conductivity" : {"value": 31, "unit": "W/(m.K)"}, - "specific_heat_capacity" : {"value": 880, "unit": "J/(kg.K)"}, - "reference_temperature" : {"value": 100.0, "unit": "oC"}, - "horizon" : {"value": 0.03, "unit": "m"}, - "final_time" : {"value": 0.01, "unit": "s"}, - "timestep" : {"value": 1.0E-7, "unit": "s"}, - "output_frequency" : {"value": 100}, - "output_reference" : {"value": true} -} diff --git a/examples/thermomechanics/inputs/thermal_deformation_cube.json b/examples/thermomechanics/inputs/thermal_deformation_heat_transfer.json similarity index 67% rename from examples/thermomechanics/inputs/thermal_deformation_cube.json rename to examples/thermomechanics/inputs/thermal_deformation_heat_transfer.json index 9d78a057..c3b68c5d 100644 --- a/examples/thermomechanics/inputs/thermal_deformation_cube.json +++ b/examples/thermomechanics/inputs/thermal_deformation_heat_transfer.json @@ -1,5 +1,5 @@ { - "num_cells" : {"value": [40, 40, 40]}, + "num_cells" : {"value": [49, 49, 49]}, "system_size" : {"value": [0.1, 0.1, 0.1], "unit": "m"}, "density" : {"value": 8915, "unit": "kg/m^3"}, "elastic_modulus" : {"value": 115e+9, "unit": "Pa"}, @@ -7,10 +7,10 @@ "thermal_conductivity" : {"value": 387, "unit": "W/(m.K)"}, "specific_heat_capacity" : {"value": 385, "unit": "J/(kg.K)"}, "reference_temperature" : {"value": 100.0, "unit": "oC"}, - "horizon" : {"value": 0.0075, "unit": "m"}, + "horizon" : {"value": 0.00615, "unit": "m"}, "final_time" : {"value": 8, "unit": "s"}, - "timestep" : {"value": 5.1e-8, "unit": "s"}, - "thermal_subcycle_steps": {"value": 100}, - "output_frequency" : {"value": 1}, + "timestep" : {"value": 4.2e-07, "unit": "s"}, + "thermal_subcycle_steps" : {"value": 5e+4}, + "output_frequency" : {"value": 10}, "output_reference" : {"value": true} } diff --git a/examples/thermomechanics/thermal_crack.cpp b/examples/thermomechanics/thermal_crack.cpp index 171eafb7..c29c16ef 100644 --- a/examples/thermomechanics/thermal_crack.cpp +++ b/examples/thermomechanics/thermal_crack.cpp @@ -43,7 +43,7 @@ void thermalCrackExample( const std::string filename ) double K = E / ( 3 * ( 1 - 2 * nu ) ); double G0 = inputs["fracture_energy"]; double delta = inputs["horizon"]; - double alpha = inputs["thermal_coefficient"]; + double alpha = inputs["thermal_expansion_coeff"]; // Problem parameters double temp0 = inputs["reference_temperature"]; diff --git a/examples/thermomechanics/thermal_deformation_heat_transfer.cpp b/examples/thermomechanics/thermal_deformation_heat_transfer.cpp index 1259beba..021ac974 100644 --- a/examples/thermomechanics/thermal_deformation_heat_transfer.cpp +++ b/examples/thermomechanics/thermal_deformation_heat_transfer.cpp @@ -18,7 +18,7 @@ #include -// Simulate thermally-induced deformation in a rectangular plate. +// Simulate heat transfer in a pseudo-1d cube. void thermalDeformationHeatTransferExample( const std::string filename ) { // ==================================================== @@ -104,89 +104,29 @@ void thermalDeformationHeatTransferExample( const std::string filename ) // ==================================================== // Boundary condition // ==================================================== - // EXAMPLE 1: Temperature profile imposed over entire domain - using plane_type = CabanaPD::RegionBoundary; - /* - plane_type plane( low_corner[0], high_corner[0], - low_corner[1], high_corner[1], - low_corner[2], high_corner[2] ); - - std::vector planes = { plane }; - */ - - // EXAMPLE 2: Temperature profile imposed on top, bottom, left, and right - // surfaces - double dx = particles->dx[0]; + // Temperature profile imposed on top and bottom surfaces double dy = particles->dx[1]; - double dz = particles->dx[2]; - - // Left surface: x-direction - plane_type plane1( low_corner[0] - dx, low_corner[0] + dx, low_corner[1], - high_corner[1], low_corner[2], high_corner[2] ); - - // Right surface: x-direction - plane_type plane2( high_corner[0] - dx, high_corner[0] + dx, low_corner[1], - high_corner[1], low_corner[2], high_corner[2] ); + using plane_type = CabanaPD::RegionBoundary; - // Top surface: y-direction - plane_type plane3( low_corner[0], high_corner[0], high_corner[1] - dy, + // Top surface + plane_type plane1( low_corner[0], high_corner[0], high_corner[1] - dy, high_corner[1] + dy, low_corner[2], high_corner[2] ); - // Bottom surface: y-direction - plane_type plane4( low_corner[0], high_corner[0], low_corner[1] - dy, + // Bottom surface + plane_type plane2( low_corner[0], high_corner[0], low_corner[1] - dy, low_corner[1] + dy, low_corner[2], high_corner[2] ); - // Front surface: z-direction - plane_type plane5( low_corner[0], high_corner[0], low_corner[1], - high_corner[1], low_corner[2] - dz, low_corner[2] + dz ); - - // Back surface: z-direction - plane_type plane6( low_corner[0], high_corner[0], low_corner[1], - high_corner[1], high_corner[2] - dz, - high_corner[2] + dz ); - - // std::vector planes = { plane1, plane2, plane3, plane4 }; - // std::vector planes = { plane1, plane2, plane3, plane4, - // plane5, plane6 }; std::vector planes = { plane1, plane2 }; - std::vector planes = { plane3, plane4 }; - /* - // EXAMPLE 3: Temperature profile imposed on top, bottom, left, and - right - // nonlocal boundaries (width delta) Top surface - plane_type plane1( - low_corner[0], high_corner[0], high_corner[1] - delta, - high_corner[1] + delta, low_corner[2], high_corner[2] ); - - // Bottom surface - plane_type plane2( - low_corner[0], high_corner[0], low_corner[1] - delta, - low_corner[1] + delta, low_corner[2], high_corner[2] ); - - // Left surface - plane_type plane3( - low_corner[0] - delta, low_corner[0] + delta, low_corner[1], - high_corner[1], low_corner[2], high_corner[2] ); - - // Right surface - plane_type plane4( - high_corner[0] - delta, high_corner[0] + delta, low_corner[1], - high_corner[1], low_corner[2], high_corner[2] ); - - std::vector planes = { plane1, plane2, plane3, - plane4 }; - */ - - auto x = particles->sliceReferencePosition(); - const double low_corner_y = low_corner[1]; // This is purposely delayed until after solver init so that ghosted // particles are correctly taken into account for lambda capture here. - auto temp_func = KOKKOS_LAMBDA( const int pid, const double t ) + temp = particles->sliceTemperature(); + auto temp_bc = KOKKOS_LAMBDA( const int pid, const double ) { - // temp( pid ) = temp0 + 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; temp( pid ) = 0.0; }; + + std::vector planes = { plane1, plane2 }; auto bc = CabanaPD::createBoundaryCondition( - temp_func, exec_space{}, *particles, planes, false, 1.0 ); + temp_bc, exec_space{}, *particles, planes, false, 1.0 ); // ==================================================== // Simulation run @@ -197,38 +137,12 @@ void thermalDeformationHeatTransferExample( const std::string filename ) // ==================================================== // Outputs // ==================================================== - // Output temperature along the y-axis - auto temp_output = particles->sliceTemperature(); - auto value = KOKKOS_LAMBDA( const int pid ) { return temp_output( pid ); }; - int profile_dim = 1; + auto value = KOKKOS_LAMBDA( const int pid ) { return temp( pid ); }; std::string file_name = "temperature_yaxis_profile.txt"; createOutputProfile( MPI_COMM_WORLD, num_cells[1], profile_dim, file_name, *particles, value ); - - /* - // Output y-displacement along the x-axis - createDisplacementProfile( MPI_COMM_WORLD, - "ydisplacement_xaxis_profile.txt", *particles, - num_cells[0], 0, 1 ); - - // Output y-displacement along the y-axis - createDisplacementProfile( MPI_COMM_WORLD, - "ydisplacement_yaxis_profile.txt", *particles, - num_cells[1], 1, 1 ); - - // Output displacement magnitude along the x-axis - createDisplacementMagnitudeProfile( - MPI_COMM_WORLD, "displacement_magnitude_xaxis_profile.txt", *particles, - num_cells[0], 0 ); - - // Output displacement magnitude along the y-axis - createDisplacementMagnitudeProfile( - MPI_COMM_WORLD, "displacement_magnitude_yaxis_profile.txt", *particles, - num_cells[1], 1 ); - - */ } // Initialize MPI+Kokkos. From 0e4cc00dde4bd48d0d0e39752d807e4954cc7ab6 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 18 Oct 2024 11:39:06 -0400 Subject: [PATCH 8/9] Add type checkers for temperature --- src/CabanaPD_Solver.hpp | 36 ++++++++++++++++++------------------ src/CabanaPD_Types.hpp | 22 ++++++++++++++++++++++ 2 files changed, 40 insertions(+), 18 deletions(-) diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index a12e173f..86e552fd 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -134,13 +134,13 @@ class SolverElastic comm = std::make_shared( *particles ); // Update temperature ghost size if needed. - if constexpr ( std::is_same::value ) + if constexpr ( is_temperature_dependent< + typename force_model_type::thermal_type>::value ) force_model.update( particles->sliceTemperature() ); // Create heat transfer if needed. - if constexpr ( std::is_same::value ) + if constexpr ( is_heat_transfer< + typename force_model_type::thermal_type>::value ) { thermal_subcycle_steps = inputs["thermal_subcycle_steps"]; heat_transfer = std::make_shared( @@ -243,8 +243,8 @@ class SolverElastic boundary_condition.apply( exec_space(), *particles, 0.0 ); // Communicate temperature. - if constexpr ( std::is_same::value ) + if constexpr ( is_temperature_dependent< + typename force_model_type::thermal_type>::value ) comm->gatherTemperature(); // Force init without particle output. @@ -274,15 +274,15 @@ class SolverElastic // Update ghost particles. comm->gatherDisplacement(); - if constexpr ( std::is_same::value ) + if constexpr ( is_heat_transfer< + typename force_model_type::thermal_type>::value ) { if ( step % thermal_subcycle_steps == 0 ) computeHeatTransfer( *heat_transfer, *particles, *neighbors, neigh_iter_tag{}, dt ); } - if constexpr ( std::is_same::value ) + if constexpr ( is_temperature_dependent< + typename force_model_type::thermal_type>::value ) comm->gatherTemperature(); // Add non-force boundary condition. @@ -324,8 +324,8 @@ class SolverElastic // Compute internal forces. updateForce(); - if constexpr ( std::is_same::value ) + if constexpr ( is_temperature_dependent< + typename force_model_type::thermal_type>::value ) comm->gatherTemperature(); // Integrate - velocity Verlet second half. @@ -548,8 +548,8 @@ class SolverFracture boundary_condition.apply( exec_space(), *particles, 0.0 ); // Communicate temperature. - if constexpr ( std::is_same::value ) + if constexpr ( is_temperature_dependent< + typename force_model_type::thermal_type>::value ) comm->gatherTemperature(); // Force init without particle output. @@ -580,8 +580,8 @@ class SolverFracture if ( !boundary_condition.forceUpdate() ) boundary_condition.apply( exec_space(), *particles, step * dt ); - if constexpr ( std::is_same::value ) + if constexpr ( is_temperature_dependent< + typename force_model_type::thermal_type>::value ) comm->gatherTemperature(); // Update ghost particles. @@ -616,8 +616,8 @@ class SolverFracture // Integrate - velocity Verlet first half. integrator->initialHalfStep( *particles ); - if constexpr ( std::is_same::value ) + if constexpr ( is_temperature_dependent< + typename force_model_type::thermal_type>::value ) comm->gatherTemperature(); // Update ghost particles. diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index b1fa4399..2e47a9c9 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -34,6 +34,28 @@ struct DynamicTemperature : public TemperatureDependent using base_type = TemperatureDependent; }; +//! Static type checkers. +template +struct is_temperature_dependent : public std::false_type +{ +}; +template <> +struct is_temperature_dependent : public std::true_type +{ +}; +template <> +struct is_temperature_dependent : public std::true_type +{ +}; +template +struct is_heat_transfer : public std::false_type +{ +}; +template <> +struct is_heat_transfer : public std::true_type +{ +}; + // Model types. struct PMB { From 0be8643f31682c50f2cd5898945840848eb6e960 Mon Sep 17 00:00:00 2001 From: pabloseleson Date: Wed, 6 Nov 2024 11:50:52 -0500 Subject: [PATCH 9/9] Fixed solver by changing order of temp b.c. and temp gather --- src/CabanaPD_Solver.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 86e552fd..44c7ab91 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -281,14 +281,15 @@ class SolverElastic computeHeatTransfer( *heat_transfer, *particles, *neighbors, neigh_iter_tag{}, dt ); } - if constexpr ( is_temperature_dependent< - typename force_model_type::thermal_type>::value ) - comm->gatherTemperature(); // Add non-force boundary condition. if ( !boundary_condition.forceUpdate() ) boundary_condition.apply( exec_space(), *particles, step * dt ); + if constexpr ( is_temperature_dependent< + typename force_model_type::thermal_type>::value ) + comm->gatherTemperature(); + // Compute internal forces. updateForce();