diff --git a/examples/thermomechanics/CMakeLists.txt b/examples/thermomechanics/CMakeLists.txt index b2633bfa..19c5a103 100644 --- a/examples/thermomechanics/CMakeLists.txt +++ b/examples/thermomechanics/CMakeLists.txt @@ -10,4 +10,11 @@ target_link_libraries(ThermalDeformationHeatTransfer LINK_PUBLIC CabanaPD) add_executable(ThermalDeformationHeatTransferPrenotched thermal_deformation_heat_transfer_prenotched.cpp) target_link_libraries(ThermalDeformationHeatTransferPrenotched LINK_PUBLIC CabanaPD) -install(TARGETS ThermalDeformation ThermalCrack ThermalDeformationHeatTransfer ThermalDeformationHeatTransferPrenotched DESTINATION ${CMAKE_INSTALL_BINDIR}) +add_executable(HIPCylinder hip_cylinder.cpp) +target_link_libraries(HIPCylinder LINK_PUBLIC CabanaPD) + +add_executable(HIPRveCube hip_rve_cube.cpp) +target_link_libraries(HIPRveCube LINK_PUBLIC CabanaPD) + + +install(TARGETS ThermalDeformation ThermalCrack ThermalDeformationHeatTransfer ThermalDeformationHeatTransferPrenotched HIPCylinder HIPRveCube DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/examples/thermomechanics/hip_cylinder.cpp b/examples/thermomechanics/hip_cylinder.cpp new file mode 100644 index 00000000..644772d2 --- /dev/null +++ b/examples/thermomechanics/hip_cylinder.cpp @@ -0,0 +1,305 @@ +/**************************************************************************** + * Copyright (c) 2022-2023 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of CabanaPD. CabanaPD 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 * + ****************************************************************************/ + +#include +#include + +#include "mpi.h" + +#include + +#include + +// Simulate a cylinder under hot isostatic pressing (HIP). +void HIPCylinderExample( const std::string filename ) +{ + // ==================================================== + // Choose Kokkos spaces + // ==================================================== + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = typename exec_space::memory_space; + + // ==================================================== + // Read inputs + // ==================================================== + CabanaPD::Inputs inputs( filename ); + + // ==================================================== + // Material parameters + // ==================================================== + double rho0 = inputs["density"]; + double E = inputs["elastic_modulus"]; + double nu = 0.25; // Use bond-based model + double K = E / ( 3 * ( 1 - 2 * nu ) ); + double G0 = inputs["fracture_energy"]; + double sigma_y = inputs["yield_stress"]; + double delta = inputs["horizon"]; + delta += 1e-10; + double alpha = inputs["thermal_expansion_coeff"]; + double temp0 = inputs["reference_temperature"]; + + // ==================================================== + // Discretization + // ==================================================== + std::array low_corner = inputs["low_corner"]; + std::array high_corner = inputs["high_corner"]; + std::array num_cells = inputs["num_cells"]; + int m = std::floor( delta / + ( ( high_corner[0] - low_corner[0] ) / num_cells[0] ) ); + int halo_width = m + 1; // Just to be safe. + + // ==================================================== + // Force model type + // ==================================================== + using model_type = CabanaPD::PMB; + using mechanics_type = CabanaPD::ElasticPerfectlyPlastic; + using thermal_type = CabanaPD::TemperatureDependent; + using density_type = CabanaPD::DynamicDensity; + + // ==================================================== + // Custom particle generation and initialization + // ==================================================== + double x_center = 0.5 * ( low_corner[0] + high_corner[0] ); + double y_center = 0.5 * ( low_corner[1] + high_corner[1] ); + double z_center = 0.5 * ( low_corner[2] + high_corner[2] ); + double Rout = inputs["cylinder_outer_radius"]; + double Rin = inputs["cylinder_inner_radius"]; + double H = inputs["cylinder_height"]; + + // Do not create particles outside given cylindrical region + auto init_op = KOKKOS_LAMBDA( const int, const double x[3] ) + { + double rsq = ( x[0] - x_center ) * ( x[0] - x_center ) + + ( x[1] - y_center ) * ( x[1] - y_center ); + if ( rsq < Rin * Rin || rsq > Rout * Rout || + x[2] > z_center + 0.5 * H || x[2] < z_center - 0.5 * H ) + return false; + return true; + }; + + CabanaPD::Particles particles( + memory_space{}, model_type{}, thermal_type{}, density_type{}, + low_corner, high_corner, num_cells, halo_width, init_op, exec_space{} ); + + // Impose separate density values for powder and container particles. + double W = inputs["wall_thickness"]; + double D0 = inputs["powder_initial_relative_density"]; + auto rho = particles.sliceDensity(); + auto temp = particles.sliceTemperature(); + auto x = particles.sliceReferencePosition(); + auto nofail = particles.sliceNoFail(); + + using pool_type = Kokkos::Random_XorShift64_Pool; + using random_type = Kokkos::Random_XorShift64; + pool_type pool; + int seed = 456854; + pool.init( seed, particles.numLocal() ); + + // Use time to seed random number generator + // std::srand( std::time( nullptr ) ); + double rho_perturb_factor = 0.02; + + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) + + ( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center ); + if ( rsq > ( Rin + W ) * ( Rin + W ) && + rsq < ( Rout - W ) * ( Rout - W ) && + x( pid, 2 ) < z_center + 0.5 * H - W && + x( pid, 2 ) > z_center - 0.5 * H + W ) + { + // Powder density + rho( pid ) = D0 * rho0; + // Perturb powder density + auto gen = pool.get_state(); + auto rand = + Kokkos::rand::draw( gen, 0.0, 1.0 ); + double factor = ( 1 + ( 2.0 * rand - 1.0 ) * rho_perturb_factor ); + + // Check perturbed density does not exceed maximum density + if ( D0 * factor < 1 ) + { + rho( pid ) *= factor; + } + else + { + rho( pid ) = rho0; + } + + // Free the state after drawing + pool.free_state( gen ); + } + else + { + // Container density + rho( pid ) = rho0; + }; + + // Temperature + temp( pid ) = temp0; + + // No fail + nofail( pid ) = 1; + }; + particles.updateParticles( exec_space{}, init_functor ); + + // ==================================================== + // Boundary conditions planes + // ==================================================== + CabanaPD::RegionBoundary plane( + low_corner[0], high_corner[0], low_corner[1], high_corner[1], + low_corner[2], high_corner[2] ); + + // ==================================================== + // Force model + // ==================================================== + rho = particles.sliceDensity(); + auto rho_current = particles.sliceCurrentDensity(); + temp = particles.sliceTemperature(); + CabanaPD::ForceDensityModel force_model( + model_type{}, mechanics_type{}, rho, rho_current, delta, K, G0, sigma_y, + rho0, temp, alpha, temp0 ); + + // ==================================================== + // Create solver + // ==================================================== + CabanaPD::Solver solver( inputs, particles, force_model ); + /* + using contact_type = CabanaPD::NormalRepulsionModel; + auto dx = particles.dx; + + // Use contact radius and extension relative to particle spacing. + double r_c = inputs["contact_horizon_factor"]; + double r_extend = inputs["contact_horizon_extend_factor"]; + // NOTE: dx/2 is when particles first touch. + r_c *= dx[0] / 2.0; + r_extend *= dx[0]; + + contact_type contact_model( delta, r_c, r_extend, K ); + + CabanaPD::Solver solver( inputs, particles, force_model, + contact_model ); + */ + + // ==================================================== + // Impose field + // ==================================================== + // Create BC last to ensure ghost particles are included. + double Pmax = inputs["maximum_pressure"]; + double trampP = inputs["ramp_pressure_time"]; + double Tmax = inputs["maximum_temperature"]; + double dx = solver.particles.dx[0]; + double dz = solver.particles.dx[2]; + double b0max = Pmax / dx; + x = solver.particles.sliceReferencePosition(); + auto f = solver.particles.sliceForce(); + auto u = solver.particles.sliceDisplacement(); + temp = solver.particles.sliceTemperature(); + + // This is purposely delayed until after solver init so that ghosted + // particles are correctly taken into account for lambda capture here. + auto force_temp_func = KOKKOS_LAMBDA( const int pid, const double t ) + { + double b0; + + // Pressure ramping + if ( t < trampP ) + { + b0 = b0max * t / trampP; + } + else + { + b0 = b0max; + } + + // ----------------------- + // Isostatic pressure BC + // ----------------------- + double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) + + ( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center ); + double theta = + Kokkos::atan2( x( pid, 1 ) - y_center, x( pid, 0 ) - x_center ); + + // BC on outer boundary + if ( rsq > ( Rout - W ) * ( Rout - W ) ) + { + f( pid, 0 ) += -b0 * Kokkos::cos( theta ); + f( pid, 1 ) += -b0 * Kokkos::sin( theta ); + } + // BC on inner boundary + else if ( rsq < ( Rin + W ) * ( Rin + W ) ) + { + f( pid, 0 ) += b0 * Kokkos::cos( theta ); + f( pid, 1 ) += b0 * Kokkos::sin( theta ); + }; + + // BC on top boundary + if ( x( pid, 2 ) > z_center + 0.5 * H - W ) + { + f( pid, 2 ) += -b0; + } + // BC on bottom boundary + else if ( x( pid, 2 ) < z_center - 0.5 * H + W ) + { + f( pid, 2 ) += b0; + }; + + // ----------------------- + // Temperature BC + // ----------------------- + temp( pid ) = Tmax; + + // ----------------------- + // Constrain displacements + // ----------------------- + // Fix x- and y-displacement on top surface: only enable motion in + // z-direction + if ( x( pid, 2 ) > z_center + 0.5 * H - W ) + { + u( pid, 0 ) = 0.0; + u( pid, 1 ) = 0.0; + } + + // Fix x- and y-displacement on bottom surface: only enable motion in + // z-direction + if ( x( pid, 2 ) < z_center - 0.5 * H + W ) + { + u( pid, 0 ) = 0.0; + u( pid, 1 ) = 0.0; + } + + // Fix z-displacement on mid surface + if ( x( pid, 2 ) > z_center - dz && x( pid, 2 ) < z_center + dz ) + { + u( pid, 2 ) = 0.0; + } + }; + CabanaPD::BodyTerm body_term( force_temp_func, solver.particles.size(), + true ); + + // ==================================================== + // Simulation run + // ==================================================== + solver.init( body_term ); + solver.run( body_term ); +} + +// Initialize MPI+Kokkos. +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + HIPCylinderExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} diff --git a/examples/thermomechanics/hip_rve_cube.cpp b/examples/thermomechanics/hip_rve_cube.cpp new file mode 100644 index 00000000..2cfb4cbe --- /dev/null +++ b/examples/thermomechanics/hip_rve_cube.cpp @@ -0,0 +1,376 @@ +/**************************************************************************** + * Copyright (c) 2022-2023 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of CabanaPD. CabanaPD 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 * + ****************************************************************************/ + +#include +#include + +#include "mpi.h" + +#include + +#include + +// Simulate a cylinder under hot isostatic pressing (HIP). +void HIPCylinderExample( const std::string filename ) +{ + // ==================================================== + // Choose Kokkos spaces + // ==================================================== + using exec_space = Kokkos::DefaultExecutionSpace; + using memory_space = typename exec_space::memory_space; + + // ==================================================== + // Read inputs + // ==================================================== + CabanaPD::Inputs inputs( filename ); + + // ==================================================== + // Material parameters + // ==================================================== + double rho0 = inputs["density"]; + double E = inputs["elastic_modulus"]; + double nu = 0.25; // Use bond-based model + double K = E / ( 3 * ( 1 - 2 * nu ) ); + double G0 = inputs["fracture_energy"]; + double sigma_y = inputs["yield_stress"]; + double delta = inputs["horizon"]; + delta += 1e-10; + double alpha = inputs["thermal_expansion_coeff"]; + double temp0 = inputs["reference_temperature"]; + + // ==================================================== + // Discretization + // ==================================================== + std::array low_corner = inputs["low_corner"]; + std::array high_corner = inputs["high_corner"]; + std::array num_cells = inputs["num_cells"]; + int m = std::floor( delta / + ( ( high_corner[0] - low_corner[0] ) / num_cells[0] ) ); + int halo_width = m + 1; // Just to be safe. + + // ==================================================== + // Force model type + // ==================================================== + using model_type = CabanaPD::PMB; + using mechanics_type = CabanaPD::ElasticPerfectlyPlastic; + using thermal_type = CabanaPD::TemperatureDependent; + using density_type = CabanaPD::DynamicDensity; + + // ==================================================== + // Custom particle generation and initialization + // ==================================================== + double x_center = 0.5 * ( low_corner[0] + high_corner[0] ); + double y_center = 0.5 * ( low_corner[1] + high_corner[1] ); + double z_center = 0.5 * ( low_corner[2] + high_corner[2] ); + // double Rout = inputs["cylinder_outer_radius"]; + // double Rin = inputs["cylinder_inner_radius"]; + // double H = inputs["cylinder_height"]; + double L = inputs["cube_side"]; + + // // Do not create particles outside given cylindrical region + // Do not create particles outside given cube region + auto init_op = KOKKOS_LAMBDA( const int, const double x[3] ) + { + /* + double rsq = ( x[0] - x_center ) * ( x[0] - x_center ) + + ( x[1] - y_center ) * ( x[1] - y_center ); + if ( rsq < Rin * Rin || rsq > Rout * Rout || + x[2] > z_center + 0.5 * H || x[2] < z_center - 0.5 * H ) + return false; + */ + + if ( Kokkos::abs( x[0] - x_center ) > 0.5 * L || + Kokkos::abs( x[1] - y_center ) > 0.5 * L || + Kokkos::abs( x[2] - z_center ) > 0.5 * L ) + return false; + return true; + }; + + CabanaPD::Particles particles( + memory_space{}, model_type{}, thermal_type{}, density_type{}, + low_corner, high_corner, num_cells, halo_width, init_op, exec_space{} ); + + // Impose separate density values for powder and container particles. + double W = inputs["wall_thickness"]; + double D0 = inputs["powder_initial_relative_density"]; + auto rho = particles.sliceDensity(); + auto temp = particles.sliceTemperature(); + auto x = particles.sliceReferencePosition(); + auto nofail = particles.sliceNoFail(); + + using pool_type = Kokkos::Random_XorShift64_Pool; + using random_type = Kokkos::Random_XorShift64; + pool_type pool; + int seed = 456854; + pool.init( seed, particles.numLocal() ); + + // Use time to seed random number generator + // std::srand( std::time( nullptr ) ); + double rho_perturb_factor = 0.02; + + auto init_functor = KOKKOS_LAMBDA( const int pid ) + { + // double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) + // + + // ( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center ); + // if ( rsq > ( Rin + W ) * ( Rin + W ) && + // rsq < ( Rout - W ) * ( Rout - W ) && + // x( pid, 2 ) < z_center + 0.5 * H - W && + // x( pid, 2 ) > z_center - 0.5 * H + W ) + + if ( Kokkos::abs( x( pid, 0 ) - x_center ) < 0.5 * L - W && + Kokkos::abs( x( pid, 1 ) - y_center ) < 0.5 * L - W && + Kokkos::abs( x( pid, 2 ) - z_center ) < 0.5 * L - W ) + { + // Powder density + rho( pid ) = D0 * rho0; + // Perturb powder density + auto gen = pool.get_state(); + auto rand = + Kokkos::rand::draw( gen, 0.0, 1.0 ); + double factor = ( 1 + ( 2.0 * rand - 1.0 ) * rho_perturb_factor ); + + // Check perturbed density does not exceed maximum density + if ( D0 * factor < 1 ) + { + rho( pid ) *= factor; + } + else + { + rho( pid ) = rho0; + } + + // Free the state after drawing + pool.free_state( gen ); + } + else + { + // Container density + rho( pid ) = rho0; + }; + + // Temperature + temp( pid ) = temp0; + + // No fail + nofail( pid ) = 1; + }; + particles.updateParticles( exec_space{}, init_functor ); + + // ==================================================== + // Boundary conditions planes + // ==================================================== + CabanaPD::RegionBoundary plane( + low_corner[0], high_corner[0], low_corner[1], high_corner[1], + low_corner[2], high_corner[2] ); + + // ==================================================== + // Force model + // ==================================================== + rho = particles.sliceDensity(); + auto rho_current = particles.sliceCurrentDensity(); + temp = particles.sliceTemperature(); + CabanaPD::ForceDensityModel force_model( + model_type{}, mechanics_type{}, rho, rho_current, delta, K, G0, sigma_y, + rho0, temp, alpha, temp0 ); + + // ==================================================== + // Create solver + // ==================================================== + CabanaPD::Solver solver( inputs, particles, force_model ); + /* + using contact_type = CabanaPD::NormalRepulsionModel; + auto dx = particles.dx; + + // Use contact radius and extension relative to particle spacing. + double r_c = inputs["contact_horizon_factor"]; + double r_extend = inputs["contact_horizon_extend_factor"]; + // NOTE: dx/2 is when particles first touch. + r_c *= dx[0] / 2.0; + r_extend *= dx[0]; + + contact_type contact_model( delta, r_c, r_extend, K ); + + CabanaPD::Solver solver( inputs, particles, force_model, + contact_model ); + */ + + // ==================================================== + // Impose field + // ==================================================== + // Create BC last to ensure ghost particles are included. + double Pmax = inputs["maximum_pressure"]; + double trampP = inputs["ramp_pressure_time"]; + double Tmax = inputs["maximum_temperature"]; + double dx = solver.particles.dx[0]; + double dz = solver.particles.dx[2]; + // double b0max = Pmax / dx; + double b0max = Pmax / W; + x = solver.particles.sliceReferencePosition(); + auto f = solver.particles.sliceForce(); + auto u = solver.particles.sliceDisplacement(); + temp = solver.particles.sliceTemperature(); + + // This is purposely delayed until after solver init so that ghosted + // particles are correctly taken into account for lambda capture here. + auto force_temp_func = KOKKOS_LAMBDA( const int pid, const double t ) + { + double b0; + + // Pressure ramping + if ( t < trampP ) + { + b0 = b0max * t / trampP; + } + else + { + b0 = b0max; + } + + // ----------------------- + // Isostatic pressure BC + // ----------------------- + // double rsq = ( x( pid, 0 ) - x_center ) * ( x( pid, 0 ) - x_center ) + // + + // ( x( pid, 1 ) - y_center ) * ( x( pid, 1 ) - y_center ); + // double theta = + // Kokkos::atan2( x( pid, 1 ) - y_center, x( pid, 0 ) - x_center ); + + // if ( kokkos::abs(x[0] - x_center) < 0.5*L && kokkos::abs(x[1] - + // y_center) < 0.5*L && kokkos::abs(x[2] - z_center) < 0.5*L ) + + // Pressure in x-direction + if ( x( pid, 0 ) > x_center + 0.5 * L - W ) + f( pid, 0 ) += -b0; + + if ( x( pid, 0 ) < x_center - 0.5 * L + W ) + f( pid, 0 ) += b0; + + // Pressure in y-direction + if ( x( pid, 1 ) > y_center + 0.5 * L - W ) + f( pid, 1 ) += -b0; + + if ( x( pid, 1 ) < y_center - 0.5 * L + W ) + f( pid, 1 ) += b0; + + // Pressure in z-direction + if ( x( pid, 2 ) > z_center + 0.5 * L - W ) + f( pid, 2 ) += -b0; + + if ( x( pid, 2 ) < z_center - 0.5 * L + W ) + f( pid, 2 ) += b0; + + /* + // BC on outer boundary + if ( rsq > ( Rout - W ) * ( Rout - W ) ) + { + f( pid, 0 ) += -b0 * Kokkos::cos( theta ); + f( pid, 1 ) += -b0 * Kokkos::sin( theta ); + } + // BC on inner boundary + else if ( rsq < ( Rin + W ) * ( Rin + W ) ) + { + f( pid, 0 ) += b0 * Kokkos::cos( theta ); + f( pid, 1 ) += b0 * Kokkos::sin( theta ); + }; + + // BC on top boundary + if ( x( pid, 2 ) > z_center + 0.5 * H - W ) + { + f( pid, 2 ) += -b0; + } + // BC on bottom boundary + else if ( x( pid, 2 ) < z_center - 0.5 * H + W ) + { + f( pid, 2 ) += b0; + }; + */ + + // ----------------------- + // Temperature BC + // ----------------------- + temp( pid ) = Tmax; + + // ----------------------- + // Constrain displacements + // ----------------------- + + // Only enable motion in x-direction on surfaces with x-normal + if ( x( pid, 0 ) > x_center + 0.5 * L - W || + x( pid, 0 ) < x_center - 0.5 * L + W ) + { + u( pid, 1 ) = 0.0; + u( pid, 2 ) = 0.0; + } + + // Only enable motion in y-direction on surfaces with y-normal + if ( x( pid, 1 ) > y_center + 0.5 * L - W || + x( pid, 1 ) < y_center - 0.5 * L + W ) + { + u( pid, 0 ) = 0.0; + u( pid, 2 ) = 0.0; + } + + // Only enable motion in z-direction on surfaces with z-normal + if ( x( pid, 2 ) > z_center + 0.5 * L - W || + x( pid, 2 ) < z_center - 0.5 * L + W ) + { + u( pid, 0 ) = 0.0; + u( pid, 1 ) = 0.0; + } + + /* + // Fix x- and y-displacement on top surface: only enable motion in + // z-direction + // if ( x( pid, 2 ) > z_center + 0.5 * H - W ) + if ( x( pid, 2 ) > z_center + 0.5 * L - W ) + { + u( pid, 0 ) = 0.0; + u( pid, 1 ) = 0.0; + } + + // Fix x- and y-displacement on bottom surface: only enable motion in + // z-direction + // if ( x( pid, 2 ) < z_center - 0.5 * H + W ) + if ( x( pid, 2 ) < z_center - 0.5 * L + W ) + { + u( pid, 0 ) = 0.0; + u( pid, 1 ) = 0.0; + } + + // Fix z-displacement on mid surface + if ( x( pid, 2 ) > z_center - dz && x( pid, 2 ) < z_center + dz ) + { + u( pid, 2 ) = 0.0; + } + */ + }; + CabanaPD::BodyTerm body_term( force_temp_func, solver.particles.size(), + true ); + + // ==================================================== + // Simulation run + // ==================================================== + solver.init( body_term ); + solver.run( body_term ); +} + +// Initialize MPI+Kokkos. +int main( int argc, char* argv[] ) +{ + MPI_Init( &argc, &argv ); + Kokkos::initialize( argc, argv ); + + HIPCylinderExample( argv[1] ); + + Kokkos::finalize(); + MPI_Finalize(); +} diff --git a/examples/thermomechanics/inputs/hip_cylinder.json b/examples/thermomechanics/inputs/hip_cylinder.json new file mode 100644 index 00000000..9844b009 --- /dev/null +++ b/examples/thermomechanics/inputs/hip_cylinder.json @@ -0,0 +1,27 @@ +{ + "num_cells" : {"value": [100, 100, 140]}, + "system_size" : {"value": [0.25, 0.25, 0.35], "unit": "m"}, + "density" : {"value": 7950, "unit": "kg/m^3"}, + "powder_initial_relative_density" : {"value": 0.85}, + "elastic_modulus" : {"value": 195.6e+9, "unit": "Pa"}, + "fracture_energy" : {"value": 1550e+6, "unit": "J/m^2"}, + "horizon" : {"value": 0.0075, "unit": "m"}, + "yield_stress_ref" : {"value": 170e6, "unit": "Pa"}, + "yield_stress" : {"value": 170e8, "unit": "Pa"}, + "thermal_expansion_coeff" : {"value": 14.56E-6, "unit": "oC^{-1}"}, + "thermal_conductivity" : {"value": 14.12, "unit": "W/(m.K)"}, + "specific_heat_capacity" : {"value": 492, "unit": "J/(kg.K)"}, + "maximum_pressure" : {"value": 110e6, "unit": "Pa"}, + "reference_temperature" : {"value": 1125, "unit": "oC"}, + "maximum_temperature" : {"value": 1125, "unit": "oC"}, + "cylinder_outer_radius" : {"value": 0.1187, "unit": "m"}, + "cylinder_inner_radius" : {"value": 0.0679, "unit": "m"}, + "cylinder_height" : {"value": 0.3048, "unit": "m"}, + "wall_thickness" : {"value": 0.0045, "unit": "m"}, + "ramp_pressure_time" : {"value": 2.5e-1, "unit": "s"}, + "final_time" : {"value": 2.5e-2, "unit": "s"}, + "timestep" : {"value": 3.0e-07, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.70}, + "output_frequency" : {"value": 100}, + "output_reference" : {"value": false} +} diff --git a/examples/thermomechanics/inputs/hip_cylinder_temp.json b/examples/thermomechanics/inputs/hip_cylinder_temp.json new file mode 100644 index 00000000..f06853ae --- /dev/null +++ b/examples/thermomechanics/inputs/hip_cylinder_temp.json @@ -0,0 +1,25 @@ +{ + "num_cells" : {"value": [100, 100, 140]}, + "system_size" : {"value": [0.25, 0.25, 0.35], "unit": "m"}, + "density" : {"value": 7950, "unit": "kg/m^3"}, + "elastic_modulus" : {"value": 195.6e+8, "unit": "Pa"}, + "fracture_energy" : {"value": 1550e+6, "unit": "J/m^2"}, + "horizon" : {"value": 0.0075, "unit": "m"}, + "yield_stress_ref" : {"value": 170e7, "unit": "Pa"}, + "yield_stress" : {"value": 500e6, "unit": "Pa"}, + "thermal_expansion_coeff" : {"value": 14.56E-6, "unit": "oC^{-1}"}, + "thermal_conductivity" : {"value": 14.12, "unit": "W/(m.K)"}, + "specific_heat_capacity" : {"value": 492, "unit": "J/(kg.K)"}, + "maximum_pressure" : {"value": 110e6, "unit": "Pa"}, + "reference_temperature" : {"value": 1125, "unit": "oC"}, + "maximum_temperature" : {"value": 1125, "unit": "oC"}, + "cylinder_outer_radius" : {"value": 0.1187, "unit": "m"}, + "cylinder_inner_radius" : {"value": 0.0679, "unit": "m"}, + "cylinder_height" : {"value": 0.3048, "unit": "m"}, + "wall_thickness" : {"value": 0.0067, "unit": "m"}, + "final_time" : {"value": 1e-3, "unit": "s"}, + "timestep" : {"value": 3.0e-07, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.70}, + "output_frequency" : {"value": 50}, + "output_reference" : {"value": false} +} diff --git a/examples/thermomechanics/inputs/hip_rve_cube.json b/examples/thermomechanics/inputs/hip_rve_cube.json new file mode 100644 index 00000000..b20e409c --- /dev/null +++ b/examples/thermomechanics/inputs/hip_rve_cube.json @@ -0,0 +1,24 @@ +{ + "num_cells" : {"value": [80, 80, 80]}, + "system_size" : {"value": [0.1, 0.1, 0.1], "unit": "m"}, + "density" : {"value": 7950, "unit": "kg/m^3"}, + "powder_initial_relative_density" : {"value": 0.85}, + "elastic_modulus" : {"value": 195.6e+9, "unit": "Pa"}, + "fracture_energy" : {"value": 1550e+6, "unit": "J/m^2"}, + "horizon" : {"value": 0.00375, "unit": "m"}, + "yield_stress" : {"value": 170e8, "unit": "Pa"}, + "thermal_expansion_coeff" : {"value": 14.56E-6, "unit": "oC^{-1}"}, + "thermal_conductivity" : {"value": 14.12, "unit": "W/(m.K)"}, + "specific_heat_capacity" : {"value": 492, "unit": "J/(kg.K)"}, + "maximum_pressure" : {"value": 700e7, "unit": "Pa"}, + "reference_temperature" : {"value": 1125, "unit": "oC"}, + "maximum_temperature" : {"value": 1125, "unit": "oC"}, + "cube_side" : {"value": 0.1, "unit": "m"}, + "wall_thickness" : {"value": 0.0045, "unit": "m"}, + "ramp_pressure_time" : {"value": 6.0e-3, "unit": "s"}, + "final_time" : {"value": 0.0027, "unit": "s"}, + "timestep" : {"value": 1.5e-07, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.70}, + "output_frequency" : {"value": 250}, + "output_reference" : {"value": false} +} diff --git a/examples/thermomechanics/inputs/hip_rve_cube_ref.json b/examples/thermomechanics/inputs/hip_rve_cube_ref.json new file mode 100644 index 00000000..4a3b4e5e --- /dev/null +++ b/examples/thermomechanics/inputs/hip_rve_cube_ref.json @@ -0,0 +1,35 @@ +{ + "num_cells_OLD" : {"value": [100, 100, 140]}, + "num_cells" : {"value": [20, 20, 20]}, + "system_size_OLD" : {"value": [0.25, 0.25, 0.35], "unit": "m"}, + "system_size" : {"value": [0.1, 0.1, 0.1], "unit": "m"}, + "density" : {"value": 7950, "unit": "kg/m^3"}, + "powder_initial_relative_density" : {"value": 0.85}, + "elastic_modulus" : {"value": 195.6e+9, "unit": "Pa"}, + "fracture_energy" : {"value": 1550e+6, "unit": "J/m^2"}, + "horizon_OLD" : {"value": 0.0075, "unit": "m"}, + "horizon" : {"value": 0.015, "unit": "m"}, + "yield_stress_ref" : {"value": 170e6, "unit": "Pa"}, + "yield_stress" : {"value": 170e8, "unit": "Pa"}, + "thermal_expansion_coeff" : {"value": 14.56E-6, "unit": "oC^{-1}"}, + "thermal_conductivity" : {"value": 14.12, "unit": "W/(m.K)"}, + "specific_heat_capacity" : {"value": 492, "unit": "J/(kg.K)"}, + "maximum_pressure_OLD" : {"value": 110e6, "unit": "Pa"}, + "maximum_pressure" : {"value": 550e7, "unit": "Pa"}, + "reference_temperature" : {"value": 1125, "unit": "oC"}, + "maximum_temperature" : {"value": 1125, "unit": "oC"}, + "cube_side" : {"value": 0.1, "unit": "m"}, + "cylinder_outer_radius_OLD" : {"value": 0.1187, "unit": "m"}, + "cylinder_inner_radius_OLD" : {"value": 0.0679, "unit": "m"}, + "cylinder_heigh_OLD" : {"value": 0.3048, "unit": "m"}, + "wall_thickness" : {"value": 0.0045, "unit": "m"}, + "ramp_pressure_time_OLD" : {"value": 2.5e-1, "unit": "s"}, + "ramp_pressure_time" : {"value": 6.0e-3, "unit": "s"}, + "final_time_OLD" : {"value": 2.5e-2, "unit": "s"}, + "final_time" : {"value": 0.0054, "unit": "s"}, + "timestep_OLD" : {"value": 3.0e-07, "unit": "s"}, + "timestep" : {"value": 6.0e-07, "unit": "s"}, + "timestep_safety_factor" : {"value": 0.70}, + "output_frequency" : {"value": 50}, + "output_reference" : {"value": false} +} diff --git a/src/CabanaPD_Force.hpp b/src/CabanaPD_Force.hpp index 28703df6..47589a72 100644 --- a/src/CabanaPD_Force.hpp +++ b/src/CabanaPD_Force.hpp @@ -63,6 +63,7 @@ #include #include +#include #include namespace CabanaPD @@ -137,165 +138,229 @@ getLinearizedDistance( const PosType& x, const PosType& u, const int i, } // Forward declaration. -template +template class Force; template -class Force +class BaseForce { - public: - using neighbor_list_type = - Cabana::VerletList; - protected: - bool _half_neigh; - neighbor_list_type _neigh_list; - Timer _timer; Timer _energy_timer; Timer _stress_timer; public: - // Primary constructor: use positions and construct neighbors. - template - Force( const bool half_neigh, const double delta, - const ParticleType& particles, const double tol = 1e-14 ) - : _half_neigh( half_neigh ) - , _neigh_list( neighbor_list_type( - particles.sliceReferencePosition(), particles.frozenOffset(), - particles.localOffset(), delta + tol, 1.0, - particles.ghost_mesh_lo, particles.ghost_mesh_hi ) ) + // Default to no-op. + template + void computeWeightedVolume( ParticleType&, const NeighborType ) const { } - - // General constructor (necessary for contact, but could be used by any - // force routine). - template - Force( const bool half_neigh, const double delta, - const PositionType& positions, const std::size_t frozen_offset, - const std::size_t local_offset, const double mesh_min[3], - const double mesh_max[3], const double tol = 1e-14 ) - : _half_neigh( half_neigh ) - , _neigh_list( neighbor_list_type( positions, frozen_offset, - local_offset, delta + tol, 1.0, - mesh_min, mesh_max ) ) + template + void computeDilatation( ParticleType&, const NeighborType ) const { } - // Constructor which stores existing neighbors. - template - Force( const bool half_neigh, const NeighborListType& neighbors ) - : _half_neigh( half_neigh ) - , _neigh_list( neighbors ) - { - } + auto time() { return _timer.time(); }; + auto timeEnergy() { return _energy_timer.time(); }; +}; + +/****************************************************************************** + Dilatation. +******************************************************************************/ +template +class Dilatation; + +template +class Dilatation +{ + protected: + using model_type = ModelType; + model_type _model; - // FIXME: should it be possible to update this list? - template - void update( const ParticleType&, const double, const bool = false ) + Timer _timer; + + public: + // Using the default exec_space. + using exec_space = typename MemorySpace::execution_space; + + Dilatation( const model_type model ) + : _model( model ) { } - unsigned getMaxLocalNeighbors() + template + void computeWeightedVolume( ParticleType& particles, + const NeighborType& neighbor ) { - auto neigh = _neigh_list; - unsigned local_max_neighbors; - auto neigh_max = KOKKOS_LAMBDA( const int, unsigned& max_n ) + _timer.start(); + + auto x = particles.sliceReferencePosition(); + auto u = particles.sliceDisplacement(); + const auto vol = particles.sliceVolume(); + auto m = particles.sliceWeightedVolume(); + Cabana::deep_copy( m, 0.0 ); + auto model = _model; + + auto weighted_volume = KOKKOS_LAMBDA( const int i, const int j ) { - max_n = - Cabana::NeighborList::maxNeighbor( neigh ); + // Get the reference positions and displacements. + double xi, r, s; + getDistance( x, u, i, j, xi, r, s ); + m( i ) += model.weightedVolume( xi, vol( j ) ); }; - using exec_space = typename MemorySpace::execution_space; - Kokkos::RangePolicy policy( 0, 1 ); - Kokkos::parallel_reduce( policy, neigh_max, local_max_neighbors ); - Kokkos::fence(); - return local_max_neighbors; + neighbor.iterate( exec_space{}, weighted_volume, particles, + "CabanaPD::Dilatation::computeWeightedVolume" ); + _timer.stop(); } - void getNeighborStatistics( unsigned& max_neighbors, - unsigned long long& total_neighbors ) + template + void computeDilatation( ParticleType& particles, + const NeighborType neighbor ) { - auto neigh = _neigh_list; - unsigned local_max_neighbors; - unsigned long long local_total_neighbors; - auto neigh_stats = KOKKOS_LAMBDA( const int, unsigned& max_n, - unsigned long long& total_n ) + _timer.start(); + + const auto x = particles.sliceReferencePosition(); + auto u = particles.sliceDisplacement(); + const auto vol = particles.sliceVolume(); + auto m = particles.sliceWeightedVolume(); + auto theta = particles.sliceDilatation(); + auto model = _model; + Cabana::deep_copy( theta, 0.0 ); + + auto dilatation = KOKKOS_LAMBDA( const int i, const int j ) { - max_n = - Cabana::NeighborList::maxNeighbor( neigh ); - total_n = Cabana::NeighborList::totalNeighbor( - neigh ); + // Get the bond distance, displacement, and stretch. + double xi, r, s; + getDistance( x, u, i, j, xi, r, s ); + theta( i ) += model.dilatation( i, s, xi, vol( j ), m( i ) ); }; - using exec_space = typename MemorySpace::execution_space; - Kokkos::RangePolicy policy( 0, 1 ); - Kokkos::parallel_reduce( policy, neigh_stats, local_max_neighbors, - local_total_neighbors ); - Kokkos::fence(); - MPI_Reduce( &local_max_neighbors, &max_neighbors, 1, MPI_UNSIGNED, - MPI_MAX, 0, MPI_COMM_WORLD ); - MPI_Reduce( &local_total_neighbors, &total_neighbors, 1, - MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD ); - } - // Default to no-op for pair models. - template - void computeWeightedVolume( ParticleType&, const ParallelType ) const - { - } - template - void computeDilatation( ParticleType&, const ParallelType ) const - { - } + neighbor.iterate( exec_space{}, dilatation, particles, + "CabanaPD::Dilatation::compute" ); - auto getNeighbors() const { return _neigh_list; } + _timer.stop(); + } auto time() { return _timer.time(); }; - auto timeEnergy() { return _energy_timer.time(); }; - auto timeNeighbor() { return 0.0; }; }; -template -class BaseFracture +template +class Dilatation { protected: - using memory_space = MemorySpace; - using NeighborView = typename Kokkos::View; - NeighborView _mu; + using model_type = ModelType; + model_type _model; + + Timer _timer; public: - BaseFracture( const int local_particles, const int max_neighbors ) + // Using the default exec_space. + using exec_space = typename MemorySpace::execution_space; + + Dilatation( const model_type model ) + : _model( model ) { - // Create View to track broken bonds. - // TODO: this could be optimized to ignore frozen particle bonds. - _mu = NeighborView( - Kokkos::ViewAllocateWithoutInitializing( "broken_bonds" ), - local_particles, max_neighbors ); - Kokkos::deep_copy( _mu, 1 ); } - BaseFracture( NeighborView mu ) - : _mu( mu ) + template + void computeWeightedVolume( ParticleType& particles, + const NeighborType& neighbor ) { + _timer.start(); + + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); + + auto x = particles.sliceReferencePosition(); + auto u = particles.sliceDisplacement(); + const auto vol = particles.sliceVolume(); + auto m = particles.sliceWeightedVolume(); + Cabana::deep_copy( m, 0.0 ); + auto model = _model; + + auto weighted_volume = KOKKOS_LAMBDA( const int i ) + { + std::size_t num_neighbors = + Cabana::NeighborList::numNeighbor( + neigh_list, i ); + for ( std::size_t n = 0; n < num_neighbors; n++ ) + { + std::size_t j = + Cabana::NeighborList::getNeighbor( + neigh_list, i, n ); + + // Get the reference positions and displacements. + double xi, r, s; + getDistance( x, u, i, j, xi, r, s ); + // mu is included to account for bond breaking. + m( i ) += mu( i, n ) * model.weightedVolume( xi, vol( j ) ); + } + }; + + neighbor.iterateLinear( exec_space{}, weighted_volume, particles, + "CabanaPD::Dilatation::computeWeightedVolume" ); + + _timer.stop(); } - template - void prenotch( ExecSpace exec_space, const ParticleType& particles, - PrenotchType& prenotch, NeighborList& neigh_list ) + template + void computeDilatation( ParticleType& particles, + const NeighborType& neighbor ) { - prenotch.create( exec_space, _mu, particles, neigh_list ); - } + _timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); + + const auto x = particles.sliceReferencePosition(); + auto u = particles.sliceDisplacement(); + const auto vol = particles.sliceVolume(); + auto m = particles.sliceWeightedVolume(); + auto theta = particles.sliceDilatation(); + auto model = _model; + + Cabana::deep_copy( theta, 0.0 ); - auto getBrokenBonds() const { return _mu; } + auto dilatation = KOKKOS_LAMBDA( const int i ) + { + std::size_t num_neighbors = + Cabana::NeighborList::numNeighbor( + neigh_list, i ); + for ( std::size_t n = 0; n < num_neighbors; n++ ) + { + std::size_t j = + Cabana::NeighborList::getNeighbor( + neigh_list, i, n ); + + // Get the bond distance, displacement, and stretch. + double xi, r, s; + getDistance( x, u, i, j, xi, r, s ); + + // Check if all bonds are broken (m=0) to avoid dividing by + // zero. Alternatively, one could check if this bond mu(i,n) is + // broken, because m=0 only occurs when all bonds are broken. + // mu is still included to account for individual bond breaking. + if ( m( i ) > 0 ) + theta( i ) += mu( i, n ) * model.dilatation( + i, s, xi, vol( j ), m( i ) ); + } + }; + + neighbor.iterateLinear( exec_space{}, dilatation, particles, + "CabanaPD::Dilatation::compute" ); + + _timer.stop(); + } }; /****************************************************************************** Force free functions. ******************************************************************************/ -template +template void computeForce( ForceType& force, ParticleType& particles, - const ParallelType& neigh_op_tag, const bool reset = true ) + NeighborType& neighbor, const bool reset = true ) { auto x = particles.sliceReferencePosition(); auto u = particles.sliceDisplacement(); @@ -312,16 +377,17 @@ void computeForce( ForceType& force, ParticleType& particles, // neigh_op_tag ); // Forces only atomic if using team threading. - if ( std::is_same::value ) - force.computeForceFull( f_a, x, u, particles, neigh_op_tag ); + if constexpr ( std::is_same::value ) + force.computeForceFull( f_a, x, u, particles, neighbor ); else - force.computeForceFull( f, x, u, particles, neigh_op_tag ); + force.computeForceFull( f, x, u, particles, neighbor ); Kokkos::fence(); } -template +template double computeEnergy( ForceType& force, ParticleType& particles, - const ParallelType& neigh_op_tag ) + const NeighborType& neighbor ) { double energy = 0.0; if constexpr ( is_energy_output::value ) @@ -339,7 +405,7 @@ double computeEnergy( ForceType& force, ParticleType& particles, // energy = computeEnergy_half( force, x, u, // neigh_op_tag ); // else - energy = force.computeEnergyFull( W, x, u, particles, neigh_op_tag ); + energy = force.computeEnergyFull( W, x, u, particles, neighbor ); Kokkos::fence(); } return energy; diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp index 6717cf05..6f5b7306 100644 --- a/src/CabanaPD_ForceModels.hpp +++ b/src/CabanaPD_ForceModels.hpp @@ -17,28 +17,14 @@ namespace CabanaPD { -struct BaseForceModel -{ - double delta; - - BaseForceModel( const double _delta ) - : delta( _delta ){}; - - // Only needed for models which store bond properties. - void updateBonds( const int, const int ) {} - - // No-op for temperature. - KOKKOS_INLINE_FUNCTION - void thermalStretch( double&, const int, const int ) const {} -}; template class BasePlasticity { protected: using memory_space = MemorySpace; - using NeighborView = typename Kokkos::View; - NeighborView _s_p; + using neighbor_view = typename Kokkos::View; + neighbor_view _s_p; public: // Must update later because number of neighbors not known at construction. @@ -72,7 +58,11 @@ struct BaseTemperatureModel temperature = model.temperature; } - void update( const TemperatureType _temp ) { temperature = _temp; } + template + void update( const ParticleType& particles ) + { + temperature = particles.sliceTemperature(); + } // Update stretch with temperature effects. KOKKOS_INLINE_FUNCTION @@ -121,6 +111,11 @@ template struct ForceModel; +template +struct ForceDensityModel; + } // namespace CabanaPD #endif diff --git a/src/CabanaPD_HeatTransfer.hpp b/src/CabanaPD_HeatTransfer.hpp index 5111fa49..7835224f 100644 --- a/src/CabanaPD_HeatTransfer.hpp +++ b/src/CabanaPD_HeatTransfer.hpp @@ -13,6 +13,7 @@ #define HEATTRANSFER_H #include +#include #include #include @@ -27,39 +28,31 @@ class HeatTransfer; template class HeatTransfer> - : public Force { public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; using model_type = ForceModel; - using base_type = Force; - using neighbor_list_type = typename base_type::neighbor_list_type; protected: - using base_type::_half_neigh; - using base_type::_neigh_list; - using base_type::_timer; - Timer _euler_timer = base_type::_energy_timer; + Timer _timer; + Timer _euler_timer; model_type _model; public: // Running with mechanics as well; no reason to rebuild neighbors. - template - HeatTransfer( const bool half_neigh, const ForceType& force, - const model_type model ) - : base_type( half_neigh, force.getNeighbors() ) - , _model( model ) + HeatTransfer( const model_type model ) + : _model( model ) { } template + class NeighborType> void computeHeatTransferFull( TemperatureType& conduction, const PosType& x, const PosType& u, const ParticleType& particles, - ParallelType& neigh_op_tag ) + NeighborType& neighbor ) { _timer.start(); @@ -77,12 +70,8 @@ class HeatTransfer policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, temp_func, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::HeatTransfer::computeFull" ); - + neighbor.iterate( exec_space{}, temp_func, particles, + "CabanaPD::HeatTransfer::computeFull" ); _timer.stop(); } @@ -111,9 +100,7 @@ class HeatTransfer> : public HeatTransfer>, - BaseFracture - + DynamicTemperature, ModelParams...>> { public: // Using the default exec_space. @@ -124,47 +111,39 @@ class HeatTransfer>; - using neighbor_list_type = typename base_type::neighbor_list_type; protected: using base_type::_euler_timer; - using base_type::_half_neigh; - using base_type::_neigh_list; using base_type::_timer; model_type _model; - using fracture_type = BaseFracture; - using fracture_type::_mu; - public: // Explicit base model construction is necessary because of the indirect // model inheritance. This could be avoided with a BaseHeatTransfer object. - template - HeatTransfer( const bool half_neigh, const ForceType& force, - const model_type model ) - : base_type( half_neigh, force, - typename base_type::model_type( - PMB{}, NoFracture{}, model.delta, model.K, - model.temperature, model.kappa, model.cp, model.alpha, - model.temp0, model.constant_microconductivity ) ) - , fracture_type( force.getBrokenBonds() ) + HeatTransfer( const model_type model ) + : base_type( typename base_type::model_type( + PMB{}, NoFracture{}, model.delta, model.K, model.temperature, + model.kappa, model.cp, model.alpha, model.temp0, + model.constant_microconductivity ) ) , _model( model ) { } template + class NeighborType> void computeHeatTransferFull( TemperatureType& conduction, const PosType& x, const PosType& u, - const ParticleType& particles, ParallelType& ) + const ParticleType& particles, + NeighborType& neighbor ) { _timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); auto model = _model; - const auto neigh_list = _neigh_list; const auto vol = particles.sliceVolume(); const auto temp = particles.sliceTemperature(); - const auto mu = _mu; auto temp_func = KOKKOS_LAMBDA( const int i ) { @@ -191,19 +170,17 @@ class HeatTransfer policy( particles.frozenOffset(), - particles.localOffset() ); - Kokkos::parallel_for( "CabanaPD::HeatTransfer::computeFull", policy, - temp_func ); + neighbor.iterateLinear( exec_space{}, temp_func, particles, + "CabanaPD::HeatTransfer::computeFull" ); _timer.stop(); } }; // Heat transfer free functions. -template +template void computeHeatTransfer( HeatTransferType& heat_transfer, - ParticleType& particles, - const ParallelType& neigh_op_tag, const double dt ) + ParticleType& particles, const NeighborType& neighbor, + const double dt ) { auto x = particles.sliceReferencePosition(); auto u = particles.sliceDisplacement(); @@ -214,12 +191,12 @@ void computeHeatTransfer( HeatTransferType& heat_transfer, Cabana::deep_copy( conduction, 0.0 ); // Temperature only needs to be atomic if using team threading. - if ( std::is_same::value ) + if ( std::is_same::value ) heat_transfer.computeHeatTransferFull( conduction_a, x, u, particles, - neigh_op_tag ); + neighbor ); else heat_transfer.computeHeatTransferFull( conduction, x, u, particles, - neigh_op_tag ); + neighbor ); Kokkos::fence(); heat_transfer.forwardEuler( particles, dt ); diff --git a/src/CabanaPD_Neighbor.hpp b/src/CabanaPD_Neighbor.hpp new file mode 100644 index 00000000..9c4efc41 --- /dev/null +++ b/src/CabanaPD_Neighbor.hpp @@ -0,0 +1,225 @@ +/**************************************************************************** + * Copyright (c) 2022 by Oak Ridge National Laboratory * + * All rights reserved. * + * * + * This file is part of CabanaPD. CabanaPD 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 NEIGHBOR_H +#define NEIGHBOR_H + +#include + +#include + +namespace CabanaPD +{ +template +class Neighbor; + +template +class Neighbor +{ + public: + using list_type = + Cabana::VerletList; + using Tag = Cabana::SerialOpTag; + Tag tag; + + protected: + using exec_space = typename MemorySpace::execution_space; + bool _half_neigh; + list_type _neigh_list; + double mesh_max[3]; + double mesh_min[3]; + + Timer _timer; + + public: + // Primary constructor: use positions and construct neighbors. + template + Neighbor( const bool half_neigh, const ModelType& model, + const ParticleType& particles, const double tol = 1e-14 ) + : _half_neigh( half_neigh ) + { + _timer.start(); + if constexpr ( !is_contact::value ) + _neigh_list = list_type( + particles.sliceReferencePosition(), particles.frozenOffset(), + particles.localOffset(), model.cutoff() + tol, 1.0, + particles.ghost_mesh_lo, particles.ghost_mesh_hi ); + else + _neigh_list = list_type( + particles.sliceCurrentPosition(), particles.frozenOffset(), + particles.localOffset(), model.cutoff() + tol, 1.0, + particles.ghost_mesh_lo, particles.ghost_mesh_hi ); + _timer.stop(); + } + + unsigned getMaxLocal() + { + auto neigh = _neigh_list; + unsigned local_max_neighbors; + auto neigh_max = KOKKOS_LAMBDA( const int, unsigned& max_n ) + { + max_n = Cabana::NeighborList::maxNeighbor( neigh ); + }; + Kokkos::RangePolicy policy( 0, 1 ); + Kokkos::parallel_reduce( policy, neigh_max, local_max_neighbors ); + Kokkos::fence(); + return local_max_neighbors; + } + + void getStatistics( unsigned& max_neighbors, + unsigned long long& total_neighbors ) + { + auto neigh = _neigh_list; + unsigned local_max_neighbors; + unsigned long long local_total_neighbors; + auto neigh_stats = KOKKOS_LAMBDA( const int, unsigned& max_n, + unsigned long long& total_n ) + { + max_n = Cabana::NeighborList::maxNeighbor( neigh ); + total_n = Cabana::NeighborList::totalNeighbor( neigh ); + }; + Kokkos::RangePolicy policy( 0, 1 ); + Kokkos::parallel_reduce( policy, neigh_stats, local_max_neighbors, + local_total_neighbors ); + Kokkos::fence(); + MPI_Reduce( &local_max_neighbors, &max_neighbors, 1, MPI_UNSIGNED, + MPI_MAX, 0, MPI_COMM_WORLD ); + MPI_Reduce( &local_total_neighbors, &total_neighbors, 1, + MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD ); + } + + auto list() const { return _neigh_list; } + + template + void update( ParticleType& particles, const double cutoff ) + { + const auto y = particles.sliceCurrentPosition(); + const int n_frozen = particles.frozenOffset(); + const int n_local = particles.localOffset(); + + _neigh_list.build( exec_space{}, y, n_frozen, n_local, cutoff, 1.0, + particles.ghost_mesh_lo, particles.ghost_mesh_hi ); + } + + // Only rebuild neighbor list as needed. + template + void update( const ParticleType& particles, const double search_radius, + const double radius_extend, const bool require_update = false ) + { + double max_displacement = particles.getMaxDisplacement(); + if ( max_displacement > radius_extend || require_update ) + { + _timer.start(); + const auto y = particles.sliceCurrentPosition(); + _neigh_list.build( y, particles.frozenOffset(), + particles.localOffset(), search_radius, 1.0, + mesh_min, mesh_max ); + // Reset neighbor update displacement. + const auto u = particles.sliceDisplacement(); + auto u_neigh = particles.sliceDisplacementNeighborBuild(); + Cabana::deep_copy( u_neigh, u ); + _timer.stop(); + } + } + + template + void iterate( ExecSpace, const FunctorType functor, + const ParticleType& particles, const std::string label ) const + { + auto policy = makePolicy( particles ); + Cabana::neighbor_parallel_for( policy, functor, _neigh_list, + Cabana::FirstNeighborsTag(), tag, + label ); + } + + template + void iterateLinear( ExecSpace, const FunctorType functor, + const ParticleType& particles, + const std::string label ) const + { + auto policy = makePolicy( particles ); + Kokkos::parallel_for( label, policy, functor ); + } + + template + auto reduce( ExecSpace, const FunctorType functor, + const ParticleType& particles, const std::string label ) const + { + auto policy = makePolicy( particles ); + double total = 0.0; + Cabana::neighbor_parallel_reduce( policy, functor, _neigh_list, + Cabana::FirstNeighborsTag(), tag, + total, label ); + return total; + } + + template + auto reduceLinear( ExecSpace, const FunctorType functor, + const ParticleType& particles, + const std::string label ) const + { + double total = 0.0; + auto policy = makePolicy( particles ); + Kokkos::parallel_reduce( label, policy, functor, total ); + return total; + } + + template + auto makePolicy( const ParticleType& particles ) const + { + return Kokkos::RangePolicy( particles.frozenOffset(), + particles.localOffset() ); + } + + auto time() { return _timer.time(); }; +}; + +template +class Neighbor : public Neighbor +{ + public: + using memory_space = MemorySpace; + using neighbor_view = typename Kokkos::View; + + protected: + using base_type = Neighbor; + using base_type::_neigh_list; + + neighbor_view _mu; + + public: + template + Neighbor( const bool half_neigh, const ModelType& model, + const ParticleType& particles, const double tol = 1e-14 ) + : base_type( half_neigh, model, particles, tol ) + { + // Create View to track broken bonds. + // TODO: this could be optimized to ignore frozen particle bonds. + _mu = neighbor_view( + Kokkos::ViewAllocateWithoutInitializing( "broken_bonds" ), + particles.localOffset(), base_type::getMaxLocal() ); + Kokkos::deep_copy( _mu, 1 ); + } + + template + void prenotch( ExecSpace exec_space, const ParticleType& particles, + PrenotchType& prenotch ) + { + prenotch.create( exec_space, _mu, particles, _neigh_list ); + } + + auto brokenBonds() const { return _mu; } +}; + +} // namespace CabanaPD + +#endif diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp index 19811cab..362acbb6 100644 --- a/src/CabanaPD_Particles.hpp +++ b/src/CabanaPD_Particles.hpp @@ -79,15 +79,17 @@ namespace CabanaPD { template + class OutputType = BaseOutput, class DensityType = StaticDensity, + int Dimension = 3> class Particles; template -class Particles +class Particles { public: using self_type = Particles; + BaseOutput, StaticDensity, Dimension>; using thermal_type = TemperatureIndependent; using output_type = BaseOutput; using memory_space = MemorySpace; @@ -744,7 +746,7 @@ class Particles Cabana::Experimental::HDF5ParticleOutput::writeTimeStep( h5_config, "particles", MPI_COMM_WORLD, output_step, output_time, localOffset(), getPosition( use_reference ), sliceForce(), - sliceDisplacement(), sliceVelocity(), + sliceDisplacement(), sliceVelocity(), sliceDensity(), std::forward( other )... ); #else #ifdef Cabana_ENABLE_SILO @@ -752,7 +754,7 @@ class Particles writePartialRangeTimeStep( "particles", local_grid->globalGrid(), output_step, output_time, 0, localOffset(), getPosition( use_reference ), sliceForce(), - sliceDisplacement(), sliceVelocity(), + sliceDisplacement(), sliceVelocity(), sliceDensity(), std::forward( other )... ); #else @@ -790,15 +792,16 @@ class Particles }; template -class Particles +class Particles : public Particles + StaticDensity, Dimension> { public: using self_type = Particles; + BaseOutput, StaticDensity, Dimension>; using base_type = Particles; + BaseOutput, StaticDensity, Dimension>; using output_type = typename base_type::output_type; using thermal_type = TemperatureIndependent; using memory_space = typename base_type::memory_space; @@ -909,15 +912,15 @@ class Particles template class Particles + StaticDensity, Dimension> : public Particles + BaseOutput, StaticDensity, Dimension> { public: using self_type = Particles; + BaseOutput, StaticDensity, Dimension>; using base_type = Particles; + BaseOutput, StaticDensity, Dimension>; using thermal_type = TemperatureDependent; using output_type = typename base_type::output_type; using memory_space = typename base_type::memory_space; @@ -1022,17 +1025,19 @@ class Particles -class Particles - : public Particles +class Particles + : public Particles { // Note: no overloaded output() since there are very few cases where this // is a desired output field. public: - using self_type = - Particles; - using base_type = - Particles; + using self_type = Particles; + using base_type = Particles; using thermal_type = typename base_type::thermal_type; using output_type = typename base_type::output_type; using memory_space = typename base_type::memory_space; @@ -1112,15 +1117,16 @@ class Particles }; template -class Particles +class Particles : public Particles + StaticDensity, Dimension> { public: - using self_type = - Particles; - using base_type = - Particles; + using self_type = Particles; + using base_type = Particles; using thermal_type = typename base_type::thermal_type; using output_type = EnergyOutput; using memory_space = typename base_type::memory_space; @@ -1234,17 +1240,18 @@ class Particles aosoa_output_type _aosoa_output; }; -template +template class Particles + DensityType, Dimension> : public Particles + DensityType, Dimension> { public: using self_type = Particles; - using base_type = - Particles; + EnergyStressOutput, DensityType, Dimension>; + using base_type = Particles; using thermal_type = typename base_type::thermal_type; using output_type = EnergyStressOutput; using memory_space = typename base_type::memory_space; @@ -1330,6 +1337,86 @@ class Particles +class Particles + : public Particles +{ + public: + using self_type = Particles; + using base_type = Particles; + using thermal_type = typename base_type::thermal_type; + using output_type = OutputType; + using memory_space = typename base_type::memory_space; + using base_type::dim; + + using aosoa_density_type = + Cabana::AoSoA, memory_space, 1>; + + template + Particles( MemorySpace space, ModelType, TemperatureDependent temp, + DynamicDensity, Args&&... args ) + : base_type( space, LPS{}, temp, std::forward( args )... ) + { + _aosoa_density = aosoa_density_type( "Particle Output Fields", + base_type::localOffset() ); + init(); + } + + template + void createParticles( Args&&... args ) + { + // Forward arguments to standard or custom particle creation. + base_type::createParticles( std::forward( args )... ); + _aosoa_density.resize( base_type::localOffset() ); + } + + auto sliceCurrentDensity() + { + return Cabana::slice<0>( _aosoa_density, "current_density" ); + } + auto sliceCurrentDensity() const + { + return Cabana::slice<0>( _aosoa_density, "current_density" ); + } + + template + void resize( Args&&... args ) + { + base_type::resize( std::forward( args )... ); + _aosoa_density.resize( base_type::localOffset() ); + } + + template + void output( const int output_step, const double output_time, + const bool use_reference, OtherFields&&... other ) + { + base_type::output( output_step, output_time, use_reference, + sliceCurrentDensity(), + std::forward( other )... ); + } + + friend class Comm; + friend class Comm; + friend class Comm; + friend class Comm; + + protected: + void init() + { + auto reference = base_type::sliceDensity(); + auto current = sliceCurrentDensity(); + Cabana::deep_copy( current, reference ); + } + + aosoa_density_type _aosoa_density; +}; + /****************************************************************************** Template deduction guides. ******************************************************************************/ @@ -1362,6 +1449,12 @@ Particles( MemorySpace, ModelType, ThermalType, OutputType, -> Particles; +template +Particles( MemorySpace, ModelType, ThermalType, DensityType, Args... ) + -> Particles; + template Particles( MemorySpace, ModelType, OutputType, std::array, diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index ee68329d..970799f8 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -82,6 +82,7 @@ #include #include #include +#include namespace CabanaPD { @@ -95,15 +96,19 @@ class Solver // Core module types - required for all problems. using force_model_type = ForceModelType; - using force_type = Force; + using force_type = Force; using comm_type = Comm; - using neigh_iter_tag = Cabana::SerialOpTag; + using neighbor_type = + Neighbor; // Optional module types. using heat_transfer_type = HeatTransfer; - using contact_type = Force; + using contact_type = Force; using contact_model_type = ContactModelType; // Flexible module types. @@ -130,10 +135,7 @@ class Solver { setup( force_model ); - _neighbor_timer.start(); - contact = std::make_shared( inputs["half_neigh"], - particles, contact_model ); - _neighbor_timer.stop(); + contact = std::make_shared( contact_model ); } void setup( force_model_type force_model ) @@ -163,18 +165,21 @@ class Solver // Update temperature ghost size if needed. if constexpr ( is_temperature_dependent< typename force_model_type::thermal_type>::value ) - force_model.update( particles.sliceTemperature() ); + force_model.update( particles ); + + neighbor = std::make_shared( inputs["half_neigh"], + force_model, particles ); + // Plastic models need correct size for the bond array. + force_model.updateBonds( particles.localOffset(), + neighbor->getMaxLocal() ); - _neighbor_timer.start(); // This will either be PD or DEM forces. - force = std::make_shared( inputs["half_neigh"], particles, - force_model ); - _neighbor_timer.stop(); + force = std::make_shared( force_model ); _init_timer.start(); unsigned max_neighbors; unsigned long long total_neighbors; - force->getNeighborStatistics( max_neighbors, total_neighbors ); + neighbor->getStatistics( max_neighbors, total_neighbors ); // Create heat transfer if needed, using the same neighbor list as // the mechanics. @@ -182,8 +187,7 @@ class Solver typename force_model_type::thermal_type>::value ) { thermal_subcycle_steps = inputs["thermal_subcycle_steps"]; - heat_transfer = std::make_shared( - inputs["half_neigh"], *force, force_model ); + heat_transfer = std::make_shared( force_model ); } print = print_rank(); @@ -221,13 +225,13 @@ class Solver if constexpr ( !is_fracture< typename force_model_type::fracture_type>::value ) { - force->computeWeightedVolume( particles, neigh_iter_tag{} ); + force->computeWeightedVolume( particles, *neighbor ); comm->gatherWeightedVolume(); } // Compute initial internal forces and energy. updateForce(); - computeEnergy( *force, particles, neigh_iter_tag() ); - computeStress( *force, particles, neigh_iter_tag() ); + computeEnergy( *force, particles, *neighbor ); + computeStress( *force, particles, *neighbor ); if ( initial_output ) particles.output( 0, 0.0, output_reference ); @@ -305,7 +309,7 @@ class Solver // FIXME: Will need to rebuild ghosts. } - void updateNeighbors() { force->update( *particles, 0.0, true ); } + void updateNeighbors() { force->update( particles, 0.0, true ); } template void run( BoundaryType boundary_condition ) @@ -327,8 +331,7 @@ class Solver typename force_model_type::thermal_type>::value ) { if ( step % thermal_subcycle_steps == 0 ) - computeHeatTransfer( *heat_transfer, particles, - neigh_iter_tag{}, + computeHeatTransfer( *heat_transfer, particles, *neighbor, thermal_subcycle_steps * dt ); } @@ -344,7 +347,7 @@ class Solver updateForce(); if constexpr ( is_contact::value ) - computeForce( *contact, particles, neigh_iter_tag{}, false ); + computeForce( *contact, particles, *neighbor, false ); // Add force boundary condition. if ( boundary_condition.forceUpdate() ) @@ -379,7 +382,7 @@ class Solver updateForce(); if constexpr ( is_contact::value ) - computeForce( *contact, particles, neigh_iter_tag{}, false ); + computeForce( *contact, particles, *neighbor, false ); if constexpr ( is_temperature_dependent< typename force_model_type::thermal_type>::value ) @@ -404,15 +407,15 @@ class Solver if constexpr ( is_fracture< typename force_model_type::fracture_type>::value ) { - force->computeWeightedVolume( particles, neigh_iter_tag{} ); + force->computeWeightedVolume( particles, *neighbor ); comm->gatherWeightedVolume(); } // Compute and communicate dilatation for LPS (does nothing for PMB). - force->computeDilatation( particles, neigh_iter_tag{} ); + force->computeDilatation( particles, *neighbor ); comm->gatherDilatation(); // Compute internal forces. - computeForce( *force, particles, neigh_iter_tag{} ); + computeForce( *force, particles, *neighbor ); } void output( const int step ) @@ -420,8 +423,8 @@ class Solver // Print output. if ( step % output_frequency == 0 ) { - auto W = computeEnergy( *force, particles, neigh_iter_tag() ); - computeStress( *force, particles, neigh_iter_tag() ); + auto W = computeEnergy( *force, particles, *neighbor ); + computeStress( *force, particles, *neighbor ); particles.output( step / output_frequency, step * dt, output_reference ); @@ -438,11 +441,11 @@ class Solver { // Output after construction and initial forces. std::ofstream out( output_file, std::ofstream::app ); - _init_time += _init_timer.time() + _neighbor_timer.time() + + _init_time += _init_timer.time() + neighbor->time() + particles.timeInit() + comm->timeInit() + integrator->timeInit() + boundary_init_time; log( out, "Init-Time(s): ", _init_time ); - log( out, "Init-Neighbor-Time(s): ", _neighbor_timer.time(), "\n" ); + log( out, "Init-Neighbor-Time(s): ", neighbor->time(), "\n" ); log( out, "#Timestep/Total-steps Simulation-time Total-strain-energy " "Step-Time(s) Force-Time(s) Comm-Time(s) Integrate-Time(s) " "Energy-Time(s) Output-Time(s) Particle*steps/s" ); @@ -461,7 +464,7 @@ class Solver double integrate_time = integrator->time(); double force_time = force->time(); double energy_time = force->timeEnergy(); - double neigh_time = force->timeNeighbor(); + double neigh_time = neighbor->time(); double output_time = particles.timeOutput(); _total_time += step_time; // Instantaneous rate. @@ -489,7 +492,7 @@ class Solver double force_time = force->time(); double energy_time = force->timeEnergy(); double output_time = particles.timeOutput(); - double neighbor_time = _neighbor_timer.time(); + double neighbor_time = neighbor->time(); _total_time = _init_time + comm_time + integrate_time + force_time + energy_time + output_time + particles.time(); @@ -534,7 +537,7 @@ class Solver "Cannot create prenotch in system without fracture." ); // Create prenotch. - force->prenotch( exec_space{}, particles, prenotch ); + neighbor->prenotch( exec_space{}, particles, prenotch ); _init_time += prenotch.time(); } @@ -542,6 +545,7 @@ class Solver Inputs inputs; std::shared_ptr comm; std::shared_ptr integrator; + std::shared_ptr neighbor; std::shared_ptr force; // Optional modules. std::shared_ptr heat_transfer; @@ -554,7 +558,6 @@ class Solver // Note: init_time is combined from many class timers. double _init_time; Timer _init_timer; - Timer _neighbor_timer; Timer _step_timer; double _total_time; bool print; diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp index 4e47fbc4..3da76b36 100644 --- a/src/CabanaPD_Types.hpp +++ b/src/CabanaPD_Types.hpp @@ -35,9 +35,15 @@ struct is_fracture : public std::true_type // Mechanics tags. struct Elastic { + using base_type = Elastic; +}; +struct Plastic +{ + using base_type = Plastic; }; struct ElasticPerfectlyPlastic { + using base_type = Plastic; }; // Model category tags. @@ -207,5 +213,12 @@ struct is_stress_output : public std::true_type { }; +struct StaticDensity +{ +}; +struct DynamicDensity +{ +}; + } // namespace CabanaPD #endif diff --git a/src/force/CabanaPD_Contact.hpp b/src/force/CabanaPD_Contact.hpp index 3cfec5a9..7f0bca8e 100644 --- a/src/force/CabanaPD_Contact.hpp +++ b/src/force/CabanaPD_Contact.hpp @@ -38,98 +38,30 @@ KOKKOS_INLINE_FUNCTION void getRelativeNormalVelocityComponents( vn /= r; }; -// Contact forces base class. -template -class BaseForceContact : public Force -{ - public: - using base_type = Force; - using neighbor_list_type = typename base_type::neighbor_list_type; - - // NOTE: using 2x radius to find neighbors when particles first touch. - template - BaseForceContact( const bool half_neigh, const ParticleType& particles, - const ModelType model ) - : base_type( half_neigh, 2.0 * model.radius + model.radius_extend, - particles.sliceCurrentPosition(), particles.frozenOffset(), - particles.localOffset(), particles.ghost_mesh_lo, - particles.ghost_mesh_hi ) - , search_radius( 2.0 * model.radius + model.radius_extend ) - , radius_extend( model.radius_extend ) - { - for ( int d = 0; d < particles.dim; d++ ) - { - mesh_min[d] = particles.ghost_mesh_lo[d]; - mesh_max[d] = particles.ghost_mesh_hi[d]; - } - } - - // Only rebuild neighbor list as needed. - template - void update( const ParticleType& particles, - const bool require_update = false ) - { - double max_displacement = particles.getMaxDisplacement(); - if ( max_displacement > radius_extend || require_update ) - { - _neigh_timer.start(); - const auto y = particles.sliceCurrentPosition(); - _neigh_list.build( y, particles.frozenOffset(), - particles.localOffset(), search_radius, 1.0, - mesh_min, mesh_max ); - // Reset neighbor update displacement. - const auto u = particles.sliceDisplacement(); - auto u_neigh = particles.sliceDisplacementNeighborBuild(); - Cabana::deep_copy( u_neigh, u ); - _neigh_timer.stop(); - } - } - - auto timeNeighbor() { return _neigh_timer.time(); }; - - protected: - double search_radius; - double radius_extend; - Timer _neigh_timer; - - using base_type::_half_neigh; - using base_type::_neigh_list; - double mesh_max[3]; - double mesh_min[3]; -}; - /****************************************************************************** Normal repulsion forces ******************************************************************************/ template -class Force - : public BaseForceContact +class Force + : public BaseForce { public: - using base_type = BaseForceContact; - using neighbor_list_type = typename base_type::neighbor_list_type; - - template - Force( const bool half_neigh, const ParticleType& particles, - const NormalRepulsionModel model ) - : base_type( half_neigh, particles, model ) - , _model( model ) + using base_type = BaseForce; + + Force( const NormalRepulsionModel model ) + : _model( model ) { } template + class NeighborType> void computeForceFull( ForceType& fc, const PosType& x, const PosType& u, - const ParticleType& particles, - ParallelType& neigh_op_tag ) + ParticleType& particles, NeighborType& neighbor ) { auto model = _model; const auto vol = particles.sliceVolume(); - const auto y = particles.sliceCurrentPosition(); - const int n_frozen = particles.frozenOffset(); - const int n_local = particles.localOffset(); - base_type::update( particles, particles.getMaxDisplacement() ); + neighbor.update( particles, model.cutoff(), model.extend() ); auto contact_full = KOKKOS_LAMBDA( const int i, const int j ) { @@ -157,32 +89,24 @@ class Force // FIXME: using default space for now. using exec_space = typename MemorySpace::execution_space; - Kokkos::RangePolicy policy( n_frozen, n_local ); - Cabana::neighbor_parallel_for( - policy, contact_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::Contact::compute_full" ); + neighbor.iterate( exec_space{}, contact_full, particles, + "CabanaPD::Contact::compute_full" ); _timer.stop(); } // FIXME: implement energy template + class NeighborType> double computeEnergyFull( WType&, const PosType&, const PosType&, - ParticleType&, const int, ParallelType& ) + ParticleType&, const int, NeighborType& ) { return 0.0; } protected: NormalRepulsionModel _model; - using base_type::_half_neigh; - using base_type::_neigh_list; - using base_type::_neigh_timer; using base_type::_timer; - - double mesh_max[3]; - double mesh_min[3]; }; } // namespace CabanaPD diff --git a/src/force/CabanaPD_Hertzian.hpp b/src/force/CabanaPD_Hertzian.hpp index 8011e1d7..7abeed1d 100644 --- a/src/force/CabanaPD_Hertzian.hpp +++ b/src/force/CabanaPD_Hertzian.hpp @@ -25,35 +25,28 @@ namespace CabanaPD Normal repulsion forces ******************************************************************************/ template -class Force : public BaseForceContact +class Force + : public BaseForce { public: - using base_type = BaseForceContact; - using neighbor_list_type = typename base_type::neighbor_list_type; - - template - Force( const bool half_neigh, const ParticleType& particles, - const HertzianModel model ) - : base_type( half_neigh, particles, model ) - , _model( model ) + using base_type = BaseForce; + + Force( const HertzianModel model ) + : _model( model ) { } template + class NeighborType> void computeForceFull( ForceType& fc, const PosType& x, const PosType& u, - const ParticleType& particles, - ParallelType& neigh_op_tag ) + ParticleType& particles, NeighborType& neighbor ) { - const int n_frozen = particles.frozenOffset(); - const int n_local = particles.localOffset(); - auto model = _model; const auto vol = particles.sliceVolume(); const auto rho = particles.sliceDensity(); const auto vel = particles.sliceVelocity(); - base_type::update( particles, particles.getMaxDisplacement() ); + neighbor.update( particles, model.cutoff(), model.extend() ); auto contact_full = KOKKOS_LAMBDA( const int i, const int j ) { @@ -76,11 +69,8 @@ class Force : public BaseForceContact // FIXME: using default space for now. using exec_space = typename MemorySpace::execution_space; - Kokkos::RangePolicy policy( n_frozen, n_local ); - Cabana::neighbor_parallel_for( - policy, contact_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::Contact::compute_full" ); - + neighbor.iterate( exec_space{}, contact_full, particles, + "CabanaPD::Hertzian::compute" ); _timer.stop(); } @@ -95,9 +85,6 @@ class Force : public BaseForceContact protected: HertzianModel _model; - using base_type::_half_neigh; - using base_type::_neigh_list; - using base_type::_neigh_timer; using base_type::_timer; }; diff --git a/src/force/CabanaPD_LPS.hpp b/src/force/CabanaPD_LPS.hpp index a46b4188..9bfe5f36 100644 --- a/src/force/CabanaPD_LPS.hpp +++ b/src/force/CabanaPD_LPS.hpp @@ -69,100 +69,33 @@ namespace CabanaPD { -template -class Force> - : public Force +template +class Force + : public Dilatation { protected: - using base_type = Force; - using base_type::_half_neigh; - using model_type = ForceModel; - model_type _model; + using base_type = Dilatation; + using model_type = typename base_type::model_type; + using base_type::_model; - using base_type::_energy_timer; - using base_type::_stress_timer; using base_type::_timer; + Timer _energy_timer; + Timer _stress_timer; public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; - using neighbor_list_type = typename base_type::neighbor_list_type; - using base_type::_neigh_list; - template - Force( const bool half_neigh, ParticleType& particles, - const model_type model ) - : base_type( half_neigh, model.delta, particles ) - , _model( model ) + Force( const model_type model ) + : base_type( model ) { } - template - void computeWeightedVolume( ParticleType& particles, - const ParallelType neigh_op_tag ) - { - _timer.start(); - - auto x = particles.sliceReferencePosition(); - auto u = particles.sliceDisplacement(); - const auto vol = particles.sliceVolume(); - auto m = particles.sliceWeightedVolume(); - Cabana::deep_copy( m, 0.0 ); - auto model = _model; - - auto weighted_volume = KOKKOS_LAMBDA( const int i, const int j ) - { - // Get the reference positions and displacements. - double xi, r, s; - getDistance( x, u, i, j, xi, r, s ); - m( i ) += model.weightedVolume( xi, vol( j ) ); - }; - - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, weighted_volume, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLPS::computeWeightedVolume" ); - - _timer.stop(); - } - - template - void computeDilatation( ParticleType& particles, - const ParallelType neigh_op_tag ) - { - _timer.start(); - - const auto x = particles.sliceReferencePosition(); - auto u = particles.sliceDisplacement(); - const auto vol = particles.sliceVolume(); - auto m = particles.sliceWeightedVolume(); - auto theta = particles.sliceDilatation(); - auto model = _model; - Cabana::deep_copy( theta, 0.0 ); - - auto dilatation = KOKKOS_LAMBDA( const int i, const int j ) - { - // Get the bond distance, displacement, and stretch. - double xi, r, s; - getDistance( x, u, i, j, xi, r, s ); - theta( i ) += model.dilatation( s, xi, vol( j ), m( i ) ); - }; - - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, dilatation, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLPS::computeDilatation" ); - - _timer.stop(); - } - template + class NeighborType> void computeForceFull( ForceType& f, const PosType& x, const PosType& u, const ParticleType& particles, - ParallelType& neigh_op_tag ) + NeighborType& neighbor ) { _timer.start(); @@ -191,29 +124,26 @@ class Force> f( i, 2 ) += fz_i; }; - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, force_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLPS::computeFull" ); - + neighbor.iterate( exec_space{}, force_full, particles, + "CabanaPD::ForceLPS::computeFull" ); _timer.stop(); } template + class NeighborType> double computeEnergyFull( WType& W, const PosType& x, const PosType& u, const ParticleType& particles, - ParallelType& neigh_op_tag ) + NeighborType& neighbor ) { _energy_timer.start(); auto model = _model; - auto neigh_list = _neigh_list; const auto vol = particles.sliceVolume(); const auto theta = particles.sliceDilatation(); auto m = particles.sliceWeightedVolume(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); auto energy_full = KOKKOS_LAMBDA( const int i, const int j, double& Phi ) @@ -231,30 +161,22 @@ class Force> Phi += w * vol( i ); }; - double strain_energy = 0.0; - - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_reduce( - policy, energy_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, strain_energy, - "CabanaPD::ForceLPS::computeEnergyFull" ); + double strain_energy = + neighbor.reduce( exec_space{}, energy_full, particles, + "CabanaPD::ForceLPS::computeEnergyFull" ); _energy_timer.stop(); return strain_energy; } - auto time() { return _timer.time(); }; auto timeEnergy() { return _energy_timer.time(); }; - template - void computeStressFull( ParticleType& particles, - ParallelType& neigh_op_tag ) + template + void computeStressFull( ParticleType& particles, NeighborType& neighbor ) { _stress_timer.start(); auto model = _model; - auto neigh_list = _neigh_list; const auto x = particles.sliceReferencePosition(); const auto u = particles.sliceDisplacement(); const auto vol = particles.sliceVolume(); @@ -291,157 +213,47 @@ class Force> stress( i, 2, 1 ) += fz_i * xi_y; }; - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, stress_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLPSFracture::computeStressFull" ); - + neighbor.iterate( exec_space{}, stress_full, particles, + "CabanaPD::ForceLPS::computeStressFull" ); _stress_timer.stop(); } }; -template -class Force> - : public Force>, - public BaseFracture +template +class Force + : public Dilatation { protected: - using fracture_type = BaseFracture; - using fracture_type::_mu; + using base_type = Dilatation; + using model_type = typename base_type::model_type; + using base_type::_model; - using base_type = Force>; - using base_type::_half_neigh; - using model_type = ForceModel; - model_type _model; - - using base_type::_energy_timer; - using base_type::_stress_timer; using base_type::_timer; + Timer _energy_timer; + Timer _stress_timer; public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; - using neighbor_list_type = typename base_type::neighbor_list_type; - using base_type::_neigh_list; - - template - Force( const bool half_neigh, const ParticleType& particles, - const model_type model ) - : base_type( half_neigh, particles, model ) - , fracture_type( particles.localOffset(), - base_type::getMaxLocalNeighbors() ) - , _model( model ) - { - } - - template - void prenotch( ExecSpace exec_space, const ParticleType& particles, - PrenotchType& prenotch ) - { - fracture_type::prenotch( exec_space, particles, prenotch, _neigh_list ); - } - - template - void computeWeightedVolume( ParticleType& particles, ParallelType ) - { - _timer.start(); - - auto x = particles.sliceReferencePosition(); - auto u = particles.sliceDisplacement(); - const auto vol = particles.sliceVolume(); - auto m = particles.sliceWeightedVolume(); - Cabana::deep_copy( m, 0.0 ); - auto model = _model; - auto neigh_list = _neigh_list; - auto mu = _mu; - - auto weighted_volume = KOKKOS_LAMBDA( const int i ) - { - std::size_t num_neighbors = - Cabana::NeighborList::numNeighbor( - neigh_list, i ); - for ( std::size_t n = 0; n < num_neighbors; n++ ) - { - std::size_t j = - Cabana::NeighborList::getNeighbor( - neigh_list, i, n ); - - // Get the reference positions and displacements. - double xi, r, s; - getDistance( x, u, i, j, xi, r, s ); - // mu is included to account for bond breaking. - m( i ) += mu( i, n ) * model.weightedVolume( xi, vol( j ) ); - } - }; - - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Kokkos::parallel_for( "CabanaPD::ForceLPSDamage::computeWeightedVolume", - policy, weighted_volume ); - _timer.stop(); - } - - template - void computeDilatation( ParticleType& particles, ParallelType ) + Force( const model_type model ) + : base_type( model ) { - _timer.start(); - - const auto x = particles.sliceReferencePosition(); - auto u = particles.sliceDisplacement(); - const auto vol = particles.sliceVolume(); - auto m = particles.sliceWeightedVolume(); - auto theta = particles.sliceDilatation(); - auto model = _model; - auto neigh_list = _neigh_list; - auto mu = _mu; - Cabana::deep_copy( theta, 0.0 ); - - auto dilatation = KOKKOS_LAMBDA( const int i ) - { - std::size_t num_neighbors = - Cabana::NeighborList::numNeighbor( - neigh_list, i ); - for ( std::size_t n = 0; n < num_neighbors; n++ ) - { - std::size_t j = - Cabana::NeighborList::getNeighbor( - neigh_list, i, n ); - - // Get the bond distance, displacement, and stretch. - double xi, r, s; - getDistance( x, u, i, j, xi, r, s ); - - // Check if all bonds are broken (m=0) to avoid dividing by - // zero. Alternatively, one could check if this bond mu(i,n) is - // broken, because m=0 only occurs when all bonds are broken. - // mu is still included to account for individual bond breaking. - if ( m( i ) > 0 ) - theta( i ) += mu( i, n ) * - model.dilatation( s, xi, vol( j ), m( i ) ); - } - }; - - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Kokkos::parallel_for( "CabanaPD::ForceLPSDamage::computeDilatation", - policy, dilatation ); - - _timer.stop(); } template + class NeighborType> void computeForceFull( ForceType& f, const PosType& x, const PosType& u, - const ParticleType& particles, ParallelType ) + const ParticleType& particles, + NeighborType& neighbor ) { _timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); auto break_coeff = _model.bond_break_coeff; auto model = _model; - auto neigh_list = _neigh_list; - auto mu = _mu; const auto vol = particles.sliceVolume(); auto theta = particles.sliceDilatation(); @@ -493,24 +305,24 @@ class Force> } }; - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Kokkos::parallel_for( "CabanaPD::ForceLPSDamage::computeFull", policy, - force_full ); + neighbor.iterateLinear( exec_space{}, force_full, particles, + "CabanaPD::ForceLPSDamage::computeFull" ); _timer.stop(); } template + class NeighborType> double computeEnergyFull( WType& W, const PosType& x, const PosType& u, - ParticleType& particles, ParallelType& ) + ParticleType& particles, + const NeighborType& neighbor ) { _energy_timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); auto model = _model; - auto neigh_list = _neigh_list; - auto mu = _mu; const auto vol = particles.sliceVolume(); const auto theta = particles.sliceDilatation(); @@ -548,30 +360,31 @@ class Force> phi( i ) = 1 - phi_i / vol_H_i; }; - double strain_energy = 0.0; - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Kokkos::parallel_reduce( "CabanaPD::ForceLPSDamage::computeEnergyFull", - policy, energy_full, strain_energy ); + double strain_energy = neighbor.reduceLinear( + exec_space{}, energy_full, particles, + "CabanaPD::ForceLPSDamage::computeEnergyFull" ); _energy_timer.stop(); return strain_energy; } - template - void computeStressFull( ParticleType& particles, ParallelType& ) + auto timeEnergy() { return _energy_timer.time(); }; + + template + void computeStressFull( ParticleType& particles, NeighborType& neighbor ) { _stress_timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); auto model = _model; - auto neigh_list = _neigh_list; const auto x = particles.sliceReferencePosition(); const auto u = particles.sliceDisplacement(); const auto vol = particles.sliceVolume(); const auto theta = particles.sliceDilatation(); const auto m = particles.sliceWeightedVolume(); auto stress = particles.sliceStress(); - auto mu = _mu; auto stress_full = KOKKOS_LAMBDA( const int i ) { @@ -618,47 +431,40 @@ class Force> } }; - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Kokkos::parallel_for( "CabanaPD::ForceLPSFracture::computeStressFull", - policy, stress_full ); - + neighbor.iterateLinear( exec_space{}, stress_full, particles, + "CabanaPD::ForceLPS::computeStressFull" ); _stress_timer.stop(); } }; -template -class Force> - : public Force> +template +class Force + : public Dilatation { protected: - using base_type = Force>; - using model_type = ForceModel; + using base_type = Dilatation; + using model_type = typename base_type::model_type; model_type _model; - using base_type::_energy_timer; - using base_type::_stress_timer; using base_type::_timer; + Timer _energy_timer; + Timer _stress_timer; public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; - using neighbor_list_type = typename base_type::neighbor_list_type; - using base_type::_neigh_list; - template - Force( const bool half_neigh, ParticleType& particles, - const model_type model ) - : base_type( half_neigh, particles, model ) + Force( const model_type model ) + : base_type( model ) , _model( model ) { } template + class NeighborType> void computeForceFull( ForceType& f, const PosType& x, const PosType& u, const ParticleType& particles, - ParallelType& neigh_op_tag ) + NeighborType& neighbor ) { _timer.start(); @@ -691,25 +497,22 @@ class Force> f( i, 2 ) += fz_i; }; - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, force_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLPS::computeFull" ); - + neighbor.iterate( exec_space{}, force_full, particles, + "CabanaPD::ForceLPS::computeFull" ); _timer.stop(); } template + class NeighborType> double computeEnergyFull( WType& W, const PosType& x, const PosType& u, const ParticleType& particles, - ParallelType& neigh_op_tag ) + NeighborType& neighbor ) { _energy_timer.start(); auto model = _model; - auto neigh_list = _neigh_list; + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); const auto vol = particles.sliceVolume(); const auto theta = particles.sliceDilatation(); @@ -732,27 +535,22 @@ class Force> Phi += w * vol( i ); }; - double strain_energy = 0.0; - - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_reduce( - policy, energy_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, strain_energy, - "CabanaPD::ForceLPS::computeEnergyFull" ); + double strain_energy = + neighbor.reduce( exec_space{}, energy_full, particles, + "CabanaPD::ForceLPS::computeEnergyFull" ); _energy_timer.stop(); return strain_energy; } - template - void computeStressFull( ParticleType& particles, - ParallelType& neigh_op_tag ) + auto timeEnergy() { return _energy_timer.time(); }; + + template + void computeStressFull( ParticleType& particles, NeighborType& neighbor ) { _stress_timer.start(); auto model = _model; - auto neigh_list = _neigh_list; const auto x = particles.sliceReferencePosition(); const auto u = particles.sliceDisplacement(); const auto vol = particles.sliceVolume(); @@ -789,12 +587,8 @@ class Force> stress( i, 2, 1 ) += fz_i * xi_y; }; - Kokkos::RangePolicy policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, stress_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLPS::computeStressFull" ); - + neighbor.iterate( exec_space{}, stress_full, particles, + "CabanaPD::ForceLPS::computeStressFull" ); _stress_timer.stop(); } }; diff --git a/src/force/CabanaPD_PMB.hpp b/src/force/CabanaPD_PMB.hpp index 3767c405..e0dda1c1 100644 --- a/src/force/CabanaPD_PMB.hpp +++ b/src/force/CabanaPD_PMB.hpp @@ -70,22 +70,17 @@ namespace CabanaPD { -template -class Force> - : public Force +template +class Force + : public BaseForce { public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; - using model_type = - ForceModel; - using base_type = Force; - using neighbor_list_type = typename base_type::neighbor_list_type; - using base_type::_neigh_list; + using model_type = ModelType; + using base_type = BaseForce; protected: - using base_type::_half_neigh; model_type _model; using base_type::_energy_timer; @@ -93,19 +88,16 @@ class Force - Force( const bool half_neigh, const ParticleType& particles, - const model_type model ) - : base_type( half_neigh, model.delta, particles ) - , _model( model ) + Force( const model_type model ) + : _model( model ) { } template + class NeighborType> void computeForceFull( ForceType& f, const PosType& x, const PosType& u, const ParticleType& particles, - ParallelType& neigh_op_tag ) + NeighborType& neighbor ) { _timer.start(); @@ -134,20 +126,16 @@ class Force policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, force_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForcePMB::computeFull" ); - + neighbor.iterate( exec_space{}, force_full, particles, + "CabanaPD::ForcePMB::computeFull" ); _timer.stop(); } template + class NeighborType> double computeEnergyFull( WType& W, const PosType& x, const PosType& u, const ParticleType& particles, - ParallelType& neigh_op_tag ) + NeighborType& neighbor ) { _energy_timer.start(); @@ -168,21 +156,15 @@ class Force policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_reduce( - policy, energy_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, strain_energy, - "CabanaPD::ForcePMB::computeEnergyFull" ); - + double strain_energy = + neighbor.reduce( exec_space{}, energy_full, particles, + "CabanaPD::ForcePMB::computeEnergyFull" ); _energy_timer.stop(); return strain_energy; } - template - void computeStressFull( ParticleType& particles, - ParallelType& neigh_op_tag ) + template + void computeStressFull( ParticleType& particles, NeighborType& neighbor ) { _stress_timer.start(); @@ -223,36 +205,23 @@ class Force policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, stress_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForcePMB::computeStressFull" ); - + neighbor.iterate( exec_space{}, stress_full, particles, + "CabanaPD::ForcePMB::computeStressFull" ); _stress_timer.stop(); } }; -template -class Force> - : public Force, - public BaseFracture +template +class Force + : public BaseForce { public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; - using model_type = ForceModel; - using base_type = Force; - using neighbor_list_type = typename base_type::neighbor_list_type; - using base_type::_neigh_list; + using model_type = ModelType; + using base_type = BaseForce; protected: - using fracture_type = BaseFracture; - using fracture_type::_mu; - - using base_model_type = typename model_type::base_type; - using base_type::_half_neigh; model_type _model; using base_type::_energy_timer; @@ -260,36 +229,24 @@ class Force - Force( const bool half_neigh, const ParticleType& particles, - const model_type model ) - : base_type( half_neigh, model.delta, particles ) - , fracture_type( particles.localOffset(), - base_type::getMaxLocalNeighbors() ) + Force( const model_type model ) + : base_type() , _model( model ) { - // Needed only for models which store per-bond information. - _model.updateBonds( particles.localOffset(), - base_type::getMaxLocalNeighbors() ); - } - - template - void prenotch( ExecSpace exec_space, const ParticleType& particles, - PrenotchType& prenotch ) - { - fracture_type::prenotch( exec_space, particles, prenotch, _neigh_list ); } template + class NeighborType> void computeForceFull( ForceType& f, const PosType& x, const PosType& u, - const ParticleType& particles, ParallelType& ) + const ParticleType& particles, + const NeighborType& neighbor ) { _timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); auto model = _model; - auto neigh_list = _neigh_list; - auto mu = _mu; const auto vol = particles.sliceVolume(); const auto nofail = particles.sliceNoFail(); @@ -337,24 +294,22 @@ class Force policy( particles.frozenOffset(), - particles.localOffset() ); - Kokkos::parallel_for( "CabanaPD::ForcePMBDamage::computeFull", policy, - force_full ); - + neighbor.iterateLinear( exec_space{}, force_full, particles, + "CabanaPD::ForcePMBDamage::computeFull" ); _timer.stop(); } template + class NeighborType> double computeEnergyFull( WType& W, const PosType& x, const PosType& u, - ParticleType& particles, ParallelType& ) + ParticleType& particles, NeighborType& neighbor ) { _energy_timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); auto model = _model; - auto neigh_list = _neigh_list; - auto mu = _mu; const auto vol = particles.sliceVolume(); auto phi = particles.sliceDamage(); @@ -386,29 +341,28 @@ class Force policy( particles.frozenOffset(), - particles.localOffset() ); - Kokkos::parallel_reduce( "CabanaPD::ForcePMBDamage::computeEnergyFull", - policy, energy_full, strain_energy ); + double strain_energy = neighbor.reduceLinear( + exec_space{}, energy_full, particles, + "CabanaPD::ForcePMBDamage::computeEnergyFull" ); _energy_timer.stop(); return strain_energy; } - template - void computeStressFull( ParticleType& particles, ParallelType& ) + template + void computeStressFull( ParticleType& particles, NeighborType& neighbor ) { _stress_timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); auto model = _model; - auto neigh_list = _neigh_list; const auto x = particles.sliceReferencePosition(); const auto u = particles.sliceDisplacement(); const auto vol = particles.sliceVolume(); const auto f = particles.sliceForce(); auto stress = particles.sliceStress(); - auto mu = _mu; auto stress_full = KOKKOS_LAMBDA( const int i ) { @@ -451,31 +405,120 @@ class Force policy( particles.frozenOffset(), - particles.localOffset() ); - Kokkos::parallel_for( "CabanaPD::ForcePMBDamage::computeStressFull", - policy, stress_full ); - + neighbor.iterateLinear( exec_space{}, stress_full, particles, + "CabanaPD::ForcePMBDamage::computeStressFull" ); _stress_timer.stop(); } }; -template -class Force> - : public Force +template +class Force + : public Force, + public Dilatation { public: // Using the default exec_space. using exec_space = typename MemorySpace::execution_space; - using model_type = - ForceModel; - using base_type = Force; - using neighbor_list_type = typename base_type::neighbor_list_type; - using base_type::_neigh_list; + using model_type = ModelType; + using base_type = Force; + using dilatation_type = Dilatation; + + using dilatation_type::computeDilatation; + using dilatation_type::computeWeightedVolume; + + protected: + using base_type::_model; + + using base_type::_energy_timer; + using base_type::_timer; + + public: + Force( const model_type model ) + : base_type( model ) + , dilatation_type( model ) + { + } + + template + void computeForceFull( ForceType& f, const PosType& x, const PosType& u, + const ParticleType& particles, + const NeighborType& neighbor ) + { + _timer.start(); + using neighbor_list_type = typename NeighborType::list_type; + auto neigh_list = neighbor.list(); + const auto mu = neighbor.brokenBonds(); + + auto model = _model; + const auto vol = particles.sliceVolume(); + const auto nofail = particles.sliceNoFail(); + auto theta = particles.sliceDilatation(); + + auto force_full = KOKKOS_LAMBDA( const int i ) + { + std::size_t num_neighbors = + Cabana::NeighborList::numNeighbor( + neigh_list, i ); + for ( std::size_t n = 0; n < num_neighbors; n++ ) + { + double fx_i = 0.0; + double fy_i = 0.0; + double fz_i = 0.0; + + std::size_t j = + Cabana::NeighborList::getNeighbor( + neigh_list, i, n ); + + // Get the reference positions and displacements. + double xi, r, s; + double rx, ry, rz; + getDistance( x, u, i, j, xi, r, s, rx, ry, rz ); + + model.thermalStretch( s, i, j ); + + model.updateDensity( i, theta( i ) ); + + // Break if beyond critical stretch unless in no-fail zone. + if ( model.criticalStretch( i, j, r, xi ) && !nofail( i ) && + !nofail( j ) ) + { + mu( i, n ) = 0; + } + // Else if statement is only for performance. + else if ( mu( i, n ) > 0 ) + { + const double coeff = model.forceCoeff( i, n, s, vol( j ) ); + + double muij = mu( i, n ); + fx_i = muij * coeff * rx / r; + fy_i = muij * coeff * ry / r; + fz_i = muij * coeff * rz / r; + + f( i, 0 ) += fx_i; + f( i, 1 ) += fy_i; + f( i, 2 ) += fz_i; + } + } + }; + + neighbor.iterateLinear( exec_space{}, force_full, particles, + "CabanaPD::ForcePMBDamage::computeFull" ); + _timer.stop(); + } +}; + +template +class Force + : public BaseForce +{ + public: + // Using the default exec_space. + using exec_space = typename MemorySpace::execution_space; + using model_type = ModelType; + using base_type = BaseForce; protected: - using base_type::_half_neigh; model_type _model; using base_type::_energy_timer; @@ -483,18 +526,15 @@ class Force - Force( const bool half_neigh, const ParticleType& particles, - const model_type model ) - : base_type( half_neigh, model.delta, particles ) - , _model( model ) + Force( const model_type model ) + : _model( model ) { } template + class NeighborType> void computeForceFull( ForceType& f, const PosType& x, const PosType& u, - ParticleType& particles, ParallelType& neigh_op_tag ) + ParticleType& particles, NeighborType& neighbor ) { _timer.start(); @@ -524,20 +564,15 @@ class Force policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, force_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLinearPMB::computeFull" ); - + neighbor.iterate( exec_space{}, force_full, particles, + "CabanaPD::ForceLinearPMB::computeFull" ); _timer.stop(); } template + class NeighborType> double computeEnergyFull( WType& W, const PosType& x, const PosType& u, - ParticleType& particles, - ParallelType& neigh_op_tag ) + ParticleType& particles, NeighborType& neighbor ) { _energy_timer.start(); @@ -558,26 +593,19 @@ class Force policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_reduce( - policy, energy_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, strain_energy, - "CabanaPD::ForceLinearPMB::computeEnergyFull" ); - + double strain_energy = + neighbor.reduce( exec_space{}, energy_full, particles, + "CabanaPD::ForceLinearPMB::computeEnergyFull" ); _energy_timer.stop(); return strain_energy; } - template - void computeStressFull( ParticleType& particles, - ParallelType& neigh_op_tag ) + template + void computeStressFull( ParticleType& particles, NeighborType& neighbor ) { _stress_timer.start(); auto model = _model; - auto neigh_list = _neigh_list; const auto x = particles.sliceReferencePosition(); const auto u = particles.sliceDisplacement(); const auto vol = particles.sliceVolume(); @@ -613,12 +641,8 @@ class Force policy( particles.frozenOffset(), - particles.localOffset() ); - Cabana::neighbor_parallel_for( - policy, stress_full, _neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "CabanaPD::ForceLinearPMB::computeStressFull" ); - + neighbor.iterate( exec_space{}, stress_full, particles, + "CabanaPD::ForceLinearPMB::computeStressFull" ); _stress_timer.stop(); } }; diff --git a/src/force_models/CabanaPD_Contact.hpp b/src/force_models/CabanaPD_Contact.hpp index be56bfff..54f7bb60 100644 --- a/src/force_models/CabanaPD_Contact.hpp +++ b/src/force_models/CabanaPD_Contact.hpp @@ -39,12 +39,18 @@ struct ContactModel ContactModel( const double _radius, const double _radius_extend ) : radius( _radius ) , radius_extend( _radius_extend ){}; + + auto cutoff() const { return 2.0 * radius + radius_extend; } + auto extend() const { return radius_extend; } + + void updateBonds( const int, const int ) {} }; /* Normal repulsion */ struct NormalRepulsionModel : public ContactModel { using base_type = ContactModel; + using model_type = ContactModel; using base_model = base_type::base_model; using fracture_type = NoFracture; using thermal_type = TemperatureIndependent; diff --git a/src/force_models/CabanaPD_Hertzian.hpp b/src/force_models/CabanaPD_Hertzian.hpp index 7f24d664..5bdc6957 100644 --- a/src/force_models/CabanaPD_Hertzian.hpp +++ b/src/force_models/CabanaPD_Hertzian.hpp @@ -23,9 +23,11 @@ namespace CabanaPD struct HertzianModel : public ContactModel { using base_type = ContactModel; + using model_type = ContactModel; using base_model = base_type::base_model; using fracture_type = NoFracture; using thermal_type = TemperatureIndependent; + using density_type = StaticDensity; using base_type::radius; double nu; // Poisson's ratio diff --git a/src/force_models/CabanaPD_LPS.hpp b/src/force_models/CabanaPD_LPS.hpp index 84dafc84..5691746c 100644 --- a/src/force_models/CabanaPD_LPS.hpp +++ b/src/force_models/CabanaPD_LPS.hpp @@ -20,15 +20,15 @@ namespace CabanaPD { template <> -struct ForceModel : public BaseForceModel +struct ForceModel { - using base_type = BaseForceModel; + using model_type = LPS; using base_model = LPS; using fracture_type = NoFracture; using thermal_type = TemperatureIndependent; + using density_type = StaticDensity; - using base_type::delta; - + double delta; int influence_type; double K; @@ -38,7 +38,7 @@ struct ForceModel : public BaseForceModel ForceModel( LPS, NoFracture, const double _delta, const double _K, const double _G, const int _influence = 0 ) - : base_type( _delta ) + : delta( _delta ) , influence_type( _influence ) , K( _K ) , G( _G ) @@ -48,7 +48,7 @@ struct ForceModel : public BaseForceModel ForceModel( LPS, Elastic, NoFracture, const double _delta, const double _K, const double _G, const int _influence = 0 ) - : base_type( _delta ) + : delta( _delta ) , influence_type( _influence ) , K( _K ) , G( _G ) @@ -56,6 +56,15 @@ struct ForceModel : public BaseForceModel init(); } + auto cutoff() const { return delta; } + + // Only needed for models which store bond properties. + void updateBonds( const int, const int ) {} + + // No-op for temperature. + KOKKOS_INLINE_FUNCTION + void thermalStretch( double&, const int, const int ) const {} + void init() { theta_coeff = 3.0 * K - 5.0 * G; @@ -76,8 +85,8 @@ struct ForceModel : public BaseForceModel return influenceFunction( xi ) * xi * xi * vol; } - KOKKOS_INLINE_FUNCTION auto dilatation( const double s, const double xi, - const double vol, + KOKKOS_INLINE_FUNCTION auto dilatation( const int, const double s, + const double xi, const double vol, const double m_i ) const { double theta_i = influenceFunction( xi ) * s * xi * xi * vol; @@ -111,6 +120,7 @@ template <> struct ForceModel : public ForceModel { + using model_type = LPS; using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = Fracture; @@ -169,6 +179,7 @@ template <> struct ForceModel : public ForceModel { + using model_type = LinearLPS; using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; @@ -192,6 +203,7 @@ template <> struct ForceModel : public ForceModel { + using model_type = LinearLPS; using base_type = ForceModel; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; diff --git a/src/force_models/CabanaPD_PMB.hpp b/src/force_models/CabanaPD_PMB.hpp index 0cb75206..2691998c 100644 --- a/src/force_models/CabanaPD_PMB.hpp +++ b/src/force_models/CabanaPD_PMB.hpp @@ -22,27 +22,26 @@ namespace CabanaPD { template <> struct ForceModel - : public BaseForceModel { - using base_type = BaseForceModel; + using model_type = PMB; using base_model = PMB; using fracture_type = NoFracture; using thermal_type = TemperatureIndependent; + using density_type = StaticDensity; - using base_type::delta; - + double delta; double c; double K; - ForceModel( PMB, NoFracture, const double delta, const double _K ) - : base_type( delta ) + ForceModel( PMB, NoFracture, const double _delta, const double _K ) + : delta( _delta ) , K( _K ) { init(); } - ForceModel( PMB, Elastic, NoFracture, const double delta, const double _K ) - : base_type( delta ) + ForceModel( PMB, Elastic, NoFracture, const double _delta, const double _K ) + : delta( _delta ) , K( _K ) { init(); @@ -50,6 +49,15 @@ struct ForceModel void init() { c = 18.0 * K / ( pi * delta * delta * delta * delta ); } + auto cutoff() const { return delta; } + + // Only needed for models which store bond properties. + void updateBonds( const int, const int ) {} + + // No-op for temperature. + KOKKOS_INLINE_FUNCTION + void thermalStretch( double&, const int, const int ) const {} + KOKKOS_INLINE_FUNCTION auto forceCoeff( const int, const int, const double s, const double vol ) const @@ -72,6 +80,7 @@ struct ForceModel : public ForceModel { using base_type = ForceModel; + using base_type::model_type; using base_model = typename base_type::base_model; using fracture_type = Fracture; using mechanics_type = Elastic; @@ -132,6 +141,7 @@ struct ForceModel; using base_plasticity_type = BasePlasticity; + using base_type::model_type; using base_model = typename base_type::base_model; using fracture_type = Fracture; using mechanics_type = ElasticPerfectlyPlastic; @@ -208,6 +218,7 @@ struct ForceModel : public ForceModel { using base_type = ForceModel; + using model_type = LinearPMB; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using mechanics_type = Elastic; @@ -282,6 +293,7 @@ struct ForceModel; using base_temperature_type = BaseTemperatureModel; + using base_type::model_type; using base_model = PMB; using fracture_type = NoFracture; using thermal_type = TemperatureDependent; @@ -315,6 +327,7 @@ struct ForceModel using base_type = ForceModel; using base_temperature_type = BaseTemperatureModel; + using base_type::model_type; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using mechanics_type = Elastic; @@ -363,10 +376,11 @@ struct ForceModel, BaseTemperatureModel { + using memory_space = typename TemperatureType::memory_space; using base_type = ForceModel; + TemperatureIndependent, memory_space>; using base_temperature_type = BaseTemperatureModel; + using typename base_type::model_type; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using mechanics_type = ElasticPerfectlyPlastic; @@ -388,10 +402,12 @@ struct ForceModel using base_type = ForceModel; using base_temperature_type = BaseDynamicTemperatureModel; + using typename base_type::model_type; using base_model = PMB; using fracture_type = NoFracture; using thermal_type = DynamicTemperature; @@ -482,6 +499,7 @@ struct ForceModel; using base_temperature_type = BaseDynamicTemperatureModel; + using typename base_type::model_type; using base_model = typename base_type::base_model; using fracture_type = typename base_type::fracture_type; using thermal_type = DynamicTemperature; @@ -525,6 +543,140 @@ ForceModel( ModelType, const double delta, const double K, const double G0, -> ForceModel; +// Force models with evolving density. +template +struct ForceDensityModel + : public ForceModel, + ForceModel +{ + using base_type = ForceModel; + using lps_base_type = + ForceModel; + using typename base_type::model_type; + using base_model = typename base_type::base_model; + using fracture_type = typename base_type::fracture_type; + using mechanics_type = typename base_type::mechanics_type; + using thermal_type = typename base_type::thermal_type; + using density_type = DynamicDensity; + + // Does not use base_type::c because it is state dependent. + using base_type::_s_p; + using base_type::bond_break_coeff; + using base_type::delta; + using base_type::G0; + using base_type::K; + using base_type::s0; + using base_type::s_Y; + double coeff; + + double rho0; + DensityType rho; + CurrentDensityType rho_current; + + // Define which base functions to use (do not use LPS). + using base_type::cutoff; + using base_type::energy; + using base_type::thermalStretch; + using base_type::updateBonds; + + // Thermal parameters + using base_type::alpha; + using base_type::temp0; + using base_type::temperature; + + ForceDensityModel( PMB model, ElasticPerfectlyPlastic mechanics, + const DensityType& _rho, + const CurrentDensityType& _rho_c, const double delta, + const double K, const double G0, const double sigma_y, + const double _rho0, const TemperatureType _temp, + const double alpha, const double temp0 = 0.0 ) + : base_type( model, mechanics, delta, K, G0, sigma_y, _temp, alpha, + temp0 ) + , lps_base_type( LPS{}, Fracture{}, delta, K, ( 3.0 / 5.0 * K ), G0 ) + , rho0( _rho0 ) + , rho( _rho ) + , rho_current( _rho_c ) + { + coeff = 18.0 / pi / delta / delta / delta / delta; + } + + KOKKOS_INLINE_FUNCTION + auto currentC( const int i ) const + { + // Initial relative density for relative density factor (HARDCODED) + double D0 = 0.69; + // Relative density: + double D = rho_current( i ) / rho0; + // Relative density factor from C. Van Nguyen et al., Journal of + // Materials Processing Technology 226 (2015): 134-145. + double RD = Kokkos::pow( ( D - D0 ) / ( 1.0 - D0 ), 1.46 * ( 1 - D0 ) ); + + // Young's modulus from C. Van Nguyen et al., Journal of Materials + // Processing Technology 226 (2015): 134-145. + // double T = temp( i ); + // double E = ( 199510 - 65.63 * T - 0.0276 * T * T - 1.754E-6 * T * T * + // T ) * RD; + // + // Bulk modulus + // double nu = 0.25; // Use bond-based model + // double K = E / ( 3 * ( 1 - 2 * nu ) ); + // + // return coeff * K; + + return coeff * K * RD; + } + + KOKKOS_INLINE_FUNCTION + auto forceCoeff( const int i, const int, const double s, + const double vol ) const + { + auto c_current = currentC( i ); + return c_current * s * vol; + } + + // Update plastic dilatation. + KOKKOS_INLINE_FUNCTION auto dilatation( const int, const double s, + const double xi, const double vol, + const double ) const + { + return coeff * s * xi * vol; + } + + // Density update using plastic dilatation. + KOKKOS_INLINE_FUNCTION void updateDensity( const int i, + const double theta_i ) const + { + // Update density using plastic dilatation. + // Note that this assumes zero initial plastic dilatation. + rho_current( i ) = Kokkos::min( rho( i ) * Kokkos::exp( -theta_i ), + rho0 ); // exp(theta_i - theta_i_0) + } + + template + void update( const ParticleType& particles ) + { + base_type::update( particles ); + rho = particles.sliceDensity(); + rho_current = particles.sliceCurrentDensity(); + } +}; + +template +ForceDensityModel( PMB, ElasticPerfectlyPlastic, DensityType rho, + const CurrentDensityType& rho_c, const double delta, + const double K, const double G0, const double sigma_y, + const double rho0, TemperatureType temp, const double _alpha, + const double _temp0 ) + -> ForceDensityModel; + } // namespace CabanaPD #endif diff --git a/unit_test/tstForce.hpp b/unit_test/tstForce.hpp index 5f26c418..a08f4398 100644 --- a/unit_test/tstForce.hpp +++ b/unit_test/tstForce.hpp @@ -772,11 +772,12 @@ void checkAnalyticalDilatation( ModelType, QuadraticTag, const double, { } -template -void initializeForce( ForceType& force, ParticleType& particles ) +template +void initializeForce( ForceType& force, ParticleType& particles, + const NeighborType& neighbor ) { - force.computeWeightedVolume( particles, Cabana::SerialOpTag{} ); - force.computeDilatation( particles, Cabana::SerialOpTag{} ); + force.computeWeightedVolume( particles, neighbor ); + force.computeDilatation( particles, neighbor ); } template @@ -808,8 +809,11 @@ void testForce( ModelType model, const double dx, const double m, // This needs to exactly match the mesh spacing to compare with the single // particle calculation. - CabanaPD::Force force( true, particles, model ); - + CabanaPD::Force + force( model ); + CabanaPD::Neighbor + neighbor( false, model, particles ); auto x = particles.sliceReferencePosition(); auto f = particles.sliceForce(); auto W = particles.sliceStrainEnergy(); @@ -817,15 +821,11 @@ void testForce( ModelType model, const double dx, const double m, auto vol = particles.sliceVolume(); // No communication needed (as in the main solver) since this test is only // intended for one rank. - initializeForce( force, particles ); - - unsigned int max_neighbors; - unsigned long long total_neighbors; - force.getNeighborStatistics( max_neighbors, total_neighbors ); + initializeForce( force, particles, neighbor ); - computeForce( force, particles, Cabana::SerialOpTag() ); - double Phi = computeEnergy( force, particles, Cabana::SerialOpTag() ); - computeStress( force, particles, Cabana::SerialOpTag() ); + computeForce( force, particles, neighbor ); + double Phi = computeEnergy( force, particles, neighbor ); + computeStress( force, particles, neighbor ); // Make a copy of final results on the host std::size_t num_particle = x.size(); @@ -857,9 +857,41 @@ void testForce( ModelType model, const double dx, const double m, boundary_width, Phi ); } +//---------------------------------------------------------------------------// +// Neighbor test function. +//---------------------------------------------------------------------------// +void testNeighbor( const double dx, const double m, const double delta ) +{ + using model_type = CabanaPD::PMB; + CabanaPD::ForceModel model( model_type{}, CabanaPD::Elastic{}, + CabanaPD::NoFracture{}, delta, 1.0 ); + + auto particles = createParticles( model_type{}, LinearTag{}, dx, 0.0 ); + + CabanaPD::Neighbor neighbor( + false, model, particles ); + + auto expected_num_neighbors = computeReferenceNeighbors( model.delta, m ); + EXPECT_EQ( expected_num_neighbors, neighbor.getMaxLocal() ); +} + //---------------------------------------------------------------------------// // GTest tests. //---------------------------------------------------------------------------// +TEST( TEST_CATEGORY, test_neighbor ) +{ + // dx needs to be decreased for increased m: boundary particles are ignored. + double m = 3; + double dx = 2.0 / 11.0; + double delta = dx * m; + testNeighbor( dx, m, delta ); + + m = 6; + dx = 2.0 / 15.0; + delta = dx * m; + testNeighbor( dx, m, delta ); +} + TEST( TEST_CATEGORY, test_force_pmb ) { // dx needs to be decreased for increased m: boundary particles are ignored.