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