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 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/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 4f0f6908..7249c936 100644 --- a/examples/thermomechanics/inputs/thermal_deformation.json +++ b/examples/thermomechanics/inputs/thermal_deformation.json @@ -1,13 +1,16 @@ { - "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": [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"}, + "thermal_subcycle_steps": {"value": 100}, + "output_frequency" : {"value": 100}, + "output_reference" : {"value": true} } diff --git a/examples/thermomechanics/inputs/thermal_deformation_heat_transfer.json b/examples/thermomechanics/inputs/thermal_deformation_heat_transfer.json new file mode 100644 index 00000000..c3b68c5d --- /dev/null +++ b/examples/thermomechanics/inputs/thermal_deformation_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": 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.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..021ac974 --- /dev/null +++ b/examples/thermomechanics/thermal_deformation_heat_transfer.cpp @@ -0,0 +1,158 @@ +/**************************************************************************** + * 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 heat transfer in a pseudo-1d cube. +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 + // ==================================================== + // Temperature profile imposed on top and bottom surfaces + double dy = particles->dx[1]; + using plane_type = CabanaPD::RegionBoundary; + + // 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 + plane_type plane2( low_corner[0], high_corner[0], low_corner[1] - dy, + low_corner[1] + dy, low_corner[2], high_corner[2] ); + + // This is purposely delayed until after solver init so that ghosted + // particles are correctly taken into account for lambda capture here. + temp = particles->sliceTemperature(); + auto temp_bc = KOKKOS_LAMBDA( const int pid, const double ) + { + temp( pid ) = 0.0; + }; + + std::vector planes = { plane1, plane2 }; + auto bc = CabanaPD::createBoundaryCondition( + temp_bc, exec_space{}, *particles, planes, false, 1.0 ); + + // ==================================================== + // Simulation run + // ==================================================== + cabana_pd->init( bc ); + cabana_pd->run( bc ); + + // ==================================================== + // Outputs + // ==================================================== + // Output temperature along the y-axis + 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 ); +} + +// Initialize MPI+Kokkos. +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + thermalDeformationHeatTransferExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} 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..44c7ab91 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -74,6 +74,7 @@ #include #include #include +#include #include #include #include @@ -98,18 +99,22 @@ 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; 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 module types. + using heat_transfer_type = HeatTransfer; + SolverElastic( input_type _inputs, std::shared_ptr _particles, force_model_type force_model ) @@ -129,10 +134,18 @@ 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 ( is_heat_transfer< + typename force_model_type::thermal_type>::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 ); @@ -230,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. @@ -258,17 +271,25 @@ class SolverElastic // Integrate - velocity Verlet first half. integrator->initialHalfStep( *particles ); + // Update ghost particles. + comm->gatherDisplacement(); + + 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 ); + } + // Add non-force boundary condition. 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. - comm->gatherDisplacement(); - // Compute internal forces. updateForce(); @@ -304,8 +325,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. @@ -429,19 +450,24 @@ class SolverElastic int output_frequency; bool output_reference; double dt; + int thermal_subcycle_steps; 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 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; @@ -523,8 +549,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. @@ -555,8 +581,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. @@ -591,8 +617,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 62076735..2e47a9c9 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -14,19 +14,49 @@ namespace CabanaPD { +// Mechanics types. struct Elastic { }; struct Fracture { }; + +// Thermal types. +struct TemperatureIndependent +{ +}; struct TemperatureDependent { }; -struct TemperatureIndependent +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 { }; 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