From b64bd3b1a0255973f9776ab6524b4905b7f1eb1b Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 7 May 2024 16:13:43 -0400 Subject: [PATCH 01/11] Add temperature field --- src/CabanaPD_Boundary.hpp | 8 ++++---- src/CabanaPD_Particles.hpp | 35 +++++++++++++++++++++++------------ src/CabanaPD_Solver.hpp | 8 ++++---- 3 files changed, 31 insertions(+), 20 deletions(-) diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 23be020d..31e97959 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -167,7 +167,7 @@ struct BoundaryCondition } template - void apply( ExecSpace, ParticleType& ) + void apply( ExecSpace, ParticleType&, double ) { _timer.start(); auto user = _user_functor; @@ -196,7 +196,7 @@ struct BoundaryCondition } template - void apply( ExecSpace, ParticleType ) + void apply( ExecSpace, ParticleType, double ) { } @@ -226,7 +226,7 @@ struct BoundaryCondition } template - void apply( ExecSpace, ParticleType& particles ) + void apply( ExecSpace, ParticleType& particles, double ) { _timer.start(); auto f = particles.sliceForce(); @@ -268,7 +268,7 @@ struct BoundaryCondition } template - void apply( ExecSpace, ParticleType& particles ) + void apply( ExecSpace, ParticleType& particles, double ) { _timer.start(); diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 02d691c2..0dba3d33 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -101,9 +101,9 @@ class Particles using scalar_type = Cabana::MemberTypes; // no-fail. using int_type = Cabana::MemberTypes; - // type, W, v, rho, damage. + // v, W, rho, damage, temperature, type. using other_types = - Cabana::MemberTypes; + Cabana::MemberTypes; // Potentially needed later: body force (b), ID. // FIXME: add vector length. @@ -258,6 +258,7 @@ class Particles auto y = sliceCurrentPosition(); auto vol = sliceVolume(); auto nofail = sliceNoFail(); + auto temp = sliceTemperature(); // Initialize particles. auto create_functor = @@ -286,6 +287,7 @@ class Particles type( pid ) = 0; nofail( pid ) = 0; rho( pid ) = 1.0; + temp( pid ) = 0.0; return create; }; @@ -354,8 +356,8 @@ class Particles { return Cabana::slice<0>( _aosoa_vol, "volume" ); } - auto sliceType() { return Cabana::slice<0>( _aosoa_other, "type" ); } - auto sliceType() const { return Cabana::slice<0>( _aosoa_other, "type" ); } + auto sliceType() { return Cabana::slice<5>( _aosoa_other, "type" ); } + auto sliceType() const { return Cabana::slice<5>( _aosoa_other, "type" ); } auto sliceStrainEnergy() { return Cabana::slice<1>( _aosoa_other, "strain_energy" ); @@ -366,21 +368,29 @@ class Particles } auto sliceVelocity() { - return Cabana::slice<2>( _aosoa_other, "velocities" ); + return Cabana::slice<0>( _aosoa_other, "velocities" ); } auto sliceVelocity() const { - return Cabana::slice<2>( _aosoa_other, "velocities" ); + return Cabana::slice<0>( _aosoa_other, "velocities" ); } - auto sliceDensity() { return Cabana::slice<3>( _aosoa_other, "density" ); } + auto sliceDensity() { return Cabana::slice<2>( _aosoa_other, "density" ); } auto sliceDensity() const { - return Cabana::slice<3>( _aosoa_other, "density" ); + return Cabana::slice<2>( _aosoa_other, "density" ); } - auto sliceDamage() { return Cabana::slice<4>( _aosoa_other, "damage" ); } + auto sliceDamage() { return Cabana::slice<3>( _aosoa_other, "damage" ); } auto sliceDamage() const { - return Cabana::slice<4>( _aosoa_other, "damage" ); + return Cabana::slice<3>( _aosoa_other, "damage" ); + } + auto sliceTemperature() + { + return Cabana::slice<4>( _aosoa_other, "temperature" ); + } + auto sliceTemperature() const + { + return Cabana::slice<4>( _aosoa_other, "temperature" ); } auto sliceNoFail() { @@ -447,7 +457,8 @@ class Particles Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( h5_config, "particles", MPI_COMM_WORLD, output_step, output_time, n_local, getPosition( use_reference ), sliceStrainEnergy(), - sliceForce(), sliceDisplacement(), sliceVelocity(), sliceDamage() ); + sliceForce(), sliceDisplacement(), sliceVelocity(), sliceDamage(), + sliceTemperature() ); #else #ifdef Cabana_ENABLE_SILO Cabana::Grid::Experimental::SiloParticleOutput:: @@ -455,7 +466,7 @@ class Particles "particles", local_grid->globalGrid(), output_step, output_time, 0, n_local, getPosition( use_reference ), sliceStrainEnergy(), sliceForce(), sliceDisplacement(), sliceVelocity(), - sliceDamage() ); + sliceDamage(), sliceTemperature() ); #else log( std::cout, "No particle output enabled." ); #endif diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 6f46bfb5..0515f367 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -201,7 +201,7 @@ class SolverElastic computeEnergy( *force, *particles, *neighbors, neigh_iter_tag() ); // Add boundary condition. - boundary_condition.apply( exec_space(), *particles ); + boundary_condition.apply( exec_space(), *particles, 0.0 ); particles->output( 0, 0.0, output_reference ); } @@ -231,7 +231,7 @@ class SolverElastic computeForce( *force, *particles, *neighbors, neigh_iter_tag{} ); // Add boundary condition. - boundary_condition.apply( exec_space(), *particles ); + boundary_condition.apply( exec_space(), *particles, step * dt ); // Integrate - velocity Verlet second half. integrator->finalHalfStep( *particles ); @@ -430,7 +430,7 @@ class SolverFracture computeEnergy( *force, *particles, *neighbors, mu, neigh_iter_tag() ); // Add boundary condition. - boundary_condition.apply( exec_space(), *particles ); + boundary_condition.apply( exec_space(), *particles, 0 ); particles->output( 0, 0.0, output_reference ); } @@ -462,7 +462,7 @@ class SolverFracture neigh_iter_tag{} ); // Add boundary condition. - boundary_condition.apply( exec_space{}, *particles ); + boundary_condition.apply( exec_space{}, *particles, step * dt ); // Integrate - velocity Verlet second half. integrator->finalHalfStep( *particles ); From 44666213a7007b992a55119e02ad29412c7afda8 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 6 May 2024 09:21:17 -0400 Subject: [PATCH 02/11] Add temperature-based models --- src/CabanaPD_ForceModels.hpp | 45 +++++++++++++- src/CabanaPD_Types.hpp | 6 ++ src/force/CabanaPD_ForceModels_PMB.hpp | 62 +++++++++++++++++-- src/force/CabanaPD_Force_PMB.hpp | 86 ++++++++++++++++---------- 4 files changed, 161 insertions(+), 38 deletions(-) diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index af1cf283..88bf7be5 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -23,9 +23,52 @@ struct BaseForceModel : delta( _delta ){}; BaseForceModel( const BaseForceModel& model ) { delta = model.delta; } + + // No-op for temperature. + KOKKOS_INLINE_FUNCTION + void thermalStretch( double&, const int, const int ) const {} +}; + +template +struct BaseTemperatureModel +{ + double alpha; + double temp0; + + // Temperature field + TemperatureType temperature; + + BaseTemperatureModel(){}; + BaseTemperatureModel( const TemperatureType _temp, const double _alpha, + const double _temp0 ) + : alpha( _alpha ) + , temp0( _temp0 ) + , temperature( _temp ){}; + + BaseTemperatureModel( const BaseTemperatureModel& model ) + { + alpha = model.alpha; + temp0 = model.temp0; + temperature = model.temperature; + } + + void set_param( const double _alpha, const double _temp0 ) + { + alpha = _alpha; + temp0 = _temp0; + } + + // Update stretch with temperature effects. + KOKKOS_INLINE_FUNCTION + void thermalStretch( double& s, const int i, const int j ) const + { + double temp_avg = 0.5 * ( temperature( i ) + temperature( j ) ) - temp0; + s -= alpha * temp_avg; + } }; -template +template struct ForceModel; } // namespace CabanaPD diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index a036d867..fd95375a 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -20,6 +20,12 @@ struct Elastic struct Fracture { }; +struct TemperatureDependent +{ +}; +struct TemperatureIndependent +{ +}; struct PMB { diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp index 3d36740a..2b71f599 100644 --- a/src/force/CabanaPD_ForceModels_PMB.hpp +++ b/src/force/CabanaPD_ForceModels_PMB.hpp @@ -18,11 +18,12 @@ namespace CabanaPD { template <> -struct ForceModel : public BaseForceModel +struct ForceModel : public BaseForceModel { using base_type = BaseForceModel; using base_model = PMB; using fracture_type = Elastic; + using thermal_type = TemperatureIndependent; using base_type::delta; @@ -52,7 +53,8 @@ struct ForceModel : public BaseForceModel }; template <> -struct ForceModel : public ForceModel +struct ForceModel + : public ForceModel { using base_type = ForceModel; using base_model = typename base_type::base_model; @@ -90,7 +92,8 @@ struct ForceModel : public ForceModel }; template <> -struct ForceModel : public ForceModel +struct ForceModel + : public ForceModel { using base_type = ForceModel; using base_model = typename base_type::base_model; @@ -104,7 +107,8 @@ struct ForceModel : public ForceModel }; template <> -struct ForceModel : public ForceModel +struct ForceModel + : public ForceModel { using base_type = ForceModel; using base_model = typename base_type::base_model; @@ -121,6 +125,56 @@ struct ForceModel : public ForceModel using base_type::s0; }; +template +struct ForceModel + : public ForceModel, + BaseTemperatureModel +{ + using base_type = ForceModel; + using base_temperature_type = BaseTemperatureModel; + using base_model = PMB; + using fracture_type = Elastic; + using thermal_type = TemperatureDependent; + + using base_type::c; + using base_type::delta; + using base_type::K; + + // Thermal parameters + using base_temperature_type::alpha; + using base_temperature_type::temp0; + + // Explicitly use the temperature-dependent stretch. + using base_temperature_type::thermalStretch; + + // ForceModel(){}; + ForceModel( const double _delta, const double _K, + const TemperatureType _temp, const double _alpha, + const double _temp0 = 0.0 ) + : base_type( _delta, _K ) + , base_temperature_type( _temp, _alpha, _temp0 ) + { + set_param( _delta, _K, _alpha, _temp0 ); + } + + void set_param( const double _delta, const double _K, const double _alpha, + const double _temp0 ) + { + base_type::set_param( _delta, _K ); + base_temperature_type::set_param( _alpha, _temp0 ); + } +}; + +template +auto createForceModel( ParticleType particles, const double delta, + const double K, const double alpha, const double temp0 ) +{ + auto temp = particles.sliceTemperature(); + using temp_type = decltype( temp ); + return ForceModel( + delta, K, temp, alpha, temp0 ); +} } // namespace CabanaPD #endif diff --git a/src/force/CabanaPD_Force_PMB.hpp b/src/force/CabanaPD_Force_PMB.hpp index 36f1223b..54991207 100644 --- a/src/force/CabanaPD_Force_PMB.hpp +++ b/src/force/CabanaPD_Force_PMB.hpp @@ -69,20 +69,22 @@ namespace CabanaPD { -template -class Force> +template +class Force> { + public: + using exec_space = ExecutionSpace; + using model_type = ForceModel; + protected: bool _half_neigh; - ForceModel _model; + model_type _model; Timer _timer; Timer _energy_timer; public: - using exec_space = ExecutionSpace; - - Force( const bool half_neigh, const ForceModel model ) + Force( const bool half_neigh, const model_type model ) : _half_neigh( half_neigh ) , _model( model ) { @@ -108,7 +110,7 @@ class Force> { _timer.start(); - auto c = _model.c; + auto model = _model; const auto vol = particles.sliceVolume(); auto force_full = KOKKOS_LAMBDA( const int i, const int j ) @@ -120,7 +122,10 @@ class Force> double xi, r, s; double rx, ry, rz; getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); - const double coeff = c * s * vol( j ); + + model.thermalStretch( s, i, j ); + + const double coeff = model.c * s * vol( j ); fx_i = coeff * rx / r; fy_i = coeff * ry / r; fz_i = coeff * rz / r; @@ -147,7 +152,7 @@ class Force> { _energy_timer.start(); - auto c = _model.c; + auto model = _model; const auto vol = particles.sliceVolume(); auto energy_full = @@ -157,9 +162,11 @@ class Force> double xi, r, s; getDistance( x, u, i, j, xi, r, s ); + model.thermalStretch( s, i, j ); + // 0.25 factor is due to 1/2 from outside the integral and 1/2 from // the integrand (pairwise potential). - double w = 0.25 * c * s * s * xi * vol( j ); + double w = 0.25 * model.c * s * s * xi * vol( j ); W( i ) += w; Phi += w * vol( i ); }; @@ -179,22 +186,25 @@ class Force> auto timeEnergy() { return _energy_timer.time(); }; }; -template -class Force> - : public Force> +template +class Force> + : public Force> { + public: + using exec_space = ExecutionSpace; + using model_type = ForceModel; + protected: - using base_type = Force>; + using base_type = + Force>; using base_type::_half_neigh; - ForceModel _model; + model_type _model; using base_type::_energy_timer; using base_type::_timer; public: - using exec_space = ExecutionSpace; - - Force( const bool half_neigh, const ForceModel model ) + Force( const bool half_neigh, const model_type model ) : base_type( half_neigh, model ) , _model( model ) { @@ -209,7 +219,7 @@ class Force> { _timer.start(); - auto c = _model.c; + auto model = _model; auto break_coeff = _model.bond_break_coeff; const auto vol = particles.sliceVolume(); const auto nofail = particles.sliceNoFail(); @@ -234,6 +244,8 @@ class Force> double rx, ry, rz; getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); + model.thermalStretch( s, i, j ); + // Break if beyond critical stretch unless in no-fail zone. if ( r * r >= break_coeff * xi * xi && !nofail( i ) && !nofail( j ) ) @@ -243,7 +255,7 @@ class Force> // Else if statement is only for performance. else if ( mu( i, n ) > 0 ) { - const double coeff = c * s * vol( j ); + const double coeff = model.c * s * vol( j ); double muij = mu( i, n ); fx_i = muij * coeff * rx / r; fy_i = muij * coeff * ry / r; @@ -272,7 +284,7 @@ class Force> { _energy_timer.start(); - auto c = _model.c; + auto model = _model; const auto vol = particles.sliceVolume(); auto energy_full = KOKKOS_LAMBDA( const int i, double& Phi ) @@ -291,9 +303,11 @@ class Force> double xi, r, s; getDistance( x, u, i, j, xi, r, s ); + model.thermalStretch( s, i, j ); + // 0.25 factor is due to 1/2 from outside the integral and 1/2 // from the integrand (pairwise potential). - double w = mu( i, n ) * 0.25 * c * s * s * xi * vol( j ); + double w = mu( i, n ) * 0.25 * model.c * s * s * xi * vol( j ); W( i ) += w; phi_i += mu( i, n ) * vol( j ); @@ -313,22 +327,24 @@ class Force> } }; -template -class Force> - : public Force> +template +class Force> + : public Force> { + public: + using exec_space = ExecutionSpace; + using model_type = ForceModel; + protected: using base_type = Force>; using base_type::_half_neigh; - ForceModel _model; + model_type _model; using base_type::_energy_timer; using base_type::_timer; public: - using exec_space = ExecutionSpace; - - Force( const bool half_neigh, const ForceModel model ) + Force( const bool half_neigh, const model_type model ) : base_type( half_neigh, model ) , _model( model ) { @@ -343,7 +359,7 @@ class Force> { _timer.start(); - auto c = _model.c; + auto model = _model; const auto vol = particles.sliceVolume(); auto force_full = KOKKOS_LAMBDA( const int i, const int j ) @@ -358,7 +374,9 @@ class Force> getLinearizedDistanceComponents( x, u, i, j, xi, linear_s, xi_x, xi_y, xi_z ); - const double coeff = c * linear_s * vol( j ); + model.thermalStretch( linear_s, i, j ); + + const double coeff = model.c * linear_s * vol( j ); fx_i = coeff * xi_x / xi; fy_i = coeff * xi_y / xi; fz_i = coeff * xi_z / xi; @@ -385,7 +403,7 @@ class Force> { _energy_timer.start(); - auto c = _model.c; + auto model = _model; const auto vol = particles.sliceVolume(); auto energy_full = @@ -395,9 +413,11 @@ class Force> double xi, linear_s; getLinearizedDistance( x, u, i, j, xi, linear_s ); + model.thermalStretch( linear_s, i, j ); + // 0.25 factor is due to 1/2 from outside the integral and 1/2 from // the integrand (pairwise potential). - double w = 0.25 * c * linear_s * linear_s * xi * vol( j ); + double w = 0.25 * model.c * linear_s * linear_s * xi * vol( j ); W( i ) += w; Phi += w * vol( i ); }; From 716eb6096bc7927fe26be823378ebec8922c9859 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 6 May 2024 09:22:05 -0400 Subject: [PATCH 03/11] Add thermal elastic example --- examples/CMakeLists.txt | 2 + examples/thermomechanics/CMakeLists.txt | 3 + .../inputs/thermal_deformation.json | 13 +++ .../thermomechanics/thermal_deformation.cpp | 82 +++++++++++++++++++ 4 files changed, 100 insertions(+) create mode 100644 examples/thermomechanics/CMakeLists.txt create mode 100644 examples/thermomechanics/inputs/thermal_deformation.json create mode 100644 examples/thermomechanics/thermal_deformation.cpp diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index df1b8e15..cc861300 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -8,3 +8,5 @@ add_executable(CrackBranching crack_branching.cpp) target_link_libraries(CrackBranching LINK_PUBLIC CabanaPD) install(TARGETS ElasticWave KalthoffWinkler CrackBranching DESTINATION ${CMAKE_INSTALL_BINDIR}) + +add_subdirectory(thermomechanics) diff --git a/examples/thermomechanics/CMakeLists.txt b/examples/thermomechanics/CMakeLists.txt new file mode 100644 index 00000000..84788a90 --- /dev/null +++ b/examples/thermomechanics/CMakeLists.txt @@ -0,0 +1,3 @@ +add_executable(ThermalDeformation thermal_deformation.cpp) +target_link_libraries(ThermalDeformation LINK_PUBLIC CabanaPD) +install(TARGETS ThermalDeformation DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/examples/thermomechanics/inputs/thermal_deformation.json b/examples/thermomechanics/inputs/thermal_deformation.json new file mode 100644 index 00000000..e6bc5205 --- /dev/null +++ b/examples/thermomechanics/inputs/thermal_deformation.json @@ -0,0 +1,13 @@ +{ + "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_coefficient" : {"value": 7.5E-6, "unit": "oC^{-1}"}, + "reference_temperature" : {"value": 0.0, "unit": "oC"}, + "horizon" : {"value": 0.03, "unit": "m"}, + "final_time" : {"value": 0.0093, "unit": "s"}, + "timestep" : {"value": 7.5E-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 new file mode 100644 index 00000000..576a2dc1 --- /dev/null +++ b/examples/thermomechanics/thermal_deformation.cpp @@ -0,0 +1,82 @@ +/**************************************************************************** + * 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 * + ****************************************************************************/ + +#include +#include + +#include "mpi.h" + +#include + +#include + +void thermalDeformationExample( const std::string filename ) +{ + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = typename exec_space::memory_space; + + CabanaPD::Inputs inputs( filename ); + double E = inputs["elastic_modulus"]; + double rho0 = inputs["density"]; + double nu = 0.25; // unitless + double K = E / ( 3 * ( 1 - 2 * nu ) ); + double delta = inputs["horizon"]; + + double alpha = inputs["thermal_coefficient"]; + double temp0 = inputs["reference_temperature"]; + + 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. + + // Choose force model type. + using model_type = CabanaPD::PMB; + + // Create particles from mesh. + auto particles = + std::make_shared>( + exec_space(), low_corner, high_corner, num_cells, halo_width ); + + auto temp = particles->sliceTemperature(); + auto x = particles->sliceReferencePosition(); + auto temp_func = KOKKOS_LAMBDA( const int pid, const double t ) + { + temp( pid ) = 5000.0 * ( x( pid, 1 ) - ( -0.15 ) ) * t; + }; + auto body_term = CabanaPD::createBodyTerm( temp_func ); + + auto rho = particles->sliceDensity(); + auto init_functor = KOKKOS_LAMBDA( const int pid ) { rho( pid ) = rho0; }; + particles->updateParticles( exec_space{}, init_functor ); + + auto force_model = + CabanaPD::createForceModel( + *particles, delta, K, alpha, temp0 ); + auto cabana_pd = CabanaPD::createSolverElastic( + inputs, particles, force_model, body_term ); + cabana_pd->init_force(); + cabana_pd->run(); +} + +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + thermalDeformationExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} From f196b2f53bf757b349485e91ec92110ba6c3815d Mon Sep 17 00:00:00 2001 From: pabloseleson Date: Wed, 8 May 2024 16:55:16 -0400 Subject: [PATCH 04/11] Formating changes to thermalDeformationExample --- .../thermomechanics/thermal_deformation.cpp | 50 ++++++++++++++++--- 1 file changed, 43 insertions(+), 7 deletions(-) diff --git a/examples/thermomechanics/thermal_deformation.cpp b/examples/thermomechanics/thermal_deformation.cpp index 576a2dc1..7d4fab1e 100644 --- a/examples/thermomechanics/thermal_deformation.cpp +++ b/examples/thermomechanics/thermal_deformation.cpp @@ -18,21 +18,38 @@ #include +// Simulate thermally-induced deformation in a rectangular plate. void thermalDeformationExample( 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 ); - double E = inputs["elastic_modulus"]; + + // ==================================================== + // Material and problem parameters + // ==================================================== + // Material parameters double rho0 = inputs["density"]; - double nu = 0.25; // unitless + double E = inputs["elastic_modulus"]; + double nu = 0.25; double K = E / ( 3 * ( 1 - 2 * nu ) ); double delta = inputs["horizon"]; - double alpha = inputs["thermal_coefficient"]; + + // 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"]; @@ -40,36 +57,55 @@ void thermalDeformationExample( const std::string filename ) ( ( high_corner[0] - low_corner[0] ) / num_cells[0] ) ); int halo_width = m + 1; // Just to be safe. - // Choose force model type. + // ==================================================== + // Force model type + // ==================================================== using model_type = CabanaPD::PMB; - // Create particles from mesh. + // ==================================================== + // Particle generation + // ==================================================== + // Does not set displacements, velocities, etc. auto particles = std::make_shared>( exec_space(), low_corner, high_corner, num_cells, halo_width ); - auto temp = particles->sliceTemperature(); + // ==================================================== + // Imposed field + // ==================================================== auto x = particles->sliceReferencePosition(); + auto temp = particles->sliceTemperature(); auto temp_func = KOKKOS_LAMBDA( const int pid, const double t ) { - temp( pid ) = 5000.0 * ( x( pid, 1 ) - ( -0.15 ) ) * t; + temp( pid ) = 5000.0 * ( x( pid, 1 ) - low_corner[1] ) * t; }; auto body_term = CabanaPD::createBodyTerm( temp_func ); + // ==================================================== + // Custom particle initialization + // ==================================================== auto rho = particles->sliceDensity(); auto init_functor = KOKKOS_LAMBDA( const int pid ) { rho( pid ) = rho0; }; particles->updateParticles( exec_space{}, init_functor ); + // ==================================================== + // Force model + // ==================================================== auto force_model = CabanaPD::createForceModel( *particles, delta, K, alpha, temp0 ); + + // ==================================================== + // Simulation run + // ==================================================== auto cabana_pd = CabanaPD::createSolverElastic( inputs, particles, force_model, body_term ); cabana_pd->init_force(); cabana_pd->run(); } +// Initialize MPI+Kokkos. int main( int argc, char* argv[] ) { MPI_Init( &argc, &argv ); From 7aa7344d1f3cadb157e4bdc8a9f9a7828f462442 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 17 May 2024 12:07:30 -0400 Subject: [PATCH 05/11] Propagate thermal type through Particles, Solver, Comm --- examples/crack_branching.cpp | 3 +- examples/elastic_wave.cpp | 3 +- examples/kalthoff_winkler.cpp | 3 +- .../thermomechanics/thermal_deformation.cpp | 13 +- src/CabanaPD_Comm.hpp | 35 +++- src/CabanaPD_ForceModels.hpp | 4 + src/CabanaPD_Particles.hpp | 175 +++++++++++++++--- src/CabanaPD_Solver.hpp | 23 ++- src/force/CabanaPD_ForceModels_LPS.hpp | 4 + src/force/CabanaPD_ForceModels_PMB.hpp | 3 + unit_test/tstComm.hpp | 8 +- unit_test/tstForce.hpp | 12 +- unit_test/tstIntegrator.hpp | 5 +- 13 files changed, 238 insertions(+), 53 deletions(-) diff --git a/examples/crack_branching.cpp b/examples/crack_branching.cpp index 40c2bad6..fd29321c 100644 --- a/examples/crack_branching.cpp +++ b/examples/crack_branching.cpp @@ -79,7 +79,8 @@ void crackBranchingExample( const std::string filename ) // ==================================================== // Does not set displacements, velocities, etc. auto particles = std::make_shared< - CabanaPD::Particles>( + CabanaPD::Particles>( exec_space(), low_corner, high_corner, num_cells, halo_width ); // ==================================================== diff --git a/examples/elastic_wave.cpp b/examples/elastic_wave.cpp index d8644abb..d70904a0 100644 --- a/examples/elastic_wave.cpp +++ b/examples/elastic_wave.cpp @@ -64,7 +64,8 @@ void elasticWaveExample( const std::string filename ) // ==================================================== // Does not set displacements, velocities, etc. auto particles = std::make_shared< - CabanaPD::Particles>( + CabanaPD::Particles>( exec_space(), low_corner, high_corner, num_cells, halo_width ); // ==================================================== diff --git a/examples/kalthoff_winkler.cpp b/examples/kalthoff_winkler.cpp index 5a7f189c..a3b2b904 100644 --- a/examples/kalthoff_winkler.cpp +++ b/examples/kalthoff_winkler.cpp @@ -89,7 +89,8 @@ void kalthoffWinklerExample( const std::string filename ) // ==================================================== // Does not set displacements, velocities, etc. auto particles = std::make_shared< - CabanaPD::Particles>( + CabanaPD::Particles>( exec_space(), low_corner, high_corner, num_cells, halo_width ); // ==================================================== diff --git a/examples/thermomechanics/thermal_deformation.cpp b/examples/thermomechanics/thermal_deformation.cpp index 7d4fab1e..e25687d6 100644 --- a/examples/thermomechanics/thermal_deformation.cpp +++ b/examples/thermomechanics/thermal_deformation.cpp @@ -61,23 +61,25 @@ void thermalDeformationExample( const std::string filename ) // Force model type // ==================================================== using model_type = CabanaPD::PMB; + using thermal_type = CabanaPD::TemperatureDependent; // ==================================================== // Particle generation // ==================================================== // Does not set displacements, velocities, etc. - auto particles = - std::make_shared>( - exec_space(), low_corner, high_corner, num_cells, halo_width ); + auto particles = std::make_shared< + CabanaPD::Particles>( + exec_space(), low_corner, high_corner, num_cells, halo_width ); // ==================================================== // Imposed field // ==================================================== auto x = particles->sliceReferencePosition(); auto temp = particles->sliceTemperature(); + const double low_corner_y = low_corner[1]; auto temp_func = KOKKOS_LAMBDA( const int pid, const double t ) { - temp( pid ) = 5000.0 * ( x( pid, 1 ) - low_corner[1] ) * t; + temp( pid ) = 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; }; auto body_term = CabanaPD::createBodyTerm( temp_func ); @@ -92,8 +94,7 @@ void thermalDeformationExample( const std::string filename ) // Force model // ==================================================== auto force_model = - CabanaPD::createForceModel( + CabanaPD::createForceModel( *particles, delta, K, alpha, temp0 ); // ==================================================== diff --git a/src/CabanaPD_Comm.hpp b/src/CabanaPD_Comm.hpp index a5f8b698..3abca2e4 100644 --- a/src/CabanaPD_Comm.hpp +++ b/src/CabanaPD_Comm.hpp @@ -230,12 +230,12 @@ struct HaloIds } }; -template +template class Comm; // FIXME: extract model from ParticleType instead. template -class Comm +class Comm { public: int mpi_size = -1; @@ -320,10 +320,11 @@ class Comm }; template -class Comm : public Comm +class Comm + : public Comm { public: - using base_type = Comm; + using base_type = Comm; using memory_space = typename base_type::memory_space; using halo_type = typename base_type::halo_type; using base_type::gather_u; @@ -348,6 +349,7 @@ class Comm : public Comm gather_theta = std::make_shared( *halo, particles._aosoa_theta ); + particles.resize( halo->numLocal(), halo->numGhost() ); _init_timer.stop(); } ~Comm() {} @@ -366,6 +368,31 @@ class Comm : public Comm } }; +template +class Comm + : public Comm +{ + public: + using base_type = Comm; + using memory_space = typename base_type::memory_space; + using halo_type = typename base_type::halo_type; + using base_type::halo; + + using gather_temp_type = + Cabana::Gather; + std::shared_ptr gather_temp; + + Comm( ParticleType& particles, int max_export_guess = 100 ) + : base_type( particles, max_export_guess ) + { + gather_temp = + std::make_shared( *halo, particles._aosoa_temp ); + particles.resize( halo->numLocal(), halo->numGhost() ); + } + + void gatherTemperature() { gather_temp->apply(); } +}; + } // namespace CabanaPD #endif diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index 88bf7be5..4a4b03d4 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -12,6 +12,8 @@ #ifndef FORCE_MODELS_H #define FORCE_MODELS_H +#include + namespace CabanaPD { struct BaseForceModel @@ -58,6 +60,8 @@ struct BaseTemperatureModel temp0 = _temp0; } + void update( const TemperatureType _temp ) { temperature = _temp; } + // Update stretch with temperature effects. KOKKOS_INLINE_FUNCTION void thermalStretch( double& s, const int i, const int j ) const diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 0dba3d33..839652a0 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -77,14 +77,16 @@ namespace CabanaPD { -template +template class Particles; template -class Particles +class Particles { public: - using self_type = Particles; + using self_type = + Particles; using memory_space = MemorySpace; using execution_space = typename memory_space::execution_space; static constexpr int dim = Dimension; @@ -101,9 +103,9 @@ class Particles using scalar_type = Cabana::MemberTypes; // no-fail. using int_type = Cabana::MemberTypes; - // v, W, rho, damage, temperature, type. + // v, W, rho, damage, type. using other_types = - Cabana::MemberTypes; + Cabana::MemberTypes; // Potentially needed later: body force (b), ID. // FIXME: add vector length. @@ -258,7 +260,6 @@ class Particles auto y = sliceCurrentPosition(); auto vol = sliceVolume(); auto nofail = sliceNoFail(); - auto temp = sliceTemperature(); // Initialize particles. auto create_functor = @@ -287,7 +288,6 @@ class Particles type( pid ) = 0; nofail( pid ) = 0; rho( pid ) = 1.0; - temp( pid ) = 0.0; return create; }; @@ -356,8 +356,8 @@ class Particles { return Cabana::slice<0>( _aosoa_vol, "volume" ); } - auto sliceType() { return Cabana::slice<5>( _aosoa_other, "type" ); } - auto sliceType() const { return Cabana::slice<5>( _aosoa_other, "type" ); } + auto sliceType() { return Cabana::slice<4>( _aosoa_other, "type" ); } + auto sliceType() const { return Cabana::slice<4>( _aosoa_other, "type" ); } auto sliceStrainEnergy() { return Cabana::slice<1>( _aosoa_other, "strain_energy" ); @@ -384,14 +384,6 @@ class Particles { return Cabana::slice<3>( _aosoa_other, "damage" ); } - auto sliceTemperature() - { - return Cabana::slice<4>( _aosoa_other, "temperature" ); - } - auto sliceTemperature() const - { - return Cabana::slice<4>( _aosoa_other, "temperature" ); - } auto sliceNoFail() { return Cabana::slice<0>( _aosoa_nofail, "no_fail_region" ); @@ -457,8 +449,7 @@ class Particles Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( h5_config, "particles", MPI_COMM_WORLD, output_step, output_time, n_local, getPosition( use_reference ), sliceStrainEnergy(), - sliceForce(), sliceDisplacement(), sliceVelocity(), sliceDamage(), - sliceTemperature() ); + sliceForce(), sliceDisplacement(), sliceVelocity(), sliceDamage() ); #else #ifdef Cabana_ENABLE_SILO Cabana::Grid::Experimental::SiloParticleOutput:: @@ -466,7 +457,7 @@ class Particles "particles", local_grid->globalGrid(), output_step, output_time, 0, n_local, getPosition( use_reference ), sliceStrainEnergy(), sliceForce(), sliceDisplacement(), sliceVelocity(), - sliceDamage(), sliceTemperature() ); + sliceDamage() ); #else log( std::cout, "No particle output enabled." ); #endif @@ -479,7 +470,8 @@ class Particles auto timeOutput() { return _output_timer.time(); }; auto time() { return _timer.time(); }; - friend class Comm; + friend class Comm; + friend class Comm; protected: aosoa_u_type _aosoa_u; @@ -501,12 +493,14 @@ class Particles }; template -class Particles - : public Particles +class Particles + : public Particles { public: - using self_type = Particles; - using base_type = Particles; + using self_type = + Particles; + using base_type = + Particles; using memory_space = typename base_type::memory_space; using base_type::dim; @@ -633,8 +627,8 @@ class Particles _output_timer.stop(); } - friend class Comm; - friend class Comm; + friend class Comm; + friend class Comm; protected: void init_lps() @@ -657,6 +651,133 @@ class Particles using base_type::_timer; }; +template +class Particles + : public Particles +{ + public: + using self_type = + Particles; + using base_type = + Particles; + using memory_space = typename base_type::memory_space; + using base_type::dim; + + // Per particle. + using base_type::n_ghost; + using base_type::n_global; + using base_type::n_local; + using base_type::size; + + // 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; + + // Per type. + using base_type::n_types; + + // Simulation total domain. + using base_type::global_mesh_ext; + + // Simulation sub domain (single MPI rank). + using base_type::ghost_mesh_hi; + using base_type::ghost_mesh_lo; + using base_type::local_mesh_ext; + using base_type::local_mesh_hi; + using base_type::local_mesh_lo; + + using base_type::dx; + using base_type::local_grid; + + using base_type::halo_width; + + // Default constructor. + Particles() + : base_type() + { + _aosoa_temp = aosoa_temp_type( "Particle Temperature", 0 ); + } + + // Constructor which initializes particles on regular grid. + template + Particles( const ExecSpace& exec_space, std::array low_corner, + std::array high_corner, + const std::array num_cells, const int max_halo_width ) + : base_type( exec_space, low_corner, high_corner, num_cells, + max_halo_width ) + { + _aosoa_temp = aosoa_temp_type( "Particle Temperature", n_local ); + init_temp(); + } + + template + void createParticles( const ExecSpace& exec_space ) + { + base_type::createParticles( exec_space ); + _aosoa_temp.resize( 0 ); + } + + auto sliceTemperature() + { + return Cabana::slice<0>( _aosoa_temp, "temperature" ); + } + auto sliceTemperature() const + { + return Cabana::slice<0>( _aosoa_temp, "temperature" ); + } + + void resize( int new_local, int new_ghost ) + { + base_type::resize( new_local, new_ghost ); + _aosoa_temp.resize( new_local + new_ghost ); + } + + void output( [[maybe_unused]] const int output_step, + [[maybe_unused]] const double output_time, + [[maybe_unused]] const bool use_reference = true ) + { +#ifdef Cabana_ENABLE_HDF5 + Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( + h5_config, "particles", MPI_COMM_WORLD, output_step, output_time, + n_local, base_type::getPosition( use_reference ), + base_type::sliceStrainEnergy(), base_type::sliceForce(), + base_type::sliceDisplacement(), base_type::sliceVelocity(), + base_type::sliceDamage(), sliceTemperature() ); +#else +#ifdef Cabana_ENABLE_SILO + Cabana::Grid::Experimental::SiloParticleOutput:: + writePartialRangeTimeStep( + "particles", local_grid->globalGrid(), output_step, output_time, + 0, n_local, base_type::getPosition( use_reference ), + base_type::sliceStrainEnergy(), base_type::sliceForce(), + base_type::sliceDisplacement(), base_type::sliceVelocity(), + base_type::sliceDamage(), sliceTemperature() ); +#else + log( std::cout, "No particle output enabled." ); +#endif +#endif + } + + friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; + + protected: + void init_temp() + { + auto temp = sliceTemperature(); + Cabana::deep_copy( temp, 0.0 ); + } + + aosoa_temp_type _aosoa_temp; + +#ifdef Cabana_ENABLE_HDF5 + using base_type::h5_config; +#endif +}; + } // namespace CabanaPD #endif diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 0515f367..1f9b5253 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -102,8 +102,8 @@ class SolverElastic using integrator_type = Integrator; using force_model_type = ForceModel; using force_type = Force; - using comm_type = - Comm; + using comm_type = Comm; using neighbor_type = Cabana::VerletList; @@ -130,6 +130,14 @@ class SolverElastic // Add ghosts from other MPI ranks. comm = std::make_shared( *particles ); + // Update temperature ghost size if needed. + if constexpr ( std::is_same::value ) + force_model.update( particles->sliceTemperature() ); + + force = + std::make_shared( inputs["half_neigh"], force_model ); + // Create the neighbor list. _neighbor_timer.start(); double mesh_min[3] = { particles->ghost_mesh_lo[0], @@ -156,9 +164,6 @@ class SolverElastic MPI_Reduce( &total_local_neighbors, &total_neighbors, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD ); - force = - std::make_shared( inputs["half_neigh"], force_model ); - print = print_rank(); if ( print ) { @@ -188,6 +193,11 @@ class SolverElastic void init_force() { + // Communicate temperature. + if constexpr ( std::is_same::value ) + comm->gatherTemperature(); + // Compute/communicate LPS weighted volume (does nothing for PMB). force->computeWeightedVolume( *particles, *neighbors, neigh_iter_tag{} ); @@ -214,6 +224,9 @@ class SolverElastic for ( int step = 1; step <= num_steps; step++ ) { _step_timer.start(); + if constexpr ( std::is_same::value ) + comm->gatherTemperature(); // Integrate - velocity Verlet first half. integrator->initialHalfStep( *particles ); diff --git a/src/force/CabanaPD_ForceModels_LPS.hpp b/src/force/CabanaPD_ForceModels_LPS.hpp index 3c0cc507..ed8d763a 100644 --- a/src/force/CabanaPD_ForceModels_LPS.hpp +++ b/src/force/CabanaPD_ForceModels_LPS.hpp @@ -23,6 +23,7 @@ struct ForceModel : public BaseForceModel using base_type = BaseForceModel; using base_model = LPS; using fracture_type = Elastic; + using thermal_type = TemperatureIndependent; using base_type::delta; @@ -67,6 +68,7 @@ struct ForceModel : public ForceModel using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = Fracture; + using thermal_type = base_type::thermal_type; using base_type::delta; using base_type::G; @@ -109,6 +111,7 @@ struct ForceModel : public ForceModel using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; + using thermal_type = base_type::thermal_type; using base_type::base_type; @@ -126,6 +129,7 @@ struct ForceModel : public ForceModel using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; + using thermal_type = base_type::thermal_type; using base_type::base_type; diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp index 2b71f599..3380de37 100644 --- a/src/force/CabanaPD_ForceModels_PMB.hpp +++ b/src/force/CabanaPD_ForceModels_PMB.hpp @@ -59,6 +59,7 @@ struct ForceModel using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = Fracture; + using thermal_type = base_type::thermal_type; using base_type::c; using base_type::delta; @@ -98,6 +99,7 @@ struct ForceModel using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; + using thermal_type = base_type::thermal_type; using base_type::base_type; @@ -113,6 +115,7 @@ struct ForceModel using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; + using thermal_type = base_type::thermal_type; using base_type::base_type; diff --git a/unit_test/tstComm.hpp b/unit_test/tstComm.hpp index afd7f907..87a73740 100644 --- a/unit_test/tstComm.hpp +++ b/unit_test/tstComm.hpp @@ -48,7 +48,9 @@ void testHalo() int halo_width = 2; // FIXME: This is for m = 1; should be calculated from m int expected_n = 6; - using particles_type = CabanaPD::Particles; + using particles_type = + CabanaPD::Particles; particles_type particles( exec_space(), box_min, box_max, num_cells, halo_width ); // Set ID equal to MPI rank. @@ -75,7 +77,9 @@ void testHalo() Cabana::deep_copy( rank_init_host, rank ); // A gather is performed on construction. - CabanaPD::Comm comm( particles ); + CabanaPD::Comm + comm( particles ); HostAoSoA aosoa_host( "host_aosoa", particles.size ); x = particles.sliceReferencePosition(); diff --git a/unit_test/tstForce.hpp b/unit_test/tstForce.hpp index be48fca7..df2c22b5 100644 --- a/unit_test/tstForce.hpp +++ b/unit_test/tstForce.hpp @@ -378,7 +378,8 @@ computeReferenceForceX( QuadraticTag, // System creation. //---------------------------------------------------------------------------// template -CabanaPD::Particles +CabanaPD::Particles createParticles( ModelType, LinearTag, const double dx, const double s0 ) { std::array box_min = { -1.0, -1.0, -1.0 }; @@ -388,7 +389,8 @@ createParticles( ModelType, LinearTag, const double dx, const double s0 ) // Create particles based on the mesh. using ptype = - CabanaPD::Particles; + CabanaPD::Particles; ptype particles( TEST_EXECSPACE{}, box_min, box_max, num_cells, 0 ); auto x = particles.sliceReferencePosition(); @@ -407,7 +409,8 @@ createParticles( ModelType, LinearTag, const double dx, const double s0 ) } template -CabanaPD::Particles +CabanaPD::Particles createParticles( ModelType, QuadraticTag, const double dx, const double s0 ) { std::array box_min = { -1.0, -1.0, -1.0 }; @@ -417,7 +420,8 @@ createParticles( ModelType, QuadraticTag, const double dx, const double s0 ) // Create particles based on the mesh. using ptype = - CabanaPD::Particles; + CabanaPD::Particles; ptype particles( TEST_EXECSPACE{}, box_min, box_max, num_cells, 0 ); auto x = particles.sliceReferencePosition(); auto u = particles.sliceDisplacement(); diff --git a/unit_test/tstIntegrator.hpp b/unit_test/tstIntegrator.hpp index 8fb41fb4..028a236b 100644 --- a/unit_test/tstIntegrator.hpp +++ b/unit_test/tstIntegrator.hpp @@ -43,8 +43,9 @@ void testIntegratorReversibility( int steps ) std::array box_max = { 1.0, 1.0, 1.0 }; std::array num_cells = { 10, 10, 10 }; - CabanaPD::Particles particles( - exec_space(), box_min, box_max, num_cells, 0 ); + CabanaPD::Particles + particles( exec_space(), box_min, box_max, num_cells, 0 ); auto x = particles.sliceReferencePosition(); std::size_t num_particle = x.size(); From 85d3bdfaeb9ec08e05eb45a3b8fe7e4ac4f713a8 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 23 May 2024 10:03:55 -0400 Subject: [PATCH 06/11] Create mechanics examples subdirectory --- examples/CMakeLists.txt | 12 +----------- examples/mechanics/CMakeLists.txt | 10 ++++++++++ examples/{ => mechanics}/crack_branching.cpp | 0 examples/{ => mechanics}/elastic_wave.cpp | 0 examples/{ => mechanics}/inputs/crack_branching.json | 0 examples/{ => mechanics}/inputs/elastic_wave.json | 0 .../{ => mechanics}/inputs/kalthoff_winkler.json | 0 examples/{ => mechanics}/kalthoff_winkler.cpp | 0 8 files changed, 11 insertions(+), 11 deletions(-) create mode 100644 examples/mechanics/CMakeLists.txt rename examples/{ => mechanics}/crack_branching.cpp (100%) rename examples/{ => mechanics}/elastic_wave.cpp (100%) rename examples/{ => mechanics}/inputs/crack_branching.json (100%) rename examples/{ => mechanics}/inputs/elastic_wave.json (100%) rename examples/{ => mechanics}/inputs/kalthoff_winkler.json (100%) rename examples/{ => mechanics}/kalthoff_winkler.cpp (100%) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index cc861300..27c4790a 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,12 +1,2 @@ -add_executable(ElasticWave elastic_wave.cpp) -target_link_libraries(ElasticWave LINK_PUBLIC CabanaPD) - -add_executable(KalthoffWinkler kalthoff_winkler.cpp) -target_link_libraries(KalthoffWinkler LINK_PUBLIC CabanaPD) - -add_executable(CrackBranching crack_branching.cpp) -target_link_libraries(CrackBranching LINK_PUBLIC CabanaPD) - -install(TARGETS ElasticWave KalthoffWinkler CrackBranching DESTINATION ${CMAKE_INSTALL_BINDIR}) - +add_subdirectory(mechanics) add_subdirectory(thermomechanics) diff --git a/examples/mechanics/CMakeLists.txt b/examples/mechanics/CMakeLists.txt new file mode 100644 index 00000000..df1b8e15 --- /dev/null +++ b/examples/mechanics/CMakeLists.txt @@ -0,0 +1,10 @@ +add_executable(ElasticWave elastic_wave.cpp) +target_link_libraries(ElasticWave LINK_PUBLIC CabanaPD) + +add_executable(KalthoffWinkler kalthoff_winkler.cpp) +target_link_libraries(KalthoffWinkler LINK_PUBLIC CabanaPD) + +add_executable(CrackBranching crack_branching.cpp) +target_link_libraries(CrackBranching LINK_PUBLIC CabanaPD) + +install(TARGETS ElasticWave KalthoffWinkler CrackBranching DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/examples/crack_branching.cpp b/examples/mechanics/crack_branching.cpp similarity index 100% rename from examples/crack_branching.cpp rename to examples/mechanics/crack_branching.cpp diff --git a/examples/elastic_wave.cpp b/examples/mechanics/elastic_wave.cpp similarity index 100% rename from examples/elastic_wave.cpp rename to examples/mechanics/elastic_wave.cpp diff --git a/examples/inputs/crack_branching.json b/examples/mechanics/inputs/crack_branching.json similarity index 100% rename from examples/inputs/crack_branching.json rename to examples/mechanics/inputs/crack_branching.json diff --git a/examples/inputs/elastic_wave.json b/examples/mechanics/inputs/elastic_wave.json similarity index 100% rename from examples/inputs/elastic_wave.json rename to examples/mechanics/inputs/elastic_wave.json diff --git a/examples/inputs/kalthoff_winkler.json b/examples/mechanics/inputs/kalthoff_winkler.json similarity index 100% rename from examples/inputs/kalthoff_winkler.json rename to examples/mechanics/inputs/kalthoff_winkler.json diff --git a/examples/kalthoff_winkler.cpp b/examples/mechanics/kalthoff_winkler.cpp similarity index 100% rename from examples/kalthoff_winkler.cpp rename to examples/mechanics/kalthoff_winkler.cpp From ef5301c01c52cf620add3625b08285576169e635 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 23 May 2024 10:19:25 -0400 Subject: [PATCH 07/11] Update readme for examples --- README.md | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index bfc983c0..168b629b 100644 --- a/README.md +++ b/README.md @@ -107,8 +107,13 @@ ctest ## Examples -Once built and installed, CabanaPD examples can be run. Timing and energy -information is output to file and particle output is written to files (if enabled within Cabana) that can be visualized with Paraview and similar applications. The first example is an elastic wave propagating through a cube from an initial Gaussian radial displacement profile from [1]. Assuming the build paths above, the example can be run with: +Once built and installed, CabanaPD `examples/` can be run. Timing and energy +information is output to file and particle output is written to files (if enabled within Cabana) that can be visualized with Paraview and similar applications. +New examples can be created by using any of the current cases as a template. All inputs are specified in the example JSON files within the relevant `inputs/` subdirectory. + +### Mechanics +Examples which only include mechanics and fracture are with `examples/mechanics`. +The first example is an elastic wave propagating through a cube from an initial Gaussian radial displacement profile from [1]. Assuming the build paths above, the example can be run with: ``` ./CabanaPD/build/install/bin/ElasticWave CabanaPD/examples/inputs/elastic_wave.json @@ -129,7 +134,14 @@ run with: ./CabanaPD/build/install/bin/CrackBranching CabanaPD/examples/inputs/crack_branching.json ``` -New examples can be created by using any of the current cases as a template. All inputs are currently specified in the example source files themselves. +### Thermomechanics +Examples which demonstrate temperature-dependent mechanics and fracture are with `examples/thermomechanics`. + +The first example demonstrates a thermoelastic problem without fracture: + +``` +./CabanaPD/build/install/bin/ThermalDeformation CabanaPD/examples/thermomechanics/thermal_deformation.json +``` ## References From f62d1c06d3dad4c333e55ac179af4af1c987c8d3 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 24 May 2024 16:07:36 -0400 Subject: [PATCH 08/11] Move temperature comm after BC update --- src/CabanaPD_Solver.hpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 1f9b5253..539a6d82 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -193,11 +193,6 @@ class SolverElastic void init_force() { - // Communicate temperature. - if constexpr ( std::is_same::value ) - comm->gatherTemperature(); - // Compute/communicate LPS weighted volume (does nothing for PMB). force->computeWeightedVolume( *particles, *neighbors, neigh_iter_tag{} ); @@ -211,7 +206,11 @@ class SolverElastic computeEnergy( *force, *particles, *neighbors, neigh_iter_tag() ); // Add boundary condition. - boundary_condition.apply( exec_space(), *particles, 0.0 ); + // boundary_condition.apply( exec_space(), *particles, 0.0 ); + // Communicate temperature. + if constexpr ( std::is_same::value ) + comm->gatherTemperature(); particles->output( 0, 0.0, output_reference ); } @@ -224,9 +223,6 @@ class SolverElastic for ( int step = 1; step <= num_steps; step++ ) { _step_timer.start(); - if constexpr ( std::is_same::value ) - comm->gatherTemperature(); // Integrate - velocity Verlet first half. integrator->initialHalfStep( *particles ); @@ -246,6 +242,10 @@ class SolverElastic // Add boundary condition. boundary_condition.apply( exec_space(), *particles, step * dt ); + if constexpr ( std::is_same::value ) + comm->gatherTemperature(); + // Integrate - velocity Verlet second half. integrator->finalHalfStep( *particles ); From 418cd7b5f5883efde13d1ef145cf0db187f06277 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Tue, 28 May 2024 10:33:24 -0400 Subject: [PATCH 09/11] Rearrange boundary conditions after solver init; Enable BC before or after force update; remove zero BC (allow no BC instead) --- examples/mechanics/crack_branching.cpp | 10 +- examples/mechanics/elastic_wave.cpp | 10 +- examples/mechanics/kalthoff_winkler.cpp | 6 +- .../thermomechanics/thermal_deformation.cpp | 38 ++-- src/CabanaPD_BodyTerm.hpp | 10 +- src/CabanaPD_Boundary.hpp | 49 ++-- src/CabanaPD_Solver.hpp | 213 ++++++++++++------ 7 files changed, 204 insertions(+), 132 deletions(-) diff --git a/examples/mechanics/crack_branching.cpp b/examples/mechanics/crack_branching.cpp index fd29321c..7e24b82d 100644 --- a/examples/mechanics/crack_branching.cpp +++ b/examples/mechanics/crack_branching.cpp @@ -110,8 +110,8 @@ void crackBranchingExample( const std::string filename ) auto sign = std::abs( ypos ) / ypos; Cabana::get( p_f, CabanaPD::Field::Force(), 1 ) += b0 * sign; }; - auto bc = - createBoundaryCondition( bc_op, exec_space{}, *particles, planes ); + auto bc = createBoundaryCondition( bc_op, exec_space{}, *particles, planes, + true ); // ==================================================== // Custom particle initialization @@ -137,9 +137,9 @@ void crackBranchingExample( const std::string filename ) // Simulation run // ==================================================== auto cabana_pd = CabanaPD::createSolverFracture( - inputs, particles, force_model, bc, prenotch ); - cabana_pd->init_force(); - cabana_pd->run(); + inputs, particles, force_model, prenotch ); + cabana_pd->init(); + cabana_pd->run( bc ); } // Initialize MPI+Kokkos. diff --git a/examples/mechanics/elastic_wave.cpp b/examples/mechanics/elastic_wave.cpp index d70904a0..27fe55a9 100644 --- a/examples/mechanics/elastic_wave.cpp +++ b/examples/mechanics/elastic_wave.cpp @@ -68,12 +68,6 @@ void elasticWaveExample( const std::string filename ) typename model_type::thermal_type>>( exec_space(), low_corner, high_corner, num_cells, halo_width ); - // ==================================================== - // Boundary conditions - // ==================================================== - auto bc = CabanaPD::createBoundaryCondition( - CabanaPD::ZeroBCTag{} ); - // ==================================================== // Custom particle initialization // ==================================================== @@ -111,8 +105,8 @@ void elasticWaveExample( const std::string filename ) // Simulation run // ==================================================== auto cabana_pd = CabanaPD::createSolverElastic( - inputs, particles, force_model, bc ); - cabana_pd->init_force(); + inputs, particles, force_model ); + cabana_pd->init(); cabana_pd->run(); // ==================================================== diff --git a/examples/mechanics/kalthoff_winkler.cpp b/examples/mechanics/kalthoff_winkler.cpp index a3b2b904..efb00ec8 100644 --- a/examples/mechanics/kalthoff_winkler.cpp +++ b/examples/mechanics/kalthoff_winkler.cpp @@ -129,9 +129,9 @@ void kalthoffWinklerExample( const std::string filename ) // Simulation run // ==================================================== auto cabana_pd = CabanaPD::createSolverFracture( - inputs, particles, force_model, bc, prenotch ); - cabana_pd->init_force(); - cabana_pd->run(); + inputs, particles, force_model, prenotch ); + cabana_pd->init(); + cabana_pd->run( bc ); } // Initialize MPI+Kokkos. diff --git a/examples/thermomechanics/thermal_deformation.cpp b/examples/thermomechanics/thermal_deformation.cpp index e25687d6..2c9cb8be 100644 --- a/examples/thermomechanics/thermal_deformation.cpp +++ b/examples/thermomechanics/thermal_deformation.cpp @@ -71,18 +71,6 @@ void thermalDeformationExample( const std::string filename ) CabanaPD::Particles>( exec_space(), low_corner, high_corner, num_cells, halo_width ); - // ==================================================== - // Imposed field - // ==================================================== - auto x = particles->sliceReferencePosition(); - auto temp = particles->sliceTemperature(); - const double low_corner_y = low_corner[1]; - auto temp_func = KOKKOS_LAMBDA( const int pid, const double t ) - { - temp( pid ) = 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; - }; - auto body_term = CabanaPD::createBodyTerm( temp_func ); - // ==================================================== // Custom particle initialization // ==================================================== @@ -98,12 +86,30 @@ void thermalDeformationExample( const std::string filename ) *particles, delta, K, alpha, temp0 ); // ==================================================== - // Simulation run + // Create solver // ==================================================== auto cabana_pd = CabanaPD::createSolverElastic( - inputs, particles, force_model, body_term ); - cabana_pd->init_force(); - cabana_pd->run(); + inputs, particles, force_model ); + + // ==================================================== + // Imposed field + // ==================================================== + auto x = particles->sliceReferencePosition(); + auto temp = particles->sliceTemperature(); + 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 ) = 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; + }; + auto body_term = CabanaPD::createBodyTerm( temp_func, false ); + + // ==================================================== + // Simulation run + // ==================================================== + cabana_pd->init( body_term ); + cabana_pd->run( body_term ); } // Initialize MPI+Kokkos. diff --git a/src/CabanaPD_BodyTerm.hpp b/src/CabanaPD_BodyTerm.hpp index ef1678c1..f156b603 100644 --- a/src/CabanaPD_BodyTerm.hpp +++ b/src/CabanaPD_BodyTerm.hpp @@ -25,11 +25,13 @@ template struct BodyTerm { UserFunctor _user_functor; + bool _force_update; Timer _timer; - BodyTerm( UserFunctor user ) + BodyTerm( UserFunctor user, const bool force ) : _user_functor( user ) + , _force_update( force ) { } @@ -47,14 +49,16 @@ struct BodyTerm _timer.stop(); } + auto forceUpdate() { return _force_update; } + auto time() { return _timer.time(); }; auto timeInit() { return 0.0; }; }; template -auto createBodyTerm( UserFunctor user_functor ) +auto createBodyTerm( UserFunctor user_functor, const bool force_update ) { - return BodyTerm( user_functor ); + return BodyTerm( user_functor, force_update ); } } // namespace CabanaPD diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp index 31e97959..48b1e466 100644 --- a/src/CabanaPD_Boundary.hpp +++ b/src/CabanaPD_Boundary.hpp @@ -141,21 +141,21 @@ struct ForceUpdateBCTag { }; -struct ZeroBCTag -{ -}; - +// Custom boundary condition. template struct BoundaryCondition { BCIndexSpace _index_space; UserFunctor _user_functor; + bool _force_update; Timer _timer; - BoundaryCondition( BCIndexSpace bc_index_space, UserFunctor user ) + BoundaryCondition( BCIndexSpace bc_index_space, UserFunctor user, + const bool force ) : _index_space( bc_index_space ) , _user_functor( user ) + , _force_update( force ) { } @@ -181,27 +181,10 @@ struct BoundaryCondition _timer.stop(); } - auto time() { return _timer.time(); }; - auto timeInit() { return _index_space.time(); }; -}; - -template -struct BoundaryCondition -{ - Timer _timer; - - template - void update( ExecSpace, Particles, RegionBoundary ) - { - } - - template - void apply( ExecSpace, ParticleType, double ) - { - } + auto forceUpdate() { return _force_update; } auto time() { return _timer.time(); }; - auto timeInit() { return 0.0; }; + auto timeInit() { return _index_space.time(); }; }; template @@ -209,6 +192,7 @@ struct BoundaryCondition { double _value; BCIndexSpace _index_space; + const bool _force_update = true; Timer _timer; @@ -242,6 +226,8 @@ struct BoundaryCondition _timer.stop(); } + auto forceUpdate() { return _force_update; } + auto time() { return _timer.time(); }; auto timeInit() { return _index_space.time(); }; }; @@ -251,6 +237,7 @@ struct BoundaryCondition { double _value; BCIndexSpace _index_space; + const bool _force_update = true; Timer _timer; @@ -286,6 +273,8 @@ struct BoundaryCondition _timer.stop(); } + auto forceUpdate() { return _force_update; } + auto time() { return _timer.time(); }; auto timeInit() { return _index_space.time(); }; }; @@ -310,21 +299,15 @@ template planes, + const bool force_update, const double initial_guess = 0.5 ) { using memory_space = typename Particles::memory_space; using bc_index_type = BoundaryIndexSpace; bc_index_type bc_indices = createBoundaryIndexSpace( exec_space, particles, planes, initial_guess ); - return BoundaryCondition( bc_indices, - user_functor ); -} - -template -auto createBoundaryCondition( ZeroBCTag ) -{ - using bc_index_type = BoundaryIndexSpace; - return BoundaryCondition(); + return BoundaryCondition( + bc_indices, user_functor, force_update ); } } // namespace CabanaPD diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 539a6d82..2e60ad8d 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -91,7 +91,7 @@ class SolverBase }; template + class ForceModel> class SolverElastic { public: @@ -109,14 +109,12 @@ class SolverElastic Cabana::VerletLayout2D, Cabana::TeamOpTag>; using neigh_iter_tag = Cabana::SerialOpTag; using input_type = InputType; - using bc_type = BoundaryCondition; SolverElastic( input_type _inputs, std::shared_ptr _particles, - force_model_type force_model, bc_type bc ) + force_model_type force_model ) : inputs( _inputs ) , particles( _particles ) - , boundary_condition( bc ) , _init_time( 0.0 ) { num_steps = inputs["num_steps"]; @@ -191,7 +189,7 @@ class SolverElastic _init_timer.stop(); } - void init_force() + void init( const bool initial_output = true ) { // Compute/communicate LPS weighted volume (does nothing for PMB). force->computeWeightedVolume( *particles, *neighbors, @@ -205,14 +203,84 @@ class SolverElastic computeForce( *force, *particles, *neighbors, neigh_iter_tag{} ); computeEnergy( *force, *particles, *neighbors, neigh_iter_tag() ); + if ( initial_output ) + particles->output( 0, 0.0, output_reference ); + } + + template + void init( BoundaryType boundary_condition, + const bool initial_output = true ) + { + if ( !boundary_condition.forceUpdate() ) + { + // Add boundary condition. + boundary_condition.apply( exec_space(), *particles, 0.0 ); + + // Communicate temperature. + if constexpr ( std::is_same::value ) + comm->gatherTemperature(); + } + + // Force init without particle output. + init( false ); + // Add boundary condition. - // boundary_condition.apply( exec_space(), *particles, 0.0 ); - // Communicate temperature. - if constexpr ( std::is_same::value ) - comm->gatherTemperature(); + if ( boundary_condition.forceUpdate() ) + boundary_condition.apply( exec_space(), *particles, 0.0 ); + + if ( initial_output ) + particles->output( 0, 0.0, output_reference ); + } + + template + void run( BoundaryType boundary_condition ) + { + init_output( boundary_condition.timeInit() ); + + // Main timestep loop. + for ( int step = 1; step <= num_steps; step++ ) + { + _step_timer.start(); + + // 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 ); + + if constexpr ( std::is_same< + typename force_model_type::thermal_type, + TemperatureDependent>::value ) + comm->gatherTemperature(); + } + + // Update ghost particles. + comm->gatherDisplacement(); + + // Do not need to recompute LPS weighted volume here without damage. + // Compute/communicate LPS dilatation (does nothing for PMB). + force->computeDilatation( *particles, *neighbors, + neigh_iter_tag{} ); + comm->gatherDilatation(); + + // Compute internal forces. + computeForce( *force, *particles, *neighbors, neigh_iter_tag{} ); - particles->output( 0, 0.0, output_reference ); + // Add force boundary condition. + if ( boundary_condition.forceUpdate() ) + boundary_condition.apply( exec_space(), *particles, step * dt ); + + // Integrate - velocity Verlet second half. + integrator->finalHalfStep( *particles ); + + output( step ); + } + + // Final output and timings. + final_output(); } void run() @@ -239,9 +307,6 @@ class SolverElastic // Compute internal forces. computeForce( *force, *particles, *neighbors, neigh_iter_tag{} ); - // Add boundary condition. - boundary_condition.apply( exec_space(), *particles, step * dt ); - if constexpr ( std::is_same::value ) comm->gatherTemperature(); @@ -249,34 +314,39 @@ class SolverElastic // Integrate - velocity Verlet second half. integrator->finalHalfStep( *particles ); - // Print output. - if ( step % output_frequency == 0 ) - { - auto W = computeEnergy( *force, *particles, *neighbors, - neigh_iter_tag() ); - - particles->output( step / output_frequency, step * dt, - output_reference ); - _step_timer.stop(); - step_output( step, W ); - } - else - { - _step_timer.stop(); - } + output( step ); } // Final output and timings. final_output(); } - void init_output() + void output( const int step ) + { + // Print output. + if ( step % output_frequency == 0 ) + { + auto W = computeEnergy( *force, *particles, *neighbors, + neigh_iter_tag() ); + + particles->output( step / output_frequency, step * dt, + output_reference ); + _step_timer.stop(); + step_output( step, W ); + } + else + { + _step_timer.stop(); + } + } + + void init_output( double boundary_init_time = 0.0 ) { // Output after construction and initial forces. std::ofstream out( output_file, std::ofstream::app ); _init_time += _init_timer.time() + _neighbor_timer.time() + - particles->timeInit() + boundary_condition.timeInit() + - comm->timeInit() + integrator->timeInit(); + particles->timeInit() + comm->timeInit() + + integrator->timeInit() + boundary_init_time; log( out, "Init-Time(s): ", _init_time ); log( out, "Init-Neighbor-Time(s): ", _neighbor_timer.time(), "\n" ); log( out, "#Timestep/Total-steps Simulation-time Total-strain-energy " @@ -358,7 +428,6 @@ class SolverElastic std::shared_ptr integrator; std::shared_ptr force; std::shared_ptr neighbors; - bc_type boundary_condition; std::string output_file; std::string error_file; @@ -373,14 +442,13 @@ class SolverElastic }; template + class ForceModel, class PrenotchType> class SolverFracture - : public SolverElastic + : public SolverElastic { public: - using base_type = SolverElastic; + using base_type = + SolverElastic; using exec_space = typename base_type::exec_space; using memory_space = typename base_type::memory_space; @@ -391,15 +459,13 @@ class SolverFracture using force_model_type = ForceModel; using force_type = typename base_type::force_type; using neigh_iter_tag = Cabana::SerialOpTag; - using bc_type = BoundaryCondition; using prenotch_type = PrenotchType; using input_type = typename base_type::input_type; SolverFracture( input_type _inputs, std::shared_ptr _particles, - force_model_type force_model, bc_type bc, - prenotch_type prenotch ) - : base_type( _inputs, _particles, force_model, bc ) + force_model_type force_model, prenotch_type prenotch ) + : base_type( _inputs, _particles, force_model ) { init_mu(); @@ -410,8 +476,8 @@ class SolverFracture SolverFracture( input_type _inputs, std::shared_ptr _particles, - force_model_type force_model, bc_type bc ) - : base_type( _inputs, _particles, force_model, bc ) + force_model_type force_model ) + : base_type( _inputs, _particles, force_model ) { init_mu(); } @@ -429,7 +495,7 @@ class SolverFracture _init_timer.stop(); } - void init_force() + void init( const bool initial_output = true ) { // Compute/communicate weighted volume for LPS (does nothing for PMB). force->computeWeightedVolume( *particles, *neighbors, mu ); @@ -442,15 +508,31 @@ class SolverFracture computeForce( *force, *particles, *neighbors, mu, neigh_iter_tag{} ); computeEnergy( *force, *particles, *neighbors, mu, neigh_iter_tag() ); - // Add boundary condition. - boundary_condition.apply( exec_space(), *particles, 0 ); + if ( initial_output ) + particles->output( 0, 0.0, output_reference ); + } + + template + void init( BoundaryType boundary_condition, + const bool initial_output = true ) + { + if ( !boundary_condition.forceUpdate() ) + boundary_condition.apply( exec_space(), *particles, 0.0 ); + + // Force init without particle output. + init( false ); + + if ( boundary_condition.forceUpdate() ) + boundary_condition.apply( exec_space(), *particles, 0.0 ); - particles->output( 0, 0.0, output_reference ); + if ( initial_output ) + particles->output( 0, 0.0, output_reference ); } - void run() + template + void run( BoundaryType boundary_condition ) { - this->init_output(); + this->init_output( boundary_condition.timeInit() ); // Main timestep loop. for ( int step = 1; step <= num_steps; step++ ) @@ -460,6 +542,10 @@ class SolverFracture // 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(); @@ -474,8 +560,9 @@ class SolverFracture computeForce( *force, *particles, *neighbors, mu, neigh_iter_tag{} ); - // Add boundary condition. - boundary_condition.apply( exec_space{}, *particles, step * dt ); + // Add force boundary condition. + if ( boundary_condition.forceUpdate() ) + boundary_condition.apply( exec_space{}, *particles, step * dt ); // Integrate - velocity Verlet second half. integrator->finalHalfStep( *particles ); @@ -507,7 +594,6 @@ class SolverFracture using base_type::output_reference; protected: - using base_type::boundary_condition; using base_type::comm; using base_type::force; using base_type::inputs; @@ -525,26 +611,25 @@ class SolverFracture }; template + class ForceModel> auto createSolverElastic( InputsType inputs, std::shared_ptr particles, - ForceModel model, BCType bc ) + ForceModel model ) { - return std::make_shared>( - inputs, particles, model, bc ); + return std::make_shared< + SolverElastic>( + inputs, particles, model ); } template + class ForceModel, class PrenotchType> auto createSolverFracture( InputsType inputs, std::shared_ptr particles, - ForceModel model, BCType bc, PrenotchType prenotch ) + ForceModel model, PrenotchType prenotch ) { - return std::make_shared< - SolverFracture>( inputs, particles, model, bc, - prenotch ); + return std::make_shared>( + inputs, particles, model, prenotch ); } } // namespace CabanaPD From c501fcdf32bc183361c451523bae1f9574ba7049 Mon Sep 17 00:00:00 2001 From: pabloseleson Date: Tue, 28 May 2024 11:53:41 -0400 Subject: [PATCH 10/11] Added reference to thermal deformation problem --- README.md | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 168b629b..f8aecea3 100644 --- a/README.md +++ b/README.md @@ -113,6 +113,7 @@ New examples can be created by using any of the current cases as a template. All ### Mechanics Examples which only include mechanics and fracture are with `examples/mechanics`. + The first example is an elastic wave propagating through a cube from an initial Gaussian radial displacement profile from [1]. Assuming the build paths above, the example can be run with: ``` @@ -137,7 +138,8 @@ run with: ### Thermomechanics Examples which demonstrate temperature-dependent mechanics and fracture are with `examples/thermomechanics`. -The first example demonstrates a thermoelastic problem without fracture: +The first example demonstrates a thermoelastic problem without fracture with a homogeneous plate +under linear thermal loading [4]. The example can be run with: ``` ./CabanaPD/build/install/bin/ThermalDeformation CabanaPD/examples/thermomechanics/thermal_deformation.json @@ -157,6 +159,8 @@ Kunze, and L.W. Meyer, eds., Vol 1, DGM Informationsgesellschaft Verlag (1988) [3] F. Bobaru and G. Zhang, Why do cracks branch? A peridynamic investigation of dynamic brittle fracture, International Journal of Fracture 196 (2015): 59–98. +[4] D. He, D. Huang, and D. Jiang, Modeling and studies of fracture in functionally graded materials under thermal shock loading using peridynamics, Theoretical and Applied Fracture Mechanics 111 (2021): 102852. + ## Contributing We encourage you to contribute to CabanaPD! Please check the From e5c395f906e84d669b33b10e9fcc97cfa6917a46 Mon Sep 17 00:00:00 2001 From: pabloseleson Date: Tue, 28 May 2024 18:42:14 -0400 Subject: [PATCH 11/11] Changes to thermal deformation example --- .../inputs/thermal_deformation.json | 6 +++--- examples/thermomechanics/thermal_deformation.cpp | 15 +++++++++++++-- 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/examples/thermomechanics/inputs/thermal_deformation.json b/examples/thermomechanics/inputs/thermal_deformation.json index e6bc5205..4f0f6908 100644 --- a/examples/thermomechanics/inputs/thermal_deformation.json +++ b/examples/thermomechanics/inputs/thermal_deformation.json @@ -1,12 +1,12 @@ { - "num_cells" : {"value": [101, 31, 3]}, + "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": 0.0, "unit": "oC"}, + "reference_temperature" : {"value": 20.0, "unit": "oC"}, "horizon" : {"value": 0.03, "unit": "m"}, - "final_time" : {"value": 0.0093, "unit": "s"}, + "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/thermal_deformation.cpp b/examples/thermomechanics/thermal_deformation.cpp index 2c9cb8be..4dea6829 100644 --- a/examples/thermomechanics/thermal_deformation.cpp +++ b/examples/thermomechanics/thermal_deformation.cpp @@ -22,7 +22,7 @@ void thermalDeformationExample( const std::string filename ) { // ==================================================== - // Use default Kokkos spaces + // Use default Kokkos spaces // ==================================================== using exec_space = Kokkos::DefaultExecutionSpace; using memory_space = typename exec_space::memory_space; @@ -101,7 +101,7 @@ void thermalDeformationExample( const std::string filename ) // particles are correctly taken into account for lambda capture here. auto temp_func = KOKKOS_LAMBDA( const int pid, const double t ) { - temp( pid ) = 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; + temp( pid ) = temp0 + 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; }; auto body_term = CabanaPD::createBodyTerm( temp_func, false ); @@ -110,6 +110,17 @@ void thermalDeformationExample( const std::string filename ) // ==================================================== cabana_pd->init( body_term ); cabana_pd->run( body_term ); + + // ==================================================== + // Outputs + // ==================================================== + // Output displacement along the x-axis + createDisplacementProfile( MPI_COMM_WORLD, num_cells[0], 0, + "xdisplacement_profile.txt", *particles ); + + // Output displacement along the y-axis + createDisplacementProfile( MPI_COMM_WORLD, num_cells[1], 1, + "ydisplacement_profile.txt", *particles ); } // Initialize MPI+Kokkos.