From 43423253521202c70b48902fefbdeac7ef3065f1 Mon Sep 17 00:00:00 2001 From: kwitaechong Date: Sun, 1 Sep 2024 19:12:23 -0400 Subject: [PATCH 01/22] initial setup dam break problem --- CMakeLists.txt | 2 +- examples/CMakeLists.txt | 12 + examples/dam_break.cpp | 122 ++++ examples/inputs/particle_init_test.json | 10 + examples/sources/CMakeLists.txt | 0 examples/sources/Example_Sources.hpp | 4 + examples/sources/Particle_Field.hpp | 87 +++ examples/sources/Particle_Init.hpp | 60 ++ .../sources/Picasso_BoundaryCondition.hpp | 52 ++ .../Picasso_ExplicitMomentumUpdate.hpp | 681 ++++++++++++++++++ examples/sources/Picasso_Output.hpp | 65 ++ 11 files changed, 1094 insertions(+), 1 deletion(-) create mode 100644 examples/dam_break.cpp create mode 100644 examples/inputs/particle_init_test.json create mode 100644 examples/sources/CMakeLists.txt create mode 100644 examples/sources/Example_Sources.hpp create mode 100644 examples/sources/Particle_Field.hpp create mode 100644 examples/sources/Particle_Init.hpp create mode 100644 examples/sources/Picasso_BoundaryCondition.hpp create mode 100644 examples/sources/Picasso_ExplicitMomentumUpdate.hpp create mode 100644 examples/sources/Picasso_Output.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 3f45e6d7..267dfe18 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -110,7 +110,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 example/*.[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..9795e0e6 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -0,0 +1,12 @@ +add_subdirectory(sources) + +configure_file( + ${CMAKE_CURRENT_SOURCE_DIR}/inputs/particle_init_test.json + ${CMAKE_CURRENT_BINARY_DIR}/particle_init_test.json + COPYONLY) + +add_executable( DamBreak dam_break.cpp ) +target_link_libraries( DamBreak PRIVATE 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..dc011b85 --- /dev/null +++ b/examples/dam_break.cpp @@ -0,0 +1,122 @@ +#include + +#include "sources/Particle_Init.hpp" +#include "sources/Picasso_BoundaryCondition.hpp" +#include "sources/Picasso_ExplicitMomentumUpdate.hpp" +#include "sources/Picasso_Output.hpp" + +#include +#include + +#include + +#include + +using namespace Picasso; + +//---------------------------------------------------------------------------// +// DamBreak example +template +void DamBreak() +{ + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = exec_space::memory_space; + + // Global bounding box. + double cell_size = 0.05; + std::array global_num_cell = { 20, 20, 20 }; + std::array global_low_corner = { 0.0, 0.0, 0.0 }; + std::array global_high_corner = { + global_low_corner[0] + cell_size * global_num_cell[0], + global_low_corner[1] + cell_size * global_num_cell[1], + global_low_corner[2] + cell_size * global_num_cell[2] }; + + // Get inputs for mesh. + InputParser parser( "dam_break.json", "json" ); + const auto& pt = parser.propertyTree(); + + Kokkos::Array global_box = { + global_low_corner[0], global_low_corner[1], global_low_corner[2], + global_high_corner[0], global_high_corner[1], global_high_corner[2] }; + int minimum_halo_size = 0; + + // Make mesh. + auto mesh = createUniformMesh( memory_space(), pt, global_box, + minimum_halo_size, MPI_COMM_WORLD ); + + // Make a particle list. + Cabana::ParticleTraits + fields; + auto particles = Cabana::Grid::createParticleList( + "test_particles", fields ); + + // Initialize particles + Kokkos::Array block = { 0.0, 0.0, 0.0, 0.4, 0.4, 0.4 }; + double density = 1e3; + + auto momentum_init_functor = + createParticleInitFunc( particles, ParticleVelocity(), block, density ); + + int ppc = 2; + auto local_grid = mesh->localGrid(); + Cabana::Grid::createParticles( Cabana::InitUniform(), exec_space(), + momentum_init_functor, particles, ppc, + *local_grid ); + + // Boundary index space + auto bc_index = local_grid->boundaryIndexSpace( + Cabana::Grid::Own(), Cabana::Grid::Node(), -1, 1, 0 ); + using bc_index_type = decltype( bc_index ); + + BoundaryCondition bc{ bc_index }; + + // Properties + double gamma = 7.0; + double bulk_modulus = 1.0e+5; + Kokkos::Array gravity = { 0.0, 0.0, -9.8 }; + Picasso::Properties props( gamma, bulk_modulus, gravity ); + + // Time integragor + auto time_integrator = Picasso::createExplicitMomentumIntegrator( + mesh, InterpolationType(), ParticleVelocity(), props, 1.0e-4 ); + auto fm = Picasso::createFieldManager( mesh ); + time_integrator.setup( *fm ); + + // steps + while ( time_integrator.totalTime() < 1.0 ) + { + // Write particle fields. + Picasso::Output::outputParticles( + MPI_COMM_WORLD, exec_space(), ParticleVelocity(), + time_integrator.totalSteps(), 10, time_integrator.totalTime(), + particles ); + + printf( "aaa\n" ); + + // Step. + time_integrator.step( exec_space(), *fm, particles, *local_grid, bc ); + } +} + +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]; + + DamBreak(); + + Kokkos::finalize(); + MPI_Finalize(); + + return 0; +} diff --git a/examples/inputs/particle_init_test.json b/examples/inputs/particle_init_test.json new file mode 100644 index 00000000..a21c2c34 --- /dev/null +++ b/examples/inputs/particle_init_test.json @@ -0,0 +1,10 @@ +{ + "mesh": { + "cell_size": 0.2, + "periodic": [false, false, false], + "partitioner": { + "type": "uniform_dim" + }, + "halo_cell_width": 0 + } +} diff --git a/examples/sources/CMakeLists.txt b/examples/sources/CMakeLists.txt new file mode 100644 index 00000000..e69de29b diff --git a/examples/sources/Example_Sources.hpp b/examples/sources/Example_Sources.hpp new file mode 100644 index 00000000..bdcadbec --- /dev/null +++ b/examples/sources/Example_Sources.hpp @@ -0,0 +1,4 @@ +#include "Particle_Field.hpp" +#include "Particle_Init.hpp" +#include "Picasso_BoundaryCondition.hpp" +#include "Picasso_ExplicitMomentumUpdate.hpp" diff --git a/examples/sources/Particle_Field.hpp b/examples/sources/Particle_Field.hpp new file mode 100644 index 00000000..0e4c93b3 --- /dev/null +++ b/examples/sources/Particle_Field.hpp @@ -0,0 +1,87 @@ +#include + +#include +#include + +#include + +#include + +namespace Picasso +{ + +namespace Field +{ + +struct Mass : Field::Scalar +{ + static std::string label() { return "Mass"; } +}; + +struct Pressure : Field::Scalar +{ + static std::string label() { return "Pressure"; } +}; + +struct Volume : Field::Scalar +{ + static std::string label() { return "Volume"; } +}; + +struct Velocity : Field::Vector +{ + static std::string label() { return "Velocity"; } +}; + +struct OldU : Field::Vector +{ + static std::string label() { return "Old_Velocity"; } +}; + +struct DeltaUGravity : Field::Vector +{ + static std::string label() { return "velocity_change_from_gravity"; } +}; + +struct Stress : Field::Matrix +{ + static std::string label() { return "stress"; } +}; + +struct DeltaUStress : Field::Vector +{ + static std::string label() { return "velocity_change_from_stress"; } +}; + +struct DetDefGrad : Field::Scalar +{ + static std::string label() { return "Det_deformation_gradient"; } +}; + +using Position = Picasso::Field::LogicalPosition<3>; + +} // namespace Field + +namespace APIC +{ + +struct APicTag +{ +}; + +namespace Field +{ + +struct Velocity : Picasso::Field::Matrix +{ + static std::string label() { return "velocity"; } +}; + +} // namespace Field +} // namespace APIC + +struct FlipTag +{ +}; + +} // namespace Picasso diff --git a/examples/sources/Particle_Init.hpp b/examples/sources/Particle_Init.hpp new file mode 100644 index 00000000..77bd6ac7 --- /dev/null +++ b/examples/sources/Particle_Init.hpp @@ -0,0 +1,60 @@ +#include + +#include "Particle_Field.hpp" + +#include +#include + +#include + +#include + +namespace Picasso +{ + +template +struct ParticleInitFunc +{ + Kokkos::Array block; + double density; + + ParticleInitFunc( const Kokkos::Array& _block, + const double _density ) + : block( _block ) + , density( _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() ) = 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::Position(), d ) = x[d]; + return true; + } + + return false; + } +}; + +template +auto createParticleInitFunc( const ParticleType&, const VelocityType&, + const Kokkos::Array& block, + const double density ) +{ + return ParticleInitFunc{ block, density }; +} + +} // namespace Picasso diff --git a/examples/sources/Picasso_BoundaryCondition.hpp b/examples/sources/Picasso_BoundaryCondition.hpp new file mode 100644 index 00000000..a745b240 --- /dev/null +++ b/examples/sources/Picasso_BoundaryCondition.hpp @@ -0,0 +1,52 @@ +#ifndef PICASSO_BOUNDARYCONDITION_HPP +#define PICASSO_BOUNDARYCONDITION_HPP + +#include +#include + +namespace Picasso +{ + +template +struct BoundaryCondition +{ + BoundaryIndexType bc_index; + + template + KOKKOS_INLINE_FUNCTION void apply( ViewType view, const int i, const int j, + const int k ) const + { + // x faces + if ( bc_index.min( Cabana::Grid::Dim::I ) <= i && + i < bc_index.max( Cabana::Grid::Dim::I ) ) + view( i, j, k, 0 ) = 0.0; + // y faces + if ( bc_index.min( Cabana::Grid::Dim::J ) <= j && + j < bc_index.max( Cabana::Grid::Dim::J ) ) + view( i, j, k, 1 ) = 0.0; + // z faces + if ( bc_index.min( Cabana::Grid::Dim::K ) <= k && + k < bc_index.max( Cabana::Grid::Dim::K ) ) + view( i, j, k, 2 ) = 0.0; + } +}; + +struct Properties +{ + double gamma; + double bulk_modulus; + Kokkos::Array gravity; + + KOKKOS_INLINE_FUNCTION + Properties( const double _gamma, const double _bulk_modulus, + const Kokkos::Array _gravity ) + : gamma( _gamma ) + , bulk_modulus( _bulk_modulus ) + , gravity( _gravity ) + { + } +}; + +} // end namespace Picasso + +#endif // end PICASSO_BOUNDARYCONDITION_HPP diff --git a/examples/sources/Picasso_ExplicitMomentumUpdate.hpp b/examples/sources/Picasso_ExplicitMomentumUpdate.hpp new file mode 100644 index 00000000..911c2a6b --- /dev/null +++ b/examples/sources/Picasso_ExplicitMomentumUpdate.hpp @@ -0,0 +1,681 @@ +#include + +#include "Picasso_BoundaryCondition.hpp" + +#include +#include + +#include + +#include + +namespace Picasso +{ + +//---------------------------------------------------------------------------// +// Compute nodal velocity from mass-weighted momentum. +//---------------------------------------------------------------------------// +struct ComputeGridVelocity +{ + template + 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(), Field::Mass() ); + auto u_i = + local_deps.get( Picasso::FieldLocation::Node(), Field::Velocity() ); + auto old_u_i = + local_deps.get( Picasso::FieldLocation::Node(), Field::OldU() ); + + // 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 ); + } + } +}; + +//---------------------------------------------------------------------------// +// Update particle stress. +//---------------------------------------------------------------------------// +template +struct ComputeParticlePressure +{ + double dt; + + // Stress property functions + StressProperties properties; + + template + 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, Field::Position() ); + auto& J_p = Picasso::get( particle, Field::DetDefGrad() ); + auto& p_p = Picasso::get( particle, Field::Pressure() ); + auto s_p = Picasso::get( particle, Field::Stress() ); + auto v_p = Picasso::get( particle, Field::Volume() ); + auto m_p = Picasso::get( particle, Field::Mass() ); + + // Get the gather dependencies. + auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), + Field::Velocity() ); + + // update strain rate + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::InterpolationOrder(), local_mesh, x_p, + Picasso::SplineValue(), Picasso::SplineGradient() ); + Picasso::Mat3 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 = properties.bulk_modulus * + ( Kokkos::pow( J_p, -properties.gamma ) - 1.0 ); + + Picasso::Mat3 I; + Picasso::LinearAlgebra::identity( I ); + + s_p = -p_p * I; + + printf( "particle %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e\n", J_p, + p_p, v_p, m_p, x_p( 0 ), x_p( 1 ), x_p( 2 ) ); + } +}; + +//---------------------------------------------------------------------------// +// Grid momentum change due to stress +//---------------------------------------------------------------------------// +template +struct ComputeGridVelocityChangeStress +{ + template + KOKKOS_INLINE_FUNCTION void + operator()( const LocalMeshType& local_mesh, const GatherDependencies&, + const ScatterDependencies& scatter_deps, + const LocalDependencies&, ParticleViewType& particle ) const + { + // Get particle data. + auto v_p = Picasso::get( particle, Field::Volume() ); + auto x_p = Picasso::get( particle, Field::Position() ); + auto s_p = Picasso::get( particle, Field::Stress() ); + + // Get the scatter dependencies. + auto delta_u_s_i = scatter_deps.get( Picasso::FieldLocation::Node(), + Field::DeltaUStress() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::InterpolationOrder(), local_mesh, x_p, + Picasso::SplineValue(), Picasso::SplineGradient() ); + + // Compute velocity update. + Picasso::P2G::divergence( spline, -v_p * s_p, delta_u_s_i ); + } +}; + +//---------------------------------------------------------------------------// +// Grid momentum change due to gravity +//---------------------------------------------------------------------------// +template +struct ComputeGridVelocityChangeGravity +{ + // Stress property functions + StressProperties properties; + + template + 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, Field::Mass() ); + auto x_p = Picasso::get( particle, Field::Position() ); + + // Get the scatter dependencies. + auto delta_u_g_i = scatter_deps.get( Picasso::FieldLocation::Node(), + Field::DeltaUGravity() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::InterpolationOrder(), local_mesh, x_p, + Picasso::SplineValue(), Picasso::SplineGradient() ); + + // Compute velocity update. + Picasso::Vec3 g_p = { properties.gravity[0], + properties.gravity[1], + properties.gravity[2] }; + + Picasso::P2G::value( spline, m_p * g_p, delta_u_g_i ); + } +}; + +//---------------------------------------------------------------------------// +// Update nodal momentum n+1 with stress, gravity +//---------------------------------------------------------------------------// +struct UpdateGridVelocity +{ + // Explicit time step size. + double dt; + + template + 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 local dependencies. + auto m_i = + gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), + Field::Velocity() ); + auto delta_u_s_i = gather_deps.get( Picasso::FieldLocation::Node(), + Field::DeltaUStress() ); + + auto delta_u_g_i = gather_deps.get( Picasso::FieldLocation::Node(), + Field::DeltaUGravity() ); + + // Compute velocity. + Picasso::Vec3 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; + + printf( "%d %d %d %8.5e %8.5e %8.5e %8.5e\n", i, j, k, m_i( i, j, k ), + delta_u_s_i( i, j, k, 2 ), delta_u_g_i( i, j, k, 2 ), + u_i( i, j, k, 2 ) ); + } +}; + +template +struct Grid2ParticleVelocity; + +//---------------------------------------------------------------------------// +// Update particle state. APIC variant. +//---------------------------------------------------------------------------// +template +struct Grid2ParticleVelocity +{ + // FIXME: only needed for FLIP/PIC + double beta; + + // Explicit time step. + double dt; + + template + 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() ); + auto x_p = Picasso::get( particle, Field::Position() ); + + // Get the gather dependencies. + auto m_i = + gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), + Field::Velocity() ); + // Get the local dependencies for getting physcial location of node + auto x_i = local_deps.get( Picasso::FieldLocation::Node(), + Picasso::Field::PhysicalPosition<3>() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::InterpolationOrder(), local_mesh, x_p, + Picasso::SplineValue(), Picasso::SplineDistance(), + Picasso::SplineGradient() ); + + // Update particle velocity using a PolyPIC 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 ) + { + printf( "g2p %d %d %d %8.5e %8.5e\n", i, j, k, x_i( i, j, k, 2 ), + u_i( i, j, k, 2 ) ); + 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 +struct Grid2ParticleVelocity +{ + double beta = 0.01; + + // Explicit time step. + double dt; + + template + 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() ); + auto x_p = Picasso::get( particle, Field::Position() ); + + // Get the gather dependencies. + auto m_i = + gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), + Field::Velocity() ); + auto old_u_i = + gather_deps.get( Picasso::FieldLocation::Node(), Field::OldU() ); + + // Get the local dependencies for getting physcial location of node + auto x_i = local_deps.get( Picasso::FieldLocation::Node(), + Picasso::Field::PhysicalPosition<3>() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::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 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 d_u_p; + Picasso::G2P::value( spline, d_u_i, d_u_p ); + + Picasso::Vec3 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 ); + } +}; + +//---------------------------------------------------------------------------// +// Compute boundary condition on grid. +//---------------------------------------------------------------------------// +template +struct ApplyBoundaryCondition +{ + BCType bc; + + template + 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 f_i = + local_deps.get( Picasso::FieldLocation::Node(), FieldType() ); + + bc.apply( f_i, i, j, k ); + } +}; + +template +struct Particle2Grid; + +//---------------------------------------------------------------------------// +// Project particle enthalpy/momentum to grid. APIC variant +//---------------------------------------------------------------------------// +template +struct Particle2Grid +{ + // Explicit time step. + double dt; + + template + 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, Field::Mass() ); + auto x_p = Picasso::get( particle, Field::Position() ); + + // Get the scatter dependencies. + auto m_i = + scatter_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + auto f_i = + scatter_deps.get( Picasso::FieldLocation::Node(), OldFieldType() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::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 +struct Particle2Grid +{ + // Explicit time step. + double dt; + + template + 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, Field::Mass() ); + auto x_p = Picasso::get( particle, Field::Position() ); + + // Get the scatter dependencies. + auto m_i = + scatter_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + auto f_i = + scatter_deps.get( Picasso::FieldLocation::Node(), OldFieldType() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::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 ); + } +}; + +//---------------------------------------------------------------------------// +// Explicit momentum time stepper +//---------------------------------------------------------------------------// + +template +class ExplicitMomentumIntegrator +{ + public: + using mass_type = + Picasso::FieldLayout; + using velocity_type = + Picasso::FieldLayout; + using old_u_type = + Picasso::FieldLayout; + using delta_u_s_type = + Picasso::FieldLayout; + using delta_u_g_type = Picasso::FieldLayout; + using grid_type = Picasso::FieldLayout>; + + // Momentum P2G operator types. + using p2g_scatter_deps = + Picasso::ScatterDependencies; + using p2g_u_op = Picasso::GridOperator; + + // Grid velocity operator types. + using compute_u_local_deps = + Picasso::LocalDependencies; + using compute_u_op = Picasso::GridOperator; + + // Update particle stress operator types + using compute_s_gather_deps = Picasso::GatherDependencies; + using compute_s_op = Picasso::GridOperator; + + // Grid boundary condition operator types. + using bc_u_local_deps = + Picasso::LocalDependencies; + using bc_u_op = Picasso::GridOperator; + + // Compute velocity change due to stress operator types. + using compute_du_s_scatter_deps = + Picasso::ScatterDependencies; + using compute_du_s_op = + Picasso::GridOperator; + + // Compute velocity change due to gravity operator types. + using compute_du_g_scatter_deps = + Picasso::ScatterDependencies; + using compute_du_g_op = + Picasso::GridOperator; + + // Grid velocity time update operator types. + using update_u_gather_deps = + Picasso::GatherDependencies; + using update_u_op = Picasso::GridOperator; + + // Particle update G2P types. + using g2p_u_gather_deps = + Picasso::GatherDependencies; + using g2p_u_local_deps = Picasso::LocalDependencies; + using g2p_u_op = + Picasso::GridOperator; + + using property_type = StressProperties; + using interpolation_type = InterpolationType; + using interpolation_variable = ParticleVelocity; + + public: + // Constructor. + ExplicitMomentumIntegrator( const std::shared_ptr& mesh, + StressProperties properties, + const double max_dt, const double cfl_number, + const double beta, const bool fixed_dt ) + : _max_dt( max_dt ) + , _cfl_number( cfl_number ) + , _beta( beta ) + , _fixed_dt( fixed_dt ) + , _p2g_momentum( mesh ) + , _compute_velocity( mesh ) + , _compute_stress( mesh ) + , _compute_delta_u_s( mesh ) + , _compute_delta_u_g( mesh ) + , _apply_bc_momentum( mesh ) + , _update_velocity( mesh ) + , _g2p_momentum( mesh ) + , _properties( properties ) + { + _total_steps = 0; + _total_time = 0.0; + _dt = _max_dt; + _cell_size = mesh->cellSize(); + } + + // Populate the field manager. + void setup( Picasso::FieldManager& fm ) + { + _p2g_momentum.setup( fm ); + _compute_velocity.setup( fm ); + _compute_stress.setup( fm ); + _compute_delta_u_s.setup( fm ); + _compute_delta_u_g.setup( fm ); + _apply_bc_momentum.setup( fm ); + _update_velocity.setup( fm ); + _g2p_momentum.setup( fm ); + } + + // Do a time step. + // Note that the boundary condition is passed here rather than saved as a + // member so that it can be initialized with the FieldManager. + template + void step( const ExecutionSpace& exec_space, + const Picasso::FieldManager& fm, ParticleList& particles, + const LocalGridType& local_grid, const BCType bc ) + { + // Spline interpolation order. + const int spline_order = 1; + + // P2G + Particle2Grid + p2g_func{ _dt }; + _p2g_momentum.apply( "Picasso::p2g_U", + Picasso::FieldLocation::Particle(), exec_space, fm, + particles, p2g_func ); + + // Compute grid velocity. + ComputeGridVelocity compute_u_func; + _compute_velocity.apply( "Picasso::grid_U", + Picasso::FieldLocation::Node(), exec_space, fm, + compute_u_func ); + + // Update particle stress. + ComputeParticlePressure compute_s_func{ + _dt, _properties }; + _compute_stress.apply( "Picasso::particle_S", + Picasso::FieldLocation::Particle(), exec_space, + fm, particles, compute_s_func ); + + // Compute grid velocity change due to stress + ComputeGridVelocityChangeStress compute_du_s_func; + _compute_delta_u_s.apply( + "Picasso::div_S", Picasso::FieldLocation::Particle(), exec_space, + fm, particles, compute_du_s_func ); + + // Compute grid velocity change due to gravity + ComputeGridVelocityChangeGravity + compute_du_g_func{ _properties }; + _compute_delta_u_g.apply( + "Picasso::rhoG", Picasso::FieldLocation::Particle(), exec_space, fm, + particles, compute_du_g_func ); + + // Compute next grid velocity. + UpdateGridVelocity update_u_func{ _dt }; + _update_velocity.apply( "Picasso::update_U", + Picasso::FieldLocation::Node(), exec_space, fm, + update_u_func ); + + // Apply boundary condition. + ApplyBoundaryCondition apply_bc{ bc }; + _apply_bc_momentum.apply( "Picasso::BC_U", + Picasso::FieldLocation::Node(), exec_space, + fm, apply_bc ); + + // G2P + Grid2ParticleVelocity g2p_func{ _beta, + _dt }; + _g2p_momentum.apply( "Picasso::g2p_U", + Picasso::FieldLocation::Particle(), exec_space, fm, + particles, g2p_func ); + + // Do not force particles redistribution + particles.redistribute( local_grid, Field::Position() ); + + _total_time += _dt; + _total_steps += 1; + } + + double dt() { return _dt; } + double totalTime() { return _total_time; } + int totalSteps() { return _total_steps; } + + protected: + double _max_dt; + double _cfl_number; + double _beta; + bool _fixed_dt; + int _total_steps; + double _total_time; + double _dt; + double _cell_size; + + p2g_u_op _p2g_momentum; + compute_u_op _compute_velocity; + compute_s_op _compute_stress; + compute_du_s_op _compute_delta_u_s; + compute_du_g_op _compute_delta_u_g; + bc_u_op _apply_bc_momentum; + update_u_op _update_velocity; + g2p_u_op _g2p_momentum; + + property_type _properties; +}; + +template +auto createExplicitMomentumIntegrator( + const std::shared_ptr& mesh, InterpolationType, ParticleVelocity, + StressProperties properties, const double max_dt, const double cfl = 0.5, + const double beta = 1.0, const bool fixed_dt = false ) +{ + return ExplicitMomentumIntegrator( + mesh, properties, max_dt, cfl, beta, fixed_dt ); +} +//---------------------------------------------------------------------------// + +} // namespace Picasso diff --git a/examples/sources/Picasso_Output.hpp b/examples/sources/Picasso_Output.hpp new file mode 100644 index 00000000..10051a96 --- /dev/null +++ b/examples/sources/Picasso_Output.hpp @@ -0,0 +1,65 @@ +#ifndef PICASSO_OUTPUT_HPP +#define PICASSO_OUTPUT_HPP + +#include + +#include +#include + +namespace Picasso +{ +namespace Output +{ + +namespace Impl +{ + +Cabana::Experimental::HDF5ParticleOutput::HDF5Config setupHDF5Config() +{ + Cabana::Experimental::HDF5ParticleOutput::HDF5Config h5_config; + h5_config.collective = true; + return h5_config; +} + +} // namespace Impl +// Momentum-only output. +template +void outputParticles( MPI_Comm comm, ExecutionSpace, ParticleVelocity, + const int n, const int write_frequency, const double time, + const ParticleList particles ) +{ + if ( write_frequency > 0 && n % write_frequency == 0 ) + { + auto h5_config = Impl::setupHDF5Config(); + + auto x = particles.slice( Picasso::Field::Position() ); + auto p_p = particles.slice( Picasso::Field::Pressure() ); + auto vel_p = particles.slice( ParticleVelocity() ); + auto m_p = particles.slice( Picasso::Field::Mass() ); + auto v_p = particles.slice( Picasso::Field::Volume() ); + + std::string prefix = "particles"; + Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( + h5_config, prefix, comm, n / write_frequency, time, x.size(), x, + p_p, vel_p, m_p, v_p ); + } +} + +void outputTimestep( MPI_Comm comm, const int n, const int write_frequency, + const double total_time, const double dt ) +{ + if ( write_frequency > 0 && n % write_frequency == 0 ) + { + int rank; + MPI_Comm_rank( comm, &rank ); + if ( rank == 0 ) + std::cout << "Current step: " << n + << "; Current time: " << total_time + << "; Current dt: " << dt << std::endl; + } +} + +} // namespace Output +} // namespace Picasso + +#endif From 04792ba55d5e8ed2479a45098d2819d29d0dcc49 Mon Sep 17 00:00:00 2001 From: kwitaechong Date: Tue, 3 Sep 2024 15:53:25 -0400 Subject: [PATCH 02/22] change namespace Field to Example --- examples/dam_break.cpp | 6 +- examples/sources/Particle_Field.hpp | 4 +- examples/sources/Particle_Init.hpp | 12 +-- .../Picasso_ExplicitMomentumUpdate.hpp | 92 +++++++++---------- examples/sources/Picasso_Output.hpp | 8 +- 5 files changed, 61 insertions(+), 61 deletions(-) diff --git a/examples/dam_break.cpp b/examples/dam_break.cpp index dc011b85..7a47978c 100644 --- a/examples/dam_break.cpp +++ b/examples/dam_break.cpp @@ -45,9 +45,9 @@ void DamBreak() minimum_halo_size, MPI_COMM_WORLD ); // Make a particle list. - Cabana::ParticleTraits + Cabana::ParticleTraits fields; auto particles = Cabana::Grid::createParticleList( "test_particles", fields ); diff --git a/examples/sources/Particle_Field.hpp b/examples/sources/Particle_Field.hpp index 0e4c93b3..18685f78 100644 --- a/examples/sources/Particle_Field.hpp +++ b/examples/sources/Particle_Field.hpp @@ -10,7 +10,7 @@ namespace Picasso { -namespace Field +namespace Example { struct Mass : Field::Scalar @@ -60,7 +60,7 @@ struct DetDefGrad : Field::Scalar using Position = Picasso::Field::LogicalPosition<3>; -} // namespace Field +} // namespace Example namespace APIC { diff --git a/examples/sources/Particle_Init.hpp b/examples/sources/Particle_Init.hpp index 77bd6ac7..10c52fcf 100644 --- a/examples/sources/Particle_Init.hpp +++ b/examples/sources/Particle_Init.hpp @@ -33,15 +33,15 @@ struct ParticleInitFunc x[1] <= block[4] && block[2] <= x[2] && x[2] <= block[5] ) { - Picasso::get( p, Picasso::Field::Stress() ) = 0.0; + Picasso::get( p, Picasso::Example::Stress() ) = 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; + Picasso::get( p, Picasso::Example::DetDefGrad() ) = 1.0; + Picasso::get( p, Picasso::Example::Mass() ) = pv * density; + Picasso::get( p, Picasso::Example::Volume() ) = pv; + Picasso::get( p, Picasso::Example::Pressure() ) = 0.0; for ( int d = 0; d < 3; ++d ) - Picasso::get( p, Picasso::Field::Position(), d ) = x[d]; + Picasso::get( p, Picasso::Example::Position(), d ) = x[d]; return true; } diff --git a/examples/sources/Picasso_ExplicitMomentumUpdate.hpp b/examples/sources/Picasso_ExplicitMomentumUpdate.hpp index 911c2a6b..20ce0f64 100644 --- a/examples/sources/Picasso_ExplicitMomentumUpdate.hpp +++ b/examples/sources/Picasso_ExplicitMomentumUpdate.hpp @@ -26,11 +26,11 @@ struct ComputeGridVelocity { // Get the local dependencies. auto m_i = - local_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); - auto u_i = - local_deps.get( Picasso::FieldLocation::Node(), Field::Velocity() ); + local_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); + auto u_i = local_deps.get( Picasso::FieldLocation::Node(), + Example::Velocity() ); auto old_u_i = - local_deps.get( Picasso::FieldLocation::Node(), Field::OldU() ); + local_deps.get( Picasso::FieldLocation::Node(), Example::OldU() ); // Compute velocity. for ( int d = 0; d < 3; ++d ) @@ -64,16 +64,16 @@ struct ComputeParticlePressure ParticleViewType& particle ) const { // Get particle data. - auto x_p = Picasso::get( particle, Field::Position() ); - auto& J_p = Picasso::get( particle, Field::DetDefGrad() ); - auto& p_p = Picasso::get( particle, Field::Pressure() ); - auto s_p = Picasso::get( particle, Field::Stress() ); - auto v_p = Picasso::get( particle, Field::Volume() ); - auto m_p = Picasso::get( particle, Field::Mass() ); + auto x_p = Picasso::get( particle, Example::Position() ); + auto& J_p = Picasso::get( particle, Example::DetDefGrad() ); + auto& p_p = Picasso::get( particle, Example::Pressure() ); + auto s_p = Picasso::get( particle, Example::Stress() ); + auto v_p = Picasso::get( particle, Example::Volume() ); + auto m_p = Picasso::get( particle, Example::Mass() ); // Get the gather dependencies. auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Field::Velocity() ); + Example::Velocity() ); // update strain rate auto spline = Picasso::createSpline( @@ -114,13 +114,13 @@ struct ComputeGridVelocityChangeStress const LocalDependencies&, ParticleViewType& particle ) const { // Get particle data. - auto v_p = Picasso::get( particle, Field::Volume() ); - auto x_p = Picasso::get( particle, Field::Position() ); - auto s_p = Picasso::get( particle, Field::Stress() ); + auto v_p = Picasso::get( particle, Example::Volume() ); + auto x_p = Picasso::get( particle, Example::Position() ); + auto s_p = Picasso::get( particle, Example::Stress() ); // Get the scatter dependencies. auto delta_u_s_i = scatter_deps.get( Picasso::FieldLocation::Node(), - Field::DeltaUStress() ); + Example::DeltaUStress() ); // Node interpolant. auto spline = Picasso::createSpline( @@ -151,12 +151,12 @@ struct ComputeGridVelocityChangeGravity const LocalDependencies&, ParticleViewType& particle ) const { // Get particle data. - auto m_p = Picasso::get( particle, Field::Mass() ); - auto x_p = Picasso::get( particle, Field::Position() ); + auto m_p = Picasso::get( particle, Example::Mass() ); + auto x_p = Picasso::get( particle, Example::Position() ); // Get the scatter dependencies. auto delta_u_g_i = scatter_deps.get( Picasso::FieldLocation::Node(), - Field::DeltaUGravity() ); + Example::DeltaUGravity() ); // Node interpolant. auto spline = Picasso::createSpline( @@ -190,14 +190,14 @@ struct UpdateGridVelocity { // Get the local dependencies. auto m_i = - gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + gather_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Field::Velocity() ); + Example::Velocity() ); auto delta_u_s_i = gather_deps.get( Picasso::FieldLocation::Node(), - Field::DeltaUStress() ); + Example::DeltaUStress() ); auto delta_u_g_i = gather_deps.get( Picasso::FieldLocation::Node(), - Field::DeltaUGravity() ); + Example::DeltaUGravity() ); // Compute velocity. Picasso::Vec3 zeros = { 0.0, 0.0, 0.0 }; @@ -239,13 +239,13 @@ struct Grid2ParticleVelocity { // Get particle data. auto f_p = Picasso::get( particle, APIC::Field::Velocity() ); - auto x_p = Picasso::get( particle, Field::Position() ); + auto x_p = Picasso::get( particle, Example::Position() ); // Get the gather dependencies. auto m_i = - gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + gather_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Field::Velocity() ); + Example::Velocity() ); // Get the local dependencies for getting physcial location of node auto x_i = local_deps.get( Picasso::FieldLocation::Node(), Picasso::Field::PhysicalPosition<3>() ); @@ -257,7 +257,7 @@ struct Grid2ParticleVelocity Picasso::SplineValue(), Picasso::SplineDistance(), Picasso::SplineGradient() ); - // Update particle velocity using a PolyPIC update. + // Update particle velocity using a APIC update. Picasso::APIC::g2p( u_i, f_p, spline ); // Update particle position. @@ -293,16 +293,16 @@ struct Grid2ParticleVelocity ParticleViewType& particle ) const { // Get particle data. - auto u_p = Picasso::get( particle, Field::Velocity() ); - auto x_p = Picasso::get( particle, Field::Position() ); + auto u_p = Picasso::get( particle, Example::Velocity() ); + auto x_p = Picasso::get( particle, Example::Position() ); // Get the gather dependencies. auto m_i = - gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + gather_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Field::Velocity() ); + Example::Velocity() ); auto old_u_i = - gather_deps.get( Picasso::FieldLocation::Node(), Field::OldU() ); + gather_deps.get( Picasso::FieldLocation::Node(), Example::OldU() ); // Get the local dependencies for getting physcial location of node auto x_i = local_deps.get( Picasso::FieldLocation::Node(), @@ -394,12 +394,12 @@ struct Particle2Grid; + Picasso::FieldLayout; using velocity_type = - Picasso::FieldLayout; + Picasso::FieldLayout; using old_u_type = - Picasso::FieldLayout; - using delta_u_s_type = - Picasso::FieldLayout; + Picasso::FieldLayout; + using delta_u_s_type = Picasso::FieldLayout; using delta_u_g_type = Picasso::FieldLayout; + Example::DeltaUGravity>; using grid_type = Picasso::FieldLayout>; @@ -580,7 +580,7 @@ class ExplicitMomentumIntegrator const int spline_order = 1; // P2G - Particle2Grid p2g_func{ _dt }; _p2g_momentum.apply( "Picasso::p2g_U", @@ -620,7 +620,7 @@ class ExplicitMomentumIntegrator update_u_func ); // Apply boundary condition. - ApplyBoundaryCondition apply_bc{ bc }; + ApplyBoundaryCondition apply_bc{ bc }; _apply_bc_momentum.apply( "Picasso::BC_U", Picasso::FieldLocation::Node(), exec_space, fm, apply_bc ); @@ -633,7 +633,7 @@ class ExplicitMomentumIntegrator particles, g2p_func ); // Do not force particles redistribution - particles.redistribute( local_grid, Field::Position() ); + particles.redistribute( local_grid, Example::Position() ); _total_time += _dt; _total_steps += 1; diff --git a/examples/sources/Picasso_Output.hpp b/examples/sources/Picasso_Output.hpp index 10051a96..4129901c 100644 --- a/examples/sources/Picasso_Output.hpp +++ b/examples/sources/Picasso_Output.hpp @@ -32,11 +32,11 @@ void outputParticles( MPI_Comm comm, ExecutionSpace, ParticleVelocity, { auto h5_config = Impl::setupHDF5Config(); - auto x = particles.slice( Picasso::Field::Position() ); - auto p_p = particles.slice( Picasso::Field::Pressure() ); + auto x = particles.slice( Picasso::Example::Position() ); + auto p_p = particles.slice( Picasso::Example::Pressure() ); auto vel_p = particles.slice( ParticleVelocity() ); - auto m_p = particles.slice( Picasso::Field::Mass() ); - auto v_p = particles.slice( Picasso::Field::Volume() ); + auto m_p = particles.slice( Picasso::Example::Mass() ); + auto v_p = particles.slice( Picasso::Example::Volume() ); std::string prefix = "particles"; Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( From 126ed550c46df4d589e8d9e0114ec32baa967fcb Mon Sep 17 00:00:00 2001 From: kwitaechong Date: Thu, 5 Sep 2024 10:50:12 -0400 Subject: [PATCH 03/22] use InputParser with dam_break.json --- CMakeLists.txt | 5 ++- examples/CMakeLists.txt | 4 +- examples/dam_break.cpp | 56 ++++++++++++++----------- examples/inputs/dam_break.json | 20 +++++++++ examples/inputs/particle_init_test.json | 10 ----- 5 files changed, 58 insertions(+), 37 deletions(-) create mode 100644 examples/inputs/dam_break.json delete mode 100644 examples/inputs/particle_init_test.json diff --git a/CMakeLists.txt b/CMakeLists.txt index 267dfe18..37c99eac 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" diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 9795e0e6..cb970489 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,8 +1,8 @@ add_subdirectory(sources) configure_file( - ${CMAKE_CURRENT_SOURCE_DIR}/inputs/particle_init_test.json - ${CMAKE_CURRENT_BINARY_DIR}/particle_init_test.json + ${CMAKE_CURRENT_SOURCE_DIR}/inputs/dam_break.json + ${CMAKE_CURRENT_BINARY_DIR}/dam_break.json COPYONLY) add_executable( DamBreak dam_break.cpp ) diff --git a/examples/dam_break.cpp b/examples/dam_break.cpp index 7a47978c..53d915c3 100644 --- a/examples/dam_break.cpp +++ b/examples/dam_break.cpp @@ -14,6 +14,21 @@ using namespace Picasso; +template +Kokkos::Array parserGetArray( InputParser parser, + const std::string name ) +{ + const auto& pt = parser.propertyTree(); + Kokkos::Array array; + + int d = 0; + for ( auto& element : pt.get_child( name ) ) + { + array[d] = element.second.get_value(); + ++d; + } + return array; +} //---------------------------------------------------------------------------// // DamBreak example template @@ -22,22 +37,12 @@ void DamBreak() using exec_space = Kokkos::DefaultExecutionSpace; using memory_space = exec_space::memory_space; - // Global bounding box. - double cell_size = 0.05; - std::array global_num_cell = { 20, 20, 20 }; - std::array global_low_corner = { 0.0, 0.0, 0.0 }; - std::array global_high_corner = { - global_low_corner[0] + cell_size * global_num_cell[0], - global_low_corner[1] + cell_size * global_num_cell[1], - global_low_corner[2] + cell_size * global_num_cell[2] }; - // Get inputs for mesh. InputParser parser( "dam_break.json", "json" ); const auto& pt = parser.propertyTree(); - Kokkos::Array global_box = { - global_low_corner[0], global_low_corner[1], global_low_corner[2], - global_high_corner[0], global_high_corner[1], global_high_corner[2] }; + // Global bounding box. + auto global_box = parserGetArray( parser, "global_box" ); int minimum_halo_size = 0; // Make mesh. @@ -53,13 +58,13 @@ void DamBreak() "test_particles", fields ); // Initialize particles - Kokkos::Array block = { 0.0, 0.0, 0.0, 0.4, 0.4, 0.4 }; - double density = 1e3; + auto particle_box = parserGetArray( parser, "particle_box" ); + auto density = pt.get( "density" ); - auto momentum_init_functor = - createParticleInitFunc( particles, ParticleVelocity(), block, density ); + auto momentum_init_functor = createParticleInitFunc( + particles, ParticleVelocity(), particle_box, density ); - int ppc = 2; + auto ppc = pt.get( "ppc" ); auto local_grid = mesh->localGrid(); Cabana::Grid::createParticles( Cabana::InitUniform(), exec_space(), momentum_init_functor, particles, ppc, @@ -73,25 +78,28 @@ void DamBreak() BoundaryCondition bc{ bc_index }; // Properties - double gamma = 7.0; - double bulk_modulus = 1.0e+5; - Kokkos::Array gravity = { 0.0, 0.0, -9.8 }; + auto gamma = pt.get( "gamma" ); + auto bulk_modulus = pt.get( "bulk_modulus" ); + auto gravity = parserGetArray( parser, "gravity" ); Picasso::Properties props( gamma, bulk_modulus, gravity ); // Time integragor + auto dt = pt.get( "dt" ); auto time_integrator = Picasso::createExplicitMomentumIntegrator( - mesh, InterpolationType(), ParticleVelocity(), props, 1.0e-4 ); + mesh, InterpolationType(), ParticleVelocity(), props, dt ); auto fm = Picasso::createFieldManager( mesh ); time_integrator.setup( *fm ); // steps - while ( time_integrator.totalTime() < 1.0 ) + auto final_time = pt.get( "final_time" ); + auto write_frequency = pt.get( "write_frequency" ); + while ( time_integrator.totalTime() < final_time ) { // Write particle fields. Picasso::Output::outputParticles( MPI_COMM_WORLD, exec_space(), ParticleVelocity(), - time_integrator.totalSteps(), 10, time_integrator.totalTime(), - particles ); + time_integrator.totalSteps(), write_frequency, + time_integrator.totalTime(), particles ); printf( "aaa\n" ); diff --git a/examples/inputs/dam_break.json b/examples/inputs/dam_break.json new file mode 100644 index 00000000..b83847ce --- /dev/null +++ b/examples/inputs/dam_break.json @@ -0,0 +1,20 @@ +{ + "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, + "final_time": 1.0, + "write_frequency" : 10 +} diff --git a/examples/inputs/particle_init_test.json b/examples/inputs/particle_init_test.json deleted file mode 100644 index a21c2c34..00000000 --- a/examples/inputs/particle_init_test.json +++ /dev/null @@ -1,10 +0,0 @@ -{ - "mesh": { - "cell_size": 0.2, - "periodic": [false, false, false], - "partitioner": { - "type": "uniform_dim" - }, - "halo_cell_width": 0 - } -} From 7c32948b48383fba95024572693ccddd66dd4c4e Mon Sep 17 00:00:00 2001 From: kwitaechong Date: Thu, 5 Sep 2024 12:23:41 -0400 Subject: [PATCH 04/22] rebase master and use nlohmann::json --- examples/dam_break.cpp | 43 ++++++++++++++++--------------------- src/Picasso_InputParser.hpp | 11 ++++++++++ 2 files changed, 29 insertions(+), 25 deletions(-) diff --git a/examples/dam_break.cpp b/examples/dam_break.cpp index 53d915c3..3ed480e0 100644 --- a/examples/dam_break.cpp +++ b/examples/dam_break.cpp @@ -14,20 +14,14 @@ using namespace Picasso; -template -Kokkos::Array parserGetArray( InputParser parser, - const std::string name ) +template +Kokkos::Array getArray( nlohmann::json input_data, + const std::string name ) { - const auto& pt = parser.propertyTree(); - Kokkos::Array array; + std::array property_std = input_data[name]; + auto property = copy( property_std ); - int d = 0; - for ( auto& element : pt.get_child( name ) ) - { - array[d] = element.second.get_value(); - ++d; - } - return array; + return property; } //---------------------------------------------------------------------------// // DamBreak example @@ -38,15 +32,14 @@ void DamBreak() using memory_space = exec_space::memory_space; // Get inputs for mesh. - InputParser parser( "dam_break.json", "json" ); - const auto& pt = parser.propertyTree(); + auto inputs = Picasso::parse( "dam_break.json" ); // Global bounding box. - auto global_box = parserGetArray( parser, "global_box" ); + auto global_box = getArray( inputs, "global_box" ); int minimum_halo_size = 0; // Make mesh. - auto mesh = createUniformMesh( memory_space(), pt, global_box, + auto mesh = createUniformMesh( memory_space(), inputs, global_box, minimum_halo_size, MPI_COMM_WORLD ); // Make a particle list. @@ -58,13 +51,13 @@ void DamBreak() "test_particles", fields ); // Initialize particles - auto particle_box = parserGetArray( parser, "particle_box" ); - auto density = pt.get( "density" ); + auto particle_box = getArray( inputs, "particle_box" ); + double density = inputs["density"]; auto momentum_init_functor = createParticleInitFunc( particles, ParticleVelocity(), particle_box, density ); - auto ppc = pt.get( "ppc" ); + double ppc = inputs["ppc"]; auto local_grid = mesh->localGrid(); Cabana::Grid::createParticles( Cabana::InitUniform(), exec_space(), momentum_init_functor, particles, ppc, @@ -78,21 +71,21 @@ void DamBreak() BoundaryCondition bc{ bc_index }; // Properties - auto gamma = pt.get( "gamma" ); - auto bulk_modulus = pt.get( "bulk_modulus" ); - auto gravity = parserGetArray( parser, "gravity" ); + auto gamma = inputs["gamma"]; + auto bulk_modulus = inputs["bulk_modulus"]; + auto gravity = getArray( inputs, "gravity" ); Picasso::Properties props( gamma, bulk_modulus, gravity ); // Time integragor - auto dt = pt.get( "dt" ); + auto dt = inputs["dt"]; auto time_integrator = Picasso::createExplicitMomentumIntegrator( mesh, InterpolationType(), ParticleVelocity(), props, dt ); auto fm = Picasso::createFieldManager( mesh ); time_integrator.setup( *fm ); // steps - auto final_time = pt.get( "final_time" ); - auto write_frequency = pt.get( "write_frequency" ); + auto final_time = inputs["final_time"]; + auto write_frequency = inputs["write_frequency"]; while ( time_integrator.totalTime() < final_time ) { // Write particle fields. 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 #include +#include + namespace Picasso { +template +Kokkos::Array copy( const std::array in ) +{ + Kokkos::Array 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 ); From 3ee18a97574f05f2875c4a1efaca9827caf24eb4 Mon Sep 17 00:00:00 2001 From: kwitaechong Date: Wed, 11 Sep 2024 13:26:22 -0400 Subject: [PATCH 05/22] move p2g/g2p inside src/Picasso_ParticleInterpolation.hpp --- examples/dam_break.cpp | 2 +- examples/sources/Example_Sources.hpp | 1 - examples/sources/Particle_Field.hpp | 87 ------- examples/sources/Particle_Init.hpp | 2 - .../Picasso_ExplicitMomentumUpdate.hpp | 221 ------------------ src/Picasso_FieldTypes.hpp | 75 ++++++ src/Picasso_ParticleInterpolation.hpp | 221 ++++++++++++++++++ 7 files changed, 297 insertions(+), 312 deletions(-) delete mode 100644 examples/sources/Particle_Field.hpp diff --git a/examples/dam_break.cpp b/examples/dam_break.cpp index 3ed480e0..e5b015e2 100644 --- a/examples/dam_break.cpp +++ b/examples/dam_break.cpp @@ -114,7 +114,7 @@ int main( int argc, char* argv[] ) $/: ./DamBreak inputs/dam_break.json\n" ); std::string filename = argv[1]; - DamBreak(); + DamBreak(); Kokkos::finalize(); MPI_Finalize(); diff --git a/examples/sources/Example_Sources.hpp b/examples/sources/Example_Sources.hpp index bdcadbec..bcf38d00 100644 --- a/examples/sources/Example_Sources.hpp +++ b/examples/sources/Example_Sources.hpp @@ -1,4 +1,3 @@ -#include "Particle_Field.hpp" #include "Particle_Init.hpp" #include "Picasso_BoundaryCondition.hpp" #include "Picasso_ExplicitMomentumUpdate.hpp" diff --git a/examples/sources/Particle_Field.hpp b/examples/sources/Particle_Field.hpp deleted file mode 100644 index 18685f78..00000000 --- a/examples/sources/Particle_Field.hpp +++ /dev/null @@ -1,87 +0,0 @@ -#include - -#include -#include - -#include - -#include - -namespace Picasso -{ - -namespace Example -{ - -struct Mass : Field::Scalar -{ - static std::string label() { return "Mass"; } -}; - -struct Pressure : Field::Scalar -{ - static std::string label() { return "Pressure"; } -}; - -struct Volume : Field::Scalar -{ - static std::string label() { return "Volume"; } -}; - -struct Velocity : Field::Vector -{ - static std::string label() { return "Velocity"; } -}; - -struct OldU : Field::Vector -{ - static std::string label() { return "Old_Velocity"; } -}; - -struct DeltaUGravity : Field::Vector -{ - static std::string label() { return "velocity_change_from_gravity"; } -}; - -struct Stress : Field::Matrix -{ - static std::string label() { return "stress"; } -}; - -struct DeltaUStress : Field::Vector -{ - static std::string label() { return "velocity_change_from_stress"; } -}; - -struct DetDefGrad : Field::Scalar -{ - static std::string label() { return "Det_deformation_gradient"; } -}; - -using Position = Picasso::Field::LogicalPosition<3>; - -} // namespace Example - -namespace APIC -{ - -struct APicTag -{ -}; - -namespace Field -{ - -struct Velocity : Picasso::Field::Matrix -{ - static std::string label() { return "velocity"; } -}; - -} // namespace Field -} // namespace APIC - -struct FlipTag -{ -}; - -} // namespace Picasso diff --git a/examples/sources/Particle_Init.hpp b/examples/sources/Particle_Init.hpp index 10c52fcf..2a5818d3 100644 --- a/examples/sources/Particle_Init.hpp +++ b/examples/sources/Particle_Init.hpp @@ -1,7 +1,5 @@ #include -#include "Particle_Field.hpp" - #include #include diff --git a/examples/sources/Picasso_ExplicitMomentumUpdate.hpp b/examples/sources/Picasso_ExplicitMomentumUpdate.hpp index 20ce0f64..67c9c82e 100644 --- a/examples/sources/Picasso_ExplicitMomentumUpdate.hpp +++ b/examples/sources/Picasso_ExplicitMomentumUpdate.hpp @@ -213,140 +213,6 @@ struct UpdateGridVelocity } }; -template -struct Grid2ParticleVelocity; - -//---------------------------------------------------------------------------// -// Update particle state. APIC variant. -//---------------------------------------------------------------------------// -template -struct Grid2ParticleVelocity -{ - // FIXME: only needed for FLIP/PIC - double beta; - - // Explicit time step. - double dt; - - template - 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() ); - auto x_p = Picasso::get( particle, Example::Position() ); - - // Get the gather dependencies. - auto m_i = - gather_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); - auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Example::Velocity() ); - // Get the local dependencies for getting physcial location of node - auto x_i = local_deps.get( Picasso::FieldLocation::Node(), - Picasso::Field::PhysicalPosition<3>() ); - - // Node interpolant. - auto spline = Picasso::createSpline( - Picasso::FieldLocation::Node(), - Picasso::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 ) - { - printf( "g2p %d %d %d %8.5e %8.5e\n", i, j, k, x_i( i, j, k, 2 ), - u_i( i, j, k, 2 ) ); - 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 -struct Grid2ParticleVelocity -{ - double beta = 0.01; - - // Explicit time step. - double dt; - - template - 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, Example::Velocity() ); - auto x_p = Picasso::get( particle, Example::Position() ); - - // Get the gather dependencies. - auto m_i = - gather_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); - auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Example::Velocity() ); - auto old_u_i = - gather_deps.get( Picasso::FieldLocation::Node(), Example::OldU() ); - - // Get the local dependencies for getting physcial location of node - auto x_i = local_deps.get( Picasso::FieldLocation::Node(), - Picasso::Field::PhysicalPosition<3>() ); - - // Node interpolant. - auto spline = Picasso::createSpline( - Picasso::FieldLocation::Node(), - Picasso::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 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 d_u_p; - Picasso::G2P::value( spline, d_u_i, d_u_p ); - - Picasso::Vec3 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 ); - } -}; - //---------------------------------------------------------------------------// // Compute boundary condition on grid. //---------------------------------------------------------------------------// @@ -370,93 +236,6 @@ struct ApplyBoundaryCondition } }; -template -struct Particle2Grid; - -//---------------------------------------------------------------------------// -// Project particle enthalpy/momentum to grid. APIC variant -//---------------------------------------------------------------------------// -template -struct Particle2Grid -{ - // Explicit time step. - double dt; - - template - 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, Example::Mass() ); - auto x_p = Picasso::get( particle, Example::Position() ); - - // Get the scatter dependencies. - auto m_i = - scatter_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); - auto f_i = - scatter_deps.get( Picasso::FieldLocation::Node(), OldFieldType() ); - - // Node interpolant. - auto spline = Picasso::createSpline( - Picasso::FieldLocation::Node(), - Picasso::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 -struct Particle2Grid -{ - // Explicit time step. - double dt; - - template - 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, Example::Mass() ); - auto x_p = Picasso::get( particle, Example::Position() ); - - // Get the scatter dependencies. - auto m_i = - scatter_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); - auto f_i = - scatter_deps.get( Picasso::FieldLocation::Node(), OldFieldType() ); - - // Node interpolant. - auto spline = Picasso::createSpline( - Picasso::FieldLocation::Node(), - Picasso::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 ); - } -}; - //---------------------------------------------------------------------------// // Explicit momentum time stepper //---------------------------------------------------------------------------// diff --git a/src/Picasso_FieldTypes.hpp b/src/Picasso_FieldTypes.hpp index 11bc6c01..465fdc52 100644 --- a/src/Picasso_FieldTypes.hpp +++ b/src/Picasso_FieldTypes.hpp @@ -969,6 +969,81 @@ struct CommRank : Scalar //---------------------------------------------------------------------------// } // end namespace Field + +namespace Example +{ + +struct Mass : Field::Scalar +{ + static std::string label() { return "Mass"; } +}; + +struct Pressure : Field::Scalar +{ + static std::string label() { return "Pressure"; } +}; + +struct Volume : Field::Scalar +{ + static std::string label() { return "Volume"; } +}; + +struct Velocity : Field::Vector +{ + static std::string label() { return "Velocity"; } +}; + +struct OldU : Field::Vector +{ + static std::string label() { return "Old_Velocity"; } +}; + +struct DeltaUGravity : Field::Vector +{ + static std::string label() { return "velocity_change_from_gravity"; } +}; + +struct Stress : Field::Matrix +{ + static std::string label() { return "stress"; } +}; + +struct DeltaUStress : Field::Vector +{ + static std::string label() { return "velocity_change_from_stress"; } +}; + +struct DetDefGrad : Field::Scalar +{ + static std::string label() { return "Det_deformation_gradient"; } +}; + +using Position = Picasso::Field::LogicalPosition<3>; + +} // namespace Example + +namespace APIC +{ + +namespace Field +{ + +struct Velocity : Picasso::Field::Matrix +{ + static std::string label() { return "velocity"; } +}; + +} // namespace Field +} // namespace APIC + +struct APicTag +{ +}; + +struct FlipTag +{ +}; + } // end namespace Picasso #endif // PICASSO_FIELDTYPES_HPP diff --git a/src/Picasso_ParticleInterpolation.hpp b/src/Picasso_ParticleInterpolation.hpp index fa41a0c6..88aa38fe 100644 --- a/src/Picasso_ParticleInterpolation.hpp +++ b/src/Picasso_ParticleInterpolation.hpp @@ -12,6 +12,7 @@ #ifndef PICASSO_PARTICLEINTERPOLATION_HPP #define PICASSO_PARTICLEINTERPOLATION_HPP +#include #include #include @@ -279,6 +280,226 @@ KOKKOS_INLINE_FUNCTION void divergence( const SplineDataType& sd, } // end namespace P2G //---------------------------------------------------------------------------// +template +struct Particle2Grid; + +//---------------------------------------------------------------------------// +// Project particle enthalpy/momentum to grid. APIC variant +//---------------------------------------------------------------------------// +template +struct Particle2Grid +{ + // Explicit time step. + double dt; + + template + 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, Example::Mass() ); + auto x_p = Picasso::get( particle, Example::Position() ); + + // Get the scatter dependencies. + auto m_i = + scatter_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); + auto f_i = + scatter_deps.get( Picasso::FieldLocation::Node(), OldFieldType() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::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 +struct Particle2Grid +{ + // Explicit time step. + double dt; + + template + 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, Example::Mass() ); + auto x_p = Picasso::get( particle, Example::Position() ); + + // Get the scatter dependencies. + auto m_i = + scatter_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); + auto f_i = + scatter_deps.get( Picasso::FieldLocation::Node(), OldFieldType() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::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 ); + } +}; + +template +struct Grid2ParticleVelocity; + +//---------------------------------------------------------------------------// +// Update particle state. APIC variant. +//---------------------------------------------------------------------------// +template +struct Grid2ParticleVelocity +{ + // FIXME: only needed for FLIP/PIC + double beta; + + // Explicit time step. + double dt; + + template + 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() ); + auto x_p = Picasso::get( particle, Example::Position() ); + + // Get the gather dependencies. + auto m_i = + gather_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); + auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), + Example::Velocity() ); + // Get the local dependencies for getting physcial location of node + auto x_i = local_deps.get( Picasso::FieldLocation::Node(), + Picasso::Field::PhysicalPosition<3>() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::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 ) + { + printf( "g2p %d %d %d %8.5e %8.5e\n", i, j, k, x_i( i, j, k, 2 ), + u_i( i, j, k, 2 ) ); + 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 +struct Grid2ParticleVelocity +{ + double beta = 0.01; + + // Explicit time step. + double dt; + + template + 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, Example::Velocity() ); + auto x_p = Picasso::get( particle, Example::Position() ); + + // Get the gather dependencies. + auto m_i = + gather_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); + auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), + Example::Velocity() ); + auto old_u_i = + gather_deps.get( Picasso::FieldLocation::Node(), Example::OldU() ); + + // Get the local dependencies for getting physcial location of node + auto x_i = local_deps.get( Picasso::FieldLocation::Node(), + Picasso::Field::PhysicalPosition<3>() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::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 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 d_u_p; + Picasso::G2P::value( spline, d_u_i, d_u_p ); + + Picasso::Vec3 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 ); + } +}; } // end namespace Picasso From 4566eff61929c7676410d92f33af7689a4e343c4 Mon Sep 17 00:00:00 2001 From: kwitaechong Date: Wed, 11 Sep 2024 14:03:40 -0400 Subject: [PATCH 06/22] add PolyPIC --- src/Picasso_FieldTypes.hpp | 16 +++++ src/Picasso_ParticleInterpolation.hpp | 98 +++++++++++++++++++++++++++ 2 files changed, 114 insertions(+) diff --git a/src/Picasso_FieldTypes.hpp b/src/Picasso_FieldTypes.hpp index 465fdc52..00011c27 100644 --- a/src/Picasso_FieldTypes.hpp +++ b/src/Picasso_FieldTypes.hpp @@ -1022,6 +1022,18 @@ using Position = Picasso::Field::LogicalPosition<3>; } // namespace Example +namespace PolyPIC +{ +namespace Field +{ + +struct Velocity : Picasso::Field::Matrix +{ + static std::string label() { return "velocity"; } +}; +} // namespace Field +} // namespace PolyPIC + namespace APIC { @@ -1036,6 +1048,10 @@ struct Velocity : Picasso::Field::Matrix } // namespace Field } // namespace APIC +struct PolyPicTag +{ +}; + struct APicTag { }; diff --git a/src/Picasso_ParticleInterpolation.hpp b/src/Picasso_ParticleInterpolation.hpp index 88aa38fe..73384865 100644 --- a/src/Picasso_ParticleInterpolation.hpp +++ b/src/Picasso_ParticleInterpolation.hpp @@ -15,6 +15,7 @@ #include #include #include +#include #include @@ -284,6 +285,48 @@ template struct Particle2Grid; +//---------------------------------------------------------------------------// +// Project particle enthalpy/momentum to grid. PolyPIC variant +//---------------------------------------------------------------------------// +template +struct Particle2Grid +{ + // Explicit time step. + double dt; + + template + 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 v_p = Picasso::get( particle, PolyPIC::Field::Velocity() ); + auto m_p = Picasso::get( particle, Example::Mass() ); + auto x_p = Picasso::get( particle, Example::Position() ); + + // Get the scatter dependencies. + auto m_i = + scatter_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); + auto f_i = + scatter_deps.get( Picasso::FieldLocation::Node(), OldFieldType() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::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 //---------------------------------------------------------------------------// @@ -370,6 +413,61 @@ struct Particle2Grid struct Grid2ParticleVelocity; +//---------------------------------------------------------------------------// +// Update particle state. PolyPIC variant. +//---------------------------------------------------------------------------// +template +struct Grid2ParticleVelocity +{ + // FIXME: only needed for FLIP/PIC + double beta; + + // Explicit time step. + double dt; + + template + 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() ); + auto x_p = Picasso::get( particle, Example::Position() ); + + // Get the gather dependencies. + auto m_i = + gather_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); + auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), + Example::Velocity() ); + // Get the local dependencies for getting physcial location of node + auto x_i = local_deps.get( Picasso::FieldLocation::Node(), + Picasso::Field::PhysicalPosition<3>() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::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 ) + { + printf( "g2p %d %d %d %8.5e %8.5e\n", i, j, k, x_i( i, j, k, 2 ), + u_i( i, j, k, 2 ) ); + 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. //---------------------------------------------------------------------------// From 58c194003f7122a5bdbfad0a4eff2caaca36bf04 Mon Sep 17 00:00:00 2001 From: kwitaechong Date: Wed, 11 Sep 2024 23:00:26 -0400 Subject: [PATCH 07/22] clean up code --- examples/dam_break.cpp | 34 ++++---- examples/sources/Particle_Init.hpp | 12 +-- .../Picasso_ExplicitMomentumUpdate.hpp | 80 +++++++++---------- examples/sources/Picasso_Output.hpp | 65 --------------- src/Picasso_FieldTypes.hpp | 5 -- src/Picasso_ParticleInterpolation.hpp | 52 ++++++------ 6 files changed, 80 insertions(+), 168 deletions(-) delete mode 100644 examples/sources/Picasso_Output.hpp diff --git a/examples/dam_break.cpp b/examples/dam_break.cpp index e5b015e2..5830438c 100644 --- a/examples/dam_break.cpp +++ b/examples/dam_break.cpp @@ -3,7 +3,6 @@ #include "sources/Particle_Init.hpp" #include "sources/Picasso_BoundaryCondition.hpp" #include "sources/Picasso_ExplicitMomentumUpdate.hpp" -#include "sources/Picasso_Output.hpp" #include #include @@ -14,15 +13,6 @@ using namespace Picasso; -template -Kokkos::Array getArray( nlohmann::json input_data, - const std::string name ) -{ - std::array property_std = input_data[name]; - auto property = copy( property_std ); - - return property; -} //---------------------------------------------------------------------------// // DamBreak example template @@ -35,7 +25,7 @@ void DamBreak() auto inputs = Picasso::parse( "dam_break.json" ); // Global bounding box. - auto global_box = getArray( inputs, "global_box" ); + auto global_box = copy( inputs["global_box"] ); int minimum_halo_size = 0; // Make mesh. @@ -43,15 +33,14 @@ void DamBreak() minimum_halo_size, MPI_COMM_WORLD ); // Make a particle list. - Cabana::ParticleTraits + Cabana::ParticleTraits fields; auto particles = Cabana::Grid::createParticleList( "test_particles", fields ); // Initialize particles - auto particle_box = getArray( inputs, "particle_box" ); + auto particle_box = copy( inputs["particle_box"] ); double density = inputs["density"]; auto momentum_init_functor = createParticleInitFunc( @@ -73,7 +62,7 @@ void DamBreak() // Properties auto gamma = inputs["gamma"]; auto bulk_modulus = inputs["bulk_modulus"]; - auto gravity = getArray( inputs, "gravity" ); + auto gravity = copy( inputs["gravity"] ); Picasso::Properties props( gamma, bulk_modulus, gravity ); // Time integragor @@ -89,10 +78,15 @@ void DamBreak() while ( time_integrator.totalTime() < final_time ) { // Write particle fields. - Picasso::Output::outputParticles( - MPI_COMM_WORLD, exec_space(), ParticleVelocity(), - time_integrator.totalSteps(), write_frequency, - time_integrator.totalTime(), particles ); + Cabana::Experimental::HDF5ParticleOutput::HDF5Config h5_config; + Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( + h5_config, "particles", MPI_COMM_WORLD, + time_integrator.totalSteps(), time_integrator.totalTime(), + particles.size(), particles.slice( Picasso::Position() ), + particles.slice( Picasso::Pressure() ), + particles.slice( ParticleVelocity() ), + particles.slice( Picasso::Mass() ), + particles.slice( Picasso::Volume() ) ); printf( "aaa\n" ); diff --git a/examples/sources/Particle_Init.hpp b/examples/sources/Particle_Init.hpp index 2a5818d3..c6bd1178 100644 --- a/examples/sources/Particle_Init.hpp +++ b/examples/sources/Particle_Init.hpp @@ -31,15 +31,15 @@ struct ParticleInitFunc x[1] <= block[4] && block[2] <= x[2] && x[2] <= block[5] ) { - Picasso::get( p, Picasso::Example::Stress() ) = 0.0; + Picasso::get( p, Picasso::Stress() ) = 0.0; Picasso::get( p, VelocityType() ) = 0.0; - Picasso::get( p, Picasso::Example::DetDefGrad() ) = 1.0; - Picasso::get( p, Picasso::Example::Mass() ) = pv * density; - Picasso::get( p, Picasso::Example::Volume() ) = pv; - Picasso::get( p, Picasso::Example::Pressure() ) = 0.0; + Picasso::get( p, Picasso::DetDefGrad() ) = 1.0; + Picasso::get( p, Picasso::Mass() ) = pv * density; + Picasso::get( p, Picasso::Volume() ) = pv; + Picasso::get( p, Picasso::Pressure() ) = 0.0; for ( int d = 0; d < 3; ++d ) - Picasso::get( p, Picasso::Example::Position(), d ) = x[d]; + Picasso::get( p, Picasso::Position(), d ) = x[d]; return true; } diff --git a/examples/sources/Picasso_ExplicitMomentumUpdate.hpp b/examples/sources/Picasso_ExplicitMomentumUpdate.hpp index 67c9c82e..e895e371 100644 --- a/examples/sources/Picasso_ExplicitMomentumUpdate.hpp +++ b/examples/sources/Picasso_ExplicitMomentumUpdate.hpp @@ -25,12 +25,9 @@ struct ComputeGridVelocity const int i, const int j, const int k ) const { // Get the local dependencies. - auto m_i = - local_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); - auto u_i = local_deps.get( Picasso::FieldLocation::Node(), - Example::Velocity() ); - auto old_u_i = - local_deps.get( Picasso::FieldLocation::Node(), Example::OldU() ); + auto m_i = local_deps.get( Picasso::FieldLocation::Node(), Mass() ); + auto u_i = local_deps.get( Picasso::FieldLocation::Node(), Velocity() ); + auto old_u_i = local_deps.get( Picasso::FieldLocation::Node(), OldU() ); // Compute velocity. for ( int d = 0; d < 3; ++d ) @@ -64,16 +61,16 @@ struct ComputeParticlePressure ParticleViewType& particle ) const { // Get particle data. - auto x_p = Picasso::get( particle, Example::Position() ); - auto& J_p = Picasso::get( particle, Example::DetDefGrad() ); - auto& p_p = Picasso::get( particle, Example::Pressure() ); - auto s_p = Picasso::get( particle, Example::Stress() ); - auto v_p = Picasso::get( particle, Example::Volume() ); - auto m_p = Picasso::get( particle, Example::Mass() ); + auto x_p = Picasso::get( particle, Position() ); + auto& J_p = Picasso::get( particle, DetDefGrad() ); + auto& p_p = Picasso::get( particle, Pressure() ); + auto s_p = Picasso::get( particle, Stress() ); + auto v_p = Picasso::get( particle, Volume() ); + auto m_p = Picasso::get( particle, Mass() ); // Get the gather dependencies. - auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Example::Velocity() ); + auto u_i = + gather_deps.get( Picasso::FieldLocation::Node(), Velocity() ); // update strain rate auto spline = Picasso::createSpline( @@ -114,13 +111,13 @@ struct ComputeGridVelocityChangeStress const LocalDependencies&, ParticleViewType& particle ) const { // Get particle data. - auto v_p = Picasso::get( particle, Example::Volume() ); - auto x_p = Picasso::get( particle, Example::Position() ); - auto s_p = Picasso::get( particle, Example::Stress() ); + auto v_p = Picasso::get( particle, Volume() ); + auto x_p = Picasso::get( particle, Position() ); + auto s_p = Picasso::get( particle, Stress() ); // Get the scatter dependencies. - auto delta_u_s_i = scatter_deps.get( Picasso::FieldLocation::Node(), - Example::DeltaUStress() ); + auto delta_u_s_i = + scatter_deps.get( Picasso::FieldLocation::Node(), DeltaUStress() ); // Node interpolant. auto spline = Picasso::createSpline( @@ -151,12 +148,12 @@ struct ComputeGridVelocityChangeGravity const LocalDependencies&, ParticleViewType& particle ) const { // Get particle data. - auto m_p = Picasso::get( particle, Example::Mass() ); - auto x_p = Picasso::get( particle, Example::Position() ); + auto m_p = Picasso::get( particle, Mass() ); + auto x_p = Picasso::get( particle, Position() ); // Get the scatter dependencies. - auto delta_u_g_i = scatter_deps.get( Picasso::FieldLocation::Node(), - Example::DeltaUGravity() ); + auto delta_u_g_i = + scatter_deps.get( Picasso::FieldLocation::Node(), DeltaUGravity() ); // Node interpolant. auto spline = Picasso::createSpline( @@ -189,15 +186,14 @@ struct UpdateGridVelocity const int i, const int j, const int k ) const { // Get the local dependencies. - auto m_i = - gather_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); - auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Example::Velocity() ); - auto delta_u_s_i = gather_deps.get( Picasso::FieldLocation::Node(), - Example::DeltaUStress() ); + auto m_i = gather_deps.get( Picasso::FieldLocation::Node(), Mass() ); + auto u_i = + gather_deps.get( Picasso::FieldLocation::Node(), Velocity() ); + auto delta_u_s_i = + gather_deps.get( Picasso::FieldLocation::Node(), DeltaUStress() ); - auto delta_u_g_i = gather_deps.get( Picasso::FieldLocation::Node(), - Example::DeltaUGravity() ); + auto delta_u_g_i = + gather_deps.get( Picasso::FieldLocation::Node(), DeltaUGravity() ); // Compute velocity. Picasso::Vec3 zeros = { 0.0, 0.0, 0.0 }; @@ -245,16 +241,14 @@ template ; + using mass_type = Picasso::FieldLayout; using velocity_type = - Picasso::FieldLayout; - using old_u_type = - Picasso::FieldLayout; - using delta_u_s_type = Picasso::FieldLayout; - using delta_u_g_type = Picasso::FieldLayout; + Picasso::FieldLayout; + using old_u_type = Picasso::FieldLayout; + using delta_u_s_type = + Picasso::FieldLayout; + using delta_u_g_type = + Picasso::FieldLayout; using grid_type = Picasso::FieldLayout>; @@ -359,7 +353,7 @@ class ExplicitMomentumIntegrator const int spline_order = 1; // P2G - Particle2Grid p2g_func{ _dt }; _p2g_momentum.apply( "Picasso::p2g_U", @@ -399,7 +393,7 @@ class ExplicitMomentumIntegrator update_u_func ); // Apply boundary condition. - ApplyBoundaryCondition apply_bc{ bc }; + ApplyBoundaryCondition apply_bc{ bc }; _apply_bc_momentum.apply( "Picasso::BC_U", Picasso::FieldLocation::Node(), exec_space, fm, apply_bc ); @@ -412,7 +406,7 @@ class ExplicitMomentumIntegrator particles, g2p_func ); // Do not force particles redistribution - particles.redistribute( local_grid, Example::Position() ); + particles.redistribute( local_grid, Position() ); _total_time += _dt; _total_steps += 1; diff --git a/examples/sources/Picasso_Output.hpp b/examples/sources/Picasso_Output.hpp deleted file mode 100644 index 4129901c..00000000 --- a/examples/sources/Picasso_Output.hpp +++ /dev/null @@ -1,65 +0,0 @@ -#ifndef PICASSO_OUTPUT_HPP -#define PICASSO_OUTPUT_HPP - -#include - -#include -#include - -namespace Picasso -{ -namespace Output -{ - -namespace Impl -{ - -Cabana::Experimental::HDF5ParticleOutput::HDF5Config setupHDF5Config() -{ - Cabana::Experimental::HDF5ParticleOutput::HDF5Config h5_config; - h5_config.collective = true; - return h5_config; -} - -} // namespace Impl -// Momentum-only output. -template -void outputParticles( MPI_Comm comm, ExecutionSpace, ParticleVelocity, - const int n, const int write_frequency, const double time, - const ParticleList particles ) -{ - if ( write_frequency > 0 && n % write_frequency == 0 ) - { - auto h5_config = Impl::setupHDF5Config(); - - auto x = particles.slice( Picasso::Example::Position() ); - auto p_p = particles.slice( Picasso::Example::Pressure() ); - auto vel_p = particles.slice( ParticleVelocity() ); - auto m_p = particles.slice( Picasso::Example::Mass() ); - auto v_p = particles.slice( Picasso::Example::Volume() ); - - std::string prefix = "particles"; - Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( - h5_config, prefix, comm, n / write_frequency, time, x.size(), x, - p_p, vel_p, m_p, v_p ); - } -} - -void outputTimestep( MPI_Comm comm, const int n, const int write_frequency, - const double total_time, const double dt ) -{ - if ( write_frequency > 0 && n % write_frequency == 0 ) - { - int rank; - MPI_Comm_rank( comm, &rank ); - if ( rank == 0 ) - std::cout << "Current step: " << n - << "; Current time: " << total_time - << "; Current dt: " << dt << std::endl; - } -} - -} // namespace Output -} // namespace Picasso - -#endif diff --git a/src/Picasso_FieldTypes.hpp b/src/Picasso_FieldTypes.hpp index 00011c27..6a01bb2c 100644 --- a/src/Picasso_FieldTypes.hpp +++ b/src/Picasso_FieldTypes.hpp @@ -970,9 +970,6 @@ struct CommRank : Scalar } // end namespace Field -namespace Example -{ - struct Mass : Field::Scalar { static std::string label() { return "Mass"; } @@ -1020,8 +1017,6 @@ struct DetDefGrad : Field::Scalar using Position = Picasso::Field::LogicalPosition<3>; -} // namespace Example - namespace PolyPIC { namespace Field diff --git a/src/Picasso_ParticleInterpolation.hpp b/src/Picasso_ParticleInterpolation.hpp index 73384865..d3d99a2b 100644 --- a/src/Picasso_ParticleInterpolation.hpp +++ b/src/Picasso_ParticleInterpolation.hpp @@ -306,12 +306,11 @@ struct Particle2Grid { // Get particle data. auto u_p = Picasso::get( particle, PolyPIC::Field::Velocity() ); - auto x_p = Picasso::get( particle, Example::Position() ); + auto x_p = Picasso::get( particle, Position() ); // Get the gather dependencies. - auto m_i = - gather_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); - auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Example::Velocity() ); + auto m_i = gather_deps.get( Picasso::FieldLocation::Node(), Mass() ); + auto u_i = + gather_deps.get( Picasso::FieldLocation::Node(), Velocity() ); // Get the local dependencies for getting physcial location of node auto x_i = local_deps.get( Picasso::FieldLocation::Node(), Picasso::Field::PhysicalPosition<3>() ); @@ -491,13 +487,12 @@ struct Grid2ParticleVelocity { // Get particle data. auto f_p = Picasso::get( particle, APIC::Field::Velocity() ); - auto x_p = Picasso::get( particle, Example::Position() ); + auto x_p = Picasso::get( particle, Position() ); // Get the gather dependencies. - auto m_i = - gather_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); - auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Example::Velocity() ); + auto m_i = gather_deps.get( Picasso::FieldLocation::Node(), Mass() ); + auto u_i = + gather_deps.get( Picasso::FieldLocation::Node(), Velocity() ); // Get the local dependencies for getting physcial location of node auto x_i = local_deps.get( Picasso::FieldLocation::Node(), Picasso::Field::PhysicalPosition<3>() ); @@ -545,16 +540,15 @@ struct Grid2ParticleVelocity ParticleViewType& particle ) const { // Get particle data. - auto u_p = Picasso::get( particle, Example::Velocity() ); - auto x_p = Picasso::get( particle, Example::Position() ); + auto u_p = Picasso::get( particle, Velocity() ); + auto x_p = Picasso::get( particle, Position() ); // Get the gather dependencies. - auto m_i = - gather_deps.get( Picasso::FieldLocation::Node(), Example::Mass() ); - auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Example::Velocity() ); + auto m_i = gather_deps.get( Picasso::FieldLocation::Node(), Mass() ); + auto u_i = + gather_deps.get( Picasso::FieldLocation::Node(), Velocity() ); auto old_u_i = - gather_deps.get( Picasso::FieldLocation::Node(), Example::OldU() ); + gather_deps.get( Picasso::FieldLocation::Node(), OldU() ); // Get the local dependencies for getting physcial location of node auto x_i = local_deps.get( Picasso::FieldLocation::Node(), From 35e9aa4df261761065c5d89eb4eb4410d23f8ab4 Mon Sep 17 00:00:00 2001 From: "Chong, Kwitae" Date: Tue, 1 Oct 2024 22:06:41 -0400 Subject: [PATCH 08/22] fix boundary condition --- examples/dam_break.cpp | 22 +++++---- .../sources/Picasso_BoundaryCondition.hpp | 46 +++++++++++++------ .../Picasso_ExplicitMomentumUpdate.hpp | 7 --- src/Picasso_ParticleInterpolation.hpp | 12 +---- 4 files changed, 47 insertions(+), 40 deletions(-) diff --git a/examples/dam_break.cpp b/examples/dam_break.cpp index 5830438c..4ac8e30c 100644 --- a/examples/dam_break.cpp +++ b/examples/dam_break.cpp @@ -53,11 +53,11 @@ void DamBreak() *local_grid ); // Boundary index space - auto bc_index = local_grid->boundaryIndexSpace( - Cabana::Grid::Own(), Cabana::Grid::Node(), -1, 1, 0 ); - using bc_index_type = decltype( bc_index ); + //auto bc_index = local_grid->boundaryIndexSpace( + // Cabana::Grid::Own(), Cabana::Grid::Node(), -1, 1, 0 ); + using local_grid_type = decltype( *local_grid ); - BoundaryCondition bc{ bc_index }; + BoundaryCondition bc{ *local_grid }; // Properties auto gamma = inputs["gamma"]; @@ -74,24 +74,26 @@ void DamBreak() // steps auto final_time = inputs["final_time"]; - auto write_frequency = inputs["write_frequency"]; + int write_frequency = inputs["write_frequency"]; + int step = 0; while ( time_integrator.totalTime() < final_time ) { // Write particle fields. Cabana::Experimental::HDF5ParticleOutput::HDF5Config h5_config; + if ( step % write_frequency == 0 ) Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( h5_config, "particles", MPI_COMM_WORLD, - time_integrator.totalSteps(), time_integrator.totalTime(), + step / write_frequency, time_integrator.totalTime(), particles.size(), particles.slice( Picasso::Position() ), particles.slice( Picasso::Pressure() ), particles.slice( ParticleVelocity() ), particles.slice( Picasso::Mass() ), particles.slice( Picasso::Volume() ) ); - printf( "aaa\n" ); - // Step. time_integrator.step( exec_space(), *fm, particles, *local_grid, bc ); + + step++; } } @@ -108,7 +110,9 @@ int main( int argc, char* argv[] ) $/: ./DamBreak inputs/dam_break.json\n" ); std::string filename = argv[1]; - DamBreak(); + //DamBreak(); + //DamBreak(); + DamBreak(); Kokkos::finalize(); MPI_Finalize(); diff --git a/examples/sources/Picasso_BoundaryCondition.hpp b/examples/sources/Picasso_BoundaryCondition.hpp index a745b240..926463a3 100644 --- a/examples/sources/Picasso_BoundaryCondition.hpp +++ b/examples/sources/Picasso_BoundaryCondition.hpp @@ -7,27 +7,45 @@ namespace Picasso { -template +template struct BoundaryCondition { - BoundaryIndexType bc_index; + LocalGridType local_grid; + // Free slip boundary condition template KOKKOS_INLINE_FUNCTION void apply( ViewType view, const int i, const int j, const int k ) const { - // x faces - if ( bc_index.min( Cabana::Grid::Dim::I ) <= i && - i < bc_index.max( Cabana::Grid::Dim::I ) ) - view( i, j, k, 0 ) = 0.0; - // y faces - if ( bc_index.min( Cabana::Grid::Dim::J ) <= j && - j < bc_index.max( Cabana::Grid::Dim::J ) ) - view( i, j, k, 1 ) = 0.0; - // z faces - if ( bc_index.min( Cabana::Grid::Dim::K ) <= k && - k < bc_index.max( Cabana::Grid::Dim::K ) ) - view( i, j, k, 2 ) = 0.0; + + auto index_space = local_grid.indexSpace( Cabana::Grid::Own(), Cabana::Grid::Node(), Cabana::Grid::Local() ); + + auto global_grid = local_grid.globalGrid(); + + // -x face + if ( i == index_space.min( Cabana::Grid::Dim::I ) && + global_grid.onLowBoundary( Cabana::Grid::Dim::I ) ) + view( i, j, k, 0 ) = 0.0; + // +x face + if ( i == index_space.min( Cabana::Grid::Dim::I ) && + global_grid.onHighBoundary( Cabana::Grid::Dim::I ) ) + view( i, j, k, 0 ) = 0.0; + // -y face + if ( j == index_space.min( Cabana::Grid::Dim::J ) && + global_grid.onLowBoundary( Cabana::Grid::Dim::J ) ) + view( i, j, k, 1 ) = 0.0; + // +y face + if ( j == index_space.min( Cabana::Grid::Dim::J ) && + global_grid.onHighBoundary( Cabana::Grid::Dim::J ) ) + view( i, j, k, 1 ) = 0.0; + // -z face + if ( k == index_space.min( Cabana::Grid::Dim::K ) && + global_grid.onLowBoundary( Cabana::Grid::Dim::K ) ) + view( i, j, k, 2 ) = 0.0; + // +z face + if ( k == index_space.min( Cabana::Grid::Dim::K ) && + global_grid.onHighBoundary( Cabana::Grid::Dim::K ) ) + view( i, j, k, 2 ) = 0.0; } }; diff --git a/examples/sources/Picasso_ExplicitMomentumUpdate.hpp b/examples/sources/Picasso_ExplicitMomentumUpdate.hpp index e895e371..f34b8bd0 100644 --- a/examples/sources/Picasso_ExplicitMomentumUpdate.hpp +++ b/examples/sources/Picasso_ExplicitMomentumUpdate.hpp @@ -90,9 +90,6 @@ struct ComputeParticlePressure Picasso::LinearAlgebra::identity( I ); s_p = -p_p * I; - - printf( "particle %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e %8.5e\n", J_p, - p_p, v_p, m_p, x_p( 0 ), x_p( 1 ), x_p( 2 ) ); } }; @@ -202,10 +199,6 @@ struct UpdateGridVelocity ? dt * ( delta_u_s_i( i, j, k ) + delta_u_g_i( i, j, k ) ) / m_i( i, j, k ) : zeros; - - printf( "%d %d %d %8.5e %8.5e %8.5e %8.5e\n", i, j, k, m_i( i, j, k ), - delta_u_s_i( i, j, k, 2 ), delta_u_g_i( i, j, k, 2 ), - u_i( i, j, k, 2 ) ); } }; diff --git a/src/Picasso_ParticleInterpolation.hpp b/src/Picasso_ParticleInterpolation.hpp index d3d99a2b..3f1e0a4e 100644 --- a/src/Picasso_ParticleInterpolation.hpp +++ b/src/Picasso_ParticleInterpolation.hpp @@ -455,11 +455,7 @@ struct Grid2ParticleVelocity // Update particle position. auto x_i_updated = [=]( const int i, const int j, const int k, const int d ) - { - printf( "g2p %d %d %d %8.5e %8.5e\n", i, j, k, x_i( i, j, k, 2 ), - u_i( i, j, k, 2 ) ); - return x_i( i, j, k, d ) + dt * u_i( i, j, k, d ); - }; + { return x_i( i, j, k, d ) + dt * u_i( i, j, k, d ); }; Picasso::G2P::value( spline, x_i_updated, x_p ); } }; @@ -510,11 +506,7 @@ struct Grid2ParticleVelocity // Update particle position. auto x_i_updated = [=]( const int i, const int j, const int k, const int d ) - { - printf( "g2p %d %d %d %8.5e %8.5e\n", i, j, k, x_i( i, j, k, 2 ), - u_i( i, j, k, 2 ) ); - return x_i( i, j, k, d ) + dt * u_i( i, j, k, d ); - }; + { return x_i( i, j, k, d ) + dt * u_i( i, j, k, d ); }; Picasso::G2P::value( spline, x_i_updated, x_p ); } }; From 898277d815ad0d17d165b48b768f7e3234d083ca Mon Sep 17 00:00:00 2001 From: "Chong, Kwitae" Date: Wed, 2 Oct 2024 11:24:48 -0400 Subject: [PATCH 09/22] fix on gpu --- examples/dam_break.cpp | 25 +++++++++++++---- .../sources/Picasso_BoundaryCondition.hpp | 27 ++++++------------- 2 files changed, 28 insertions(+), 24 deletions(-) diff --git a/examples/dam_break.cpp b/examples/dam_break.cpp index 4ac8e30c..72ad732d 100644 --- a/examples/dam_break.cpp +++ b/examples/dam_break.cpp @@ -53,11 +53,26 @@ void DamBreak() *local_grid ); // Boundary index space - //auto bc_index = local_grid->boundaryIndexSpace( - // Cabana::Grid::Own(), Cabana::Grid::Node(), -1, 1, 0 ); - using local_grid_type = decltype( *local_grid ); - - BoundaryCondition bc{ *local_grid }; + auto index_space = local_grid->indexSpace( Cabana::Grid::Own(), Cabana::Grid::Node(), Cabana::Grid::Local() ); + auto global_grid = local_grid->globalGrid(); + + Kokkos::Array 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 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 auto gamma = inputs["gamma"]; diff --git a/examples/sources/Picasso_BoundaryCondition.hpp b/examples/sources/Picasso_BoundaryCondition.hpp index 926463a3..769665a1 100644 --- a/examples/sources/Picasso_BoundaryCondition.hpp +++ b/examples/sources/Picasso_BoundaryCondition.hpp @@ -7,44 +7,33 @@ namespace Picasso { -template struct BoundaryCondition { - LocalGridType local_grid; + Kokkos::Array bc_index_space; + Kokkos::Array on_boundary; // Free slip boundary condition template KOKKOS_INLINE_FUNCTION void apply( ViewType view, const int i, const int j, const int k ) const { - - auto index_space = local_grid.indexSpace( Cabana::Grid::Own(), Cabana::Grid::Node(), Cabana::Grid::Local() ); - - auto global_grid = local_grid.globalGrid(); - // -x face - if ( i == index_space.min( Cabana::Grid::Dim::I ) && - global_grid.onLowBoundary( Cabana::Grid::Dim::I ) ) + if ( i == bc_index_space[0] && on_boundary[0] ) view( i, j, k, 0 ) = 0.0; // +x face - if ( i == index_space.min( Cabana::Grid::Dim::I ) && - global_grid.onHighBoundary( Cabana::Grid::Dim::I ) ) + if ( i == bc_index_space[3] && on_boundary[3] ) view( i, j, k, 0 ) = 0.0; // -y face - if ( j == index_space.min( Cabana::Grid::Dim::J ) && - global_grid.onLowBoundary( Cabana::Grid::Dim::J ) ) + if ( j == bc_index_space[1] && on_boundary[1] ) view( i, j, k, 1 ) = 0.0; // +y face - if ( j == index_space.min( Cabana::Grid::Dim::J ) && - global_grid.onHighBoundary( Cabana::Grid::Dim::J ) ) + if ( j == bc_index_space[4] && on_boundary[4] ) view( i, j, k, 1 ) = 0.0; // -z face - if ( k == index_space.min( Cabana::Grid::Dim::K ) && - global_grid.onLowBoundary( Cabana::Grid::Dim::K ) ) + if ( k == bc_index_space[2] && on_boundary[2] ) view( i, j, k, 2 ) = 0.0; // +z face - if ( k == index_space.min( Cabana::Grid::Dim::K ) && - global_grid.onHighBoundary( Cabana::Grid::Dim::K ) ) + if ( k == bc_index_space[5] && on_boundary[5] ) view( i, j, k, 2 ) = 0.0; } }; From 709404b85759b9ff68a7a53e42230a692214422b Mon Sep 17 00:00:00 2001 From: kwitaechong Date: Thu, 3 Oct 2024 10:35:25 -0400 Subject: [PATCH 10/22] format --- examples/dam_break.cpp | 59 ++++++++++--------- .../sources/Picasso_BoundaryCondition.hpp | 14 ++--- 2 files changed, 37 insertions(+), 36 deletions(-) diff --git a/examples/dam_break.cpp b/examples/dam_break.cpp index 72ad732d..f46f5e38 100644 --- a/examples/dam_break.cpp +++ b/examples/dam_break.cpp @@ -53,26 +53,27 @@ void DamBreak() *local_grid ); // Boundary index space - auto index_space = local_grid->indexSpace( Cabana::Grid::Own(), Cabana::Grid::Node(), Cabana::Grid::Local() ); + auto index_space = local_grid->indexSpace( + Cabana::Grid::Own(), Cabana::Grid::Node(), Cabana::Grid::Local() ); auto global_grid = local_grid->globalGrid(); - - Kokkos::Array 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 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 }; + + Kokkos::Array 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 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 auto gamma = inputs["gamma"]; @@ -96,14 +97,14 @@ index_space.max( Cabana::Grid::Dim::K ) - 1}; // Write particle fields. Cabana::Experimental::HDF5ParticleOutput::HDF5Config h5_config; if ( step % write_frequency == 0 ) - Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( - h5_config, "particles", MPI_COMM_WORLD, - step / write_frequency, time_integrator.totalTime(), - particles.size(), particles.slice( Picasso::Position() ), - particles.slice( Picasso::Pressure() ), - particles.slice( ParticleVelocity() ), - particles.slice( Picasso::Mass() ), - particles.slice( Picasso::Volume() ) ); + Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( + h5_config, "particles", MPI_COMM_WORLD, step / write_frequency, + time_integrator.totalTime(), particles.size(), + particles.slice( Picasso::Position() ), + particles.slice( Picasso::Pressure() ), + particles.slice( ParticleVelocity() ), + particles.slice( Picasso::Mass() ), + particles.slice( Picasso::Volume() ) ); // Step. time_integrator.step( exec_space(), *fm, particles, *local_grid, bc ); @@ -125,8 +126,8 @@ int main( int argc, char* argv[] ) $/: ./DamBreak inputs/dam_break.json\n" ); std::string filename = argv[1]; - //DamBreak(); - //DamBreak(); + // DamBreak(); + // DamBreak(); DamBreak(); Kokkos::finalize(); diff --git a/examples/sources/Picasso_BoundaryCondition.hpp b/examples/sources/Picasso_BoundaryCondition.hpp index 769665a1..34305b96 100644 --- a/examples/sources/Picasso_BoundaryCondition.hpp +++ b/examples/sources/Picasso_BoundaryCondition.hpp @@ -10,7 +10,7 @@ namespace Picasso struct BoundaryCondition { Kokkos::Array bc_index_space; - Kokkos::Array on_boundary; + Kokkos::Array on_boundary; // Free slip boundary condition template @@ -19,22 +19,22 @@ struct BoundaryCondition { // -x face if ( i == bc_index_space[0] && on_boundary[0] ) - view( i, j, k, 0 ) = 0.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; + 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; + 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; + 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; + 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; + view( i, j, k, 2 ) = 0.0; } }; From 104d061e933e8c084d10a277b65cc457862684a8 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Sat, 5 Oct 2024 13:54:54 -0400 Subject: [PATCH 11/22] Rearrange example, src --- examples/CMakeLists.txt | 4 +- examples/dam_break.cpp | 369 +++++++++++++-- examples/inputs/dam_break.json | 1 + examples/sources/CMakeLists.txt | 0 examples/sources/Example_Sources.hpp | 3 - examples/sources/Particle_Init.hpp | 58 --- .../sources/Picasso_BoundaryCondition.hpp | 59 --- .../Picasso_ExplicitMomentumUpdate.hpp | 447 ------------------ src/CMakeLists.txt | 1 + src/Picasso.hpp | 1 + src/Picasso_FieldTypes.hpp | 18 +- src/Picasso_GridUpdate.hpp | 57 +++ src/Picasso_ParticleInterpolation.hpp | 61 +-- 13 files changed, 426 insertions(+), 653 deletions(-) delete mode 100644 examples/sources/CMakeLists.txt delete mode 100644 examples/sources/Example_Sources.hpp delete mode 100644 examples/sources/Particle_Init.hpp delete mode 100644 examples/sources/Picasso_BoundaryCondition.hpp delete mode 100644 examples/sources/Picasso_ExplicitMomentumUpdate.hpp create mode 100644 src/Picasso_GridUpdate.hpp diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index cb970489..c06f315d 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,12 +1,10 @@ -add_subdirectory(sources) - 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 ) +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 index f46f5e38..36fa163a 100644 --- a/examples/dam_break.cpp +++ b/examples/dam_break.cpp @@ -1,9 +1,5 @@ #include -#include "sources/Particle_Init.hpp" -#include "sources/Picasso_BoundaryCondition.hpp" -#include "sources/Picasso_ExplicitMomentumUpdate.hpp" - #include #include @@ -13,28 +9,242 @@ using namespace Picasso; +// Dam break particle creation and initialization. +template +struct ParticleInitFunc +{ + Kokkos::Array 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() ) = 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::Position(), d ) = x[d]; + return true; + } + + return false; + } +}; + +// Dam break system boundaries. +struct BoundaryCondition +{ + Kokkos::Array bc_index_space; + Kokkos::Array on_boundary; + + // Free slip boundary condition + template + 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 +{ + static std::string label() { return "velocity_change_from_gravity"; } +}; +struct DeltaUStress : Field::Vector +{ + static std::string label() { return "velocity_change_from_stress"; } +}; + +//---------------------------------------------------------------------------// +// Update particle stress. +//---------------------------------------------------------------------------// +template +struct ComputeParticlePressure +{ + double dt; + double bulk_modulus; + double gamma; + + template + 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::Position() ); + 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() ); + 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() ); + + // update strain rate + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::InterpolationOrder(), local_mesh, x_p, + Picasso::SplineValue(), Picasso::SplineGradient() ); + Picasso::Mat3 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 I; + Picasso::LinearAlgebra::identity( I ); + + s_p = -p_p * I; + } +}; + +//---------------------------------------------------------------------------// +// Grid momentum change due to stress +//---------------------------------------------------------------------------// +template +struct ComputeGridVelocityChange +{ + Picasso::Vec3 gravity; + + ComputeGridVelocityChange( std::array g ) + { + for ( int d = 0; d < 3; ++d ) + gravity( d ) = g[d]; + } + + template + 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::Position() ); + auto s_p = Picasso::get( particle, Picasso::Field::Stress() ); + + // 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(), 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 + 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() ); + 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 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 -void DamBreak() +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( "dam_break.json" ); + auto inputs = Picasso::parse( filename ); // Global bounding box. auto global_box = copy( inputs["global_box"] ); int minimum_halo_size = 0; // Make mesh. - auto mesh = createUniformMesh( memory_space(), inputs, global_box, - minimum_halo_size, MPI_COMM_WORLD ); + using mesh_type = Picasso::UniformMesh; + std::shared_ptr mesh = Picasso::createUniformMesh( + memory_space(), inputs, global_box, minimum_halo_size, MPI_COMM_WORLD ); // Make a particle list. - Cabana::ParticleTraits + Cabana::ParticleTraits fields; auto particles = Cabana::Grid::createParticleList( "test_particles", fields ); @@ -42,13 +252,12 @@ void DamBreak() // Initialize particles auto particle_box = copy( inputs["particle_box"] ); double density = inputs["density"]; - - auto momentum_init_functor = createParticleInitFunc( - particles, ParticleVelocity(), particle_box, density ); + ParticleInitFunc + momentum_init_functor{ particle_box, density }; double ppc = inputs["ppc"]; auto local_grid = mesh->localGrid(); - Cabana::Grid::createParticles( Cabana::InitUniform(), exec_space(), + Cabana::Grid::createParticles( Cabana::InitUniform(), exec_space{}, momentum_init_functor, particles, ppc, *local_grid ); @@ -76,40 +285,121 @@ void DamBreak() BoundaryCondition bc{ bc_index_space, on_boundary }; // Properties - auto gamma = inputs["gamma"]; - auto bulk_modulus = inputs["bulk_modulus"]; - auto gravity = copy( inputs["gravity"] ); - Picasso::Properties props( gamma, bulk_modulus, gravity ); - - // Time integragor - auto dt = inputs["dt"]; - auto time_integrator = Picasso::createExplicitMomentumIntegrator( - mesh, InterpolationType(), ParticleVelocity(), props, dt ); + double gamma = inputs["gamma"]; + double bulk_modulus = inputs["bulk_modulus"]; + auto gravity = 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 ); - time_integrator.setup( *fm ); - // steps + // Field locations and variables. + using node_type = Picasso::FieldLocation::Node; + using grid_type = + Picasso::FieldLayout>; + using particle_type = Picasso::FieldLocation::Particle; + using mass_type = Picasso::FieldLayout; + using velocity_type = + Picasso::FieldLayout; + using old_u_type = Picasso::FieldLayout; + using delta_u_s_type = Picasso::FieldLayout; + using delta_u_g_type = Picasso::FieldLayout; + + // Define particle/grid parallel data dependencies. + Picasso::GridOperator> + p2g( mesh ); + p2g.setup( *fm ); + + Picasso::GridOperator> + compute_velocity( mesh ); + compute_velocity.setup( *fm ); + + Picasso::GridOperator> + compute_stress( mesh ); + compute_stress.setup( *fm ); + + Picasso::GridOperator< + mesh_type, Picasso::ScatterDependencies> + compute_velocity_change( mesh ); + compute_velocity_change.setup( *fm ); + + Picasso::GridOperator< + mesh_type, Picasso::GatherDependencies> + update_velocity( mesh ); + update_velocity.setup( *fm ); + + Picasso::GridOperator< + mesh_type, + Picasso::GatherDependencies, + Picasso::LocalDependencies> + g2p( mesh ); + g2p.setup( *fm ); + + // Timestep loop. + double time = 0.0; auto final_time = inputs["final_time"]; int write_frequency = inputs["write_frequency"]; - int step = 0; - while ( time_integrator.totalTime() < final_time ) + int steps = 0; + while ( time < final_time ) { + // Particle interpolation (Picasso built-in). + Picasso::Particle2Grid + p2g_func{ 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 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 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). + Picasso::Grid2ParticleVelocity + g2p_func{ beta, dt }; + g2p.apply( "g2p_U", particle_type{}, exec_space{}, *fm, particles, + g2p_func ); + + // Redistribution of particles across MPI subdomains. + particles.redistribute( *local_grid, Picasso::Field::Position() ); + // Write particle fields. Cabana::Experimental::HDF5ParticleOutput::HDF5Config h5_config; - if ( step % write_frequency == 0 ) + if ( steps % write_frequency == 0 ) Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( - h5_config, "particles", MPI_COMM_WORLD, step / write_frequency, - time_integrator.totalTime(), particles.size(), - particles.slice( Picasso::Position() ), - particles.slice( Picasso::Pressure() ), + h5_config, "particles", MPI_COMM_WORLD, steps / write_frequency, + time, particles.size(), + particles.slice( Picasso::Field::Position() ), + particles.slice( Picasso::Field::Pressure() ), particles.slice( ParticleVelocity() ), - particles.slice( Picasso::Mass() ), - particles.slice( Picasso::Volume() ) ); - - // Step. - time_integrator.step( exec_space(), *fm, particles, *local_grid, bc ); + particles.slice( Picasso::Field::Mass() ), + particles.slice( Picasso::Field::Volume() ) ); - step++; + time += dt; + steps++; } } @@ -126,9 +416,10 @@ int main( int argc, char* argv[] ) $/: ./DamBreak inputs/dam_break.json\n" ); std::string filename = argv[1]; + // Problem can run with any interpolation scheme. // DamBreak(); // DamBreak(); - DamBreak(); + DamBreak( filename ); Kokkos::finalize(); MPI_Finalize(); diff --git a/examples/inputs/dam_break.json b/examples/inputs/dam_break.json index b83847ce..31de2a87 100644 --- a/examples/inputs/dam_break.json +++ b/examples/inputs/dam_break.json @@ -15,6 +15,7 @@ "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/examples/sources/CMakeLists.txt b/examples/sources/CMakeLists.txt deleted file mode 100644 index e69de29b..00000000 diff --git a/examples/sources/Example_Sources.hpp b/examples/sources/Example_Sources.hpp deleted file mode 100644 index bcf38d00..00000000 --- a/examples/sources/Example_Sources.hpp +++ /dev/null @@ -1,3 +0,0 @@ -#include "Particle_Init.hpp" -#include "Picasso_BoundaryCondition.hpp" -#include "Picasso_ExplicitMomentumUpdate.hpp" diff --git a/examples/sources/Particle_Init.hpp b/examples/sources/Particle_Init.hpp deleted file mode 100644 index c6bd1178..00000000 --- a/examples/sources/Particle_Init.hpp +++ /dev/null @@ -1,58 +0,0 @@ -#include - -#include -#include - -#include - -#include - -namespace Picasso -{ - -template -struct ParticleInitFunc -{ - Kokkos::Array block; - double density; - - ParticleInitFunc( const Kokkos::Array& _block, - const double _density ) - : block( _block ) - , density( _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::Stress() ) = 0.0; - Picasso::get( p, VelocityType() ) = 0.0; - Picasso::get( p, Picasso::DetDefGrad() ) = 1.0; - Picasso::get( p, Picasso::Mass() ) = pv * density; - Picasso::get( p, Picasso::Volume() ) = pv; - Picasso::get( p, Picasso::Pressure() ) = 0.0; - - for ( int d = 0; d < 3; ++d ) - Picasso::get( p, Picasso::Position(), d ) = x[d]; - return true; - } - - return false; - } -}; - -template -auto createParticleInitFunc( const ParticleType&, const VelocityType&, - const Kokkos::Array& block, - const double density ) -{ - return ParticleInitFunc{ block, density }; -} - -} // namespace Picasso diff --git a/examples/sources/Picasso_BoundaryCondition.hpp b/examples/sources/Picasso_BoundaryCondition.hpp deleted file mode 100644 index 34305b96..00000000 --- a/examples/sources/Picasso_BoundaryCondition.hpp +++ /dev/null @@ -1,59 +0,0 @@ -#ifndef PICASSO_BOUNDARYCONDITION_HPP -#define PICASSO_BOUNDARYCONDITION_HPP - -#include -#include - -namespace Picasso -{ - -struct BoundaryCondition -{ - Kokkos::Array bc_index_space; - Kokkos::Array on_boundary; - - // Free slip boundary condition - template - 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; - } -}; - -struct Properties -{ - double gamma; - double bulk_modulus; - Kokkos::Array gravity; - - KOKKOS_INLINE_FUNCTION - Properties( const double _gamma, const double _bulk_modulus, - const Kokkos::Array _gravity ) - : gamma( _gamma ) - , bulk_modulus( _bulk_modulus ) - , gravity( _gravity ) - { - } -}; - -} // end namespace Picasso - -#endif // end PICASSO_BOUNDARYCONDITION_HPP diff --git a/examples/sources/Picasso_ExplicitMomentumUpdate.hpp b/examples/sources/Picasso_ExplicitMomentumUpdate.hpp deleted file mode 100644 index f34b8bd0..00000000 --- a/examples/sources/Picasso_ExplicitMomentumUpdate.hpp +++ /dev/null @@ -1,447 +0,0 @@ -#include - -#include "Picasso_BoundaryCondition.hpp" - -#include -#include - -#include - -#include - -namespace Picasso -{ - -//---------------------------------------------------------------------------// -// Compute nodal velocity from mass-weighted momentum. -//---------------------------------------------------------------------------// -struct ComputeGridVelocity -{ - template - 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(), Mass() ); - auto u_i = local_deps.get( Picasso::FieldLocation::Node(), Velocity() ); - auto old_u_i = local_deps.get( Picasso::FieldLocation::Node(), OldU() ); - - // 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 ); - } - } -}; - -//---------------------------------------------------------------------------// -// Update particle stress. -//---------------------------------------------------------------------------// -template -struct ComputeParticlePressure -{ - double dt; - - // Stress property functions - StressProperties properties; - - template - 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, Position() ); - auto& J_p = Picasso::get( particle, DetDefGrad() ); - auto& p_p = Picasso::get( particle, Pressure() ); - auto s_p = Picasso::get( particle, Stress() ); - auto v_p = Picasso::get( particle, Volume() ); - auto m_p = Picasso::get( particle, Mass() ); - - // Get the gather dependencies. - auto u_i = - gather_deps.get( Picasso::FieldLocation::Node(), Velocity() ); - - // update strain rate - auto spline = Picasso::createSpline( - Picasso::FieldLocation::Node(), - Picasso::InterpolationOrder(), local_mesh, x_p, - Picasso::SplineValue(), Picasso::SplineGradient() ); - Picasso::Mat3 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 = properties.bulk_modulus * - ( Kokkos::pow( J_p, -properties.gamma ) - 1.0 ); - - Picasso::Mat3 I; - Picasso::LinearAlgebra::identity( I ); - - s_p = -p_p * I; - } -}; - -//---------------------------------------------------------------------------// -// Grid momentum change due to stress -//---------------------------------------------------------------------------// -template -struct ComputeGridVelocityChangeStress -{ - template - KOKKOS_INLINE_FUNCTION void - operator()( const LocalMeshType& local_mesh, const GatherDependencies&, - const ScatterDependencies& scatter_deps, - const LocalDependencies&, ParticleViewType& particle ) const - { - // Get particle data. - auto v_p = Picasso::get( particle, Volume() ); - auto x_p = Picasso::get( particle, Position() ); - auto s_p = Picasso::get( particle, Stress() ); - - // Get the scatter dependencies. - auto delta_u_s_i = - scatter_deps.get( Picasso::FieldLocation::Node(), DeltaUStress() ); - - // Node interpolant. - auto spline = Picasso::createSpline( - Picasso::FieldLocation::Node(), - Picasso::InterpolationOrder(), local_mesh, x_p, - Picasso::SplineValue(), Picasso::SplineGradient() ); - - // Compute velocity update. - Picasso::P2G::divergence( spline, -v_p * s_p, delta_u_s_i ); - } -}; - -//---------------------------------------------------------------------------// -// Grid momentum change due to gravity -//---------------------------------------------------------------------------// -template -struct ComputeGridVelocityChangeGravity -{ - // Stress property functions - StressProperties properties; - - template - 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, Mass() ); - auto x_p = Picasso::get( particle, Position() ); - - // Get the scatter dependencies. - auto delta_u_g_i = - scatter_deps.get( Picasso::FieldLocation::Node(), DeltaUGravity() ); - - // Node interpolant. - auto spline = Picasso::createSpline( - Picasso::FieldLocation::Node(), - Picasso::InterpolationOrder(), local_mesh, x_p, - Picasso::SplineValue(), Picasso::SplineGradient() ); - - // Compute velocity update. - Picasso::Vec3 g_p = { properties.gravity[0], - properties.gravity[1], - properties.gravity[2] }; - - Picasso::P2G::value( spline, m_p * g_p, delta_u_g_i ); - } -}; - -//---------------------------------------------------------------------------// -// Update nodal momentum n+1 with stress, gravity -//---------------------------------------------------------------------------// -struct UpdateGridVelocity -{ - // Explicit time step size. - double dt; - - template - 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 local dependencies. - auto m_i = gather_deps.get( Picasso::FieldLocation::Node(), Mass() ); - auto u_i = - gather_deps.get( Picasso::FieldLocation::Node(), Velocity() ); - 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 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; - } -}; - -//---------------------------------------------------------------------------// -// Compute boundary condition on grid. -//---------------------------------------------------------------------------// -template -struct ApplyBoundaryCondition -{ - BCType bc; - - template - 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 f_i = - local_deps.get( Picasso::FieldLocation::Node(), FieldType() ); - - bc.apply( f_i, i, j, k ); - } -}; - -//---------------------------------------------------------------------------// -// Explicit momentum time stepper -//---------------------------------------------------------------------------// - -template -class ExplicitMomentumIntegrator -{ - public: - using mass_type = Picasso::FieldLayout; - using velocity_type = - Picasso::FieldLayout; - using old_u_type = Picasso::FieldLayout; - using delta_u_s_type = - Picasso::FieldLayout; - using delta_u_g_type = - Picasso::FieldLayout; - using grid_type = Picasso::FieldLayout>; - - // Momentum P2G operator types. - using p2g_scatter_deps = - Picasso::ScatterDependencies; - using p2g_u_op = Picasso::GridOperator; - - // Grid velocity operator types. - using compute_u_local_deps = - Picasso::LocalDependencies; - using compute_u_op = Picasso::GridOperator; - - // Update particle stress operator types - using compute_s_gather_deps = Picasso::GatherDependencies; - using compute_s_op = Picasso::GridOperator; - - // Grid boundary condition operator types. - using bc_u_local_deps = - Picasso::LocalDependencies; - using bc_u_op = Picasso::GridOperator; - - // Compute velocity change due to stress operator types. - using compute_du_s_scatter_deps = - Picasso::ScatterDependencies; - using compute_du_s_op = - Picasso::GridOperator; - - // Compute velocity change due to gravity operator types. - using compute_du_g_scatter_deps = - Picasso::ScatterDependencies; - using compute_du_g_op = - Picasso::GridOperator; - - // Grid velocity time update operator types. - using update_u_gather_deps = - Picasso::GatherDependencies; - using update_u_op = Picasso::GridOperator; - - // Particle update G2P types. - using g2p_u_gather_deps = - Picasso::GatherDependencies; - using g2p_u_local_deps = Picasso::LocalDependencies; - using g2p_u_op = - Picasso::GridOperator; - - using property_type = StressProperties; - using interpolation_type = InterpolationType; - using interpolation_variable = ParticleVelocity; - - public: - // Constructor. - ExplicitMomentumIntegrator( const std::shared_ptr& mesh, - StressProperties properties, - const double max_dt, const double cfl_number, - const double beta, const bool fixed_dt ) - : _max_dt( max_dt ) - , _cfl_number( cfl_number ) - , _beta( beta ) - , _fixed_dt( fixed_dt ) - , _p2g_momentum( mesh ) - , _compute_velocity( mesh ) - , _compute_stress( mesh ) - , _compute_delta_u_s( mesh ) - , _compute_delta_u_g( mesh ) - , _apply_bc_momentum( mesh ) - , _update_velocity( mesh ) - , _g2p_momentum( mesh ) - , _properties( properties ) - { - _total_steps = 0; - _total_time = 0.0; - _dt = _max_dt; - _cell_size = mesh->cellSize(); - } - - // Populate the field manager. - void setup( Picasso::FieldManager& fm ) - { - _p2g_momentum.setup( fm ); - _compute_velocity.setup( fm ); - _compute_stress.setup( fm ); - _compute_delta_u_s.setup( fm ); - _compute_delta_u_g.setup( fm ); - _apply_bc_momentum.setup( fm ); - _update_velocity.setup( fm ); - _g2p_momentum.setup( fm ); - } - - // Do a time step. - // Note that the boundary condition is passed here rather than saved as a - // member so that it can be initialized with the FieldManager. - template - void step( const ExecutionSpace& exec_space, - const Picasso::FieldManager& fm, ParticleList& particles, - const LocalGridType& local_grid, const BCType bc ) - { - // Spline interpolation order. - const int spline_order = 1; - - // P2G - Particle2Grid - p2g_func{ _dt }; - _p2g_momentum.apply( "Picasso::p2g_U", - Picasso::FieldLocation::Particle(), exec_space, fm, - particles, p2g_func ); - - // Compute grid velocity. - ComputeGridVelocity compute_u_func; - _compute_velocity.apply( "Picasso::grid_U", - Picasso::FieldLocation::Node(), exec_space, fm, - compute_u_func ); - - // Update particle stress. - ComputeParticlePressure compute_s_func{ - _dt, _properties }; - _compute_stress.apply( "Picasso::particle_S", - Picasso::FieldLocation::Particle(), exec_space, - fm, particles, compute_s_func ); - - // Compute grid velocity change due to stress - ComputeGridVelocityChangeStress compute_du_s_func; - _compute_delta_u_s.apply( - "Picasso::div_S", Picasso::FieldLocation::Particle(), exec_space, - fm, particles, compute_du_s_func ); - - // Compute grid velocity change due to gravity - ComputeGridVelocityChangeGravity - compute_du_g_func{ _properties }; - _compute_delta_u_g.apply( - "Picasso::rhoG", Picasso::FieldLocation::Particle(), exec_space, fm, - particles, compute_du_g_func ); - - // Compute next grid velocity. - UpdateGridVelocity update_u_func{ _dt }; - _update_velocity.apply( "Picasso::update_U", - Picasso::FieldLocation::Node(), exec_space, fm, - update_u_func ); - - // Apply boundary condition. - ApplyBoundaryCondition apply_bc{ bc }; - _apply_bc_momentum.apply( "Picasso::BC_U", - Picasso::FieldLocation::Node(), exec_space, - fm, apply_bc ); - - // G2P - Grid2ParticleVelocity g2p_func{ _beta, - _dt }; - _g2p_momentum.apply( "Picasso::g2p_U", - Picasso::FieldLocation::Particle(), exec_space, fm, - particles, g2p_func ); - - // Do not force particles redistribution - particles.redistribute( local_grid, Position() ); - - _total_time += _dt; - _total_steps += 1; - } - - double dt() { return _dt; } - double totalTime() { return _total_time; } - int totalSteps() { return _total_steps; } - - protected: - double _max_dt; - double _cfl_number; - double _beta; - bool _fixed_dt; - int _total_steps; - double _total_time; - double _dt; - double _cell_size; - - p2g_u_op _p2g_momentum; - compute_u_op _compute_velocity; - compute_s_op _compute_stress; - compute_du_s_op _compute_delta_u_s; - compute_du_g_op _compute_delta_u_g; - bc_u_op _apply_bc_momentum; - update_u_op _update_velocity; - g2p_u_op _g2p_momentum; - - property_type _properties; -}; - -template -auto createExplicitMomentumIntegrator( - const std::shared_ptr& mesh, InterpolationType, ParticleVelocity, - StressProperties properties, const double max_dt, const double cfl = 0.5, - const double beta = 1.0, const bool fixed_dt = false ) -{ - return ExplicitMomentumIntegrator( - mesh, properties, max_dt, cfl, beta, fixed_dt ); -} -//---------------------------------------------------------------------------// - -} // namespace Picasso diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 34b8f23f..b536a1fe 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -12,6 +12,7 @@ set(HEADERS Picasso_FieldManager.hpp Picasso_FieldTypes.hpp Picasso_GridOperator.hpp + Picasso_GridUpdate.hpp Picasso_InputParser.hpp Picasso_LevelSet.hpp Picasso_LevelSetRedistance.hpp diff --git a/src/Picasso.hpp b/src/Picasso.hpp index a5104c29..df803d1e 100644 --- a/src/Picasso.hpp +++ b/src/Picasso.hpp @@ -23,6 +23,7 @@ #include #include #include +#include #include #include #include diff --git a/src/Picasso_FieldTypes.hpp b/src/Picasso_FieldTypes.hpp index 6a01bb2c..bdf89bd1 100644 --- a/src/Picasso_FieldTypes.hpp +++ b/src/Picasso_FieldTypes.hpp @@ -966,10 +966,6 @@ struct CommRank : Scalar static std::string label() { return "comm_rank"; } }; -//---------------------------------------------------------------------------// - -} // end namespace Field - struct Mass : Field::Scalar { static std::string label() { return "Mass"; } @@ -995,21 +991,11 @@ struct OldU : Field::Vector static std::string label() { return "Old_Velocity"; } }; -struct DeltaUGravity : Field::Vector -{ - static std::string label() { return "velocity_change_from_gravity"; } -}; - struct Stress : Field::Matrix { static std::string label() { return "stress"; } }; -struct DeltaUStress : Field::Vector -{ - static std::string label() { return "velocity_change_from_stress"; } -}; - struct DetDefGrad : Field::Scalar { static std::string label() { return "Det_deformation_gradient"; } @@ -1017,6 +1003,10 @@ struct DetDefGrad : Field::Scalar using Position = Picasso::Field::LogicalPosition<3>; +//---------------------------------------------------------------------------// + +} // end namespace Field + namespace PolyPIC { namespace Field diff --git a/src/Picasso_GridUpdate.hpp b/src/Picasso_GridUpdate.hpp new file mode 100644 index 00000000..5ae64e28 --- /dev/null +++ b/src/Picasso_GridUpdate.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_GRIDUPDATE_HPP +#define PICASSO_GRIDUPDATE_HPP + +#include +#include + +#include + +#include + +namespace Picasso +{ +//---------------------------------------------------------------------------// +// Compute nodal velocity from mass-weighted momentum. +//---------------------------------------------------------------------------// +template +struct ComputeGridVelocity +{ + template + 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_GRIDUPDATE_HPP diff --git a/src/Picasso_ParticleInterpolation.hpp b/src/Picasso_ParticleInterpolation.hpp index 3f1e0a4e..2877efa2 100644 --- a/src/Picasso_ParticleInterpolation.hpp +++ b/src/Picasso_ParticleInterpolation.hpp @@ -306,11 +306,12 @@ struct Particle2Grid struct Grid2ParticleVelocity { - // FIXME: only needed for FLIP/PIC - double beta; - // Explicit time step. double dt; @@ -433,12 +433,13 @@ struct Grid2ParticleVelocity { // Get particle data. auto u_p = Picasso::get( particle, PolyPIC::Field::Velocity() ); - auto x_p = Picasso::get( particle, Position() ); + auto x_p = Picasso::get( particle, Field::Position() ); // Get the gather dependencies. - auto m_i = gather_deps.get( Picasso::FieldLocation::Node(), Mass() ); - auto u_i = - gather_deps.get( Picasso::FieldLocation::Node(), Velocity() ); + auto m_i = + gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), + Field::Velocity() ); // Get the local dependencies for getting physcial location of node auto x_i = local_deps.get( Picasso::FieldLocation::Node(), Picasso::Field::PhysicalPosition<3>() ); @@ -466,9 +467,6 @@ struct Grid2ParticleVelocity template struct Grid2ParticleVelocity { - // FIXME: only needed for FLIP/PIC - double beta; - // Explicit time step. double dt; @@ -483,12 +481,13 @@ struct Grid2ParticleVelocity { // Get particle data. auto f_p = Picasso::get( particle, APIC::Field::Velocity() ); - auto x_p = Picasso::get( particle, Position() ); + auto x_p = Picasso::get( particle, Field::Position() ); // Get the gather dependencies. - auto m_i = gather_deps.get( Picasso::FieldLocation::Node(), Mass() ); - auto u_i = - gather_deps.get( Picasso::FieldLocation::Node(), Velocity() ); + auto m_i = + gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), + Field::Velocity() ); // Get the local dependencies for getting physcial location of node auto x_i = local_deps.get( Picasso::FieldLocation::Node(), Picasso::Field::PhysicalPosition<3>() ); @@ -524,7 +523,8 @@ struct Grid2ParticleVelocity template + class ParticleViewType, class PositionType = Field::Position, + class VelocityType = Field::Velocity> KOKKOS_INLINE_FUNCTION void operator()( const LocalMeshType& local_mesh, const GatherDependencies& gather_deps, @@ -532,15 +532,16 @@ struct Grid2ParticleVelocity ParticleViewType& particle ) const { // Get particle data. - auto u_p = Picasso::get( particle, Velocity() ); - auto x_p = Picasso::get( particle, Position() ); + auto u_p = Picasso::get( particle, VelocityType{} ); + auto x_p = Picasso::get( particle, PositionType{} ); // Get the gather dependencies. - auto m_i = gather_deps.get( Picasso::FieldLocation::Node(), Mass() ); - auto u_i = - gather_deps.get( Picasso::FieldLocation::Node(), Velocity() ); + auto m_i = + gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), + Field::Velocity() ); auto old_u_i = - gather_deps.get( Picasso::FieldLocation::Node(), OldU() ); + gather_deps.get( Picasso::FieldLocation::Node(), Field::OldU() ); // Get the local dependencies for getting physcial location of node auto x_i = local_deps.get( Picasso::FieldLocation::Node(), From b27b83bff95286e1ccf4e99d3d34150427b8aab4 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 14 Oct 2024 17:04:52 -0400 Subject: [PATCH 12/22] fixup: constructors for p2g/g2p functors --- examples/dam_break.cpp | 2 +- src/Picasso_ParticleInterpolation.hpp | 81 ++++++++++++++++++++++----- 2 files changed, 69 insertions(+), 14 deletions(-) diff --git a/examples/dam_break.cpp b/examples/dam_break.cpp index 36fa163a..29dca7d6 100644 --- a/examples/dam_break.cpp +++ b/examples/dam_break.cpp @@ -379,7 +379,7 @@ void DamBreak( std::string filename ) // Grid interpolation (Picasso built-in). Picasso::Grid2ParticleVelocity - g2p_func{ beta, dt }; + g2p_func{ dt, beta }; g2p.apply( "g2p_U", particle_type{}, exec_space{}, *fm, particles, g2p_func ); diff --git a/src/Picasso_ParticleInterpolation.hpp b/src/Picasso_ParticleInterpolation.hpp index 2877efa2..205c9c10 100644 --- a/src/Picasso_ParticleInterpolation.hpp +++ b/src/Picasso_ParticleInterpolation.hpp @@ -293,7 +293,13 @@ struct Particle2Grid { // Explicit time step. - double dt; + double _dt; + + // Constructor + Particle2Grid( const double dt ) + : _dt( dt ) + { + } template { // Explicit time step. - double dt; + double _dt; + + // Constructor + Particle2Grid( const double dt ) + : _dt( dt ) + { + } template { // Explicit time step. - double dt; + double _dt; + + // Constructor + Particle2Grid( const double dt ) + : _dt( dt ) + { + } template struct Grid2ParticleVelocity { // Explicit time step. - double dt; + 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 // 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 ); }; + { return x_i( i, j, k, d ) + _dt * u_i( i, j, k, d ); }; Picasso::G2P::value( spline, x_i_updated, x_p ); } }; @@ -468,7 +498,19 @@ template struct Grid2ParticleVelocity { // Explicit time step. - double dt; + 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 // 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 ); }; + { return x_i( i, j, k, d ) + _dt * u_i( i, j, k, d ); }; Picasso::G2P::value( spline, x_i_updated, x_p ); } }; @@ -516,10 +558,23 @@ struct Grid2ParticleVelocity template struct Grid2ParticleVelocity { - double beta = 0.01; - // Explicit time step. - double dt; + 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 u_p_flip = u_p + d_u_p; // Update particle velocity. - u_p = beta * u_p_flip + ( 1.0 - beta ) * u_p_pic; + 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 ); }; + { return x_i( i, j, k, d ) + _dt * u_i( i, j, k, d ); }; Picasso::G2P::value( spline, x_i_updated, x_p ); } }; From 8a04c8efbaad908c6ebae537fc5fc164fbc9da46 Mon Sep 17 00:00:00 2001 From: "Chong, Kwitae" Date: Tue, 15 Oct 2024 10:59:35 -0400 Subject: [PATCH 13/22] add globalConservation --- examples/dam_break.cpp | 9 ++- src/Picasso.hpp | 1 + src/Picasso_FieldTypes.hpp | 31 ++++++++ src/Picasso_PostProcess.hpp | 140 ++++++++++++++++++++++++++++++++++++ 4 files changed, 179 insertions(+), 2 deletions(-) create mode 100644 src/Picasso_PostProcess.hpp diff --git a/examples/dam_break.cpp b/examples/dam_break.cpp index 29dca7d6..15fe54ae 100644 --- a/examples/dam_break.cpp +++ b/examples/dam_break.cpp @@ -398,6 +398,11 @@ void DamBreak( std::string filename ) particles.slice( Picasso::Field::Mass() ), particles.slice( Picasso::Field::Volume() ) ); + // Write conservation sums + Picasso::PostProcess::globalConservation( + MPI_COMM_WORLD, exec_space(), mesh, *fm, + steps, write_frequency, particles, ParticleVelocity() ); + time += dt; steps++; } @@ -417,8 +422,8 @@ int main( int argc, char* argv[] ) std::string filename = argv[1]; // Problem can run with any interpolation scheme. - // DamBreak(); - // DamBreak(); + // DamBreak( filename); + // DamBreak( filename ); DamBreak( filename ); Kokkos::finalize(); diff --git a/src/Picasso.hpp b/src/Picasso.hpp index df803d1e..d7d063dc 100644 --- a/src/Picasso.hpp +++ b/src/Picasso.hpp @@ -34,6 +34,7 @@ #endif #include #include +#include #include #include #include diff --git a/src/Picasso_FieldTypes.hpp b/src/Picasso_FieldTypes.hpp index bdf89bd1..d7d496a2 100644 --- a/src/Picasso_FieldTypes.hpp +++ b/src/Picasso_FieldTypes.hpp @@ -1045,6 +1045,37 @@ struct FlipTag { }; +// Particle indexing for PIC/FLIP vs PolyPIC scalar field +KOKKOS_INLINE_FUNCTION auto getField( const double field ) { return field; } + +template +KOKKOS_INLINE_FUNCTION + std::enable_if_t::value, + typename FieldType::value_type> + getField( const FieldType& field ) +{ + return field( 0, 0 ); +} + +// Particle indexing for PIC/FLIP vs PolyPIC vector field +template +KOKKOS_INLINE_FUNCTION + std::enable_if_t::value, + typename FieldType::value_type> + getField( const FieldType& field, const int dir ) +{ + return field( 0, dir ); +} + +template +KOKKOS_INLINE_FUNCTION + std::enable_if_t::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_PostProcess.hpp b/src/Picasso_PostProcess.hpp new file mode 100644 index 00000000..dcc35c2a --- /dev/null +++ b/src/Picasso_PostProcess.hpp @@ -0,0 +1,140 @@ +/**************************************************************************** + * 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_POSTPROCESS_HPP +#define PICASSO_VPOSTPROCESSHPP + +#include + +#include + + +namespace Picasso +{ +namespace PostProcess +{ + +// Output routine for integral conservation sums of particle quantities +// (momentum-only) +template +void particleConservation( + MPI_Comm comm, ExecutionSpace, ParticleVelocity, const int n, + const int write_frequency, const ParticleList particles ) +{ + if ( write_frequency > 0 && n % write_frequency == 0 ) + { + std::size_t num_p = particles.size(); + auto aosoa = particles.aosoa(); + + double mass_integral = 0; + double ke_integral = 0; + Kokkos::parallel_reduce( + "Picasso::Particle Conservation", + Kokkos::RangePolicy( 0, num_p ), + KOKKOS_LAMBDA( const int pn, double& mass_cont, double& ke_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 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; + }, + mass_integral, ke_integral ); + + double global_mass = 0.0; + double global_ke = 0.0; + MPI_Allreduce( &mass_integral, &global_mass, 1, MPI_DOUBLE, MPI_SUM, comm ); + MPI_Allreduce( &ke_integral, &global_ke, 1, MPI_DOUBLE, MPI_SUM, comm ); + + int rank; + MPI_Comm_rank( comm, &rank ); + if ( rank == 0 ) + std::cout << "Current Particle Integrals for Step: " << n + << " Mass:" << global_mass + << " Kinetic Energy:" << global_ke << std::endl; + } +} + +// Output for integral conservation sums for grid quantities (momentum-only) +template +void gridConservation( MPI_Comm comm, ExecutionSpace exec_space, + const std::shared_ptr& mesh, + Picasso::FieldManager& fm, const int n, + const int write_frequency, ParticleVelocity ) +{ + if ( write_frequency > 0 && n % write_frequency == 0 ) + { + 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() ); + double mass_integral = 0; + double ke_integral = 0; + + // FIX ME: 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 ); + + double global_mass = 0.0; + double global_ke = 0.0; + MPI_Allreduce( &mass_integral, &global_mass, 1, MPI_DOUBLE, MPI_SUM, comm ); + MPI_Allreduce( &ke_integral, &global_ke, 1, MPI_DOUBLE, MPI_SUM, comm ); + + int rank; + MPI_Comm_rank( comm, &rank ); + if ( rank == 0 ) + std::cout << "Current Grid Integrals for Step: " << n + << " Mass:" << global_mass + << " Kinetic Energy:" << global_ke << std::endl; + } +} + +// Global conservation output for momentum problem +template +void globalConservation( MPI_Comm comm, ExecutionSpace exec_space, + const std::shared_ptr& mesh, + Picasso::FieldManager& fm, const int n, + const int write_frequency, + const ParticleList particles, ParticleVelocity p_v ) +{ + particleConservation( comm, exec_space, p_v, n, write_frequency, + particles ); + + gridConservation( comm, exec_space, mesh, fm, n, write_frequency, p_v ); +} + +} // end namespace PostProcess +} // end namespace Picasso + +#endif // end PICASSO_POSTPROCESS_HPP From 292ab9dd3c707ba1b125e289b4988dbe2c8d46c2 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 21 Oct 2024 11:32:00 -0400 Subject: [PATCH 14/22] Rearrange conservation files --- examples/dam_break.cpp | 27 +++++-- src/Picasso.hpp | 2 +- src/Picasso_Conservation.hpp | 94 +++++++++++++++++++++++ src/Picasso_PostProcess.hpp | 140 ----------------------------------- 4 files changed, 115 insertions(+), 148 deletions(-) create mode 100644 src/Picasso_Conservation.hpp delete mode 100644 src/Picasso_PostProcess.hpp diff --git a/examples/dam_break.cpp b/examples/dam_break.cpp index 15fe54ae..c6e6c595 100644 --- a/examples/dam_break.cpp +++ b/examples/dam_break.cpp @@ -389,20 +389,33 @@ void DamBreak( std::string filename ) // Write particle fields. Cabana::Experimental::HDF5ParticleOutput::HDF5Config h5_config; if ( steps % write_frequency == 0 ) + { Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( - h5_config, "particles", MPI_COMM_WORLD, steps / write_frequency, - time, particles.size(), + h5_config, "particles", global_grid.comm(), + steps / write_frequency, time, particles.size(), particles.slice( Picasso::Field::Position() ), particles.slice( Picasso::Field::Pressure() ), particles.slice( ParticleVelocity() ), particles.slice( Picasso::Field::Mass() ), particles.slice( Picasso::Field::Volume() ) ); - // Write conservation sums - Picasso::PostProcess::globalConservation( - MPI_COMM_WORLD, exec_space(), mesh, *fm, - steps, write_frequency, particles, ParticleVelocity() ); - + // Calculate conservation sums + double mass_particles = 0.0; + double ke_particles = 0.0; + Picasso::particleConservation( global_grid.comm(), exec_space(), + particles, ParticleVelocity(), + mass_particles, ke_particles ); + double mass_grid = 0.0; + double ke_grid = 0.0; + Picasso::gridConservation( global_grid.comm(), exec_space(), mesh, + *fm, mass_grid, ke_grid ); + + if ( global_grid.blockId() == 0 ) + std::cout << "Particle/Grid Mass: " << mass_particles << " / " + << mass_grid << "\n" + << "Particle/Grid Kinetic Energy: " << ke_particles + << " / " << ke_grid << "\n\n"; + } time += dt; steps++; } diff --git a/src/Picasso.hpp b/src/Picasso.hpp index d7d063dc..9a8de3b9 100644 --- a/src/Picasso.hpp +++ b/src/Picasso.hpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -34,7 +35,6 @@ #endif #include #include -#include #include #include #include diff --git a/src/Picasso_Conservation.hpp b/src/Picasso_Conservation.hpp new file mode 100644 index 00000000..476dbf9c --- /dev/null +++ b/src/Picasso_Conservation.hpp @@ -0,0 +1,94 @@ +/**************************************************************************** + * 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 + +#include + +namespace Picasso +{ +// Output routine for integral conservation sums of particle quantities +// (momentum-only) +template +void particleConservation( MPI_Comm comm, ExecutionSpace, + const ParticleList particles, ParticleVelocity, + double& global_mass, double& global_ke ) +{ + std::size_t num_p = particles.size(); + auto aosoa = particles.aosoa(); + + double mass_integral = 0; + double ke_integral = 0; + Kokkos::parallel_reduce( + "Picasso::Particle Conservation", + Kokkos::RangePolicy( 0, num_p ), + KOKKOS_LAMBDA( const int pn, double& mass_cont, double& ke_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 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; + }, + mass_integral, ke_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 ); +} + +// Output for integral conservation sums for grid quantities (momentum-only) +template +void gridConservation( MPI_Comm comm, ExecutionSpace exec_space, + const std::shared_ptr& mesh, + Picasso::FieldManager& fm, double& global_mass, + double& global_ke ) +{ + 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() ); + double mass_integral = 0; + double ke_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 ); + + MPI_Allreduce( &mass_integral, &global_mass, 1, MPI_DOUBLE, MPI_SUM, comm ); + MPI_Allreduce( &ke_integral, &global_ke, 1, MPI_DOUBLE, MPI_SUM, comm ); +} + +} // end namespace Picasso + +#endif // end PICASSO_CONSERVATION_HPP diff --git a/src/Picasso_PostProcess.hpp b/src/Picasso_PostProcess.hpp deleted file mode 100644 index dcc35c2a..00000000 --- a/src/Picasso_PostProcess.hpp +++ /dev/null @@ -1,140 +0,0 @@ -/**************************************************************************** - * 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_POSTPROCESS_HPP -#define PICASSO_VPOSTPROCESSHPP - -#include - -#include - - -namespace Picasso -{ -namespace PostProcess -{ - -// Output routine for integral conservation sums of particle quantities -// (momentum-only) -template -void particleConservation( - MPI_Comm comm, ExecutionSpace, ParticleVelocity, const int n, - const int write_frequency, const ParticleList particles ) -{ - if ( write_frequency > 0 && n % write_frequency == 0 ) - { - std::size_t num_p = particles.size(); - auto aosoa = particles.aosoa(); - - double mass_integral = 0; - double ke_integral = 0; - Kokkos::parallel_reduce( - "Picasso::Particle Conservation", - Kokkos::RangePolicy( 0, num_p ), - KOKKOS_LAMBDA( const int pn, double& mass_cont, double& ke_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 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; - }, - mass_integral, ke_integral ); - - double global_mass = 0.0; - double global_ke = 0.0; - MPI_Allreduce( &mass_integral, &global_mass, 1, MPI_DOUBLE, MPI_SUM, comm ); - MPI_Allreduce( &ke_integral, &global_ke, 1, MPI_DOUBLE, MPI_SUM, comm ); - - int rank; - MPI_Comm_rank( comm, &rank ); - if ( rank == 0 ) - std::cout << "Current Particle Integrals for Step: " << n - << " Mass:" << global_mass - << " Kinetic Energy:" << global_ke << std::endl; - } -} - -// Output for integral conservation sums for grid quantities (momentum-only) -template -void gridConservation( MPI_Comm comm, ExecutionSpace exec_space, - const std::shared_ptr& mesh, - Picasso::FieldManager& fm, const int n, - const int write_frequency, ParticleVelocity ) -{ - if ( write_frequency > 0 && n % write_frequency == 0 ) - { - 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() ); - double mass_integral = 0; - double ke_integral = 0; - - // FIX ME: 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 ); - - double global_mass = 0.0; - double global_ke = 0.0; - MPI_Allreduce( &mass_integral, &global_mass, 1, MPI_DOUBLE, MPI_SUM, comm ); - MPI_Allreduce( &ke_integral, &global_ke, 1, MPI_DOUBLE, MPI_SUM, comm ); - - int rank; - MPI_Comm_rank( comm, &rank ); - if ( rank == 0 ) - std::cout << "Current Grid Integrals for Step: " << n - << " Mass:" << global_mass - << " Kinetic Energy:" << global_ke << std::endl; - } -} - -// Global conservation output for momentum problem -template -void globalConservation( MPI_Comm comm, ExecutionSpace exec_space, - const std::shared_ptr& mesh, - Picasso::FieldManager& fm, const int n, - const int write_frequency, - const ParticleList particles, ParticleVelocity p_v ) -{ - particleConservation( comm, exec_space, p_v, n, write_frequency, - particles ); - - gridConservation( comm, exec_space, mesh, fm, n, write_frequency, p_v ); -} - -} // end namespace PostProcess -} // end namespace Picasso - -#endif // end PICASSO_POSTPROCESS_HPP From 76e85d565ba881e38bf350e3307f1fd1d7f9e039 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 21 Oct 2024 11:44:04 -0400 Subject: [PATCH 15/22] Separate step kernels --- src/Picasso.hpp | 3 +- ...GridUpdate.hpp => Picasso_GridKernels.hpp} | 6 +- src/Picasso_InterpolationKernels.hpp | 389 ++++++++++++++++++ src/Picasso_ParticleInterpolation.hpp | 363 ---------------- 4 files changed, 394 insertions(+), 367 deletions(-) rename src/{Picasso_GridUpdate.hpp => Picasso_GridKernels.hpp} (95%) create mode 100644 src/Picasso_InterpolationKernels.hpp diff --git a/src/Picasso.hpp b/src/Picasso.hpp index 9a8de3b9..9a21c2bf 100644 --- a/src/Picasso.hpp +++ b/src/Picasso.hpp @@ -23,9 +23,10 @@ #include #include #include +#include #include -#include #include +#include #include #include #include diff --git a/src/Picasso_GridUpdate.hpp b/src/Picasso_GridKernels.hpp similarity index 95% rename from src/Picasso_GridUpdate.hpp rename to src/Picasso_GridKernels.hpp index 5ae64e28..76c52306 100644 --- a/src/Picasso_GridUpdate.hpp +++ b/src/Picasso_GridKernels.hpp @@ -9,8 +9,8 @@ * SPDX-License-Identifier: BSD-3-Clause * ****************************************************************************/ -#ifndef PICASSO_GRIDUPDATE_HPP -#define PICASSO_GRIDUPDATE_HPP +#ifndef PICASSO_GRIDKERNELS_HPP +#define PICASSO_GRIDKERNELS_HPP #include #include @@ -54,4 +54,4 @@ struct ComputeGridVelocity }; } // namespace Picasso -#endif // end PICASSO_GRIDUPDATE_HPP +#endif // end PICASSO_GRIDKERNELS_HPP diff --git a/src/Picasso_InterpolationKernels.hpp b/src/Picasso_InterpolationKernels.hpp new file mode 100644 index 00000000..854755ed --- /dev/null +++ b/src/Picasso_InterpolationKernels.hpp @@ -0,0 +1,389 @@ +/**************************************************************************** + * 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 +#include +#include +#include + +#include + +#include + +namespace Picasso +{ +//---------------------------------------------------------------------------// +template +struct Particle2Grid; + +//---------------------------------------------------------------------------// +// Project particle enthalpy/momentum to grid. PolyPIC variant +//---------------------------------------------------------------------------// +template +struct Particle2Grid +{ + // Explicit time step. + double _dt; + + // Constructor + Particle2Grid( const double dt ) + : _dt( dt ) + { + } + + template + 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 v_p = Picasso::get( particle, PolyPIC::Field::Velocity() ); + auto m_p = Picasso::get( particle, Field::Mass() ); + auto x_p = Picasso::get( particle, Field::Position() ); + + // Get the scatter dependencies. + auto m_i = + scatter_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + auto f_i = + scatter_deps.get( Picasso::FieldLocation::Node(), OldFieldType() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::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 +struct Particle2Grid +{ + // Explicit time step. + double _dt; + + // Constructor + Particle2Grid( const double dt ) + : _dt( dt ) + { + } + + template + 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, Field::Mass() ); + auto x_p = Picasso::get( particle, Field::Position() ); + + // Get the scatter dependencies. + auto m_i = + scatter_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + auto f_i = + scatter_deps.get( Picasso::FieldLocation::Node(), OldFieldType() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::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 +struct Particle2Grid +{ + // Explicit time step. + double _dt; + + // Constructor + Particle2Grid( const double dt ) + : _dt( dt ) + { + } + + template + 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, Field::Mass() ); + auto x_p = Picasso::get( particle, Field::Position() ); + + // Get the scatter dependencies. + auto m_i = + scatter_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + auto f_i = + scatter_deps.get( Picasso::FieldLocation::Node(), OldFieldType() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::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 ); + } +}; + +template +struct Grid2ParticleVelocity; + +//---------------------------------------------------------------------------// +// Update particle state. PolyPIC variant. +//---------------------------------------------------------------------------// +template +struct Grid2ParticleVelocity +{ + // 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 + 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() ); + auto x_p = Picasso::get( particle, Field::Position() ); + + // Get the gather dependencies. + auto m_i = + gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), + Field::Velocity() ); + // Get the local dependencies for getting physcial location of node + auto x_i = local_deps.get( Picasso::FieldLocation::Node(), + Picasso::Field::PhysicalPosition<3>() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::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 +struct Grid2ParticleVelocity +{ + // 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 + 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() ); + auto x_p = Picasso::get( particle, Field::Position() ); + + // Get the gather dependencies. + auto m_i = + gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), + Field::Velocity() ); + // Get the local dependencies for getting physcial location of node + auto x_i = local_deps.get( Picasso::FieldLocation::Node(), + Picasso::Field::PhysicalPosition<3>() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::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 +struct Grid2ParticleVelocity +{ + // 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 + 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, VelocityType{} ); + auto x_p = Picasso::get( particle, PositionType{} ); + + // Get the gather dependencies. + auto m_i = + gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); + auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), + Field::Velocity() ); + auto old_u_i = + gather_deps.get( Picasso::FieldLocation::Node(), Field::OldU() ); + + // Get the local dependencies for getting physcial location of node + auto x_i = local_deps.get( Picasso::FieldLocation::Node(), + Picasso::Field::PhysicalPosition<3>() ); + + // Node interpolant. + auto spline = Picasso::createSpline( + Picasso::FieldLocation::Node(), + Picasso::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 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 d_u_p; + Picasso::G2P::value( spline, d_u_i, d_u_p ); + + Picasso::Vec3 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 ); + } +}; + +} // end namespace Picasso + +#endif // end PICASSO_INTERPOLATIONKERNELS_HPP diff --git a/src/Picasso_ParticleInterpolation.hpp b/src/Picasso_ParticleInterpolation.hpp index 205c9c10..b8cf5a8c 100644 --- a/src/Picasso_ParticleInterpolation.hpp +++ b/src/Picasso_ParticleInterpolation.hpp @@ -12,10 +12,8 @@ #ifndef PICASSO_PARTICLEINTERPOLATION_HPP #define PICASSO_PARTICLEINTERPOLATION_HPP -#include #include #include -#include #include @@ -280,367 +278,6 @@ KOKKOS_INLINE_FUNCTION void divergence( const SplineDataType& sd, } // end namespace P2G -//---------------------------------------------------------------------------// -template -struct Particle2Grid; - -//---------------------------------------------------------------------------// -// Project particle enthalpy/momentum to grid. PolyPIC variant -//---------------------------------------------------------------------------// -template -struct Particle2Grid -{ - // Explicit time step. - double _dt; - - // Constructor - Particle2Grid( const double dt ) - : _dt( dt ) - { - } - - template - 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 v_p = Picasso::get( particle, PolyPIC::Field::Velocity() ); - auto m_p = Picasso::get( particle, Field::Mass() ); - auto x_p = Picasso::get( particle, Field::Position() ); - - // Get the scatter dependencies. - auto m_i = - scatter_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); - auto f_i = - scatter_deps.get( Picasso::FieldLocation::Node(), OldFieldType() ); - - // Node interpolant. - auto spline = Picasso::createSpline( - Picasso::FieldLocation::Node(), - Picasso::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 -struct Particle2Grid -{ - // Explicit time step. - double _dt; - - // Constructor - Particle2Grid( const double dt ) - : _dt( dt ) - { - } - - template - 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, Field::Mass() ); - auto x_p = Picasso::get( particle, Field::Position() ); - - // Get the scatter dependencies. - auto m_i = - scatter_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); - auto f_i = - scatter_deps.get( Picasso::FieldLocation::Node(), OldFieldType() ); - - // Node interpolant. - auto spline = Picasso::createSpline( - Picasso::FieldLocation::Node(), - Picasso::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 -struct Particle2Grid -{ - // Explicit time step. - double _dt; - - // Constructor - Particle2Grid( const double dt ) - : _dt( dt ) - { - } - - template - 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, Field::Mass() ); - auto x_p = Picasso::get( particle, Field::Position() ); - - // Get the scatter dependencies. - auto m_i = - scatter_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); - auto f_i = - scatter_deps.get( Picasso::FieldLocation::Node(), OldFieldType() ); - - // Node interpolant. - auto spline = Picasso::createSpline( - Picasso::FieldLocation::Node(), - Picasso::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 ); - } -}; - -template -struct Grid2ParticleVelocity; - -//---------------------------------------------------------------------------// -// Update particle state. PolyPIC variant. -//---------------------------------------------------------------------------// -template -struct Grid2ParticleVelocity -{ - // 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 - 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() ); - auto x_p = Picasso::get( particle, Field::Position() ); - - // Get the gather dependencies. - auto m_i = - gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); - auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Field::Velocity() ); - // Get the local dependencies for getting physcial location of node - auto x_i = local_deps.get( Picasso::FieldLocation::Node(), - Picasso::Field::PhysicalPosition<3>() ); - - // Node interpolant. - auto spline = Picasso::createSpline( - Picasso::FieldLocation::Node(), - Picasso::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 -struct Grid2ParticleVelocity -{ - // 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 - 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() ); - auto x_p = Picasso::get( particle, Field::Position() ); - - // Get the gather dependencies. - auto m_i = - gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); - auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Field::Velocity() ); - // Get the local dependencies for getting physcial location of node - auto x_i = local_deps.get( Picasso::FieldLocation::Node(), - Picasso::Field::PhysicalPosition<3>() ); - - // Node interpolant. - auto spline = Picasso::createSpline( - Picasso::FieldLocation::Node(), - Picasso::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 -struct Grid2ParticleVelocity -{ - // 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 - 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, VelocityType{} ); - auto x_p = Picasso::get( particle, PositionType{} ); - - // Get the gather dependencies. - auto m_i = - gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); - auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Field::Velocity() ); - auto old_u_i = - gather_deps.get( Picasso::FieldLocation::Node(), Field::OldU() ); - - // Get the local dependencies for getting physcial location of node - auto x_i = local_deps.get( Picasso::FieldLocation::Node(), - Picasso::Field::PhysicalPosition<3>() ); - - // Node interpolant. - auto spline = Picasso::createSpline( - Picasso::FieldLocation::Node(), - Picasso::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 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 d_u_p; - Picasso::G2P::value( spline, d_u_i, d_u_p ); - - Picasso::Vec3 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 ); - } -}; - } // end namespace Picasso #endif // end PICASSO_PARTICLEINTERPOLATION_HPP From 282a6a3ef5479eff9c1b920d215adda76ef38dd9 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 21 Oct 2024 12:18:07 -0400 Subject: [PATCH 16/22] fixup: field names --- examples/dam_break.cpp | 3 ++- src/Picasso_FieldTypes.hpp | 16 +++++++--------- src/Picasso_GridKernels.hpp | 2 +- src/Picasso_InterpolationKernels.hpp | 4 ++-- 4 files changed, 12 insertions(+), 13 deletions(-) diff --git a/examples/dam_break.cpp b/examples/dam_break.cpp index c6e6c595..c4413586 100644 --- a/examples/dam_break.cpp +++ b/examples/dam_break.cpp @@ -305,7 +305,8 @@ void DamBreak( std::string filename ) using mass_type = Picasso::FieldLayout; using velocity_type = Picasso::FieldLayout; - using old_u_type = Picasso::FieldLayout; + using old_u_type = + Picasso::FieldLayout; using delta_u_s_type = Picasso::FieldLayout; using delta_u_g_type = Picasso::FieldLayout; diff --git a/src/Picasso_FieldTypes.hpp b/src/Picasso_FieldTypes.hpp index d7d496a2..55400d4b 100644 --- a/src/Picasso_FieldTypes.hpp +++ b/src/Picasso_FieldTypes.hpp @@ -968,27 +968,27 @@ struct CommRank : Scalar struct Mass : Field::Scalar { - static std::string label() { return "Mass"; } + static std::string label() { return "mass"; } }; struct Pressure : Field::Scalar { - static std::string label() { return "Pressure"; } + static std::string label() { return "pressure"; } }; struct Volume : Field::Scalar { - static std::string label() { return "Volume"; } + static std::string label() { return "volume"; } }; struct Velocity : Field::Vector { - static std::string label() { return "Velocity"; } + static std::string label() { return "velocity"; } }; -struct OldU : Field::Vector +struct OldVelocity : Field::Vector { - static std::string label() { return "Old_Velocity"; } + static std::string label() { return "old_velocity"; } }; struct Stress : Field::Matrix @@ -998,11 +998,9 @@ struct Stress : Field::Matrix struct DetDefGrad : Field::Scalar { - static std::string label() { return "Det_deformation_gradient"; } + static std::string label() { return "determinant_deformation_gradient"; } }; -using Position = Picasso::Field::LogicalPosition<3>; - //---------------------------------------------------------------------------// } // end namespace Field diff --git a/src/Picasso_GridKernels.hpp b/src/Picasso_GridKernels.hpp index 76c52306..05921d3e 100644 --- a/src/Picasso_GridKernels.hpp +++ b/src/Picasso_GridKernels.hpp @@ -25,7 +25,7 @@ namespace Picasso // Compute nodal velocity from mass-weighted momentum. //---------------------------------------------------------------------------// template + class OldVelType = Field::OldVelocity> struct ComputeGridVelocity { template gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), Field::Velocity() ); - auto old_u_i = - gather_deps.get( Picasso::FieldLocation::Node(), Field::OldU() ); + auto old_u_i = gather_deps.get( Picasso::FieldLocation::Node(), + Field::OldVelocity() ); // Get the local dependencies for getting physcial location of node auto x_i = local_deps.get( Picasso::FieldLocation::Node(), From 05e168131e4e4dea67ce3a75bec27f17a7a8623a Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 21 Oct 2024 13:16:50 -0400 Subject: [PATCH 17/22] fixup: defaults for interpolation kernels --- src/Picasso_InterpolationKernels.hpp | 132 ++++++++++++++++++--------- 1 file changed, 89 insertions(+), 43 deletions(-) diff --git a/src/Picasso_InterpolationKernels.hpp b/src/Picasso_InterpolationKernels.hpp index 304be7a7..a18488bb 100644 --- a/src/Picasso_InterpolationKernels.hpp +++ b/src/Picasso_InterpolationKernels.hpp @@ -15,6 +15,7 @@ #include #include #include +#include #include #include @@ -23,17 +24,22 @@ namespace Picasso { +//---------------------------------------------------------------------------// +// P2G +//---------------------------------------------------------------------------// + //---------------------------------------------------------------------------// template + class MassType, class PositionType, class InterpolationType> struct Particle2Grid; //---------------------------------------------------------------------------// // Project particle enthalpy/momentum to grid. PolyPIC variant //---------------------------------------------------------------------------// -template +template struct Particle2Grid + MassType, PositionType, PolyPicTag> { // Explicit time step. double _dt; @@ -54,13 +60,14 @@ struct Particle2Grid +template struct Particle2Grid + MassType, PositionType, APicTag> { // Explicit time step. double _dt; @@ -102,12 +110,12 @@ struct Particle2Grid +template struct Particle2Grid + MassType, PositionType, FlipTag> { // Explicit time step. double _dt; @@ -150,12 +159,12 @@ struct Particle2Grid +// P2G creation function (used instead of direct class defaults because of +// partial specializations). +template > +auto createParticle2Grid( const double dt ) +{ + return Particle2Grid( dt ); +} + +//---------------------------------------------------------------------------// +// G2P +//---------------------------------------------------------------------------// + +template struct Grid2ParticleVelocity; //---------------------------------------------------------------------------// // Update particle state. PolyPIC variant. //---------------------------------------------------------------------------// -template -struct Grid2ParticleVelocity +template +struct Grid2ParticleVelocity { // Explicit time step. double _dt; @@ -206,16 +235,16 @@ struct Grid2ParticleVelocity { // Get particle data. auto u_p = Picasso::get( particle, PolyPIC::Field::Velocity() ); - auto x_p = Picasso::get( particle, Field::Position() ); + auto x_p = Picasso::get( particle, PositionType() ); // Get the gather dependencies. auto m_i = - gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); - auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Field::Velocity() ); + 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(), - Picasso::Field::PhysicalPosition<3>() ); + auto x_i = + local_deps.get( Picasso::FieldLocation::Node(), GridPosition() ); // Node interpolant. auto spline = Picasso::createSpline( @@ -237,8 +266,10 @@ struct Grid2ParticleVelocity //---------------------------------------------------------------------------// // Update particle state. APIC variant. //---------------------------------------------------------------------------// -template -struct Grid2ParticleVelocity +template +struct Grid2ParticleVelocity { // Explicit time step. double _dt; @@ -266,16 +297,16 @@ struct Grid2ParticleVelocity { // Get particle data. auto f_p = Picasso::get( particle, APIC::Field::Velocity() ); - auto x_p = Picasso::get( particle, Field::Position() ); + auto x_p = Picasso::get( particle, PositionType() ); // Get the gather dependencies. auto m_i = - gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); - auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Field::Velocity() ); + 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(), - Picasso::Field::PhysicalPosition<3>() ); + auto x_i = + local_deps.get( Picasso::FieldLocation::Node(), GridPosition() ); // Node interpolant. auto spline = Picasso::createSpline( @@ -298,8 +329,10 @@ struct Grid2ParticleVelocity //---------------------------------------------------------------------------// // Update particle state. FLIP/PIC variant. //---------------------------------------------------------------------------// -template -struct Grid2ParticleVelocity +template +struct Grid2ParticleVelocity { // Explicit time step. double _dt; @@ -321,8 +354,7 @@ struct Grid2ParticleVelocity template + class ParticleViewType> KOKKOS_INLINE_FUNCTION void operator()( const LocalMeshType& local_mesh, const GatherDependencies& gather_deps, @@ -330,20 +362,20 @@ struct Grid2ParticleVelocity ParticleViewType& particle ) const { // Get particle data. - auto u_p = Picasso::get( particle, VelocityType{} ); + auto u_p = Picasso::get( particle, Field::Velocity{} ); auto x_p = Picasso::get( particle, PositionType{} ); // Get the gather dependencies. auto m_i = - gather_deps.get( Picasso::FieldLocation::Node(), Field::Mass() ); - auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Field::Velocity() ); + 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() ); // Get the local dependencies for getting physcial location of node - auto x_i = local_deps.get( Picasso::FieldLocation::Node(), - Picasso::Field::PhysicalPosition<3>() ); + auto x_i = + local_deps.get( Picasso::FieldLocation::Node(), GridPosition() ); // Node interpolant. auto spline = Picasso::createSpline( @@ -384,6 +416,20 @@ struct Grid2ParticleVelocity } }; +// G2P creation function (used instead of direct class defaults because of +// partial specializations). +template , + class GridVelocity = Picasso::Field::Velocity, + class GridPosition = Picasso::Field::LogicalPosition<3>> +auto createGrid2ParticleVelocity( const double dt, const double beta ) +{ + return Grid2ParticleVelocity( + dt, beta ); +} + } // end namespace Picasso #endif // end PICASSO_INTERPOLATIONKERNELS_HPP From bd6b3828874499e7ae5b00c122048da6d621741f Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 21 Oct 2024 13:17:15 -0400 Subject: [PATCH 18/22] fixup: build --- src/CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b536a1fe..8e3866c0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -12,8 +12,9 @@ set(HEADERS Picasso_FieldManager.hpp Picasso_FieldTypes.hpp Picasso_GridOperator.hpp - Picasso_GridUpdate.hpp + Picasso_GridKernels.hpp Picasso_InputParser.hpp + Picasso_InterpolationKernels.hpp Picasso_LevelSet.hpp Picasso_LevelSetRedistance.hpp Picasso_MarchingCubes.hpp From e8f588732bac024d36e00d7c2e06ebfca6c40ecd Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Mon, 21 Oct 2024 13:19:19 -0400 Subject: [PATCH 19/22] fixup: template dimension for all fields --- examples/dam_break.cpp | 53 +++++++++++++++------------- src/Picasso_Conservation.hpp | 2 +- src/Picasso_FieldTypes.hpp | 15 +++++--- src/Picasso_GridKernels.hpp | 4 +-- src/Picasso_InterpolationKernels.hpp | 18 +++++----- 5 files changed, 51 insertions(+), 41 deletions(-) diff --git a/examples/dam_break.cpp b/examples/dam_break.cpp index c4413586..cb967755 100644 --- a/examples/dam_break.cpp +++ b/examples/dam_break.cpp @@ -24,7 +24,7 @@ struct ParticleInitFunc x[1] <= block[4] && block[2] <= x[2] && x[2] <= block[5] ) { - Picasso::get( p, Picasso::Field::Stress() ) = 0.0; + 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; @@ -32,7 +32,8 @@ struct ParticleInitFunc Picasso::get( p, Picasso::Field::Pressure() ) = 0.0; for ( int d = 0; d < 3; ++d ) - Picasso::get( p, Picasso::Field::Position(), d ) = x[d]; + Picasso::get( p, Picasso::Field::LogicalPosition<3>(), d ) = + x[d]; return true; } @@ -102,16 +103,17 @@ struct ComputeParticlePressure ParticleViewType& particle ) const { // Get particle data. - auto x_p = Picasso::get( particle, Picasso::Field::Position() ); + 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() ); + 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() ); + Picasso::Field::Velocity<3>() ); // update strain rate auto spline = Picasso::createSpline( @@ -158,8 +160,9 @@ struct ComputeGridVelocityChange // 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::Position() ); - auto s_p = Picasso::get( particle, Picasso::Field::Stress() ); + 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 = @@ -200,7 +203,7 @@ struct UpdateGridVelocity auto m_i = gather_deps.get( Picasso::FieldLocation::Node(), Picasso::Field::Mass() ); auto u_i = gather_deps.get( Picasso::FieldLocation::Node(), - Picasso::Field::Velocity() ); + Picasso::Field::Velocity<3>() ); auto delta_u_s_i = gather_deps.get( Picasso::FieldLocation::Node(), DeltaUStress() ); @@ -241,10 +244,10 @@ void DamBreak( std::string filename ) memory_space(), inputs, global_box, minimum_halo_size, MPI_COMM_WORLD ); // Make a particle list. - Cabana::ParticleTraits + Cabana::ParticleTraits, ParticleVelocity, + Picasso::Field::LogicalPosition<3>, + Picasso::Field::Mass, Picasso::Field::Pressure, + Picasso::Field::Volume, Picasso::Field::DetDefGrad> fields; auto particles = Cabana::Grid::createParticleList( "test_particles", fields ); @@ -304,9 +307,9 @@ void DamBreak( std::string filename ) using particle_type = Picasso::FieldLocation::Particle; using mass_type = Picasso::FieldLayout; using velocity_type = - Picasso::FieldLayout; + Picasso::FieldLayout>; using old_u_type = - Picasso::FieldLayout; + Picasso::FieldLayout>; using delta_u_s_type = Picasso::FieldLayout; using delta_u_g_type = Picasso::FieldLayout; @@ -351,9 +354,9 @@ void DamBreak( std::string filename ) while ( time < final_time ) { // Particle interpolation (Picasso built-in). - Picasso::Particle2Grid - p2g_func{ dt }; + auto p2g_func = + Picasso::createParticle2Grid( dt ); p2g.apply( "p2g_u", particle_type{}, exec_space{}, *fm, particles, p2g_func ); @@ -379,13 +382,15 @@ void DamBreak( std::string filename ) update_u_func ); // Grid interpolation (Picasso built-in). - Picasso::Grid2ParticleVelocity - g2p_func{ dt, beta }; + auto g2p_func = + Picasso::createGrid2ParticleVelocity( 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::Position() ); + particles.redistribute( *local_grid, + Picasso::Field::LogicalPosition<3>() ); // Write particle fields. Cabana::Experimental::HDF5ParticleOutput::HDF5Config h5_config; @@ -394,7 +399,7 @@ void DamBreak( std::string filename ) Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( h5_config, "particles", global_grid.comm(), steps / write_frequency, time, particles.size(), - particles.slice( Picasso::Field::Position() ), + particles.slice( Picasso::Field::LogicalPosition<3>() ), particles.slice( Picasso::Field::Pressure() ), particles.slice( ParticleVelocity() ), particles.slice( Picasso::Field::Mass() ), @@ -436,9 +441,9 @@ int main( int argc, char* argv[] ) std::string filename = argv[1]; // Problem can run with any interpolation scheme. - // DamBreak( filename); - // DamBreak( filename ); - DamBreak( filename ); + // DamBreak>( filename); + // DamBreak>( filename ); + DamBreak>( filename ); Kokkos::finalize(); MPI_Finalize(); diff --git a/src/Picasso_Conservation.hpp b/src/Picasso_Conservation.hpp index 476dbf9c..6191ac72 100644 --- a/src/Picasso_Conservation.hpp +++ b/src/Picasso_Conservation.hpp @@ -62,7 +62,7 @@ void gridConservation( MPI_Comm comm, ExecutionSpace exec_space, auto m_i = fm.view( Picasso::FieldLocation::Node(), Picasso::Field::Mass() ); auto v_i = - fm.view( Picasso::FieldLocation::Node(), Picasso::Field::Velocity() ); + fm.view( Picasso::FieldLocation::Node(), Picasso::Field::Velocity<3>() ); double mass_integral = 0; double ke_integral = 0; diff --git a/src/Picasso_FieldTypes.hpp b/src/Picasso_FieldTypes.hpp index 55400d4b..333afad2 100644 --- a/src/Picasso_FieldTypes.hpp +++ b/src/Picasso_FieldTypes.hpp @@ -981,17 +981,20 @@ struct Volume : Field::Scalar static std::string label() { return "volume"; } }; -struct Velocity : Field::Vector +template +struct Velocity : Field::Vector { static std::string label() { return "velocity"; } }; -struct OldVelocity : Field::Vector +template +struct OldVelocity : Field::Vector { static std::string label() { return "old_velocity"; } }; -struct Stress : Field::Matrix +template +struct Stress : Field::Matrix { static std::string label() { return "stress"; } }; @@ -1010,7 +1013,8 @@ namespace PolyPIC namespace Field { -struct Velocity : Picasso::Field::Matrix +template +struct Velocity : Picasso::Field::Matrix { static std::string label() { return "velocity"; } }; @@ -1023,7 +1027,8 @@ namespace APIC namespace Field { -struct Velocity : Picasso::Field::Matrix +template +struct Velocity : Picasso::Field::Matrix { static std::string label() { return "velocity"; } }; diff --git a/src/Picasso_GridKernels.hpp b/src/Picasso_GridKernels.hpp index 05921d3e..7f5e8a46 100644 --- a/src/Picasso_GridKernels.hpp +++ b/src/Picasso_GridKernels.hpp @@ -24,8 +24,8 @@ namespace Picasso //---------------------------------------------------------------------------// // Compute nodal velocity from mass-weighted momentum. //---------------------------------------------------------------------------// -template +template , + class OldVelType = Field::OldVelocity<3>> struct ComputeGridVelocity { template () ); auto m_p = Picasso::get( particle, MassType() ); auto x_p = Picasso::get( particle, PositionType() ); @@ -184,7 +184,7 @@ struct Particle2Grid, class MassType = Picasso::Field::Mass, class PositionType = Picasso::Field::LogicalPosition<3>> auto createParticle2Grid( const double dt ) @@ -234,7 +234,7 @@ struct Grid2ParticleVelocity() ); auto x_p = Picasso::get( particle, PositionType() ); // Get the gather dependencies. @@ -296,7 +296,7 @@ struct Grid2ParticleVelocity() ); auto x_p = Picasso::get( particle, PositionType() ); // Get the gather dependencies. @@ -362,7 +362,7 @@ struct Grid2ParticleVelocity{} ); auto x_p = Picasso::get( particle, PositionType{} ); // Get the gather dependencies. @@ -371,7 +371,7 @@ struct Grid2ParticleVelocity() ); // Get the local dependencies for getting physcial location of node auto x_i = @@ -420,9 +420,9 @@ struct Grid2ParticleVelocity, - class GridVelocity = Picasso::Field::Velocity, - class GridPosition = Picasso::Field::LogicalPosition<3>> + 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 Date: Mon, 21 Oct 2024 13:31:13 -0400 Subject: [PATCH 20/22] fixup: format --- src/Picasso_Conservation.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Picasso_Conservation.hpp b/src/Picasso_Conservation.hpp index 6191ac72..93ff152e 100644 --- a/src/Picasso_Conservation.hpp +++ b/src/Picasso_Conservation.hpp @@ -61,8 +61,8 @@ void gridConservation( MPI_Comm comm, ExecutionSpace exec_space, 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 v_i = fm.view( Picasso::FieldLocation::Node(), + Picasso::Field::Velocity<3>() ); double mass_integral = 0; double ke_integral = 0; From 271155db3184e23b73a94dadf13aa2b1645feca7 Mon Sep 17 00:00:00 2001 From: "Chong, Kwitae" Date: Mon, 21 Oct 2024 21:51:20 -0400 Subject: [PATCH 21/22] fix wrong folder for clang-format --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 37c99eac..091f9bda 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -113,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 example/*.[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}) From a05937b39f8d3be8430dd1fdf5d49d747180f769 Mon Sep 17 00:00:00 2001 From: "Chong, Kwitae" Date: Mon, 21 Oct 2024 21:51:47 -0400 Subject: [PATCH 22/22] add potential energy --- examples/dam_break.cpp | 23 ++++++++++++++--------- src/Picasso_Conservation.hpp | 32 ++++++++++++++++++++++++++++---- 2 files changed, 42 insertions(+), 13 deletions(-) diff --git a/examples/dam_break.cpp b/examples/dam_break.cpp index cb967755..a414c52c 100644 --- a/examples/dam_break.cpp +++ b/examples/dam_break.cpp @@ -143,7 +143,7 @@ struct ComputeGridVelocityChange { Picasso::Vec3 gravity; - ComputeGridVelocityChange( std::array g ) + ComputeGridVelocityChange( Kokkos::Array g ) { for ( int d = 0; d < 3; ++d ) gravity( d ) = g[d]; @@ -290,7 +290,7 @@ void DamBreak( std::string filename ) // Properties double gamma = inputs["gamma"]; double bulk_modulus = inputs["bulk_modulus"]; - auto gravity = inputs["gravity"]; + auto gravity = copy( inputs["gravity"] ); // Time integragor inputs. double dt = inputs["dt"]; @@ -408,19 +408,24 @@ void DamBreak( std::string filename ) // Calculate conservation sums double mass_particles = 0.0; double ke_particles = 0.0; - Picasso::particleConservation( global_grid.comm(), exec_space(), - particles, ParticleVelocity(), - mass_particles, ke_particles ); + 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; - Picasso::gridConservation( global_grid.comm(), exec_space(), mesh, - *fm, mass_grid, ke_grid ); + 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 Kinetic Energy: " << ke_particles - << " / " << ke_grid << "\n\n"; + << "Particle/Grid Total Energy: " + << ke_particles + pe_particles << " / " + << ke_grid + pe_grid << "\n\n"; } time += dt; steps++; diff --git a/src/Picasso_Conservation.hpp b/src/Picasso_Conservation.hpp index 93ff152e..11c1bf6e 100644 --- a/src/Picasso_Conservation.hpp +++ b/src/Picasso_Conservation.hpp @@ -23,31 +23,40 @@ namespace Picasso template void particleConservation( MPI_Comm comm, ExecutionSpace, const ParticleList particles, ParticleVelocity, - double& global_mass, double& global_ke ) + double& global_mass, double& global_ke, + double& global_pe, const int pe_dir, + const double zero_ground, + const Kokkos::Array 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( 0, num_p ), - KOKKOS_LAMBDA( const int pn, double& mass_cont, double& ke_cont ) { + 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 ); + 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) @@ -55,7 +64,9 @@ template void gridConservation( MPI_Comm comm, ExecutionSpace exec_space, const std::shared_ptr& mesh, Picasso::FieldManager& fm, double& global_mass, - double& global_ke ) + double& global_ke, double& global_pe, const int pe_dir, + const double zero_ground, + const Kokkos::Array gravity ) { auto own_index_space = mesh->localGrid()->indexSpace( Cabana::Grid::Own(), Cabana::Grid::Node(), Cabana::Grid::Local() ); @@ -63,8 +74,11 @@ void gridConservation( MPI_Comm comm, ExecutionSpace exec_space, 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 @@ -85,8 +99,18 @@ void gridConservation( MPI_Comm comm, ExecutionSpace exec_space, }, 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