Add interpolation to powder fill example

examples/dem/powder_fill.cpp

@@ -143,7 +143,25 @@ void powderSettlingExample( const std::string filename )
     // ====================================================
     //                   Simulation run
     // ====================================================
+    solver.updateNeighbors();
     solver.run( gravity );
+
+    // ====================================================
+    // Interpolate for consolidation
+    // ====================================================
+    // Do this first to get the coarser grid.
+    std::array<int, 3> coarse_cells{ num_cells[0] / 2, num_cells[1] / 2,
+                                     num_cells[2] / 2 };
+    CabanaPD::Particles pd_particles(
+        memory_space{}, model_type{}, CabanaPD::BaseOutput{}, low_corner,
+        high_corner, coarse_cells, halo_width, create_powder, exec_space{} );
+    // FIXME: missing container.
+
+    CabanaPD::interpolate( solver.particles, pd_particles );
+
+    // Output
+    std::string name = "particles_pd";
+    pd_particles.output( name, 0, 0.0, true );
 }
 
 // Initialize MPI+Kokkos.

src/CabanaPD.hpp

@@ -22,6 +22,7 @@
 #include <CabanaPD_Geometry.hpp>
 #include <CabanaPD_Input.hpp>
 #include <CabanaPD_Integrate.hpp>
+#include <CabanaPD_Interpolation.hpp>
 #include <CabanaPD_Output.hpp>
 #include <CabanaPD_OutputProfiles.hpp>
 #include <CabanaPD_Particles.hpp>

src/CabanaPD_Force.hpp

@@ -194,7 +194,7 @@ class BaseForce
 
     // FIXME: should it be possible to update this list?
     template <class ParticleType>
-    void update( const ParticleType&, const double, const bool = false )
+    void update( const ParticleType&, const bool = false )
     {
     }
 

src/CabanaPD_Interpolation.hpp

@@ -0,0 +1,38 @@
+
+namespace CabanaPD
+{
+template <typename ParticleType, typename NewParticleType>
+void interpolate( const ParticleType& particles,
+                  NewParticleType& new_particles )
+{
+    using memory_space = typename NewParticleType::memory_space;
+
+    // Interpolate to the background grid.
+    auto powder_positions = particles.sliceCurrentPosition();
+    auto density = particles.sliceDensity();
+    auto num_particles = particles.localOffset();
+    auto scalar_p2g = Cabana::Grid::createScalarValueP2G( density, 0.125 );
+    auto scalar_layout = Cabana::Grid::createArrayLayout(
+        new_particles.local_grid, 1, Cabana::Grid::Node() );
+    auto scalar_grid_field = Cabana::Grid::createArray<double, memory_space>(
+        "scalar_grid_field", scalar_layout );
+    auto scalar_halo = Cabana::Grid::createHalo(
+        Cabana::Grid::NodeHaloPattern<3>(), new_particles.halo_width,
+        *scalar_grid_field );
+    Cabana::Grid::p2g( scalar_p2g, powder_positions, num_particles,
+                       Cabana::Grid::Spline<0>(), *scalar_halo,
+                       *scalar_grid_field );
+    Cabana::Grid::Experimental::BovWriter::writeTimeStep(
+        "grid_density", 0, 0.0, *scalar_grid_field );
+
+    // Now interpolate to what would be the consolidation particles.
+    auto consolidation_positions = new_particles.sliceCurrentPosition();
+    auto consolidation_density = new_particles.sliceDensity();
+    auto pd_num_particles = new_particles.localOffset();
+    auto scalar_value_g2p =
+        Cabana::Grid::createScalarValueG2P( consolidation_density, 1.0 );
+    Cabana::Grid::g2p( *scalar_grid_field, *scalar_halo,
+                       consolidation_positions, pd_num_particles,
+                       Cabana::Grid::Spline<0>(), scalar_value_g2p );
+}
+} // namespace CabanaPD

src/CabanaPD_Particles.hpp

@@ -479,7 +479,7 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
                     v( pid, d ) = 0.0;
                     f( pid, d ) = 0.0;
                 }
-                type( pid ) = 0;
+                type( pid ) = 1;
                 nofail( pid ) = 0;
                 rho( pid ) = 1.0;
             } );
@@ -735,14 +735,26 @@ class Particles<MemorySpace, PMB, TemperatureIndependent, BaseOutput, Dimension>
                  [[maybe_unused]] const double output_time,
                  [[maybe_unused]] const bool use_reference,
                  [[maybe_unused]] OtherFields&&... other )
+    {
+        output( "particles", output_step, output_time, use_reference,
+                other... );
+    }
+
+    // TODO: enable ignoring frozen particles.
+    template <typename... OtherFields>
+    void output( [[maybe_unused]] const std::string name,
+                 [[maybe_unused]] const int output_step,
+                 [[maybe_unused]] const double output_time,
+                 [[maybe_unused]] const bool use_reference,
+                 [[maybe_unused]] OtherFields&&... other )
     {
         _output_timer.start();
 
 #ifdef Cabana_ENABLE_HDF5
         Cabana::Experimental::HDF5ParticleOutput::writeTimeStep(
-            h5_config, "particles", MPI_COMM_WORLD, output_step, output_time,
+            h5_config, name, MPI_COMM_WORLD, output_step, output_time,
             localOffset(), getPosition( use_reference ), sliceForce(),
-            sliceDisplacement(), sliceVelocity(),
+            sliceDisplacement(), sliceVelocity(), sliceDensity(), sliceType(),
             std::forward<OtherFields>( other )... );
 #else
 #ifdef Cabana_ENABLE_SILO
@@ -880,6 +892,27 @@ class Particles<MemorySpace, LPS, TemperatureIndependent, BaseOutput, Dimension>
         _timer.stop();
     }
 
+    template <typename KeepType>
+    void remove( const int num_keep, const KeepType& keep )
+    {
+        base_type::remove( num_keep, keep );
+        _timer.start();
+        Cabana::remove( typename base_type::execution_space(), num_keep, keep,
+                        _aosoa_theta, base_type::numFrozen() );
+        Cabana::remove( typename base_type::execution_space(), num_keep, keep,
+                        _aosoa_m, base_type::numFrozen() );
+        _timer.stop();
+    }
+
+    template <typename... OtherFields>
+    void output( const std::string name, const int output_step,
+                 const double output_time, const bool use_reference,
+                 OtherFields&&... other )
+    {
+        base_type::output( name, output_step, output_time, use_reference,
+                           sliceWeightedVolume(), sliceDilatation(),
+                           std::forward<OtherFields>( other )... );
+    }
     template <typename... OtherFields>
     void output( const int output_step, const double output_time,
                  const bool use_reference, OtherFields&&... other )
@@ -998,6 +1031,25 @@ class Particles<MemorySpace, ModelType, TemperatureDependent, BaseOutput,
         _aosoa_temp.resize( base_type::referenceOffset() );
     }
 
+    template <typename KeepType>
+    void remove( const int num_keep, const KeepType& keep )
+    {
+        base_type::remove( num_keep, keep );
+        _timer.start();
+        Cabana::remove( typename base_type::execution_space(), num_keep, keep,
+                        _aosoa_temp, base_type::numFrozen() );
+        _timer.stop();
+    }
+
+    template <typename... OtherFields>
+    void output( const std::string name, const int output_step,
+                 const double output_time, const bool use_reference,
+                 OtherFields&&... other )
+    {
+        base_type::output( name, output_step, output_time, use_reference,
+                           sliceTemperature(),
+                           std::forward<OtherFields>( other )... );
+    }
     template <typename... OtherFields>
     void output( const int output_step, const double output_time,
                  const bool use_reference, OtherFields&&... other )
@@ -1020,6 +1072,7 @@ class Particles<MemorySpace, ModelType, TemperatureDependent, BaseOutput,
     }
 
     aosoa_temp_type _aosoa_temp;
+    using base_type::_timer;
 };
 
 template <class MemorySpace, class ThermalType, int Dimension>
@@ -1091,6 +1144,16 @@ class Particles<MemorySpace, Contact, ThermalType, BaseOutput, Dimension>
         _aosoa_u_neigh.resize( base_type::localOffset() );
     }
 
+    template <typename KeepType>
+    void remove( const int num_keep, const KeepType& keep )
+    {
+        base_type::remove( num_keep, keep );
+        _timer.start();
+        Cabana::remove( typename base_type::execution_space(), num_keep, keep,
+                        _aosoa_u_neigh, base_type::numFrozen() );
+        _timer.stop();
+    }
+
     void setMaxDisplacement( double new_max ) { _max_displacement = new_max; }
     double getMaxDisplacement() const { return _max_displacement; }
 
@@ -1110,6 +1173,7 @@ class Particles<MemorySpace, Contact, ThermalType, BaseOutput, Dimension>
     double _max_displacement;
 
     aosoa_u_neigh_type _aosoa_u_neigh;
+    using base_type::_timer;
 };
 
 template <class MemorySpace, class ModelType, class ThermalType, int Dimension>
@@ -1209,6 +1273,26 @@ class Particles<MemorySpace, ModelType, ThermalType, EnergyOutput, Dimension>
         _aosoa_output.resize( base_type::localOffset() );
     }
 
+    template <typename KeepType>
+    void remove( const int num_keep, const KeepType& keep )
+    {
+        base_type::remove( num_keep, keep );
+        _timer.start();
+        Cabana::remove( typename base_type::execution_space(), num_keep, keep,
+                        _aosoa_output, base_type::numFrozen() );
+        _timer.stop();
+    }
+
+    template <typename... OtherFields>
+    void output( const std::string name, const int output_step,
+                 const double output_time, const bool use_reference,
+                 OtherFields&&... other )
+    {
+        base_type::output( name, output_step, output_time, use_reference,
+                           sliceStrainEnergy(), sliceDamage(),
+                           std::forward<OtherFields>( other )... );
+    }
+
     template <typename... OtherFields>
     void output( const int output_step, const double output_time,
                  const bool use_reference, OtherFields&&... other )
@@ -1233,6 +1317,7 @@ class Particles<MemorySpace, ModelType, ThermalType, EnergyOutput, Dimension>
     }
 
     aosoa_output_type _aosoa_output;
+    using base_type::_timer;
 };
 
 template <class MemorySpace, class ModelType, class ThermalType, int Dimension>
@@ -1307,6 +1392,25 @@ class Particles<MemorySpace, ModelType, ThermalType, EnergyStressOutput,
         _aosoa_stress.resize( new_local + new_ghost );
     }
 
+    template <typename KeepType>
+    void remove( const int num_keep, const KeepType& keep )
+    {
+        base_type::remove( num_keep, keep );
+        _timer.start();
+        Cabana::remove( typename base_type::execution_space(), num_keep, keep,
+                        _aosoa_stress, base_type::numFrozen() );
+        _timer.stop();
+    }
+
+    template <typename... OtherFields>
+    void output( const std::string name, const int output_step,
+                 const double output_time, const bool use_reference,
+                 OtherFields&&... other )
+    {
+        base_type::output( name, output_step, output_time, use_reference,
+                           sliceStress(),
+                           std::forward<OtherFields>( other )... );
+    }
     template <typename... OtherFields>
     void output( const int output_step, const double output_time,
                  const bool use_reference, OtherFields&&... other )
@@ -1329,6 +1433,7 @@ class Particles<MemorySpace, ModelType, ThermalType, EnergyStressOutput,
     }
 
     aosoa_stress_type _aosoa_stress;
+    using base_type::_timer;
 };
 
 /******************************************************************************

src/CabanaPD_Solver.hpp

@@ -312,7 +312,7 @@ class Solver
         // FIXME: Will need to rebuild ghosts.
     }
 
-    void updateNeighbors() { force->update( particles, 0.0, true ); }
+    void updateNeighbors() { force->update( particles, true ); }
 
     template <typename BoundaryType>
     void runStep( const int step, BoundaryType boundary_condition )