Virial stress

.gitignore

@@ -3,3 +3,5 @@
 build/
 *.out
 *.silo
+run/
+.vscode/

examples/mechanics/crack_branching.cpp

@@ -75,7 +75,7 @@ void crackBranchingExample( const std::string filename )
     // ====================================================
     // Note that individual inputs can be passed instead (see other examples).
     auto particles = CabanaPD::createParticles<memory_space, model_type>(
-        exec_space{}, inputs );
+        exec_space{}, inputs, CabanaPD::EnergyStressOutput{} );
 
     // ====================================================
     //                Boundary conditions planes

src/CabanaPD_Force.hpp

@@ -137,6 +137,7 @@ class Force<MemorySpace, BaseForceModel>
 
     Timer _timer;
     Timer _energy_timer;
+    Timer _stress_timer;
 
   public:
     using neighbor_list_type =
@@ -276,6 +277,22 @@ double computeEnergy( ForceType& force, ParticleType& particles,
     return energy;
 }
 
+template <class ForceType, class ParticleType, class ParallelType>
+void computeStress( ForceType& force, ParticleType& particles,
+                    const ParallelType& neigh_op_tag )
+{
+    if constexpr ( is_stress_output<typename ParticleType::output_type>::value )
+    {
+        auto stress = particles.sliceStress();
+
+        // Reset stress.
+        Cabana::deep_copy( stress, 0.0 );
+
+        force.computeStressFull( particles, neigh_op_tag );
+        Kokkos::fence();
+    }
+}
+
 // Forces with bond breaking.
 template <class ForceType, class ParticleType, class NeighborView,
           class ParallelType>
@@ -334,6 +351,22 @@ double computeEnergy( ForceType& force, ParticleType& particles,
     return energy;
 }
 
+template <class ForceType, class ParticleType, class NeighborView,
+          class ParallelType>
+void computeStress( ForceType& force, ParticleType& particles, NeighborView& mu,
+                    const ParallelType& neigh_op_tag )
+{
+    if constexpr ( is_stress_output<typename ParticleType::output_type>::value )
+    {
+        auto stress = particles.sliceStress();
+
+        // Reset stress.
+        Cabana::deep_copy( stress, 0.0 );
+
+        force.computeStressFull( particles, mu, neigh_op_tag );
+        Kokkos::fence();
+    }
+}
 } // namespace CabanaPD
 
 #endif

src/CabanaPD_Particles.hpp

@@ -346,6 +346,12 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
         return Cabana::slice<0>( _aosoa_u, "displacements" );
     }
     auto sliceForce() { return _plist_f.slice( CabanaPD::Field::Force() ); }
+
+    auto sliceForce() const
+    {
+        return _plist_f.slice( CabanaPD::Field::Force() );
+    }
+
     auto sliceForceAtomic()
     {
         auto f = sliceForce();
@@ -864,6 +870,112 @@ class Particles<MemorySpace, ModelType, ThermalType, EnergyOutput, Dimension>
     aosoa_output_type _aosoa_output;
 };
 
+template <class MemorySpace, class ModelType, class ThermalType, int Dimension>
+class Particles<MemorySpace, ModelType, ThermalType, EnergyStressOutput,
+                Dimension>
+    : public Particles<MemorySpace, ModelType, ThermalType, EnergyOutput,
+                       Dimension>
+{
+  public:
+    using self_type = Particles<MemorySpace, ModelType, ThermalType,
+                                EnergyStressOutput, Dimension>;
+    using base_type =
+        Particles<MemorySpace, ModelType, ThermalType, EnergyOutput, Dimension>;
+    using thermal_type = typename base_type::thermal_type;
+    using output_type = EnergyStressOutput;
+    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;
+
+    using output_types = Cabana::MemberTypes<double[dim][dim]>;
+    using aosoa_output_type = Cabana::AoSoA<output_types, memory_space, 1>;
+
+    // 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_output = aosoa_output_type( "Particle Stress Output", 0 );
+    }
+
+    // Constructor which initializes particles on regular grid.
+    template <class ExecSpace>
+    Particles( const ExecSpace& exec_space, std::array<double, dim> low_corner,
+               std::array<double, dim> high_corner,
+               const std::array<int, dim> num_cells, const int max_halo_width )
+        : base_type( exec_space, low_corner, high_corner, num_cells,
+                     max_halo_width )
+    {
+        _aosoa_output = aosoa_output_type( "Particle Stress Output", n_local );
+        init_output();
+    }
+
+    template <typename... Args>
+    void createParticles( Args&&... args )
+    {
+        // Forward arguments to standard or custom particle creation.
+        base_type::createParticles( std::forward<Args>( args )... );
+        _aosoa_output.resize( n_local );
+    }
+
+    auto sliceStress() { return Cabana::slice<0>( _aosoa_output, "stress" ); }
+    auto sliceStress() const
+    {
+        return Cabana::slice<0>( _aosoa_output, "stress" );
+    }
+
+    void resize( int new_local, int new_ghost )
+    {
+        base_type::resize( new_local, new_ghost );
+        _aosoa_output.resize( new_local + new_ghost );
+    }
+
+    template <typename... OtherFields>
+    void output( const int output_step, const double output_time,
+                 const bool use_reference, OtherFields&&... other )
+    {
+        base_type::output( output_step, output_time, use_reference,
+                           sliceStress(),
+                           std::forward<OtherFields>( other )... );
+    }
+
+    friend class Comm<self_type, PMB, TemperatureIndependent>;
+    friend class Comm<self_type, LPS, TemperatureIndependent>;
+    friend class Comm<self_type, PMB, TemperatureDependent>;
+    friend class Comm<self_type, LPS, TemperatureDependent>;
+
+  protected:
+    void init_output()
+    {
+        auto stress = sliceStress();
+        Cabana::deep_copy( stress, 0.0 );
+    }
+
+    aosoa_output_type _aosoa_output;
+};
+
 template <typename MemorySpace, typename ModelType, typename ExecSpace,
           typename OutputType>
 auto createParticles( ExecSpace exec_space, CabanaPD::Inputs inputs,

src/CabanaPD_Solver.hpp

@@ -191,6 +191,7 @@ class SolverElastic
         // Compute initial internal forces and energy.
         updateForce();
         computeEnergy( *force, *particles, neigh_iter_tag() );
+        computeStress( *force, *particles, neigh_iter_tag() );
 
         if ( initial_output )
             particles->output( 0, 0.0, output_reference );
@@ -320,6 +321,7 @@ class SolverElastic
         if ( step % output_frequency == 0 )
         {
             auto W = computeEnergy( *force, *particles, neigh_iter_tag() );
+            computeStress( *force, *particles, neigh_iter_tag() );
 
             particles->output( step / output_frequency, step * dt,
                                output_reference );
@@ -494,6 +496,7 @@ class SolverFracture
         // Compute initial internal forces and energy.
         updateForce();
         computeEnergy( *force, *particles, mu, neigh_iter_tag() );
+        computeStress( *force, *particles, mu, neigh_iter_tag() );
 
         if ( initial_output )
             particles->output( 0, 0.0, output_reference );
@@ -619,6 +622,7 @@ class SolverFracture
         if ( step % output_frequency == 0 )
         {
             auto W = computeEnergy( *force, *particles, mu, neigh_iter_tag() );
+            computeStress( *force, *particles, mu, neigh_iter_tag() );
 
             particles->output( step / output_frequency, step * dt,
                                output_reference );

src/CabanaPD_Types.hpp

@@ -78,6 +78,9 @@ struct BaseOutput
 struct EnergyOutput
 {
 };
+struct EnergyStressOutput
+{
+};
 
 template <class>
 struct is_energy_output : public std::false_type
@@ -87,6 +90,19 @@ template <>
 struct is_energy_output<EnergyOutput> : public std::true_type
 {
 };
+template <>
+struct is_energy_output<EnergyStressOutput> : public std::true_type
+{
+};
+
+template <class>
+struct is_stress_output : public std::false_type
+{
+};
+template <>
+struct is_stress_output<EnergyStressOutput> : public std::true_type
+{
+};
 
 } // namespace CabanaPD
 #endif