diff --git a/examples/mechanics/CMakeLists.txt b/examples/mechanics/CMakeLists.txt index a10b05f1..2671d574 100644 --- a/examples/mechanics/CMakeLists.txt +++ b/examples/mechanics/CMakeLists.txt @@ -10,4 +10,7 @@ target_link_libraries(CrackBranching LINK_PUBLIC CabanaPD) add_executable(FragmentingCylinder fragmenting_cylinder.cpp) target_link_libraries(FragmentingCylinder LINK_PUBLIC CabanaPD) -install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder DESTINATION ${CMAKE_INSTALL_BINDIR}) +add_executable(FiberReinforcedComposite fiber_reinforced_composite.cpp) +target_link_libraries(FiberReinforcedComposite LINK_PUBLIC CabanaPD) + +install(TARGETS ElasticWave KalthoffWinkler CrackBranching FragmentingCylinder FiberReinforcedComposite DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/examples/mechanics/fiber_reinforced_composite.cpp b/examples/mechanics/fiber_reinforced_composite.cpp new file mode 100644 index 00000000..f7a8737c --- /dev/null +++ b/examples/mechanics/fiber_reinforced_composite.cpp @@ -0,0 +1,233 @@ +/**************************************************************************** + * Copyright (c) 2022 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of CabanaPD. CabanaPD is distributed under a * + * BSD 3-clause license. For the licensing terms see the LICENSE file in * + * the top-level directory. * + * * + * SPDX-License-Identifier: BSD-3-Clause * + ****************************************************************************/ + +#include +#include + +#include "mpi.h" + +#include + +#include + +#include + +// Generate a unidirectional fiber-reinforced composite geometry. +void fiberReinforcedCompositeExample( 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 + // ==================================================== + // Matrix material + double rho0_m = inputs["density"][0]; + double E_m = inputs["elastic_modulus"][0]; + double nu_m = 1.0 / 3.0; + double K_m = E_m / ( 3.0 * ( 1.0 - 2.0 * nu_m ) ); + double G0_m = inputs["fracture_energy_matrix"]; + // double G_m = E_m / ( 2.0 * ( 1.0 + nu_m ) ); // Only for LPS. + + // Read horizon first because it is needed to compute the fiber fracture + // energy. + double delta = inputs["horizon"]; + delta += 1e-10; + + // Fiber material + double rho0_f = inputs["density"][1]; + double E_f = inputs["elastic_modulus"][1]; + double nu_f = 1.0 / 3.0; + double K_f = E_f / ( 3.0 * ( 1.0 - 2.0 * nu_f ) ); + double sc_f = inputs["critical_stretch_fiber"]; + double G0_f = 9 * K_f * delta * ( sc_f * sc_f ) / 5; + // double G_f = E_f / ( 2.0 * ( 1.0 + nu_f ) ); // Only for LPS. + + // ==================================================== + // 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 + // ==================================================== + // Matrix material + using model_matrix = CabanaPD::ForceModel; + model_matrix force_model_matrix( delta, K_m, G0_m ); + + // Fiber material + using model_fiber = CabanaPD::ForceModel; + model_fiber force_model_fiber( delta, K_f, G0_f ); + + // ==================================================== + // Particle generation + // ==================================================== + // Does not set displacements, velocities, etc. + auto particles = CabanaPD::createParticles( + exec_space(), low_corner, high_corner, num_cells, halo_width ); + + // ==================================================== + // Custom particle initialization + // ==================================================== + + std::array system_size = inputs["system_size"]; + + auto rho = particles->sliceDensity(); + auto x = particles->sliceReferencePosition(); + auto type = particles->sliceType(); + + // Fiber-reinforced composite geometry parameters + double Vf = inputs["fiber_volume_fraction"]; + double Df = inputs["fiber_diameter"]; + std::vector stacking_vector = inputs["stacking_sequence"]; + Kokkos::View stacking_sequence( + stacking_vector.data(), stacking_vector.size() ); + + // Fiber radius + double Rf = 0.5 * Df; + + // System sizes + double Lx = system_size[0]; + double Ly = system_size[1]; + double Lz = system_size[2]; + + // Number of plies + auto Nplies = stacking_sequence.size(); + // Ply thickness (in z-direction) + double dzply = Lz / Nplies; + + // Single-fiber volume (assume a 0° fiber orientation) + double Vfs = CabanaPD::pi * Rf * Rf * Lx; + // Domain volume + double Vd = Lx * Ly * Lz; + // Total fiber volume + double Vftotal = Vf * Vd; + // Total number of fibers + int Nf = std::floor( Vftotal / Vfs ); + // Cross section corresponding to a single fiber in the YZ-plane + // (assume all plies have 0° fiber orientation) + double Af = Ly * Lz / Nf; + // Number of fibers in y-direction (assume Af is a square area) + int Nfy = std::round( Ly / std::sqrt( Af ) ); + // Ensure Nfy is even. + if ( Nfy % 2 == 1 ) + Nfy = Nfy + 1; + + // Number of fibers in z-direction + int Nfz = std::round( Nf / Nfy ); + // Ensure number of fibers in z-direction within each ply is even. + int nfz = std::round( Nfz / Nplies ); + if ( nfz % 2 == 0 ) + { + Nfz = nfz * Nplies; + } + else + { + Nfz = ( nfz + 1 ) * Nplies; + }; + + // Fiber grid spacings (assume all plies have 0° fiber orientation) + double dyf = Ly / Nfy; + double dzf = Lz / Nfz; + + // Domain center x- and y-coordinates + double Xc = 0.5 * ( low_corner[0] + high_corner[0] ); + double Yc = 0.5 * ( low_corner[1] + high_corner[1] ); + + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + // Particle position + double xi = x( pid, 0 ); + double yi = x( pid, 1 ); + double zi = x( pid, 2 ); + + // Find ply number of particle (counting from 0). + int nply = Kokkos::floor( ( zi - low_corner[2] ) / dzply ); + + // Ply fiber orientation (in radians) + double theta = stacking_sequence( nply ) * CabanaPD::pi / 180; + + // Translate then rotate (clockwise) y-coordinate of particle in + // XY-plane. + double yinew = -Kokkos::sin( theta ) * ( xi - Xc ) + + Kokkos::cos( theta ) * ( yi - Yc ); + + // Find center of ply in z-direction (first ply has nply = 0). + double Zply_bot = low_corner[2] + nply * dzply; + double Zply_top = Zply_bot + dzply; + double Zcply = 0.5 * ( Zply_bot + Zply_top ); + + // Translate point in z-direction. + double zinew = zi - Zcply; + + // Find nearest fiber grid center point in YZ plane. + double Iyf = Kokkos::floor( yinew / dyf ); + double Izf = Kokkos::floor( zinew / dzf ); + double YI = 0.5 * dyf + dyf * Iyf; + double ZI = 0.5 * dzf + dzf * Izf; + + // Check if point belongs to fiber + if ( ( yinew - YI ) * ( yinew - YI ) + ( zinew - ZI ) * ( zinew - ZI ) < + Rf * Rf + 1e-8 ) + { + // Material type: 1 = fiber (default is 0 = matrix) + type( pid ) = 1; + // Density (fiber) + rho( pid ) = rho0_f; + } + else + { + // Density (matrix) + rho( pid ) = rho0_m; + } + }; + particles->updateParticles( exec_space{}, init_functor ); + + // ==================================================== + // Create solver + // ==================================================== + auto models = CabanaPD::createMultiForceModel( + *particles, CabanaPD::AverageTag{}, force_model_matrix, + force_model_fiber ); + auto cabana_pd = CabanaPD::createSolverFracture( + inputs, particles, models ); + + // ==================================================== + // Simulation run + // ==================================================== + cabana_pd->init(); + cabana_pd->run(); +} + +// Initialize MPI+Kokkos. +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + fiberReinforcedCompositeExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} diff --git a/examples/mechanics/inputs/fiber_reinforced_composite.json b/examples/mechanics/inputs/fiber_reinforced_composite.json new file mode 100644 index 00000000..835629a6 --- /dev/null +++ b/examples/mechanics/inputs/fiber_reinforced_composite.json @@ -0,0 +1,17 @@ +{ + "num_cells" : {"value": [201, 201, 101]}, + "system_size" : {"value": [0.25, 0.25, 0.125], "unit": "m"}, + "density" : {"value": [1200.3, 1780], "unit": "kg/m^3"}, + "elastic_modulus" : {"value": [2.97e+9, 276e+9], "unit": "Pa"}, + "critical_stretch_fiber" : {"value": 0.017}, + "fracture_energy_matrix" : {"value": 138.11, "unit": "J/m^2"}, + "horizon" : {"value": 0.003732, "unit": "m"}, + "fiber_volume_fraction" : {"value": 0.30}, + "fiber_diameter" : {"value": 0.005, "unit": "m"}, + "stacking_sequence" : {"value": [0, 45, 90]}, + "final_time" : {"value": 1.16e-6, "unit": "s"}, + "timestep" : {"value": 5.8e-7, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.85}, + "output_frequency" : {"value": 1}, + "output_reference" : {"value": true} +} diff --git a/examples/mechanics/kalthoff_winkler.cpp b/examples/mechanics/kalthoff_winkler.cpp index aed4f4e6..0e8a26e9 100644 --- a/examples/mechanics/kalthoff_winkler.cpp +++ b/examples/mechanics/kalthoff_winkler.cpp @@ -78,17 +78,16 @@ void kalthoffWinklerExample( const std::string filename ) // ==================================================== // 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 ); + using model_type1 = CabanaPD::ForceModel; + model_type1 force_model1( delta, K, G0 ); + using model_type2 = CabanaPD::ForceModel; + model_type2 force_model2( delta, K / 10.0, G0 / 10.0 ); // ==================================================== // Particle generation // ==================================================== // Does not set displacements, velocities, etc. - auto particles = CabanaPD::createParticles( + auto particles = CabanaPD::createParticles( exec_space(), low_corner, high_corner, num_cells, halo_width ); // ==================================================== @@ -98,6 +97,7 @@ void kalthoffWinklerExample( const std::string filename ) auto x = particles->sliceReferencePosition(); auto v = particles->sliceVelocity(); auto f = particles->sliceForce(); + auto type = particles->sliceType(); double dx = particles->dx[0]; double v0 = inputs["impactor_velocity"]; @@ -109,14 +109,19 @@ void kalthoffWinklerExample( const std::string filename ) if ( x( pid, 1 ) > y_prenotch1 && x( pid, 1 ) < y_prenotch2 && x( pid, 0 ) < -0.5 * height + dx ) v( pid, 0 ) = v0; + + if ( x( pid, 1 ) > 0.0 ) + type( pid ) = 1; }; particles->updateParticles( exec_space{}, init_functor ); // ==================================================== // Create solver // ==================================================== + auto models = CabanaPD::createMultiForceModel( + *particles, CabanaPD::AverageTag{}, force_model1, force_model2 ); auto cabana_pd = CabanaPD::createSolverFracture( - inputs, particles, force_model ); + inputs, particles, models ); // ==================================================== // Boundary conditions diff --git a/src/CabanaPD_Comm.hpp b/src/CabanaPD_Comm.hpp index 1a7269f0..c4a3d5cc 100644 --- a/src/CabanaPD_Comm.hpp +++ b/src/CabanaPD_Comm.hpp @@ -230,12 +230,13 @@ struct HaloIds } }; -template +template class Comm; // FIXME: extract model from ParticleType instead. template -class Comm +class Comm { public: int mpi_size = -1; @@ -322,11 +323,12 @@ class Comm }; template -class Comm - : public Comm +class Comm + : public Comm { public: - using base_type = Comm; + using base_type = + Comm; using memory_space = typename base_type::memory_space; using halo_type = typename base_type::halo_type; using base_type::gather_u; @@ -370,12 +372,51 @@ class Comm } }; -template -class Comm - : public Comm +template +class Comm + : public Comm +{ + public: + using base_type = + Comm; + using memory_space = typename base_type::memory_space; + using halo_type = typename base_type::halo_type; + using base_type::halo; + + using base_type::_init_timer; + using base_type::_timer; + + using gather_material_type = + Cabana::Gather; + std::shared_ptr gather_material; + + Comm( ParticleType& particles, int max_export_guess = 100 ) + : base_type( particles, max_export_guess ) + { + _init_timer.start(); + + gather_material = std::make_shared( + *halo, particles._aosoa_material ); + + particles.resize( halo->numLocal(), halo->numGhost() ); + _init_timer.stop(); + } + + void gatherMaterial() + { + _timer.start(); + gather_material->apply(); + _timer.stop(); + } +}; + +template +class Comm + : public Comm { public: - using base_type = Comm; + using base_type = + Comm; using memory_space = typename base_type::memory_space; using halo_type = typename base_type::halo_type; using base_type::halo; diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index 0e87bd59..9eca64aa 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -126,11 +126,12 @@ getLinearizedDistance( const PosType& x, const PosType& u, const int i, } // Forward declaration. -template +template class Force; template -class Force +class BaseForce { protected: bool _half_neigh; @@ -146,8 +147,8 @@ class Force // Primary constructor: use positions and construct neighbors. template - Force( const bool half_neigh, const double delta, - const ParticleType& particles, const double tol = 1e-14 ) + BaseForce( 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(), @@ -159,10 +160,10 @@ class Force // 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 ) + BaseForce( 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, @@ -172,7 +173,7 @@ class Force // Constructor which stores existing neighbors. template - Force( const bool half_neigh, const NeighborListType& neighbors ) + BaseForce( const bool half_neigh, const NeighborListType& neighbors ) : _half_neigh( half_neigh ) , _neigh_list( neighbors ) { diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index 26f600a8..ada80260 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -19,19 +19,161 @@ namespace CabanaPD { struct BaseForceModel { + using material_type = SingleMaterial; double delta; - BaseForceModel(){}; BaseForceModel( const double _delta ) : delta( _delta ){}; - BaseForceModel( const BaseForceModel& model ) { delta = model.delta; } - // No-op for temperature. KOKKOS_INLINE_FUNCTION - void thermalStretch( double&, const int, const int ) const {} + void thermalStretch( const int, const int, double& ) const {} +}; + +struct AverageTag +{ }; +// Wrap multiple models in a single object. +// TODO: this currently only supports bi-material systems. +template +struct ForceModels +{ + using material_type = MultiMaterial; + + using tuple_type = std::tuple; + using first_model = typename std::tuple_element<0, tuple_type>::type; + using model_type = typename first_model::model_type; + using base_model = typename first_model::base_model; + using thermal_type = typename first_model::thermal_type; + using fracture_type = typename first_model::fracture_type; + + ForceModels( MaterialType t, const ModelType... m ) + : delta( 0.0 ) + , type( t ) + , models( std::make_tuple( m... ) ) + { + setHorizon(); + } + + ForceModels( MaterialType t, const tuple_type m ) + : delta( 0.0 ) + , type( t ) + , models( m ) + { + setHorizon(); + } + + void setHorizon() + { + std::apply( [this]( auto&&... m ) { ( this->maxDelta( m ), ... ); }, + models ); + } + + template + auto maxDelta( Model m ) + { + if ( m.delta > delta ) + delta = m.delta; + } + + template + KOKKOS_INLINE_FUNCTION auto forceCoeff( const int i, const int j, + Args... args ) const + { + const int type_i = type( i ); + const int type_j = type( j ); + if ( type_i == type_j ) + if ( type_i == 0 ) + return std::get<0>( models ).forceCoeff( i, j, args... ); + else + return std::get<1>( models ).forceCoeff( i, j, args... ); + else + return std::get<2>( models ).forceCoeff( i, j, args... ); + } + + template + KOKKOS_INLINE_FUNCTION auto energy( const int i, const int j, + Args... args ) const + { + const int type_i = type( i ); + const int type_j = type( j ); + if ( type_i == type_j ) + if ( type_i == 0 ) + return std::get<0>( models ).energy( i, j, args... ); + else + return std::get<1>( models ).energy( i, j, args... ); + else + return std::get<2>( models ).energy( i, j, args... ); + } + + template + KOKKOS_INLINE_FUNCTION auto thermalStretch( const int i, const int j, + Args... args ) const + { + const int type_i = type( i ); + const int type_j = type( j ); + if ( type_i == type_j ) + if ( type_i == 0 ) + return std::get<0>( models ).thermalStretch( i, j, args... ); + else + return std::get<1>( models ).thermalStretch( i, j, args... ); + else + return std::get<2>( models ).thermalStretch( i, j, args... ); + } + + template + KOKKOS_INLINE_FUNCTION auto criticalStretch( const int i, const int j, + Args... args ) const + { + const int type_i = type( i ); + const int type_j = type( j ); + if ( type_i == type_j ) + if ( type_i == 0 ) + return std::get<0>( models ).criticalStretch( i, j, args... ); + else + return std::get<1>( models ).criticalStretch( i, j, args... ); + else + return std::get<2>( models ).criticalStretch( i, j, args... ); + } + + auto horizon( const int ) { return delta; } + auto maxHorizon() { return delta; } + + void update( const MaterialType _type ) { type = _type; } + + double delta; + MaterialType type; + tuple_type models; +}; + +template +auto createMultiForceModel( ParticleType particles, ModelType... models ) +{ + auto type = particles.sliceType(); + using material_type = decltype( type ); + return ForceModels( type, models... ); +} + +template +auto createMultiForceModel( ParticleType particles, AverageTag, + ModelType... models ) +{ + auto tuple = std::make_tuple( models... ); + auto m1 = std::get<0>( tuple ); + auto m2 = std::get<1>( tuple ); + + using first_model = + typename std::tuple_element<0, std::tuple>::type; + auto m12 = std::make_tuple( first_model( m1, m2 ) ); + auto all_models = std::tuple_cat( tuple, m12 ); + + auto type = particles.sliceType(); + using material_type = decltype( type ); + return ForceModels( type, + all_models ); +} + template struct BaseTemperatureModel { @@ -41,31 +183,25 @@ struct BaseTemperatureModel // Temperature field TemperatureType temperature; - BaseTemperatureModel(){}; BaseTemperatureModel( const TemperatureType _temp, const double _alpha, const double _temp0 ) : alpha( _alpha ) , temp0( _temp0 ) , temperature( _temp ){}; - BaseTemperatureModel( const BaseTemperatureModel& model ) - { - alpha = model.alpha; - temp0 = model.temp0; - temperature = model.temperature; - } - - void set_param( const double _alpha, const double _temp0 ) + // Average from existing models. + BaseTemperatureModel( const BaseTemperatureModel& model1, + const BaseTemperatureModel& model2 ) { - alpha = _alpha; - temp0 = _temp0; + alpha = ( model1.alpha + model2.alpha ) / 2.0; + temp0 = ( model1.temp0 + model2.temp0 ) / 2.0; } void update( const TemperatureType _temp ) { temperature = _temp; } // Update stretch with temperature effects. KOKKOS_INLINE_FUNCTION - void thermalStretch( double& s, const int i, const int j ) const + void thermalStretch( const int i, const int j, double& s ) const { double temp_avg = 0.5 * ( temperature( i ) + temperature( j ) ) - temp0; s -= alpha * temp_avg; @@ -87,12 +223,6 @@ struct BaseDynamicTemperatureModel BaseDynamicTemperatureModel( const double _delta, const double _kappa, const double _cp, const bool _constant_microconductivity = true ) - { - set_param( _delta, _kappa, _cp, _constant_microconductivity ); - } - - void set_param( const double _delta, const double _kappa, const double _cp, - const bool _constant_microconductivity ) { delta = _delta; kappa = _kappa; @@ -102,6 +232,16 @@ struct BaseDynamicTemperatureModel constant_microconductivity = _constant_microconductivity; } + // Average from existing models. + BaseDynamicTemperatureModel( const BaseDynamicTemperatureModel& model1, + const BaseDynamicTemperatureModel& model2 ) + { + delta = ( model1.delta + model2.delta ) / 2.0; + kappa = ( model1.kappa + model2.kappa ) / 2.0; + cp = ( model1.cp + model2.cp ) / 2.0; + constant_microconductivity = model1.constant_microconductivity; + } + KOKKOS_INLINE_FUNCTION double microconductivity_function( double r ) const { if ( constant_microconductivity ) diff --git a/src/CabanaPD_HeatTransfer.hpp b/src/CabanaPD_HeatTransfer.hpp index f34b0bf9..15c75527 100644 --- a/src/CabanaPD_HeatTransfer.hpp +++ b/src/CabanaPD_HeatTransfer.hpp @@ -20,10 +20,10 @@ namespace CabanaPD // Peridynamic heat transfer with forward-Euler time integration. // Inherits only because this is a similar neighbor-based kernel. template -class HeatTransfer : public Force +class HeatTransfer : public BaseForce { protected: - using base_type = Force; + using base_type = BaseForce; using base_type::_half_neigh; using base_type::_timer; diff --git a/src/CabanaPD_Input.hpp b/src/CabanaPD_Input.hpp index b91fe698..f4efbbcd 100644 --- a/src/CabanaPD_Input.hpp +++ b/src/CabanaPD_Input.hpp @@ -185,25 +185,45 @@ class Inputs double sum_ht = 0; // Estimate the bulk modulus if needed. - double K; + std::size_t min_index; + std::vector K; if ( inputs.contains( "bulk_modulus" ) ) { - K = inputs["bulk_modulus"]["value"]; + // Necessary because JSON can't determine the type otherwise. + std::vector K_vec = inputs["bulk_modulus"]["value"]; + K = K_vec; } else { - double E = inputs["elastic_modulus"]["value"]; + std::vector E = inputs["elastic_modulus"]["value"]; // This is only exact for bond-based (PMB). double nu = 0.25; - K = E / ( 3 * ( 1 - 2 * nu ) ); + + // Resize K to match the size of E + K.resize( E.size() ); + + for ( std::size_t i = 0; i < E.size(); i++ ) + K[i] = E[i] / ( 3.0 * ( 1.0 - 2.0 * nu ) ); } + // Support for multi-material: find the minimum density-bulk modulus + // ratio for the correct critical timestep. + std::vector rho_over_K; + std::vector rho = inputs["density"]["value"]; + for ( std::size_t i = 0; i < rho.size(); i++ ) + rho_over_K[i] = rho[i] / K[i]; + min_index = std::distance( std::begin( rho_over_K ), + std::min_element( std::begin( rho_over_K ), + std::end( rho_over_K ) ) ); + auto min_K = K[min_index]; + auto min_rho = rho[min_index]; + // Run over the neighborhood of a point in the bulk of a body (at the // origin). int m = inputs["m"]["value"]; double delta = inputs["horizon"]["value"]; // FIXME: this is copied from the forces - double c = 18.0 * K / ( pi * delta * delta * delta * delta ); + double c = 18.0 * min_K / ( pi * delta * delta * delta * delta ); for ( int i = -( m + 1 ); i < m + 2; i++ ) { @@ -250,8 +270,8 @@ class Inputs } double dt = inputs["timestep"]["value"]; - double rho = inputs["density"]["value"]; - double dt_crit = std::sqrt( 2.0 * rho / sum ); + // This supports multi-material. + double dt_crit = std::sqrt( 2.0 * min_rho / sum ); compareCriticalTimeStep( "mechanics", dt, dt_crit ); // Heat transfer timestep. @@ -262,7 +282,7 @@ class Inputs dt_ht *= dt; double cp = inputs["specific_heat_capacity"]["value"]; - double dt_ht_crit = rho * cp / sum_ht; + double dt_ht_crit = min_rho * cp / sum_ht; compareCriticalTimeStep( "heat_transfer", dt_ht, dt_ht_crit ); } } diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 08d69a7b..a63f54a4 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -105,10 +105,10 @@ class Particles using vector_type = Cabana::MemberTypes; // volume, dilatation, weighted_volume. using scalar_type = Cabana::MemberTypes; - // no-fail. + // no-fail, type. using int_type = Cabana::MemberTypes; - // v, rho, type. - using other_types = Cabana::MemberTypes; + // v, rho. + using other_types = Cabana::MemberTypes; // FIXME: add vector length. // FIXME: enable variable aosoa. @@ -116,6 +116,7 @@ class Particles using aosoa_y_type = Cabana::AoSoA; using aosoa_vol_type = Cabana::AoSoA; using aosoa_nofail_type = Cabana::AoSoA; + using aosoa_material_type = Cabana::AoSoA; using aosoa_other_type = Cabana::AoSoA; // Using grid here for the particle init. using plist_x_type = @@ -461,6 +462,11 @@ class Particles { return Cabana::slice<0>( _aosoa_vol, "volume" ); } + auto sliceType() { return Cabana::slice<0>( _aosoa_material, "type" ); } + auto sliceType() const + { + return Cabana::slice<0>( _aosoa_material, "type" ); + } auto sliceVelocity() { return Cabana::slice<0>( _aosoa_other, "velocities" ); @@ -474,8 +480,6 @@ class Particles { return Cabana::slice<1>( _aosoa_other, "density" ); } - auto sliceType() { return Cabana::slice<2>( _aosoa_other, "type" ); } - auto sliceType() const { return Cabana::slice<2>( _aosoa_other, "type" ); } auto sliceNoFail() { @@ -524,6 +528,7 @@ class Particles _plist_f.aosoa().resize( localOffset() ); _aosoa_other.resize( localOffset() ); _aosoa_nofail.resize( referenceOffset() ); + _aosoa_material.resize( referenceOffset() ); size = _plist_x.size(); _timer.stop(); }; @@ -549,7 +554,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(), sliceType(), std::forward( other )... ); #else #ifdef Cabana_ENABLE_SILO @@ -557,7 +562,7 @@ class Particles writePartialRangeTimeStep( "particles", local_grid->globalGrid(), output_step, output_time, 0, localOffset(), getPosition( use_reference ), sliceForce(), - sliceDisplacement(), sliceVelocity(), + sliceDisplacement(), sliceVelocity(), sliceType(), std::forward( other )... ); #else @@ -572,14 +577,17 @@ class Particles auto timeOutput() { return _output_timer.time(); }; auto time() { return _timer.time(); }; - friend class Comm; - friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; protected: aosoa_u_type _aosoa_u; aosoa_y_type _aosoa_y; aosoa_vol_type _aosoa_vol; aosoa_nofail_type _aosoa_nofail; + aosoa_material_type _aosoa_material; aosoa_other_type _aosoa_other; plist_x_type _plist_x; @@ -706,8 +714,10 @@ class Particles std::forward( other )... ); } - friend class Comm; - friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; protected: void init_lps() @@ -831,10 +841,10 @@ class Particles std::forward( other )... ); } - friend class Comm; - friend class Comm; - friend class Comm; - friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; protected: void init_temp() @@ -937,10 +947,14 @@ class Particles std::forward( other )... ); } - friend class Comm; - friend class Comm; - friend class Comm; - friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; protected: void init_output() diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index 6c091589..9e7834b6 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -103,16 +103,26 @@ class SolverNoFracture using particle_type = ParticleType; using integrator_type = Integrator; using force_model_type = ForceModelType; - using force_type = Force; - using comm_type = Comm; + using force_model_tag = typename force_model_type::model_type; + using force_fracture_type = typename force_model_type::fracture_type; + using force_type = Force; + using force_thermal_type = + typename force_model_type::thermal_type::base_type; + using comm_type = + Comm; using neigh_iter_tag = Cabana::SerialOpTag; using input_type = InputType; // Optional module types. using heat_transfer_type = HeatTransfer; - using contact_type = Force; using contact_model_type = ContactModelType; + using contact_model_tag = typename contact_model_type::model_type; + using contact_base_model_type = typename contact_model_type::base_model; + using contact_fracture_type = typename contact_model_type::fracture_type; + using contact_type = Force; SolverNoFracture( input_type _inputs, std::shared_ptr _particles, @@ -162,6 +172,10 @@ class SolverNoFracture "Contact with MPI is currently disabled." ); } + // Update optional property ghost sizes if needed. + if constexpr ( std::is_same::value ) + force_model.update( particles->sliceType() ); // Update temperature ghost size if needed. if constexpr ( is_temperature_dependent< typename force_model_type::thermal_type>::value ) @@ -238,7 +252,10 @@ class SolverNoFracture if ( !boundary_condition.forceUpdate() ) boundary_condition.apply( exec_space(), *particles, 0.0 ); - // Communicate temperature. + // Communicate optional properties. + if constexpr ( std::is_same::value ) + comm->gatherMaterial(); if constexpr ( is_temperature_dependent< typename force_model_type::thermal_type>::value ) comm->gatherTemperature(); @@ -270,6 +287,12 @@ class SolverNoFracture // Update ghost particles. comm->gatherDisplacement(); + // Communicate optional type. + if constexpr ( std::is_same< + typename force_model_type::material_type, + MultiMaterial>::value ) + comm->gatherMaterial(); + if constexpr ( is_heat_transfer< typename force_model_type::thermal_type>::value ) { @@ -322,6 +345,12 @@ class SolverNoFracture // Update ghost particles. comm->gatherDisplacement(); + // Communicate optional properties. + if constexpr ( std::is_same< + typename force_model_type::material_type, + MultiMaterial>::value ) + comm->gatherMaterial(); + // Compute internal forces. updateForce(); @@ -464,7 +493,6 @@ class SolverNoFracture // Optional modules. std::shared_ptr heat_transfer; std::shared_ptr contact; - contact_model_type contact_model; // Output files. std::string output_file; @@ -565,7 +593,10 @@ class SolverFracture if ( !boundary_condition.forceUpdate() ) boundary_condition.apply( exec_space(), *particles, 0.0 ); - // Communicate temperature. + // Communicate optional properties. + if constexpr ( std::is_same::value ) + comm->gatherMaterial(); if constexpr ( is_temperature_dependent< typename force_model_type::thermal_type>::value ) comm->gatherTemperature(); @@ -613,6 +644,12 @@ class SolverFracture // Update ghost particles. comm->gatherDisplacement(); + // Communicate optional type. + if constexpr ( std::is_same< + typename force_model_type::material_type, + MultiMaterial>::value ) + comm->gatherMaterial(); + // Compute internal forces. updateForce(); @@ -651,6 +688,11 @@ class SolverFracture // Update ghost particles. comm->gatherDisplacement(); + // Communicate optional properties. + if constexpr ( std::is_same< + typename force_model_type::material_type, + MultiMaterial>::value ) + comm->gatherMaterial(); // Compute internal forces. updateForce(); diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index 1fc303c7..f70a30e7 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -12,6 +12,8 @@ #ifndef TYPES_H #define TYPES_H +#include + namespace CabanaPD { // Fracture tags. @@ -27,12 +29,11 @@ struct Elastic { }; -// Contact and DEM (contact without PD) tags. -struct NoContact +// Material option tags. +struct SingleMaterial { }; -template -struct is_contact : public std::false_type +struct MultiMaterial { }; @@ -73,6 +74,19 @@ struct is_heat_transfer : public std::true_type { }; +// Contact and DEM (contact without PD) tags. +struct NoContact +{ + using base_model = std::false_type; + using model_type = std::false_type; + using thermal_type = TemperatureIndependent; + using fracture_type = NoFracture; +}; +template +struct is_contact : public std::false_type +{ +}; + // Force model tags. struct PMB { diff --git a/src/force/CabanaPD_ContactModels.hpp b/src/force/CabanaPD_ContactModels.hpp index 3afe3ae3..d6e9426d 100644 --- a/src/force/CabanaPD_ContactModels.hpp +++ b/src/force/CabanaPD_ContactModels.hpp @@ -25,10 +25,11 @@ namespace CabanaPD ******************************************************************************/ struct ContactModel { + using material_type = SingleMaterial; + double delta; double Rc; - ContactModel(){}; // PD horizon // Contact radius ContactModel( const double _delta, const double _Rc ) @@ -41,6 +42,7 @@ struct ContactModel struct NormalRepulsionModel : public ContactModel { // FIXME: This is for use as the primary force model. + using model_type = NormalRepulsionModel; using base_model = PMB; using fracture_type = NoFracture; using thermal_type = TemperatureIndependent; @@ -51,22 +53,23 @@ struct NormalRepulsionModel : public ContactModel double c; double K; - NormalRepulsionModel(){}; NormalRepulsionModel( const double delta, const double Rc, const double _K ) : ContactModel( delta, Rc ) , K( _K ) { - set_param( delta, Rc, K ); - } - - void set_param( const double _delta, const double _Rc, const double _K ) - { - delta = _delta; - Rc = _Rc; K = _K; // This could inherit from PMB (same c) c = 18.0 * K / ( pi * delta * delta * delta * delta ); } + + KOKKOS_INLINE_FUNCTION + auto forceCoeff( const double r, const double vol ) const + { + // Contact "stretch" + const double sc = ( r - Rc ) / delta; + // Normal repulsion uses a 15 factor compared to the PMB force + return 15.0 * c * sc * vol; + } }; template <> diff --git a/src/force/CabanaPD_ForceModels_LPS.hpp b/src/force/CabanaPD_ForceModels_LPS.hpp index 2317b267..f7966227 100644 --- a/src/force/CabanaPD_ForceModels_LPS.hpp +++ b/src/force/CabanaPD_ForceModels_LPS.hpp @@ -23,9 +23,11 @@ template <> struct ForceModel : public BaseForceModel { using base_type = BaseForceModel; + using model_type = LPS; using base_model = LPS; using fracture_type = NoFracture; using thermal_type = TemperatureIndependent; + using material_type = typename base_type::material_type; using base_type::delta; @@ -36,31 +38,64 @@ struct ForceModel : public BaseForceModel double theta_coeff; double s_coeff; + double c; + ForceModel( const double _delta, const double _K, const double _G, const int _influence = 0 ) : base_type( _delta ) , influence_type( _influence ) + , K( _K ) + , G( _G ) { - set_param( _delta, _K, _G ); - } - - void set_param( const double _delta, const double _K, const double _G ) - { - delta = _delta; - K = _K; - G = _G; - theta_coeff = 3.0 * K - 5.0 * G; s_coeff = 15.0 * G; + + // Equivalent PMB inputs. + c = 18.0 * K / ( pi * delta * delta * delta * delta ); } - KOKKOS_INLINE_FUNCTION double influence_function( double xi ) const + KOKKOS_INLINE_FUNCTION double influenceFunction( double xi ) const { if ( influence_type == 1 ) return 1.0 / xi; else return 1.0; } + + KOKKOS_INLINE_FUNCTION auto weightedVolume( const double xi, + const double vol ) const + { + return influenceFunction( xi ) * xi * xi * vol; + } + + KOKKOS_INLINE_FUNCTION auto dilatation( const double s, const double xi, + const double vol, + const double m_i ) const + { + double theta_i = influenceFunction( xi ) * s * xi * xi * vol; + return 3.0 * theta_i / m_i; + } + + KOKKOS_INLINE_FUNCTION auto + forceCoeff( const int, const int, const double s, const double xi, + const double vol, const double m_i, const double m_j, + const double theta_i, const double theta_j ) const + { + return ( theta_coeff * ( theta_i / m_i + theta_j / m_j ) + + s_coeff * s * ( 1.0 / m_i + 1.0 / m_j ) ) * + influenceFunction( xi ) * xi * vol; + } + + KOKKOS_INLINE_FUNCTION + auto energy( const int, const int, const double s, const double xi, + const double vol, const double m_i, const double theta_i, + const double num_bonds ) const + { + return 1.0 / num_bonds * 0.5 * theta_coeff / 3.0 * + ( theta_i * theta_i ) + + 0.5 * ( s_coeff / m_i ) * influenceFunction( xi ) * s * s * xi * + xi * vol; + } }; template <> @@ -68,6 +103,7 @@ struct ForceModel : public ForceModel { using base_type = ForceModel; + using model_type = LPS; using base_model = typename base_type::base_model; using fracture_type = Fracture; using thermal_type = base_type::thermal_type; @@ -82,18 +118,13 @@ struct ForceModel double s0; double bond_break_coeff; + using base_type::c; + ForceModel( const double _delta, const double _K, const double _G, const double _G0, const int _influence = 0 ) : base_type( _delta, _K, _G, _influence ) + , G0( _G0 ) { - set_param( _delta, _K, _G, _G0 ); - } - - void set_param( const double _delta, const double _K, const double _G, - const double _G0 ) - { - base_type::set_param( _delta, _K, _G ); - G0 = _G0; if ( influence_type == 1 ) { s0 = Kokkos::sqrt( 5.0 * G0 / 9.0 / K / delta ); // 1/xi @@ -111,6 +142,7 @@ struct ForceModel : public ForceModel { using base_type = ForceModel; + using model_type = LinearLPS; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using thermal_type = base_type::thermal_type; @@ -130,6 +162,7 @@ struct ForceModel : public ForceModel { using base_type = ForceModel; + using model_type = LinearLPS; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using thermal_type = base_type::thermal_type; diff --git a/src/force/CabanaPD_ForceModels_PMB.hpp b/src/force/CabanaPD_ForceModels_PMB.hpp index 3b6180b8..d7cbc760 100644 --- a/src/force/CabanaPD_ForceModels_PMB.hpp +++ b/src/force/CabanaPD_ForceModels_PMB.hpp @@ -26,32 +26,46 @@ struct ForceModel { using base_type = BaseForceModel; using base_model = PMB; + using model_type = PMB; using fracture_type = NoFracture; using thermal_type = TemperatureIndependent; + using material_type = typename base_type::material_type; using base_type::delta; double c; double K; - ForceModel( const double delta, const double K ) + ForceModel( const double delta, const double _K ) : base_type( delta ) + , K( _K ) { - set_param( delta, K ); + c = 18.0 * K / ( pi * delta * delta * delta * delta ); } - ForceModel( const ForceModel& model ) - : base_type( model ) + // Average from existing models. + template + ForceModel( const ModelType1& model1, const ModelType2& model2 ) + : base_type( model1.delta ) { - c = model.c; - K = model.K; + K = ( model1.K + model2.K ) / 2.0; + c = ( model1.c + model2.c ) / 2.0; } - void set_param( const double _delta, const double _K ) + KOKKOS_INLINE_FUNCTION + auto forceCoeff( const int, const int, const double s, + const double vol ) const { - delta = _delta; - K = _K; - c = 18.0 * K / ( pi * delta * delta * delta * delta ); + return c * s * vol; + } + + KOKKOS_INLINE_FUNCTION + auto energy( const int, const int, const double s, const double xi, + const double vol ) const + { + // 0.25 factor is due to 1/2 from outside the integral and 1/2 from + // the integrand (pairwise potential). + return 0.25 * c * s * s * xi * vol; } }; @@ -60,6 +74,7 @@ struct ForceModel : public ForceModel { using base_type = ForceModel; + using model_type = PMB; using base_model = typename base_type::base_model; using fracture_type = Fracture; using thermal_type = base_type::thermal_type; @@ -67,29 +82,28 @@ struct ForceModel using base_type::c; using base_type::delta; using base_type::K; + double G0; double s0; double bond_break_coeff; - ForceModel( const double delta, const double K, const double G0 ) - : base_type( delta, K ) - { - set_param( delta, K, G0 ); - } - - ForceModel( const ForceModel& model ) - : base_type( model ) + ForceModel( const double delta, const double _K, const double _G0 ) + : base_type( delta, _K ) + , G0( _G0 ) { - G0 = model.G0; - s0 = model.s0; - bond_break_coeff = model.bond_break_coeff; + s0 = Kokkos::sqrt( 5.0 * G0 / 9.0 / K / delta ); + bond_break_coeff = ( 1.0 + s0 ) * ( 1.0 + s0 ); } - void set_param( const double _delta, const double _K, const double _G0 ) + // Average from existing models. + template + ForceModel( const ModelType1& model1, const ModelType2& model2 ) + : base_type( model1, model2 ) { - base_type::set_param( _delta, _K ); - G0 = _G0; - s0 = Kokkos::sqrt( 5.0 * G0 / 9.0 / K / delta ); + G0 = ( model1.G0 + model2.G0 ) / 2.0; + s0 = Kokkos::sqrt( ( model1.s0 * model1.s0 * model1.c + + model2.s0 * model2.s0 * model2.c ) / + ( model1.c + model2.c ) ); bond_break_coeff = ( 1.0 + s0 ) * ( 1.0 + s0 ); } @@ -106,9 +120,10 @@ 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 thermal_type = base_type::thermal_type; + using thermal_type = typename base_type::thermal_type; using base_type::base_type; @@ -122,9 +137,10 @@ 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 thermal_type = base_type::thermal_type; + using thermal_type = typename base_type::thermal_type; using base_type::base_type; @@ -146,6 +162,7 @@ struct ForceModel; using base_temperature_type = BaseTemperatureModel; + using model_type = PMB; using base_model = PMB; using fracture_type = NoFracture; using thermal_type = TemperatureDependent; @@ -167,14 +184,6 @@ struct ForceModel using base_type = ForceModel; using base_temperature_type = BaseTemperatureModel; + using model_type = PMB; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using thermal_type = TemperatureDependent; @@ -232,14 +242,13 @@ struct ForceModel : base_type( _delta, _K, _G0 ) , base_temperature_type( _temp, _alpha, _temp0 ) { - set_param( _delta, _K, _G0, _alpha, _temp0 ); } - void set_param( const double _delta, const double _K, const double _G0, - const double _alpha, const double _temp0 ) + // Average from existing models. + ForceModel( const ForceModel& model1, const ForceModel& model2 ) + : base_type( model1, model2 ) + , base_temperature_type( model1, model2 ) { - base_type::set_param( _delta, _K, _G0 ); - base_temperature_type::set_param( _alpha, _temp0 ); } KOKKOS_INLINE_FUNCTION @@ -273,6 +282,7 @@ struct ForceModel using base_type = ForceModel; using base_temperature_type = BaseDynamicTemperatureModel; + using model_type = PMB; using base_model = PMB; using fracture_type = NoFracture; using thermal_type = DynamicTemperature; @@ -301,17 +311,6 @@ struct ForceModel , base_temperature_type( _delta, _kappa, _cp, _constant_microconductivity ) { - set_param( _delta, _K, _kappa, _cp, _alpha, _temp0, - _constant_microconductivity ); - } - - void set_param( const double _delta, const double _K, const double _kappa, - const double _cp, const double _alpha, const double _temp0, - const bool _constant_microconductivity ) - { - base_type::set_param( _delta, _K, _alpha, _temp0 ); - base_temperature_type::set_param( _delta, _kappa, _cp, - _constant_microconductivity ); } }; @@ -335,6 +334,7 @@ struct ForceModel using base_type = ForceModel; using base_temperature_type = BaseDynamicTemperatureModel; + using model_type = PMB; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using thermal_type = DynamicTemperature; @@ -364,20 +364,9 @@ struct ForceModel const double _temp0 = 0.0, const bool _constant_microconductivity = true ) : base_type( _delta, _K, _G0, _temp, _alpha, _temp0 ) + , base_temperature_type( _delta, _kappa, _cp, + _constant_microconductivity ) { - set_param( _delta, _K, _G0, _kappa, _cp, _alpha, _temp0 ); - base_temperature_type::set_param( _delta, _kappa, _cp, - _constant_microconductivity ); - } - - void set_param( const double _delta, const double _K, const double _G0, - const double _kappa, const double _cp, const double _alpha, - const double _temp0, - const bool _constant_microconductivity ) - { - base_type::set_param( _delta, _K, _G0, _alpha, _temp0 ); - base_temperature_type::set_param( _delta, _kappa, _cp, - _constant_microconductivity ); } }; diff --git a/src/force/CabanaPD_Force_Contact.hpp b/src/force/CabanaPD_Force_Contact.hpp index ba143635..828b8366 100644 --- a/src/force/CabanaPD_Force_Contact.hpp +++ b/src/force/CabanaPD_Force_Contact.hpp @@ -41,11 +41,11 @@ KOKKOS_INLINE_FUNCTION void getRelativeNormalVelocityComponents( Normal repulsion forces ******************************************************************************/ template -class Force - : public Force +class Force + : public BaseForce { public: - using base_type = Force; + using base_type = BaseForce; using neighbor_list_type = typename base_type::neighbor_list_type; template @@ -69,16 +69,15 @@ class Force const ParticleType& particles, ParallelType& neigh_op_tag ) { - auto delta = _model.delta; - auto Rc = _model.Rc; - auto c = _model.c; + 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(); _neigh_timer.start(); - _neigh_list.build( y, n_frozen, n_local, Rc, 1.0, mesh_min, mesh_max ); + _neigh_list.build( y, n_frozen, n_local, model.Rc, 1.0, mesh_min, + mesh_max ); _neigh_timer.stop(); auto contact_full = KOKKOS_LAMBDA( const int i, const int j ) @@ -91,11 +90,7 @@ class Force double rx, ry, rz; getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); - // Contact "stretch" - const double sc = ( r - Rc ) / delta; - - // Normal repulsion uses a 15 factor compared to the PMB force - const double coeff = 15 * c * sc * vol( j ); + const double coeff = model.forceCoeff( r, vol( j ) ); fcx_i = coeff * rx / r; fcy_i = coeff * ry / r; fcz_i = coeff * rz / r; diff --git a/src/force/CabanaPD_Force_HertzianContact.hpp b/src/force/CabanaPD_Force_HertzianContact.hpp index 688636b1..7b9ef12e 100644 --- a/src/force/CabanaPD_Force_HertzianContact.hpp +++ b/src/force/CabanaPD_Force_HertzianContact.hpp @@ -26,11 +26,11 @@ namespace CabanaPD Normal repulsion forces ******************************************************************************/ template -class Force - : public Force +class Force + : public BaseForce { public: - using base_type = Force; + using base_type = BaseForce; using neighbor_list_type = typename base_type::neighbor_list_type; template @@ -54,75 +54,44 @@ class Force const ParticleType& particles, ParallelType& neigh_op_tag ) { - auto Rc = _model.Rc; - auto radius = _model.radius; - auto Es = _model.Es; - auto Rs = _model.Rs; - auto beta = _model.beta; const int n_frozen = particles.frozenOffset(); const int n_local = particles.localOffset(); - const double coeff_h_n = 4.0 / 3.0 * Es * Kokkos::sqrt( Rs ); - const double coeff_h_d = -2.0 * Kokkos::sqrt( 5.0 / 6.0 ) * beta; - + auto model = _model; const auto vol = particles.sliceVolume(); const auto rho = particles.sliceDensity(); const auto y = particles.sliceCurrentPosition(); const auto vel = particles.sliceVelocity(); _neigh_timer.start(); - _neigh_list.build( y, n_frozen, n_local, Rc, 1.0, mesh_min, mesh_max ); + _neigh_list.build( y, n_frozen, n_local, model.Rc, 1.0, mesh_min, + mesh_max ); _neigh_timer.stop(); auto contact_full = KOKKOS_LAMBDA( const int i, const int j ) { - double fcx_i = 0.0; - double fcy_i = 0.0; - double fcz_i = 0.0; - double xi, r, s; double rx, ry, rz; getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); - // Contact "overlap" - const double delta_n = ( r - 2.0 * radius ); - - // Hertz normal force coefficient - double coeff = 0.0; - double Sn = 0.0; - if ( delta_n < 0.0 ) - { - coeff = Kokkos::min( - 0.0, -coeff_h_n * - Kokkos::pow( Kokkos::abs( delta_n ), 3.0 / 2.0 ) ); - Sn = 2.0 * Es * Kokkos::sqrt( Rs * Kokkos::abs( delta_n ) ); - } + const double delta_n = ( r - 2.0 * model.radius ); - coeff /= vol( i ); - - fcx_i = coeff * rx / r; - fcy_i = coeff * ry / r; - fcz_i = coeff * rz / r; - - fc( i, 0 ) += fcx_i; - fc( i, 1 ) += fcy_i; - fc( i, 2 ) += fcz_i; + const double coeff_normal = + model.normalForceCoeff( delta_n, vol( i ) ); + fc( i, 0 ) += coeff_normal * rx / r; + fc( i, 1 ) += coeff_normal * ry / r; + fc( i, 2 ) += coeff_normal * rz / r; // Hertz normal force damping component double vx, vy, vz, vn; getRelativeNormalVelocityComponents( vel, i, j, rx, ry, rz, r, vx, vy, vz, vn ); - double ms = ( rho( i ) * vol( i ) ) / 2.0; - double fnd = coeff_h_d * Kokkos::sqrt( Sn * ms ) * vn / vol( i ); - - fcx_i = fnd * rx / r; - fcy_i = fnd * ry / r; - fcz_i = fnd * rz / r; - - fc( i, 0 ) += fcx_i; - fc( i, 1 ) += fcy_i; - fc( i, 2 ) += fcz_i; + const double coeff_damp = + model.dampingForceCoeff( delta_n, vn, vol( i ), rho( i ) ); + fc( i, 0 ) += coeff_damp * rx / r; + fc( i, 1 ) += coeff_damp * ry / r; + fc( i, 2 ) += coeff_damp * rz / r; }; _timer.start(); diff --git a/src/force/CabanaPD_Force_LPS.hpp b/src/force/CabanaPD_Force_LPS.hpp index 8262e625..09cee8ab 100644 --- a/src/force/CabanaPD_Force_LPS.hpp +++ b/src/force/CabanaPD_Force_LPS.hpp @@ -69,14 +69,14 @@ namespace CabanaPD { -template -class Force> - : public Force +template +class Force + : public BaseForce { protected: - using base_type = Force; + using base_type = BaseForce; using base_type::_half_neigh; - using model_type = ForceModel; + using model_type = ModelType; model_type _model; using base_type::_energy_timer; @@ -114,8 +114,7 @@ class Force> // Get the reference positions and displacements. double xi, r, s; getDistance( x, u, i, j, xi, r, s ); - double m_j = model.influence_function( xi ) * xi * xi * vol( j ); - m( i ) += m_j; + m( i ) += model.weightedVolume( xi, vol( j ) ); }; Kokkos::RangePolicy policy( particles.frozenOffset(), @@ -146,9 +145,7 @@ class Force> // Get the bond distance, displacement, and stretch. double xi, r, s; getDistance( x, u, i, j, xi, r, s ); - double theta_i = - model.influence_function( xi ) * s * xi * xi * vol( j ); - theta( i ) += 3.0 * theta_i / m( i ); + theta( i ) += model.dilatation( s, xi, vol( j ), m( i ) ); }; Kokkos::RangePolicy policy( particles.frozenOffset(), @@ -168,10 +165,7 @@ class Force> { _timer.start(); - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; auto model = _model; - const auto vol = particles.sliceVolume(); auto theta = particles.sliceDilatation(); auto m = particles.sliceWeightedVolume(); @@ -184,10 +178,9 @@ class Force> double xi, r, s; double rx, ry, rz; getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); - const double coeff = - ( theta_coeff * ( theta( i ) / m( i ) + theta( j ) / m( j ) ) + - s_coeff * s * ( 1.0 / m( i ) + 1.0 / m( j ) ) ) * - model.influence_function( xi ) * xi * vol( j ); + + const double coeff = model.forceCoeff( + i, j, s, xi, vol( j ), m( i ), m( j ), theta( i ), theta( j ) ); fx_i = coeff * rx / r; fy_i = coeff * ry / r; fz_i = coeff * rz / r; @@ -214,8 +207,6 @@ class Force> { _energy_timer.start(); - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; auto model = _model; auto neigh_list = _neigh_list; @@ -233,11 +224,8 @@ class Force> Cabana::NeighborList::numNeighbor( neigh_list, i ) ); - double w = ( 1.0 / num_neighbors ) * 0.5 * theta_coeff / 3.0 * - ( theta( i ) * theta( i ) ) + - 0.5 * ( s_coeff / m( i ) ) * - model.influence_function( xi ) * s * s * xi * xi * - vol( j ); + double w = model.energy( i, j, s, xi, vol( j ), m( i ), theta( i ), + num_neighbors ); W( i ) += w; Phi += w * vol( i ); }; @@ -254,19 +242,16 @@ class Force> _energy_timer.stop(); return strain_energy; } - - auto time() { return _timer.time(); }; - auto timeEnergy() { return _energy_timer.time(); }; }; -template -class Force> - : public Force> +template +class Force + : public Force { protected: - using base_type = Force>; + using base_type = Force; using base_type::_half_neigh; - using model_type = ForceModel; + using model_type = ModelType; model_type _model; using base_type::_energy_timer; @@ -314,9 +299,7 @@ class Force> double xi, r, s; getDistance( x, u, i, j, xi, r, s ); // mu is included to account for bond breaking. - double m_j = mu( i, n ) * model.influence_function( xi ) * xi * - xi * vol( j ); - m( i ) += m_j; + m( i ) += mu( i, n ) * model.weightedVolume( xi, vol( j ) ); } }; @@ -356,14 +339,14 @@ class Force> // Get the bond distance, displacement, and stretch. double xi, r, s; getDistance( x, u, i, j, xi, r, s ); - // mu is included to account for bond breaking. - double theta_i = mu( i, n ) * model.influence_function( xi ) * - s * xi * xi * vol( j ); + // 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 ) += 3.0 * theta_i / m( i ); + theta( i ) += mu( i, n ) * + model.dilatation( s, xi, vol( j ), m( i ) ); } }; @@ -384,8 +367,6 @@ class Force> _timer.start(); auto break_coeff = _model.bond_break_coeff; - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; auto model = _model; auto neigh_list = _neigh_list; @@ -425,10 +406,8 @@ class Force> else if ( mu( i, n ) > 0 ) { const double coeff = - ( theta_coeff * - ( theta( i ) / m( i ) + theta( j ) / m( j ) ) + - s_coeff * s * ( 1.0 / m( i ) + 1.0 / m( j ) ) ) * - model.influence_function( xi ) * xi * vol( j ); + model.forceCoeff( i, j, s, xi, vol( j ), m( i ), m( j ), + theta( i ), theta( j ) ); double muij = mu( i, n ); fx_i = muij * coeff * rx / r; fy_i = muij * coeff * ry / r; @@ -457,8 +436,6 @@ class Force> { _energy_timer.start(); - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; auto model = _model; auto neigh_list = _neigh_list; @@ -487,11 +464,8 @@ class Force> getDistance( x, u, i, j, xi, r, s ); double w = - mu( i, n ) * ( ( 1.0 / num_bonds ) * 0.5 * theta_coeff / - 3.0 * ( theta( i ) * theta( i ) ) + - 0.5 * ( s_coeff / m( i ) ) * - model.influence_function( xi ) * s * s * - xi * xi * vol( j ) ); + mu( i, n ) * model.energy( i, j, s, xi, vol( j ), m( i ), + theta( i ), num_bonds ); W( i ) += w; phi_i += mu( i, n ) * vol( j ); @@ -512,13 +486,13 @@ class Force> } }; -template -class Force> - : public Force> +template +class Force + : public Force { protected: - using base_type = Force>; - using model_type = ForceModel; + using base_type = Force; + using model_type = ModelType; model_type _model; using base_type::_energy_timer; @@ -546,12 +520,9 @@ class Force> { _timer.start(); - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; auto model = _model; - const auto vol = particles.sliceVolume(); - auto linear_theta = particles.sliceDilatation(); + auto theta = particles.sliceDilatation(); // Using weighted volume from base LPS class. auto m = particles.sliceWeightedVolume(); @@ -568,10 +539,8 @@ class Force> xi_y, xi_z ); const double coeff = - ( theta_coeff * ( linear_theta( i ) / m( i ) + - linear_theta( j ) / m( j ) ) + - s_coeff * linear_s * ( 1.0 / m( i ) + 1.0 / m( j ) ) ) * - model.influence_function( xi ) * xi * vol( j ); + model.forceCoeff( i, j, linear_s, xi, vol( j ), m( i ), m( j ), + theta( i ), theta( j ) ); fx_i = coeff * xi_x / xi; fy_i = coeff * xi_y / xi; fz_i = coeff * xi_z / xi; @@ -598,20 +567,17 @@ class Force> { _energy_timer.start(); - auto theta_coeff = _model.theta_coeff; - auto s_coeff = _model.s_coeff; auto model = _model; auto neigh_list = _neigh_list; const auto vol = particles.sliceVolume(); - const auto linear_theta = particles.sliceDilatation(); + const auto theta = particles.sliceDilatation(); auto m = particles.sliceWeightedVolume(); auto energy_full = KOKKOS_LAMBDA( const int i, const int j, double& Phi ) { // Do we need to recompute linear_theta_i? - double xi, linear_s; getLinearizedDistance( x, u, i, j, xi, linear_s ); @@ -619,11 +585,8 @@ class Force> Cabana::NeighborList::numNeighbor( neigh_list, i ) ); - double w = ( 1.0 / num_neighbors ) * 0.5 * theta_coeff / 3.0 * - ( linear_theta( i ) * linear_theta( i ) ) + - 0.5 * ( s_coeff / m( i ) ) * - model.influence_function( xi ) * linear_s * - linear_s * xi * xi * vol( j ); + double w = model.energy( i, j, linear_s, xi, vol( j ), m( i ), + theta( i ), num_neighbors ); W( i ) += w; Phi += w * vol( i ); }; diff --git a/src/force/CabanaPD_Force_PMB.hpp b/src/force/CabanaPD_Force_PMB.hpp index 8a099639..cb2e73b0 100644 --- a/src/force/CabanaPD_Force_PMB.hpp +++ b/src/force/CabanaPD_Force_PMB.hpp @@ -70,15 +70,15 @@ namespace CabanaPD { -template -class Force> - : public Force +template +class 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 model_type = ModelType; + using base_type = BaseForce; using neighbor_list_type = typename base_type::neighbor_list_type; using base_type::_neigh_list; @@ -119,9 +119,9 @@ class Force> double rx, ry, rz; getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); - model.thermalStretch( s, i, j ); + model.thermalStretch( i, j, s ); - const double coeff = model.c * s * vol( j ); + const double coeff = model.forceCoeff( i, j, s, vol( j ) ); fx_i = coeff * rx / r; fy_i = coeff * ry / r; fz_i = coeff * rz / r; @@ -158,11 +158,9 @@ class Force> double xi, r, s; getDistance( x, u, i, j, xi, r, s ); - model.thermalStretch( s, i, j ); + model.thermalStretch( i, j, s ); - // 0.25 factor is due to 1/2 from outside the integral and 1/2 from - // the integrand (pairwise potential). - double w = 0.25 * model.c * s * s * xi * vol( j ); + double w = model.energy( i, j, s, xi, vol( j ) ); W( i ) += w; Phi += w * vol( i ); }; @@ -180,20 +178,19 @@ class Force> } }; -template -class Force> - : public Force +template +class 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 model_type = ModelType; + using base_type = BaseForce; using neighbor_list_type = typename base_type::neighbor_list_type; using base_type::_neigh_list; protected: - using base_model_type = typename model_type::base_type; using base_type::_half_neigh; model_type _model; @@ -242,7 +239,7 @@ class Force> double rx, ry, rz; getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz ); - model.thermalStretch( s, i, j ); + model.thermalStretch( i, j, s ); // Break if beyond critical stretch unless in no-fail zone. if ( model.criticalStretch( i, j, r, xi ) && !nofail( i ) && @@ -253,7 +250,8 @@ class Force> // Else if statement is only for performance. else if ( mu( i, n ) > 0 ) { - const double coeff = model.c * s * vol( j ); + const double coeff = model.forceCoeff( i, j, s, vol( j ) ); + double muij = mu( i, n ); fx_i = muij * coeff * rx / r; fy_i = muij * coeff * ry / r; @@ -302,11 +300,9 @@ class Force> double xi, r, s; getDistance( x, u, i, j, xi, r, s ); - model.thermalStretch( s, i, j ); + model.thermalStretch( i, j, s ); - // 0.25 factor is due to 1/2 from outside the integral and 1/2 - // from the integrand (pairwise potential). - double w = mu( i, n ) * 0.25 * model.c * s * s * xi * vol( j ); + double w = mu( i, n ) * model.energy( i, j, s, xi, vol( j ) ); W( i ) += w; phi_i += mu( i, n ) * vol( j ); @@ -327,17 +323,15 @@ class Force> } }; -template -class Force> - : public Force +template +class 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 model_type = ModelType; + using base_type = BaseForce; using neighbor_list_type = typename base_type::neighbor_list_type; using base_type::_neigh_list; @@ -379,9 +373,9 @@ class Force 2*radius) @@ -36,15 +38,12 @@ struct HertzianModel : public ContactModel double e; // Coefficient of restitution double beta; // Damping coefficient + double coeff_h_n; + double coeff_h_d; + HertzianModel( const double _Rc, const double _radius, const double _nu, const double _E, const double _e ) : ContactModel( 1.0, _Rc ) - { - set_param( _radius, _nu, _E, _e ); - } - - void set_param( const double _radius, const double _nu, const double _E, - const double _e ) { nu = _nu; radius = _radius; @@ -54,6 +53,38 @@ struct HertzianModel : public ContactModel double ln_e = Kokkos::log( e ); beta = -ln_e / Kokkos::sqrt( Kokkos::pow( ln_e, 2.0 ) + Kokkos::pow( pi, 2.0 ) ); + + // Derived constants. + coeff_h_n = 4.0 / 3.0 * Es * Kokkos::sqrt( Rs ); + coeff_h_d = -2.0 * Kokkos::sqrt( 5.0 / 6.0 ) * beta; + } + + KOKKOS_INLINE_FUNCTION + auto normalForceCoeff( const double delta_n, const double vol ) const + { + // Hertz normal force coefficient + double coeff = 0.0; + if ( delta_n < 0.0 ) + { + coeff = Kokkos::min( + 0.0, + -coeff_h_n * Kokkos::pow( Kokkos::abs( delta_n ), 3.0 / 2.0 ) ); + } + + coeff /= vol; + return coeff; + } + + KOKKOS_INLINE_FUNCTION + auto dampingForceCoeff( const double delta_n, const double vn, + const double vol, const double rho ) const + { + double Sn = 0.0; + if ( delta_n < 0.0 ) + Sn = 2.0 * Es * Kokkos::sqrt( Rs * Kokkos::abs( delta_n ) ); + double ms = ( rho * vol ) / 2.0; + double coeff = coeff_h_d * Kokkos::sqrt( Sn * ms ) * vn / vol; + return coeff; } }; diff --git a/unit_test/tstComm.hpp b/unit_test/tstComm.hpp index 95760534..d1781eae 100644 --- a/unit_test/tstComm.hpp +++ b/unit_test/tstComm.hpp @@ -77,7 +77,7 @@ void testHalo() Cabana::deep_copy( rank_init_host, rank ); // A gather is performed on construction. - CabanaPD::Comm comm( particles ); diff --git a/unit_test/tstForce.hpp b/unit_test/tstForce.hpp index 036329e2..c2580471 100644 --- a/unit_test/tstForce.hpp +++ b/unit_test/tstForce.hpp @@ -174,7 +174,7 @@ double computeReferenceWeightedVolume( ModelType model, const int m, double xi = sqrt( xi_x * xi_x + xi_y * xi_y + xi_z * xi_z ); if ( xi > 0.0 && xi < model.delta + 1e-14 ) weighted_volume += - model.influence_function( xi ) * xi * xi * vol; + model.influenceFunction( xi ) * xi * xi * vol; } return weighted_volume; } @@ -198,8 +198,7 @@ double computeReferenceDilatation( ModelType model, const int m, if ( xi > 0.0 && xi < model.delta + 1e-14 ) theta += 3.0 / weighted_volume * - model.influence_function( xi ) * s0 * xi * xi * - vol; + model.influenceFunction( xi ) * s0 * xi * xi * vol; } return theta; } @@ -231,7 +230,7 @@ double computeReferenceDilatation( ModelType model, const int m, if ( xi > 0.0 && xi < model.delta + 1e-14 ) theta += 3.0 / weighted_volume * - model.influence_function( xi ) * s * xi * xi * vol; + model.influenceFunction( xi ) * s * xi * xi * vol; } return theta; } @@ -287,8 +286,8 @@ double computeReferenceStrainEnergyDensity( W += ( 1.0 / num_neighbors ) * 0.5 * model.theta_coeff / 3.0 * ( theta * theta ) + 0.5 * ( model.s_coeff / weighted_volume ) * - model.influence_function( xi ) * s0 * s0 * xi * - xi * vol; + model.influenceFunction( xi ) * s0 * s0 * xi * xi * + vol; } } return W; @@ -335,7 +334,7 @@ double computeReferenceStrainEnergyDensity( W += ( 1.0 / num_neighbors ) * 0.5 * model.theta_coeff / 3.0 * ( theta_i * theta_j ) + 0.5 * ( model.s_coeff / weighted_volume ) * - model.influence_function( xi ) * s * s * xi * xi * + model.influenceFunction( xi ) * s * s * xi * xi * vol; } } @@ -384,7 +383,7 @@ double computeReferenceForceX( model.s_coeff * s * ( 1.0 / weighted_volume + 1.0 / weighted_volume ) ) * - model.influence_function( xi ) * xi * vol * rx / r; + model.influenceFunction( xi ) * xi * vol * rx / r; } } return fx; @@ -728,7 +727,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( true, particles, model ); + CabanaPD::Force + force( true, particles, model ); auto x = particles.sliceReferencePosition(); auto f = particles.sliceForce();