diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3f45e6d7..091f9bda 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -91,7 +91,10 @@ if(Picasso_ENABLE_TESTING)
 endif()
 
 # examples
-add_subdirectory(examples)
+option(Picasso_ENABLE_EXAMPLES "Build examples" OFF)
+if(Picasso_ENABLE_EXAMPLES)
+  add_subdirectory(examples)
+endif()
 
 # Package Configuration
 write_basic_package_version_file("PicassoConfigVersion.cmake"
@@ -110,7 +113,7 @@ install(FILES
 
 # Clang Format
 if(CLANG_FORMAT_FOUND)
-  file(GLOB_RECURSE FORMAT_SOURCES src/*.[c,h]pp unit_test/*.[c,h]pp)
+  file(GLOB_RECURSE FORMAT_SOURCES src/*.[c,h]pp unit_test/*.[c,h]pp examples/*.[c,h]pp)
   add_custom_target(picasso-format
     COMMAND ${CLANG_FORMAT_EXECUTABLE} -i -style=file ${FORMAT_SOURCES}
     DEPENDS ${FORMAT_SOURCES})
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index e69de29b..c06f315d 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -0,0 +1,10 @@
+configure_file(
+  ${CMAKE_CURRENT_SOURCE_DIR}/inputs/dam_break.json
+  ${CMAKE_CURRENT_BINARY_DIR}/dam_break.json
+  COPYONLY)
+
+add_executable( DamBreak dam_break.cpp )
+target_link_libraries( DamBreak PRIVATE Picasso::picasso )
+target_include_directories( DamBreak PRIVATE ${CMAKE_CURRENT_BINARY_DIR})
+
+install(TARGETS DamBreak DESTINATION ${CMAKE_INSTALL_BINDIR})
diff --git a/examples/dam_break.cpp b/examples/dam_break.cpp
new file mode 100644
index 00000000..a414c52c
--- /dev/null
+++ b/examples/dam_break.cpp
@@ -0,0 +1,457 @@
+#include <Picasso.hpp>
+
+#include <Cabana_Core.hpp>
+#include <Cabana_Grid.hpp>
+
+#include <Kokkos_Core.hpp>
+
+#include <mpi.h>
+
+using namespace Picasso;
+
+// Dam break particle creation and initialization.
+template <typename ParticleType, typename VelocityType>
+struct ParticleInitFunc
+{
+    Kokkos::Array<double, 6> block;
+    double density;
+
+    KOKKOS_INLINE_FUNCTION bool
+    operator()( const int, const double x[3], const double pv,
+                typename ParticleType::particle_type& p ) const
+    {
+        if ( block[0] <= x[0] && x[0] <= block[3] && block[1] <= x[1] &&
+             x[1] <= block[4] && block[2] <= x[2] && x[2] <= block[5] )
+        {
+
+            Picasso::get( p, Picasso::Field::Stress<3>() ) = 0.0;
+            Picasso::get( p, VelocityType() ) = 0.0;
+            Picasso::get( p, Picasso::Field::DetDefGrad() ) = 1.0;
+            Picasso::get( p, Picasso::Field::Mass() ) = pv * density;
+            Picasso::get( p, Picasso::Field::Volume() ) = pv;
+            Picasso::get( p, Picasso::Field::Pressure() ) = 0.0;
+
+            for ( int d = 0; d < 3; ++d )
+                Picasso::get( p, Picasso::Field::LogicalPosition<3>(), d ) =
+                    x[d];
+            return true;
+        }
+
+        return false;
+    }
+};
+
+// Dam break system boundaries.
+struct BoundaryCondition
+{
+    Kokkos::Array<long int, 6> bc_index_space;
+    Kokkos::Array<bool, 6> on_boundary;
+
+    // Free slip boundary condition
+    template <class ViewType>
+    KOKKOS_INLINE_FUNCTION void apply( ViewType view, const int i, const int j,
+                                       const int k ) const
+    {
+        // -x face
+        if ( i == bc_index_space[0] && on_boundary[0] )
+            view( i, j, k, 0 ) = 0.0;
+        // +x face
+        if ( i == bc_index_space[3] && on_boundary[3] )
+            view( i, j, k, 0 ) = 0.0;
+        // -y face
+        if ( j == bc_index_space[1] && on_boundary[1] )
+            view( i, j, k, 1 ) = 0.0;
+        // +y face
+        if ( j == bc_index_space[4] && on_boundary[4] )
+            view( i, j, k, 1 ) = 0.0;
+        // -z face
+        if ( k == bc_index_space[2] && on_boundary[2] )
+            view( i, j, k, 2 ) = 0.0;
+        // +z face
+        if ( k == bc_index_space[5] && on_boundary[5] )
+            view( i, j, k, 2 ) = 0.0;
+    }
+};
+
+// Custom field types for stress and gravity kernels.
+struct DeltaUGravity : Field::Vector<double, 3>
+{
+    static std::string label() { return "velocity_change_from_gravity"; }
+};
+struct DeltaUStress : Field::Vector<double, 3>
+{
+    static std::string label() { return "velocity_change_from_stress"; }
+};
+
+//---------------------------------------------------------------------------//
+// Update particle stress.
+//---------------------------------------------------------------------------//
+template <int InterpolationOrder>
+struct ComputeParticlePressure
+{
+    double dt;
+    double bulk_modulus;
+    double gamma;
+
+    template <class LocalMeshType, class GatherDependencies,
+              class ScatterDependencies, class LocalDependencies,
+              class ParticleViewType>
+    KOKKOS_INLINE_FUNCTION void
+    operator()( const LocalMeshType& local_mesh,
+                const GatherDependencies& gather_deps,
+                const ScatterDependencies&, const LocalDependencies&,
+                ParticleViewType& particle ) const
+    {
+        // Get particle data.
+        auto x_p =
+            Picasso::get( particle, Picasso::Field::LogicalPosition<3>() );
+        auto& J_p = Picasso::get( particle, Picasso::Field::DetDefGrad() );
+        auto& p_p = Picasso::get( particle, Picasso::Field::Pressure() );
+        auto s_p = Picasso::get( particle, Picasso::Field::Stress<3>() );
+        auto v_p = Picasso::get( particle, Picasso::Field::Volume() );
+        auto m_p = Picasso::get( particle, Picasso::Field::Mass() );
+
+        // Get the gather dependencies.
+        auto u_i = gather_deps.get( Picasso::FieldLocation::Node(),
+                                    Picasso::Field::Velocity<3>() );
+
+        // update strain rate
+        auto spline = Picasso::createSpline(
+            Picasso::FieldLocation::Node(),
+            Picasso::InterpolationOrder<InterpolationOrder>(), local_mesh, x_p,
+            Picasso::SplineValue(), Picasso::SplineGradient() );
+        Picasso::Mat3<double> vel_grad;
+        Picasso::G2P::gradient( spline, u_i, vel_grad );
+
+        // J_p = Kokkos::abs( !F_p );
+        J_p *= Kokkos::exp( Picasso::LinearAlgebra::trace( dt * vel_grad ) );
+
+        p_p = bulk_modulus * ( Kokkos::pow( J_p, -gamma ) - 1.0 );
+
+        Picasso::Mat3<double> I;
+        Picasso::LinearAlgebra::identity( I );
+
+        s_p = -p_p * I;
+    }
+};
+
+//---------------------------------------------------------------------------//
+// Grid momentum change due to stress
+//---------------------------------------------------------------------------//
+template <int InterpolationOrder>
+struct ComputeGridVelocityChange
+{
+    Picasso::Vec3<double> gravity;
+
+    ComputeGridVelocityChange( Kokkos::Array<double, 3> g )
+    {
+        for ( int d = 0; d < 3; ++d )
+            gravity( d ) = g[d];
+    }
+
+    template <class LocalMeshType, class GatherDependencies,
+              class ScatterDependencies, class LocalDependencies,
+              class ParticleViewType>
+    KOKKOS_INLINE_FUNCTION void
+    operator()( const LocalMeshType& local_mesh, const GatherDependencies&,
+                const ScatterDependencies& scatter_deps,
+                const LocalDependencies&, ParticleViewType& particle ) const
+    {
+        // Get particle data.
+        auto m_p = Picasso::get( particle, Picasso::Field::Mass() );
+        auto v_p = Picasso::get( particle, Picasso::Field::Volume() );
+        auto x_p =
+            Picasso::get( particle, Picasso::Field::LogicalPosition<3>() );
+        auto s_p = Picasso::get( particle, Picasso::Field::Stress<3>() );
+
+        // Get the scatter dependencies.
+        auto delta_u_s_i =
+            scatter_deps.get( Picasso::FieldLocation::Node(), DeltaUStress() );
+        auto delta_u_g_i =
+            scatter_deps.get( Picasso::FieldLocation::Node(), DeltaUGravity() );
+
+        // Node interpolant.
+        auto spline = Picasso::createSpline(
+            Picasso::FieldLocation::Node(),
+            Picasso::InterpolationOrder<InterpolationOrder>(), local_mesh, x_p,
+            Picasso::SplineValue(), Picasso::SplineGradient() );
+
+        // Compute velocity update from stress.
+        Picasso::P2G::divergence( spline, -v_p * s_p, delta_u_s_i );
+        // Compute velocity update from gravity.
+        Picasso::P2G::value( spline, m_p * gravity, delta_u_g_i );
+    }
+};
+
+//---------------------------------------------------------------------------//
+// Update nodal momentum n+1 with stress, gravity, and BC
+//---------------------------------------------------------------------------//
+struct UpdateGridVelocity
+{
+    // Explicit time step size.
+    double dt;
+    BoundaryCondition bc;
+
+    template <class LocalMeshType, class GatherDependencies,
+              class ScatterDependencies, class LocalDependencies>
+    KOKKOS_INLINE_FUNCTION void
+    operator()( const LocalMeshType&, const GatherDependencies& gather_deps,
+                const ScatterDependencies&, const LocalDependencies&,
+                const int i, const int j, const int k ) const
+    {
+        // Get the gather dependencies.
+        auto m_i = gather_deps.get( Picasso::FieldLocation::Node(),
+                                    Picasso::Field::Mass() );
+        auto u_i = gather_deps.get( Picasso::FieldLocation::Node(),
+                                    Picasso::Field::Velocity<3>() );
+        auto delta_u_s_i =
+            gather_deps.get( Picasso::FieldLocation::Node(), DeltaUStress() );
+
+        auto delta_u_g_i =
+            gather_deps.get( Picasso::FieldLocation::Node(), DeltaUGravity() );
+
+        // Compute velocity.
+        Picasso::Vec3<double> zeros = { 0.0, 0.0, 0.0 };
+        u_i( i, j, k ) +=
+            ( m_i( i, j, k ) > 0.0 )
+                ? dt * ( delta_u_s_i( i, j, k ) + delta_u_g_i( i, j, k ) ) /
+                      m_i( i, j, k )
+                : zeros;
+
+        // Add boundary conditions last.
+        bc.apply( u_i, i, j, k );
+    }
+};
+
+//---------------------------------------------------------------------------//
+// DamBreak example
+template <class InterpolationType, class ParticleVelocity>
+void DamBreak( std::string filename )
+{
+    using exec_space = Kokkos::DefaultExecutionSpace;
+    using memory_space = exec_space::memory_space;
+
+    // Get inputs for mesh.
+    auto inputs = Picasso::parse( filename );
+
+    // Global bounding box.
+    auto global_box = copy<double, 6>( inputs["global_box"] );
+    int minimum_halo_size = 0;
+
+    // Make mesh.
+    using mesh_type = Picasso::UniformMesh<memory_space>;
+    std::shared_ptr<mesh_type> mesh = Picasso::createUniformMesh(
+        memory_space(), inputs, global_box, minimum_halo_size, MPI_COMM_WORLD );
+
+    // Make a particle list.
+    Cabana::ParticleTraits<Picasso::Field::Stress<3>, ParticleVelocity,
+                           Picasso::Field::LogicalPosition<3>,
+                           Picasso::Field::Mass, Picasso::Field::Pressure,
+                           Picasso::Field::Volume, Picasso::Field::DetDefGrad>
+        fields;
+    auto particles = Cabana::Grid::createParticleList<memory_space>(
+        "test_particles", fields );
+
+    // Initialize particles
+    auto particle_box = copy<double, 6>( inputs["particle_box"] );
+    double density = inputs["density"];
+    ParticleInitFunc<decltype( particles ), ParticleVelocity>
+        momentum_init_functor{ particle_box, density };
+
+    double ppc = inputs["ppc"];
+    auto local_grid = mesh->localGrid();
+    Cabana::Grid::createParticles( Cabana::InitUniform(), exec_space{},
+                                   momentum_init_functor, particles, ppc,
+                                   *local_grid );
+
+    // Boundary index space
+    auto index_space = local_grid->indexSpace(
+        Cabana::Grid::Own(), Cabana::Grid::Node(), Cabana::Grid::Local() );
+    auto global_grid = local_grid->globalGrid();
+
+    Kokkos::Array<long int, 6> bc_index_space{
+        index_space.min( Cabana::Grid::Dim::I ),
+        index_space.min( Cabana::Grid::Dim::J ),
+        index_space.min( Cabana::Grid::Dim::K ),
+        index_space.max( Cabana::Grid::Dim::I ) - 1,
+        index_space.max( Cabana::Grid::Dim::J ) - 1,
+        index_space.max( Cabana::Grid::Dim::K ) - 1 };
+
+    Kokkos::Array<bool, 6> on_boundary{
+        global_grid.onLowBoundary( Cabana::Grid::Dim::I ),
+        global_grid.onLowBoundary( Cabana::Grid::Dim::J ),
+        global_grid.onLowBoundary( Cabana::Grid::Dim::K ),
+        global_grid.onHighBoundary( Cabana::Grid::Dim::I ),
+        global_grid.onHighBoundary( Cabana::Grid::Dim::J ),
+        global_grid.onHighBoundary( Cabana::Grid::Dim::K ) };
+
+    BoundaryCondition bc{ bc_index_space, on_boundary };
+
+    // Properties
+    double gamma = inputs["gamma"];
+    double bulk_modulus = inputs["bulk_modulus"];
+    auto gravity = copy<double, 3>( inputs["gravity"] );
+
+    // Time integragor inputs.
+    double dt = inputs["dt"];
+    // Only used for FLIP.
+    double beta = inputs["beta"];
+    const int spline_order = 1;
+
+    auto fm = Picasso::createFieldManager( mesh );
+
+    // Field locations and variables.
+    using node_type = Picasso::FieldLocation::Node;
+    using grid_type =
+        Picasso::FieldLayout<node_type, Picasso::Field::PhysicalPosition<3>>;
+    using particle_type = Picasso::FieldLocation::Particle;
+    using mass_type = Picasso::FieldLayout<node_type, Picasso::Field::Mass>;
+    using velocity_type =
+        Picasso::FieldLayout<node_type, Picasso::Field::Velocity<3>>;
+    using old_u_type =
+        Picasso::FieldLayout<node_type, Picasso::Field::OldVelocity<3>>;
+    using delta_u_s_type = Picasso::FieldLayout<node_type, DeltaUStress>;
+    using delta_u_g_type = Picasso::FieldLayout<node_type, DeltaUGravity>;
+
+    // Define particle/grid parallel data dependencies.
+    Picasso::GridOperator<mesh_type,
+                          Picasso::ScatterDependencies<mass_type, old_u_type>>
+        p2g( mesh );
+    p2g.setup( *fm );
+
+    Picasso::GridOperator<mesh_type, Picasso::LocalDependencies<
+                                         mass_type, velocity_type, old_u_type>>
+        compute_velocity( mesh );
+    compute_velocity.setup( *fm );
+
+    Picasso::GridOperator<mesh_type, Picasso::GatherDependencies<velocity_type>>
+        compute_stress( mesh );
+    compute_stress.setup( *fm );
+
+    Picasso::GridOperator<
+        mesh_type, Picasso::ScatterDependencies<delta_u_s_type, delta_u_g_type>>
+        compute_velocity_change( mesh );
+    compute_velocity_change.setup( *fm );
+
+    Picasso::GridOperator<
+        mesh_type, Picasso::GatherDependencies<mass_type, velocity_type,
+                                               delta_u_s_type, delta_u_g_type>>
+        update_velocity( mesh );
+    update_velocity.setup( *fm );
+
+    Picasso::GridOperator<
+        mesh_type,
+        Picasso::GatherDependencies<mass_type, velocity_type, old_u_type>,
+        Picasso::LocalDependencies<grid_type>>
+        g2p( mesh );
+    g2p.setup( *fm );
+
+    // Timestep loop.
+    double time = 0.0;
+    auto final_time = inputs["final_time"];
+    int write_frequency = inputs["write_frequency"];
+    int steps = 0;
+    while ( time < final_time )
+    {
+        // Particle interpolation (Picasso built-in).
+        auto p2g_func =
+            Picasso::createParticle2Grid<spline_order, InterpolationType,
+                                         ParticleVelocity>( dt );
+        p2g.apply( "p2g_u", particle_type{}, exec_space{}, *fm, particles,
+                   p2g_func );
+
+        // Compute grid velocity (Picasso built-in).
+        Picasso::ComputeGridVelocity compute_u_func;
+        compute_velocity.apply( "grid_U", node_type{}, exec_space{}, *fm,
+                                compute_u_func );
+
+        // Update particle stress.
+        ComputeParticlePressure<spline_order> compute_s_func{ dt, bulk_modulus,
+                                                              gamma };
+        compute_stress.apply( "particle_S", particle_type{}, exec_space{}, *fm,
+                              particles, compute_s_func );
+
+        // Compute grid velocity change due to stress
+        ComputeGridVelocityChange<spline_order> compute_du_s_func( gravity );
+        compute_velocity_change.apply( "div_S", particle_type{}, exec_space{},
+                                       *fm, particles, compute_du_s_func );
+
+        // Compute next grid velocity for gravity/stress.
+        UpdateGridVelocity update_u_func{ dt, bc };
+        update_velocity.apply( "update_U", node_type{}, exec_space{}, *fm,
+                               update_u_func );
+
+        // Grid interpolation (Picasso built-in).
+        auto g2p_func =
+            Picasso::createGrid2ParticleVelocity<spline_order,
+                                                 InterpolationType>( dt, beta );
+        g2p.apply( "g2p_U", particle_type{}, exec_space{}, *fm, particles,
+                   g2p_func );
+
+        // Redistribution of particles across MPI subdomains.
+        particles.redistribute( *local_grid,
+                                Picasso::Field::LogicalPosition<3>() );
+
+        // Write particle fields.
+        Cabana::Experimental::HDF5ParticleOutput::HDF5Config h5_config;
+        if ( steps % write_frequency == 0 )
+        {
+            Cabana::Experimental::HDF5ParticleOutput::writeTimeStep(
+                h5_config, "particles", global_grid.comm(),
+                steps / write_frequency, time, particles.size(),
+                particles.slice( Picasso::Field::LogicalPosition<3>() ),
+                particles.slice( Picasso::Field::Pressure() ),
+                particles.slice( ParticleVelocity() ),
+                particles.slice( Picasso::Field::Mass() ),
+                particles.slice( Picasso::Field::Volume() ) );
+
+            // Calculate conservation sums
+            double mass_particles = 0.0;
+            double ke_particles = 0.0;
+            double pe_particles = 0.0;
+            Picasso::particleConservation(
+                global_grid.comm(), exec_space(), particles, ParticleVelocity(),
+                mass_particles, ke_particles, pe_particles,
+                Cabana::Grid::Dim::K, global_box[2], gravity );
+            double mass_grid = 0.0;
+            double ke_grid = 0.0;
+            double pe_grid = 0.0;
+            Picasso::gridConservation(
+                global_grid.comm(), exec_space(), mesh, *fm, mass_grid, ke_grid,
+                pe_grid, Cabana::Grid::Dim::K, global_box[2], gravity );
+
+            if ( global_grid.blockId() == 0 )
+                std::cout << "Particle/Grid Mass: " << mass_particles << " / "
+                          << mass_grid << "\n"
+                          << "Particle/Grid Total Energy: "
+                          << ke_particles + pe_particles << " / "
+                          << ke_grid + pe_grid << "\n\n";
+        }
+        time += dt;
+        steps++;
+    }
+}
+
+int main( int argc, char* argv[] )
+{
+    MPI_Init( &argc, &argv );
+    Kokkos::initialize( argc, argv );
+
+    if ( argc < 2 )
+        throw std::runtime_error( "Incorrect number of arguments. \n \
+            First argument - file name for output \n \
+            \n \
+            Example: \n \
+            $/: ./DamBreak inputs/dam_break.json\n" );
+    std::string filename = argv[1];
+
+    // Problem can run with any interpolation scheme.
+    // DamBreak<PolyPicTag, Picasso::PolyPIC::Field::Velocity<3>>( filename);
+    // DamBreak<APicTag, Picasso::APIC::Field::Velocity<3>>( filename );
+    DamBreak<FlipTag, Picasso::Field::Velocity<3>>( filename );
+
+    Kokkos::finalize();
+    MPI_Finalize();
+
+    return 0;
+}
diff --git a/examples/inputs/dam_break.json b/examples/inputs/dam_break.json
new file mode 100644
index 00000000..31de2a87
--- /dev/null
+++ b/examples/inputs/dam_break.json
@@ -0,0 +1,21 @@
+{
+    "mesh": {
+        "cell_size": 0.05,
+        "periodic": [false, false, false],
+        "partitioner": {
+            "type": "uniform_dim"
+        },
+        "halo_cell_width": 0
+    },
+    "global_box": [ 0.0, 0.0, 0.0, 1.0, 1.0, 1.0 ],
+    "particle_box": [ 0.0, 0.0, 0.0, 0.4, 0.4, 0.4 ],
+    "density": 1.0e3,
+    "ppc": 2,
+    "gamma": 7.0,
+    "bulk_modulus": 1.0e5,
+    "gravity": [ 0.0, 0.0, -9.8],
+    "dt": 1.0e-4,
+    "beta": 0.99,
+    "final_time": 1.0,
+    "write_frequency" : 10
+}
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 34b8f23f..8e3866c0 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -12,7 +12,9 @@ set(HEADERS
   Picasso_FieldManager.hpp
   Picasso_FieldTypes.hpp
   Picasso_GridOperator.hpp
+  Picasso_GridKernels.hpp
   Picasso_InputParser.hpp
+  Picasso_InterpolationKernels.hpp
   Picasso_LevelSet.hpp
   Picasso_LevelSetRedistance.hpp
   Picasso_MarchingCubes.hpp
diff --git a/src/Picasso.hpp b/src/Picasso.hpp
index a5104c29..9a21c2bf 100644
--- a/src/Picasso.hpp
+++ b/src/Picasso.hpp
@@ -18,12 +18,15 @@
 #include <Picasso_AdaptiveMesh.hpp>
 #include <Picasso_BatchedLinearAlgebra.hpp>
 #include <Picasso_BilinearMeshMapping.hpp>
+#include <Picasso_Conservation.hpp>
 #include <Picasso_CurvilinearMesh.hpp>
 #include <Picasso_FacetGeometry.hpp>
 #include <Picasso_FieldManager.hpp>
 #include <Picasso_FieldTypes.hpp>
+#include <Picasso_GridKernels.hpp>
 #include <Picasso_GridOperator.hpp>
 #include <Picasso_InputParser.hpp>
+#include <Picasso_InterpolationKernels.hpp>
 #include <Picasso_LevelSet.hpp>
 #include <Picasso_LevelSetRedistance.hpp>
 #include <Picasso_ParticleInit.hpp>
diff --git a/src/Picasso_Conservation.hpp b/src/Picasso_Conservation.hpp
new file mode 100644
index 00000000..11c1bf6e
--- /dev/null
+++ b/src/Picasso_Conservation.hpp
@@ -0,0 +1,118 @@
+/****************************************************************************
+ * Copyright (c) 2021 by the Picasso authors                                *
+ * All rights reserved.                                                     *
+ *                                                                          *
+ * This file is part of the Picasso library. Picasso 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                                    *
+ ****************************************************************************/
+
+#ifndef PICASSO_CONSERVATION_HPP
+#define PICASSO_CONSERVATION_HPP
+
+#include <Picasso.hpp>
+
+#include <mpi.h>
+
+namespace Picasso
+{
+// Output routine for integral conservation sums of particle quantities
+// (momentum-only)
+template <class ExecutionSpace, class ParticleList, class ParticleVelocity>
+void particleConservation( MPI_Comm comm, ExecutionSpace,
+                           const ParticleList particles, ParticleVelocity,
+                           double& global_mass, double& global_ke,
+                           double& global_pe, const int pe_dir,
+                           const double zero_ground,
+                           const Kokkos::Array<double, 3> gravity )
+{
+    std::size_t num_p = particles.size();
+    auto aosoa = particles.aosoa();
+
+    double mass_integral = 0;
+    double ke_integral = 0;
+    double pe_integral = 0;
+    Kokkos::parallel_reduce(
+        "Picasso::Particle Conservation",
+        Kokkos::RangePolicy<ExecutionSpace>( 0, num_p ),
+        KOKKOS_LAMBDA( const int pn, double& mass_cont, double& ke_cont,
+                       double& pe_cont ) {
+            typename ParticleList::particle_type p( aosoa.getTuple( pn ) );
+            auto m_p = Picasso::get( p, Picasso::Field::Mass() );
+            auto v_p = Picasso::get( p, ParticleVelocity() );
+            auto x_p = Picasso::get( p, Picasso::Field::LogicalPosition<3>() );
+            auto u = Picasso::getField( v_p, 0 );
+            auto v = Picasso::getField( v_p, 1 );
+            auto w = Picasso::getField( v_p, 2 );
+            mass_cont += m_p;
+            double mag_v = sqrt( u * u + v * v + w * w );
+            ke_cont += m_p * mag_v * mag_v / 2;
+            pe_cont += m_p * Kokkos::abs( gravity[pe_dir] ) *
+                       ( x_p( pe_dir ) - zero_ground );
+        },
+        mass_integral, ke_integral, pe_integral );
+
+    MPI_Allreduce( &mass_integral, &global_mass, 1, MPI_DOUBLE, MPI_SUM, comm );
+    MPI_Allreduce( &ke_integral, &global_ke, 1, MPI_DOUBLE, MPI_SUM, comm );
+    MPI_Allreduce( &pe_integral, &global_pe, 1, MPI_DOUBLE, MPI_SUM, comm );
+}
+
+// Output for integral conservation sums for grid quantities (momentum-only)
+template <class ExecutionSpace, class Mesh>
+void gridConservation( MPI_Comm comm, ExecutionSpace exec_space,
+                       const std::shared_ptr<Mesh>& mesh,
+                       Picasso::FieldManager<Mesh>& fm, double& global_mass,
+                       double& global_ke, double& global_pe, const int pe_dir,
+                       const double zero_ground,
+                       const Kokkos::Array<double, 3> gravity )
+{
+    auto own_index_space = mesh->localGrid()->indexSpace(
+        Cabana::Grid::Own(), Cabana::Grid::Node(), Cabana::Grid::Local() );
+    auto m_i =
+        fm.view( Picasso::FieldLocation::Node(), Picasso::Field::Mass() );
+    auto v_i = fm.view( Picasso::FieldLocation::Node(),
+                        Picasso::Field::Velocity<3>() );
+    auto x_i = fm.view( Picasso::FieldLocation::Node(),
+                        Picasso::Field::PhysicalPosition<3>() );
+    double mass_integral = 0;
+    double ke_integral = 0;
+    double pe_integral = 0;
+
+    // FIXME: Collapse the parallel reduce calls once that capability is
+    // available
+    Cabana::Grid::grid_parallel_reduce(
+        "Picasso:MassIntegral", exec_space, own_index_space,
+        KOKKOS_LAMBDA( const int i, const int j, const int k,
+                       double& mass_cont ) { mass_cont += m_i( i, j, k, 0 ); },
+        mass_integral );
+
+    Cabana::Grid::grid_parallel_reduce(
+        "Picasso:KineticEnergyIntegral", exec_space, own_index_space,
+        KOKKOS_LAMBDA( const int i, const int j, const int k,
+                       double& ke_cont ) {
+            double mag_v = sqrt( v_i( i, j, k, 0 ) * v_i( i, j, k, 0 ) +
+                                 v_i( i, j, k, 1 ) * v_i( i, j, k, 1 ) +
+                                 v_i( i, j, k, 2 ) * v_i( i, j, k, 2 ) );
+            ke_cont += m_i( i, j, k, 0 ) * mag_v * mag_v / 2;
+        },
+        ke_integral );
+
+    Cabana::Grid::grid_parallel_reduce(
+        "Picasso:PotentialEnergyIntegral", exec_space, own_index_space,
+        KOKKOS_LAMBDA( const int i, const int j, const int k,
+                       double& pe_cont ) {
+            pe_cont += m_i( i, j, k, 0 ) * Kokkos::abs( gravity[pe_dir] ) *
+                       ( x_i( i, j, k, pe_dir ) - zero_ground );
+        },
+        pe_integral );
+
+    MPI_Allreduce( &mass_integral, &global_mass, 1, MPI_DOUBLE, MPI_SUM, comm );
+    MPI_Allreduce( &ke_integral, &global_ke, 1, MPI_DOUBLE, MPI_SUM, comm );
+    MPI_Allreduce( &pe_integral, &global_pe, 1, MPI_DOUBLE, MPI_SUM, comm );
+}
+
+} // end namespace Picasso
+
+#endif // end PICASSO_CONSERVATION_HPP
diff --git a/src/Picasso_FieldTypes.hpp b/src/Picasso_FieldTypes.hpp
index 11bc6c01..333afad2 100644
--- a/src/Picasso_FieldTypes.hpp
+++ b/src/Picasso_FieldTypes.hpp
@@ -966,9 +966,119 @@ struct CommRank : Scalar<int>
     static std::string label() { return "comm_rank"; }
 };
 
+struct Mass : Field::Scalar<double>
+{
+    static std::string label() { return "mass"; }
+};
+
+struct Pressure : Field::Scalar<double>
+{
+    static std::string label() { return "pressure"; }
+};
+
+struct Volume : Field::Scalar<double>
+{
+    static std::string label() { return "volume"; }
+};
+
+template <std::size_t NumSpaceDim>
+struct Velocity : Field::Vector<double, NumSpaceDim>
+{
+    static std::string label() { return "velocity"; }
+};
+
+template <std::size_t NumSpaceDim>
+struct OldVelocity : Field::Vector<double, NumSpaceDim>
+{
+    static std::string label() { return "old_velocity"; }
+};
+
+template <std::size_t NumSpaceDim>
+struct Stress : Field::Matrix<double, NumSpaceDim, NumSpaceDim>
+{
+    static std::string label() { return "stress"; }
+};
+
+struct DetDefGrad : Field::Scalar<double>
+{
+    static std::string label() { return "determinant_deformation_gradient"; }
+};
+
 //---------------------------------------------------------------------------//
 
 } // end namespace Field
+
+namespace PolyPIC
+{
+namespace Field
+{
+
+template <std::size_t NumSpaceDim>
+struct Velocity : Picasso::Field::Matrix<double, 8, NumSpaceDim>
+{
+    static std::string label() { return "velocity"; }
+};
+} // namespace Field
+} // namespace PolyPIC
+
+namespace APIC
+{
+
+namespace Field
+{
+
+template <std::size_t NumSpaceDim>
+struct Velocity : Picasso::Field::Matrix<double, 4, NumSpaceDim>
+{
+    static std::string label() { return "velocity"; }
+};
+
+} // namespace Field
+} // namespace APIC
+
+struct PolyPicTag
+{
+};
+
+struct APicTag
+{
+};
+
+struct FlipTag
+{
+};
+
+// Particle indexing for PIC/FLIP vs PolyPIC scalar field
+KOKKOS_INLINE_FUNCTION auto getField( const double field ) { return field; }
+
+template <class FieldType>
+KOKKOS_INLINE_FUNCTION
+    std::enable_if_t<Picasso::LinearAlgebra::is_matrix<FieldType>::value,
+                     typename FieldType::value_type>
+    getField( const FieldType& field )
+{
+    return field( 0, 0 );
+}
+
+// Particle indexing for PIC/FLIP vs PolyPIC vector field
+template <class FieldType>
+KOKKOS_INLINE_FUNCTION
+    std::enable_if_t<Picasso::LinearAlgebra::is_matrix<FieldType>::value,
+                     typename FieldType::value_type>
+    getField( const FieldType& field, const int dir )
+{
+    return field( 0, dir );
+}
+
+template <class FieldType>
+KOKKOS_INLINE_FUNCTION
+    std::enable_if_t<Picasso::LinearAlgebra::is_vector<FieldType>::value,
+                     typename FieldType::value_type>
+    getField( const FieldType& field, const int dir )
+{
+    return field( dir );
+}
+
 } // end namespace Picasso
 
 #endif // PICASSO_FIELDTYPES_HPP
diff --git a/src/Picasso_GridKernels.hpp b/src/Picasso_GridKernels.hpp
new file mode 100644
index 00000000..7f5e8a46
--- /dev/null
+++ b/src/Picasso_GridKernels.hpp
@@ -0,0 +1,57 @@
+/****************************************************************************
+ * Copyright (c) 2021 by the Picasso authors                                *
+ * All rights reserved.                                                     *
+ *                                                                          *
+ * This file is part of the Picasso library. Picasso 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                                    *
+ ****************************************************************************/
+
+#ifndef PICASSO_GRIDKERNELS_HPP
+#define PICASSO_GRIDKERNELS_HPP
+
+#include <Picasso_BatchedLinearAlgebra.hpp>
+#include <Picasso_FieldTypes.hpp>
+
+#include <Cabana_Grid.hpp>
+
+#include <type_traits>
+
+namespace Picasso
+{
+//---------------------------------------------------------------------------//
+// Compute nodal velocity from mass-weighted momentum.
+//---------------------------------------------------------------------------//
+template <class MassType = Field::Mass, class VelocityType = Field::Velocity<3>,
+          class OldVelType = Field::OldVelocity<3>>
+struct ComputeGridVelocity
+{
+    template <class LocalMeshType, class GatherDependencies,
+              class ScatterDependencies, class LocalDependencies>
+    KOKKOS_INLINE_FUNCTION void
+    operator()( const LocalMeshType&, const GatherDependencies&,
+                const ScatterDependencies&, const LocalDependencies& local_deps,
+                const int i, const int j, const int k ) const
+    {
+        // Get the local dependencies.
+        auto m_i = local_deps.get( Picasso::FieldLocation::Node(), MassType{} );
+        auto u_i =
+            local_deps.get( Picasso::FieldLocation::Node(), VelocityType{} );
+        auto old_u_i =
+            local_deps.get( Picasso::FieldLocation::Node(), OldVelType{} );
+
+        // Compute velocity.
+        for ( int d = 0; d < 3; ++d )
+        {
+            old_u_i( i, j, k, d ) = ( m_i( i, j, k ) > 0.0 )
+                                        ? old_u_i( i, j, k, d ) / m_i( i, j, k )
+                                        : 0.0;
+            u_i( i, j, k, d ) = old_u_i( i, j, k, d );
+        }
+    }
+};
+} // namespace Picasso
+
+#endif // end PICASSO_GRIDKERNELS_HPP
diff --git a/src/Picasso_InputParser.hpp b/src/Picasso_InputParser.hpp
index 1888213e..c01214d4 100644
--- a/src/Picasso_InputParser.hpp
+++ b/src/Picasso_InputParser.hpp
@@ -17,9 +17,20 @@
 #include <fstream>
 #include <string>
 
+#include <Kokkos_Core.hpp>
+
 namespace Picasso
 {
 
+template <typename T, std::size_t N>
+Kokkos::Array<T, N> copy( const std::array<T, N> in )
+{
+    Kokkos::Array<T, N> out;
+    for ( std::size_t d = 0; d < in.size(); ++d )
+        out[d] = in[d];
+    return out;
+}
+
 inline nlohmann::json parse( const std::string& filename )
 {
     std::ifstream stream( filename );
diff --git a/src/Picasso_InterpolationKernels.hpp b/src/Picasso_InterpolationKernels.hpp
new file mode 100644
index 00000000..14500ca0
--- /dev/null
+++ b/src/Picasso_InterpolationKernels.hpp
@@ -0,0 +1,435 @@
+/****************************************************************************
+ * Copyright (c) 2021 by the Picasso authors                                *
+ * All rights reserved.                                                     *
+ *                                                                          *
+ * This file is part of the Picasso library. Picasso 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                                    *
+ ****************************************************************************/
+
+#ifndef PICASSO_INTERPOLATIONKERNELS_HPP
+#define PICASSO_INTERPOLATIONKERNELS_HPP
+
+#include <Picasso_APIC.hpp>
+#include <Picasso_BatchedLinearAlgebra.hpp>
+#include <Picasso_FieldTypes.hpp>
+#include <Picasso_ParticleInterpolation.hpp>
+#include <Picasso_PolyPIC.hpp>
+
+#include <Cabana_Grid.hpp>
+
+#include <type_traits>
+
+namespace Picasso
+{
+//---------------------------------------------------------------------------//
+// P2G
+//---------------------------------------------------------------------------//
+
+//---------------------------------------------------------------------------//
+template <int InterpolationOrder, class ParticleFieldType, class OldFieldType,
+          class MassType, class PositionType, class InterpolationType>
+struct Particle2Grid;
+
+//---------------------------------------------------------------------------//
+// Project particle enthalpy/momentum to grid. PolyPIC variant
+//---------------------------------------------------------------------------//
+template <int InterpolationOrder, class ParticleFieldType, class OldFieldType,
+          class MassType, class PositionType>
+struct Particle2Grid<InterpolationOrder, ParticleFieldType, OldFieldType,
+                     MassType, PositionType, PolyPicTag>
+{
+    // Explicit time step.
+    double _dt;
+
+    // Constructor
+    Particle2Grid( const double dt )
+        : _dt( dt )
+    {
+    }
+
+    template <class LocalMeshType, class GatherDependencies,
+              class ScatterDependencies, class LocalDependencies,
+              class ParticleViewType>
+    KOKKOS_INLINE_FUNCTION void
+    operator()( const LocalMeshType& local_mesh, const GatherDependencies&,
+                const ScatterDependencies& scatter_deps,
+                const LocalDependencies&, ParticleViewType& particle ) const
+    {
+        // Get particle data.
+        auto f_p = Picasso::get( particle, ParticleFieldType() );
+        // Require the PolyPIC specific field-type.
+        auto v_p = Picasso::get( particle, PolyPIC::Field::Velocity<3>() );
+        auto m_p = Picasso::get( particle, MassType() );
+        auto x_p = Picasso::get( particle, PositionType() );
+
+        // Get the scatter dependencies.
+        auto m_i =
+            scatter_deps.get( Picasso::FieldLocation::Node(), MassType() );
+        auto f_i =
+            scatter_deps.get( Picasso::FieldLocation::Node(), OldFieldType() );
+
+        // Node interpolant.
+        auto spline = Picasso::createSpline(
+            Picasso::FieldLocation::Node(),
+            Picasso::InterpolationOrder<InterpolationOrder>(), local_mesh, x_p,
+            Picasso::SplineValue(), Picasso::SplineDistance() );
+
+        // Interpolate mass and mass-weighted enthalpy/momentum to grid with
+        // PolyPIC.
+        Picasso::PolyPIC::p2g( m_p, v_p, f_p, f_i, m_i, _dt, spline );
+    }
+};
+
+//---------------------------------------------------------------------------//
+// Project particle enthalpy/momentum to grid. APIC variant
+//---------------------------------------------------------------------------//
+template <int InterpolationOrder, class ParticleFieldType, class OldFieldType,
+          class MassType, class PositionType>
+struct Particle2Grid<InterpolationOrder, ParticleFieldType, OldFieldType,
+                     MassType, PositionType, APicTag>
+{
+    // Explicit time step.
+    double _dt;
+
+    // Constructor
+    Particle2Grid( const double dt )
+        : _dt( dt )
+    {
+    }
+
+    template <class LocalMeshType, class GatherDependencies,
+              class ScatterDependencies, class LocalDependencies,
+              class ParticleViewType>
+    KOKKOS_INLINE_FUNCTION void
+    operator()( const LocalMeshType& local_mesh, const GatherDependencies&,
+                const ScatterDependencies& scatter_deps,
+                const LocalDependencies&, ParticleViewType& particle ) const
+    {
+        // Get particle data.
+        auto f_p = Picasso::get( particle, ParticleFieldType() );
+        auto m_p = Picasso::get( particle, MassType() );
+        auto x_p = Picasso::get( particle, PositionType() );
+
+        // Get the scatter dependencies.
+        auto m_i =
+            scatter_deps.get( Picasso::FieldLocation::Node(), MassType() );
+        auto f_i =
+            scatter_deps.get( Picasso::FieldLocation::Node(), OldFieldType() );
+
+        // Node interpolant.
+        auto spline = Picasso::createSpline(
+            Picasso::FieldLocation::Node(),
+            Picasso::InterpolationOrder<InterpolationOrder>(), local_mesh, x_p,
+            Picasso::SplineValue(), Picasso::SplineDistance(),
+            Picasso::SplineGradient() );
+
+        // Interpolate mass and mass-weighted enthalpy/momentum to grid with
+        // APIC.
+        Picasso::APIC::p2g( m_p, f_p, m_i, f_i, spline );
+    }
+};
+
+//---------------------------------------------------------------------------//
+// Project particle enthalpy/momentum to grid. FLIP/PIC variant
+//---------------------------------------------------------------------------//
+template <int InterpolationOrder, class ParticleFieldType, class OldFieldType,
+          class MassType, class PositionType>
+struct Particle2Grid<InterpolationOrder, ParticleFieldType, OldFieldType,
+                     MassType, PositionType, FlipTag>
+{
+    // Explicit time step.
+    double _dt;
+
+    // Constructor
+    Particle2Grid( const double dt )
+        : _dt( dt )
+    {
+    }
+
+    template <class LocalMeshType, class GatherDependencies,
+              class ScatterDependencies, class LocalDependencies,
+              class ParticleViewType>
+    KOKKOS_INLINE_FUNCTION void
+    operator()( const LocalMeshType& local_mesh, const GatherDependencies&,
+                const ScatterDependencies& scatter_deps,
+                const LocalDependencies&, ParticleViewType& particle ) const
+    {
+        // Get particle data.
+        auto f_p = Picasso::get( particle, ParticleFieldType() );
+        auto m_p = Picasso::get( particle, MassType() );
+        auto x_p = Picasso::get( particle, PositionType() );
+
+        // Get the scatter dependencies.
+        auto m_i =
+            scatter_deps.get( Picasso::FieldLocation::Node(), MassType() );
+        auto f_i =
+            scatter_deps.get( Picasso::FieldLocation::Node(), OldFieldType() );
+
+        // Node interpolant.
+        auto spline = Picasso::createSpline(
+            Picasso::FieldLocation::Node(),
+            Picasso::InterpolationOrder<InterpolationOrder>(), local_mesh, x_p,
+            Picasso::SplineValue(), Picasso::SplineDistance() );
+
+        // Interpolate mass and mass-weighted enthalpy/momentum to grid.
+        Picasso::P2G::value( spline, m_p, m_i );
+        Picasso::P2G::value( spline, m_p * f_p, f_i );
+    }
+};
+
+// P2G creation function (used instead of direct class defaults because of
+// partial specializations).
+template <int InterpolationOrder, class InterpolationType,
+          class ParticleFieldType,
+          class OldFieldType = Picasso::Field::OldVelocity<3>,
+          class MassType = Picasso::Field::Mass,
+          class PositionType = Picasso::Field::LogicalPosition<3>>
+auto createParticle2Grid( const double dt )
+{
+    return Particle2Grid<InterpolationOrder, ParticleFieldType, OldFieldType,
+                         MassType, PositionType, InterpolationType>( dt );
+}
+
+//---------------------------------------------------------------------------//
+// G2P
+//---------------------------------------------------------------------------//
+
+template <int InterpolationOrder, class MassType, class PositionType,
+          class GridVelocity, class GridPosition, class InterpolationType>
+struct Grid2ParticleVelocity;
+
+//---------------------------------------------------------------------------//
+// Update particle state. PolyPIC variant.
+//---------------------------------------------------------------------------//
+template <int InterpolationOrder, class MassType, class PositionType,
+          class GridVelocity, class GridPosition>
+struct Grid2ParticleVelocity<InterpolationOrder, MassType, PositionType,
+                             GridVelocity, GridPosition, PolyPicTag>
+{
+    // Explicit time step.
+    double _dt;
+
+    // Primary constructor
+    Grid2ParticleVelocity( const double dt )
+        : _dt( dt )
+    {
+    }
+
+    // Constructor for easier compatibility with FLIP
+    Grid2ParticleVelocity( const double dt, const double )
+        : _dt( dt )
+    {
+    }
+
+    template <class LocalMeshType, class GatherDependencies,
+              class ScatterDependencies, class LocalDependencies,
+              class ParticleViewType>
+    KOKKOS_INLINE_FUNCTION void
+    operator()( const LocalMeshType& local_mesh,
+                const GatherDependencies& gather_deps,
+                const ScatterDependencies&, const LocalDependencies& local_deps,
+                ParticleViewType& particle ) const
+    {
+        // Get particle data.
+        auto u_p = Picasso::get( particle, PolyPIC::Field::Velocity<3>() );
+        auto x_p = Picasso::get( particle, PositionType() );
+
+        // Get the gather dependencies.
+        auto m_i =
+            gather_deps.get( Picasso::FieldLocation::Node(), MassType() );
+        auto u_i =
+            gather_deps.get( Picasso::FieldLocation::Node(), GridVelocity() );
+        // Get the local dependencies for getting physcial location of node
+        auto x_i =
+            local_deps.get( Picasso::FieldLocation::Node(), GridPosition() );
+
+        // Node interpolant.
+        auto spline = Picasso::createSpline(
+            Picasso::FieldLocation::Node(),
+            Picasso::InterpolationOrder<InterpolationOrder>(), local_mesh, x_p,
+            Picasso::SplineValue(), Picasso::SplineGradient() );
+
+        // Update particle velocity using a PolyPIC update.
+        Picasso::PolyPIC::g2p( u_i, u_p, spline );
+
+        // Update particle position.
+        auto x_i_updated =
+            [=]( const int i, const int j, const int k, const int d )
+        { return x_i( i, j, k, d ) + _dt * u_i( i, j, k, d ); };
+        Picasso::G2P::value( spline, x_i_updated, x_p );
+    }
+};
+
+//---------------------------------------------------------------------------//
+// Update particle state. APIC variant.
+//---------------------------------------------------------------------------//
+template <int InterpolationOrder, class MassType, class PositionType,
+          class GridVelocity, class GridPosition>
+struct Grid2ParticleVelocity<InterpolationOrder, MassType, PositionType,
+                             GridVelocity, GridPosition, APicTag>
+{
+    // Explicit time step.
+    double _dt;
+
+    // Primary constructor
+    Grid2ParticleVelocity( const double dt )
+        : _dt( dt )
+    {
+    }
+
+    // Constructor for easier compatibility with FLIP
+    Grid2ParticleVelocity( const double dt, const double )
+        : _dt( dt )
+    {
+    }
+
+    template <class LocalMeshType, class GatherDependencies,
+              class ScatterDependencies, class LocalDependencies,
+              class ParticleViewType>
+    KOKKOS_INLINE_FUNCTION void
+    operator()( const LocalMeshType& local_mesh,
+                const GatherDependencies& gather_deps,
+                const ScatterDependencies&, const LocalDependencies& local_deps,
+                ParticleViewType& particle ) const
+    {
+        // Get particle data.
+        auto f_p = Picasso::get( particle, APIC::Field::Velocity<3>() );
+        auto x_p = Picasso::get( particle, PositionType() );
+
+        // Get the gather dependencies.
+        auto m_i =
+            gather_deps.get( Picasso::FieldLocation::Node(), MassType() );
+        auto u_i =
+            gather_deps.get( Picasso::FieldLocation::Node(), GridVelocity() );
+        // Get the local dependencies for getting physcial location of node
+        auto x_i =
+            local_deps.get( Picasso::FieldLocation::Node(), GridPosition() );
+
+        // Node interpolant.
+        auto spline = Picasso::createSpline(
+            Picasso::FieldLocation::Node(),
+            Picasso::InterpolationOrder<InterpolationOrder>(), local_mesh, x_p,
+            Picasso::SplineValue(), Picasso::SplineDistance(),
+            Picasso::SplineGradient() );
+
+        // Update particle velocity using a APIC update.
+        Picasso::APIC::g2p( u_i, f_p, spline );
+
+        // Update particle position.
+        auto x_i_updated =
+            [=]( const int i, const int j, const int k, const int d )
+        { return x_i( i, j, k, d ) + _dt * u_i( i, j, k, d ); };
+        Picasso::G2P::value( spline, x_i_updated, x_p );
+    }
+};
+
+//---------------------------------------------------------------------------//
+// Update particle state. FLIP/PIC variant.
+//---------------------------------------------------------------------------//
+template <int InterpolationOrder, class MassType, class PositionType,
+          class GridVelocity, class GridPosition>
+struct Grid2ParticleVelocity<InterpolationOrder, MassType, PositionType,
+                             GridVelocity, GridPosition, FlipTag>
+{
+    // Explicit time step.
+    double _dt;
+    // FLIP/PIC ratio
+    double _beta = 0.99;
+
+    // Primary constructor
+    Grid2ParticleVelocity( const double dt, const double beta )
+        : _dt( dt )
+        , _beta( beta )
+    {
+    }
+
+    // Default beta constructor
+    Grid2ParticleVelocity( const double dt )
+        : _dt( dt )
+    {
+    }
+
+    template <class LocalMeshType, class GatherDependencies,
+              class ScatterDependencies, class LocalDependencies,
+              class ParticleViewType>
+    KOKKOS_INLINE_FUNCTION void
+    operator()( const LocalMeshType& local_mesh,
+                const GatherDependencies& gather_deps,
+                const ScatterDependencies&, const LocalDependencies& local_deps,
+                ParticleViewType& particle ) const
+    {
+        // Get particle data.
+        auto u_p = Picasso::get( particle, Field::Velocity<3>{} );
+        auto x_p = Picasso::get( particle, PositionType{} );
+
+        // Get the gather dependencies.
+        auto m_i =
+            gather_deps.get( Picasso::FieldLocation::Node(), MassType() );
+        auto u_i =
+            gather_deps.get( Picasso::FieldLocation::Node(), GridVelocity() );
+        auto old_u_i = gather_deps.get( Picasso::FieldLocation::Node(),
+                                        Field::OldVelocity<3>() );
+
+        // Get the local dependencies for getting physcial location of node
+        auto x_i =
+            local_deps.get( Picasso::FieldLocation::Node(), GridPosition() );
+
+        // Node interpolant.
+        auto spline = Picasso::createSpline(
+            Picasso::FieldLocation::Node(),
+            Picasso::InterpolationOrder<InterpolationOrder>(), local_mesh, x_p,
+            Picasso::SplineValue(), Picasso::SplineGradient() );
+
+        // Update particle velocity using a hybrid PIC/FLIP update.
+        // Note the need for 4 indices here because this is passed to the Cabana
+        // grid.
+        auto d_u_i = [=]( const int i, const int j, const int k, const int d )
+        {
+            return ( m_i( i, j, k ) > 0.0 )
+                       ? u_i( i, j, k, d ) - old_u_i( i, j, k, d )
+                       : 0.0;
+        };
+
+        Picasso::Vec3<double> u_p_pic;
+        // auto u_i_pic = [=]( const int i, const int j, const int k,
+        //                    const int d ) { return u_i( i, j, k, d ); };
+        Picasso::G2P::value( spline, u_i, u_p_pic );
+
+        Picasso::Vec3<double> d_u_p;
+        Picasso::G2P::value( spline, d_u_i, d_u_p );
+
+        Picasso::Vec3<double> u_p_flip;
+
+        u_p_flip = u_p + d_u_p;
+
+        // Update particle velocity.
+        u_p = _beta * u_p_flip + ( 1.0 - _beta ) * u_p_pic;
+
+        // Update particle position.
+        auto x_i_updated =
+            [=]( const int i, const int j, const int k, const int d )
+        { return x_i( i, j, k, d ) + _dt * u_i( i, j, k, d ); };
+        Picasso::G2P::value( spline, x_i_updated, x_p );
+    }
+};
+
+// G2P creation function (used instead of direct class defaults because of
+// partial specializations).
+template <int InterpolationOrder, class InterpolationType,
+          class MassType = Picasso::Field::Mass,
+          class PositionType = Picasso::Field::LogicalPosition<3>,
+          class GridVelocity = Picasso::Field::Velocity<3>,
+          class GridPosition = Picasso::Field::PhysicalPosition<3>>
+auto createGrid2ParticleVelocity( const double dt, const double beta )
+{
+    return Grid2ParticleVelocity<InterpolationOrder, MassType, PositionType,
+                                 GridVelocity, GridPosition, InterpolationType>(
+        dt, beta );
+}
+
+} // end namespace Picasso
+
+#endif // end PICASSO_INTERPOLATIONKERNELS_HPP
diff --git a/src/Picasso_ParticleInterpolation.hpp b/src/Picasso_ParticleInterpolation.hpp
index fa41a0c6..b8cf5a8c 100644
--- a/src/Picasso_ParticleInterpolation.hpp
+++ b/src/Picasso_ParticleInterpolation.hpp
@@ -278,8 +278,6 @@ KOKKOS_INLINE_FUNCTION void divergence( const SplineDataType& sd,
 
 } // end namespace P2G
 
-//---------------------------------------------------------------------------//
-
 } // end namespace Picasso
 
 #endif // end PICASSO_PARTICLEINTERPOLATION_HPP