diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index f7c8e73..8d82df8 100644 --- a/.github/workflows/cmake.yml +++ b/.github/workflows/cmake.yml @@ -37,16 +37,16 @@ jobs: working-directory: ${{github.workspace}}/build run: make UpdatePluginDataBase && ctest --test-dir example -C --rerun-failed --output-on-failure ${{env.BUILD_TYPE}} - - name: Generate a code coverage report - run: gcovr --html-details coverage.html - - - name: Generate a code coverage report - uses: threeal/gcovr-action@v1.0.0 - with: - root: ${{github.workspace}} +# - name: Generate a code coverage report +# uses: threeal/gcovr-action@v1.0.0 +# with: +# root: ${{github.workspace}} - - name: Archive code coverage results - uses: actions/upload-artifact@v4 - with: - name: code-coverage-report - path: coverage.html +# - name: Generate a code coverage report +# run: gcovr --html-details coverage.html + +# - name: Archive code coverage results +# uses: actions/upload-artifact@v4 +# with: +# name: code-coverage-report +# path: coverage.html diff --git a/example/spheres/rigid-surface/rigid_surface.msp b/example/spheres/rigid-surface/rigid_surface.msp index 80e7d13..9692085 100644 --- a/example/spheres/rigid-surface/rigid_surface.msp +++ b/example/spheres/rigid-surface/rigid_surface.msp @@ -19,28 +19,22 @@ input_data: - set_rand_vrot_arot - update_inertia -update_nbh_friction: - rcut: 1.1 m - -+first_iteration: - - chunk_neighbors_stats - compute_force: - rigid_surface: normal: [0,0,1] - offset: -1 + offset: 0.0 kt: 80000 kn: 100000 kr : 1 - mu: 0.9 + mu: 0.5 damprate: 0.999 - gravity_force - hooke_force: - config: { rcut: 1.1 m , dncut: 1.1 m, kn: 100000, kt: 100000, kr: 0.1, fc: 0.05, mu: 0.9, damp_rate: 0.9} + config: { rcut: 1.1 m , dncut: 1.1 m, kn: 100000, kt: 100000, kr: 0.1, fc: 0.05, mu: 0.1, damp_rate: 0.9} domain: - cell_size: 2 m - periodic: [false,true,false] + cell_size: 4 m + periodic: [false,false,false] write_vtklegacy: ascii: true @@ -50,9 +44,8 @@ global: simulation_dump_frequency: -1 simulation_end_iteration: 100000 simulation_log_frequency: 1000 - simulation_paraview_frequency: 5000 - dt: 0.00005 s -# dt: 0.000005 s - rcut_inc: 0.01 m + simulation_paraview_frequency: 10000 + dt: 0.0001 s + rcut_inc: 0.4 m friction_rcut: 1.1 m diff --git a/example/spheres/rigid-surface/rigid_surface_full.msp b/example/spheres/rigid-surface/rigid_surface_full.msp new file mode 100644 index 0000000..1c293b3 --- /dev/null +++ b/example/spheres/rigid-surface/rigid_surface_full.msp @@ -0,0 +1,62 @@ +grid_flavor: grid_flavor_dem + +includes: + - config_spheres.msp + +replicate_domain: + repeat: [ 5 , 5 , 1 ] + +input_data: + - read_xyz: + file: input_file_rigid_surface.xyz + bounds_mode: FILE + enlarge_bounds: 1.0 m + - set_radius: + rad: 0.5 + - set_quaternion + - set_rand_velocity: + var: 0.1 + mean: [0.0,0.0,0.0] + - set_density: + density: 0.02 + - set_rand_vrot_arot + - update_inertia + - replicate_domain + +update_nbh_friction: + rcut: 1.1 m + ++first_iteration: + - chunk_neighbors_stats + +compute_force: + - rigid_surface: + normal: [0,0,1] + offset: 0.5 + kt: 80000 + kn: 100000 + kr : 1 + mu: 0.9 + damprate: 0.999 + - gravity_force + - hooke_force: + config: { rcut: 1.1 m , dncut: 1.1 m, kn: 100000, kt: 100000, kr: 0.1, fc: 0.05, mu: 0.9, damp_rate: 0.9} + +domain: + cell_size: 2 m + periodic: [false,true,false] + +write_vtklegacy: + ascii: true + ghost: false + +global: + simulation_dump_frequency: -1 + simulation_end_iteration: 100000 + simulation_log_frequency: 1000 + simulation_paraview_frequency: 5000 + dt: 0.00005 s +# dt: 0.0001 s + rcut_inc: 0.01 m + friction_rcut: 1.1 m + diff --git a/src/forcefield/include/exaDEM/compute_wall.h b/src/forcefield/include/exaDEM/compute_wall.h index ff021a5..04c8c67 100644 --- a/src/forcefield/include/exaDEM/compute_wall.h +++ b/src/forcefield/include/exaDEM/compute_wall.h @@ -5,59 +5,86 @@ using exanb::Vec3d; -struct RigidSurfaceFunctor +namespace exaDEM { - Vec3d m_normal; - double m_offset; - double m_vel; - double m_dt; - double m_kt; - double m_kn; - double m_kr; - double m_mu; - double m_dampRate; - - ONIKA_HOST_DEVICE_FUNC inline double operator () ( - const double a_rx, const double a_ry, const double a_rz, - const double a_vx, const double a_vy, const double a_vz, - Vec3d& a_vrot, - double a_particle_radius, - double& a_fx, double& a_fy, double& a_fz, - const double a_mass, - Vec3d& a_mom, - Vec3d& a_ft) const + struct RigidSurfaceFunctor { - Vec3d pos = {a_rx,a_ry,a_rz}; - Vec3d vel = {a_vx,a_vy,a_vz}; - - const Vec3d pos_proj = dot(pos, m_normal) * m_normal; - const Vec3d contact_position = m_offset * m_normal; - - Vec3d vec_n = pos_proj - contact_position; - double n = norm(vec_n); - vec_n = vec_n / n; - const double dn = n - a_particle_radius; - Vec3d rigid_surface_center = contact_position; - const Vec3d rigid_surface_velocity = m_normal * m_vel; - constexpr Vec3d rigid_surface_angular_velocity = {0.0,0.0,0.0}; - - Vec3d f = {0.0,0.0,0.0}; - constexpr double meff = 1; - - exaDEM::hooke_force_core_v2( - dn, vec_n, - m_dt, m_kn, m_kt, m_kr, m_mu, m_dampRate, meff, - a_ft, contact_position, pos_proj, vel, f, a_mom, a_vrot, - rigid_surface_center, rigid_surface_velocity, rigid_surface_angular_velocity - ); - - // compute forces (norm) - const double res = (-1) * exanb::dot(Vec3d{f.x-a_ft.x,f.y-a_ft.y,f.z-a_ft.z}, m_normal); - - // === update forces - a_fx += f.x ; - a_fy += f.y ; - a_fz += f.z ; - return res; - } -}; + Vec3d m_normal; + double m_offset; + double m_vel; + double m_dt; + double m_kt; + double m_kn; + double m_kr; + double m_mu; + double m_dampRate; + + void test() {std::cout << "test" << std::endl;} + + ONIKA_HOST_DEVICE_FUNC inline double operator () ( + const double a_rx, const double a_ry, const double a_rz, + const double a_vx, const double a_vy, const double a_vz, + Vec3d& a_vrot, + double a_particle_radius, + double& a_fx, double& a_fy, double& a_fz, + const double a_mass, + Vec3d& a_mom, + Vec3d& a_ft) const + { + std::cout <<" la "<< std::endl; + Vec3d pos = {a_rx,a_ry,a_rz}; + Vec3d vel = {a_vx,a_vy,a_vz}; + + const Vec3d pos_proj = dot(pos, m_normal) * m_normal; + const Vec3d rigid_surface_center = m_offset * m_normal; + + Vec3d vec_n = pos_proj - rigid_surface_center; + double n = norm(vec_n); + vec_n = vec_n / n; + const double dn = n - a_particle_radius; + if (dn < 0.0) + { + const Vec3d contact_position = pos_proj - vec_n * ( a_particle_radius + 0.5 * dn ) ; + const Vec3d rigid_surface_velocity = m_normal * m_vel; + constexpr Vec3d rigid_surface_angular_velocity = {0.0,0.0,0.0}; + + Vec3d f = {0.0,0.0,0.0}; + const double meff = a_mass; + + /* + exaDEM::hooke_force_core_v2( + dn, vec_n, + m_dt, m_kn, m_kt, m_kr, m_mu, m_dampRate, meff, + a_ft, contact_position, pos_proj, vel, f, a_mom, a_vrot, + rigid_surface_center, rigid_surface_velocity, rigid_surface_angular_velocity + ); + */ + exaDEM::hooke_force_core( + dn, vec_n, + m_dt, m_kn, m_kt, m_kr, m_mu, m_dampRate, meff, + a_ft, contact_position, pos_proj, vel, f, a_mom, a_vrot, + rigid_surface_center, rigid_surface_velocity, rigid_surface_angular_velocity + ); + // compute forces (norm) + const double res = (-1) * exanb::dot(Vec3d{f.x-a_ft.x,f.y-a_ft.y,f.z-a_ft.z}, m_normal); + + // === update forces + a_fx += f.x ; + a_fy += f.y ; + a_fz += f.z ; + return res; + } + else return 0.0; + } + }; +} + + +namespace exanb +{ + template<> struct ComputeCellParticlesTraits + { + static inline constexpr bool RequiresBlockSynchronousCall = false; + static inline constexpr bool CudaCompatible = true; + }; +} diff --git a/src/forcefield/rigid_surface.cpp b/src/forcefield/rigid_surface.cpp index af146d0..2518ab2 100644 --- a/src/forcefield/rigid_surface.cpp +++ b/src/forcefield/rigid_surface.cpp @@ -34,7 +34,7 @@ namespace exaDEM ADD_SLOT( MPI_Comm , mpi , INPUT , MPI_COMM_WORLD); ADD_SLOT( GridT , grid , INPUT_OUTPUT ); ADD_SLOT( Domain , domain , INPUT , REQUIRED ); - ADD_SLOT( Vec3d , normal , INPUT , Vec3d{1.0,0.0,0.0} , DocString{"Normal vector of the rigid surface"}); + ADD_SLOT( Vec3d , normal , INPUT , Vec3d{0.0,0.0,1.0} , DocString{"Normal vector of the rigid surface"}); ADD_SLOT( double , offset , INPUT , 0.0, DocString{"Offset from the origin (0,0,0) of the rigid surface"} ); ADD_SLOT( double , dt , INPUT , REQUIRED , DocString{"Timestep of the simulation"}); ADD_SLOT( double , kt , INPUT , REQUIRED , DocString{"Parameter of the force law used to model contact cyclinder/sphere"}); @@ -56,6 +56,7 @@ namespace exaDEM // no velocity version const double vel_null = 0.; RigidSurfaceFunctor func {*normal, *offset, vel_null, *dt, *kt, *kn, *kr, *mu, *damprate}; +func.test(); compute_cell_particles( *grid , false , func , compute_field_set , parallel_execution_context() ); } };