From 508f3976b8c506bef515060e4f38379d09275325 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 31 Jan 2025 12:40:29 -0500 Subject: [PATCH 01/30] Create Neighbor class and test; stores list, handles iteration, and broken bonds --- src/CabanaPD_Force.hpp | 157 ++---------------- src/CabanaPD_ForceModels.hpp | 6 +- src/CabanaPD_HeatTransfer.hpp | 81 ++++------ src/CabanaPD_Neighbor.hpp | 225 ++++++++++++++++++++++++++ src/CabanaPD_Solver.hpp | 57 ++++--- src/CabanaPD_Types.hpp | 6 + src/force/CabanaPD_Contact.hpp | 101 ++---------- src/force/CabanaPD_Hertzian.hpp | 34 ++-- src/force/CabanaPD_LPS.hpp | 201 +++++++++-------------- src/force/CabanaPD_PMB.hpp | 147 ++++++----------- src/force_models/CabanaPD_Contact.hpp | 5 + unit_test/tstForce.hpp | 61 +++++-- 12 files changed, 503 insertions(+), 578 deletions(-) create mode 100644 src/CabanaPD_Neighbor.hpp diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index 05aef53f..458b7334 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -63,6 +63,7 @@ #include #include +#include #include namespace CabanaPD @@ -130,160 +131,33 @@ template class Force; template -class Force +class BaseForce { - public: - using neighbor_list_type = - Cabana::VerletList; - protected: - bool _half_neigh; - neighbor_list_type _neigh_list; - Timer _timer; Timer _energy_timer; public: - // Primary constructor: use positions and construct neighbors. - template - Force( const bool half_neigh, const double delta, - const ParticleType& particles, const double tol = 1e-14 ) - : _half_neigh( half_neigh ) - , _neigh_list( neighbor_list_type( - particles.sliceReferencePosition(), particles.frozenOffset(), - particles.localOffset(), delta + tol, 1.0, - particles.ghost_mesh_lo, particles.ghost_mesh_hi ) ) - { - } - - // General constructor (necessary for contact, but could be used by any - // force routine). - template - Force( const bool half_neigh, const double delta, - const PositionType& positions, const std::size_t frozen_offset, - const std::size_t local_offset, const double mesh_min[3], - const double mesh_max[3], const double tol = 1e-14 ) - : _half_neigh( half_neigh ) - , _neigh_list( neighbor_list_type( positions, frozen_offset, - local_offset, delta + tol, 1.0, - mesh_min, mesh_max ) ) - { - } - - // Constructor which stores existing neighbors. - template - Force( const bool half_neigh, const NeighborListType& neighbors ) - : _half_neigh( half_neigh ) - , _neigh_list( neighbors ) - { - } - - // FIXME: should it be possible to update this list? - template - void update( const ParticleType&, const double, const bool = false ) - { - } - - unsigned getMaxLocalNeighbors() - { - auto neigh = _neigh_list; - unsigned local_max_neighbors; - auto neigh_max = KOKKOS_LAMBDA( const int, unsigned& max_n ) - { - max_n = - Cabana::NeighborList::maxNeighbor( neigh ); - }; - using exec_space = typename MemorySpace::execution_space; - Kokkos::RangePolicy policy( 0, 1 ); - Kokkos::parallel_reduce( policy, neigh_max, local_max_neighbors ); - Kokkos::fence(); - return local_max_neighbors; - } - - void getNeighborStatistics( unsigned& max_neighbors, - unsigned long long& total_neighbors ) - { - auto neigh = _neigh_list; - unsigned local_max_neighbors; - unsigned long long local_total_neighbors; - auto neigh_stats = KOKKOS_LAMBDA( const int, unsigned& max_n, - unsigned long long& total_n ) - { - max_n = - Cabana::NeighborList::maxNeighbor( neigh ); - total_n = Cabana::NeighborList::totalNeighbor( - neigh ); - }; - using exec_space = typename MemorySpace::execution_space; - Kokkos::RangePolicy policy( 0, 1 ); - Kokkos::parallel_reduce( policy, neigh_stats, local_max_neighbors, - local_total_neighbors ); - Kokkos::fence(); - MPI_Reduce( &local_max_neighbors, &max_neighbors, 1, MPI_UNSIGNED, - MPI_MAX, 0, MPI_COMM_WORLD ); - MPI_Reduce( &local_total_neighbors, &total_neighbors, 1, - MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD ); - } - - // Default to no-op for pair models. - template - void computeWeightedVolume( ParticleType&, const ParallelType ) const + // Default to no-op. + template + void computeWeightedVolume( ParticleType&, const NeighborType ) const { } - template - void computeDilatation( ParticleType&, const ParallelType ) const + template + void computeDilatation( ParticleType&, const NeighborType ) const { } - auto getNeighbors() const { return _neigh_list; } - auto time() { return _timer.time(); }; auto timeEnergy() { return _energy_timer.time(); }; - auto timeNeighbor() { return 0.0; }; -}; - -template -class BaseFracture -{ - protected: - using memory_space = MemorySpace; - using NeighborView = typename Kokkos::View; - NeighborView _mu; - - public: - BaseFracture( const int local_particles, const int max_neighbors ) - { - // Create View to track broken bonds. - // TODO: this could be optimized to ignore frozen particle bonds. - _mu = NeighborView( - Kokkos::ViewAllocateWithoutInitializing( "broken_bonds" ), - local_particles, max_neighbors ); - Kokkos::deep_copy( _mu, 1 ); - } - - BaseFracture( NeighborView mu ) - : _mu( mu ) - { - } - - template - void prenotch( ExecSpace exec_space, const ParticleType& particles, - PrenotchType& prenotch, NeighborList& neigh_list ) - { - prenotch.create( exec_space, _mu, particles, neigh_list ); - } - - auto getBrokenBonds() const { return _mu; } }; /****************************************************************************** Force free functions. ******************************************************************************/ -template +template void computeForce( ForceType& force, ParticleType& particles, - const ParallelType& neigh_op_tag, const bool reset = true ) + NeighborType& neighbor, const bool reset = true ) { auto x = particles.sliceReferencePosition(); auto u = particles.sliceDisplacement(); @@ -300,16 +174,17 @@ void computeForce( ForceType& force, ParticleType& particles, // neigh_op_tag ); // Forces only atomic if using team threading. - if ( std::is_same::value ) - force.computeForceFull( f_a, x, u, particles, neigh_op_tag ); + if constexpr ( std::is_same::value ) + force.computeForceFull( f_a, x, u, particles, neighbor ); else - force.computeForceFull( f, x, u, particles, neigh_op_tag ); + force.computeForceFull( f, x, u, particles, neighbor ); Kokkos::fence(); } -template +template double computeEnergy( ForceType& force, ParticleType& particles, - const ParallelType& neigh_op_tag ) + const NeighborType& neighbor ) { double energy = 0.0; if constexpr ( is_energy_output::value ) @@ -327,7 +202,7 @@ double computeEnergy( ForceType& force, ParticleType& particles, // energy = computeEnergy_half( force, x, u, // neigh_op_tag ); // else - energy = force.computeEnergyFull( W, x, u, particles, neigh_op_tag ); + energy = force.computeEnergyFull( W, x, u, particles, neighbor ); Kokkos::fence(); } return energy; diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index 6717cf05..450105a4 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -24,6 +24,8 @@ struct BaseForceModel BaseForceModel( const double _delta ) : delta( _delta ){}; + auto cutoff() const { return delta; } + // Only needed for models which store bond properties. void updateBonds( const int, const int ) {} @@ -37,8 +39,8 @@ class BasePlasticity { protected: using memory_space = MemorySpace; - using NeighborView = typename Kokkos::View; - NeighborView _s_p; + using neighbor_view = typename Kokkos::View; + neighbor_view _s_p; public: // Must update later because number of neighbors not known at construction. diff --git a/src/CabanaPD_HeatTransfer.hpp b/src/CabanaPD_HeatTransfer.hpp index 5111fa49..7835224f 100644 --- a/src/CabanaPD_HeatTransfer.hpp +++ b/src/CabanaPD_HeatTransfer.hpp @@ -13,6 +13,7 @@ #define HEATTRANSFER_H #include +#include #include #include @@ -27,39 +28,31 @@ class HeatTransfer; template class HeatTransfer> - : public Force { public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; using model_type = ForceModel; - using base_type = Force; - using neighbor_list_type = typename base_type::neighbor_list_type; protected: - using base_type::_half_neigh; - using base_type::_neigh_list; - using base_type::_timer; - Timer _euler_timer = base_type::_energy_timer; + Timer _timer; + Timer _euler_timer; model_type _model; public: // Running with mechanics as well; no reason to rebuild neighbors. - template - HeatTransfer( const bool half_neigh, const ForceType& force, - const model_type model ) - : base_type( half_neigh, force.getNeighbors() ) - , _model( model ) + HeatTransfer( const model_type model ) + : _model( model ) { } template + class NeighborType> void computeHeatTransferFull( TemperatureType& conduction, const PosType& x, const PosType& u, const ParticleType& particles, - ParallelType& neigh_op_tag ) + NeighborType& neighbor ) { _timer.start(); @@ -77,12 +70,8 @@ class HeatTransfer policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, temp_func, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::HeatTransfer::computeFull" ); - + neighbor.iterate( exec_space{}, temp_func, particles, + "CabanaPD::HeatTransfer::computeFull" ); _timer.stop(); } @@ -111,9 +100,7 @@ class HeatTransfer> : public HeatTransfer>, - BaseFracture - + DynamicTemperature, ModelParams...>> { public: // Using the default exec_space. @@ -124,47 +111,39 @@ class HeatTransfer>; - using neighbor_list_type = typename base_type::neighbor_list_type; protected: using base_type::_euler_timer; - using base_type::_half_neigh; - using base_type::_neigh_list; using base_type::_timer; model_type _model; - using fracture_type = BaseFracture; - using fracture_type::_mu; - public: // Explicit base model construction is necessary because of the indirect // model inheritance. This could be avoided with a BaseHeatTransfer object. - template - HeatTransfer( const bool half_neigh, const ForceType& force, - const model_type model ) - : base_type( half_neigh, force, - typename base_type::model_type( - PMB{}, NoFracture{}, model.delta, model.K, - model.temperature, model.kappa, model.cp, model.alpha, - model.temp0, model.constant_microconductivity ) ) - , fracture_type( force.getBrokenBonds() ) + HeatTransfer( const model_type model ) + : base_type( typename base_type::model_type( + PMB{}, NoFracture{}, model.delta, model.K, model.temperature, + model.kappa, model.cp, model.alpha, model.temp0, + model.constant_microconductivity ) ) , _model( model ) { } template + class NeighborType> void computeHeatTransferFull( TemperatureType& conduction, const PosType& x, const PosType& u, - const ParticleType& particles, ParallelType& ) + const ParticleType& particles, + NeighborType& neighbor ) { _timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); auto model = _model; - const auto neigh_list = _neigh_list; const auto vol = particles.sliceVolume(); const auto temp = particles.sliceTemperature(); - const auto mu = _mu; auto temp_func = KOKKOS_LAMBDA( const int i ) { @@ -191,19 +170,17 @@ class HeatTransfer policy( particles.frozenOffset(), - particles.localOffset() ); - Kokkos::parallel_for( "CabanaPD::HeatTransfer::computeFull", policy, - temp_func ); + neighbor.iterateLinear( exec_space{}, temp_func, particles, + "CabanaPD::HeatTransfer::computeFull" ); _timer.stop(); } }; // Heat transfer free functions. -template +template void computeHeatTransfer( HeatTransferType& heat_transfer, - ParticleType& particles, - const ParallelType& neigh_op_tag, const double dt ) + ParticleType& particles, const NeighborType& neighbor, + const double dt ) { auto x = particles.sliceReferencePosition(); auto u = particles.sliceDisplacement(); @@ -214,12 +191,12 @@ void computeHeatTransfer( HeatTransferType& heat_transfer, Cabana::deep_copy( conduction, 0.0 ); // Temperature only needs to be atomic if using team threading. - if ( std::is_same::value ) + if ( std::is_same::value ) heat_transfer.computeHeatTransferFull( conduction_a, x, u, particles, - neigh_op_tag ); + neighbor ); else heat_transfer.computeHeatTransferFull( conduction, x, u, particles, - neigh_op_tag ); + neighbor ); Kokkos::fence(); heat_transfer.forwardEuler( particles, dt ); diff --git a/src/CabanaPD_Neighbor.hpp b/src/CabanaPD_Neighbor.hpp new file mode 100644 index 00000000..9c4efc41 --- /dev/null +++ b/src/CabanaPD_Neighbor.hpp @@ -0,0 +1,225 @@ +/**************************************************************************** + * 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 * + ****************************************************************************/ + +#ifndef NEIGHBOR_H +#define NEIGHBOR_H + +#include + +#include + +namespace CabanaPD +{ +template +class Neighbor; + +template +class Neighbor +{ + public: + using list_type = + Cabana::VerletList; + using Tag = Cabana::SerialOpTag; + Tag tag; + + protected: + using exec_space = typename MemorySpace::execution_space; + bool _half_neigh; + list_type _neigh_list; + double mesh_max[3]; + double mesh_min[3]; + + Timer _timer; + + public: + // Primary constructor: use positions and construct neighbors. + template + Neighbor( const bool half_neigh, const ModelType& model, + const ParticleType& particles, const double tol = 1e-14 ) + : _half_neigh( half_neigh ) + { + _timer.start(); + if constexpr ( !is_contact::value ) + _neigh_list = list_type( + particles.sliceReferencePosition(), particles.frozenOffset(), + particles.localOffset(), model.cutoff() + tol, 1.0, + particles.ghost_mesh_lo, particles.ghost_mesh_hi ); + else + _neigh_list = list_type( + particles.sliceCurrentPosition(), particles.frozenOffset(), + particles.localOffset(), model.cutoff() + tol, 1.0, + particles.ghost_mesh_lo, particles.ghost_mesh_hi ); + _timer.stop(); + } + + unsigned getMaxLocal() + { + auto neigh = _neigh_list; + unsigned local_max_neighbors; + auto neigh_max = KOKKOS_LAMBDA( const int, unsigned& max_n ) + { + max_n = Cabana::NeighborList::maxNeighbor( neigh ); + }; + Kokkos::RangePolicy policy( 0, 1 ); + Kokkos::parallel_reduce( policy, neigh_max, local_max_neighbors ); + Kokkos::fence(); + return local_max_neighbors; + } + + void getStatistics( unsigned& max_neighbors, + unsigned long long& total_neighbors ) + { + auto neigh = _neigh_list; + unsigned local_max_neighbors; + unsigned long long local_total_neighbors; + auto neigh_stats = KOKKOS_LAMBDA( const int, unsigned& max_n, + unsigned long long& total_n ) + { + max_n = Cabana::NeighborList::maxNeighbor( neigh ); + total_n = Cabana::NeighborList::totalNeighbor( neigh ); + }; + Kokkos::RangePolicy policy( 0, 1 ); + Kokkos::parallel_reduce( policy, neigh_stats, local_max_neighbors, + local_total_neighbors ); + Kokkos::fence(); + MPI_Reduce( &local_max_neighbors, &max_neighbors, 1, MPI_UNSIGNED, + MPI_MAX, 0, MPI_COMM_WORLD ); + MPI_Reduce( &local_total_neighbors, &total_neighbors, 1, + MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD ); + } + + auto list() const { return _neigh_list; } + + template + void update( ParticleType& particles, const double cutoff ) + { + const auto y = particles.sliceCurrentPosition(); + const int n_frozen = particles.frozenOffset(); + const int n_local = particles.localOffset(); + + _neigh_list.build( exec_space{}, y, n_frozen, n_local, cutoff, 1.0, + particles.ghost_mesh_lo, particles.ghost_mesh_hi ); + } + + // Only rebuild neighbor list as needed. + template + void update( const ParticleType& particles, const double search_radius, + const double radius_extend, const bool require_update = false ) + { + double max_displacement = particles.getMaxDisplacement(); + if ( max_displacement > radius_extend || require_update ) + { + _timer.start(); + const auto y = particles.sliceCurrentPosition(); + _neigh_list.build( y, particles.frozenOffset(), + particles.localOffset(), search_radius, 1.0, + mesh_min, mesh_max ); + // Reset neighbor update displacement. + const auto u = particles.sliceDisplacement(); + auto u_neigh = particles.sliceDisplacementNeighborBuild(); + Cabana::deep_copy( u_neigh, u ); + _timer.stop(); + } + } + + template + void iterate( ExecSpace, const FunctorType functor, + const ParticleType& particles, const std::string label ) const + { + auto policy = makePolicy( particles ); + Cabana::neighbor_parallel_for( policy, functor, _neigh_list, + Cabana::FirstNeighborsTag(), tag, + label ); + } + + template + void iterateLinear( ExecSpace, const FunctorType functor, + const ParticleType& particles, + const std::string label ) const + { + auto policy = makePolicy( particles ); + Kokkos::parallel_for( label, policy, functor ); + } + + template + auto reduce( ExecSpace, const FunctorType functor, + const ParticleType& particles, const std::string label ) const + { + auto policy = makePolicy( particles ); + double total = 0.0; + Cabana::neighbor_parallel_reduce( policy, functor, _neigh_list, + Cabana::FirstNeighborsTag(), tag, + total, label ); + return total; + } + + template + auto reduceLinear( ExecSpace, const FunctorType functor, + const ParticleType& particles, + const std::string label ) const + { + double total = 0.0; + auto policy = makePolicy( particles ); + Kokkos::parallel_reduce( label, policy, functor, total ); + return total; + } + + template + auto makePolicy( const ParticleType& particles ) const + { + return Kokkos::RangePolicy( particles.frozenOffset(), + particles.localOffset() ); + } + + auto time() { return _timer.time(); }; +}; + +template +class Neighbor : public Neighbor +{ + public: + using memory_space = MemorySpace; + using neighbor_view = typename Kokkos::View; + + protected: + using base_type = Neighbor; + using base_type::_neigh_list; + + neighbor_view _mu; + + public: + template + Neighbor( const bool half_neigh, const ModelType& model, + const ParticleType& particles, const double tol = 1e-14 ) + : base_type( half_neigh, model, particles, tol ) + { + // Create View to track broken bonds. + // TODO: this could be optimized to ignore frozen particle bonds. + _mu = neighbor_view( + Kokkos::ViewAllocateWithoutInitializing( "broken_bonds" ), + particles.localOffset(), base_type::getMaxLocal() ); + Kokkos::deep_copy( _mu, 1 ); + } + + template + void prenotch( ExecSpace exec_space, const ParticleType& particles, + PrenotchType& prenotch ) + { + prenotch.create( exec_space, _mu, particles, _neigh_list ); + } + + auto brokenBonds() const { return _mu; } +}; + +} // namespace CabanaPD + +#endif diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 930c12d1..8ae80b42 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -99,7 +99,8 @@ class Solver using comm_type = Comm; - using neigh_iter_tag = Cabana::SerialOpTag; + using neighbor_type = + Neighbor; // Optional module types. using heat_transfer_type = HeatTransfer; @@ -130,10 +131,7 @@ class Solver { setup( force_model ); - _neighbor_timer.start(); - contact = std::make_shared( inputs["half_neigh"], - particles, contact_model ); - _neighbor_timer.stop(); + contact = std::make_shared( contact_model ); } void setup( force_model_type force_model ) @@ -165,16 +163,19 @@ class Solver typename force_model_type::thermal_type>::value ) force_model.update( particles.sliceTemperature() ); - _neighbor_timer.start(); + neighbor = std::make_shared( inputs["half_neigh"], + force_model, particles ); + // Plastic models need correct size for the bond array. + force_model.updateBonds( particles.localOffset(), + neighbor->getMaxLocal() ); + // This will either be PD or DEM forces. - force = std::make_shared( inputs["half_neigh"], particles, - force_model ); - _neighbor_timer.stop(); + force = std::make_shared( force_model ); _init_timer.start(); unsigned max_neighbors; unsigned long long total_neighbors; - force->getNeighborStatistics( max_neighbors, total_neighbors ); + neighbor->getStatistics( max_neighbors, total_neighbors ); // Create heat transfer if needed, using the same neighbor list as // the mechanics. @@ -182,8 +183,7 @@ class Solver typename force_model_type::thermal_type>::value ) { thermal_subcycle_steps = inputs["thermal_subcycle_steps"]; - heat_transfer = std::make_shared( - inputs["half_neigh"], *force, force_model ); + heat_transfer = std::make_shared( force_model ); } print = print_rank(); @@ -221,12 +221,12 @@ class Solver if constexpr ( !is_fracture< typename force_model_type::fracture_type>::value ) { - force->computeWeightedVolume( particles, neigh_iter_tag{} ); + force->computeWeightedVolume( particles, *neighbor ); comm->gatherWeightedVolume(); } // Compute initial internal forces and energy. updateForce(); - computeEnergy( *force, particles, neigh_iter_tag() ); + computeEnergy( *force, particles, *neighbor ); if ( initial_output ) particles.output( 0, 0.0, output_reference ); @@ -304,7 +304,7 @@ class Solver // FIXME: Will need to rebuild ghosts. } - void updateNeighbors() { force->update( *particles, 0.0, true ); } + void updateNeighbors() { force->update( particles, 0.0, true ); } template void run( BoundaryType boundary_condition ) @@ -326,8 +326,7 @@ class Solver typename force_model_type::thermal_type>::value ) { if ( step % thermal_subcycle_steps == 0 ) - computeHeatTransfer( *heat_transfer, particles, - neigh_iter_tag{}, + computeHeatTransfer( *heat_transfer, particles, *neighbor, thermal_subcycle_steps * dt ); } @@ -343,7 +342,7 @@ class Solver updateForce(); if constexpr ( is_contact::value ) - computeForce( *contact, particles, neigh_iter_tag{}, false ); + computeForce( *contact, particles, *neighbor, false ); // Add force boundary condition. if ( boundary_condition.forceUpdate() ) @@ -378,7 +377,7 @@ class Solver updateForce(); if constexpr ( is_contact::value ) - computeForce( *contact, particles, neigh_iter_tag{}, false ); + computeForce( *contact, particles, *neighbor, false ); if constexpr ( is_temperature_dependent< typename force_model_type::thermal_type>::value ) @@ -403,15 +402,15 @@ class Solver if constexpr ( is_fracture< typename force_model_type::fracture_type>::value ) { - force->computeWeightedVolume( particles, neigh_iter_tag{} ); + force->computeWeightedVolume( particles, *neighbor ); comm->gatherWeightedVolume(); } // Compute and communicate dilatation for LPS (does nothing for PMB). - force->computeDilatation( particles, neigh_iter_tag{} ); + force->computeDilatation( particles, *neighbor ); comm->gatherDilatation(); // Compute internal forces. - computeForce( *force, particles, neigh_iter_tag{} ); + computeForce( *force, particles, *neighbor ); } void output( const int step ) @@ -419,7 +418,7 @@ class Solver // Print output. if ( step % output_frequency == 0 ) { - auto W = computeEnergy( *force, particles, neigh_iter_tag() ); + auto W = computeEnergy( *force, particles, *neighbor ); particles.output( step / output_frequency, step * dt, output_reference ); @@ -436,11 +435,11 @@ class Solver { // Output after construction and initial forces. std::ofstream out( output_file, std::ofstream::app ); - _init_time += _init_timer.time() + _neighbor_timer.time() + + _init_time += _init_timer.time() + neighbor->time() + 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, "Init-Neighbor-Time(s): ", neighbor->time(), "\n" ); log( out, "#Timestep/Total-steps Simulation-time Total-strain-energy " "Step-Time(s) Force-Time(s) Comm-Time(s) Integrate-Time(s) " "Energy-Time(s) Output-Time(s) Particle*steps/s" ); @@ -459,7 +458,7 @@ class Solver double integrate_time = integrator->time(); double force_time = force->time(); double energy_time = force->timeEnergy(); - double neigh_time = force->timeNeighbor(); + double neigh_time = neighbor->time(); double output_time = particles.timeOutput(); _total_time += step_time; // Instantaneous rate. @@ -487,7 +486,7 @@ class Solver double force_time = force->time(); double energy_time = force->timeEnergy(); double output_time = particles.timeOutput(); - double neighbor_time = _neighbor_timer.time(); + double neighbor_time = neighbor->time(); _total_time = _init_time + comm_time + integrate_time + force_time + energy_time + output_time + particles.time(); @@ -532,7 +531,7 @@ class Solver "Cannot create prenotch in system without fracture." ); // Create prenotch. - force->prenotch( exec_space{}, particles, prenotch ); + neighbor->prenotch( exec_space{}, particles, prenotch ); _init_time += prenotch.time(); } @@ -540,6 +539,7 @@ class Solver Inputs inputs; std::shared_ptr comm; std::shared_ptr integrator; + std::shared_ptr neighbor; std::shared_ptr force; // Optional modules. std::shared_ptr heat_transfer; @@ -552,7 +552,6 @@ class Solver // Note: init_time is combined from many class timers. double _init_time; Timer _init_timer; - Timer _neighbor_timer; Timer _step_timer; double _total_time; bool print; diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index 1e265b3f..18a776cb 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -35,9 +35,15 @@ struct is_fracture : public std::true_type // Mechanics tags. struct Elastic { + using base_type = Elastic; +}; +struct Plastic +{ + using base_type = Plastic; }; struct ElasticPerfectlyPlastic { + using base_type = Plastic; }; // Model category tags. diff --git a/src/force/CabanaPD_Contact.hpp b/src/force/CabanaPD_Contact.hpp index b2ab8df5..f6bd49d5 100644 --- a/src/force/CabanaPD_Contact.hpp +++ b/src/force/CabanaPD_Contact.hpp @@ -38,98 +38,29 @@ KOKKOS_INLINE_FUNCTION void getRelativeNormalVelocityComponents( vn /= r; }; -// Contact forces base class. -template -class BaseForceContact : public Force -{ - public: - using base_type = Force; - using neighbor_list_type = typename base_type::neighbor_list_type; - - // NOTE: using 2x radius to find neighbors when particles first touch. - template - BaseForceContact( const bool half_neigh, const ParticleType& particles, - const ModelType model ) - : base_type( half_neigh, 2.0 * model.radius + model.radius_extend, - particles.sliceCurrentPosition(), particles.frozenOffset(), - particles.localOffset(), particles.ghost_mesh_lo, - particles.ghost_mesh_hi ) - , search_radius( 2.0 * model.radius + model.radius_extend ) - , radius_extend( model.radius_extend ) - { - for ( int d = 0; d < particles.dim; d++ ) - { - mesh_min[d] = particles.ghost_mesh_lo[d]; - mesh_max[d] = particles.ghost_mesh_hi[d]; - } - } - - // Only rebuild neighbor list as needed. - template - void update( const ParticleType& particles, - const bool require_update = false ) - { - double max_displacement = particles.getMaxDisplacement(); - if ( max_displacement > radius_extend || require_update ) - { - _neigh_timer.start(); - const auto y = particles.sliceCurrentPosition(); - _neigh_list.build( y, particles.frozenOffset(), - particles.localOffset(), search_radius, 1.0, - mesh_min, mesh_max ); - // Reset neighbor update displacement. - const auto u = particles.sliceDisplacement(); - auto u_neigh = particles.sliceDisplacementNeighborBuild(); - Cabana::deep_copy( u_neigh, u ); - _neigh_timer.stop(); - } - } - - auto timeNeighbor() { return _neigh_timer.time(); }; - - protected: - double search_radius; - double radius_extend; - Timer _neigh_timer; - - using base_type::_half_neigh; - using base_type::_neigh_list; - double mesh_max[3]; - double mesh_min[3]; -}; - /****************************************************************************** Normal repulsion forces ******************************************************************************/ template -class Force - : public BaseForceContact +class Force : public BaseForce { public: - using base_type = BaseForceContact; - using neighbor_list_type = typename base_type::neighbor_list_type; - - template - Force( const bool half_neigh, const ParticleType& particles, - const NormalRepulsionModel model ) - : base_type( half_neigh, particles, model ) - , _model( model ) + using base_type = BaseForce; + + Force( const NormalRepulsionModel model ) + : _model( model ) { } template + class NeighborType> void computeForceFull( ForceType& fc, const PosType& x, const PosType& u, - const ParticleType& particles, - ParallelType& neigh_op_tag ) + ParticleType& particles, NeighborType& neighbor ) { auto model = _model; const auto vol = particles.sliceVolume(); - const auto y = particles.sliceCurrentPosition(); - const int n_frozen = particles.frozenOffset(); - const int n_local = particles.localOffset(); - base_type::update( particles, particles.getMaxDisplacement() ); + neighbor.update( particles, model.cutoff(), model.extend() ); auto contact_full = KOKKOS_LAMBDA( const int i, const int j ) { @@ -157,32 +88,24 @@ class Force // FIXME: using default space for now. using exec_space = typename MemorySpace::execution_space; - Kokkos::RangePolicy policy( n_frozen, n_local ); - Cabana::neighbor_parallel_for( - policy, contact_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::Contact::compute_full" ); + neighbor.iterate( exec_space{}, contact_full, particles, + "CabanaPD::Contact::compute_full" ); _timer.stop(); } // FIXME: implement energy template + class NeighborType> double computeEnergyFull( WType&, const PosType&, const PosType&, - ParticleType&, const int, ParallelType& ) + ParticleType&, const int, NeighborType& ) { return 0.0; } protected: NormalRepulsionModel _model; - using base_type::_half_neigh; - using base_type::_neigh_list; - using base_type::_neigh_timer; using base_type::_timer; - - double mesh_max[3]; - double mesh_min[3]; }; } // namespace CabanaPD diff --git a/src/force/CabanaPD_Hertzian.hpp b/src/force/CabanaPD_Hertzian.hpp index 487cd686..e03a4aa6 100644 --- a/src/force/CabanaPD_Hertzian.hpp +++ b/src/force/CabanaPD_Hertzian.hpp @@ -25,35 +25,27 @@ namespace CabanaPD Normal repulsion forces ******************************************************************************/ template -class Force : public BaseForceContact +class Force : public BaseForce { public: - using base_type = BaseForceContact; - using neighbor_list_type = typename base_type::neighbor_list_type; - - template - Force( const bool half_neigh, const ParticleType& particles, - const HertzianModel model ) - : base_type( half_neigh, particles, model ) - , _model( model ) + using base_type = BaseForce; + + Force( const HertzianModel model ) + : _model( model ) { } template + class NeighborType> void computeForceFull( ForceType& fc, const PosType& x, const PosType& u, - const ParticleType& particles, - ParallelType& neigh_op_tag ) + ParticleType& particles, NeighborType& neighbor ) { - const int n_frozen = particles.frozenOffset(); - const int n_local = particles.localOffset(); - auto model = _model; const auto vol = particles.sliceVolume(); const auto rho = particles.sliceDensity(); const auto vel = particles.sliceVelocity(); - base_type::update( particles, particles.getMaxDisplacement() ); + neighbor.update( particles, model.cutoff(), model.extend() ); auto contact_full = KOKKOS_LAMBDA( const int i, const int j ) { @@ -76,11 +68,8 @@ class Force : public BaseForceContact // FIXME: using default space for now. using exec_space = typename MemorySpace::execution_space; - Kokkos::RangePolicy policy( n_frozen, n_local ); - Cabana::neighbor_parallel_for( - policy, contact_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::Contact::compute_full" ); - + neighbor.iterate( exec_space{}, contact_full, particles, + "CabanaPD::Hertzian::compute" ); _timer.stop(); } @@ -95,9 +84,6 @@ class Force : public BaseForceContact protected: HertzianModel _model; - using base_type::_half_neigh; - using base_type::_neigh_list; - using base_type::_neigh_timer; using base_type::_timer; }; diff --git a/src/force/CabanaPD_LPS.hpp b/src/force/CabanaPD_LPS.hpp index 6a8d46d9..f5d189fa 100644 --- a/src/force/CabanaPD_LPS.hpp +++ b/src/force/CabanaPD_LPS.hpp @@ -71,11 +71,10 @@ namespace CabanaPD { template class Force> - : public Force + : public BaseForce { protected: - using base_type = Force; - using base_type::_half_neigh; + using base_type = BaseForce; using model_type = ForceModel; model_type _model; @@ -85,20 +84,15 @@ class Force> public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; - using neighbor_list_type = typename base_type::neighbor_list_type; - using base_type::_neigh_list; - template - Force( const bool half_neigh, ParticleType& particles, - const model_type model ) - : base_type( half_neigh, model.delta, particles ) - , _model( model ) + Force( const model_type model ) + : _model( model ) { } - template + template void computeWeightedVolume( ParticleType& particles, - const ParallelType neigh_op_tag ) + const NeighborType& neighbor ) { _timer.start(); @@ -116,19 +110,14 @@ class Force> getDistance( x, u, i, j, xi, r, s ); m( i ) += model.weightedVolume( xi, vol( j ) ); }; - - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, weighted_volume, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLPS::computeWeightedVolume" ); - + neighbor.iterate( exec_space{}, weighted_volume, particles, + "CabanaPD::Dilatation::computeWeightedVolume" ); _timer.stop(); } - template + template void computeDilatation( ParticleType& particles, - const ParallelType neigh_op_tag ) + const NeighborType neighbor ) { _timer.start(); @@ -148,20 +137,17 @@ class Force> theta( i ) += model.dilatation( s, xi, vol( j ), m( i ) ); }; - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, dilatation, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLPS::computeDilatation" ); + neighbor.iterate( exec_space{}, dilatation, particles, + "CabanaPD::Dilatation::compute" ); _timer.stop(); } template + class NeighborType> void computeForceFull( ForceType& f, const PosType& x, const PosType& u, const ParticleType& particles, - ParallelType& neigh_op_tag ) + NeighborType& neighbor ) { _timer.start(); @@ -190,29 +176,26 @@ class Force> f( i, 2 ) += fz_i; }; - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, force_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLPS::computeFull" ); - + neighbor.iterate( exec_space{}, force_full, particles, + "CabanaPD::ForceLPS::computeFull" ); _timer.stop(); } template + class NeighborType> double computeEnergyFull( WType& W, const PosType& x, const PosType& u, const ParticleType& particles, - ParallelType& neigh_op_tag ) + NeighborType& neighbor ) { _energy_timer.start(); auto model = _model; - auto neigh_list = _neigh_list; const auto vol = particles.sliceVolume(); const auto theta = particles.sliceDilatation(); auto m = particles.sliceWeightedVolume(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); auto energy_full = KOKKOS_LAMBDA( const int i, const int j, double& Phi ) @@ -230,14 +213,9 @@ class Force> Phi += w * vol( i ); }; - double strain_energy = 0.0; - - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_reduce( - policy, energy_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, strain_energy, - "CabanaPD::ForceLPS::computeEnergyFull" ); + double strain_energy = + neighbor.reduce( exec_space{}, energy_full, particles, + "CabanaPD::ForceLPS::computeEnergyFull" ); _energy_timer.stop(); return strain_energy; @@ -249,15 +227,10 @@ class Force> template class Force> - : public Force>, - public BaseFracture + : public Force> { protected: - using fracture_type = BaseFracture; - using fracture_type::_mu; - using base_type = Force>; - using base_type::_half_neigh; using model_type = ForceModel; model_type _model; @@ -267,39 +240,29 @@ class Force> public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; - using neighbor_list_type = typename base_type::neighbor_list_type; - using base_type::_neigh_list; - - template - Force( const bool half_neigh, const ParticleType& particles, - const model_type model ) - : base_type( half_neigh, particles, model ) - , fracture_type( particles.localOffset(), - base_type::getMaxLocalNeighbors() ) - , _model( model ) - { - } - template - void prenotch( ExecSpace exec_space, const ParticleType& particles, - PrenotchType& prenotch ) + Force( const model_type model ) + : base_type( model ) + , _model( model ) { - fracture_type::prenotch( exec_space, particles, prenotch, _neigh_list ); } - template - void computeWeightedVolume( ParticleType& particles, ParallelType ) + template + void computeWeightedVolume( ParticleType& particles, + const NeighborType& neighbor ) { _timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); + auto x = particles.sliceReferencePosition(); auto u = particles.sliceDisplacement(); const auto vol = particles.sliceVolume(); auto m = particles.sliceWeightedVolume(); Cabana::deep_copy( m, 0.0 ); auto model = _model; - auto neigh_list = _neigh_list; - auto mu = _mu; auto weighted_volume = KOKKOS_LAMBDA( const int i ) { @@ -320,18 +283,21 @@ class Force> } }; - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Kokkos::parallel_for( "CabanaPD::ForceLPSDamage::computeWeightedVolume", - policy, weighted_volume ); + neighbor.iterateLinear( + exec_space{}, weighted_volume, particles, + "CabanaPD::ForceLPSDamage::computeWeightedVolume" ); _timer.stop(); } - template - void computeDilatation( ParticleType& particles, ParallelType ) + template + void computeDilatation( ParticleType& particles, + const NeighborType& neighbor ) { _timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); const auto x = particles.sliceReferencePosition(); auto u = particles.sliceDisplacement(); @@ -339,8 +305,7 @@ class Force> auto m = particles.sliceWeightedVolume(); auto theta = particles.sliceDilatation(); auto model = _model; - auto neigh_list = _neigh_list; - auto mu = _mu; + Cabana::deep_copy( theta, 0.0 ); auto dilatation = KOKKOS_LAMBDA( const int i ) @@ -368,25 +333,25 @@ class Force> } }; - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Kokkos::parallel_for( "CabanaPD::ForceLPSDamage::computeDilatation", - policy, dilatation ); + neighbor.iterateLinear( exec_space{}, dilatation, particles, + "CabanaPD::ForceLPSDamage::computeDilatation" ); _timer.stop(); } template + class NeighborType> void computeForceFull( ForceType& f, const PosType& x, const PosType& u, - const ParticleType& particles, ParallelType ) + const ParticleType& particles, + NeighborType& neighbor ) { _timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); auto break_coeff = _model.bond_break_coeff; auto model = _model; - auto neigh_list = _neigh_list; - auto mu = _mu; const auto vol = particles.sliceVolume(); auto theta = particles.sliceDilatation(); @@ -438,24 +403,24 @@ class Force> } }; - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Kokkos::parallel_for( "CabanaPD::ForceLPSDamage::computeFull", policy, - force_full ); + neighbor.iterateLinear( exec_space{}, force_full, particles, + "CabanaPD::ForceLPSDamage::computeFull" ); _timer.stop(); } template + class NeighborType> double computeEnergyFull( WType& W, const PosType& x, const PosType& u, - ParticleType& particles, ParallelType& ) + ParticleType& particles, + const NeighborType& neighbor ) { _energy_timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); auto model = _model; - auto neigh_list = _neigh_list; - auto mu = _mu; const auto vol = particles.sliceVolume(); const auto theta = particles.sliceDilatation(); @@ -493,11 +458,9 @@ class Force> phi( i ) = 1 - phi_i / vol_H_i; }; - double strain_energy = 0.0; - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Kokkos::parallel_reduce( "CabanaPD::ForceLPSDamage::computeEnergyFull", - policy, energy_full, strain_energy ); + double strain_energy = neighbor.reduceLinear( + exec_space{}, energy_full, particles, + "CabanaPD::ForceLPSDamage::computeEnergyFull" ); _energy_timer.stop(); return strain_energy; @@ -519,22 +482,18 @@ class Force> public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; - using neighbor_list_type = typename base_type::neighbor_list_type; - using base_type::_neigh_list; - template - Force( const bool half_neigh, ParticleType& particles, - const model_type model ) - : base_type( half_neigh, particles, model ) + Force( const model_type model ) + : base_type( model ) , _model( model ) { } template + class NeighborType> void computeForceFull( ForceType& f, const PosType& x, const PosType& u, const ParticleType& particles, - ParallelType& neigh_op_tag ) + NeighborType& neighbor ) { _timer.start(); @@ -568,25 +527,22 @@ class Force> f( i, 2 ) += fz_i; }; - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, force_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLPS::computeFull" ); - + neighbor.iterate( exec_space{}, force_full, particles, + "CabanaPD::ForceLPS::computeFull" ); _timer.stop(); } template + class NeighborType> double computeEnergyFull( WType& W, const PosType& x, const PosType& u, const ParticleType& particles, - ParallelType& neigh_op_tag ) + NeighborType& neighbor ) { _energy_timer.start(); auto model = _model; - auto neigh_list = _neigh_list; + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); const auto vol = particles.sliceVolume(); const auto theta = particles.sliceDilatation(); @@ -609,14 +565,9 @@ class Force> Phi += w * vol( i ); }; - double strain_energy = 0.0; - - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_reduce( - policy, energy_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, strain_energy, - "CabanaPD::ForceLPS::computeEnergyFull" ); + double strain_energy = + neighbor.reduce( exec_space{}, energy_full, particles, + "CabanaPD::ForceLPS::computeEnergyFull" ); _energy_timer.stop(); return strain_energy; diff --git a/src/force/CabanaPD_PMB.hpp b/src/force/CabanaPD_PMB.hpp index 65b9f5c5..9559591b 100644 --- a/src/force/CabanaPD_PMB.hpp +++ b/src/force/CabanaPD_PMB.hpp @@ -73,38 +73,32 @@ namespace CabanaPD template class Force> - : public Force + : public BaseForce { public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; using model_type = ForceModel; - using base_type = Force; - using neighbor_list_type = typename base_type::neighbor_list_type; - using base_type::_neigh_list; + using base_type = BaseForce; protected: - using base_type::_half_neigh; model_type _model; using base_type::_energy_timer; using base_type::_timer; public: - template - Force( const bool half_neigh, const ParticleType& particles, - const model_type model ) - : base_type( half_neigh, model.delta, particles ) - , _model( model ) + Force( const model_type model ) + : _model( model ) { } template + class NeighborType> void computeForceFull( ForceType& f, const PosType& x, const PosType& u, const ParticleType& particles, - ParallelType& neigh_op_tag ) + NeighborType& neighbor ) { _timer.start(); @@ -133,20 +127,16 @@ class Force policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, force_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForcePMB::computeFull" ); - + neighbor.iterate( exec_space{}, force_full, particles, + "CabanaPD::ForcePMB::computeFull" ); _timer.stop(); } template + class NeighborType> double computeEnergyFull( WType& W, const PosType& x, const PosType& u, const ParticleType& particles, - ParallelType& neigh_op_tag ) + NeighborType& neighbor ) { _energy_timer.start(); @@ -167,14 +157,9 @@ class Force policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_reduce( - policy, energy_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, strain_energy, - "CabanaPD::ForcePMB::computeEnergyFull" ); - + double strain_energy = + neighbor.reduce( exec_space{}, energy_full, particles, + "CabanaPD::ForcePMB::computeEnergyFull" ); _energy_timer.stop(); return strain_energy; } @@ -183,59 +168,39 @@ class Force class Force> - : public Force, - public BaseFracture + : public BaseForce { public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; using model_type = ForceModel; - using base_type = Force; - using neighbor_list_type = typename base_type::neighbor_list_type; - using base_type::_neigh_list; + using base_type = BaseForce; protected: - using fracture_type = BaseFracture; - using fracture_type::_mu; - - using base_model_type = typename model_type::base_type; - using base_type::_half_neigh; model_type _model; using base_type::_energy_timer; using base_type::_timer; public: - template - Force( const bool half_neigh, const ParticleType& particles, - const model_type model ) - : base_type( half_neigh, model.delta, particles ) - , fracture_type( particles.localOffset(), - base_type::getMaxLocalNeighbors() ) + Force( const model_type model ) + : base_type() , _model( model ) { - // Needed only for models which store per-bond information. - _model.updateBonds( particles.localOffset(), - base_type::getMaxLocalNeighbors() ); - } - - template - void prenotch( ExecSpace exec_space, const ParticleType& particles, - PrenotchType& prenotch ) - { - fracture_type::prenotch( exec_space, particles, prenotch, _neigh_list ); } template + class NeighborType> void computeForceFull( ForceType& f, const PosType& x, const PosType& u, - const ParticleType& particles, ParallelType& ) + const ParticleType& particles, + const NeighborType& neighbor ) { _timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); auto model = _model; - auto neigh_list = _neigh_list; - auto mu = _mu; const auto vol = particles.sliceVolume(); const auto nofail = particles.sliceNoFail(); @@ -284,24 +249,22 @@ class Force policy( particles.frozenOffset(), - particles.localOffset() ); - Kokkos::parallel_for( "CabanaPD::ForcePMBDamage::computeFull", policy, - force_full ); - + neighbor.iterateLinear( exec_space{}, force_full, particles, + "CabanaPD::ForcePMBDamage::computeFull" ); _timer.stop(); } template + class NeighborType> double computeEnergyFull( WType& W, const PosType& x, const PosType& u, - ParticleType& particles, ParallelType& ) + ParticleType& particles, NeighborType& neighbor ) { _energy_timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); auto model = _model; - auto neigh_list = _neigh_list; - auto mu = _mu; const auto vol = particles.sliceVolume(); auto phi = particles.sliceDamage(); @@ -333,11 +296,9 @@ class Force policy( particles.frozenOffset(), - particles.localOffset() ); - Kokkos::parallel_reduce( "CabanaPD::ForcePMBDamage::computeEnergyFull", - policy, energy_full, strain_energy ); + double strain_energy = neighbor.reduceLinear( + exec_space{}, energy_full, particles, + "CabanaPD::ForcePMBDamage::computeEnergyFull" ); _energy_timer.stop(); return strain_energy; @@ -347,37 +308,31 @@ class Force class Force> - : public Force + : public BaseForce { public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; using model_type = ForceModel; - using base_type = Force; - using neighbor_list_type = typename base_type::neighbor_list_type; - using base_type::_neigh_list; + using base_type = BaseForce; protected: - using base_type::_half_neigh; model_type _model; using base_type::_energy_timer; using base_type::_timer; public: - template - Force( const bool half_neigh, const ParticleType& particles, - const model_type model ) - : base_type( half_neigh, model.delta, particles ) - , _model( model ) + Force( const model_type model ) + : _model( model ) { } template + class NeighborType> void computeForceFull( ForceType& f, const PosType& x, const PosType& u, - ParticleType& particles, ParallelType& neigh_op_tag ) + ParticleType& particles, NeighborType& neighbor ) { _timer.start(); @@ -408,20 +363,15 @@ class Force policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, force_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLinearPMB::computeFull" ); - + neighbor.iterate( exec_space{}, force_full, particles, + "CabanaPD::ForceLinearPMB::computeFull" ); _timer.stop(); } template + class NeighborType> double computeEnergyFull( WType& W, const PosType& x, const PosType& u, - ParticleType& particles, - ParallelType& neigh_op_tag ) + ParticleType& particles, NeighborType& neighbor ) { _energy_timer.start(); @@ -442,14 +392,9 @@ class Force policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_reduce( - policy, energy_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, strain_energy, - "CabanaPD::ForceLinearPMB::computeEnergyFull" ); - + double strain_energy = + neighbor.reduce( exec_space{}, energy_full, particles, + "CabanaPD::ForceLinearPMB::computeEnergyFull" ); _energy_timer.stop(); return strain_energy; } diff --git a/src/force_models/CabanaPD_Contact.hpp b/src/force_models/CabanaPD_Contact.hpp index be56bfff..8e5d6d2a 100644 --- a/src/force_models/CabanaPD_Contact.hpp +++ b/src/force_models/CabanaPD_Contact.hpp @@ -39,6 +39,11 @@ struct ContactModel ContactModel( const double _radius, const double _radius_extend ) : radius( _radius ) , radius_extend( _radius_extend ){}; + + auto cutoff() const { return 2.0 * radius + radius_extend; } + auto extend() const { return radius_extend; } + + void updateBonds( const int, const int ) {} }; /* Normal repulsion */ diff --git a/unit_test/tstForce.hpp b/unit_test/tstForce.hpp index 416361bf..fc053966 100644 --- a/unit_test/tstForce.hpp +++ b/unit_test/tstForce.hpp @@ -648,20 +648,21 @@ void checkAnalyticalDilatation( ModelType, QuadraticTag, const double, { } -template +template double computeEnergyAndForce( ForceType force, ParticleType& particles, - const int ) + const NeighborType& neighbor ) { - computeForce( force, particles, Cabana::SerialOpTag() ); - double Phi = computeEnergy( force, particles, Cabana::SerialOpTag() ); + computeForce( force, particles, neighbor ); + double Phi = computeEnergy( force, particles, neighbor ); return Phi; } -template -void initializeForce( ForceType& force, ParticleType& particles ) +template +void initializeForce( ForceType& force, ParticleType& particles, + const NeighborType& neighbor ) { - force.computeWeightedVolume( particles, Cabana::SerialOpTag{} ); - force.computeDilatation( particles, Cabana::SerialOpTag{} ); + force.computeWeightedVolume( particles, neighbor ); + force.computeDilatation( particles, neighbor ); } template @@ -692,20 +693,18 @@ void testForce( ModelType model, const double dx, const double m, // This needs to exactly match the mesh spacing to compare with the single // particle calculation. - CabanaPD::Force force( true, particles, model ); - + CabanaPD::Force force( model ); + CabanaPD::Neighbor + neighbor( false, model, particles ); auto x = particles.sliceReferencePosition(); auto f = particles.sliceForce(); auto W = particles.sliceStrainEnergy(); auto vol = particles.sliceVolume(); // No communication needed (as in the main solver) since this test is only // intended for one rank. - initializeForce( force, particles ); + initializeForce( force, particles, neighbor ); - unsigned int max_neighbors; - unsigned long long total_neighbors; - force.getNeighborStatistics( max_neighbors, total_neighbors ); - double Phi = computeEnergyAndForce( force, particles, max_neighbors ); + double Phi = computeEnergyAndForce( force, particles, neighbor ); // Make a copy of final results on the host std::size_t num_particle = x.size(); @@ -734,9 +733,41 @@ void testForce( ModelType model, const double dx, const double m, boundary_width, Phi ); } +//---------------------------------------------------------------------------// +// Neighbor test function. +//---------------------------------------------------------------------------// +void testNeighbor( const double dx, const double m, const double delta ) +{ + using model_type = CabanaPD::PMB; + CabanaPD::ForceModel model( model_type{}, CabanaPD::Elastic{}, + CabanaPD::NoFracture{}, delta, 1.0 ); + + auto particles = createParticles( model, LinearTag{}, dx, 0.0 ); + + CabanaPD::Neighbor neighbor( + false, model, particles ); + + auto expected_num_neighbors = computeReferenceNeighbors( model.delta, m ); + EXPECT_EQ( expected_num_neighbors, neighbor.getMaxLocal() ); +} + //---------------------------------------------------------------------------// // GTest tests. //---------------------------------------------------------------------------// +TEST( TEST_CATEGORY, test_neighbor ) +{ + // dx needs to be decreased for increased m: boundary particles are ignored. + double m = 3; + double dx = 2.0 / 11.0; + double delta = dx * m; + testNeighbor( dx, m, delta ); + + m = 6; + dx = 2.0 / 15.0; + delta = dx * m; + testNeighbor( dx, m, delta ); +} + TEST( TEST_CATEGORY, test_force_pmb ) { // dx needs to be decreased for increased m: boundary particles are ignored. From 496a323611fe8e2e47b9be9c9397b672056a8a51 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 31 Jan 2025 13:00:53 -0500 Subject: [PATCH 02/30] Separate dilatation from LPS forces --- src/CabanaPD_Force.hpp | 198 +++++++++++++++++++++++++++++++++++++ src/force/CabanaPD_LPS.hpp | 184 +++++----------------------------- 2 files changed, 220 insertions(+), 162 deletions(-) diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index 458b7334..b0aac11a 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -152,6 +152,204 @@ class BaseForce auto timeEnergy() { return _energy_timer.time(); }; }; +/****************************************************************************** + Dilatation. +******************************************************************************/ +template +class Dilatation; + +template class ModelType, + class PDModelType, class MechanicsType, class... ModelParams> +class Dilatation> +{ + protected: + using model_type = + ModelType; + model_type _model; + + Timer _timer; + + public: + // Using the default exec_space. + using exec_space = typename MemorySpace::execution_space; + + Dilatation( const model_type model ) + : _model( model ) + { + } + + template + void computeWeightedVolume( ParticleType& particles, + const NeighborType& neighbor ) + { + _timer.start(); + + auto x = particles.sliceReferencePosition(); + auto u = particles.sliceDisplacement(); + const auto vol = particles.sliceVolume(); + auto m = particles.sliceWeightedVolume(); + Cabana::deep_copy( m, 0.0 ); + auto model = _model; + + auto weighted_volume = KOKKOS_LAMBDA( const int i, const int j ) + { + // Get the reference positions and displacements. + double xi, r, s; + getDistance( x, u, i, j, xi, r, s ); + m( i ) += model.weightedVolume( xi, vol( j ) ); + }; + neighbor.iterate( exec_space{}, weighted_volume, particles, + "CabanaPD::Dilatation::computeWeightedVolume" ); + _timer.stop(); + } + + template + void computeDilatation( ParticleType& particles, + const NeighborType neighbor ) + { + _timer.start(); + + const auto x = particles.sliceReferencePosition(); + auto u = particles.sliceDisplacement(); + const auto vol = particles.sliceVolume(); + auto m = particles.sliceWeightedVolume(); + auto theta = particles.sliceDilatation(); + auto model = _model; + Cabana::deep_copy( theta, 0.0 ); + + auto dilatation = KOKKOS_LAMBDA( const int i, const int j ) + { + // Get the bond distance, displacement, and stretch. + double xi, r, s; + getDistance( x, u, i, j, xi, r, s ); + theta( i ) += model.dilatation( s, xi, vol( j ), m( i ) ); + }; + + neighbor.iterate( exec_space{}, dilatation, particles, + "CabanaPD::Dilatation::compute" ); + + _timer.stop(); + } + + auto time() { return _timer.time(); }; +}; + +template class ModelType, + class PDModelType, class MechanicsType, class... ModelParams> +class Dilatation> +{ + protected: + using model_type = + ModelType; + model_type _model; + + Timer _timer; + + public: + // Using the default exec_space. + using exec_space = typename MemorySpace::execution_space; + + Dilatation( const model_type model ) + : _model( model ) + { + } + + template + void computeWeightedVolume( ParticleType& particles, + const NeighborType& neighbor ) + { + _timer.start(); + + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); + + auto x = particles.sliceReferencePosition(); + auto u = particles.sliceDisplacement(); + const auto vol = particles.sliceVolume(); + auto m = particles.sliceWeightedVolume(); + Cabana::deep_copy( m, 0.0 ); + auto model = _model; + + auto weighted_volume = KOKKOS_LAMBDA( const int i ) + { + std::size_t num_neighbors = + Cabana::NeighborList::numNeighbor( + neigh_list, i ); + for ( std::size_t n = 0; n < num_neighbors; n++ ) + { + std::size_t j = + Cabana::NeighborList::getNeighbor( + neigh_list, i, n ); + + // Get the reference positions and displacements. + double xi, r, s; + getDistance( x, u, i, j, xi, r, s ); + // mu is included to account for bond breaking. + m( i ) += mu( i, n ) * model.weightedVolume( xi, vol( j ) ); + } + }; + + neighbor.iterateLinear( exec_space{}, weighted_volume, particles, + "CabanaPD::Dilatation::computeWeightedVolume" ); + + _timer.stop(); + } + + template + void computeDilatation( ParticleType& particles, + const NeighborType& neighbor ) + { + _timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); + + const auto x = particles.sliceReferencePosition(); + auto u = particles.sliceDisplacement(); + const auto vol = particles.sliceVolume(); + auto m = particles.sliceWeightedVolume(); + auto theta = particles.sliceDilatation(); + auto model = _model; + + Cabana::deep_copy( theta, 0.0 ); + + auto dilatation = KOKKOS_LAMBDA( const int i ) + { + std::size_t num_neighbors = + Cabana::NeighborList::numNeighbor( + neigh_list, i ); + for ( std::size_t n = 0; n < num_neighbors; n++ ) + { + std::size_t j = + Cabana::NeighborList::getNeighbor( + neigh_list, i, n ); + + // Get the bond distance, displacement, and stretch. + double xi, r, s; + getDistance( x, u, i, j, xi, r, s ); + + // Check if all bonds are broken (m=0) to avoid dividing by + // zero. Alternatively, one could check if this bond mu(i,n) is + // broken, because m=0 only occurs when all bonds are broken. + // mu is still included to account for individual bond breaking. + if ( m( i ) > 0 ) + theta( i ) += mu( i, n ) * + model.dilatation( s, xi, vol( j ), m( i ) ); + } + }; + + neighbor.iterateLinear( exec_space{}, dilatation, particles, + "CabanaPD::Dilatation::compute" ); + + _timer.stop(); + } +}; + /****************************************************************************** Force free functions. ******************************************************************************/ diff --git a/src/force/CabanaPD_LPS.hpp b/src/force/CabanaPD_LPS.hpp index f5d189fa..cf4a4586 100644 --- a/src/force/CabanaPD_LPS.hpp +++ b/src/force/CabanaPD_LPS.hpp @@ -71,76 +71,24 @@ namespace CabanaPD { template class Force> - : public BaseForce + : public Dilatation> { protected: - using base_type = BaseForce; - using model_type = ForceModel; - model_type _model; + using base_type = + Dilatation>; + using model_type = typename base_type::model_type; + using base_type::_model; - using base_type::_energy_timer; using base_type::_timer; + Timer _energy_timer; public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; Force( const model_type model ) - : _model( model ) - { - } - - template - void computeWeightedVolume( ParticleType& particles, - const NeighborType& neighbor ) - { - _timer.start(); - - auto x = particles.sliceReferencePosition(); - auto u = particles.sliceDisplacement(); - const auto vol = particles.sliceVolume(); - auto m = particles.sliceWeightedVolume(); - Cabana::deep_copy( m, 0.0 ); - auto model = _model; - - auto weighted_volume = KOKKOS_LAMBDA( const int i, const int j ) - { - // Get the reference positions and displacements. - double xi, r, s; - getDistance( x, u, i, j, xi, r, s ); - m( i ) += model.weightedVolume( xi, vol( j ) ); - }; - neighbor.iterate( exec_space{}, weighted_volume, particles, - "CabanaPD::Dilatation::computeWeightedVolume" ); - _timer.stop(); - } - - template - void computeDilatation( ParticleType& particles, - const NeighborType neighbor ) + : base_type( model ) { - _timer.start(); - - const auto x = particles.sliceReferencePosition(); - auto u = particles.sliceDisplacement(); - const auto vol = particles.sliceVolume(); - auto m = particles.sliceWeightedVolume(); - auto theta = particles.sliceDilatation(); - auto model = _model; - Cabana::deep_copy( theta, 0.0 ); - - auto dilatation = KOKKOS_LAMBDA( const int i, const int j ) - { - // Get the bond distance, displacement, and stretch. - double xi, r, s; - getDistance( x, u, i, j, xi, r, s ); - theta( i ) += model.dilatation( s, xi, vol( j ), m( i ) ); - }; - - neighbor.iterate( exec_space{}, dilatation, particles, - "CabanaPD::Dilatation::compute" ); - - _timer.stop(); } template > return strain_energy; } - auto time() { return _timer.time(); }; auto timeEnergy() { return _energy_timer.time(); }; }; template class Force> - : public Force> + : public Dilatation> { protected: - using base_type = Force>; - using model_type = ForceModel; - model_type _model; + using base_type = + Dilatation>; + using model_type = typename base_type::model_type; + using base_type::_model; - using base_type::_energy_timer; using base_type::_timer; + Timer _energy_timer; public: // Using the default exec_space. @@ -243,102 +191,9 @@ class Force> Force( const model_type model ) : base_type( model ) - , _model( model ) { } - template - void computeWeightedVolume( ParticleType& particles, - const NeighborType& neighbor ) - { - _timer.start(); - - using neighbor_list_type = typename NeighborType::list_type; - auto neigh_list = neighbor.list(); - const auto mu = neighbor.brokenBonds(); - - auto x = particles.sliceReferencePosition(); - auto u = particles.sliceDisplacement(); - const auto vol = particles.sliceVolume(); - auto m = particles.sliceWeightedVolume(); - Cabana::deep_copy( m, 0.0 ); - auto model = _model; - - auto weighted_volume = KOKKOS_LAMBDA( const int i ) - { - std::size_t num_neighbors = - Cabana::NeighborList::numNeighbor( - neigh_list, i ); - for ( std::size_t n = 0; n < num_neighbors; n++ ) - { - std::size_t j = - Cabana::NeighborList::getNeighbor( - neigh_list, i, n ); - - // Get the reference positions and displacements. - double xi, r, s; - getDistance( x, u, i, j, xi, r, s ); - // mu is included to account for bond breaking. - m( i ) += mu( i, n ) * model.weightedVolume( xi, vol( j ) ); - } - }; - - neighbor.iterateLinear( - exec_space{}, weighted_volume, particles, - "CabanaPD::ForceLPSDamage::computeWeightedVolume" ); - - _timer.stop(); - } - - template - void computeDilatation( ParticleType& particles, - const NeighborType& neighbor ) - { - _timer.start(); - using neighbor_list_type = typename NeighborType::list_type; - auto neigh_list = neighbor.list(); - const auto mu = neighbor.brokenBonds(); - - const auto x = particles.sliceReferencePosition(); - auto u = particles.sliceDisplacement(); - const auto vol = particles.sliceVolume(); - auto m = particles.sliceWeightedVolume(); - auto theta = particles.sliceDilatation(); - auto model = _model; - - Cabana::deep_copy( theta, 0.0 ); - - auto dilatation = KOKKOS_LAMBDA( const int i ) - { - std::size_t num_neighbors = - Cabana::NeighborList::numNeighbor( - neigh_list, i ); - for ( std::size_t n = 0; n < num_neighbors; n++ ) - { - std::size_t j = - Cabana::NeighborList::getNeighbor( - neigh_list, i, n ); - - // Get the bond distance, displacement, and stretch. - double xi, r, s; - getDistance( x, u, i, j, xi, r, s ); - - // Check if all bonds are broken (m=0) to avoid dividing by - // zero. Alternatively, one could check if this bond mu(i,n) is - // broken, because m=0 only occurs when all bonds are broken. - // mu is still included to account for individual bond breaking. - if ( m( i ) > 0 ) - theta( i ) += mu( i, n ) * - model.dilatation( s, xi, vol( j ), m( i ) ); - } - }; - - neighbor.iterateLinear( exec_space{}, dilatation, particles, - "CabanaPD::ForceLPSDamage::computeDilatation" ); - - _timer.stop(); - } - template void computeForceFull( ForceType& f, const PosType& x, const PosType& u, @@ -465,19 +320,22 @@ class Force> _energy_timer.stop(); return strain_energy; } + + auto timeEnergy() { return _energy_timer.time(); }; }; template class Force> - : public Force> + : public Dilatation> { protected: - using base_type = Force>; - using model_type = ForceModel; + using base_type = + Dilatation>; + using model_type = typename base_type::model_type; model_type _model; - using base_type::_energy_timer; using base_type::_timer; + Timer _energy_timer; public: // Using the default exec_space. @@ -572,6 +430,8 @@ class Force> _energy_timer.stop(); return strain_energy; } + + auto timeEnergy() { return _energy_timer.time(); }; }; } // namespace CabanaPD From bdad261367962cd7640e33194fbf7b560e999482 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 31 Jan 2025 13:34:59 -0500 Subject: [PATCH 03/30] Initial draft of density update --- src/CabanaPD_Force.hpp | 6 +-- src/CabanaPD_Types.hpp | 4 ++ src/force_models/CabanaPD_LPS.hpp | 4 +- src/force_models/CabanaPD_PMB.hpp | 65 +++++++++++++++++++++++++++++++ 4 files changed, 74 insertions(+), 5 deletions(-) diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index b0aac11a..7fd70261 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -224,7 +224,7 @@ class Dilatation 0 ) - theta( i ) += mu( i, n ) * - model.dilatation( s, xi, vol( j ), m( i ) ); + theta( i ) += mu( i, n ) * model.dilatation( + i, s, xi, vol( j ), m( i ) ); } }; diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index 18a776cb..0c602427 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -194,5 +194,9 @@ struct is_energy_output : public std::true_type { }; +struct DynamicDensity +{ +}; + } // namespace CabanaPD #endif diff --git a/src/force_models/CabanaPD_LPS.hpp b/src/force_models/CabanaPD_LPS.hpp index 84dafc84..8767c714 100644 --- a/src/force_models/CabanaPD_LPS.hpp +++ b/src/force_models/CabanaPD_LPS.hpp @@ -76,8 +76,8 @@ struct ForceModel : public BaseForceModel return influenceFunction( xi ) * xi * xi * vol; } - KOKKOS_INLINE_FUNCTION auto dilatation( const double s, const double xi, - const double vol, + KOKKOS_INLINE_FUNCTION auto dilatation( const int, const double s, + const double xi, const double vol, const double m_i ) const { double theta_i = influenceFunction( xi ) * s * xi * xi * vol; diff --git a/src/force_models/CabanaPD_PMB.hpp b/src/force_models/CabanaPD_PMB.hpp index 0cb75206..f9c15147 100644 --- a/src/force_models/CabanaPD_PMB.hpp +++ b/src/force_models/CabanaPD_PMB.hpp @@ -203,6 +203,64 @@ struct ForceModel +struct ForceModel + : public ForceModel, + ForceModel +{ + using base_type = + ForceModel; + using lps_base_type = + ForceModel; + using base_model = typename base_type::base_model; + using fracture_type = typename base_type::fracture_type; + using mechanics_type = typename base_type::mechanics_type; + using thermal_type = typename base_type::thermal_type; + + using base_type::_s_p; + using base_type::bond_break_coeff; + using base_type::c; + using base_type::delta; + using base_type::G0; + using base_type::K; + using base_type::s0; + using base_type::s_Y; + + double rho0; + DensityType rho; + + // Define which base functions to use. + using base_type::energy; + using base_type::forceCoeff; + + ForceModel( const double delta, const double K, const double G0, + const double sigma_y, const double _rho0, + const DensityType _rho ) + : base_type( delta, K, G0, sigma_y ) + , rho0( _rho0 ) + , rho( _rho ) + { + } + + // Fused density update using plastic dilatation. + KOKKOS_INLINE_FUNCTION auto dilatation( const int i, const double s, + const double xi, const double vol, + const double ) const + { + double theta_i = + 3.0 / m_i * influenceFunction( xi ) * s * xi * xi * vol; + + // Update density using previous and current dilatation. + rho( i ) = rho0 * Kokkos::exp( theta_i ) / Kokkos::exp( theta_i_n ); + + return theta_i; + } +}; + template <> struct ForceModel : public ForceModel @@ -422,6 +480,13 @@ ForceModel( ModelType, const double delta, const double K, const double _G0, -> ForceModel; +template +ForceModel( PMB, ElasticPerfectlyPlastic, DensityType rho, const double delta, + const double K, const double G0, const double sigma_y, + const double rho0 ) + -> ForceModel; + template struct ForceModel : public BaseDynamicTemperatureModel, From 4b4c43098ab6a0171bc8802f62ee89d16df52b58 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 5 Feb 2025 09:55:09 -0500 Subject: [PATCH 04/30] fixup: specific form of the PMB dilatation and density update --- src/force_models/CabanaPD_PMB.hpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/force_models/CabanaPD_PMB.hpp b/src/force_models/CabanaPD_PMB.hpp index f9c15147..67a0b348 100644 --- a/src/force_models/CabanaPD_PMB.hpp +++ b/src/force_models/CabanaPD_PMB.hpp @@ -229,6 +229,7 @@ struct ForceModel Date: Wed, 2 Apr 2025 22:24:05 -0400 Subject: [PATCH 05/30] Add force for PMB with density update --- src/force/CabanaPD_PMB.hpp | 98 +++++++++++++++++++++++++++++++ src/force_models/CabanaPD_PMB.hpp | 15 +++-- 2 files changed, 107 insertions(+), 6 deletions(-) diff --git a/src/force/CabanaPD_PMB.hpp b/src/force/CabanaPD_PMB.hpp index 9559591b..8553894d 100644 --- a/src/force/CabanaPD_PMB.hpp +++ b/src/force/CabanaPD_PMB.hpp @@ -305,6 +305,104 @@ class Force +class Force> + : public Force>, + public Dilatation> +{ + public: + // Using the default exec_space. + using exec_space = typename MemorySpace::execution_space; + using model_type = ForceModel; + using base_type = BaseForce; + + protected: + model_type _model; + + using base_type::_energy_timer; + using base_type::_timer; + + public: + Force( const model_type model ) + : base_type() + , _model( model ) + { + } + + template + void computeForceFull( ForceType& f, const PosType& x, const PosType& u, + const ParticleType& particles, + const NeighborType& neighbor ) + { + _timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); + + auto model = _model; + const auto vol = particles.sliceVolume(); + const auto nofail = particles.sliceNoFail(); + auto theta = particles.sliceDilatation(); + + auto force_full = KOKKOS_LAMBDA( const int i ) + { + std::size_t num_neighbors = + Cabana::NeighborList::numNeighbor( + neigh_list, i ); + for ( std::size_t n = 0; n < num_neighbors; n++ ) + { + double fx_i = 0.0; + double fy_i = 0.0; + double fz_i = 0.0; + + std::size_t j = + Cabana::NeighborList::getNeighbor( + neigh_list, i, n ); + + // Get the reference positions and displacements. + double xi, r, s; + double rx, ry, rz; + getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); + + model.thermalStretch( s, i, j ); + + model.updateDensity( i, theta( i ) ); + + // Break if beyond critical stretch unless in no-fail zone. + if ( model.criticalStretch( i, j, r, xi ) && !nofail( i ) && + !nofail( j ) ) + { + mu( i, n ) = 0; + } + // Else if statement is only for performance. + else if ( mu( i, n ) > 0 ) + { + const double coeff = model.forceCoeff( i, n, s, vol( j ) ); + + double muij = mu( i, n ); + fx_i = muij * coeff * rx / r; + fy_i = muij * coeff * ry / r; + fz_i = muij * coeff * rz / r; + + f( i, 0 ) += fx_i; + f( i, 1 ) += fy_i; + f( i, 2 ) += fz_i; + } + } + }; + + neighbor.iterateLinear( exec_space{}, force_full, particles, + "CabanaPD::ForcePMBDamage::computeFull" ); + _timer.stop(); + } +}; + template class Force> diff --git a/src/force_models/CabanaPD_PMB.hpp b/src/force_models/CabanaPD_PMB.hpp index 67a0b348..ec200722 100644 --- a/src/force_models/CabanaPD_PMB.hpp +++ b/src/force_models/CabanaPD_PMB.hpp @@ -248,18 +248,21 @@ struct ForceModel Date: Thu, 10 Apr 2025 18:58:55 -0400 Subject: [PATCH 06/30] Adding hip cylinder example --- examples/thermomechanics/CMakeLists.txt | 5 +- examples/thermomechanics/hip_cylinder.cpp | 188 ++++++++++++++++ examples/thermomechanics/hip_cylinder_ref.cpp | 210 ++++++++++++++++++ .../thermomechanics/inputs/hip_cylinder.json | 27 +++ .../inputs/hip_cylinder_ref.json | 21 ++ 5 files changed, 450 insertions(+), 1 deletion(-) create mode 100644 examples/thermomechanics/hip_cylinder.cpp create mode 100644 examples/thermomechanics/hip_cylinder_ref.cpp create mode 100644 examples/thermomechanics/inputs/hip_cylinder.json create mode 100644 examples/thermomechanics/inputs/hip_cylinder_ref.json diff --git a/examples/thermomechanics/CMakeLists.txt b/examples/thermomechanics/CMakeLists.txt index b2633bfa..e43e60d2 100644 --- a/examples/thermomechanics/CMakeLists.txt +++ b/examples/thermomechanics/CMakeLists.txt @@ -10,4 +10,7 @@ target_link_libraries(ThermalDeformationHeatTransfer LINK_PUBLIC CabanaPD) add_executable(ThermalDeformationHeatTransferPrenotched thermal_deformation_heat_transfer_prenotched.cpp) target_link_libraries(ThermalDeformationHeatTransferPrenotched LINK_PUBLIC CabanaPD) -install(TARGETS ThermalDeformation ThermalCrack ThermalDeformationHeatTransfer ThermalDeformationHeatTransferPrenotched DESTINATION ${CMAKE_INSTALL_BINDIR}) +add_executable(HIPCylinder hip_cylinder.cpp) +target_link_libraries(HIPCylinder LINK_PUBLIC CabanaPD) + +install(TARGETS ThermalDeformation ThermalCrack ThermalDeformationHeatTransfer ThermalDeformationHeatTransferPrenotched HIPCylinder DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp new file mode 100644 index 00000000..8a4be1a6 --- /dev/null +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -0,0 +1,188 @@ +/**************************************************************************** + * 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 + +// Simulate a cylinder under hot isostatic pressing (HIP). +void HIPCylinderExample( const std::string filename ) +{ + // ==================================================== + // Choose Kokkos spaces + // ==================================================== + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = typename exec_space::memory_space; + + // ==================================================== + // Read inputs + // ==================================================== + CabanaPD::Inputs inputs( filename ); + + // ==================================================== + // Material parameters + // ==================================================== + double rho0 = inputs["density"]; + double K = inputs["bulk_modulus"]; + // double G = inputs["shear_modulus"]; // Only for LPS. + double sc = inputs["critical_stretch"]; + double delta = inputs["horizon"]; + delta += 1e-10; + // For PMB or LPS with influence_type == 1 + double G0 = 9 * K * delta * ( sc * sc ) / 5; + // For LPS with influence_type == 0 (default) + // double G0 = 15 * K * delta * ( sc * sc ) / 8; + + // ==================================================== + // Discretization + // ==================================================== + 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 + // ==================================================== + using model_type = CabanaPD::PMB; + CabanaPD::ForceModel force_model( model_type{}, delta, K, G0 ); + + // ==================================================== + // Custom particle generation and initialization + // ==================================================== + double x_center = 0.5 * ( low_corner[0] + high_corner[0] ); + double y_center = 0.5 * ( low_corner[1] + high_corner[1] ); + double z_center = 0.5 * ( low_corner[2] + high_corner[2] ); + double Rout = inputs["cylinder_outer_radius"]; + double Rin = inputs["cylinder_inner_radius"]; + double H = inputs["cylinder_height"]; + + // Do not create particles outside given cylindrical region + auto init_op = KOKKOS_LAMBDA( const int, const double x[3] ) + { + double rsq = ( x[0] - x_center ) * ( x[0] - x_center ) + + ( x[1] - y_center ) * ( x[1] - y_center ); + if ( rsq < Rin * Rin || rsq > Rout * Rout || + x[2] > z_center + 0.5 * H || x[2] < z_center - 0.5 * H ) + return false; + return true; + }; + + // ==================================================== + // Simulation run with contact physics + // ==================================================== + if ( inputs["use_contact"] ) + { + using contact_type = CabanaPD::NormalRepulsionModel; + CabanaPD::Particles particles( + memory_space{}, contact_type{}, low_corner, high_corner, num_cells, + halo_width, Cabana::InitRandom{}, init_op, exec_space{} ); + + auto rho = particles.sliceDensity(); + auto x = particles.sliceReferencePosition(); + auto v = particles.sliceVelocity(); + auto f = particles.sliceForce(); + auto dx = particles.dx; + + double vrmax = inputs["max_radial_velocity"]; + double vrmin = inputs["min_radial_velocity"]; + double vzmax = inputs["max_vertical_velocity"]; + double zmin = z_center - 0.5 * H; + + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + // Density + rho( pid ) = rho0; + + // Velocity + double zfactor = ( ( x( pid, 2 ) - zmin ) / ( 0.5 * H ) ) - 1; + double vr = vrmax - vrmin * zfactor * zfactor; + v( pid, 0 ) = + vr * Kokkos::cos( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) ); + v( pid, 1 ) = + vr * Kokkos::sin( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) ); + v( pid, 2 ) = vzmax * zfactor; + }; + particles.updateParticles( exec_space{}, init_functor ); + + // Use contact radius and extension relative to particle spacing. + double r_c = inputs["contact_horizon_factor"]; + double r_extend = inputs["contact_horizon_extend_factor"]; + // NOTE: dx/2 is when particles first touch. + r_c *= dx[0] / 2.0; + r_extend *= dx[0]; + + contact_type contact_model( delta, r_c, r_extend, K ); + + CabanaPD::Solver solver( inputs, particles, force_model, + contact_model ); + solver.init(); + solver.run(); + } + // ==================================================== + // Simulation run without contact + // ==================================================== + else + { + CabanaPD::Particles particles( + memory_space{}, model_type{}, low_corner, high_corner, num_cells, + halo_width, Cabana::InitRandom{}, init_op, exec_space{} ); + + auto rho = particles.sliceDensity(); + auto x = particles.sliceReferencePosition(); + auto v = particles.sliceVelocity(); + auto f = particles.sliceForce(); + + double vrmax = inputs["max_radial_velocity"]; + double vrmin = inputs["min_radial_velocity"]; + double vzmax = inputs["max_vertical_velocity"]; + double zmin = z_center - 0.5 * H; + + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + // Density + rho( pid ) = rho0; + + // Velocity + double zfactor = ( ( x( pid, 2 ) - zmin ) / ( 0.5 * H ) ) - 1; + double vr = vrmax - vrmin * zfactor * zfactor; + v( pid, 0 ) = + vr * Kokkos::cos( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) ); + v( pid, 1 ) = + vr * Kokkos::sin( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) ); + v( pid, 2 ) = vzmax * zfactor; + }; + particles.updateParticles( exec_space{}, init_functor ); + + CabanaPD::Solver solver( inputs, particles, force_model ); + solver.init(); + solver.run(); + } +} + +// Initialize MPI+Kokkos. +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + HIPCylinderExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} diff --git a/examples/thermomechanics/hip_cylinder_ref.cpp b/examples/thermomechanics/hip_cylinder_ref.cpp new file mode 100644 index 00000000..2604f86a --- /dev/null +++ b/examples/thermomechanics/hip_cylinder_ref.cpp @@ -0,0 +1,210 @@ +/**************************************************************************** + * 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 + +// Simulate a cylinder under hot isostatic pressing (HIP) +void fragmentingCylinderExample( 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 parameters + // ==================================================== + double rho0 = inputs["density"]; + double K = inputs["bulk_modulus"]; + // double G = inputs["shear_modulus"]; // Only for LPS. + double sc = inputs["critical_stretch"]; + double delta = inputs["horizon"]; + delta += 1e-10; + // For PMB or LPS with influence_type == 1 + double G0 = 9 * K * delta * ( sc * sc ) / 5; + // For LPS with influence_type == 0 (default) + // double G0 = 15 * K * delta * ( sc * sc ) / 8; + + // ==================================================== + // 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 + // ==================================================== + using model_type = CabanaPD::ForceModel; + model_type force_model( delta, K, G0 ); + // using model_type = + // CabanaPD::ForceModel; + // model_type force_model( delta, K, G, G0 ); + + // ==================================================== + // Particle generation + // ==================================================== + auto particles = CabanaPD::createParticles( + exec_space(), low_corner, high_corner, num_cells, halo_width ); + + // ==================================================== + // Boundary conditions planes + // ==================================================== + double dy = particles->dx[1]; + CabanaPD::RegionBoundary plane1( + low_corner[0], high_corner[0], low_corner[1] - dy, low_corner[1] + dy, + low_corner[2], high_corner[2] ); + CabanaPD::RegionBoundary plane2( + low_corner[0], high_corner[0], high_corner[1] - dy, high_corner[1] + dy, + low_corner[2], high_corner[2] ); + std::vector> planes = { + plane1, plane2 }; + + // ==================================================== + // Custom particle initialization + // ==================================================== + double x_center = 0.5 * ( low_corner[0] + high_corner[0] ); + double y_center = 0.5 * ( low_corner[1] + high_corner[1] ); + double Rout = inputs["cylinder_outer_radius"]; + double Rin = inputs["cylinder_inner_radius"]; + + // Do not create particles outside given cylindrical region + auto init_op = KOKKOS_LAMBDA( const int, const double x[3] ) + { + double rsq = ( x[0] - x_center ) * ( x[0] - x_center ) + + ( x[1] - y_center ) * ( x[1] - y_center ); + if ( rsq < Rin * Rin || rsq > Rout * Rout ) + return false; + return true; + }; + particles->createParticles( exec_space(), init_op ); + + auto rho = particles->sliceDensity(); + /* + auto x = particles->sliceReferencePosition(); + auto v = particles->sliceVelocity(); + auto f = particles->sliceForce(); + + double vrmax = inputs["max_radial_velocity"]; + double vrmin = inputs["min_radial_velocity"]; + double vzmax = inputs["max_vertical_velocity"]; + double zmin = low_corner[2]; + + auto dx = particles->dx; + double height = inputs["system_size"][2]; + double factor = inputs["grid_perturbation_factor"]; + + using pool_type = Kokkos::Random_XorShift64_Pool; + using random_type = Kokkos::Random_XorShift64; + pool_type pool; + int seed = 456854; + pool.init( seed, particles->n_local ); + */ + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + // Density + rho( pid ) = rho0; + + /* + // Perturb particle positions + auto gen = pool.get_state(); + for ( std::size_t d = 0; d < 3; d++ ) + { + auto rand = + Kokkos::rand::draw( gen, 0.0, 1.0 ); + x( pid, d ) += ( 2.0 * rand - 1.0 ) * factor * dx[d]; + } + pool.free_state( gen ); + + // Velocity + double zfactor = ( ( x( pid, 2 ) - zmin ) / ( 0.5 * height ) ) - 1; + double vr = vrmax - vrmin * zfactor * zfactor; + v( pid, 0 ) = + vr * Kokkos::cos( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) ); + v( pid, 1 ) = + vr * Kokkos::sin( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) ); + v( pid, 2 ) = vzmax * zfactor; + */ + }; + particles->updateParticles( exec_space{}, init_functor ); + + // ==================================================== + // Boundary conditions + // ==================================================== + /* + // Create BC last to ensure ghost particles are included. + double sigma0 = inputs["traction"]; + double b0 = sigma0 / dy; + f = particles->sliceForce(); + x = particles->sliceReferencePosition(); + // Create a symmetric force BC in the y-direction. + auto bc_op = KOKKOS_LAMBDA( const int pid, const double ) + { + auto ypos = x( pid, 1 ); + auto sign = std::abs( ypos ) / ypos; + f( pid, 1 ) += b0 * sign; + }; + auto bc = createBoundaryCondition( bc_op, exec_space{}, *particles, planes, + true ); + */ + + // ==================================================== + // Imposed field + // ==================================================== + auto x = particles->sliceReferencePosition(); + f = particles->sliceForce(); + // 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 ) = temp0 + 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; + }; + auto body_term = CabanaPD::createBodyTerm( temp_func, false ); + + // ==================================================== + // Simulation run + // ==================================================== + auto cabana_pd = CabanaPD::createSolverFracture( + inputs, particles, force_model ); + cabana_pd->init(); + cabana_pd->run(); +} + +// Initialize MPI+Kokkos. +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + fragmentingCylinderExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} diff --git a/examples/thermomechanics/inputs/hip_cylinder.json b/examples/thermomechanics/inputs/hip_cylinder.json new file mode 100644 index 00000000..813acf21 --- /dev/null +++ b/examples/thermomechanics/inputs/hip_cylinder.json @@ -0,0 +1,27 @@ +{ + "num_cells" : {"value": [100, 100, 140]}, + "system_size" : {"value": [0.25, 0.25, 0.35], "unit": "m"}, + "density" : {"value": 7800, "unit": "kg/m^3"}, + "bulk_modulus" : {"value": 130e+9, "unit": "Pa"}, + "shear_modulus" : {"value": 78e+9, "unit": "Pa"}, + "critical_stretch" : {"value": 0.02}, + "horizon" : {"value": 0.003, "unit": "m"}, + "use_contact" : {"value": false}, + "contact_horizon_factor" : {"value": 0.9}, + "contact_horizon_extend_factor" : {"value": 0.01}, + "cylinder_outer_radius_old" : {"value": 0.025, "unit": "m"}, + "cylinder_inner_radius_old" : {"value": 0.02, "unit": "m"}, + "cylinder_height_old" : {"value": 0.1, "unit": "m"}, + "cylinder_outer_radius" : {"value": 0.1187, "unit": "m"}, + "cylinder_inner_radius" : {"value": 0.0679, "unit": "m"}, + "cylinder_height" : {"value": 0.3048, "unit": "m"}, + "wall_thickness" : {"value": 0.3048, "unit": "m"}, + "max_radial_velocity" : {"value": 200, "unit": "m/s"}, + "min_radial_velocity" : {"value": 50, "unit": "m/s"}, + "max_vertical_velocity" : {"value": 100, "unit": "m/s"}, + "final_time" : {"value": 2.5e-4, "unit": "s"}, + "timestep" : {"value": 1.4e-07, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.70}, + "output_frequency" : {"value": 50}, + "output_reference" : {"value": false} +} diff --git a/examples/thermomechanics/inputs/hip_cylinder_ref.json b/examples/thermomechanics/inputs/hip_cylinder_ref.json new file mode 100644 index 00000000..3087e0d9 --- /dev/null +++ b/examples/thermomechanics/inputs/hip_cylinder_ref.json @@ -0,0 +1,21 @@ +{ + "num_cells" : {"value": [51, 51, 101]}, + "system_size" : {"value": [0.254, 0.254, 0.3048], "unit": "m"}, + "grid_perturbation_factor" : {"value": 0.1}, + "density" : {"value": 7800, "unit": "kg/m^3"}, + "bulk_modulus" : {"value": 130e+9, "unit": "Pa"}, + "shear_modulus" : {"value": 78e+9, "unit": "Pa"}, + "critical_stretch" : {"value": 0.02}, + "horizon" : {"value": 0.00417462, "unit": "m"}, + "pressure" : {"value": 2e6, "unit": "Pa"}, + "cylinder_outer_radius" : {"value": 0.127, "unit": "m"}, + "cylinder_inner_radius" : {"value": 0.0635, "unit": "m"}, + "max_radial_velocity" : {"value": 200, "unit": "m/s"}, + "min_radial_velocity" : {"value": 50, "unit": "m/s"}, + "max_vertical_velocity" : {"value": 100, "unit": "m/s"}, + "final_time" : {"value": 2.5e-4, "unit": "s"}, + "timestep" : {"value": 1.7e-07, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.70}, + "output_frequency" : {"value": 50}, + "output_reference" : {"value": false} +} From 3cfeeedd985736b51855e433d27a110ea876ed1e Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 11 Apr 2025 18:06:09 -0400 Subject: [PATCH 07/30] fixup: remove contact --- examples/thermomechanics/hip_cylinder.cpp | 121 ++++++---------------- 1 file changed, 32 insertions(+), 89 deletions(-) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index 8a4be1a6..e26c6ee6 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -83,96 +83,39 @@ void HIPCylinderExample( const std::string filename ) return true; }; - // ==================================================== - // Simulation run with contact physics - // ==================================================== - if ( inputs["use_contact"] ) - { - using contact_type = CabanaPD::NormalRepulsionModel; - CabanaPD::Particles particles( - memory_space{}, contact_type{}, low_corner, high_corner, num_cells, - halo_width, Cabana::InitRandom{}, init_op, exec_space{} ); - - auto rho = particles.sliceDensity(); - auto x = particles.sliceReferencePosition(); - auto v = particles.sliceVelocity(); - auto f = particles.sliceForce(); - auto dx = particles.dx; - - double vrmax = inputs["max_radial_velocity"]; - double vrmin = inputs["min_radial_velocity"]; - double vzmax = inputs["max_vertical_velocity"]; - double zmin = z_center - 0.5 * H; - - auto init_functor = KOKKOS_LAMBDA( const int pid ) - { - // Density - rho( pid ) = rho0; - - // Velocity - double zfactor = ( ( x( pid, 2 ) - zmin ) / ( 0.5 * H ) ) - 1; - double vr = vrmax - vrmin * zfactor * zfactor; - v( pid, 0 ) = - vr * Kokkos::cos( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) ); - v( pid, 1 ) = - vr * Kokkos::sin( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) ); - v( pid, 2 ) = vzmax * zfactor; - }; - particles.updateParticles( exec_space{}, init_functor ); - - // Use contact radius and extension relative to particle spacing. - double r_c = inputs["contact_horizon_factor"]; - double r_extend = inputs["contact_horizon_extend_factor"]; - // NOTE: dx/2 is when particles first touch. - r_c *= dx[0] / 2.0; - r_extend *= dx[0]; - - contact_type contact_model( delta, r_c, r_extend, K ); - - CabanaPD::Solver solver( inputs, particles, force_model, - contact_model ); - solver.init(); - solver.run(); - } - // ==================================================== - // Simulation run without contact - // ==================================================== - else + CabanaPD::Particles particles( + memory_space{}, model_type{}, low_corner, high_corner, num_cells, + halo_width, Cabana::InitRandom{}, init_op, exec_space{} ); + + auto rho = particles.sliceDensity(); + auto x = particles.sliceReferencePosition(); + auto v = particles.sliceVelocity(); + auto f = particles.sliceForce(); + + double vrmax = inputs["max_radial_velocity"]; + double vrmin = inputs["min_radial_velocity"]; + double vzmax = inputs["max_vertical_velocity"]; + double zmin = z_center - 0.5 * H; + + auto init_functor = KOKKOS_LAMBDA( const int pid ) { - CabanaPD::Particles particles( - memory_space{}, model_type{}, low_corner, high_corner, num_cells, - halo_width, Cabana::InitRandom{}, init_op, exec_space{} ); - - auto rho = particles.sliceDensity(); - auto x = particles.sliceReferencePosition(); - auto v = particles.sliceVelocity(); - auto f = particles.sliceForce(); - - double vrmax = inputs["max_radial_velocity"]; - double vrmin = inputs["min_radial_velocity"]; - double vzmax = inputs["max_vertical_velocity"]; - double zmin = z_center - 0.5 * H; - - auto init_functor = KOKKOS_LAMBDA( const int pid ) - { - // Density - rho( pid ) = rho0; - - // Velocity - double zfactor = ( ( x( pid, 2 ) - zmin ) / ( 0.5 * H ) ) - 1; - double vr = vrmax - vrmin * zfactor * zfactor; - v( pid, 0 ) = - vr * Kokkos::cos( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) ); - v( pid, 1 ) = - vr * Kokkos::sin( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) ); - v( pid, 2 ) = vzmax * zfactor; - }; - particles.updateParticles( exec_space{}, init_functor ); - - CabanaPD::Solver solver( inputs, particles, force_model ); - solver.init(); - solver.run(); - } + // Density + rho( pid ) = rho0; + + // Velocity + double zfactor = ( ( x( pid, 2 ) - zmin ) / ( 0.5 * H ) ) - 1; + double vr = vrmax - vrmin * zfactor * zfactor; + v( pid, 0 ) = + vr * Kokkos::cos( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) ); + v( pid, 1 ) = + vr * Kokkos::sin( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) ); + v( pid, 2 ) = vzmax * zfactor; + }; + particles.updateParticles( exec_space{}, init_functor ); + + CabanaPD::Solver solver( inputs, particles, force_model ); + solver.init(); + solver.run(); } // Initialize MPI+Kokkos. From e1b3481ad208926262f1d189ac5a138a9a26f4f1 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 11 Apr 2025 18:15:11 -0400 Subject: [PATCH 08/30] fixup: output density --- src/CabanaPD_Particles.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index e57c502d..598462c4 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -738,7 +738,7 @@ class Particles Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( h5_config, "particles", MPI_COMM_WORLD, output_step, output_time, localOffset(), getPosition( use_reference ), sliceForce(), - sliceDisplacement(), sliceVelocity(), + sliceDisplacement(), sliceVelocity(), sliceDensity(), std::forward( other )... ); #else #ifdef Cabana_ENABLE_SILO @@ -746,7 +746,7 @@ class Particles writePartialRangeTimeStep( "particles", local_grid->globalGrid(), output_step, output_time, 0, localOffset(), getPosition( use_reference ), sliceForce(), - sliceDisplacement(), sliceVelocity(), + sliceDisplacement(), sliceVelocity(), sliceDensity(), std::forward( other )... ); #else From b0fbe08140a9a11384afaa74696f8a1c89878dd1 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 14 Apr 2025 15:41:42 -0400 Subject: [PATCH 09/30] Major changes for thermal-density evolution --- examples/thermomechanics/hip_cylinder.cpp | 17 +- src/CabanaPD_Force.hpp | 25 +-- src/CabanaPD_ForceModels.hpp | 27 ++- src/CabanaPD_Particles.hpp | 111 +++++++++--- src/CabanaPD_Solver.hpp | 10 +- src/CabanaPD_Types.hpp | 3 + src/force/CabanaPD_Contact.hpp | 3 +- src/force/CabanaPD_Hertzian.hpp | 3 +- src/force/CabanaPD_LPS.hpp | 27 ++- src/force/CabanaPD_PMB.hpp | 50 +++--- src/force_models/CabanaPD_Contact.hpp | 1 + src/force_models/CabanaPD_Hertzian.hpp | 2 + src/force_models/CabanaPD_LPS.hpp | 24 ++- src/force_models/CabanaPD_PMB.hpp | 204 +++++++++++++--------- unit_test/tstForce.hpp | 4 +- 15 files changed, 315 insertions(+), 196 deletions(-) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index e26c6ee6..bf8f100e 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -45,6 +45,9 @@ void HIPCylinderExample( const std::string filename ) double G0 = 9 * K * delta * ( sc * sc ) / 5; // For LPS with influence_type == 0 (default) // double G0 = 15 * K * delta * ( sc * sc ) / 8; + double sigma_y = 10.0; + double alpha = inputs["thermal_expansion_coeff"]; + double temp0 = inputs["reference_temperature"]; // ==================================================== // Discretization @@ -60,7 +63,9 @@ void HIPCylinderExample( const std::string filename ) // Force model // ==================================================== using model_type = CabanaPD::PMB; - CabanaPD::ForceModel force_model( model_type{}, delta, K, G0 ); + using mechanics_type = CabanaPD::ElasticPerfectlyPlastic; + using thermal_type = CabanaPD::TemperatureDependent; + using density_type = CabanaPD::DynamicDensity; // ==================================================== // Custom particle generation and initialization @@ -84,8 +89,8 @@ void HIPCylinderExample( const std::string filename ) }; CabanaPD::Particles particles( - memory_space{}, model_type{}, low_corner, high_corner, num_cells, - halo_width, Cabana::InitRandom{}, init_op, exec_space{} ); + memory_space{}, model_type{}, thermal_type{}, density_type{}, + low_corner, high_corner, num_cells, halo_width, init_op, exec_space{} ); auto rho = particles.sliceDensity(); auto x = particles.sliceReferencePosition(); @@ -113,6 +118,12 @@ void HIPCylinderExample( const std::string filename ) }; particles.updateParticles( exec_space{}, init_functor ); + rho = particles.sliceDensity(); + auto temp = particles.sliceTemperature(); + CabanaPD::ForceDensityModel force_model( model_type{}, mechanics_type{}, + rho, delta, K, G0, sigma_y, rho0, + temp, alpha, temp0 ); + CabanaPD::Solver solver( inputs, particles, force_model ); solver.init(); solver.run(); diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index 7fd70261..00dcbdc0 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -127,7 +127,8 @@ getLinearizedDistance( const PosType& x, const PosType& u, const int i, } // Forward declaration. -template +template class Force; template @@ -155,18 +156,14 @@ class BaseForce /****************************************************************************** Dilatation. ******************************************************************************/ -template +template class Dilatation; -template class ModelType, - class PDModelType, class MechanicsType, class... ModelParams> -class Dilatation> +template +class Dilatation { protected: - using model_type = - ModelType; + using model_type = ModelType; model_type _model; Timer _timer; @@ -236,15 +233,11 @@ class Dilatation class ModelType, - class PDModelType, class MechanicsType, class... ModelParams> -class Dilatation> +template +class Dilatation { protected: - using model_type = - ModelType; + using model_type = ModelType; model_type _model; Timer _timer; diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index 450105a4..6f5b7306 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -17,22 +17,6 @@ namespace CabanaPD { -struct BaseForceModel -{ - double delta; - - BaseForceModel( const double _delta ) - : delta( _delta ){}; - - auto cutoff() const { return delta; } - - // Only needed for models which store bond properties. - void updateBonds( const int, const int ) {} - - // No-op for temperature. - KOKKOS_INLINE_FUNCTION - void thermalStretch( double&, const int, const int ) const {} -}; template class BasePlasticity @@ -74,7 +58,11 @@ struct BaseTemperatureModel temperature = model.temperature; } - void update( const TemperatureType _temp ) { temperature = _temp; } + template + void update( const ParticleType& particles ) + { + temperature = particles.sliceTemperature(); + } // Update stretch with temperature effects. KOKKOS_INLINE_FUNCTION @@ -123,6 +111,11 @@ template struct ForceModel; +template +struct ForceDensityModel; + } // namespace CabanaPD #endif diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 598462c4..2f2e1bca 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -79,15 +79,17 @@ namespace CabanaPD { template + class OutputType = BaseOutput, class DensityType = StaticDensity, + int Dimension = 3> class Particles; template -class Particles +class Particles { public: using self_type = Particles; + BaseOutput, StaticDensity, Dimension>; using thermal_type = TemperatureIndependent; using output_type = BaseOutput; using memory_space = MemorySpace; @@ -784,15 +786,16 @@ class Particles }; template -class Particles +class Particles : public Particles + StaticDensity, Dimension> { public: using self_type = Particles; + BaseOutput, StaticDensity, Dimension>; using base_type = Particles; + BaseOutput, StaticDensity, Dimension>; using output_type = typename base_type::output_type; using thermal_type = TemperatureIndependent; using memory_space = typename base_type::memory_space; @@ -903,15 +906,15 @@ class Particles template class Particles + StaticDensity, Dimension> : public Particles + BaseOutput, StaticDensity, Dimension> { public: using self_type = Particles; + BaseOutput, StaticDensity, Dimension>; using base_type = Particles; + BaseOutput, StaticDensity, Dimension>; using thermal_type = TemperatureDependent; using output_type = typename base_type::output_type; using memory_space = typename base_type::memory_space; @@ -1016,17 +1019,19 @@ class Particles -class Particles - : public Particles +class Particles + : public Particles { // Note: no overloaded output() since there are very few cases where this // is a desired output field. public: - using self_type = - Particles; - using base_type = - Particles; + using self_type = Particles; + using base_type = Particles; using thermal_type = typename base_type::thermal_type; using output_type = typename base_type::output_type; using memory_space = typename base_type::memory_space; @@ -1106,15 +1111,16 @@ class Particles }; template -class Particles +class Particles : public Particles + StaticDensity, Dimension> { public: - using self_type = - Particles; - using base_type = - Particles; + using self_type = Particles; + using base_type = Particles; using thermal_type = typename base_type::thermal_type; using output_type = EnergyOutput; using memory_space = typename base_type::memory_space; @@ -1206,6 +1212,59 @@ class Particles aosoa_output_type _aosoa_output; }; +// Requires LPS base regardless of model type. +template +class Particles + : public Particles +{ + public: + using self_type = Particles; + using base_type = Particles; + using thermal_type = typename base_type::thermal_type; + using output_type = EnergyOutput; + using memory_space = typename base_type::memory_space; + using base_type::dim; + + // energy, damage + using output_types = Cabana::MemberTypes; + using aosoa_output_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; + + template + Particles( MemorySpace space, ModelType, TemperatureDependent temp, + DynamicDensity, Args&&... args ) + : base_type( space, LPS{}, temp, std::forward( args )... ) + { + } + + friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; +}; + /****************************************************************************** Template deduction guides. ******************************************************************************/ @@ -1238,6 +1297,12 @@ Particles( MemorySpace, ModelType, ThermalType, OutputType, -> Particles; +template +Particles( MemorySpace, ModelType, ThermalType, DensityType, Args... ) + -> Particles; + template Particles( MemorySpace, ModelType, OutputType, std::array, diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 8ae80b42..3ba0738b 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -82,6 +82,7 @@ #include #include #include +#include namespace CabanaPD { @@ -95,7 +96,10 @@ class Solver // Core module types - required for all problems. using force_model_type = ForceModelType; - using force_type = Force; + using force_type = Force; using comm_type = Comm; @@ -104,7 +108,7 @@ class Solver // Optional module types. using heat_transfer_type = HeatTransfer; - using contact_type = Force; + using contact_type = Force; using contact_model_type = ContactModelType; // Flexible module types. @@ -161,7 +165,7 @@ class Solver // Update temperature ghost size if needed. if constexpr ( is_temperature_dependent< typename force_model_type::thermal_type>::value ) - force_model.update( particles.sliceTemperature() ); + force_model.update( particles ); neighbor = std::make_shared( inputs["half_neigh"], force_model, particles ); diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index 0c602427..73576838 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -194,6 +194,9 @@ struct is_energy_output : public std::true_type { }; +struct StaticDensity +{ +}; struct DynamicDensity { }; diff --git a/src/force/CabanaPD_Contact.hpp b/src/force/CabanaPD_Contact.hpp index f6bd49d5..3ecac3bf 100644 --- a/src/force/CabanaPD_Contact.hpp +++ b/src/force/CabanaPD_Contact.hpp @@ -42,7 +42,8 @@ KOKKOS_INLINE_FUNCTION void getRelativeNormalVelocityComponents( Normal repulsion forces ******************************************************************************/ template -class Force : public BaseForce +class Force + : public BaseForce { public: using base_type = BaseForce; diff --git a/src/force/CabanaPD_Hertzian.hpp b/src/force/CabanaPD_Hertzian.hpp index e03a4aa6..37efcebf 100644 --- a/src/force/CabanaPD_Hertzian.hpp +++ b/src/force/CabanaPD_Hertzian.hpp @@ -25,7 +25,8 @@ namespace CabanaPD Normal repulsion forces ******************************************************************************/ template -class Force : public BaseForce +class Force + : public BaseForce { public: using base_type = BaseForce; diff --git a/src/force/CabanaPD_LPS.hpp b/src/force/CabanaPD_LPS.hpp index cf4a4586..867c7f17 100644 --- a/src/force/CabanaPD_LPS.hpp +++ b/src/force/CabanaPD_LPS.hpp @@ -69,13 +69,12 @@ namespace CabanaPD { -template -class Force> - : public Dilatation> +template +class Force + : public Dilatation { protected: - using base_type = - Dilatation>; + using base_type = Dilatation; using model_type = typename base_type::model_type; using base_type::_model; @@ -172,13 +171,12 @@ class Force> auto timeEnergy() { return _energy_timer.time(); }; }; -template -class Force> - : public Dilatation> +template +class Force + : public Dilatation { protected: - using base_type = - Dilatation>; + using base_type = Dilatation; using model_type = typename base_type::model_type; using base_type::_model; @@ -324,13 +322,12 @@ class Force> auto timeEnergy() { return _energy_timer.time(); }; }; -template -class Force> - : public Dilatation> +template +class Force + : public Dilatation { protected: - using base_type = - Dilatation>; + using base_type = Dilatation; using model_type = typename base_type::model_type; model_type _model; diff --git a/src/force/CabanaPD_PMB.hpp b/src/force/CabanaPD_PMB.hpp index 8553894d..091da1a2 100644 --- a/src/force/CabanaPD_PMB.hpp +++ b/src/force/CabanaPD_PMB.hpp @@ -70,16 +70,14 @@ namespace CabanaPD { -template -class Force> +template +class Force : public BaseForce { public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; - using model_type = - ForceModel; + using model_type = ModelType; using base_type = BaseForce; protected: @@ -165,15 +163,14 @@ class Force -class Force> +template +class Force : public BaseForce { public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; - using model_type = ForceModel; + using model_type = ModelType; using base_type = BaseForce; protected: @@ -305,32 +302,31 @@ class Force -class Force> - : public Force>, - public Dilatation> +template +class Force + : public Force, + public Dilatation { public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; - using model_type = ForceModel; - using base_type = BaseForce; + using model_type = ModelType; + using base_type = Force; + using dilatation_type = Dilatation; + + using dilatation_type::computeDilatation; + using dilatation_type::computeWeightedVolume; protected: - model_type _model; + using base_type::_model; using base_type::_energy_timer; using base_type::_timer; public: Force( const model_type model ) - : base_type() - , _model( model ) + : base_type( model ) + , dilatation_type( model ) { } @@ -403,16 +399,14 @@ class Force -class Force> +template +class Force : public BaseForce { public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; - using model_type = - ForceModel; + using model_type = ModelType; using base_type = BaseForce; protected: diff --git a/src/force_models/CabanaPD_Contact.hpp b/src/force_models/CabanaPD_Contact.hpp index 8e5d6d2a..54f7bb60 100644 --- a/src/force_models/CabanaPD_Contact.hpp +++ b/src/force_models/CabanaPD_Contact.hpp @@ -50,6 +50,7 @@ struct ContactModel struct NormalRepulsionModel : public ContactModel { using base_type = ContactModel; + using model_type = ContactModel; using base_model = base_type::base_model; using fracture_type = NoFracture; using thermal_type = TemperatureIndependent; diff --git a/src/force_models/CabanaPD_Hertzian.hpp b/src/force_models/CabanaPD_Hertzian.hpp index 7f24d664..5bdc6957 100644 --- a/src/force_models/CabanaPD_Hertzian.hpp +++ b/src/force_models/CabanaPD_Hertzian.hpp @@ -23,9 +23,11 @@ namespace CabanaPD struct HertzianModel : public ContactModel { using base_type = ContactModel; + using model_type = ContactModel; using base_model = base_type::base_model; using fracture_type = NoFracture; using thermal_type = TemperatureIndependent; + using density_type = StaticDensity; using base_type::radius; double nu; // Poisson's ratio diff --git a/src/force_models/CabanaPD_LPS.hpp b/src/force_models/CabanaPD_LPS.hpp index 8767c714..5691746c 100644 --- a/src/force_models/CabanaPD_LPS.hpp +++ b/src/force_models/CabanaPD_LPS.hpp @@ -20,15 +20,15 @@ namespace CabanaPD { template <> -struct ForceModel : public BaseForceModel +struct ForceModel { - using base_type = BaseForceModel; + using model_type = LPS; using base_model = LPS; using fracture_type = NoFracture; using thermal_type = TemperatureIndependent; + using density_type = StaticDensity; - using base_type::delta; - + double delta; int influence_type; double K; @@ -38,7 +38,7 @@ struct ForceModel : public BaseForceModel ForceModel( LPS, NoFracture, const double _delta, const double _K, const double _G, const int _influence = 0 ) - : base_type( _delta ) + : delta( _delta ) , influence_type( _influence ) , K( _K ) , G( _G ) @@ -48,7 +48,7 @@ struct ForceModel : public BaseForceModel ForceModel( LPS, Elastic, NoFracture, const double _delta, const double _K, const double _G, const int _influence = 0 ) - : base_type( _delta ) + : delta( _delta ) , influence_type( _influence ) , K( _K ) , G( _G ) @@ -56,6 +56,15 @@ struct ForceModel : public BaseForceModel init(); } + auto cutoff() const { return delta; } + + // Only needed for models which store bond properties. + void updateBonds( const int, const int ) {} + + // No-op for temperature. + KOKKOS_INLINE_FUNCTION + void thermalStretch( double&, const int, const int ) const {} + void init() { theta_coeff = 3.0 * K - 5.0 * G; @@ -111,6 +120,7 @@ template <> struct ForceModel : public ForceModel { + using model_type = LPS; using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = Fracture; @@ -169,6 +179,7 @@ template <> struct ForceModel : public ForceModel { + using model_type = LinearLPS; using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; @@ -192,6 +203,7 @@ template <> struct ForceModel : public ForceModel { + using model_type = LinearLPS; using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; diff --git a/src/force_models/CabanaPD_PMB.hpp b/src/force_models/CabanaPD_PMB.hpp index ec200722..a8cd5fe9 100644 --- a/src/force_models/CabanaPD_PMB.hpp +++ b/src/force_models/CabanaPD_PMB.hpp @@ -22,27 +22,26 @@ namespace CabanaPD { template <> struct ForceModel - : public BaseForceModel { - using base_type = BaseForceModel; + using model_type = PMB; using base_model = PMB; using fracture_type = NoFracture; using thermal_type = TemperatureIndependent; + using density_type = StaticDensity; - using base_type::delta; - + double delta; double c; double K; - ForceModel( PMB, NoFracture, const double delta, const double _K ) - : base_type( delta ) + ForceModel( PMB, NoFracture, const double _delta, const double _K ) + : delta( _delta ) , K( _K ) { init(); } - ForceModel( PMB, Elastic, NoFracture, const double delta, const double _K ) - : base_type( delta ) + ForceModel( PMB, Elastic, NoFracture, const double _delta, const double _K ) + : delta( _delta ) , K( _K ) { init(); @@ -50,6 +49,15 @@ struct ForceModel void init() { c = 18.0 * K / ( pi * delta * delta * delta * delta ); } + auto cutoff() const { return delta; } + + // Only needed for models which store bond properties. + void updateBonds( const int, const int ) {} + + // No-op for temperature. + KOKKOS_INLINE_FUNCTION + void thermalStretch( double&, const int, const int ) const {} + KOKKOS_INLINE_FUNCTION auto forceCoeff( const int, const int, const double s, const double vol ) const @@ -72,6 +80,7 @@ struct ForceModel : public ForceModel { using base_type = ForceModel; + using base_type::model_type; using base_model = typename base_type::base_model; using fracture_type = Fracture; using mechanics_type = Elastic; @@ -132,6 +141,7 @@ struct ForceModel; using base_plasticity_type = BasePlasticity; + using base_type::model_type; using base_model = typename base_type::base_model; using fracture_type = Fracture; using mechanics_type = ElasticPerfectlyPlastic; @@ -203,74 +213,12 @@ struct ForceModel -struct ForceModel - : public ForceModel, - ForceModel -{ - using base_type = - ForceModel; - using lps_base_type = - ForceModel; - using base_model = typename base_type::base_model; - using fracture_type = typename base_type::fracture_type; - using mechanics_type = typename base_type::mechanics_type; - using thermal_type = typename base_type::thermal_type; - - using base_type::_s_p; - using base_type::bond_break_coeff; - using base_type::c; - using base_type::delta; - using base_type::G0; - using base_type::K; - using base_type::s0; - using base_type::s_Y; - double coeff; - - double rho0; - DensityType rho; - - // Define which base functions to use. - using base_type::energy; - using base_type::forceCoeff; - - ForceModel( const double delta, const double K, const double G0, - const double sigma_y, const double _rho0, - const DensityType _rho ) - : base_type( delta, K, G0, sigma_y ) - , rho0( _rho0 ) - , rho( _rho ) - { - coeff = 3.0 / pi / delta / delta / delta / delta; - } - - // Update plastic dilatation. - KOKKOS_INLINE_FUNCTION auto dilatation( const int, const double s, - const double xi, const double vol, - const double ) const - { - return coeff * s * xi * vol; - } - - // Density update using plastic dilatation. - KOKKOS_INLINE_FUNCTION void updateDensity( const int i, - const double theta_i ) - { - // Update density using plastic dilatation. - // Note that this assumes zero initial plastic dilatation. - rho( i ) = rho0 * Kokkos::exp( theta_i ); // exp(theta_i - theta_i_0) - } -}; - template <> struct ForceModel : public ForceModel { using base_type = ForceModel; + using model_type = LinearPMB; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using mechanics_type = Elastic; @@ -345,6 +293,7 @@ struct ForceModel; using base_temperature_type = BaseTemperatureModel; + using base_type::model_type; using base_model = PMB; using fracture_type = NoFracture; using thermal_type = TemperatureDependent; @@ -378,6 +327,7 @@ struct ForceModel using base_type = ForceModel; using base_temperature_type = BaseTemperatureModel; + using base_type::model_type; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using mechanics_type = Elastic; @@ -426,10 +376,11 @@ struct ForceModel, BaseTemperatureModel { + using memory_space = typename TemperatureType::memory_space; using base_type = ForceModel; + TemperatureIndependent, memory_space>; using base_temperature_type = BaseTemperatureModel; + using typename base_type::model_type; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using mechanics_type = ElasticPerfectlyPlastic; @@ -451,10 +402,12 @@ struct ForceModel ForceModel; -template -ForceModel( PMB, ElasticPerfectlyPlastic, DensityType rho, const double delta, - const double K, const double G0, const double sigma_y, - const double rho0 ) - -> ForceModel; - template struct ForceModel : public BaseDynamicTemperatureModel, @@ -501,6 +447,7 @@ struct ForceModel using base_type = ForceModel; using base_temperature_type = BaseDynamicTemperatureModel; + using typename base_type::model_type; using base_model = PMB; using fracture_type = NoFracture; using thermal_type = DynamicTemperature; @@ -552,6 +499,7 @@ struct ForceModel; using base_temperature_type = BaseDynamicTemperatureModel; + using typename base_type::model_type; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using thermal_type = DynamicTemperature; @@ -595,6 +543,98 @@ ForceModel( ModelType, const double delta, const double K, const double G0, -> ForceModel; +// Force models with evolving density. +template +struct ForceDensityModel + : public ForceModel, + ForceModel +{ + using base_type = ForceModel; + using lps_base_type = + ForceModel; + using typename base_type::model_type; + using base_model = typename base_type::base_model; + using fracture_type = typename base_type::fracture_type; + using mechanics_type = typename base_type::mechanics_type; + using thermal_type = typename base_type::thermal_type; + using density_type = DynamicDensity; + + using base_type::_s_p; + using base_type::bond_break_coeff; + using base_type::c; + using base_type::delta; + using base_type::G0; + using base_type::K; + using base_type::s0; + using base_type::s_Y; + double coeff; + + double rho0; + DensityType rho; + + // Define which base functions to use (do not use LPS). + using base_type::cutoff; + using base_type::energy; + using base_type::forceCoeff; + using base_type::thermalStretch; + using base_type::updateBonds; + + // Thermal parameters + using base_type::alpha; + using base_type::temp0; + using base_type::temperature; + + ForceDensityModel( PMB model, ElasticPerfectlyPlastic mechanics, + const DensityType& _rho, const double delta, + const double K, const double G0, const double sigma_y, + const double _rho0, const TemperatureType _temp, + const double alpha, const double temp0 = 0.0 ) + : base_type( model, mechanics, delta, K, G0, sigma_y, _temp, alpha, + temp0 ) + , lps_base_type( LPS{}, Fracture{}, delta, K, ( 3.0 / 5.0 * K ), G0 ) + , rho0( _rho0 ) + , rho( _rho ) + { + coeff = 3.0 / pi / delta / delta / delta / delta; + } + + // Update plastic dilatation. + KOKKOS_INLINE_FUNCTION auto dilatation( const int, const double s, + const double xi, const double vol, + const double ) const + { + return coeff * s * xi * vol; + } + + // Density update using plastic dilatation. + KOKKOS_INLINE_FUNCTION void updateDensity( const int i, + const double theta_i ) const + { + // Update density using plastic dilatation. + // Note that this assumes zero initial plastic dilatation. + rho( i ) = rho0 * Kokkos::exp( theta_i ); // exp(theta_i - theta_i_0) + } + + template + void update( const ParticleType& particles ) + { + base_type::update( particles ); + rho = particles.sliceDensity(); + } +}; + +template +ForceDensityModel( PMB, ElasticPerfectlyPlastic, DensityType rho, + const double delta, const double K, const double G0, + const double sigma_y, const double rho0, + TemperatureType temp, const double _alpha, + const double _temp0 ) + -> ForceDensityModel; + } // namespace CabanaPD #endif diff --git a/unit_test/tstForce.hpp b/unit_test/tstForce.hpp index fc053966..ddb5abfe 100644 --- a/unit_test/tstForce.hpp +++ b/unit_test/tstForce.hpp @@ -693,7 +693,9 @@ void testForce( ModelType model, const double dx, const double m, // This needs to exactly match the mesh spacing to compare with the single // particle calculation. - CabanaPD::Force force( model ); + CabanaPD::Force + force( model ); CabanaPD::Neighbor neighbor( false, model, particles ); auto x = particles.sliceReferencePosition(); From 2742c0531660dd84d13c0f3d7b4548566be33cb5 Mon Sep 17 00:00:00 2001 From: pabloseleson Date: Fri, 11 Apr 2025 19:20:53 -0400 Subject: [PATCH 10/30] Updated HIP cylinder example; new model form --- examples/thermomechanics/hip_cylinder.cpp | 131 ++++++++++++++---- .../thermomechanics/inputs/hip_cylinder.json | 15 +- 2 files changed, 112 insertions(+), 34 deletions(-) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index bf8f100e..1d050bc0 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -36,18 +36,20 @@ void HIPCylinderExample( const std::string filename ) // Material parameters // ==================================================== double rho0 = inputs["density"]; - double K = inputs["bulk_modulus"]; - // double G = inputs["shear_modulus"]; // Only for LPS. - double sc = inputs["critical_stretch"]; + double E = inputs["elastic_modulus"]; + double nu = 0.25; // Use bond-based model + double K = E / ( 3 * ( 1 - 2 * nu ) ); + double G0 = inputs["fracture_energy"]; + double sigma_y = inputs["yield_stress"]; double delta = inputs["horizon"]; delta += 1e-10; // For PMB or LPS with influence_type == 1 - double G0 = 9 * K * delta * ( sc * sc ) / 5; + // double G0 = 9 * K * delta * ( sc * sc ) / 5; // For LPS with influence_type == 0 (default) // double G0 = 15 * K * delta * ( sc * sc ) / 8; - double sigma_y = 10.0; double alpha = inputs["thermal_expansion_coeff"]; double temp0 = inputs["reference_temperature"]; + double sigma0 = inputs["traction"]; // ==================================================== // Discretization @@ -60,7 +62,7 @@ void HIPCylinderExample( const std::string filename ) int halo_width = m + 1; // Just to be safe. // ==================================================== - // Force model + // Force model type // ==================================================== using model_type = CabanaPD::PMB; using mechanics_type = CabanaPD::ElasticPerfectlyPlastic; @@ -94,39 +96,116 @@ void HIPCylinderExample( const std::string filename ) auto rho = particles.sliceDensity(); auto x = particles.sliceReferencePosition(); - auto v = particles.sliceVelocity(); - auto f = particles.sliceForce(); - - double vrmax = inputs["max_radial_velocity"]; - double vrmin = inputs["min_radial_velocity"]; - double vzmax = inputs["max_vertical_velocity"]; - double zmin = z_center - 0.5 * H; + double W = inputs["wall_thickness"]; auto init_functor = KOKKOS_LAMBDA( const int pid ) { - // Density - rho( pid ) = rho0; - - // Velocity - double zfactor = ( ( x( pid, 2 ) - zmin ) / ( 0.5 * H ) ) - 1; - double vr = vrmax - vrmin * zfactor * zfactor; - v( pid, 0 ) = - vr * Kokkos::cos( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) ); - v( pid, 1 ) = - vr * Kokkos::sin( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) ); - v( pid, 2 ) = vzmax * zfactor; + double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) + + ( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center ); + if ( rsq > ( Rin + W ) * ( Rin + W ) && + rsq > ( Rout - W ) * ( Rout - W ) && + x( pid, 2 ) > low_corner[2] + W ) + { // Powder density + rho( pid ) = 0.7 * rho0; + } + else + { // Container density + rho( pid ) = rho0; + }; }; particles.updateParticles( exec_space{}, init_functor ); + // ==================================================== + // Boundary conditions planes + // ==================================================== + CabanaPD::RegionBoundary plane( + low_corner[0], high_corner[0], low_corner[1], high_corner[1], + low_corner[2], high_corner[2] ); + + // ==================================================== + // Force model + // ==================================================== rho = particles.sliceDensity(); auto temp = particles.sliceTemperature(); CabanaPD::ForceDensityModel force_model( model_type{}, mechanics_type{}, rho, delta, K, G0, sigma_y, rho0, temp, alpha, temp0 ); + // ==================================================== + // Create solver + // ==================================================== CabanaPD::Solver solver( inputs, particles, force_model ); - solver.init(); - solver.run(); + + // ==================================================== + // Boundary conditions + // ==================================================== + + // Create BC last to ensure ghost particles are included. + double dx = particles.dx[0]; + // double dy = particles.dx[1]; + double dz = particles.dx[2]; + // double sigma0 = inputs["traction"]; + double b0 = sigma0 / dx; + auto f = particles.sliceForce(); + x = particles.sliceReferencePosition(); + // Create an isostatic pressure BC. + auto bc_op = KOKKOS_LAMBDA( const int pid, const double ) + { + double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) + + ( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center ); + double theta = + Kokkos::atan2( x( pid, 1 ) - y_center, x( pid, 0 ) - x_center ); + + // BC on outer boundary + if ( rsq > ( Rout - dx ) * ( Rout - dx ) ) + { + f( pid, 0 ) += -b0 * Kokkos::cos( theta ); + f( pid, 1 ) += -b0 * Kokkos::sin( theta ); + } + // BC on inner boundary + else if ( rsq < ( Rin + dx ) * ( Rin + dx ) ) + { + f( pid, 0 ) += b0 * Kokkos::cos( theta ); + f( pid, 1 ) += b0 * Kokkos::sin( theta ); + }; + + // BC on top boundary + if ( x( pid, 2 ) > z_center + 0.5 * H - dz ) + { + f( pid, 2 ) += -b0; + } + // BC on bottom boundary + else if ( x( pid, 2 ) < z_center - 0.5 * H + dz ) + { + f( pid, 2 ) += b0; + }; + }; + auto bc = + createBoundaryCondition( bc_op, exec_space{}, particles, true, plane ); + + // ==================================================== + // Imposed field + // ==================================================== + /* + 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 ) = temp0 + 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; + }; + CabanaPD::BodyTerm body_term( temp_func, particles.size(), false ); + */ + + // ==================================================== + // Simulation run + // ==================================================== + solver.init( bc ); + solver.run( bc ); + // solver.init( body_term ); + // solver.run( body_term ); } // Initialize MPI+Kokkos. diff --git a/examples/thermomechanics/inputs/hip_cylinder.json b/examples/thermomechanics/inputs/hip_cylinder.json index 813acf21..2b221ed2 100644 --- a/examples/thermomechanics/inputs/hip_cylinder.json +++ b/examples/thermomechanics/inputs/hip_cylinder.json @@ -1,14 +1,13 @@ { "num_cells" : {"value": [100, 100, 140]}, "system_size" : {"value": [0.25, 0.25, 0.35], "unit": "m"}, - "density" : {"value": 7800, "unit": "kg/m^3"}, - "bulk_modulus" : {"value": 130e+9, "unit": "Pa"}, - "shear_modulus" : {"value": 78e+9, "unit": "Pa"}, - "critical_stretch" : {"value": 0.02}, - "horizon" : {"value": 0.003, "unit": "m"}, - "use_contact" : {"value": false}, - "contact_horizon_factor" : {"value": 0.9}, - "contact_horizon_extend_factor" : {"value": 0.01}, + "density" : {"value": 1200, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 3e+9, "unit": "Pa"}, + "fracture_energy" : {"value": 3000, "unit": "J/m^2"}, + "horizon" : {"value": 0.003, "unit": "m"}, + "yield_stress" : {"value": 30e6, "unit": "Pa"}, + "traction" : {"value": 2e6, "unit": "Pa"}, + "reference_temperature" : {"value": 300.0, "unit": "oC"}, "cylinder_outer_radius_old" : {"value": 0.025, "unit": "m"}, "cylinder_inner_radius_old" : {"value": 0.02, "unit": "m"}, "cylinder_height_old" : {"value": 0.1, "unit": "m"}, From c9e11e313c0e8da2cfa3715a651b861f2af403a6 Mon Sep 17 00:00:00 2001 From: pabloseleson Date: Mon, 14 Apr 2025 17:13:16 -0400 Subject: [PATCH 11/30] Revised density update - need to replace rho( i ) by rho0( i ) --- src/force_models/CabanaPD_PMB.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/force_models/CabanaPD_PMB.hpp b/src/force_models/CabanaPD_PMB.hpp index a8cd5fe9..1e6d8db9 100644 --- a/src/force_models/CabanaPD_PMB.hpp +++ b/src/force_models/CabanaPD_PMB.hpp @@ -615,7 +615,8 @@ struct ForceDensityModel From c5ad0c2277312a397489d60cc729aba6190d8e90 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 14 Apr 2025 17:44:26 -0400 Subject: [PATCH 12/30] Add current density --- examples/thermomechanics/hip_cylinder.cpp | 7 +-- src/CabanaPD_Particles.hpp | 53 +++++++++++++++++++++-- src/force_models/CabanaPD_PMB.hpp | 28 +++++++----- 3 files changed, 71 insertions(+), 17 deletions(-) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index 1d050bc0..7cfc6cc4 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -126,10 +126,11 @@ void HIPCylinderExample( const std::string filename ) // Force model // ==================================================== rho = particles.sliceDensity(); + auto rho_current = particles.sliceCurrentDensity(); auto temp = particles.sliceTemperature(); - CabanaPD::ForceDensityModel force_model( model_type{}, mechanics_type{}, - rho, delta, K, G0, sigma_y, rho0, - temp, alpha, temp0 ); + CabanaPD::ForceDensityModel force_model( + model_type{}, mechanics_type{}, rho, rho_current, delta, K, G0, sigma_y, + rho0, temp, alpha, temp0 ); // ==================================================== // Create solver diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 2f2e1bca..40e58cb0 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -1226,13 +1226,12 @@ class Particles; using thermal_type = typename base_type::thermal_type; - using output_type = EnergyOutput; + using output_type = OutputType; using memory_space = typename base_type::memory_space; using base_type::dim; - // energy, damage - using output_types = Cabana::MemberTypes; - using aosoa_output_type = Cabana::AoSoA; + using aosoa_density_type = + Cabana::AoSoA, memory_space, 1>; // Per type. using base_type::n_types; @@ -1257,12 +1256,58 @@ class Particles( args )... ) { + _aosoa_density = aosoa_density_type( "Particle Output Fields", + base_type::localOffset() ); + init(); + } + + template + void createParticles( Args&&... args ) + { + // Forward arguments to standard or custom particle creation. + base_type::createParticles( std::forward( args )... ); + _aosoa_density.resize( base_type::localOffset() ); + } + + auto sliceCurrentDensity() + { + return Cabana::slice<0>( _aosoa_density, "current_density" ); + } + auto sliceCurrentDensity() const + { + return Cabana::slice<0>( _aosoa_density, "current_density" ); + } + + template + void resize( Args&&... args ) + { + base_type::resize( std::forward( args )... ); + _aosoa_density.resize( base_type::localOffset() ); + } + + template + void output( const int output_step, const double output_time, + const bool use_reference, OtherFields&&... other ) + { + base_type::output( output_step, output_time, use_reference, + sliceCurrentDensity(), + std::forward( other )... ); } friend class Comm; friend class Comm; friend class Comm; friend class Comm; + + protected: + void init() + { + auto reference = sliceCurrentDensity(); + auto current = sliceCurrentDensity(); + Cabana::deep_copy( current, reference ); + } + + aosoa_density_type _aosoa_density; }; /****************************************************************************** diff --git a/src/force_models/CabanaPD_PMB.hpp b/src/force_models/CabanaPD_PMB.hpp index 1e6d8db9..aee350a9 100644 --- a/src/force_models/CabanaPD_PMB.hpp +++ b/src/force_models/CabanaPD_PMB.hpp @@ -544,9 +544,11 @@ ForceModel( ModelType, const double delta, const double K, const double G0, TemperatureType>; // Force models with evolving density. -template +template struct ForceDensityModel + TemperatureDependent, TemperatureType, DensityType, + CurrentDensityType> : public ForceModel, ForceModel @@ -574,6 +576,7 @@ struct ForceDensityModel @@ -624,17 +629,20 @@ struct ForceDensityModel +template ForceDensityModel( PMB, ElasticPerfectlyPlastic, DensityType rho, - const double delta, const double K, const double G0, - const double sigma_y, const double rho0, - TemperatureType temp, const double _alpha, + const CurrentDensityType& rho_c, const double delta, + const double K, const double G0, const double sigma_y, + const double rho0, TemperatureType temp, const double _alpha, const double _temp0 ) -> ForceDensityModel; + TemperatureDependent, TemperatureType, DensityType, + CurrentDensityType>; } // namespace CabanaPD From 34a4c0431a620f04f70c54cd1cfb4ab34ff22675 Mon Sep 17 00:00:00 2001 From: pabloseleson Date: Mon, 14 Apr 2025 21:25:24 -0400 Subject: [PATCH 13/30] Changes to example --- examples/thermomechanics/hip_cylinder.cpp | 53 +++++++------------ .../thermomechanics/inputs/hip_cylinder.json | 8 +-- 2 files changed, 24 insertions(+), 37 deletions(-) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index 7cfc6cc4..ff9e9d17 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -43,10 +43,6 @@ void HIPCylinderExample( const std::string filename ) double sigma_y = inputs["yield_stress"]; double delta = inputs["horizon"]; delta += 1e-10; - // For PMB or LPS with influence_type == 1 - // double G0 = 9 * K * delta * ( sc * sc ) / 5; - // For LPS with influence_type == 0 (default) - // double G0 = 15 * K * delta * ( sc * sc ) / 8; double alpha = inputs["thermal_expansion_coeff"]; double temp0 = inputs["reference_temperature"]; double sigma0 = inputs["traction"]; @@ -94,17 +90,19 @@ void HIPCylinderExample( const std::string filename ) memory_space{}, model_type{}, thermal_type{}, density_type{}, low_corner, high_corner, num_cells, halo_width, init_op, exec_space{} ); + // Impose separate density values for powder and container particles. + double W = inputs["wall_thickness"]; auto rho = particles.sliceDensity(); auto x = particles.sliceReferencePosition(); - double W = inputs["wall_thickness"]; auto init_functor = KOKKOS_LAMBDA( const int pid ) { double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) + ( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center ); if ( rsq > ( Rin + W ) * ( Rin + W ) && - rsq > ( Rout - W ) * ( Rout - W ) && - x( pid, 2 ) > low_corner[2] + W ) + rsq < ( Rout - W ) * ( Rout - W ) && + x( pid, 2 ) < z_center + 0.5 * H - W && + x( pid, 2 ) > z_center - 0.5 * H + W ) { // Powder density rho( pid ) = 0.7 * rho0; } @@ -138,20 +136,22 @@ void HIPCylinderExample( const std::string filename ) CabanaPD::Solver solver( inputs, particles, force_model ); // ==================================================== - // Boundary conditions + // Impose field // ==================================================== - // Create BC last to ensure ghost particles are included. double dx = particles.dx[0]; - // double dy = particles.dx[1]; double dz = particles.dx[2]; - // double sigma0 = inputs["traction"]; double b0 = sigma0 / dx; auto f = particles.sliceForce(); x = particles.sliceReferencePosition(); - // Create an isostatic pressure BC. - auto bc_op = KOKKOS_LAMBDA( const int pid, const double ) + temp = particles.sliceTemperature(); + // 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 ) { + // --------------------- + // Isostatic pressure BC + // --------------------- double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) + ( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center ); double theta = @@ -180,33 +180,20 @@ void HIPCylinderExample( const std::string filename ) { f( pid, 2 ) += b0; }; - }; - auto bc = - createBoundaryCondition( bc_op, exec_space{}, particles, true, plane ); - // ==================================================== - // Imposed field - // ==================================================== - /* - 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 ) = temp0 + 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; + // --------------------- + // Temperature BC + // --------------------- + + // temp( pid ) = temp0 + 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; }; CabanaPD::BodyTerm body_term( temp_func, particles.size(), false ); - */ // ==================================================== // Simulation run // ==================================================== - solver.init( bc ); - solver.run( bc ); - // solver.init( body_term ); - // solver.run( body_term ); + solver.init( body_term ); + solver.run( body_term ); } // Initialize MPI+Kokkos. diff --git a/examples/thermomechanics/inputs/hip_cylinder.json b/examples/thermomechanics/inputs/hip_cylinder.json index 2b221ed2..a36b8b79 100644 --- a/examples/thermomechanics/inputs/hip_cylinder.json +++ b/examples/thermomechanics/inputs/hip_cylinder.json @@ -6,15 +6,15 @@ "fracture_energy" : {"value": 3000, "unit": "J/m^2"}, "horizon" : {"value": 0.003, "unit": "m"}, "yield_stress" : {"value": 30e6, "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)"}, "traction" : {"value": 2e6, "unit": "Pa"}, "reference_temperature" : {"value": 300.0, "unit": "oC"}, - "cylinder_outer_radius_old" : {"value": 0.025, "unit": "m"}, - "cylinder_inner_radius_old" : {"value": 0.02, "unit": "m"}, - "cylinder_height_old" : {"value": 0.1, "unit": "m"}, "cylinder_outer_radius" : {"value": 0.1187, "unit": "m"}, "cylinder_inner_radius" : {"value": 0.0679, "unit": "m"}, "cylinder_height" : {"value": 0.3048, "unit": "m"}, - "wall_thickness" : {"value": 0.3048, "unit": "m"}, + "wall_thickness" : {"value": 0.0045, "unit": "m"}, "max_radial_velocity" : {"value": 200, "unit": "m/s"}, "min_radial_velocity" : {"value": 50, "unit": "m/s"}, "max_vertical_velocity" : {"value": 100, "unit": "m/s"}, From 0dcd7ae35278cd2b926051eadc7d6133606ce3f9 Mon Sep 17 00:00:00 2001 From: pabloseleson Date: Tue, 15 Apr 2025 12:25:16 -0400 Subject: [PATCH 14/30] Improvements to HIP example --- examples/thermomechanics/hip_cylinder.cpp | 7 ++--- .../thermomechanics/inputs/hip_cylinder.json | 26 +++++++++---------- 2 files changed, 16 insertions(+), 17 deletions(-) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index ff9e9d17..f54d56a1 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -45,7 +45,8 @@ void HIPCylinderExample( const std::string filename ) delta += 1e-10; double alpha = inputs["thermal_expansion_coeff"]; double temp0 = inputs["reference_temperature"]; - double sigma0 = inputs["traction"]; + double Pmax = inputs["maximum_pressure"]; + double Tmax = inputs["maximum_temperature"]; // ==================================================== // Discretization @@ -141,7 +142,7 @@ void HIPCylinderExample( const std::string filename ) // Create BC last to ensure ghost particles are included. double dx = particles.dx[0]; double dz = particles.dx[2]; - double b0 = sigma0 / dx; + double b0 = Pmax / dx; auto f = particles.sliceForce(); x = particles.sliceReferencePosition(); temp = particles.sliceTemperature(); @@ -184,8 +185,8 @@ void HIPCylinderExample( const std::string filename ) // --------------------- // Temperature BC // --------------------- - // temp( pid ) = temp0 + 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; + temp( pid ) = Tmax; }; CabanaPD::BodyTerm body_term( temp_func, particles.size(), false ); diff --git a/examples/thermomechanics/inputs/hip_cylinder.json b/examples/thermomechanics/inputs/hip_cylinder.json index a36b8b79..a545c520 100644 --- a/examples/thermomechanics/inputs/hip_cylinder.json +++ b/examples/thermomechanics/inputs/hip_cylinder.json @@ -1,26 +1,24 @@ { "num_cells" : {"value": [100, 100, 140]}, "system_size" : {"value": [0.25, 0.25, 0.35], "unit": "m"}, - "density" : {"value": 1200, "unit": "kg/m^3"}, - "elastic_modulus" : {"value": 3e+9, "unit": "Pa"}, - "fracture_energy" : {"value": 3000, "unit": "J/m^2"}, - "horizon" : {"value": 0.003, "unit": "m"}, - "yield_stress" : {"value": 30e6, "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)"}, - "traction" : {"value": 2e6, "unit": "Pa"}, - "reference_temperature" : {"value": 300.0, "unit": "oC"}, + "density" : {"value": 7950, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 195.6e+9, "unit": "Pa"}, + "fracture_energy" : {"value": 1550e+6, "unit": "J/m^2"}, + "horizon" : {"value": 0.0075, "unit": "m"}, + "yield_stress" : {"value": 170e6, "unit": "Pa"}, + "thermal_expansion_coeff" : {"value": 14.56E-6, "unit": "oC^{-1}"}, + "thermal_conductivity" : {"value": 14.12, "unit": "W/(m.K)"}, + "specific_heat_capacity" : {"value": 492, "unit": "J/(kg.K)"}, + "reference_temperature" : {"value": 20.0, "unit": "oC"}, + "maximum_pressure" : {"value": 110e6, "unit": "Pa"}, + "maximum_temperature" : {"value": 1125, "unit": "oC"}, "cylinder_outer_radius" : {"value": 0.1187, "unit": "m"}, "cylinder_inner_radius" : {"value": 0.0679, "unit": "m"}, "cylinder_height" : {"value": 0.3048, "unit": "m"}, "wall_thickness" : {"value": 0.0045, "unit": "m"}, - "max_radial_velocity" : {"value": 200, "unit": "m/s"}, - "min_radial_velocity" : {"value": 50, "unit": "m/s"}, - "max_vertical_velocity" : {"value": 100, "unit": "m/s"}, "final_time" : {"value": 2.5e-4, "unit": "s"}, "timestep" : {"value": 1.4e-07, "unit": "s"}, "timestep_safety_factor" : {"value": 0.70}, - "output_frequency" : {"value": 50}, + "output_frequency" : {"value": 10}, "output_reference" : {"value": false} } From 683716d415d1dacf666cce932640052271109d0a Mon Sep 17 00:00:00 2001 From: pabloseleson Date: Tue, 15 Apr 2025 15:39:28 -0400 Subject: [PATCH 15/30] Removed old files and added random number generator for density initialization --- examples/thermomechanics/hip_cylinder.cpp | 11 +- examples/thermomechanics/hip_cylinder_ref.cpp | 210 ------------------ .../thermomechanics/inputs/hip_cylinder.json | 2 +- .../inputs/hip_cylinder_ref.json | 21 -- 4 files changed, 11 insertions(+), 233 deletions(-) delete mode 100644 examples/thermomechanics/hip_cylinder_ref.cpp delete mode 100644 examples/thermomechanics/inputs/hip_cylinder_ref.json diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index f54d56a1..884f6885 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -96,6 +96,10 @@ void HIPCylinderExample( const std::string filename ) auto rho = particles.sliceDensity(); auto x = particles.sliceReferencePosition(); + // Use time to seed random number generator + std::srand( std::time( nullptr ) ); + double rho_perturb_factor = 0.1; + auto init_functor = KOKKOS_LAMBDA( const int pid ) { double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) + @@ -106,6 +110,11 @@ void HIPCylinderExample( const std::string filename ) x( pid, 2 ) > z_center - 0.5 * H + W ) { // Powder density rho( pid ) = 0.7 * rho0; + // Perturb powder density + double factor = + ( 1 + ( -1 + 2 * ( (double)std::rand() / ( RAND_MAX ) ) ) * + rho_perturb_factor ); + rho( pid ) *= factor; } else { // Container density @@ -186,7 +195,7 @@ void HIPCylinderExample( const std::string filename ) // Temperature BC // --------------------- // temp( pid ) = temp0 + 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; - temp( pid ) = Tmax; + // temp( pid ) = Tmax; }; CabanaPD::BodyTerm body_term( temp_func, particles.size(), false ); diff --git a/examples/thermomechanics/hip_cylinder_ref.cpp b/examples/thermomechanics/hip_cylinder_ref.cpp deleted file mode 100644 index 2604f86a..00000000 --- a/examples/thermomechanics/hip_cylinder_ref.cpp +++ /dev/null @@ -1,210 +0,0 @@ -/**************************************************************************** - * 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 - -// Simulate a cylinder under hot isostatic pressing (HIP) -void fragmentingCylinderExample( 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 parameters - // ==================================================== - double rho0 = inputs["density"]; - double K = inputs["bulk_modulus"]; - // double G = inputs["shear_modulus"]; // Only for LPS. - double sc = inputs["critical_stretch"]; - double delta = inputs["horizon"]; - delta += 1e-10; - // For PMB or LPS with influence_type == 1 - double G0 = 9 * K * delta * ( sc * sc ) / 5; - // For LPS with influence_type == 0 (default) - // double G0 = 15 * K * delta * ( sc * sc ) / 8; - - // ==================================================== - // 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 - // ==================================================== - using model_type = CabanaPD::ForceModel; - model_type force_model( delta, K, G0 ); - // using model_type = - // CabanaPD::ForceModel; - // model_type force_model( delta, K, G, G0 ); - - // ==================================================== - // Particle generation - // ==================================================== - auto particles = CabanaPD::createParticles( - exec_space(), low_corner, high_corner, num_cells, halo_width ); - - // ==================================================== - // Boundary conditions planes - // ==================================================== - double dy = particles->dx[1]; - CabanaPD::RegionBoundary plane1( - low_corner[0], high_corner[0], low_corner[1] - dy, low_corner[1] + dy, - low_corner[2], high_corner[2] ); - CabanaPD::RegionBoundary plane2( - low_corner[0], high_corner[0], high_corner[1] - dy, high_corner[1] + dy, - low_corner[2], high_corner[2] ); - std::vector> planes = { - plane1, plane2 }; - - // ==================================================== - // Custom particle initialization - // ==================================================== - double x_center = 0.5 * ( low_corner[0] + high_corner[0] ); - double y_center = 0.5 * ( low_corner[1] + high_corner[1] ); - double Rout = inputs["cylinder_outer_radius"]; - double Rin = inputs["cylinder_inner_radius"]; - - // Do not create particles outside given cylindrical region - auto init_op = KOKKOS_LAMBDA( const int, const double x[3] ) - { - double rsq = ( x[0] - x_center ) * ( x[0] - x_center ) + - ( x[1] - y_center ) * ( x[1] - y_center ); - if ( rsq < Rin * Rin || rsq > Rout * Rout ) - return false; - return true; - }; - particles->createParticles( exec_space(), init_op ); - - auto rho = particles->sliceDensity(); - /* - auto x = particles->sliceReferencePosition(); - auto v = particles->sliceVelocity(); - auto f = particles->sliceForce(); - - double vrmax = inputs["max_radial_velocity"]; - double vrmin = inputs["min_radial_velocity"]; - double vzmax = inputs["max_vertical_velocity"]; - double zmin = low_corner[2]; - - auto dx = particles->dx; - double height = inputs["system_size"][2]; - double factor = inputs["grid_perturbation_factor"]; - - using pool_type = Kokkos::Random_XorShift64_Pool; - using random_type = Kokkos::Random_XorShift64; - pool_type pool; - int seed = 456854; - pool.init( seed, particles->n_local ); - */ - auto init_functor = KOKKOS_LAMBDA( const int pid ) - { - // Density - rho( pid ) = rho0; - - /* - // Perturb particle positions - auto gen = pool.get_state(); - for ( std::size_t d = 0; d < 3; d++ ) - { - auto rand = - Kokkos::rand::draw( gen, 0.0, 1.0 ); - x( pid, d ) += ( 2.0 * rand - 1.0 ) * factor * dx[d]; - } - pool.free_state( gen ); - - // Velocity - double zfactor = ( ( x( pid, 2 ) - zmin ) / ( 0.5 * height ) ) - 1; - double vr = vrmax - vrmin * zfactor * zfactor; - v( pid, 0 ) = - vr * Kokkos::cos( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) ); - v( pid, 1 ) = - vr * Kokkos::sin( Kokkos::atan2( x( pid, 1 ), x( pid, 0 ) ) ); - v( pid, 2 ) = vzmax * zfactor; - */ - }; - particles->updateParticles( exec_space{}, init_functor ); - - // ==================================================== - // Boundary conditions - // ==================================================== - /* - // Create BC last to ensure ghost particles are included. - double sigma0 = inputs["traction"]; - double b0 = sigma0 / dy; - f = particles->sliceForce(); - x = particles->sliceReferencePosition(); - // Create a symmetric force BC in the y-direction. - auto bc_op = KOKKOS_LAMBDA( const int pid, const double ) - { - auto ypos = x( pid, 1 ); - auto sign = std::abs( ypos ) / ypos; - f( pid, 1 ) += b0 * sign; - }; - auto bc = createBoundaryCondition( bc_op, exec_space{}, *particles, planes, - true ); - */ - - // ==================================================== - // Imposed field - // ==================================================== - auto x = particles->sliceReferencePosition(); - f = particles->sliceForce(); - // 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 ) = temp0 + 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; - }; - auto body_term = CabanaPD::createBodyTerm( temp_func, false ); - - // ==================================================== - // Simulation run - // ==================================================== - auto cabana_pd = CabanaPD::createSolverFracture( - inputs, particles, force_model ); - cabana_pd->init(); - cabana_pd->run(); -} - -// Initialize MPI+Kokkos. -int main( int argc, char* argv[] ) -{ - MPI_Init( &argc, &argv ); - Kokkos::initialize( argc, argv ); - - fragmentingCylinderExample( argv[1] ); - - Kokkos::finalize(); - MPI_Finalize(); -} diff --git a/examples/thermomechanics/inputs/hip_cylinder.json b/examples/thermomechanics/inputs/hip_cylinder.json index a545c520..9c360b61 100644 --- a/examples/thermomechanics/inputs/hip_cylinder.json +++ b/examples/thermomechanics/inputs/hip_cylinder.json @@ -17,7 +17,7 @@ "cylinder_height" : {"value": 0.3048, "unit": "m"}, "wall_thickness" : {"value": 0.0045, "unit": "m"}, "final_time" : {"value": 2.5e-4, "unit": "s"}, - "timestep" : {"value": 1.4e-07, "unit": "s"}, + "timestep" : {"value": 3.0e-07, "unit": "s"}, "timestep_safety_factor" : {"value": 0.70}, "output_frequency" : {"value": 10}, "output_reference" : {"value": false} diff --git a/examples/thermomechanics/inputs/hip_cylinder_ref.json b/examples/thermomechanics/inputs/hip_cylinder_ref.json deleted file mode 100644 index 3087e0d9..00000000 --- a/examples/thermomechanics/inputs/hip_cylinder_ref.json +++ /dev/null @@ -1,21 +0,0 @@ -{ - "num_cells" : {"value": [51, 51, 101]}, - "system_size" : {"value": [0.254, 0.254, 0.3048], "unit": "m"}, - "grid_perturbation_factor" : {"value": 0.1}, - "density" : {"value": 7800, "unit": "kg/m^3"}, - "bulk_modulus" : {"value": 130e+9, "unit": "Pa"}, - "shear_modulus" : {"value": 78e+9, "unit": "Pa"}, - "critical_stretch" : {"value": 0.02}, - "horizon" : {"value": 0.00417462, "unit": "m"}, - "pressure" : {"value": 2e6, "unit": "Pa"}, - "cylinder_outer_radius" : {"value": 0.127, "unit": "m"}, - "cylinder_inner_radius" : {"value": 0.0635, "unit": "m"}, - "max_radial_velocity" : {"value": 200, "unit": "m/s"}, - "min_radial_velocity" : {"value": 50, "unit": "m/s"}, - "max_vertical_velocity" : {"value": 100, "unit": "m/s"}, - "final_time" : {"value": 2.5e-4, "unit": "s"}, - "timestep" : {"value": 1.7e-07, "unit": "s"}, - "timestep_safety_factor" : {"value": 0.70}, - "output_frequency" : {"value": 50}, - "output_reference" : {"value": false} -} From 69d4b0e390ee33db5b42dfb584ed897db23d376d Mon Sep 17 00:00:00 2001 From: pabloseleson Date: Tue, 15 Apr 2025 17:37:57 -0400 Subject: [PATCH 16/30] Added initial temperature --- examples/thermomechanics/hip_cylinder.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index 884f6885..60f8612c 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -94,6 +94,7 @@ void HIPCylinderExample( const std::string filename ) // Impose separate density values for powder and container particles. double W = inputs["wall_thickness"]; auto rho = particles.sliceDensity(); + auto temp = particles.sliceTemperature(); auto x = particles.sliceReferencePosition(); // Use time to seed random number generator @@ -120,6 +121,8 @@ void HIPCylinderExample( const std::string filename ) { // Container density rho( pid ) = rho0; }; + // Temperature + temp( pid ) = temp0; }; particles.updateParticles( exec_space{}, init_functor ); @@ -135,7 +138,7 @@ void HIPCylinderExample( const std::string filename ) // ==================================================== rho = particles.sliceDensity(); auto rho_current = particles.sliceCurrentDensity(); - auto temp = particles.sliceTemperature(); + temp = particles.sliceTemperature(); CabanaPD::ForceDensityModel force_model( model_type{}, mechanics_type{}, rho, rho_current, delta, K, G0, sigma_y, rho0, temp, alpha, temp0 ); From f9006c3b585686482499bbaef6724f05225c5818 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 16 Apr 2025 13:37:38 -0400 Subject: [PATCH 17/30] fixup: wrong slice --- src/CabanaPD_Particles.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 40e58cb0..c07f079e 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -1302,7 +1302,7 @@ class Particles Date: Wed, 16 Apr 2025 14:34:00 -0400 Subject: [PATCH 18/30] Changes to HIP example --- examples/thermomechanics/hip_cylinder.cpp | 7 +++---- examples/thermomechanics/inputs/hip_cylinder.json | 2 +- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index 60f8612c..43f5bba2 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -160,7 +160,7 @@ void HIPCylinderExample( const std::string filename ) temp = particles.sliceTemperature(); // 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 ) + auto force_temp_func = KOKKOS_LAMBDA( const int pid, const double t ) { // --------------------- // Isostatic pressure BC @@ -197,10 +197,9 @@ void HIPCylinderExample( const std::string filename ) // --------------------- // Temperature BC // --------------------- - // temp( pid ) = temp0 + 5000.0 * ( x( pid, 1 ) - low_corner_y ) * t; - // temp( pid ) = Tmax; + temp( pid ) = Tmax; }; - CabanaPD::BodyTerm body_term( temp_func, particles.size(), false ); + CabanaPD::BodyTerm body_term( force_temp_func, particles.size(), false ); // ==================================================== // Simulation run diff --git a/examples/thermomechanics/inputs/hip_cylinder.json b/examples/thermomechanics/inputs/hip_cylinder.json index 9c360b61..fedc4228 100644 --- a/examples/thermomechanics/inputs/hip_cylinder.json +++ b/examples/thermomechanics/inputs/hip_cylinder.json @@ -9,8 +9,8 @@ "thermal_expansion_coeff" : {"value": 14.56E-6, "unit": "oC^{-1}"}, "thermal_conductivity" : {"value": 14.12, "unit": "W/(m.K)"}, "specific_heat_capacity" : {"value": 492, "unit": "J/(kg.K)"}, - "reference_temperature" : {"value": 20.0, "unit": "oC"}, "maximum_pressure" : {"value": 110e6, "unit": "Pa"}, + "reference_temperature" : {"value": 1125, "unit": "oC"}, "maximum_temperature" : {"value": 1125, "unit": "oC"}, "cylinder_outer_radius" : {"value": 0.1187, "unit": "m"}, "cylinder_inner_radius" : {"value": 0.0679, "unit": "m"}, From 73388504bc94a9528496a856e60bec2ca230d95d Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 16 Apr 2025 17:14:16 -0400 Subject: [PATCH 19/30] fixup: use solver particles for HIP example --- examples/thermomechanics/hip_cylinder.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index 43f5bba2..cf0f2d56 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -152,12 +152,12 @@ void HIPCylinderExample( const std::string filename ) // Impose field // ==================================================== // Create BC last to ensure ghost particles are included. - double dx = particles.dx[0]; - double dz = particles.dx[2]; + double dx = solver.particles.dx[0]; + double dz = solver.particles.dx[2]; double b0 = Pmax / dx; - auto f = particles.sliceForce(); - x = particles.sliceReferencePosition(); - temp = particles.sliceTemperature(); + auto f = solver.particles.sliceForce(); + x = solver.particles.sliceReferencePosition(); + temp = solver.particles.sliceTemperature(); // This is purposely delayed until after solver init so that ghosted // particles are correctly taken into account for lambda capture here. auto force_temp_func = KOKKOS_LAMBDA( const int pid, const double t ) @@ -199,7 +199,8 @@ void HIPCylinderExample( const std::string filename ) // --------------------- temp( pid ) = Tmax; }; - CabanaPD::BodyTerm body_term( force_temp_func, particles.size(), false ); + CabanaPD::BodyTerm body_term( force_temp_func, solver.particles.size(), + false ); // ==================================================== // Simulation run From f6e80d4cc807eae0542137ce3ba5705bfca87263 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 16 Apr 2025 17:14:24 -0400 Subject: [PATCH 20/30] fixup: unused --- examples/thermomechanics/hip_cylinder.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index cf0f2d56..572992d8 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -160,7 +160,7 @@ void HIPCylinderExample( const std::string filename ) temp = solver.particles.sliceTemperature(); // This is purposely delayed until after solver init so that ghosted // particles are correctly taken into account for lambda capture here. - auto force_temp_func = KOKKOS_LAMBDA( const int pid, const double t ) + auto force_temp_func = KOKKOS_LAMBDA( const int pid, const double ) { // --------------------- // Isostatic pressure BC From 6d084dc93d696a50639b52dc0ba2d06935c2fd40 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 16 Apr 2025 17:30:32 -0400 Subject: [PATCH 21/30] fixup: should be a force update BC; only works for now since the temp BC is not important --- examples/thermomechanics/hip_cylinder.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index 572992d8..b7eecad4 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -200,7 +200,7 @@ void HIPCylinderExample( const std::string filename ) temp( pid ) = Tmax; }; CabanaPD::BodyTerm body_term( force_temp_func, solver.particles.size(), - false ); + true ); // ==================================================== // Simulation run From 70355f1d29527c30a8478a5f87cd1a93628e5061 Mon Sep 17 00:00:00 2001 From: pabloseleson Date: Thu, 17 Apr 2025 16:09:12 -0400 Subject: [PATCH 22/30] Initial version of running HIP example --- examples/thermomechanics/hip_cylinder.cpp | 37 ++++++++++++++++--- .../inputs/hip_cylinder_temp.json | 25 +++++++++++++ src/force_models/CabanaPD_PMB.hpp | 2 +- 3 files changed, 57 insertions(+), 7 deletions(-) create mode 100644 examples/thermomechanics/inputs/hip_cylinder_temp.json diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index b7eecad4..6158332d 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -96,10 +96,12 @@ void HIPCylinderExample( const std::string filename ) auto rho = particles.sliceDensity(); auto temp = particles.sliceTemperature(); auto x = particles.sliceReferencePosition(); + auto nofail = particles.sliceNoFail(); // Use time to seed random number generator std::srand( std::time( nullptr ) ); - double rho_perturb_factor = 0.1; + // double rho_perturb_factor = 0.1; + double rho_perturb_factor = 0.02; auto init_functor = KOKKOS_LAMBDA( const int pid ) { @@ -110,7 +112,8 @@ void HIPCylinderExample( const std::string filename ) x( pid, 2 ) < z_center + 0.5 * H - W && x( pid, 2 ) > z_center - 0.5 * H + W ) { // Powder density - rho( pid ) = 0.7 * rho0; + // rho( pid ) = 0.7 * rho0; + rho( pid ) = 0.85 * rho0; // Perturb powder density double factor = ( 1 + ( -1 + 2 * ( (double)std::rand() / ( RAND_MAX ) ) ) * @@ -123,6 +126,8 @@ void HIPCylinderExample( const std::string filename ) }; // Temperature temp( pid ) = temp0; + // No fail + nofail( pid ) = 1; }; particles.updateParticles( exec_space{}, init_functor ); @@ -147,6 +152,22 @@ void HIPCylinderExample( const std::string filename ) // Create solver // ==================================================== CabanaPD::Solver solver( inputs, particles, force_model ); + /* + using contact_type = CabanaPD::NormalRepulsionModel; + auto dx = particles.dx; + + // Use contact radius and extension relative to particle spacing. + double r_c = inputs["contact_horizon_factor"]; + double r_extend = inputs["contact_horizon_extend_factor"]; + // NOTE: dx/2 is when particles first touch. + r_c *= dx[0] / 2.0; + r_extend *= dx[0]; + + contact_type contact_model( delta, r_c, r_extend, K ); + + CabanaPD::Solver solver( inputs, particles, force_model, + contact_model ); + */ // ==================================================== // Impose field @@ -171,25 +192,29 @@ void HIPCylinderExample( const std::string filename ) Kokkos::atan2( x( pid, 1 ) - y_center, x( pid, 0 ) - x_center ); // BC on outer boundary - if ( rsq > ( Rout - dx ) * ( Rout - dx ) ) + // if ( rsq > ( Rout - dx ) * ( Rout - dx ) ) + if ( rsq > ( Rout - W ) * ( Rout - W ) ) { f( pid, 0 ) += -b0 * Kokkos::cos( theta ); f( pid, 1 ) += -b0 * Kokkos::sin( theta ); } // BC on inner boundary - else if ( rsq < ( Rin + dx ) * ( Rin + dx ) ) + // else if ( rsq < ( Rin + dx ) * ( Rin + dx ) ) + else if ( rsq < ( Rin + W ) * ( Rin + W ) ) { f( pid, 0 ) += b0 * Kokkos::cos( theta ); f( pid, 1 ) += b0 * Kokkos::sin( theta ); }; // BC on top boundary - if ( x( pid, 2 ) > z_center + 0.5 * H - dz ) + // if ( x( pid, 2 ) > z_center + 0.5 * H - dz ) + if ( x( pid, 2 ) > z_center + 0.5 * H - W ) { f( pid, 2 ) += -b0; } // BC on bottom boundary - else if ( x( pid, 2 ) < z_center - 0.5 * H + dz ) + // else if ( x( pid, 2 ) < z_center - 0.5 * H + dz ) + else if ( x( pid, 2 ) < z_center - 0.5 * H + W ) { f( pid, 2 ) += b0; }; diff --git a/examples/thermomechanics/inputs/hip_cylinder_temp.json b/examples/thermomechanics/inputs/hip_cylinder_temp.json new file mode 100644 index 00000000..f06853ae --- /dev/null +++ b/examples/thermomechanics/inputs/hip_cylinder_temp.json @@ -0,0 +1,25 @@ +{ + "num_cells" : {"value": [100, 100, 140]}, + "system_size" : {"value": [0.25, 0.25, 0.35], "unit": "m"}, + "density" : {"value": 7950, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 195.6e+8, "unit": "Pa"}, + "fracture_energy" : {"value": 1550e+6, "unit": "J/m^2"}, + "horizon" : {"value": 0.0075, "unit": "m"}, + "yield_stress_ref" : {"value": 170e7, "unit": "Pa"}, + "yield_stress" : {"value": 500e6, "unit": "Pa"}, + "thermal_expansion_coeff" : {"value": 14.56E-6, "unit": "oC^{-1}"}, + "thermal_conductivity" : {"value": 14.12, "unit": "W/(m.K)"}, + "specific_heat_capacity" : {"value": 492, "unit": "J/(kg.K)"}, + "maximum_pressure" : {"value": 110e6, "unit": "Pa"}, + "reference_temperature" : {"value": 1125, "unit": "oC"}, + "maximum_temperature" : {"value": 1125, "unit": "oC"}, + "cylinder_outer_radius" : {"value": 0.1187, "unit": "m"}, + "cylinder_inner_radius" : {"value": 0.0679, "unit": "m"}, + "cylinder_height" : {"value": 0.3048, "unit": "m"}, + "wall_thickness" : {"value": 0.0067, "unit": "m"}, + "final_time" : {"value": 1e-3, "unit": "s"}, + "timestep" : {"value": 3.0e-07, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.70}, + "output_frequency" : {"value": 50}, + "output_reference" : {"value": false} +} diff --git a/src/force_models/CabanaPD_PMB.hpp b/src/force_models/CabanaPD_PMB.hpp index aee350a9..1499d1f2 100644 --- a/src/force_models/CabanaPD_PMB.hpp +++ b/src/force_models/CabanaPD_PMB.hpp @@ -620,7 +620,7 @@ struct ForceDensityModel Date: Thu, 17 Apr 2025 16:11:01 -0400 Subject: [PATCH 23/30] Add density dependent modulus --- src/force_models/CabanaPD_PMB.hpp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/force_models/CabanaPD_PMB.hpp b/src/force_models/CabanaPD_PMB.hpp index 1499d1f2..3514aeb1 100644 --- a/src/force_models/CabanaPD_PMB.hpp +++ b/src/force_models/CabanaPD_PMB.hpp @@ -581,7 +581,6 @@ struct ForceDensityModel Date: Thu, 29 May 2025 10:20:53 -0400 Subject: [PATCH 24/30] fixup: unused --- examples/thermomechanics/hip_cylinder.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index 6158332d..130d71d0 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -174,7 +174,6 @@ void HIPCylinderExample( const std::string filename ) // ==================================================== // Create BC last to ensure ghost particles are included. double dx = solver.particles.dx[0]; - double dz = solver.particles.dx[2]; double b0 = Pmax / dx; auto f = solver.particles.sliceForce(); x = solver.particles.sliceReferencePosition(); From 1c3034c3ffe7b45a732692fe95994ac34f6f3b91 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 29 May 2025 10:21:13 -0400 Subject: [PATCH 25/30] Add density dependent c --- src/force_models/CabanaPD_PMB.hpp | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/force_models/CabanaPD_PMB.hpp b/src/force_models/CabanaPD_PMB.hpp index 3514aeb1..784e8bc8 100644 --- a/src/force_models/CabanaPD_PMB.hpp +++ b/src/force_models/CabanaPD_PMB.hpp @@ -564,9 +564,9 @@ struct ForceDensityModel Date: Thu, 19 Jun 2025 18:50:26 -0400 Subject: [PATCH 26/30] Added actual density dependence of elastic modulus based on literature --- examples/thermomechanics/hip_cylinder.cpp | 9 +++---- .../thermomechanics/inputs/hip_cylinder.json | 1 + src/force_models/CabanaPD_PMB.hpp | 27 ++++++++++++++++--- 3 files changed, 28 insertions(+), 9 deletions(-) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index 130d71d0..d8ecf06b 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -45,8 +45,6 @@ void HIPCylinderExample( const std::string filename ) delta += 1e-10; double alpha = inputs["thermal_expansion_coeff"]; double temp0 = inputs["reference_temperature"]; - double Pmax = inputs["maximum_pressure"]; - double Tmax = inputs["maximum_temperature"]; // ==================================================== // Discretization @@ -93,6 +91,7 @@ void HIPCylinderExample( const std::string filename ) // Impose separate density values for powder and container particles. double W = inputs["wall_thickness"]; + double D0 = inputs["powder_initial_relative_density"]; auto rho = particles.sliceDensity(); auto temp = particles.sliceTemperature(); auto x = particles.sliceReferencePosition(); @@ -100,7 +99,6 @@ void HIPCylinderExample( const std::string filename ) // Use time to seed random number generator std::srand( std::time( nullptr ) ); - // double rho_perturb_factor = 0.1; double rho_perturb_factor = 0.02; auto init_functor = KOKKOS_LAMBDA( const int pid ) @@ -112,8 +110,7 @@ void HIPCylinderExample( const std::string filename ) x( pid, 2 ) < z_center + 0.5 * H - W && x( pid, 2 ) > z_center - 0.5 * H + W ) { // Powder density - // rho( pid ) = 0.7 * rho0; - rho( pid ) = 0.85 * rho0; + rho( pid ) = D0 * rho0; // Perturb powder density double factor = ( 1 + ( -1 + 2 * ( (double)std::rand() / ( RAND_MAX ) ) ) * @@ -173,6 +170,8 @@ void HIPCylinderExample( const std::string filename ) // Impose field // ==================================================== // Create BC last to ensure ghost particles are included. + double Pmax = inputs["maximum_pressure"]; + double Tmax = inputs["maximum_temperature"]; double dx = solver.particles.dx[0]; double b0 = Pmax / dx; auto f = solver.particles.sliceForce(); diff --git a/examples/thermomechanics/inputs/hip_cylinder.json b/examples/thermomechanics/inputs/hip_cylinder.json index fedc4228..204cbd01 100644 --- a/examples/thermomechanics/inputs/hip_cylinder.json +++ b/examples/thermomechanics/inputs/hip_cylinder.json @@ -2,6 +2,7 @@ "num_cells" : {"value": [100, 100, 140]}, "system_size" : {"value": [0.25, 0.25, 0.35], "unit": "m"}, "density" : {"value": 7950, "unit": "kg/m^3"}, + "powder_initial_relative_density" : {"value": 0.85}, "elastic_modulus" : {"value": 195.6e+9, "unit": "Pa"}, "fracture_energy" : {"value": 1550e+6, "unit": "J/m^2"}, "horizon" : {"value": 0.0075, "unit": "m"}, diff --git a/src/force_models/CabanaPD_PMB.hpp b/src/force_models/CabanaPD_PMB.hpp index 784e8bc8..2691998c 100644 --- a/src/force_models/CabanaPD_PMB.hpp +++ b/src/force_models/CabanaPD_PMB.hpp @@ -602,21 +602,40 @@ struct ForceDensityModel Date: Wed, 25 Jun 2025 18:36:09 -0400 Subject: [PATCH 27/30] Fixed issue with initial density profile when running on GPU --- examples/thermomechanics/hip_cylinder.cpp | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index d8ecf06b..34c03fad 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -97,8 +97,14 @@ void HIPCylinderExample( const std::string filename ) auto x = particles.sliceReferencePosition(); auto nofail = particles.sliceNoFail(); + using pool_type = Kokkos::Random_XorShift64_Pool; + using random_type = Kokkos::Random_XorShift64; + pool_type pool; + int seed = 456854; + pool.init( seed, particles.numLocal() ); + // Use time to seed random number generator - std::srand( std::time( nullptr ) ); + //std::srand( std::time( nullptr ) ); double rho_perturb_factor = 0.02; auto init_functor = KOKKOS_LAMBDA( const int pid ) @@ -112,9 +118,16 @@ void HIPCylinderExample( const std::string filename ) { // Powder density rho( pid ) = D0 * rho0; // Perturb powder density - double factor = - ( 1 + ( -1 + 2 * ( (double)std::rand() / ( RAND_MAX ) ) ) * - rho_perturb_factor ); + auto gen = pool.get_state(); + auto rand = + Kokkos::rand::draw( gen, 0.0, 1.0 ); + double factor = + ( 1 + ( 2.0 * rand - 1.0 ) * + rho_perturb_factor ); + + // double factor = + // ( 1 + ( -1 + 2 * ( (double)std::rand() / ( RAND_MAX ) ) ) * + // rho_perturb_factor ); rho( pid ) *= factor; } else From 5d54b54da250ff0d65047c541e4be28e4ef95ba0 Mon Sep 17 00:00:00 2001 From: pabloseleson Date: Tue, 15 Jul 2025 19:05:26 -0400 Subject: [PATCH 28/30] Added displacement constraints to avoid translation/rotation --- examples/thermomechanics/hip_cylinder.cpp | 46 ++++++++++++------- .../thermomechanics/inputs/hip_cylinder.json | 4 +- 2 files changed, 32 insertions(+), 18 deletions(-) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index 34c03fad..0216f57d 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -104,7 +104,7 @@ void HIPCylinderExample( const std::string filename ) pool.init( seed, particles.numLocal() ); // Use time to seed random number generator - //std::srand( std::time( nullptr ) ); + // std::srand( std::time( nullptr ) ); double rho_perturb_factor = 0.02; auto init_functor = KOKKOS_LAMBDA( const int pid ) @@ -119,15 +119,13 @@ void HIPCylinderExample( const std::string filename ) rho( pid ) = D0 * rho0; // Perturb powder density auto gen = pool.get_state(); - auto rand = + auto rand = Kokkos::rand::draw( gen, 0.0, 1.0 ); - double factor = - ( 1 + ( 2.0 * rand - 1.0 ) * - rho_perturb_factor ); - - // double factor = - // ( 1 + ( -1 + 2 * ( (double)std::rand() / ( RAND_MAX ) ) ) * - // rho_perturb_factor ); + double factor = ( 1 + ( 2.0 * rand - 1.0 ) * rho_perturb_factor ); + + // double factor = + // ( 1 + ( -1 + 2 * ( (double)std::rand() / ( RAND_MAX ) ) ) * + // rho_perturb_factor ); rho( pid ) *= factor; } else @@ -186,17 +184,19 @@ void HIPCylinderExample( const std::string filename ) double Pmax = inputs["maximum_pressure"]; double Tmax = inputs["maximum_temperature"]; double dx = solver.particles.dx[0]; + double dz = solver.particles.dx[2]; double b0 = Pmax / dx; - auto f = solver.particles.sliceForce(); x = solver.particles.sliceReferencePosition(); + auto f = solver.particles.sliceForce(); + auto u = solver.particles.sliceDisplacement(); temp = solver.particles.sliceTemperature(); // This is purposely delayed until after solver init so that ghosted // particles are correctly taken into account for lambda capture here. auto force_temp_func = KOKKOS_LAMBDA( const int pid, const double ) { - // --------------------- - // Isostatic pressure BC - // --------------------- + // ----------------------- + // Isostatic pressure BC + // ----------------------- double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) + ( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center ); double theta = @@ -230,10 +230,24 @@ void HIPCylinderExample( const std::string filename ) f( pid, 2 ) += b0; }; - // --------------------- - // Temperature BC - // --------------------- + // ----------------------- + // Temperature BC + // ----------------------- temp( pid ) = Tmax; + + // ----------------------- + // Constrain displacements + // ----------------------- + double Rmid = 0.5 * ( Rin + Rout ); + + if ( x( pid, 2 ) > z_center - dz && x( pid, 2 ) < z_center + dz && + rsq > ( Rmid - dx ) * ( Rmid - dx ) && + rsq < ( Rmid + dx ) * ( Rmid + dx ) ) + { + u( pid, 0 ) = 0.0; + u( pid, 1 ) = 0.0; + u( pid, 2 ) = 0.0; + } }; CabanaPD::BodyTerm body_term( force_temp_func, solver.particles.size(), true ); diff --git a/examples/thermomechanics/inputs/hip_cylinder.json b/examples/thermomechanics/inputs/hip_cylinder.json index 204cbd01..cebb7f04 100644 --- a/examples/thermomechanics/inputs/hip_cylinder.json +++ b/examples/thermomechanics/inputs/hip_cylinder.json @@ -17,9 +17,9 @@ "cylinder_inner_radius" : {"value": 0.0679, "unit": "m"}, "cylinder_height" : {"value": 0.3048, "unit": "m"}, "wall_thickness" : {"value": 0.0045, "unit": "m"}, - "final_time" : {"value": 2.5e-4, "unit": "s"}, + "final_time" : {"value": 2.5e-2, "unit": "s"}, "timestep" : {"value": 3.0e-07, "unit": "s"}, "timestep_safety_factor" : {"value": 0.70}, - "output_frequency" : {"value": 10}, + "output_frequency" : {"value": 50}, "output_reference" : {"value": false} } From 68829155988c891310bf4f3b3f949b05360ba5b9 Mon Sep 17 00:00:00 2001 From: pabloseleson Date: Mon, 21 Jul 2025 00:14:47 -0400 Subject: [PATCH 29/30] Changes to HIP example --- examples/thermomechanics/hip_cylinder.cpp | 58 ++++++++++++++++++- .../thermomechanics/inputs/hip_cylinder.json | 2 +- 2 files changed, 58 insertions(+), 2 deletions(-) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index 0216f57d..d7e695c9 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -105,7 +105,9 @@ void HIPCylinderExample( const std::string filename ) // Use time to seed random number generator // std::srand( std::time( nullptr ) ); - double rho_perturb_factor = 0.02; + + // double rho_perturb_factor = 0.02; + double rho_perturb_factor = 0.2; auto init_functor = KOKKOS_LAMBDA( const int pid ) { @@ -127,6 +129,8 @@ void HIPCylinderExample( const std::string filename ) // ( 1 + ( -1 + 2 * ( (double)std::rand() / ( RAND_MAX ) ) ) * // rho_perturb_factor ); rho( pid ) *= factor; + // std::cout << factor << std::endl; + std::cout << "rand: " << rand << std::endl; } else { // Container density @@ -240,6 +244,43 @@ void HIPCylinderExample( const std::string filename ) // ----------------------- double Rmid = 0.5 * ( Rin + Rout ); + // Fix x- and y-displacement on top surface: only enable motion in + // z-direction if ( x( pid, 2 ) > high_corner[2] - dz && x( pid, 2 ) < + // high_corner[2] + dz ) + if ( x( pid, 2 ) > z_center + 0.5 * H - W ) + { + u( pid, 0 ) = 0.0; + u( pid, 1 ) = 0.0; + } + + // Fix x- and y-displacement on bottom surface: only enable motion in + // z-direction + if ( x( pid, 2 ) < z_center - 0.5 * H + W ) + // if ( x( pid, 2 ) > low_corner[2] - dz && x( pid, 2 ) < low_corner[2] + // + dz ) + { + u( pid, 0 ) = 0.0; + u( pid, 1 ) = 0.0; + } + + // Fix z-displacement on mid surface + if ( x( pid, 2 ) > z_center - dz && x( pid, 2 ) < z_center + dz ) + { + u( pid, 2 ) = 0.0; + } + + /* + if ( x( pid, 2 ) > low_corner[0] - dz && x( pid, 2 ) < low_corner[0] + + dz && rsq > ( Rmid - dx ) * ( Rmid - dx ) && rsq < ( Rmid + dx ) * ( + Rmid + dx ) ) + { + u( pid, 0 ) = 0.0; + u( pid, 1 ) = 0.0; + u( pid, 2 ) = 0.0; + } + */ + + /* if ( x( pid, 2 ) > z_center - dz && x( pid, 2 ) < z_center + dz && rsq > ( Rmid - dx ) * ( Rmid - dx ) && rsq < ( Rmid + dx ) * ( Rmid + dx ) ) @@ -248,6 +289,21 @@ void HIPCylinderExample( const std::string filename ) u( pid, 1 ) = 0.0; u( pid, 2 ) = 0.0; } + */ + + /* + if ( x( pid, 2 ) > z_center - dz && x( pid, 2 ) < z_center + dz ) + { + u( pid, 2 ) = 0.0; + + if ( rsq > ( Rmid - dx ) * ( Rmid - dx ) && rsq < ( Rmid + dx ) * ( + Rmid + dx ) ) + { + u( pid, 0 ) = 0.0; + u( pid, 1 ) = 0.0; + } + } + */ }; CabanaPD::BodyTerm body_term( force_temp_func, solver.particles.size(), true ); diff --git a/examples/thermomechanics/inputs/hip_cylinder.json b/examples/thermomechanics/inputs/hip_cylinder.json index cebb7f04..3e1ab066 100644 --- a/examples/thermomechanics/inputs/hip_cylinder.json +++ b/examples/thermomechanics/inputs/hip_cylinder.json @@ -20,6 +20,6 @@ "final_time" : {"value": 2.5e-2, "unit": "s"}, "timestep" : {"value": 3.0e-07, "unit": "s"}, "timestep_safety_factor" : {"value": 0.70}, - "output_frequency" : {"value": 50}, + "output_frequency" : {"value": 1}, "output_reference" : {"value": false} } From 482e2bd43591448e83b0ecc6d2f38fd99e2abf6d Mon Sep 17 00:00:00 2001 From: pabloseleson Date: Mon, 28 Jul 2025 19:14:06 -0400 Subject: [PATCH 30/30] Added demonstration example of consolidation in a cubic RVE --- examples/thermomechanics/CMakeLists.txt | 6 +- examples/thermomechanics/hip_cylinder.cpp | 95 ++--- examples/thermomechanics/hip_rve_cube.cpp | 376 ++++++++++++++++++ .../thermomechanics/inputs/hip_cylinder.json | 6 +- .../thermomechanics/inputs/hip_rve_cube.json | 24 ++ .../inputs/hip_rve_cube_ref.json | 35 ++ 6 files changed, 480 insertions(+), 62 deletions(-) create mode 100644 examples/thermomechanics/hip_rve_cube.cpp create mode 100644 examples/thermomechanics/inputs/hip_rve_cube.json create mode 100644 examples/thermomechanics/inputs/hip_rve_cube_ref.json diff --git a/examples/thermomechanics/CMakeLists.txt b/examples/thermomechanics/CMakeLists.txt index e43e60d2..19c5a103 100644 --- a/examples/thermomechanics/CMakeLists.txt +++ b/examples/thermomechanics/CMakeLists.txt @@ -13,4 +13,8 @@ target_link_libraries(ThermalDeformationHeatTransferPrenotched LINK_PUBLIC Caban add_executable(HIPCylinder hip_cylinder.cpp) target_link_libraries(HIPCylinder LINK_PUBLIC CabanaPD) -install(TARGETS ThermalDeformation ThermalCrack ThermalDeformationHeatTransfer ThermalDeformationHeatTransferPrenotched HIPCylinder DESTINATION ${CMAKE_INSTALL_BINDIR}) +add_executable(HIPRveCube hip_rve_cube.cpp) +target_link_libraries(HIPRveCube LINK_PUBLIC CabanaPD) + + +install(TARGETS ThermalDeformation ThermalCrack ThermalDeformationHeatTransfer ThermalDeformationHeatTransferPrenotched HIPCylinder HIPRveCube DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp index d7e695c9..644772d2 100644 --- a/examples/thermomechanics/hip_cylinder.cpp +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -105,9 +105,7 @@ void HIPCylinderExample( const std::string filename ) // Use time to seed random number generator // std::srand( std::time( nullptr ) ); - - // double rho_perturb_factor = 0.02; - double rho_perturb_factor = 0.2; + double rho_perturb_factor = 0.02; auto init_functor = KOKKOS_LAMBDA( const int pid ) { @@ -117,7 +115,8 @@ void HIPCylinderExample( const std::string filename ) rsq < ( Rout - W ) * ( Rout - W ) && x( pid, 2 ) < z_center + 0.5 * H - W && x( pid, 2 ) > z_center - 0.5 * H + W ) - { // Powder density + { + // Powder density rho( pid ) = D0 * rho0; // Perturb powder density auto gen = pool.get_state(); @@ -125,19 +124,28 @@ void HIPCylinderExample( const std::string filename ) Kokkos::rand::draw( gen, 0.0, 1.0 ); double factor = ( 1 + ( 2.0 * rand - 1.0 ) * rho_perturb_factor ); - // double factor = - // ( 1 + ( -1 + 2 * ( (double)std::rand() / ( RAND_MAX ) ) ) * - // rho_perturb_factor ); - rho( pid ) *= factor; - // std::cout << factor << std::endl; - std::cout << "rand: " << rand << std::endl; + // Check perturbed density does not exceed maximum density + if ( D0 * factor < 1 ) + { + rho( pid ) *= factor; + } + else + { + rho( pid ) = rho0; + } + + // Free the state after drawing + pool.free_state( gen ); } else - { // Container density + { + // Container density rho( pid ) = rho0; }; + // Temperature temp( pid ) = temp0; + // No fail nofail( pid ) = 1; }; @@ -186,18 +194,32 @@ void HIPCylinderExample( const std::string filename ) // ==================================================== // Create BC last to ensure ghost particles are included. double Pmax = inputs["maximum_pressure"]; + double trampP = inputs["ramp_pressure_time"]; double Tmax = inputs["maximum_temperature"]; double dx = solver.particles.dx[0]; double dz = solver.particles.dx[2]; - double b0 = Pmax / dx; + double b0max = Pmax / dx; x = solver.particles.sliceReferencePosition(); auto f = solver.particles.sliceForce(); auto u = solver.particles.sliceDisplacement(); temp = solver.particles.sliceTemperature(); + // This is purposely delayed until after solver init so that ghosted // particles are correctly taken into account for lambda capture here. - auto force_temp_func = KOKKOS_LAMBDA( const int pid, const double ) + auto force_temp_func = KOKKOS_LAMBDA( const int pid, const double t ) { + double b0; + + // Pressure ramping + if ( t < trampP ) + { + b0 = b0max * t / trampP; + } + else + { + b0 = b0max; + } + // ----------------------- // Isostatic pressure BC // ----------------------- @@ -207,14 +229,12 @@ void HIPCylinderExample( const std::string filename ) Kokkos::atan2( x( pid, 1 ) - y_center, x( pid, 0 ) - x_center ); // BC on outer boundary - // if ( rsq > ( Rout - dx ) * ( Rout - dx ) ) if ( rsq > ( Rout - W ) * ( Rout - W ) ) { f( pid, 0 ) += -b0 * Kokkos::cos( theta ); f( pid, 1 ) += -b0 * Kokkos::sin( theta ); } // BC on inner boundary - // else if ( rsq < ( Rin + dx ) * ( Rin + dx ) ) else if ( rsq < ( Rin + W ) * ( Rin + W ) ) { f( pid, 0 ) += b0 * Kokkos::cos( theta ); @@ -222,13 +242,11 @@ void HIPCylinderExample( const std::string filename ) }; // BC on top boundary - // if ( x( pid, 2 ) > z_center + 0.5 * H - dz ) if ( x( pid, 2 ) > z_center + 0.5 * H - W ) { f( pid, 2 ) += -b0; } // BC on bottom boundary - // else if ( x( pid, 2 ) < z_center - 0.5 * H + dz ) else if ( x( pid, 2 ) < z_center - 0.5 * H + W ) { f( pid, 2 ) += b0; @@ -242,11 +260,8 @@ void HIPCylinderExample( const std::string filename ) // ----------------------- // Constrain displacements // ----------------------- - double Rmid = 0.5 * ( Rin + Rout ); - // Fix x- and y-displacement on top surface: only enable motion in - // z-direction if ( x( pid, 2 ) > high_corner[2] - dz && x( pid, 2 ) < - // high_corner[2] + dz ) + // z-direction if ( x( pid, 2 ) > z_center + 0.5 * H - W ) { u( pid, 0 ) = 0.0; @@ -256,8 +271,6 @@ void HIPCylinderExample( const std::string filename ) // Fix x- and y-displacement on bottom surface: only enable motion in // z-direction if ( x( pid, 2 ) < z_center - 0.5 * H + W ) - // if ( x( pid, 2 ) > low_corner[2] - dz && x( pid, 2 ) < low_corner[2] - // + dz ) { u( pid, 0 ) = 0.0; u( pid, 1 ) = 0.0; @@ -268,42 +281,6 @@ void HIPCylinderExample( const std::string filename ) { u( pid, 2 ) = 0.0; } - - /* - if ( x( pid, 2 ) > low_corner[0] - dz && x( pid, 2 ) < low_corner[0] + - dz && rsq > ( Rmid - dx ) * ( Rmid - dx ) && rsq < ( Rmid + dx ) * ( - Rmid + dx ) ) - { - u( pid, 0 ) = 0.0; - u( pid, 1 ) = 0.0; - u( pid, 2 ) = 0.0; - } - */ - - /* - if ( x( pid, 2 ) > z_center - dz && x( pid, 2 ) < z_center + dz && - rsq > ( Rmid - dx ) * ( Rmid - dx ) && - rsq < ( Rmid + dx ) * ( Rmid + dx ) ) - { - u( pid, 0 ) = 0.0; - u( pid, 1 ) = 0.0; - u( pid, 2 ) = 0.0; - } - */ - - /* - if ( x( pid, 2 ) > z_center - dz && x( pid, 2 ) < z_center + dz ) - { - u( pid, 2 ) = 0.0; - - if ( rsq > ( Rmid - dx ) * ( Rmid - dx ) && rsq < ( Rmid + dx ) * ( - Rmid + dx ) ) - { - u( pid, 0 ) = 0.0; - u( pid, 1 ) = 0.0; - } - } - */ }; CabanaPD::BodyTerm body_term( force_temp_func, solver.particles.size(), true ); diff --git a/examples/thermomechanics/hip_rve_cube.cpp b/examples/thermomechanics/hip_rve_cube.cpp new file mode 100644 index 00000000..2cfb4cbe --- /dev/null +++ b/examples/thermomechanics/hip_rve_cube.cpp @@ -0,0 +1,376 @@ +/**************************************************************************** + * 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 + +// Simulate a cylinder under hot isostatic pressing (HIP). +void HIPCylinderExample( const std::string filename ) +{ + // ==================================================== + // Choose Kokkos spaces + // ==================================================== + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = typename exec_space::memory_space; + + // ==================================================== + // Read inputs + // ==================================================== + CabanaPD::Inputs inputs( filename ); + + // ==================================================== + // Material parameters + // ==================================================== + double rho0 = inputs["density"]; + double E = inputs["elastic_modulus"]; + double nu = 0.25; // Use bond-based model + double K = E / ( 3 * ( 1 - 2 * nu ) ); + double G0 = inputs["fracture_energy"]; + double sigma_y = inputs["yield_stress"]; + double delta = inputs["horizon"]; + delta += 1e-10; + double alpha = inputs["thermal_expansion_coeff"]; + double temp0 = inputs["reference_temperature"]; + + // ==================================================== + // Discretization + // ==================================================== + 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 mechanics_type = CabanaPD::ElasticPerfectlyPlastic; + using thermal_type = CabanaPD::TemperatureDependent; + using density_type = CabanaPD::DynamicDensity; + + // ==================================================== + // Custom particle generation and initialization + // ==================================================== + double x_center = 0.5 * ( low_corner[0] + high_corner[0] ); + double y_center = 0.5 * ( low_corner[1] + high_corner[1] ); + double z_center = 0.5 * ( low_corner[2] + high_corner[2] ); + // double Rout = inputs["cylinder_outer_radius"]; + // double Rin = inputs["cylinder_inner_radius"]; + // double H = inputs["cylinder_height"]; + double L = inputs["cube_side"]; + + // // Do not create particles outside given cylindrical region + // Do not create particles outside given cube region + auto init_op = KOKKOS_LAMBDA( const int, const double x[3] ) + { + /* + double rsq = ( x[0] - x_center ) * ( x[0] - x_center ) + + ( x[1] - y_center ) * ( x[1] - y_center ); + if ( rsq < Rin * Rin || rsq > Rout * Rout || + x[2] > z_center + 0.5 * H || x[2] < z_center - 0.5 * H ) + return false; + */ + + if ( Kokkos::abs( x[0] - x_center ) > 0.5 * L || + Kokkos::abs( x[1] - y_center ) > 0.5 * L || + Kokkos::abs( x[2] - z_center ) > 0.5 * L ) + return false; + return true; + }; + + CabanaPD::Particles particles( + memory_space{}, model_type{}, thermal_type{}, density_type{}, + low_corner, high_corner, num_cells, halo_width, init_op, exec_space{} ); + + // Impose separate density values for powder and container particles. + double W = inputs["wall_thickness"]; + double D0 = inputs["powder_initial_relative_density"]; + auto rho = particles.sliceDensity(); + auto temp = particles.sliceTemperature(); + auto x = particles.sliceReferencePosition(); + auto nofail = particles.sliceNoFail(); + + using pool_type = Kokkos::Random_XorShift64_Pool; + using random_type = Kokkos::Random_XorShift64; + pool_type pool; + int seed = 456854; + pool.init( seed, particles.numLocal() ); + + // Use time to seed random number generator + // std::srand( std::time( nullptr ) ); + double rho_perturb_factor = 0.02; + + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + // double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) + // + + // ( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center ); + // if ( rsq > ( Rin + W ) * ( Rin + W ) && + // rsq < ( Rout - W ) * ( Rout - W ) && + // x( pid, 2 ) < z_center + 0.5 * H - W && + // x( pid, 2 ) > z_center - 0.5 * H + W ) + + if ( Kokkos::abs( x( pid, 0 ) - x_center ) < 0.5 * L - W && + Kokkos::abs( x( pid, 1 ) - y_center ) < 0.5 * L - W && + Kokkos::abs( x( pid, 2 ) - z_center ) < 0.5 * L - W ) + { + // Powder density + rho( pid ) = D0 * rho0; + // Perturb powder density + auto gen = pool.get_state(); + auto rand = + Kokkos::rand::draw( gen, 0.0, 1.0 ); + double factor = ( 1 + ( 2.0 * rand - 1.0 ) * rho_perturb_factor ); + + // Check perturbed density does not exceed maximum density + if ( D0 * factor < 1 ) + { + rho( pid ) *= factor; + } + else + { + rho( pid ) = rho0; + } + + // Free the state after drawing + pool.free_state( gen ); + } + else + { + // Container density + rho( pid ) = rho0; + }; + + // Temperature + temp( pid ) = temp0; + + // No fail + nofail( pid ) = 1; + }; + particles.updateParticles( exec_space{}, init_functor ); + + // ==================================================== + // Boundary conditions planes + // ==================================================== + CabanaPD::RegionBoundary plane( + low_corner[0], high_corner[0], low_corner[1], high_corner[1], + low_corner[2], high_corner[2] ); + + // ==================================================== + // Force model + // ==================================================== + rho = particles.sliceDensity(); + auto rho_current = particles.sliceCurrentDensity(); + temp = particles.sliceTemperature(); + CabanaPD::ForceDensityModel force_model( + model_type{}, mechanics_type{}, rho, rho_current, delta, K, G0, sigma_y, + rho0, temp, alpha, temp0 ); + + // ==================================================== + // Create solver + // ==================================================== + CabanaPD::Solver solver( inputs, particles, force_model ); + /* + using contact_type = CabanaPD::NormalRepulsionModel; + auto dx = particles.dx; + + // Use contact radius and extension relative to particle spacing. + double r_c = inputs["contact_horizon_factor"]; + double r_extend = inputs["contact_horizon_extend_factor"]; + // NOTE: dx/2 is when particles first touch. + r_c *= dx[0] / 2.0; + r_extend *= dx[0]; + + contact_type contact_model( delta, r_c, r_extend, K ); + + CabanaPD::Solver solver( inputs, particles, force_model, + contact_model ); + */ + + // ==================================================== + // Impose field + // ==================================================== + // Create BC last to ensure ghost particles are included. + double Pmax = inputs["maximum_pressure"]; + double trampP = inputs["ramp_pressure_time"]; + double Tmax = inputs["maximum_temperature"]; + double dx = solver.particles.dx[0]; + double dz = solver.particles.dx[2]; + // double b0max = Pmax / dx; + double b0max = Pmax / W; + x = solver.particles.sliceReferencePosition(); + auto f = solver.particles.sliceForce(); + auto u = solver.particles.sliceDisplacement(); + temp = solver.particles.sliceTemperature(); + + // This is purposely delayed until after solver init so that ghosted + // particles are correctly taken into account for lambda capture here. + auto force_temp_func = KOKKOS_LAMBDA( const int pid, const double t ) + { + double b0; + + // Pressure ramping + if ( t < trampP ) + { + b0 = b0max * t / trampP; + } + else + { + b0 = b0max; + } + + // ----------------------- + // Isostatic pressure BC + // ----------------------- + // double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) + // + + // ( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center ); + // double theta = + // Kokkos::atan2( x( pid, 1 ) - y_center, x( pid, 0 ) - x_center ); + + // if ( kokkos::abs(x[0] - x_center) < 0.5*L && kokkos::abs(x[1] - + // y_center) < 0.5*L && kokkos::abs(x[2] - z_center) < 0.5*L ) + + // Pressure in x-direction + if ( x( pid, 0 ) > x_center + 0.5 * L - W ) + f( pid, 0 ) += -b0; + + if ( x( pid, 0 ) < x_center - 0.5 * L + W ) + f( pid, 0 ) += b0; + + // Pressure in y-direction + if ( x( pid, 1 ) > y_center + 0.5 * L - W ) + f( pid, 1 ) += -b0; + + if ( x( pid, 1 ) < y_center - 0.5 * L + W ) + f( pid, 1 ) += b0; + + // Pressure in z-direction + if ( x( pid, 2 ) > z_center + 0.5 * L - W ) + f( pid, 2 ) += -b0; + + if ( x( pid, 2 ) < z_center - 0.5 * L + W ) + f( pid, 2 ) += b0; + + /* + // BC on outer boundary + if ( rsq > ( Rout - W ) * ( Rout - W ) ) + { + f( pid, 0 ) += -b0 * Kokkos::cos( theta ); + f( pid, 1 ) += -b0 * Kokkos::sin( theta ); + } + // BC on inner boundary + else if ( rsq < ( Rin + W ) * ( Rin + W ) ) + { + f( pid, 0 ) += b0 * Kokkos::cos( theta ); + f( pid, 1 ) += b0 * Kokkos::sin( theta ); + }; + + // BC on top boundary + if ( x( pid, 2 ) > z_center + 0.5 * H - W ) + { + f( pid, 2 ) += -b0; + } + // BC on bottom boundary + else if ( x( pid, 2 ) < z_center - 0.5 * H + W ) + { + f( pid, 2 ) += b0; + }; + */ + + // ----------------------- + // Temperature BC + // ----------------------- + temp( pid ) = Tmax; + + // ----------------------- + // Constrain displacements + // ----------------------- + + // Only enable motion in x-direction on surfaces with x-normal + if ( x( pid, 0 ) > x_center + 0.5 * L - W || + x( pid, 0 ) < x_center - 0.5 * L + W ) + { + u( pid, 1 ) = 0.0; + u( pid, 2 ) = 0.0; + } + + // Only enable motion in y-direction on surfaces with y-normal + if ( x( pid, 1 ) > y_center + 0.5 * L - W || + x( pid, 1 ) < y_center - 0.5 * L + W ) + { + u( pid, 0 ) = 0.0; + u( pid, 2 ) = 0.0; + } + + // Only enable motion in z-direction on surfaces with z-normal + if ( x( pid, 2 ) > z_center + 0.5 * L - W || + x( pid, 2 ) < z_center - 0.5 * L + W ) + { + u( pid, 0 ) = 0.0; + u( pid, 1 ) = 0.0; + } + + /* + // Fix x- and y-displacement on top surface: only enable motion in + // z-direction + // if ( x( pid, 2 ) > z_center + 0.5 * H - W ) + if ( x( pid, 2 ) > z_center + 0.5 * L - W ) + { + u( pid, 0 ) = 0.0; + u( pid, 1 ) = 0.0; + } + + // Fix x- and y-displacement on bottom surface: only enable motion in + // z-direction + // if ( x( pid, 2 ) < z_center - 0.5 * H + W ) + if ( x( pid, 2 ) < z_center - 0.5 * L + W ) + { + u( pid, 0 ) = 0.0; + u( pid, 1 ) = 0.0; + } + + // Fix z-displacement on mid surface + if ( x( pid, 2 ) > z_center - dz && x( pid, 2 ) < z_center + dz ) + { + u( pid, 2 ) = 0.0; + } + */ + }; + CabanaPD::BodyTerm body_term( force_temp_func, solver.particles.size(), + true ); + + // ==================================================== + // Simulation run + // ==================================================== + solver.init( body_term ); + solver.run( body_term ); +} + +// Initialize MPI+Kokkos. +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + HIPCylinderExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} diff --git a/examples/thermomechanics/inputs/hip_cylinder.json b/examples/thermomechanics/inputs/hip_cylinder.json index 3e1ab066..9844b009 100644 --- a/examples/thermomechanics/inputs/hip_cylinder.json +++ b/examples/thermomechanics/inputs/hip_cylinder.json @@ -6,7 +6,8 @@ "elastic_modulus" : {"value": 195.6e+9, "unit": "Pa"}, "fracture_energy" : {"value": 1550e+6, "unit": "J/m^2"}, "horizon" : {"value": 0.0075, "unit": "m"}, - "yield_stress" : {"value": 170e6, "unit": "Pa"}, + "yield_stress_ref" : {"value": 170e6, "unit": "Pa"}, + "yield_stress" : {"value": 170e8, "unit": "Pa"}, "thermal_expansion_coeff" : {"value": 14.56E-6, "unit": "oC^{-1}"}, "thermal_conductivity" : {"value": 14.12, "unit": "W/(m.K)"}, "specific_heat_capacity" : {"value": 492, "unit": "J/(kg.K)"}, @@ -17,9 +18,10 @@ "cylinder_inner_radius" : {"value": 0.0679, "unit": "m"}, "cylinder_height" : {"value": 0.3048, "unit": "m"}, "wall_thickness" : {"value": 0.0045, "unit": "m"}, + "ramp_pressure_time" : {"value": 2.5e-1, "unit": "s"}, "final_time" : {"value": 2.5e-2, "unit": "s"}, "timestep" : {"value": 3.0e-07, "unit": "s"}, "timestep_safety_factor" : {"value": 0.70}, - "output_frequency" : {"value": 1}, + "output_frequency" : {"value": 100}, "output_reference" : {"value": false} } diff --git a/examples/thermomechanics/inputs/hip_rve_cube.json b/examples/thermomechanics/inputs/hip_rve_cube.json new file mode 100644 index 00000000..b20e409c --- /dev/null +++ b/examples/thermomechanics/inputs/hip_rve_cube.json @@ -0,0 +1,24 @@ +{ + "num_cells" : {"value": [80, 80, 80]}, + "system_size" : {"value": [0.1, 0.1, 0.1], "unit": "m"}, + "density" : {"value": 7950, "unit": "kg/m^3"}, + "powder_initial_relative_density" : {"value": 0.85}, + "elastic_modulus" : {"value": 195.6e+9, "unit": "Pa"}, + "fracture_energy" : {"value": 1550e+6, "unit": "J/m^2"}, + "horizon" : {"value": 0.00375, "unit": "m"}, + "yield_stress" : {"value": 170e8, "unit": "Pa"}, + "thermal_expansion_coeff" : {"value": 14.56E-6, "unit": "oC^{-1}"}, + "thermal_conductivity" : {"value": 14.12, "unit": "W/(m.K)"}, + "specific_heat_capacity" : {"value": 492, "unit": "J/(kg.K)"}, + "maximum_pressure" : {"value": 700e7, "unit": "Pa"}, + "reference_temperature" : {"value": 1125, "unit": "oC"}, + "maximum_temperature" : {"value": 1125, "unit": "oC"}, + "cube_side" : {"value": 0.1, "unit": "m"}, + "wall_thickness" : {"value": 0.0045, "unit": "m"}, + "ramp_pressure_time" : {"value": 6.0e-3, "unit": "s"}, + "final_time" : {"value": 0.0027, "unit": "s"}, + "timestep" : {"value": 1.5e-07, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.70}, + "output_frequency" : {"value": 250}, + "output_reference" : {"value": false} +} diff --git a/examples/thermomechanics/inputs/hip_rve_cube_ref.json b/examples/thermomechanics/inputs/hip_rve_cube_ref.json new file mode 100644 index 00000000..4a3b4e5e --- /dev/null +++ b/examples/thermomechanics/inputs/hip_rve_cube_ref.json @@ -0,0 +1,35 @@ +{ + "num_cells_OLD" : {"value": [100, 100, 140]}, + "num_cells" : {"value": [20, 20, 20]}, + "system_size_OLD" : {"value": [0.25, 0.25, 0.35], "unit": "m"}, + "system_size" : {"value": [0.1, 0.1, 0.1], "unit": "m"}, + "density" : {"value": 7950, "unit": "kg/m^3"}, + "powder_initial_relative_density" : {"value": 0.85}, + "elastic_modulus" : {"value": 195.6e+9, "unit": "Pa"}, + "fracture_energy" : {"value": 1550e+6, "unit": "J/m^2"}, + "horizon_OLD" : {"value": 0.0075, "unit": "m"}, + "horizon" : {"value": 0.015, "unit": "m"}, + "yield_stress_ref" : {"value": 170e6, "unit": "Pa"}, + "yield_stress" : {"value": 170e8, "unit": "Pa"}, + "thermal_expansion_coeff" : {"value": 14.56E-6, "unit": "oC^{-1}"}, + "thermal_conductivity" : {"value": 14.12, "unit": "W/(m.K)"}, + "specific_heat_capacity" : {"value": 492, "unit": "J/(kg.K)"}, + "maximum_pressure_OLD" : {"value": 110e6, "unit": "Pa"}, + "maximum_pressure" : {"value": 550e7, "unit": "Pa"}, + "reference_temperature" : {"value": 1125, "unit": "oC"}, + "maximum_temperature" : {"value": 1125, "unit": "oC"}, + "cube_side" : {"value": 0.1, "unit": "m"}, + "cylinder_outer_radius_OLD" : {"value": 0.1187, "unit": "m"}, + "cylinder_inner_radius_OLD" : {"value": 0.0679, "unit": "m"}, + "cylinder_heigh_OLD" : {"value": 0.3048, "unit": "m"}, + "wall_thickness" : {"value": 0.0045, "unit": "m"}, + "ramp_pressure_time_OLD" : {"value": 2.5e-1, "unit": "s"}, + "ramp_pressure_time" : {"value": 6.0e-3, "unit": "s"}, + "final_time_OLD" : {"value": 2.5e-2, "unit": "s"}, + "final_time" : {"value": 0.0054, "unit": "s"}, + "timestep_OLD" : {"value": 3.0e-07, "unit": "s"}, + "timestep" : {"value": 6.0e-07, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.70}, + "output_frequency" : {"value": 50}, + "output_reference" : {"value": false} +}