Virial stress

.gitignore

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

.vscode/settings.json

@@ -0,0 +1,25 @@
+{
+    "files.associations": {
+        "array": "cpp",
+        "chrono": "cpp",
+        "cwchar": "cpp",
+        "forward_list": "cpp",
+        "vector": "cpp",
+        "exception": "cpp",
+        "fstream": "cpp",
+        "functional": "cpp",
+        "initializer_list": "cpp",
+        "iosfwd": "cpp",
+        "istream": "cpp",
+        "new": "cpp",
+        "ostream": "cpp",
+        "ratio": "cpp",
+        "sstream": "cpp",
+        "stdexcept": "cpp",
+        "streambuf": "cpp",
+        "tuple": "cpp",
+        "type_traits": "cpp",
+        "utility": "cpp",
+        "valarray": "cpp"
+    }
+}

src/CabanaPD_Particles.hpp

@@ -106,7 +106,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
     using int_type = Cabana::MemberTypes<int>;
     // v, W, rho, damage,  type.
     using other_types =
-        Cabana::MemberTypes<double[dim], double, double, double, int>;
+        Cabana::MemberTypes<double[dim * 2],double[dim], double, double, double, int>; //added double[dim * 2] for virial stress
     // Potentially needed later: body force (b), ID.
 
     // FIXME: add vector length.
@@ -344,6 +344,9 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, 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() ); } //added const version for slicing in energy calculation
+
     auto sliceForceAtomic()
     {
         auto f = sliceForce();
@@ -361,29 +364,38 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
     auto sliceType() const { return Cabana::slice<4>( _aosoa_other, "type" ); }
     auto sliceStrainEnergy()
     {
-        return Cabana::slice<1>( _aosoa_other, "strain_energy" );
+        return Cabana::slice<2>( _aosoa_other, "strain_energy" );
     }
     auto sliceStrainEnergy() const
     {
-        return Cabana::slice<1>( _aosoa_other, "strain_energy" );
+        return Cabana::slice<2>( _aosoa_other, "strain_energy" );
+    }
+    auto sliceVirialStress()
+    {
+        return Cabana::slice<0>(_aosoa_other, "virial_stress"); //added virial stress
     }
+    auto sliceVirialStress() const
+    {
+        return Cabana::slice<0>(_aosoa_other, "virial_stress"); //added virial stress const
+    }
+
     auto sliceVelocity()
     {
-        return Cabana::slice<0>( _aosoa_other, "velocities" );
+        return Cabana::slice<1>( _aosoa_other, "velocities" );
     }
     auto sliceVelocity() const
     {
-        return Cabana::slice<0>( _aosoa_other, "velocities" );
+        return Cabana::slice<1>( _aosoa_other, "velocities" );
     }
-    auto sliceDensity() { return Cabana::slice<2>( _aosoa_other, "density" ); }
+    auto sliceDensity() { return Cabana::slice<3>( _aosoa_other, "density" ); }
     auto sliceDensity() const
     {
-        return Cabana::slice<2>( _aosoa_other, "density" );
+        return Cabana::slice<3>( _aosoa_other, "density" );
     }
-    auto sliceDamage() { return Cabana::slice<3>( _aosoa_other, "damage" ); }
+    auto sliceDamage() { return Cabana::slice<4>( _aosoa_other, "damage" ); }
     auto sliceDamage() const
     {
-        return Cabana::slice<3>( _aosoa_other, "damage" );
+        return Cabana::slice<4>( _aosoa_other, "damage" );
     }
     auto sliceNoFail()
     {
@@ -450,7 +462,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
         Cabana::Experimental::HDF5ParticleOutput::writeTimeStep(
             h5_config, "particles", MPI_COMM_WORLD, output_step, output_time,
             n_local, getPosition( use_reference ), sliceStrainEnergy(),
-            sliceForce(), sliceDisplacement(), sliceVelocity(), sliceDamage() );
+            sliceForce(), sliceDisplacement(), sliceVelocity(), sliceDamage(), sliceVirialStress() );
 #else
 #ifdef Cabana_ENABLE_SILO
         Cabana::Grid::Experimental::SiloParticleOutput::

src/force/CabanaPD_Force_PMB.hpp

@@ -146,13 +146,17 @@ class Force<ExecutionSpace, ForceModel<PMB, Elastic, ModelParams...>>
 
         auto model = _model;
         const auto vol = particles.sliceVolume();
+        auto virial_stress = particles.sliceVirialStress(); // Moved outside of kernel
+        const auto f = particles.sliceForce(); //Slice forces outside of Kernel
+
 
         auto energy_full =
             KOKKOS_LAMBDA( const int i, const int j, double& Phi )
         {
             // Get the bond distance, displacement, and stretch.
             double xi, r, s;
-            getDistance( x, u, i, j, xi, r, s );
+            double rx, ry, rz;
+            getDistanceComponents( x, u, i, j, xi, r, s, rx, ry, rz );
 
             model.thermalStretch( s, i, j );
 
@@ -161,6 +165,19 @@ class Force<ExecutionSpace, ForceModel<PMB, Elastic, ModelParams...>>
             double w = 0.25 * model.c * s * s * xi * vol( j );
             W( i ) += w;
             Phi += w * vol( i );
+
+
+
+            // Use precomputed force for virial stress
+            virial_stress(i, 0) += (rx * f(i, 0)) / vol( j );  // σ_xx
+            virial_stress(i, 1) += (ry * f(i, 1)) / vol( j );  // σ_yy
+            virial_stress(i, 2) += (rz * f(i, 2)) / vol( j );  // σ_zz
+
+            virial_stress(i, 3) += (rx * f(i, 1)) / vol( j );  // σ_xy
+            virial_stress(i, 4) += (rx * f(i, 2)) / vol( j );  // σ_xz
+            virial_stress(i, 5) += (ry * f(i, 2)) / vol( j );  // σ_yz
+
+
         };
 
         double strain_energy = 0.0;