From 63f64ab49cd2d949682a8ad11b22d4210bf56c24 Mon Sep 17 00:00:00 2001 From: Audrey Collard-Daigneault <71884806+acdaigneault@users.noreply.github.com> Date: Tue, 30 Jul 2024 14:01:42 -0400 Subject: [PATCH] Introduce the DEM action manager (#1206) Description Introducing singleton design pattern with the DEM action manager. In DEM, we often check if an action has to be performed (insertion, load balancing, contact detection, etc), and those action result in resetting containers and/or do a particle contact search. The way it was handled was turning into a spaghetti code with the addition of new features. Before, when an event was performed, a flag was set to true. And then a few places in the code was checking the flag to know if the current action has to be executed. It was getting weird because some flags where nested into others if they triggered the same actions, but the flag has lost its meaning (ex: checkpoint in load balance flag, same for the solid object mapping). In order to overcome this,, the action manager is introduced in the DEM and CFD-DEM to manage what event triggered what. How does it works? Any event that needs some following actions are handled by the action manager. The event communicates to the action manager that it happened and the action manager will change the trigger flags to the actions to do. When this action is about to be performed, it checks with the action manager if it should be trigger. Many events can triggered the same actions, but when the action is about to be performed, it only checks if it has to be performed, not what event have happened. Many checks in the DEM and CFDDEM solvers were moved into the called function, even though the feature might not be implemented. The first thing in the function is a check with a return if so. Co-authored-by: Bruno Blais --- .../dynamic_contact_search.mpirun=1.output | 180 +-- ...rt-gas-solid-fluidized-bed.mpirun=2.output | 2 +- ..._balancing_mobility_status.mpirun=2.output | 174 ++- ...oad_balancing_solid_object.mpirun=2.output | 126 +- include/dem/adaptive_sparse_contacts.h | 119 +- include/dem/dem.h | 488 +++++-- include/dem/dem_action_manager.h | 470 ++++++ include/dem/dem_contact_manager.h | 91 +- include/dem/find_contact_detection_step.h | 5 +- include/dem/grid_motion.h | 30 +- include/dem/integrator.h | 8 +- include/dem/load_balancing.h | 24 +- include/dem/periodic_boundaries_manipulator.h | 112 +- include/dem/uniform_insertion.h | 84 -- .../dem/update_local_particle_containers.h | 3 +- include/fem-dem/cfd_dem_coupling.h | 206 +-- ...nchmark_simulation_time_load_balancing.dat | 37 +- .../packing_10kparticles_dynamic.prm | 4 +- ...ng_10kparticles_dynamic_load_balancing.prm | 4 +- source/dem/CMakeLists.txt | 3 + source/dem/adaptive_sparse_contacts.cc | 33 +- source/dem/dem.cc | 982 ++++++------- source/dem/dem_action_manager.cc | 12 + source/dem/dem_contact_manager.cc | 269 ++-- source/dem/explicit_euler_integrator.cc | 24 +- source/dem/find_boundary_cells_information.cc | 5 + source/dem/find_contact_detection_step.cc | 48 +- source/dem/gear3_integrator.cc | 26 +- source/dem/grid_motion.cc | 53 +- source/dem/input_parameter_inspection.cc | 2 +- source/dem/load_balancing.cc | 38 +- source/dem/particle_wall_contact_force.cc | 2 +- source/dem/periodic_boundaries_manipulator.cc | 28 +- source/dem/read_checkpoint.cc | 4 + .../dem/update_local_particle_containers.cc | 71 +- source/dem/velocity_verlet_integrator.cc | 14 + source/fem-dem/cfd_dem_coupling.cc | 1264 ++++++++--------- tests/dem/find_contact_pairs.cc | 6 +- .../particle_particle_contact_force_linear.cc | 5 +- ...rticle_particle_contact_force_nonlinear.cc | 7 +- ...icle_particle_contact_on_two_processors.cc | 9 +- tests/dem/particle_particle_fine_search.cc | 5 +- ...wo_particles_multiple_contacts_parallel.cc | 11 +- 43 files changed, 2719 insertions(+), 2369 deletions(-) create mode 100644 include/dem/dem_action_manager.h delete mode 100644 include/dem/uniform_insertion.h create mode 100644 source/dem/dem_action_manager.cc diff --git a/applications_tests/lethe-fluid-particles/dynamic_contact_search.mpirun=1.output b/applications_tests/lethe-fluid-particles/dynamic_contact_search.mpirun=1.output index be8e138dc5..f694952468 100644 --- a/applications_tests/lethe-fluid-particles/dynamic_contact_search.mpirun=1.output +++ b/applications_tests/lethe-fluid-particles/dynamic_contact_search.mpirun=1.output @@ -95,7 +95,6 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 17 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -117,7 +116,6 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 55 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -139,7 +137,6 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 65 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -161,7 +158,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 61 +DEM contact search at dem step 92 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -183,7 +180,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 49 +DEM contact search at dem step 84 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -205,7 +202,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 31 +DEM contact search at dem step 79 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -227,8 +224,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 9 -DEM contact search at dem step 84 +DEM contact search at dem step 74 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -250,7 +246,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 57 +DEM contact search at dem step 71 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -272,8 +268,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 28 -DEM contact search at dem step 98 +DEM contact search at dem step 69 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -295,7 +290,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 66 +DEM contact search at dem step 67 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -317,7 +312,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 33 +DEM contact search at dem step 66 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -339,8 +334,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 0 -DEM contact search at dem step 66 +DEM contact search at dem step 65 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -362,8 +356,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 31 -DEM contact search at dem step 96 +DEM contact search at dem step 64 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -385,7 +378,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 61 +DEM contact search at dem step 63 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -407,8 +400,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 25 -DEM contact search at dem step 89 +DEM contact search at dem step 63 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -430,7 +422,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 53 +DEM contact search at dem step 63 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -452,8 +444,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 16 -DEM contact search at dem step 79 +DEM contact search at dem step 62 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -475,7 +466,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 42 +DEM contact search at dem step 62 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -497,8 +488,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 5 -DEM contact search at dem step 68 +DEM contact search at dem step 62 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -520,8 +510,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 31 -DEM contact search at dem step 94 +DEM contact search at dem step 62 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -543,7 +532,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 57 +DEM contact search at dem step 62 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -565,8 +554,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 20 -DEM contact search at dem step 83 +DEM contact search at dem step 62 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -588,7 +576,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 46 +DEM contact search at dem step 62 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -610,8 +598,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 9 -DEM contact search at dem step 71 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -633,8 +620,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 33 -DEM contact search at dem step 95 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -656,7 +642,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 57 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -678,8 +664,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 19 -DEM contact search at dem step 81 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -701,7 +686,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 43 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -723,8 +708,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 5 -DEM contact search at dem step 67 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -746,8 +730,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 29 -DEM contact search at dem step 91 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -769,7 +752,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 53 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -791,8 +774,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 15 -DEM contact search at dem step 77 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -814,7 +796,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 39 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -836,8 +818,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 1 -DEM contact search at dem step 63 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -859,8 +840,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 25 -DEM contact search at dem step 87 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -882,7 +862,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 49 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -904,8 +884,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 11 -DEM contact search at dem step 73 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -927,8 +906,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 35 -DEM contact search at dem step 97 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -950,7 +928,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 59 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -972,8 +950,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 21 -DEM contact search at dem step 83 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -995,7 +972,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 45 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1017,8 +994,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 7 -DEM contact search at dem step 69 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1040,8 +1016,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 31 -DEM contact search at dem step 93 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1063,7 +1038,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 55 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1085,8 +1060,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 17 -DEM contact search at dem step 79 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1108,7 +1082,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 41 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1130,8 +1104,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 3 -DEM contact search at dem step 65 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1153,8 +1126,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 27 -DEM contact search at dem step 89 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1176,7 +1148,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 51 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1198,8 +1170,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 13 -DEM contact search at dem step 75 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1221,7 +1192,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 37 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1265,8 +1236,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 23 -DEM contact search at dem step 85 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1288,7 +1258,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 47 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1310,8 +1280,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 9 -DEM contact search at dem step 71 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1333,8 +1302,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 33 -DEM contact search at dem step 95 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1356,7 +1324,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 57 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1378,8 +1346,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 19 -DEM contact search at dem step 81 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1401,7 +1368,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 43 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1423,8 +1390,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 5 -DEM contact search at dem step 67 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1446,8 +1412,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 29 -DEM contact search at dem step 91 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1469,7 +1434,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 53 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1491,8 +1456,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 15 -DEM contact search at dem step 77 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1514,7 +1478,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 39 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1536,8 +1500,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 1 -DEM contact search at dem step 63 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1559,8 +1522,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 25 -DEM contact search at dem step 87 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1582,7 +1544,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 49 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1604,8 +1566,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 11 -DEM contact search at dem step 73 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1627,8 +1588,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 35 -DEM contact search at dem step 97 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1650,7 +1610,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 59 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1672,8 +1632,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 21 -DEM contact search at dem step 83 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1695,7 +1654,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 45 +DEM contact search at dem step 61 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1717,7 +1676,7 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 7 +DEM contact search at dem step 91 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- @@ -1739,7 +1698,6 @@ id, x, y, z, v_x, v_y, v_z ---- DEM ---- -DEM contact search at dem step 36 DEM contact search at dem step 99 Finished 100 DEM iterations --------------------------------------------------------------- diff --git a/applications_tests/lethe-fluid-particles/restart-gas-solid-fluidized-bed.mpirun=2.output b/applications_tests/lethe-fluid-particles/restart-gas-solid-fluidized-bed.mpirun=2.output index 2ae1de26f4..2278af98ac 100644 --- a/applications_tests/lethe-fluid-particles/restart-gas-solid-fluidized-bed.mpirun=2.output +++ b/applications_tests/lethe-fluid-particles/restart-gas-solid-fluidized-bed.mpirun=2.output @@ -28,7 +28,7 @@ Volume-Averaged Fluid Dynamics DEM ---- DEM contact search at dem step 0 -DEM contact search at dem step 32 +DEM contact search at dem step 33 DEM contact search at dem step 49 Finished 50 DEM iterations --------------------------------------------------------------- diff --git a/applications_tests/lethe-particles/load_balancing_mobility_status.mpirun=2.output b/applications_tests/lethe-particles/load_balancing_mobility_status.mpirun=2.output index 34626655fb..ac0febde94 100644 --- a/applications_tests/lethe-particles/load_balancing_mobility_status.mpirun=2.output +++ b/applications_tests/lethe-particles/load_balancing_mobility_status.mpirun=2.output @@ -72,13 +72,13 @@ Transient iteration: 600 Time: 0.6 Time step: 0.001 ***************************************************************** | Variable | Min | Max | Average | Total | | Contact list generation | 0.0000e+00 | 1.0000e+02 | 9.9833e+01 | 5.9900e+02 | -| Velocity magnitude | 1.5208e-03 | 3.9426e-01 | 3.2227e-02 | 4.2540e+00 | -| Angular velocity magnitude | 5.1062e-02 | 1.7206e+01 | 1.5901e+00 | 2.0990e+02 | -| Translational kinetic energy | 9.6882e-08 | 6.5110e-03 | 1.9059e-04 | 2.5158e-02 | -| Rotational kinetic energy | 3.4949e-08 | 3.9684e-03 | 1.4500e-04 | 1.9140e-02 | +| Velocity magnitude | 7.1469e-04 | 3.9434e-01 | 3.1661e-02 | 4.1792e+00 | +| Angular velocity magnitude | 6.1211e-02 | 1.7212e+01 | 1.5850e+00 | 2.0922e+02 | +| Translational kinetic energy | 2.1396e-08 | 6.5137e-03 | 1.8957e-04 | 2.5023e-02 | +| Rotational kinetic energy | 5.0222e-08 | 3.9711e-03 | 1.4555e-04 | 1.9212e-02 | -->Repartitionning triangulation Load balance finished -Average, minimum and maximum number of particles on the processors are 66 , 34 and 98 +Average, minimum and maximum number of particles on the processors are 66 , 44 and 88 Minimum and maximum number of cells owned by the processors are 48 and 48 ***************************************************************** @@ -86,46 +86,42 @@ Transient iteration: 700 Time: 0.7 Time step: 0.001 ***************************************************************** | Variable | Min | Max | Average | Total | | Contact list generation | 0.0000e+00 | 1.0000e+02 | 9.9857e+01 | 6.9900e+02 | -| Velocity magnitude | 1.4196e-03 | 6.0595e-01 | 3.5198e-02 | 4.6461e+00 | -| Angular velocity magnitude | 5.9106e-02 | 1.8663e+01 | 1.4933e+00 | 1.9712e+02 | -| Translational kinetic energy | 8.4416e-08 | 1.5380e-02 | 2.8638e-04 | 3.7803e-02 | -| Rotational kinetic energy | 4.6827e-08 | 4.6690e-03 | 1.5516e-04 | 2.0481e-02 | --->Repartitionning triangulation -Load balance finished -Average, minimum and maximum number of particles on the processors are 66 , 61 and 71 -Minimum and maximum number of cells owned by the processors are 48 and 48 +| Velocity magnitude | 1.3741e-03 | 5.9858e-01 | 3.1911e-02 | 4.2123e+00 | +| Angular velocity magnitude | 4.3883e-02 | 1.8674e+01 | 1.4823e+00 | 1.9566e+02 | +| Translational kinetic energy | 7.9090e-08 | 1.5008e-02 | 2.7342e-04 | 3.6091e-02 | +| Rotational kinetic energy | 2.5812e-08 | 4.6744e-03 | 1.5938e-04 | 2.1038e-02 | ***************************************************************** Transient iteration: 800 Time: 0.8 Time step: 0.001 ***************************************************************** | Variable | Min | Max | Average | Total | | Contact list generation | 0.0000e+00 | 1.0000e+02 | 9.9875e+01 | 7.9900e+02 | -| Velocity magnitude | 1.2520e-03 | 1.6787e-01 | 2.8552e-02 | 3.7689e+00 | -| Angular velocity magnitude | 4.0283e-02 | 6.0529e+00 | 8.9381e-01 | 1.1798e+02 | -| Translational kinetic energy | 6.5663e-08 | 1.1805e-03 | 6.2609e-05 | 8.2644e-03 | -| Rotational kinetic energy | 2.1751e-08 | 4.9110e-04 | 2.5481e-05 | 3.3635e-03 | +| Velocity magnitude | 9.2879e-04 | 1.9633e-01 | 1.3990e-02 | 1.8466e+00 | +| Angular velocity magnitude | 6.5169e-02 | 6.1099e+00 | 7.6817e-01 | 1.0140e+02 | +| Translational kinetic energy | 3.6135e-08 | 1.6146e-03 | 2.4703e-05 | 3.2608e-03 | +| Rotational kinetic energy | 5.6928e-08 | 5.0038e-04 | 2.4742e-05 | 3.2659e-03 | ***************************************************************** Transient iteration: 900 Time: 0.9 Time step: 0.001 ***************************************************************** | Variable | Min | Max | Average | Total | | Contact list generation | 0.0000e+00 | 1.0000e+02 | 9.9889e+01 | 8.9900e+02 | -| Velocity magnitude | 4.1095e-04 | 2.6264e-02 | 1.0525e-02 | 1.3893e+00 | -| Angular velocity magnitude | 5.0877e-02 | 2.2496e+00 | 5.0410e-01 | 6.6541e+01 | -| Translational kinetic energy | 7.0742e-09 | 2.8894e-05 | 6.3690e-06 | 8.4071e-04 | -| Rotational kinetic energy | 3.4696e-08 | 6.7831e-05 | 6.0337e-06 | 7.9645e-04 | +| Velocity magnitude | 1.1500e-03 | 3.5821e-02 | 1.0673e-02 | 1.4089e+00 | +| Angular velocity magnitude | 8.8647e-02 | 2.5515e+00 | 4.7555e-01 | 6.2773e+01 | +| Translational kinetic energy | 5.5395e-08 | 5.3748e-05 | 6.9573e-06 | 9.1836e-04 | +| Rotational kinetic energy | 1.0533e-07 | 8.7260e-05 | 5.1835e-06 | 6.8422e-04 | ***************************************************************** Transient iteration: 1000 Time: 1 Time step: 0.001 ***************************************************************** | Variable | Min | Max | Average | Total | | Contact list generation | 0.0000e+00 | 1.0000e+02 | 9.9900e+01 | 9.9900e+02 | -| Velocity magnitude | 2.0828e-04 | 2.6264e-02 | 9.5592e-03 | 1.2618e+00 | -| Angular velocity magnitude | 5.0289e-02 | 2.2496e+00 | 4.7891e-01 | 6.3215e+01 | -| Translational kinetic energy | 1.8172e-09 | 2.8894e-05 | 5.8651e-06 | 7.7420e-04 | -| Rotational kinetic energy | 3.3899e-08 | 6.7831e-05 | 5.5538e-06 | 7.3311e-04 | +| Velocity magnitude | 2.1002e-04 | 3.5821e-02 | 9.6097e-03 | 1.2685e+00 | +| Angular velocity magnitude | 1.2169e-02 | 2.5515e+00 | 4.2320e-01 | 5.5862e+01 | +| Translational kinetic energy | 1.8477e-09 | 5.3748e-05 | 6.3275e-06 | 8.3522e-04 | +| Rotational kinetic energy | 1.9848e-09 | 8.7260e-05 | 4.5398e-06 | 5.9925e-04 | print_from_processor_0 -20 +24 [deal.II intermediate Patch<3,3>] 7 0 0 0 0.0416667 0 0 0 0.0625 0 0.0416667 0.0625 0 0 0 0.1 0.0416667 0 0.1 0 0.0625 0.1 0.0416667 0.0625 0.1 @@ -299,7 +295,7 @@ print_from_processor_0 [deal.II intermediate Patch<3,3>] 7 0.0833333 0.3125 0 0.125 0.3125 0 0.0833333 0.375 0 0.125 0.375 0 0.0833333 0.3125 0.1 0.125 0.3125 0.1 0.0833333 0.375 0.1 0.125 0.375 0.1 -16 4294967295 14 4294967295 4294967295 4294967295 +16 4294967295 14 20 4294967295 4294967295 17 1 0 1 8 @@ -309,7 +305,7 @@ print_from_processor_0 [deal.II intermediate Patch<3,3>] 7 0 0.375 0 0.0416667 0.375 0 0 0.4375 0 0.0416667 0.4375 0 0 0.375 0.1 0.0416667 0.375 0.1 0 0.4375 0.1 0.0416667 0.4375 0.1 -4294967295 19 15 4294967295 4294967295 4294967295 +4294967295 19 15 21 4294967295 4294967295 18 1 0 1 8 @@ -319,22 +315,18 @@ print_from_processor_0 [deal.II intermediate Patch<3,3>] 7 0.0416667 0.375 0 0.0833333 0.375 0 0.0416667 0.4375 0 0.0833333 0.4375 0 0.0416667 0.375 0.1 0.0833333 0.375 0.1 0.0416667 0.4375 0.1 0.0833333 0.4375 0.1 -18 4294967295 16 4294967295 4294967295 4294967295 +18 20 16 22 4294967295 4294967295 19 1 0 1 8 0 0 0 0 0 0 0 0 -0 - -print_from_processor_1 -28 [deal.II intermediate Patch<3,3>] 7 0.0833333 0.375 0 0.125 0.375 0 0.0833333 0.4375 0 0.125 0.4375 0 0.0833333 0.375 0.1 0.125 0.375 0.1 0.0833333 0.4375 0.1 0.125 0.4375 0.1 -4294967295 4294967295 4294967295 3 4294967295 4294967295 -0 1 +19 4294967295 17 23 4294967295 4294967295 +20 1 0 1 8 0 0 0 0 0 0 0 0 @@ -343,8 +335,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0 0.4375 0 0.0416667 0.4375 0 0 0.5 0 0.0416667 0.5 0 0 0.4375 0.1 0.0416667 0.4375 0.1 0 0.5 0.1 0.0416667 0.5 0.1 -4294967295 2 4294967295 4 4294967295 4294967295 -1 1 +4294967295 22 18 4294967295 4294967295 4294967295 +21 1 0 1 8 0 0 0 0 0 0 0 0 @@ -353,8 +345,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0416667 0.4375 0 0.0833333 0.4375 0 0.0416667 0.5 0 0.0833333 0.5 0 0.0416667 0.4375 0.1 0.0833333 0.4375 0.1 0.0416667 0.5 0.1 0.0833333 0.5 0.1 -1 3 4294967295 5 4294967295 4294967295 -2 1 +21 23 19 4294967295 4294967295 4294967295 +22 1 0 1 8 0 0 0 0 0 0 0 0 @@ -363,18 +355,22 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0833333 0.4375 0 0.125 0.4375 0 0.0833333 0.5 0 0.125 0.5 0 0.0833333 0.4375 0.1 0.125 0.4375 0.1 0.0833333 0.5 0.1 0.125 0.5 0.1 -2 4294967295 0 6 4294967295 4294967295 -3 1 +22 4294967295 20 4294967295 4294967295 4294967295 +23 1 0 1 8 0 0 0 0 0 0 0 0 +0 + +print_from_processor_1 +24 [deal.II intermediate Patch<3,3>] 7 0 0.5 0 0.0416667 0.5 0 0 0.5625 0 0.0416667 0.5625 0 0 0.5 0.1 0.0416667 0.5 0.1 0 0.5625 0.1 0.0416667 0.5625 0.1 -4294967295 5 1 7 4294967295 4294967295 -4 1 +4294967295 1 4294967295 3 4294967295 4294967295 +0 1 0 1 8 0 0 0 0 0 0 0 0 @@ -383,8 +379,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0416667 0.5 0 0.0833333 0.5 0 0.0416667 0.5625 0 0.0833333 0.5625 0 0.0416667 0.5 0.1 0.0833333 0.5 0.1 0.0416667 0.5625 0.1 0.0833333 0.5625 0.1 -4 6 2 8 4294967295 4294967295 -5 1 +0 2 4294967295 4 4294967295 4294967295 +1 1 0 1 8 0 0 0 0 0 0 0 0 @@ -393,8 +389,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0833333 0.5 0 0.125 0.5 0 0.0833333 0.5625 0 0.125 0.5625 0 0.0833333 0.5 0.1 0.125 0.5 0.1 0.0833333 0.5625 0.1 0.125 0.5625 0.1 -5 4294967295 3 9 4294967295 4294967295 -6 1 +1 4294967295 4294967295 5 4294967295 4294967295 +2 1 0 1 8 0 0 0 0 0 0 0 0 @@ -403,8 +399,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0 0.5625 0 0.0416667 0.5625 0 0 0.625 0 0.0416667 0.625 0 0 0.5625 0.1 0.0416667 0.5625 0.1 0 0.625 0.1 0.0416667 0.625 0.1 -4294967295 8 4 10 4294967295 4294967295 -7 1 +4294967295 4 0 6 4294967295 4294967295 +3 1 0 1 8 1 1 1 1 1 1 1 1 @@ -413,8 +409,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0416667 0.5625 0 0.0833333 0.5625 0 0.0416667 0.625 0 0.0833333 0.625 0 0.0416667 0.5625 0.1 0.0833333 0.5625 0.1 0.0416667 0.625 0.1 0.0833333 0.625 0.1 -7 9 5 11 4294967295 4294967295 -8 1 +3 5 1 7 4294967295 4294967295 +4 1 0 1 8 1 1 1 1 1 1 1 1 @@ -423,8 +419,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0833333 0.5625 0 0.125 0.5625 0 0.0833333 0.625 0 0.125 0.625 0 0.0833333 0.5625 0.1 0.125 0.5625 0.1 0.0833333 0.625 0.1 0.125 0.625 0.1 -8 4294967295 6 12 4294967295 4294967295 -9 1 +4 4294967295 2 8 4294967295 4294967295 +5 1 0 1 8 1 1 1 1 1 1 1 1 @@ -433,8 +429,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0 0.625 0 0.0416667 0.625 0 0 0.6875 0 0.0416667 0.6875 0 0 0.625 0.1 0.0416667 0.625 0.1 0 0.6875 0.1 0.0416667 0.6875 0.1 -4294967295 11 7 13 4294967295 4294967295 -10 1 +4294967295 7 3 9 4294967295 4294967295 +6 1 0 1 8 4 4 4 4 4 4 4 4 @@ -443,8 +439,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0416667 0.625 0 0.0833333 0.625 0 0.0416667 0.6875 0 0.0833333 0.6875 0 0.0416667 0.625 0.1 0.0833333 0.625 0.1 0.0416667 0.6875 0.1 0.0833333 0.6875 0.1 -10 12 8 14 4294967295 4294967295 -11 1 +6 8 4 10 4294967295 4294967295 +7 1 0 1 8 4 4 4 4 4 4 4 4 @@ -453,8 +449,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0833333 0.625 0 0.125 0.625 0 0.0833333 0.6875 0 0.125 0.6875 0 0.0833333 0.625 0.1 0.125 0.625 0.1 0.0833333 0.6875 0.1 0.125 0.6875 0.1 -11 4294967295 9 15 4294967295 4294967295 -12 1 +7 4294967295 5 11 4294967295 4294967295 +8 1 0 1 8 4 4 4 4 4 4 4 4 @@ -463,8 +459,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0 0.6875 0 0.0416667 0.6875 0 0 0.75 0 0.0416667 0.75 0 0 0.6875 0.1 0.0416667 0.6875 0.1 0 0.75 0.1 0.0416667 0.75 0.1 -4294967295 14 10 16 4294967295 4294967295 -13 1 +4294967295 10 6 12 4294967295 4294967295 +9 1 0 1 8 4 4 4 4 4 4 4 4 @@ -473,8 +469,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0416667 0.6875 0 0.0833333 0.6875 0 0.0416667 0.75 0 0.0833333 0.75 0 0.0416667 0.6875 0.1 0.0833333 0.6875 0.1 0.0416667 0.75 0.1 0.0833333 0.75 0.1 -13 15 11 17 4294967295 4294967295 -14 1 +9 11 7 13 4294967295 4294967295 +10 1 0 1 8 4 4 4 4 4 4 4 4 @@ -483,8 +479,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0833333 0.6875 0 0.125 0.6875 0 0.0833333 0.75 0 0.125 0.75 0 0.0833333 0.6875 0.1 0.125 0.6875 0.1 0.0833333 0.75 0.1 0.125 0.75 0.1 -14 4294967295 12 18 4294967295 4294967295 -15 1 +10 4294967295 8 14 4294967295 4294967295 +11 1 0 1 8 4 4 4 4 4 4 4 4 @@ -493,8 +489,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0 0.75 0 0.0416667 0.75 0 0 0.8125 0 0.0416667 0.8125 0 0 0.75 0.1 0.0416667 0.75 0.1 0 0.8125 0.1 0.0416667 0.8125 0.1 -4294967295 17 13 19 4294967295 4294967295 -16 1 +4294967295 13 9 15 4294967295 4294967295 +12 1 0 1 8 0 0 0 0 0 0 0 0 @@ -503,8 +499,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0416667 0.75 0 0.0833333 0.75 0 0.0416667 0.8125 0 0.0833333 0.8125 0 0.0416667 0.75 0.1 0.0833333 0.75 0.1 0.0416667 0.8125 0.1 0.0833333 0.8125 0.1 -16 18 14 20 4294967295 4294967295 -17 1 +12 14 10 16 4294967295 4294967295 +13 1 0 1 8 0 0 0 0 0 0 0 0 @@ -513,8 +509,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0833333 0.75 0 0.125 0.75 0 0.0833333 0.8125 0 0.125 0.8125 0 0.0833333 0.75 0.1 0.125 0.75 0.1 0.0833333 0.8125 0.1 0.125 0.8125 0.1 -17 4294967295 15 21 4294967295 4294967295 -18 1 +13 4294967295 11 17 4294967295 4294967295 +14 1 0 1 8 0 0 0 0 0 0 0 0 @@ -523,8 +519,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0 0.8125 0 0.0416667 0.8125 0 0 0.875 0 0.0416667 0.875 0 0 0.8125 0.1 0.0416667 0.8125 0.1 0 0.875 0.1 0.0416667 0.875 0.1 -4294967295 20 16 22 4294967295 4294967295 -19 1 +4294967295 16 12 18 4294967295 4294967295 +15 1 0 1 8 0 0 0 0 0 0 0 0 @@ -533,8 +529,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0416667 0.8125 0 0.0833333 0.8125 0 0.0416667 0.875 0 0.0833333 0.875 0 0.0416667 0.8125 0.1 0.0833333 0.8125 0.1 0.0416667 0.875 0.1 0.0833333 0.875 0.1 -19 21 17 23 4294967295 4294967295 -20 1 +15 17 13 19 4294967295 4294967295 +16 1 0 1 8 0 0 0 0 0 0 0 0 @@ -543,8 +539,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0833333 0.8125 0 0.125 0.8125 0 0.0833333 0.875 0 0.125 0.875 0 0.0833333 0.8125 0.1 0.125 0.8125 0.1 0.0833333 0.875 0.1 0.125 0.875 0.1 -20 4294967295 18 24 4294967295 4294967295 -21 1 +16 4294967295 14 20 4294967295 4294967295 +17 1 0 1 8 0 0 0 0 0 0 0 0 @@ -553,8 +549,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0 0.875 0 0.0416667 0.875 0 0 0.9375 0 0.0416667 0.9375 0 0 0.875 0.1 0.0416667 0.875 0.1 0 0.9375 0.1 0.0416667 0.9375 0.1 -4294967295 23 19 25 4294967295 4294967295 -22 1 +4294967295 19 15 21 4294967295 4294967295 +18 1 0 1 8 0 0 0 0 0 0 0 0 @@ -563,8 +559,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0416667 0.875 0 0.0833333 0.875 0 0.0416667 0.9375 0 0.0833333 0.9375 0 0.0416667 0.875 0.1 0.0833333 0.875 0.1 0.0416667 0.9375 0.1 0.0833333 0.9375 0.1 -22 24 20 26 4294967295 4294967295 -23 1 +18 20 16 22 4294967295 4294967295 +19 1 0 1 8 0 0 0 0 0 0 0 0 @@ -573,8 +569,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0833333 0.875 0 0.125 0.875 0 0.0833333 0.9375 0 0.125 0.9375 0 0.0833333 0.875 0.1 0.125 0.875 0.1 0.0833333 0.9375 0.1 0.125 0.9375 0.1 -23 4294967295 21 27 4294967295 4294967295 -24 1 +19 4294967295 17 23 4294967295 4294967295 +20 1 0 1 8 0 0 0 0 0 0 0 0 @@ -583,8 +579,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0 0.9375 0 0.0416667 0.9375 0 0 1 0 0.0416667 1 0 0 0.9375 0.1 0.0416667 0.9375 0.1 0 1 0.1 0.0416667 1 0.1 -4294967295 26 22 4294967295 4294967295 4294967295 -25 1 +4294967295 22 18 4294967295 4294967295 4294967295 +21 1 0 1 8 0 0 0 0 0 0 0 0 @@ -593,8 +589,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0416667 0.9375 0 0.0833333 0.9375 0 0.0416667 1 0 0.0833333 1 0 0.0416667 0.9375 0.1 0.0833333 0.9375 0.1 0.0416667 1 0.1 0.0833333 1 0.1 -25 27 23 4294967295 4294967295 4294967295 -26 1 +21 23 19 4294967295 4294967295 4294967295 +22 1 0 1 8 0 0 0 0 0 0 0 0 @@ -603,8 +599,8 @@ print_from_processor_1 [deal.II intermediate Patch<3,3>] 7 0.0833333 0.9375 0 0.125 0.9375 0 0.0833333 1 0 0.125 1 0 0.0833333 0.9375 0.1 0.125 0.9375 0.1 0.0833333 1 0.1 0.125 1 0.1 -26 4294967295 24 4294967295 4294967295 4294967295 -27 1 +22 4294967295 20 4294967295 4294967295 4294967295 +23 1 0 1 8 0 0 0 0 0 0 0 0 diff --git a/applications_tests/lethe-particles/load_balancing_solid_object.mpirun=2.output b/applications_tests/lethe-particles/load_balancing_solid_object.mpirun=2.output index 6854189f22..8a78de8704 100644 --- a/applications_tests/lethe-particles/load_balancing_solid_object.mpirun=2.output +++ b/applications_tests/lethe-particles/load_balancing_solid_object.mpirun=2.output @@ -102,9 +102,9 @@ Transient iteration: 7000 Time: 0.07 Time step: 1e-05 | Variable | Min | Max | Average | Total | | Contact list generation | 0.0000e+00 | 1.8000e+01 | 1.0000e+01 | 7.0000e+01 | | Velocity magnitude | 2.3855e-02 | 6.8660e-01 | 5.7428e-01 | 5.7428e+02 | -| Angular velocity magnitude | 0.0000e+00 | 2.6094e+02 | 3.2025e+01 | 3.2025e+04 | +| Angular velocity magnitude | 0.0000e+00 | 2.6094e+02 | 3.2027e+01 | 3.2027e+04 | | Translational kinetic energy | 8.0447e-09 | 6.6646e-06 | 4.9470e-06 | 4.9470e-03 | -| Rotational kinetic energy | 0.0000e+00 | 1.7327e-06 | 7.3872e-08 | 7.3872e-05 | +| Rotational kinetic energy | 0.0000e+00 | 1.7327e-06 | 7.3875e-08 | 7.3875e-05 | -->Repartitionning triangulation Load balance finished Average, minimum and maximum number of particles on the processors are 500 , 298 and 702 @@ -116,7 +116,7 @@ Transient iteration: 8000 Time: 0.08 Time step: 1e-05 | Variable | Min | Max | Average | Total | | Contact list generation | 0.0000e+00 | 2.1000e+01 | 1.1375e+01 | 9.1000e+01 | | Velocity magnitude | 4.4251e-02 | 7.8470e-01 | 6.3095e-01 | 6.3095e+02 | -| Angular velocity magnitude | 0.0000e+00 | 2.7191e+02 | 3.9396e+01 | 3.9396e+04 | +| Angular velocity magnitude | 0.0000e+00 | 2.7191e+02 | 3.9397e+01 | 3.9397e+04 | | Translational kinetic energy | 2.7683e-08 | 8.7051e-06 | 6.1625e-06 | 6.1625e-03 | | Rotational kinetic energy | 0.0000e+00 | 1.8814e-06 | 1.0361e-07 | 1.0361e-04 | -->Repartitionning triangulation @@ -130,9 +130,9 @@ Transient iteration: 9000 Time: 0.09 Time step: 1e-05 | Variable | Min | Max | Average | Total | | Contact list generation | 0.0000e+00 | 2.3000e+01 | 1.2667e+01 | 1.1400e+02 | | Velocity magnitude | 3.4220e-02 | 8.8280e-01 | 6.7348e-01 | 6.7348e+02 | -| Angular velocity magnitude | 0.0000e+00 | 2.9721e+02 | 4.9548e+01 | 4.9548e+04 | -| Translational kinetic energy | 1.6555e-08 | 1.1018e-05 | 7.2862e-06 | 7.2862e-03 | -| Rotational kinetic energy | 0.0000e+00 | 2.2478e-06 | 1.5431e-07 | 1.5431e-04 | +| Angular velocity magnitude | 0.0000e+00 | 2.9721e+02 | 4.9565e+01 | 4.9565e+04 | +| Translational kinetic energy | 1.6555e-08 | 1.1018e-05 | 7.2861e-06 | 7.2861e-03 | +| Rotational kinetic energy | 0.0000e+00 | 2.2478e-06 | 1.5445e-07 | 1.5445e-04 | -->Repartitionning triangulation Load balance finished Average, minimum and maximum number of particles on the processors are 500 , 277 and 723 @@ -143,10 +143,10 @@ Transient iteration: 10000 Time: 0.1 Time step: 1e-05 ***************************************************************** | Variable | Min | Max | Average | Total | | Contact list generation | 0.0000e+00 | 2.6000e+01 | 1.4000e+01 | 1.4000e+02 | -| Velocity magnitude | 1.7027e-02 | 9.7575e-01 | 7.0033e-01 | 7.0033e+02 | -| Angular velocity magnitude | 4.0468e-04 | 3.4646e+02 | 5.7410e+01 | 5.7410e+04 | -| Translational kinetic energy | 4.0987e-09 | 1.3460e-05 | 8.2239e-06 | 8.2239e-03 | -| Rotational kinetic energy | 4.1673e-18 | 3.0546e-06 | 1.9232e-07 | 1.9232e-04 | +| Velocity magnitude | 1.6797e-02 | 9.7575e-01 | 7.0037e-01 | 7.0037e+02 | +| Angular velocity magnitude | 4.0468e-04 | 3.4629e+02 | 5.7392e+01 | 5.7392e+04 | +| Translational kinetic energy | 3.9884e-09 | 1.3460e-05 | 8.2241e-06 | 8.2241e-03 | +| Rotational kinetic energy | 4.1673e-18 | 3.0515e-06 | 1.9224e-07 | 1.9224e-04 | -->Repartitionning triangulation Load balance finished Average, minimum and maximum number of particles on the processors are 500 , 246 and 754 @@ -157,10 +157,10 @@ Transient iteration: 11000 Time: 0.11 Time step: 1e-05 ***************************************************************** | Variable | Min | Max | Average | Total | | Contact list generation | 0.0000e+00 | 3.0000e+01 | 1.5455e+01 | 1.7000e+02 | -| Velocity magnitude | 2.1886e-02 | 1.0703e+00 | 7.0773e-01 | 7.0773e+02 | -| Angular velocity magnitude | 8.4763e-04 | 3.2120e+02 | 7.0106e+01 | 7.0106e+04 | -| Translational kinetic energy | 6.7720e-09 | 1.6195e-05 | 8.7622e-06 | 8.7622e-03 | -| Rotational kinetic energy | 1.8283e-17 | 2.6253e-06 | 2.5114e-07 | 2.5114e-04 | +| Velocity magnitude | 2.2057e-02 | 1.0703e+00 | 7.0757e-01 | 7.0757e+02 | +| Angular velocity magnitude | 8.4763e-04 | 3.1912e+02 | 7.0827e+01 | 7.0827e+04 | +| Translational kinetic energy | 6.8778e-09 | 1.6195e-05 | 8.7638e-06 | 8.7638e-03 | +| Rotational kinetic energy | 1.8283e-17 | 2.5915e-06 | 2.5632e-07 | 2.5632e-04 | -->Repartitionning triangulation Load balance finished Average, minimum and maximum number of particles on the processors are 500 , 275 and 725 @@ -171,13 +171,13 @@ Transient iteration: 12000 Time: 0.12 Time step: 1e-05 ***************************************************************** | Variable | Min | Max | Average | Total | | Contact list generation | 0.0000e+00 | 3.1000e+01 | 1.6750e+01 | 2.0100e+02 | -| Velocity magnitude | 2.6218e-02 | 1.1851e+00 | 6.9068e-01 | 6.9068e+02 | -| Angular velocity magnitude | 8.4763e-04 | 3.3083e+02 | 8.4993e+01 | 8.4993e+04 | -| Translational kinetic energy | 9.7178e-09 | 1.9856e-05 | 8.8361e-06 | 8.8361e-03 | -| Rotational kinetic energy | 1.8283e-17 | 2.7851e-06 | 3.2713e-07 | 3.2713e-04 | +| Velocity magnitude | 1.6864e-02 | 1.1681e+00 | 6.9464e-01 | 6.9464e+02 | +| Angular velocity magnitude | 8.4763e-04 | 4.2365e+02 | 8.2901e+01 | 8.2901e+04 | +| Translational kinetic energy | 4.0206e-09 | 1.9288e-05 | 8.8829e-06 | 8.8829e-03 | +| Rotational kinetic energy | 1.8283e-17 | 4.5672e-06 | 3.1785e-07 | 3.1785e-04 | -->Repartitionning triangulation Load balance finished -Average, minimum and maximum number of particles on the processors are 500 , 340 and 660 +Average, minimum and maximum number of particles on the processors are 500 , 339 and 661 Minimum and maximum number of cells owned by the processors are 332 and 5113 ***************************************************************** @@ -185,10 +185,10 @@ Transient iteration: 13000 Time: 0.13 Time step: 1e-05 ***************************************************************** | Variable | Min | Max | Average | Total | | Contact list generation | 0.0000e+00 | 3.4000e+01 | 1.8077e+01 | 2.3500e+02 | -| Velocity magnitude | 3.0764e-02 | 1.2852e+00 | 6.5137e-01 | 6.5137e+02 | -| Angular velocity magnitude | 8.4763e-04 | 3.7422e+02 | 9.7562e+01 | 9.7562e+04 | -| Translational kinetic energy | 1.3380e-08 | 2.3349e-05 | 8.1601e-06 | 8.1601e-03 | -| Rotational kinetic energy | 1.8283e-17 | 3.5636e-06 | 3.8522e-07 | 3.8522e-04 | +| Velocity magnitude | 2.3651e-02 | 1.2663e+00 | 6.5452e-01 | 6.5452e+02 | +| Angular velocity magnitude | 8.4763e-04 | 4.3430e+02 | 1.0014e+02 | 1.0014e+05 | +| Translational kinetic energy | 7.9081e-09 | 2.2668e-05 | 8.1847e-06 | 8.1847e-03 | +| Rotational kinetic energy | 1.8283e-17 | 4.7997e-06 | 4.0369e-07 | 4.0369e-04 | -->Repartitionning triangulation Load balance finished Average, minimum and maximum number of particles on the processors are 500 , 292 and 708 @@ -199,13 +199,13 @@ Transient iteration: 14000 Time: 0.14 Time step: 1e-05 ***************************************************************** | Variable | Min | Max | Average | Total | | Contact list generation | 0.0000e+00 | 3.7000e+01 | 1.9429e+01 | 2.7200e+02 | -| Velocity magnitude | 2.4705e-02 | 1.3679e+00 | 5.9963e-01 | 5.9963e+02 | -| Angular velocity magnitude | 8.4763e-04 | 4.2556e+02 | 1.1511e+02 | 1.1511e+05 | -| Translational kinetic energy | 8.6285e-09 | 2.6454e-05 | 6.8543e-06 | 6.8543e-03 | -| Rotational kinetic energy | 1.8283e-17 | 4.6084e-06 | 4.6909e-07 | 4.6909e-04 | +| Velocity magnitude | 2.1204e-02 | 1.3853e+00 | 5.9335e-01 | 5.9335e+02 | +| Angular velocity magnitude | 8.4763e-04 | 4.0346e+02 | 1.1614e+02 | 1.1614e+05 | +| Translational kinetic energy | 6.3563e-09 | 2.7129e-05 | 6.8245e-06 | 6.8245e-03 | +| Rotational kinetic energy | 1.8283e-17 | 4.1423e-06 | 4.8052e-07 | 4.8052e-04 | -->Repartitionning triangulation Load balance finished -Average, minimum and maximum number of particles on the processors are 500 , 288 and 712 +Average, minimum and maximum number of particles on the processors are 500 , 286 and 714 Minimum and maximum number of cells owned by the processors are 332 and 5113 ***************************************************************** @@ -213,13 +213,13 @@ Transient iteration: 15000 Time: 0.15 Time step: 1e-05 ***************************************************************** | Variable | Min | Max | Average | Total | | Contact list generation | 0.0000e+00 | 3.9000e+01 | 2.0733e+01 | 3.1100e+02 | -| Velocity magnitude | 3.9610e-02 | 1.3991e+00 | 5.2027e-01 | 5.2027e+02 | -| Angular velocity magnitude | 8.4763e-04 | 4.1235e+02 | 1.2589e+02 | 1.2589e+05 | -| Translational kinetic energy | 2.2180e-08 | 2.7673e-05 | 4.9278e-06 | 4.9278e-03 | -| Rotational kinetic energy | 1.8283e-17 | 4.3269e-06 | 5.1795e-07 | 5.1795e-04 | +| Velocity magnitude | 2.7973e-02 | 1.4042e+00 | 5.1935e-01 | 5.1935e+02 | +| Angular velocity magnitude | 8.4763e-04 | 3.6824e+02 | 1.3071e+02 | 1.3071e+05 | +| Translational kinetic energy | 1.1062e-08 | 2.7875e-05 | 4.9312e-06 | 4.9312e-03 | +| Rotational kinetic energy | 1.8283e-17 | 3.4506e-06 | 5.5743e-07 | 5.5743e-04 | -->Repartitionning triangulation Load balance finished -Average, minimum and maximum number of particles on the processors are 500 , 292 and 708 +Average, minimum and maximum number of particles on the processors are 500 , 295 and 705 Minimum and maximum number of cells owned by the processors are 332 and 5113 ***************************************************************** @@ -227,66 +227,66 @@ Transient iteration: 16000 Time: 0.16 Time step: 1e-05 ***************************************************************** | Variable | Min | Max | Average | Total | | Contact list generation | 0.0000e+00 | 3.9000e+01 | 2.1875e+01 | 3.5000e+02 | -| Velocity magnitude | 6.9846e-03 | 1.2844e+00 | 4.4613e-01 | 4.4613e+02 | -| Angular velocity magnitude | 9.9629e+00 | 4.1550e+02 | 1.2936e+02 | 1.2936e+05 | -| Translational kinetic energy | 6.8969e-10 | 2.3322e-05 | 3.3045e-06 | 3.3045e-03 | -| Rotational kinetic energy | 2.5258e-09 | 4.3931e-06 | 5.2586e-07 | 5.2586e-04 | +| Velocity magnitude | 3.4935e-02 | 1.4057e+00 | 4.4759e-01 | 4.4759e+02 | +| Angular velocity magnitude | 6.1475e+00 | 3.7461e+02 | 1.2585e+02 | 1.2585e+05 | +| Translational kinetic energy | 1.7253e-08 | 2.7936e-05 | 3.3245e-06 | 3.3245e-03 | +| Rotational kinetic energy | 9.6167e-10 | 3.5710e-06 | 5.0046e-07 | 5.0046e-04 | -->Repartitionning triangulation Load balance finished -Average, minimum and maximum number of particles on the processors are 500 , 330 and 670 +Average, minimum and maximum number of particles on the processors are 500 , 332 and 668 Minimum and maximum number of cells owned by the processors are 332 and 5113 ***************************************************************** Transient iteration: 17000 Time: 0.17 Time step: 1e-05 ***************************************************************** | Variable | Min | Max | Average | Total | -| Contact list generation | 0.0000e+00 | 3.9000e+01 | 2.2235e+01 | 3.7800e+02 | -| Velocity magnitude | 1.0104e-02 | 9.0157e-01 | 4.0160e-01 | 4.0160e+02 | -| Angular velocity magnitude | 3.2229e+00 | 3.3639e+02 | 1.1417e+02 | 1.1417e+05 | -| Translational kinetic energy | 1.4432e-09 | 1.1491e-05 | 2.7331e-06 | 2.7331e-03 | -| Rotational kinetic energy | 2.6433e-10 | 2.8795e-06 | 4.1552e-07 | 4.1552e-04 | +| Contact list generation | 0.0000e+00 | 3.9000e+01 | 2.2118e+01 | 3.7600e+02 | +| Velocity magnitude | 1.6010e-02 | 8.0070e-01 | 4.0575e-01 | 4.0575e+02 | +| Angular velocity magnitude | 8.0177e+00 | 3.7461e+02 | 1.1360e+02 | 1.1360e+05 | +| Translational kinetic energy | 3.6238e-09 | 9.0637e-06 | 2.7583e-06 | 2.7583e-03 | +| Rotational kinetic energy | 1.6358e-09 | 3.5710e-06 | 4.1614e-07 | 4.1614e-04 | -->Repartitionning triangulation Load balance finished -Average, minimum and maximum number of particles on the processors are 500 , 250 and 750 -Minimum and maximum number of cells owned by the processors are 367 and 5113 +Average, minimum and maximum number of particles on the processors are 500 , 256 and 744 +Minimum and maximum number of cells owned by the processors are 395 and 5113 ***************************************************************** Transient iteration: 18000 Time: 0.18 Time step: 1e-05 ***************************************************************** | Variable | Min | Max | Average | Total | -| Contact list generation | 0.0000e+00 | 3.9000e+01 | 2.2333e+01 | 4.0200e+02 | -| Velocity magnitude | 1.4412e-02 | 8.5655e-01 | 3.8150e-01 | 3.8150e+02 | -| Angular velocity magnitude | 7.9116e+00 | 2.9213e+02 | 1.0815e+02 | 1.0815e+05 | -| Translational kinetic energy | 2.9363e-09 | 1.0372e-05 | 2.5150e-06 | 2.5150e-03 | -| Rotational kinetic energy | 1.5928e-09 | 2.1716e-06 | 3.8217e-07 | 3.8217e-04 | +| Contact list generation | 0.0000e+00 | 3.9000e+01 | 2.2167e+01 | 3.9900e+02 | +| Velocity magnitude | 1.6299e-02 | 7.7962e-01 | 3.8385e-01 | 3.8385e+02 | +| Angular velocity magnitude | 4.8333e+00 | 2.9344e+02 | 1.0554e+02 | 1.0554e+05 | +| Translational kinetic energy | 3.7557e-09 | 8.5926e-06 | 2.5155e-06 | 2.5155e-03 | +| Rotational kinetic energy | 5.9445e-10 | 2.1911e-06 | 3.6853e-07 | 3.6853e-04 | -->Repartitionning triangulation Load balance finished -Average, minimum and maximum number of particles on the processors are 500 , 194 and 806 -Minimum and maximum number of cells owned by the processors are 409 and 5113 +Average, minimum and maximum number of particles on the processors are 500 , 317 and 683 +Minimum and maximum number of cells owned by the processors are 395 and 5113 ***************************************************************** Transient iteration: 19000 Time: 0.19 Time step: 1e-05 ***************************************************************** | Variable | Min | Max | Average | Total | -| Contact list generation | 0.0000e+00 | 3.9000e+01 | 2.2421e+01 | 4.2600e+02 | -| Velocity magnitude | 2.5351e-03 | 7.9437e-01 | 3.6784e-01 | 3.6784e+02 | -| Angular velocity magnitude | 3.7865e+00 | 3.1022e+02 | 1.0685e+02 | 1.0685e+05 | -| Translational kinetic energy | 9.0854e-11 | 8.9209e-06 | 2.3438e-06 | 2.3438e-03 | -| Rotational kinetic energy | 3.6484e-10 | 2.4489e-06 | 3.8105e-07 | 3.8105e-04 | +| Contact list generation | 0.0000e+00 | 3.9000e+01 | 2.2158e+01 | 4.2100e+02 | +| Velocity magnitude | 1.2672e-02 | 7.7803e-01 | 3.7292e-01 | 3.7292e+02 | +| Angular velocity magnitude | 8.5584e+00 | 2.9797e+02 | 1.0390e+02 | 1.0390e+05 | +| Translational kinetic energy | 2.2701e-09 | 8.5576e-06 | 2.3636e-06 | 2.3636e-03 | +| Rotational kinetic energy | 1.8639e-09 | 2.2593e-06 | 3.6722e-07 | 3.6722e-04 | -->Repartitionning triangulation Load balance finished -Average, minimum and maximum number of particles on the processors are 500 , 276 and 724 +Average, minimum and maximum number of particles on the processors are 500 , 261 and 739 Minimum and maximum number of cells owned by the processors are 409 and 5113 ***************************************************************** Transient iteration: 20000 Time: 0.2 Time step: 1e-05 ***************************************************************** | Variable | Min | Max | Average | Total | -| Contact list generation | 0.0000e+00 | 3.9000e+01 | 2.2400e+01 | 4.4800e+02 | -| Velocity magnitude | 8.1812e-03 | 8.0106e-01 | 3.5873e-01 | 3.5873e+02 | -| Angular velocity magnitude | 3.8120e+00 | 3.2775e+02 | 1.0566e+02 | 1.0566e+05 | -| Translational kinetic energy | 9.4623e-10 | 9.0717e-06 | 2.2170e-06 | 2.2170e-03 | -| Rotational kinetic energy | 3.6978e-10 | 2.7336e-06 | 3.7822e-07 | 3.7822e-04 | +| Contact list generation | 0.0000e+00 | 3.9000e+01 | 2.2150e+01 | 4.4300e+02 | +| Velocity magnitude | 9.9901e-03 | 7.8873e-01 | 3.6033e-01 | 3.6033e+02 | +| Angular velocity magnitude | 4.1493e+00 | 3.1002e+02 | 1.0293e+02 | 1.0293e+05 | +| Translational kinetic energy | 1.4109e-09 | 8.7947e-06 | 2.2137e-06 | 2.2137e-03 | +| Rotational kinetic energy | 4.3811e-10 | 2.4458e-06 | 3.6528e-07 | 3.6528e-04 | -->Repartitionning triangulation Load balance finished Average, minimum and maximum number of particles on the processors are 500 , 236 and 764 diff --git a/include/dem/adaptive_sparse_contacts.h b/include/dem/adaptive_sparse_contacts.h index fe6d32fad5..19c1b37421 100644 --- a/include/dem/adaptive_sparse_contacts.h +++ b/include/dem/adaptive_sparse_contacts.h @@ -12,7 +12,6 @@ * the top level of the Lethe distribution. * * --------------------------------------------------------------------- - * */ #ifndef lethe_adaptive_sparse_contacts_h @@ -21,6 +20,7 @@ #include #include +#include #include @@ -112,7 +112,6 @@ template class LinearAlgebra::distributed::Vector; * mobility status at nodes to check the status of the neighboring cells. * * @tparam dim Spatial dimension - * */ template class AdaptiveSparseContacts @@ -164,7 +163,6 @@ class AdaptiveSparseContacts * Without this identification of the empty cells, we cannot identify the cell * that have an empty neighbor cell, which is critical for simulations using * floating walls or mesh. - * */ enum mobility_status : unsigned int @@ -177,6 +175,34 @@ class AdaptiveSparseContacts empty = 5 // used for nodes only }; + + /** + * @brief Set the threshold values for the mobile status criteria (granular + * temperature and solid fraction) and the flag for the advection of + * particles. + * + * @param[in] granular_temperature The threshold value for the granular + * temperature. + * @param[in] solid_fraction The threshold value for the volumic solid + * fraction. + * @param[in] advect_particles The flag for the advection of particles. + */ + void + set_parameters(const double granular_temperature, + const double solid_fraction, + const double advect_particles) + { + // If the function is reached, the adaptive sparse contacts is enabled. + sparse_contacts_enabled = true; + + // Communicate to the action manager that the sparse contacts is enabled + DEMActionManager::get_action_manager()->set_sparse_contacts_enabled(); + + granular_temperature_threshold = granular_temperature; + solid_fraction_threshold = solid_fraction; + advect_particles_enabled = advect_particles; + } + /** * @brief Create or update a set of the active and ghost cells so that * there is no loop over all the cells of the triangulation for the granular @@ -187,7 +213,6 @@ class AdaptiveSparseContacts * redistributed among processors. * * @param[in] background_dh The DoFHandler of the background grid. - * */ void update_local_and_ghost_cell_set(const DoFHandler &background_dh); @@ -210,7 +235,6 @@ class AdaptiveSparseContacts * sources). * @param[in] force The contact or hydrodynamic forces. * @param[in] dt The DEM time step. - * */ void update_average_velocities_acceleration( @@ -248,7 +272,6 @@ class AdaptiveSparseContacts * particles. * @param[in] n_active_cells The number of active cells in triangulation. * @param[in] mpi_communicator The MPI communicator. - * */ void identify_mobility_status( @@ -257,47 +280,6 @@ class AdaptiveSparseContacts const unsigned int n_active_cells, MPI_Comm mpi_communicator); - /** - * @brief Identify of the mobility status of each cell through a node-based - * identification and check. Only the active and ghost cells are processed. - * For CFD-DEM, the node status have to be reset at each first DEM iteration - * of a CFD time step as mobile in order to update all forces. - * - * @param[in] background_dh The dof handler of the background grid. - * @param[in] particle_handler The particle handler that contains all the - * particles. - * @param[in] n_active_cells The number of active cells in triangulation. - * @param[in] mpi_communicator The MPI communicator. - * @param[in] counter The counter of the DEM iteration. - * - */ - inline void - identify_mobility_status( - const DoFHandler &background_dh, - const Particles::ParticleHandler &particle_handler, - const unsigned int n_active_cells, - MPI_Comm mpi_communicator, - const unsigned int counter) - { - if (counter > 0) - { - identify_mobility_status(background_dh, - particle_handler, - n_active_cells, - mpi_communicator); - } - else - { - cell_mobility_status.clear(); - - for (auto cell : local_and_ghost_cells) - { - assign_mobility_status(cell->active_cell_index(), - mobility_status::mobile); - } - } - } - /** * @brief Map the periodic nodes pairs of the triangulation using the * constraints. It allows to compare the mobility status of the nodes with the @@ -306,7 +288,6 @@ class AdaptiveSparseContacts * Note: this might not be the efficient way to pair the periodic nodes. * * @param constraints[in] The background constraints of triangulation. - * */ inline void map_periodic_nodes(const AffineConstraints &constraints) @@ -331,7 +312,6 @@ class AdaptiveSparseContacts * @brief Find the mobility status of a cell. * * @param[in] cell The iterator of the cell that needs mobility evaluation. - * */ inline unsigned int check_cell_mobility( @@ -347,7 +327,6 @@ class AdaptiveSparseContacts * index. * * @param[out] status The initiated vector for the conversion. - * */ void get_mobility_status_vector(Vector &status) @@ -358,31 +337,8 @@ class AdaptiveSparseContacts } } - /** - * @brief Set the threshold values for the mobile status criteria (granular - * temperature and solid fraction) and the flag for the advection of - * particles. - * - * @param[in] granular_temperature The threshold value for the granular - * temperature. - * @param[in] solid_fraction The threshold value for the solid fraction - * (volume of particles in the cell). - * @param[in] advect_particles The flag for the advection of particles. - * - */ - void - set_parameters(const double granular_temperature, - const double solid_fraction, - const double advect_particles) - { - granular_temperature_threshold = granular_temperature; - solid_fraction_threshold = solid_fraction; - advect_particles_enabled = advect_particles; - } - /** * @brief Give the map of the mobility status of the cells. - * */ typename DEM::dem_data_structures::cell_index_int_map & get_mobility_status() @@ -392,7 +348,6 @@ class AdaptiveSparseContacts /** * @brief Give the map of the cell-averaged velocities and accelerations * dt. - * */ std::map::active_cell_iterator, std::pair, Tensor<1, 3>>> & @@ -403,7 +358,6 @@ class AdaptiveSparseContacts /** * @brief Give the advected particles flag. - * */ bool has_advected_particles() const @@ -423,7 +377,6 @@ class AdaptiveSparseContacts * @param[out] cell_granular_temperature The empty vector of granular * temperature. * @param[out] cell_solid_fraction The empty vector of solid fraction. - * */ void calculate_granular_temperature_and_solid_fraction( @@ -445,7 +398,6 @@ class AdaptiveSparseContacts * @param[in] dof_indices The vector of the DoF indices of the cell. * @param[in] cell_status The mobility status of cell to assign. * @param[in] node_status The current mobility status of the node. - * */ inline void assign_mobility_status(unsigned int cell_id, @@ -481,7 +433,6 @@ class AdaptiveSparseContacts * * @param[in] cell_id The current cell index. * @param[in] cell_status The mobility status of cell to assign. - * */ inline void assign_mobility_status(unsigned int cell_id, const int cell_status) @@ -493,14 +444,12 @@ class AdaptiveSparseContacts * @brief Set of locally owned and ghost cells: * Used to loop over only the locally owned and ghost cells without looping * over all the cells in the triangulation numerous times. - * */ std::set::active_cell_iterator> local_and_ghost_cells; /** * @brief Map of cell mobility status: - * */ typename DEM::dem_data_structures::cell_index_int_map cell_mobility_status; @@ -509,34 +458,34 @@ class AdaptiveSparseContacts * @brief Vector of mobility status at nodes: [mobility status] * Used to check the value at node to determine the mobility status of the * cell, this type of vector is used to allow update values in parallel. - * */ LinearAlgebra::distributed::Vector mobility_at_nodes; /** * @brief Map of periodic nodes: - * */ std::unordered_map periodic_node_ids; /** * @brief Map of cell velocities and accelerations * dt: * > - * */ std::map::active_cell_iterator, std::pair, Tensor<1, 3>>> cell_velocities_accelerations; + /** + * @brief Flag for adaptive sparse contacts enabling. Useful to exit functions. + */ + bool sparse_contacts_enabled; + /** * @brief Flag for the advection of particles. - * */ bool advect_particles_enabled; /** * @brief Threshold value for granular temperature. - * */ double granular_temperature_threshold; diff --git a/include/dem/dem.h b/include/dem/dem.h index a2119a936d..b30165785a 100644 --- a/include/dem/dem.h +++ b/include/dem/dem.h @@ -24,9 +24,11 @@ #include #include +#include #include #include #include +#include #include #include #include @@ -56,8 +58,8 @@ using namespace DEM; /** - * The DEM class which initializes all the required parameters and iterates over - * the DEM iterator + * @brief Solver using the soft-sphere model of the discrete element method + * (DEM) to simulate granular systems. */ template class DEMSolver @@ -66,251 +68,465 @@ class DEMSolver DEMSolver(DEMSolverParameters dem_parameters); /** - * Initializes all the required parameters and iterates over the DEM iterator - * (DEM engine). + * @brief Engine of the DEM solver. Calls all the necessary functions to set + * parameters, solve the simulation, and finish the simulation. */ void solve(); private: /** - * @brief Manages the call to the load balance by first identifying if - * load balancing is required and then performing the load balance. + * @brief Set the parameters for the DEM simulation */ void - load_balance(); + setup_parameters(); /** - * @brief Sets the right iteration check function according to the chosen contact detection method. - * - * @return Return a function. This function returns a bool indicating if the contact search should be carried out in the current iteration. + * @brief Initialize the distribution type for the particles, and sets the + * maximum particle diameter and the neighborhood threshold squared in the + * process. */ - inline std::function - set_contact_search_iteration_function(); + void + setup_distribution_type(); /** - * @brief Establish if this is a contact detection iteration using the constant contact detection frequency. - * If the iteration number is a multiple of the frequency, this iteration is - * considered to be a contact detection iteration. - * - * @return bool indicating if the contact search should be carried out in the current iteration. + * @brief Set up the solid object containers with serial objects. + * TODO: Modify this set up once the solid volumes are fully implemented. */ - inline bool - check_contact_search_iteration_constant(); + void + setup_solid_objects(); /** - * @brief Establish if this is a contact detection iteration using the maximal displacement of the particles. - * If this particle displacement surpasses a threshold, this iteration is a - * contact detection iteration. + * @brief Set up the pointers to the functions and classes according to the + * parameters. + */ + void + setup_functions_and_pointers(); + + /** + * @brief Set the iteration check function according to the chosen + * contact detection method. * - * @return bool indicating if the contact search should be carried out in the current iteration. + * @return Contact detection method function. */ - inline bool - check_contact_search_iteration_dynamic(); + inline std::function + set_contact_search_iteration_function(); /** - * @brief Manages the call to the particle insertion. Returns true if - * particles were inserted + * @brief Set the insertion method. * + * @return The pointer to the insertion object */ - bool - insert_particles(); + std::shared_ptr> + set_insertion_type(); /** - * @brief Updates moment of inertia container after sorting particles - * into subdomains + * @brief Set the integration method. * + * @return The pointer to the integration object */ - void - update_moment_of_inertia( - dealii::Particles::ParticleHandler &particle_handler, - std::vector &MOI); + std::shared_ptr> + set_integrator_type(); + /** + * @brief Set up the triangulation dependent parameters after the reading of + * the mesh and/or the simulation restart. + */ + void + setup_triangulation_dependent_parameters(); /** - * @brief Calculates particles-wall contact forces - * + * @brief Set up the background degree of freedom used for parallel grid + * output or mapping of periodic dofs when using periodic boundaries. */ void - particle_wall_contact_force(); + setup_background_dofs(); /** - * @brief finish_simulation - * Finishes the simulation by calling all - * the post-processing elements that are required + * @brief Check if a load balancing is required according to the load + * balancing method and perform it if necessary. */ void - finish_simulation(); + load_balance(); /** - * Sets the chosen insertion method in the parameter handler file - * - * @param dem_parameters DEM parameters - * @return A pointer to the insertion object + * @brief Establish if this is a contact detection iteration using the + * constant contact detection frequency. + * If the iteration number is a multiple of the frequency, this iteration is + * considered to be a contact detection iteration. */ - std::shared_ptr> - set_insertion_type(const DEMSolverParameters &dem_parameters); + inline void + check_contact_search_iteration_constant(); /** - * Sets the chosen integration method in the parameter handler file - * - * @param dem_parameters DEM parameters - * @return A pointer to the integration object + * @brief Establish if this is a contact detection iteration using the + * maximal displacement of the particles. If this particle displacement + * surpasses a threshold, this iteration is a contact detection iteration. */ - std::shared_ptr> - set_integrator_type(const DEMSolverParameters &dem_parameters); + inline void + check_contact_search_iteration_dynamic(); /** - * Sets the chosen particle-particle contact force model in the parameter - * handler file - * - * @param dem_parameters DEM parameters - * @return A pointer to the particle-particle contact force object + * @brief Check if particles have to be inserted at this iteration and + * perform it if necessary. */ - std::shared_ptr> - set_particle_particle_contact_force( - const DEMSolverParameters &dem_parameters); + void + insert_particles(); /** - * Sets the chosen particle-wall contact force model in the parameter handler - * file - * - * @param dem_parameters DEM parameters - * @return A pointer to the particle-wall contact force object + * @brief Calculates particles-wall contact forces. */ - std::shared_ptr> - set_particle_wall_contact_force( - const DEMSolverParameters &dem_parameters); + void + particle_wall_contact_force(); /** - * Sets the background degree of freedom used for parallel grid output - * + * @brief Move all solid objects. */ void - setup_background_dofs(); + move_solid_objects(); /** - * @brief write_output_results - * Generates VTU file with particles information for visualization and - * post-processing - * Post-processing as parallel VTU files + * @brief Execute the last post-processing at the end of the simulation and + * output test results if necessary. + */ + void + finish_simulation(); + + /** + * @brief Generate VTU file with particles information for visualization. */ void write_output_results(); /** - * @brief post_process_results Calculates average velocity and other post-processed quantities that need to be outputted to files + * @brief Calculate average velocity and other post-processed quantities that + * need to be outputted to files. */ void post_process_results(); - /** - * @brief reports_statistics Calculates statistics on the particles and report them to the terminal. This function is notably used to calculate min/max/avg/total values of the linear kinetic energy, angular kinetic energy, linear velocity, angular velocity + * @brief Calculate statistics on the particles and report them to the + * terminal. This function is notably used to monitor the time min, max and + * total performed contact searches, and the instant min, max, avg and total + * values of the linear velocity, angular velocity, linear kinetic energy and + * angular kinetic. */ void report_statistics(); - MPI_Comm mpi_communicator; - const unsigned int n_mpi_processes; - const unsigned int this_mpi_process; - ConditionalOStream pcout; - DEMSolverParameters parameters; + /** + * @brief Execute the sorting of particle into subdomains and cells, and + * reinitialize the containers dependent on the local particle ids. + */ + void + sort_particles_into_subdomains_and_cells(); + + /** + * @brief The MPI communicator used for the parallel simulation. + */ + MPI_Comm mpi_communicator; + + /** + * @brief The number of MPI processes used for the parallel simulation. + */ + const unsigned int n_mpi_processes; + + /** + * @brief The rank of the current MPI process. + */ + const unsigned int this_mpi_process; + + /** + * @brief The output stream used for the parallel simulation, only used by + * the process 0. + */ + ConditionalOStream pcout; + + /** + * @brief The parameters of the DEM simulation. + */ + DEMSolverParameters parameters; + + /** + * @brief The distributed triangulation used for the DEM simulation. + */ parallel::distributed::Triangulation triangulation; - MappingQGeneric mapping; - bool particles_insertion_step; - unsigned int contact_build_number; - TimerOutput computing_timer; - double smallest_contact_search_criterion; - double smallest_floating_mesh_mapping_criterion; + /** + * @brief The polynomial mapping. The mapping is first order in DEM + * simulations. + */ + MappingQGeneric mapping; + + /** + * @brief The particle handler that manages the particles. + */ Particles::ParticleHandler particle_handler; - bool contact_detection_step; - bool load_balance_step; - bool checkpoint_step; - Tensor<1, 3> g; - - DEM::DEMProperties properties_class; - std::vector> properties = - properties_class.get_properties_name(); - double neighborhood_threshold_squared; - double maximum_particle_diameter; + + /** + * @brief The action manager that manages the actions triggered by events. + */ + DEMActionManager *action_manager; + + /** + * @brief The timer that keeps track of the time spent in some functions. + * Currently theses functions are: load balancing and VTU output. + */ + TimerOutput computing_timer; + + /** + * @brief The properties of the DEM simulation. + */ + DEM::DEMProperties properties_class; + + /** + * @brief The acceleration acting on the particles. + */ + Tensor<1, 3> g; + + /** + * @brief The contact search frequency criterion for the dynamic contact + * detection method. + * The value is the smallest value between those two quantities: + * - The diameter of the smallest active cell of a triangulation minus the + * largest particle radius. + * \f$D_{c,min} - D_{p,max} / 2\f$ + * - A security factor on the neighboring threshold times the largest particle + * radius. + * \f$factor * (neighborhood_threshold - 1) * D_{p,max} / 2\f$ + */ + double smallest_contact_search_criterion; + + /** + * @brief The smallest solid object mapping criterion. + * The value is \f$2^-0.5 * D_{c,min}\f$ + */ + double smallest_solid_object_mapping_criterion; + + /** + * @brief The neighborhood threshold squared. + * \f$(\text{neighborhood_threshold} * D_{p,max})^2\f$ + */ + double neighborhood_threshold_squared; + + /** + * @brief The maximum particle diameter. + */ + double maximum_particle_diameter; + + /** + * @brief The frequency of the contact searches. + */ const unsigned int contact_detection_frequency; + + /** + * @brief The frequency of the insertion of particles + */ const unsigned int insertion_frequency; - // Initialization of classes and building objects - DEMContactManager contact_manager; - LagrangianLoadBalancing load_balancing; + /** + * @brief The counter to keep track of the number of contact search. + */ + unsigned int contact_build_number; + + /** + * @brief The manager of all the contact search operations. + */ + DEMContactManager contact_manager; + + /** + * @brief The load balancing handler. + */ + LagrangianLoadBalancing load_balancing; + + /** + * @brief The simulation control (DEM Transient). + */ std::shared_ptr simulation_control; - BoundaryCellsInformation boundary_cell_object; - std::shared_ptr> grid_motion_object; - ParticlePointLineForce particle_point_line_contact_force_object; - std::shared_ptr> integrator_object; - std::shared_ptr> insertion_object; + + /** + * @brief The boundary cells object. + */ + BoundaryCellsInformation boundary_cell_object; + + /** + * @brief The grid motion object. + */ + std::shared_ptr> grid_motion_object; + + /** + * @brief The particle-particle contact force object. + */ std::shared_ptr> particle_particle_contact_force_object; - std::shared_ptr> particles_force_chains_object; + + /** + * @brief The particle-wall contact force object. + */ std::shared_ptr> - particle_wall_contact_force_object; - Visualization visualization_object; + particle_wall_contact_force_object; + + /** + * @brief The particle-point-line contact force object. + */ + ParticlePointLineForce particle_point_line_contact_force_object; + + /** + * @brief The particle force chains post-processing object. + */ + std::shared_ptr> particles_force_chains_object; + + /** + * @brief The post-processing object. + */ LagrangianPostProcessing post_processing_object; - PVDHandler particles_pvdhandler; - PVDHandler particles_pvdhandler_force_chains; + /** + * @brief The integrator object. + */ + std::shared_ptr> integrator_object; + + /** + * @brief The insertion object. + */ + std::shared_ptr> insertion_object; + + /** + * @brief The visualization object. + */ + Visualization visualization_object; + + /** + * @brief The particle PVD handler. + */ + PVDHandler particles_pvdhandler; + + /** + * @brief The force chains PVD handler. + */ + PVDHandler particles_pvdhandler_force_chains; + + /** + * @brief The vector of torque of particles. + */ std::vector> torque; + + /** + * @brief The vector of force of particles. + */ std::vector> force; - std::vector displacement; - std::vector MOI; - // Mesh and boundary information + /** + * @brief The displacement tracking of particles for the dynamic contact + * detection. + */ + std::vector displacement; + + /** + * @brief The moment of inertia of particles. + */ + std::vector MOI; + + /** + * @brief The vector of vector of pairs of the mapping of the background mesh + * to the solid surface. + */ typename dem_data_structures::solid_surfaces_mesh_information solid_surfaces_mesh_info; + + /** + * @brief The vector of vector of pairs of the mapping of the background mesh + * to the solid volume. + */ typename dem_data_structures::solid_volumes_mesh_info solid_volumes_mesh_info; + + /** + * @brief The map of the point and normal vector of grid. + * Only used with grid motion. + */ typename dem_data_structures::boundary_points_and_normal_vectors updated_boundary_points_and_normal_vectors; + + /** + * @brief The boundary information of the forces. + */ typename dem_data_structures::vector_on_boundary forces_boundary_information; + + /** + * @brief The torques boundary information. + */ typename dem_data_structures::vector_on_boundary torques_boundary_information; + + /** + * @brief The periodic boundaries cells information map. + * Only used with a simulation with periodic boundaries. + */ typename DEM::dem_data_structures::periodic_boundaries_cells_info periodic_boundaries_cells_information; - // Information for periodic boundaries + /** + * @brief The periodic boundaries manipulator. + * Only used with a simulation with periodic boundaries. + */ PeriodicBoundariesManipulator periodic_boundaries_object; - Tensor<1, dim> periodic_offset; - bool has_periodic_boundaries; - // Information for parallel grid processing + /** + * @brief The background degree of freedom handler uses for parallel grid + * processing. + */ DoFHandler background_dh; - PVDHandler grid_pvdhandler; - bool has_solid_objects; - // Storage of statistics about time and contact lists + /** + * @brief The grid PVD handler. + */ + PVDHandler grid_pvdhandler; + + /** + * @brief The storage of statistics about time and contact lists + */ statistics contact_list; - statistics simulation_time; - // Solid DEM objects + /** + * @brief The container of the solid surfaces. + */ std::vector>> solid_surfaces; - std::vector>> solid_volumes; - // Distribution objects + /** + * @brief The container of the solid volumes. + */ + std::vector>> solid_volumes; + + /** + * @brief The container of the distribution objects. + */ std::vector> size_distribution_object_container; - // Adaptive sparce contacts (ASC) in cells object + /** + * @brief The object handling the adaptive sparse contacts (ASC). + * Only used with the ASC method. + */ AdaptiveSparseContacts sparse_contacts_object; - // Flag to indicate if sparse contacts are enabled - bool has_sparse_contacts; - - // Contraints for the background grid needed for ASC with PBC + /** + * @brief The constraints for the background grid needed for the adaptive sparse. + */ AffineConstraints background_constraints; - // Load balancing iteration check function + /** + * @brief The pointer to the function that checks if a load balancing is + * required. + */ std::function load_balance_iteration_check_function; - // Contact detection iteration check function - std::function contact_detection_iteration_check_function; + /** + * @brief The pointer to the contact detection function that checks if a + * contact search is required. + */ + std::function contact_detection_iteration_check_function; }; #endif diff --git a/include/dem/dem_action_manager.h b/include/dem/dem_action_manager.h new file mode 100644 index 0000000000..00f65285f4 --- /dev/null +++ b/include/dem/dem_action_manager.h @@ -0,0 +1,470 @@ +#ifndef dem_action_manager_h +#define dem_action_manager_h + + +/** + * @brief Manage the actions to be performed in the DEM solver and on the DEM + * part of the coupling CFD-DEM solver. + * + * Any event that needs some following actions is handled by the action + * manager. The event communicates to the action manager that it happened and + * the action manager will change the trigger flags to the actions to be taken. + * When this action is about to be performed, it checks with the action manager + * if it should be triggered. + * Many events can trigger the same actions, but when the action is about to + * be performed, it only checks if it has to be performed, not what event have + * happened. + * + * If you implement a new feature that needs resetting or contact search after + * (or other actions), you should add a new "step" function to the action + * manager that will set the trigger to the action to be performed. + * If you implement new actions, you should add a new trigger flag and a new + * check function to the action manager. Obviously, you should also have to + * initialize the trigger in the constructor, and reset it in the reset + * function called after a DEM iteration. + * If the action depends on a enabling of a feature, you should add a new + * function to set the enabled flag. + */ +class DEMActionManager +{ +public: + /** + * @brief Copy constructor as a delete function to make sure it can not be + * copied. It will never be used. + * + * @param copy The object to be copied + */ + DEMActionManager(const DEMActionManager ©) = delete; + + /** + * @brief Copy constructor as a delete function to make sure it can not be + * assigned. It will never be used. + * + * @param copy The object to be copied + */ + DEMActionManager & + operator=(const DEMActionManager ©) = delete; + + /** + * @brief Getter of the unique instance of the DEMActionManager. + */ + static DEMActionManager * + get_action_manager(); + + /** + * @brief Reset all triggers to false. + * TODO: modify the handling of the contact search at counter = 1 because of + * the resetting of the mobility status. + */ + inline void + reset_triggers() + { + // First mobility status identification of the CFD time step (from the + // velocity computed at the first DEM time step (counter = 0) of the CFD + // time step) The contact search is executed to make sure the mobility + // status of the cell matches the particles that are in. + load_balance_trigger = false; + contact_search_trigger = mobility_status_reset_trigger; + clear_tangential_overlap_trigger = false; + solid_object_search_trigger = false; + sparse_contacts_cells_update_trigger = false; + read_checkpoint_trigger = false; + mobility_status_reset_trigger = false; + } + + /** + * @brief Flag that there are periodic boundary in the simulation. + */ + inline void + set_periodic_boundaries_enabled() + { + this->periodic_boundaries_enabled = true; + } + + /** + * @brief Check if the periodic boundaries are enabled to perform some actions + * that are not handled by triggers of the action manager. + * + * @return True if the periodic boundaries are enabled. + */ + inline bool + check_periodic_boundaries_enabled() + { + return periodic_boundaries_enabled; + } + + /** + * @brief Flag that there are solid objects in the simulation. Already trigger + * a solid object search in the background map. + */ + inline void + set_solid_objects_enabled() + { + this->solid_objects_enabled = true; + + // Allowing the first contact search of solid objects + solid_object_search_trigger = true; + } + + /** + * @brief Check if the solid objects are enabled to perform some actions + * that are not handled by triggers of the action manager. + * + * @return True if the solid objects are enabled. + */ + inline bool + check_solid_objects_enabled() + { + return solid_objects_enabled; + } + + /** + * @brief Flag that the sparse contacts is enabled and trigger a + * mobility status reset to initialize the containers. + */ + inline void + set_sparse_contacts_enabled() + { + this->sparse_contacts_enabled = true; + + // Allowing the full broad search without mobility status at first iteration + mobility_status_reset_trigger = true; + } + + /** + * @brief Check if the sparse contacts is enabled to perform some actions + * that are not handled by triggers of the action manager. + * + * @return True if the sparse contacts is enabled. + */ + inline bool + check_sparse_contacts_enabled() + { + return sparse_contacts_enabled; + } + + /** + * @brief Flag that the grid will motion in the simulation. + */ + inline void + set_grid_motion_enabled() + { + this->grid_motion_enabled = true; + } + + /** + * @brief Check if the solid objects are enabled to perform some actions + * that are not handled by triggers of the action manager. + * + * @return True if the solid objects are enabled. + */ + inline bool + check_grid_motion_enabled() + { + return grid_motion_enabled; + } + + /** + * @brief Set triggers for actions to be performed in the current time step + * because the simulation starts from restart files (previously called + * checkpoint step). + * + * It triggers the reading of checkpoints, the contact search, + * clearing the contact structures, the solid object mapping (if enabled), and + * the update of the sparse contacts cells (if enabled). + */ + inline void + restart_simulation() + { + read_checkpoint_trigger = true; + contact_search_trigger = true; + clear_tangential_overlap_trigger = true; + solid_object_search_trigger = solid_objects_enabled ? true : false; + sparse_contacts_cells_update_trigger = + sparse_contacts_enabled ? true : false; + } + + /** + * @brief Set triggers for actions to be performed in the current iteration + * because it is a load balancing step. + * + * It triggers: + * - load balancing: the used method has determined that the load needs to be + * balanced + * - contact search: cells and particles will not be on the same processor, + * so all the contact search needs to be performed. + * - clearing the tangential overlap: the tangential overlap history of the + * particles pairs (or particle-wall) are + * lost when particles are handled by + * another processor. For consistency + * reasons, all the tangential overlap + * history is reset for all particles. + * - resizing the containers: the containers need to be resized to the new + * number of local particles. + * - solid object search (if enabled): the solid objects need to be mapped to + * the rebalanced cells. + * - update the sparse contacts cells (if enabled): the sparse contacts cells + * need to be updated since they only + * contain local and ghost cells. + */ + inline void + load_balance_step() + { + load_balance_trigger = true; + contact_search_trigger = true; + clear_tangential_overlap_trigger = true; + resize_containers_trigger = true; + solid_object_search_trigger = solid_objects_enabled ? true : false; + sparse_contacts_cells_update_trigger = + sparse_contacts_enabled ? true : false; + } + + /** + * @brief Set trigger for the contact search to be performed in the current + * time step because of a particle insertion step. + * + * It triggers: + * - contact search: since there are new particles in the simulation, contact + * search needs to be performed for all particles. + */ + inline void + particle_insertion_step() + { + contact_search_trigger = true; + } + + /** + * @brief Set trigger for the contact search to be performed in the current + * time step because of a contact detection step. + * + * It triggers: + * - contact search: might be from a selected frequency for the contact search + * (frequent method) or from the displacement evaluation + * (dynamic) + */ + inline void + contact_detection_step() + { + contact_search_trigger = true; + } + + /** + * @brief Set triggers for actions to be performed in the current time step + * because of a solid object search step. + * + * It triggers the solid object search and the contact search. + */ + inline void + solid_objects_search_step() + { + solid_object_search_trigger = true; + contact_search_trigger = true; + } + + /** + * @brief Set trigger that the mobility status need to be reset to mobile at + * the first DEM time step if the sparse contacts are enabled. + * + * It trigger: + * - mobility status reset: all the mobility status of the cells need to be + * mobile since the velocity of the particles are + * not yet computed and the status can not be + * evaluated. In the case of CFD-DEM, we force the + * full computation of the particles at the first DEM + * time step since hydrodynamic forces, which have + * just been calculated, may change the behavior of + * particles. + */ + inline void + first_dem_of_cfddem_iteration_step() + { + // Only triggered if the sparse contacts are enabled + mobility_status_reset_trigger = sparse_contacts_enabled ? true : false; + } + + /** + * @brief Set trigger for the contact search to be performed at the last DEM + * iteration of the CFD iteration. + * + * It triggers: + * - contact search: the contact search is performed in order to call the + * sorting of particles in cells and for the update of the + * particle reference locations prior of the void fraction + * calculation. + */ + inline void + last_dem_of_cfddem_iteration_step() + { + contact_search_trigger = true; + } + + /** + * @brief Check if checkpoint needs to be read. + */ + inline bool + check_restart_simulation() + { + return read_checkpoint_trigger; + } + + /** + * @brief Check if the repartitioning (load balancing) needs to be performed. + */ + inline bool + check_load_balance() const + { + return load_balance_trigger; + } + + /** + * @brief Check if the contact search needs to be performed. + */ + inline bool + check_contact_search() + { + return contact_search_trigger; + } + + /** + * @brief Check if the sparse contacts cells need to be updated. + */ + inline bool + check_update_sparse_contacts_cells() + { + return sparse_contacts_cells_update_trigger; + } + + /** + * @brief Check if the solid object has to be searched. + */ + inline bool + check_solid_object_search() + { + return solid_object_search_trigger; + } + + /** + * @brief Check if the tangential overlap history needs to be cleared. + */ + inline bool + check_clear_tangential_overlap() + { + return clear_tangential_overlap_trigger; + } + + /** + * @brief Check if the containers need to be resized. + */ + inline bool + check_resize_containers() + { + return resize_containers_trigger; + } + + /** + * @brief Check if the mobility status need to be reset to mobile. + */ + inline bool + check_mobility_status_reset() + { + return mobility_status_reset_trigger; + } + + /** + * @brief Check if the default broad search functions should be used. It + * depends if adaptive sparse contacts is enabled, if so, it is a step where + * the mobility status are reset. + */ + inline bool + use_default_broad_search_functions() + { + return !sparse_contacts_enabled || mobility_status_reset_trigger; + } + +private: + /** + * @brief Constructor of the DEMActionManager. + */ + DEMActionManager() + : periodic_boundaries_enabled(false) + , solid_objects_enabled(false) + , sparse_contacts_enabled(false) + , grid_motion_enabled(false) + , read_checkpoint_trigger(false) + , load_balance_trigger(false) + , clear_tangential_overlap_trigger(false) + , resize_containers_trigger(true) + , contact_search_trigger(true) + , solid_object_search_trigger(false) + , sparse_contacts_cells_update_trigger(false) + , mobility_status_reset_trigger(false) + {} + + /** + * @brief Pointer to the unique of the DEMActionManager. + */ + static DEMActionManager *instance; + + /** + * @brief Flag for periodic boundaries in the simulation. + */ + bool periodic_boundaries_enabled; + + /** + * @brief Flag for solid objects in the simulation. + */ + bool solid_objects_enabled; + + /** + * @brief Flag for the enabling of the sparse contacts. + */ + bool sparse_contacts_enabled; + + /** + * @brief Flag for motion of grid in the simulation. + */ + bool grid_motion_enabled; + + /** + * @brief Flag of the trigger for reading the checkpoint. + */ + bool read_checkpoint_trigger; + + /** + * @brief Flag of the trigger for load balancing execution. + */ + bool load_balance_trigger; + + /** + * @brief Flag of the trigger for clearing the tangential overlap history. + */ + bool clear_tangential_overlap_trigger; + + /** + * @brief Flag of the trigger for resize the vector dependant of the local + * particles. + */ + bool resize_containers_trigger; + + /** + * @brief Flag of the trigger for the contact search. + */ + bool contact_search_trigger; + + /** + * @brief Flag of the trigger for the solid object search, mapping the solid + * with the background triangulation. + */ + bool solid_object_search_trigger; + + /** + * @brief Flag of the trigger for the update of the sparse contacts cells. + */ + bool sparse_contacts_cells_update_trigger; + + /** + * @brief Flag of the trigger for the mobility status reset to mobile status. + */ + bool mobility_status_reset_trigger; +}; +#endif diff --git a/include/dem/dem_contact_manager.h b/include/dem/dem_contact_manager.h index 0a1ffab3cd..0bfd15ae23 100644 --- a/include/dem/dem_contact_manager.h +++ b/include/dem/dem_contact_manager.h @@ -15,6 +15,9 @@ * */ +#ifndef lethe_dem_contact_manager_h +#define lethe_dem_contact_manager_h + #include #include @@ -35,9 +38,6 @@ using namespace DEM; -#ifndef lethe_dem_contact_manager_h -# define lethe_dem_contact_manager_h - /** * @brief Manage the numerous of contact detection search operations in the DEM. * @@ -72,9 +72,7 @@ class DEMContactManager execute_cell_neighbors_search( const parallel::distributed::Triangulation &triangulation, const typename DEM::dem_data_structures::periodic_boundaries_cells_info - periodic_boundaries_cells_information, - const bool has_periodic_boundaries = false, - const bool has_solid_objects = false); + periodic_boundaries_cells_information); /** * @brief Execute functions to clear and update the neighbor lists of all the @@ -93,9 +91,7 @@ class DEMContactManager update_cell_neighbors( const parallel::distributed::Triangulation &triangulation, const typename DEM::dem_data_structures::periodic_boundaries_cells_info - periodic_boundaries_cells_information, - const bool has_periodic_boundaries = false, - const bool has_solid_objects = false) + periodic_boundaries_cells_information) { cells_local_neighbor_list.clear(); cells_ghost_neighbor_list.clear(); @@ -105,9 +101,7 @@ class DEMContactManager total_neighbor_list.clear(); execute_cell_neighbors_search(triangulation, - periodic_boundaries_cells_information, - has_periodic_boundaries, - has_solid_objects); + periodic_boundaries_cells_information); } /** @@ -125,7 +119,7 @@ class DEMContactManager * containers if required. */ void - update_contacts(const bool has_periodic_boundaries = false); + update_contacts(); /** * @breif Execute functions to update the particle iterators in local-local @@ -143,27 +137,7 @@ class DEMContactManager */ void update_local_particles_in_cells( - const Particles::ParticleHandler &particle_handler, - const bool clear_contact_structures, - const bool has_periodic_boundaries = false); - - /** - * @brief Execute the particle-particle broad searches. - * - * Calls proper functions to find the candidates of local and ghost - * particle-particle contact pairs and the periodic particle-particle contacts - * if required. These contact pairs will be used in the fine search step to - * investigate if they are in contact. - * - * @param[in] particle_handler Storage of particles and their accessor - * functions. search[in] - * @param has_periodic_boundaries Allow manipulations of periodic containers - * if required. - */ - void - execute_particle_particle_broad_search( - dealii::Particles::ParticleHandler &particle_handler, - const bool has_periodic_boundaries = false); + const Particles::ParticleHandler &particle_handler); /** * @brief Execute the particle-particle broad searches with adaptive sparse @@ -172,8 +146,9 @@ class DEMContactManager * Calls proper functions to find the candidates of local and ghost * particle-particle contact pairs and the periodic particle-particle contacts * if required. These contact pairs will be used in the fine search step to - * investigate if they are in contact. This version of the function is used - * when adaptive sparse contacts regards mobility is enabled. + * investigate if they are in contact. + * It checks if the adaptive sparse contacts is enabled and use proper + * functions. * * @param[in,out] particle_handler Storage of particles and their accessor * functions. @@ -185,35 +160,7 @@ class DEMContactManager void execute_particle_particle_broad_search( dealii::Particles::ParticleHandler &particle_handler, - const AdaptiveSparseContacts &sparse_particle_contact_object, - const bool has_periodic_boundaries = false); - - /** - * @brief Execute the particle-wall broad searches. - * - * Calls proper functions to find the candidates of particle-wall contact - * pairs and the particle-floating wall or particle-floating mesh if required. - * These contact pairs will be used in the fine search step to investigate if - * they are in contact. - * - * @param[in] particle_handler Storage of particles and their accessor - * functions. - * @param[in] boundary_cells_object Information of the boundary cells and - * faces. - * @param[in] solid_surfaces_mesh_info Mapping of solid surfaces meshes. - * @param[in] floating_wall Properties of the floating walls. - * @param[in] simulation_time Current simulation time. - * @param[in] has_solid_objects Allow dealing with floating mesh neighbors. - */ - void - execute_particle_wall_broad_search( - const Particles::ParticleHandler &particle_handler, - BoundaryCellsInformation &boundary_cell_object, - const typename dem_data_structures::solid_surfaces_mesh_information - solid_surfaces_mesh_info, - const Parameters::Lagrangian::FloatingWalls &floating_walls, - const double simulation_time, - const bool has_solid_objects = false); + const AdaptiveSparseContacts &sparse_particle_contact_object); /** * @brief Execute the particle-wall broad searches with adaptive sparse @@ -222,8 +169,9 @@ class DEMContactManager * Calls proper functions to find the candidates of particle-wall contact * pairs and the particle-floating wall or particle-floating mesh if required. * These contact pairs will be used in the fine search step to investigate if - * they are in contact. This version of the function is used when adaptive - * sparse contacts regards mobility is enabled. + * they are in contact. + * It checks if the adaptive sparse contacts is enabled and use proper + * functions. * * @param[in] particle_handler Storage of particles and their accessor * functions. @@ -244,8 +192,7 @@ class DEMContactManager solid_surfaces_mesh_info, const Parameters::Lagrangian::FloatingWalls &floating_walls, const double simulation_time, - const AdaptiveSparseContacts &sparse_particle_contact_object, - const bool has_solid_objects = false); + const AdaptiveSparseContacts &sparse_particle_contact_object); /** * @brief Execute the particle-particles fine searches. @@ -262,8 +209,7 @@ class DEMContactManager void execute_particle_particle_fine_search( const double neighborhood_threshold, - const bool has_periodic_boundaries = false, - const Tensor<1, dim> periodic_offset = Tensor<1, dim>()); + const Tensor<1, dim> periodic_offset = Tensor<1, dim>()); /** * @brief Execute the particle-wall fine searches. @@ -280,8 +226,7 @@ class DEMContactManager execute_particle_wall_fine_search( const Parameters::Lagrangian::FloatingWalls &floating_walls, const double simulation_time, - const double neighborhood_threshold, - const bool has_solid_objects = false); + const double neighborhood_threshold); // Container with the iterators to all local and ghost particles diff --git a/include/dem/find_contact_detection_step.h b/include/dem/find_contact_detection_step.h index 2830fa3c12..3574ac1b8c 100644 --- a/include/dem/find_contact_detection_step.h +++ b/include/dem/find_contact_detection_step.h @@ -50,13 +50,12 @@ using namespace dealii; */ template -bool +void find_particle_contact_detection_step( Particles::ParticleHandler &particle_handler, const double dt, const double smallest_contact_search_criterion, MPI_Comm &mpi_communicator, - const bool sorting_in_subdomains_step, std::vector &displacement, const bool parallel_update = true); @@ -72,7 +71,7 @@ find_particle_contact_detection_step( */ template -bool +void find_floating_mesh_mapping_step( const double smallest_contact_search_criterion, std::vector>> solids); diff --git a/include/dem/grid_motion.h b/include/dem/grid_motion.h index 1c89d1fde9..73239d442f 100644 --- a/include/dem/grid_motion.h +++ b/include/dem/grid_motion.h @@ -15,6 +15,9 @@ * */ +#ifndef grid_motion_h +#define grid_motion_h + #include #include #include @@ -23,16 +26,13 @@ using namespace dealii; -#ifndef grid_motion_h -# define grid_motion_h - /** - * Carries out motion of the triangulation, including rotational, translational + * @brief Carries out motion of the triangulation, including rotational, translational * and rotational-translational motions. * * @note Grid motion with a ParsedFunction should be added to the code as a future * improvement - * @param triangulation Triagulation + * @param triangulation Triangulation * @param dem_parameters Input DEM parameters in the parameter handler file * */ @@ -46,7 +46,7 @@ class GridMotion public: /** - * The constructor sets up the grid motion function defined in the parameter + * @brief The constructor sets up the grid motion function defined in the parameter * handler, and calculates the rotation angle for rotational and * translational-rotational motions. * @@ -58,7 +58,7 @@ class GridMotion const double dem_time_step); /** - * Calls the desired grid motion. + * @brief Calls the desired grid motion. * * @param triangulation Triangulation */ @@ -69,7 +69,7 @@ class GridMotion } /** - * Carries out updating the boundary points and normal vectors in the + * @brief Carries out updating the boundary points and normal vectors in the * particle-wall contact list. * * @param particle_wall_pairs_in_contact The particle-wall contact list container. @@ -91,7 +91,7 @@ class GridMotion private: /** - * Carries out rotational motion of the triangulation + * @brief Carries out rotational motion of the triangulation * * @param triangulation Triangulation */ @@ -106,6 +106,18 @@ class GridMotion void move_grid_translational(Triangulation &triangulation); + /** + * @brief Dummy function for no motion for grid motion pointer to point to. + * + * @param triangulation Triangulation + */ + void + no_motion(Triangulation & /* triangulation */) + { + return; + } + + // Since the DEM time-step and rotational speed are constant, we calculate the // rotation angle at each time-step once in the constructor and define it as a // member variable. diff --git a/include/dem/integrator.h b/include/dem/integrator.h index 98a544dba2..40d61d5975 100644 --- a/include/dem/integrator.h +++ b/include/dem/integrator.h @@ -15,17 +15,19 @@ * */ + +#ifndef integrator_h +#define integrator_h + #include #include +#include #include #include using namespace dealii; -#ifndef integration_h -# define integration_h - /** * Base interface for classes that carry out the integration of the velocity and * position of particles with inertia diff --git a/include/dem/load_balancing.h b/include/dem/load_balancing.h index 4de72023a4..cdaf31e4e7 100644 --- a/include/dem/load_balancing.h +++ b/include/dem/load_balancing.h @@ -106,7 +106,7 @@ class LagrangianLoadBalancing * * @return bool indicating if this is a load balance iteration. */ - bool + inline void check_load_balance_iteration() { return iteration_check_function(); @@ -121,7 +121,7 @@ class LagrangianLoadBalancing * @return Return a function that returns a bool indicating if the current * time step is a load balance iteration. */ - inline std::function + inline std::function set_iteration_check_function() { using namespace Parameters::Lagrangian; @@ -129,15 +129,15 @@ class LagrangianLoadBalancing switch (load_balance_method) { case ModelParameters::LoadBalanceMethod::once: - return [&] { return check_load_balance_once(); }; + return [&] { check_load_balance_once(); }; case ModelParameters::LoadBalanceMethod::frequent: - return [&] { return check_load_balance_frequent(); }; + return [&] { check_load_balance_frequent(); }; case ModelParameters::LoadBalanceMethod::dynamic: - return [&] { return check_load_balance_dynamic(); }; + return [&] { check_load_balance_dynamic(); }; case ModelParameters::LoadBalanceMethod::dynamic_with_sparse_contacts: - return [&] { return check_load_balance_with_sparse_contacts(); }; + return [&] { check_load_balance_with_sparse_contacts(); }; default: // Default is no load balance (none) - return [&]() { return false; }; + return [&]() { return; }; } } @@ -148,7 +148,7 @@ class LagrangianLoadBalancing * @return bool indicating if the current time step is a load balance * iteration according to the load balancing `step`. */ - bool + void check_load_balance_once(); /** @@ -158,7 +158,7 @@ class LagrangianLoadBalancing * @return bool indicating if the current time step is a load balance * iteration according to the load balancing `frequency`. */ - bool + void check_load_balance_frequent(); /** @@ -169,7 +169,7 @@ class LagrangianLoadBalancing * iteration according to the load balancing `dynamic check frequency`and * `threshold`. */ - bool + void check_load_balance_dynamic(); /** @@ -182,7 +182,7 @@ class LagrangianLoadBalancing * iteration according to the load balancing `dynamic check frequency`, * `threshold` and load factors. */ - bool + void check_load_balance_with_sparse_contacts(); /** @@ -367,7 +367,7 @@ class LagrangianLoadBalancing * @brief The load balancing iteration check function according to the load * balancing method. */ - std::function iteration_check_function; + std::function iteration_check_function; /** * @brief Default hard-coded load weight of a cell. diff --git a/include/dem/periodic_boundaries_manipulator.h b/include/dem/periodic_boundaries_manipulator.h index 79505f0d32..e62dded81f 100644 --- a/include/dem/periodic_boundaries_manipulator.h +++ b/include/dem/periodic_boundaries_manipulator.h @@ -15,8 +15,13 @@ * */ + +#ifndef periodic_boundaries_manipulator_h +#define periodic_boundaries_manipulator_h + #include #include +#include #include #include @@ -31,9 +36,6 @@ using namespace dealii; -#ifndef periodic_boundaries_manipulator_h -# define periodic_boundaries_manipulator_h - /** * This class corresponds to a manipulator of the particles crossing periodic * cells in DEM. It maps cell information of the pairs of periodic cells @@ -54,30 +56,33 @@ class PeriodicBoundariesManipulator PeriodicBoundariesManipulator(); /** - * @brief Sets the periodic boundaries parameters + * @brief Sets the periodic boundaries parameters. * - * @param periodic_boundary_id_0 Id of the first periodic boundary - * @param periodic_boundary_id_1 Id of the second periodic boundary - * @param periodic_direction Perpendicular axis of PB + * @param[in] periodic_boundary_id_0 Id of the first periodic boundary. + * @param[in] periodic_direction Perpendicular axis of PB. */ void set_periodic_boundaries_information( const types::boundary_id periodic_boundary_id_0, - const types::boundary_id periodic_boundary_id_1, const unsigned int periodic_direction) { + // If function is reached, there are periodic boundaries in the simulation + periodic_boundaries_enabled = true; + + // Communicate to the action manager that there are periodic boundaries + DEMActionManager::get_action_manager()->set_periodic_boundaries_enabled(); + periodic_boundary_0 = periodic_boundary_id_0; - periodic_boundary_1 = periodic_boundary_id_1; direction = periodic_direction; } /** * @brief Execute the mapping of the cells on periodic boundaries and store - * information in periodic_boundaries_cells_information + * information in periodic_boundaries_cells_information. * - * @param triangulation Triangulation of mesh - * @param periodic_boundaries_cells_information Map of information of the - * pair of cells on periodic boundaries + * @param[in] triangulation Triangulation. + * @param[out] periodic_boundaries_cells_information Map of information of the + * pair of cells on periodic boundaries. */ void map_periodic_cells( @@ -87,24 +92,29 @@ class PeriodicBoundariesManipulator /** * @brief Moves particles crossing periodic boundaries (any side) - * particle_handle doesn't allow automated particle displacement since it is - * not linked to triangulation and its periodic mapping + * particle_handler doesn't allow automated particle displacement since it is + * not linked to triangulation and its periodic mapping. + * + * @param[in,out] particle_handler Particle handler of particles located in + * boundary cells. + * @param[in] periodic_boundaries_cells_information Map of information of the + * pair of cells on periodic boundaries. * - * @param particle_handler Particle handler of particles located in boundary - * cells - * @param periodic_boundaries_cells_information Map of information of the - * pair of cells on periodic boundaries + * @return Flag if at least one particle has been moved to the other periodic + * cell. */ - void + bool execute_particles_displacement( const Particles::ParticleHandler &particle_handler, const typename DEM::dem_data_structures::periodic_boundaries_cells_info &periodic_boundaries_cells_information); /** - * @brief Return the periodic offset distance, it is calculated from the first pair of - * cells on periodic boundaries, all pair of cells are assumed to have the - * same offset + * @brief Return the periodic offset distance, it is calculated from the first + * pair of cells on periodic boundaries, all pair of cells are assumed to + * have the same offset. + * + * @return Offset distance between periodic boundaries. */ inline Tensor<1, dim> get_periodic_offset_distance() @@ -116,12 +126,12 @@ class PeriodicBoundariesManipulator /** * @brief Gets boundaries information related to the face at periodic * boundaries 0 and 1 and stores in periodic_boundaries_cells_info_struct - * object + * object. * - * @param cell Current cell on boundary - * @param face_id Face located on boundary - * @param boundaries_information Reference to the object with periodic - * boundaries information + * @param[in] cell Current cell on boundary. + * @param[in] face_id Face located on boundary. + * @param[out] boundaries_information Reference to the object with periodic + * boundaries information. */ void get_periodic_boundaries_info( @@ -132,27 +142,49 @@ class PeriodicBoundariesManipulator /** * @brief Checks if particle is outside of the cell, if so, modifies the * location of the particle with the distance (offset) between the periodic - * faces + * faces. * - * @param boundaries_cells_content Reference to the object with periodic - * boundaries information - * @param particles_in_cell Iterator to the particles in cell - * @param particles_in_pb0_cell If the particles are linked to a cell on the - * periodic boundary 0 (true) or the periodic boundary 1 (false) + * @param[in] boundaries_cells_content Reference to the object with periodic + * boundaries information. + * @param[in] particles_in_pb0_cell If the particles are linked to a cell on + * the periodic boundary 0 (true) or the periodic boundary 1 (false). + * @param[in,out] particles_in_cell Iterator to the particles in cell. + * @param[out] particle_has_been_moved Flag the particle has been moved to + * the other periodic cell. */ void check_and_move_particles( const periodic_boundaries_cells_info_struct &boundaries_cells_content, + const bool &particles_in_pb0_cell, typename Particles::ParticleHandler::particle_iterator_range &particles_in_cell, - bool &particles_in_pb0_cell); + bool &particle_has_been_moved); + /** + * @brief Flag for periodic boundary conditions in simulation. Useful to + * exit function when there are no periodic boundaries. + */ + bool periodic_boundaries_enabled; + + /** + * @brief Id of the first periodic boundary. No needs to store the second one + * since there are linked on the triangulation, and accessible through + * functions on cells on the boundary condition 0. + */ types::boundary_id periodic_boundary_0; - types::boundary_id periodic_boundary_1; - unsigned int direction; - Tensor<1, dim> constant_periodic_offset; -}; + /** + * @brief Direction of the periodic boundary, it is the perpendicular axis of + * the periodic boundaries. + */ + unsigned int direction; + /** + * @brief Offset distance between periodic boundaries, it is calculated from + * the first pair of cells on periodic boundaries, all pair of cells are + * assumed to have the same offset. + */ + Tensor<1, dim> constant_periodic_offset; +}; -#endif /* periodic_boundaries_manipulator_h */ +#endif diff --git a/include/dem/uniform_insertion.h b/include/dem/uniform_insertion.h deleted file mode 100644 index c43dff61ed..0000000000 --- a/include/dem/uniform_insertion.h +++ /dev/null @@ -1,84 +0,0 @@ -/* --------------------------------------------------------------------- - * - * Copyright (C) 2019 - 2024 by the Lethe authors - * - * This file is part of the Lethe library - * - * The Lethe library is free software; you can use it, redistribute - * it, and/or modify it under the terms of the GNU Lesser General - * Public License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * The full text of the license can be found in the file LICENSE at - * the top level of the Lethe distribution. - * - * --------------------------------------------------------------------- - * - */ - -#include - -#include -#include - -#include - -#include - -#ifndef uniform_insertion_h -# define uniform_insertion_h - -/** - * Uniform insertion of particles in a rectangular box - * - * @note - * - * @author Shahab Golshan, Bruno Blais, Polytechnique Montreal 2019- - */ - -template -class UniformInsertion : public Insertion -{ -public: - UniformInsertion(const DEMSolverParameters &dem_parameters, - const double maximum_particle_diameter); - - /** - * Carries out the insertion of particles by discretizing and looping over the - * insertion box, finding the initial position of particles in the insertion - * box, their id and other initial properties of the particles. This virtual - * function that serves as a template for all insertion functions. - * - * @param particle_handler The particle handler of particles which are being - * inserted - * @param triangulation Triangulation to access the cells in which the - * particles are inserted - * @param dem_parameters DEM parameters declared in the .prm file - */ - virtual void - insert(Particles::ParticleHandler &particle_handler, - const parallel::distributed::Triangulation &triangulation, - const DEMSolverParameters &dem_parameters) override; - -private: - /** - * Converts id of particles to uniform insertion location - * - * @param insertion_location Insertion location of the particle - * @param id Particle_id - * @param insertion_information DEM insertion parameters declared in the .prm - * file - */ - void - find_insertion_location_uniform( - Point &insertion_location, - const unsigned int &id, - const Parameters::Lagrangian::InsertionInfo &insertion_information); - - // Number of remained particles of each type that should be inserted in the - // upcoming insertion steps - unsigned int particles_of_each_type_remaining; - - unsigned int current_inserting_particle_type; -}; - -#endif /* uniform_insertion_h */ diff --git a/include/dem/update_local_particle_containers.h b/include/dem/update_local_particle_containers.h index cac6e1a668..5c036cf17d 100644 --- a/include/dem/update_local_particle_containers.h +++ b/include/dem/update_local_particle_containers.h @@ -62,8 +62,7 @@ void update_contact_container_iterators( pairs_structure &pairs_in_contact, const typename DEM::dem_data_structures::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures = false); + &particle_container); #endif /* locate_local_particles_h */ diff --git a/include/fem-dem/cfd_dem_coupling.h b/include/fem-dem/cfd_dem_coupling.h index 7f04ceba7b..8a12cf5095 100644 --- a/include/fem-dem/cfd_dem_coupling.h +++ b/include/fem-dem/cfd_dem_coupling.h @@ -23,6 +23,8 @@ #include #include +#include +#include #include #include #include @@ -45,221 +47,169 @@ using namespace dealii; +/** + * @brief Solver using the GLS Volume-averaged Navier-Stokes (VANS) solver for + * fluid and the soft-sphere model of the discrete element method (DEM) to + * simulate solid-fluid flow with two-way coupling. + */ template class CFDDEMSolver : public GLSVANSSolver { - using FuncPtrType = bool (CFDDEMSolver::*)(const unsigned int &counter); - FuncPtrType check_contact_search_step; - public: CFDDEMSolver(CFDDEMSimulationParameters &nsparam); ~CFDDEMSolver(); - virtual void - solve() override; - /** - * @brief Manages the call to the load balancing. Returns true if - * load balancing is performed - * + * @brief Engine of the CFD-DEM solver. Calls all the necessary functions to + * set parameters, solve the simulation, and finish the simulation. */ - void - load_balance(); + virtual void + solve() override; -protected: private: /** - * @brief Carries out the DEM calculations in the DEM_CFD solver. Particle-particle and particle-wall contact force calculations, integration and update_ghost + * @brief Set up the various parameters related to the DEM. */ void - dem_iterator(unsigned int counter); + dem_setup_parameters(); /** - * @brief Carries out the particle-particle and particle-wall contact searches, sort_particles_into_subdomains_and_cells and exchange_ghost + * @brief Initialize the distribution type for the particles, and sets the + * maximum particle diameter and the neighborhood threshold squared in the + * process. */ void - dem_contact_build(unsigned int counter); + setup_distribution_type(); /** - * @brief Sets up the various parameters related to the DEM contacts + * @brief Set the integration method. + * + * @return The pointer to the integration object */ - void - dem_setup_contact_parameters(); + std::shared_ptr> + set_integrator_type(); /** - * @brief Carries out the initialization of DEM parameters + * @brief Initialize some DEM parameters. */ void initialize_dem_parameters(); /** - * @brief write DEM_output_results - * Post-processing as parallel VTU files + * @brief Read the DEM restart files from a DEM simulation. */ void - write_DEM_output_results(); + read_dem(); /** - * @brief Carries out the fine particled-wall contact detection - * + * @brief Write the CFD-DEM restart files. */ void - particle_wall_fine_search(); + write_checkpoint() override; /** - * @brief Calculates particles-wall contact forces - * + * @brief Read the CFD-DEM restart files. */ void - particle_wall_contact_force(); + read_checkpoint() override; /** - * @brief Updates moment of inertia container after sorting particles - * into subdomains + * @brief Execute the contact detection method. * + * @param[in] counter The DEM iteration. */ void - update_moment_of_inertia( - dealii::Particles::ParticleHandler &particle_handler, - std::vector &MOI); + check_contact_detection_method(unsigned int counter); /** - * Sets the chosen integration method in the parameter handler file - * - * @return A pointer to the integration object + * @brief Check if a load balancing is required according to the load + * balancing method and perform it if necessary. */ - std::shared_ptr> - set_integrator_type(); + void + load_balance(); /** - * Adds fluid-particle interaction force to the "force" container - * + * @brief Add fluid-particle interaction force to the "force" container. */ void add_fluid_particle_interaction_force(); /** - * Adds fluid-particle interaction torque to the "torque" container - * + * @brief Add fluid-particle interaction torque to the "torque" container. */ void add_fluid_particle_interaction_torque(); /** - * Sets the chosen particle-particle contact force model in the parameter - * handler file - * - * @return A pointer to the particle-particle contact force object + * @brief Calculate particles-wall contact forces. */ - std::shared_ptr> - set_particle_particle_contact_force(); + void + particle_wall_contact_force(); /** - * Sets the chosen particle-wall contact force model in the parameter handler - * file - * - * @return A pointer to the particle-wall contact force object + * @brief Post-processing as parallel VTU files. */ - std::shared_ptr> - set_particle_wall_contact_force(); - void - read_dem(); + write_dem_output_results(); + /** + * @brief Calculate statistics on the particles and report them to the + * terminal. This function is notably used to monitor the time min, max and + * total performed contact searches, and the instant min, max, avg and total + * values of the linear velocity, angular velocity, linear kinetic energy and + * angular kinetic. + */ void - write_checkpoint() override; + report_particle_statistics(); + /** + * @brief Print the final summary of the particles. + */ void - read_checkpoint() override; + print_particles_summary(); + /** + * @brief Post-process fluid dynamics after an iteration. + */ void - print_particles_summary(); + postprocess_fd(bool first_iteration) override; /** - * @brief dem_post_process_results + * @brief Post-process cfd-dem after an iteration. */ void - dem_post_process_results(); + postprocess_cfd_dem(); /** - * @brief postprocess - * Post-process fluid dynamics after an iteration + * @brief Dynamic flow control calculation that take into account the void + * fraction for the average velocity calculation. */ void - postprocess_fd(bool first_iteration) override; + dynamic_flow_control() override; /** - * @brief postprocess for cfd-dem - * Post-process cfd-dem after an iteration + * @brief Execute the sorting of particle into subdomains and cells, and + * reinitialize the containers dependent on the local particle ids. */ void - postprocess_cfd_dem(); + sort_particles_into_subdomains_and_cells(); /** - * @brief dynamic_flow_control - * Dynamic flow control calculation that take into account the void fraction - * for the average velocity calculation + * @brief Execute a complete DEM iteration, including the particle-particle + * and particle-wall contact force calculations, integration. */ void - dynamic_flow_control() override; + dem_iterator(unsigned int counter); /** - * @brief Checks all the conditions that require a contact search step. The - * check of conditions is done in order of suspected frequency occurrence. - * - * @param counter The counter of DEM iterations in a CFD iteration. + * @brief Check if a contact search has to be performed and execute it if so. */ - inline bool - contact_search_step(const unsigned int counter) - { - if (contact_detection_step) - { - // Contact search step according to the contact detection method - return true; - } - else if (counter == (coupling_frequency - 1)) - { - // Execute the contact search at the last DEM iteration of the CFD - // iteration. It updates the reference location of particles with - // the last calculated locations prior the void fraction calculation. - return true; - } - else if (has_sparse_contacts && (counter == 1)) - { - // First mobility status identification of the CFD time step (from the - // velocity computed at the first DEM time step (counter = 0) of the CFD - // time step) The contact search is executed to make sure the mobility - // status of cell match the particles that are in. - return true; - } - else if (load_balance_step) - { - // Needs to update contacts since particles/cells may have been - // distributed to a different subdomain - return true; - } - else if ((this->simulation_control->is_at_start() && (counter == 0)) || - checkpoint_step) - { - // First contact search of the simulation - return true; - } - else if (has_periodic_boundaries && particle_displaced_in_pbc) - { - // Particles have been displaced in periodic boundaries - return true; - } - else - { - return false; - } - } + void + dem_contact_build(unsigned int counter); + unsigned int coupling_frequency; - bool contact_detection_step; - bool checkpoint_step; - bool load_balance_step; Tensor<1, 3> g; std::vector> torque; std::vector> force; @@ -297,15 +247,10 @@ class CFDDEMSolver : public GLSVANSSolver // Information for periodic boundaries PeriodicBoundariesManipulator periodic_boundaries_object; Tensor<1, dim> periodic_offset; - bool has_periodic_boundaries; - bool particle_displaced_in_pbc; // Object handling the sparse contacts AdaptiveSparseContacts sparse_contacts_object; - // Flag to indicate if sparse contacts are used - bool has_sparse_contacts; - // Counter of contact searches in a CFD iteration unsigned int contact_search_counter; @@ -314,7 +259,6 @@ class CFDDEMSolver : public GLSVANSSolver // Storage of statistics about time and contact lists statistics contact_list; - statistics simulation_time; DEM::DEMProperties properties_class; @@ -331,5 +275,7 @@ class CFDDEMSolver : public GLSVANSSolver /// Post-processing variables to output total fluid volume and total particles /// volume TableHandler table_phase_volumes; + + DEMActionManager *dem_action_manager; }; #endif diff --git a/performance_analyses/dem/short_packing_10k_particles/benchmark_simulation_time_load_balancing.dat b/performance_analyses/dem/short_packing_10k_particles/benchmark_simulation_time_load_balancing.dat index aaef0ff424..b4070bd247 100755 --- a/performance_analyses/dem/short_packing_10k_particles/benchmark_simulation_time_load_balancing.dat +++ b/performance_analyses/dem/short_packing_10k_particles/benchmark_simulation_time_load_balancing.dat @@ -1,33 +1,32 @@ Processor: AMD Ryzen 9 3950X 16-Core Processor _____________________________________________________________________________________________________________________________________________________________________________________ -|Branch | SHA | date | time_1p (s) | time_4p (s) | time_8p (s) | time_12p (s) | time_16p (s) -|master | 402d5805085080ef0034141a60c466217c48916f | 24/06/2022 | 77.5 | 42.8 | 28.8 | 24.2 | 22 -|emplace_back_fine_search | 1ddb2e31dcf5ff44a14273a8f282020eadc24307 | 11/12/2022 | 55.8 | 35.1 | 25.2 | 21.0 | 20 -|master | c51887b58f365ab571a9abd20215d123d032f34c | 03/03/2024 | 53.6 | 29.8 | 22.6 | 21.0 | 21.2 +|Branch | SHA | date | time_1p (s) | time_4p (s) | time_8p (s) | time_12p (s) | time_16p (s) +|master | 402d5805085080ef0034141a60c466217c48916f | 24/06/2022 | 77.5 | 42.8 | 28.8 | 24.2 | 22 +|emplace_back_fine_search | 1ddb2e31dcf5ff44a14273a8f282020eadc24307 | 11/12/2022 | 55.8 | 35.1 | 25.2 | 21.0 | 20 +|master | c51887b58f365ab571a9abd20215d123d032f34c | 03/03/2024 | 53.6 | 29.8 | 22.6 | 21.0 | 21.2 Processor: 12th Gen Intel(R) Core(TM) i9-12900K ___________________________________________________________________________________________________________________________________________________________________________________________ -|Branch | SHA | date | time_1p (s) | time_4p (s) | time_8p (s) | time_12p (s) | time_16p (s) -|master | 2212c29fcea13a3aeb2daae9dbc190ca83ed9340 | 18/10/2022 | 50.4 | 26.8 | 16.63 | 20.3 | 19.5 -|dem_optimize_find_cell_neighbors | abdb2998346093680fa7df8c507b9ed47ab93c8e | 08/12/2022 | 48.3 | 25.0 | 17.2 | 19.3 | 17.4 -|dem_refactor_containers_for_contact_properties | 4daeeaa4179fe3302f3109bf76361a3e6dc16ada | 08/12/2022 | 48.8 | 25.4 | 15.3 | 16.7 | 15.5 -|master | d3b06b889e54356dd1addcf2cb7db867566f96ff | 16/12/2022 | 52.6 | 23.1 | 14.6 | 18.1 | 16.4 -|dem_bcs_optimize | 927ca632a9f8d6e73b29e4e2bf5be0695b2bd4f4 | 21/12/2022 | 49.3 | 21.3 | 13.7 | 17.3 | 15.1 -|master | 0652b8c75da540931f3ec035ae4fd417bf8064ae | 30/03/2023 | 45.3 | 19.3 | 11.5 | 15.1 | 14.1 -|master | 58e9f4f1ddf7287954687d196a98cf0b597e79b1 | 24/10/2023 | 41.8 | 21.4 | 16.0 #comment : MPI Issues on Manjaro -|dem_contact_storage | 58e9f4f1ddf7287954687d196a98cf0b597e79b1 | 24/10/2023 | 41.1 | 21 | 14.3 #comment : MPI Issues on Manjaro +|Branch | SHA | date | time_1p (s) | time_4p (s) | time_8p (s) | time_12p (s) | time_16p (s) +|master | 2212c29fcea13a3aeb2daae9dbc190ca83ed9340 | 18/10/2022 | 50.4 | 26.8 | 16.63 | 20.3 | 19.5 +|dem_optimize_find_cell_neighbors | abdb2998346093680fa7df8c507b9ed47ab93c8e | 08/12/2022 | 48.3 | 25.0 | 17.2 | 19.3 | 17.4 +|dem_refactor_containers_for_contact_properties | 4daeeaa4179fe3302f3109bf76361a3e6dc16ada | 08/12/2022 | 48.8 | 25.4 | 15.3 | 16.7 | 15.5 +|master | d3b06b889e54356dd1addcf2cb7db867566f96ff | 16/12/2022 | 52.6 | 23.1 | 14.6 | 18.1 | 16.4 +|dem_bcs_optimize | 927ca632a9f8d6e73b29e4e2bf5be0695b2bd4f4 | 21/12/2022 | 49.3 | 21.3 | 13.7 | 17.3 | 15.1 +|master | 0652b8c75da540931f3ec035ae4fd417bf8064ae | 30/03/2023 | 45.3 | 19.3 | 11.5 | 15.1 | 14.1 +|master | 58e9f4f1ddf7287954687d196a98cf0b597e79b1 | 24/10/2023 | 41.8 | 21.4 | 16.0 #comment : MPI Issues on Manjaro +|dem_contact_storage | 58e9f4f1ddf7287954687d196a98cf0b597e79b1 | 24/10/2023 | 41.1 | 21 | 14.3 #comment : MPI Issues on Manjaro Processor: AMD Ryzen 9 7950X 16-Core Processor ___________________________________________________________________________________________________________________________________________________________________________________________ -|Branch | SHA | date | time_1p (s) | time_4p (s) | time_8p (s) | time_12p (s) | time_16p (s) -|master | 2212c29fcea13a3aeb2daae9dbc190ca83ed9340 | 18/10/2022 | 31.3 | 17.6 | 12.5 | 11.5 | 11.2 -|fix_DEM_memory_growth | d01db9579180768ccf51d5af5cbdb73749961127 | 27/02/2024 | 32.3 | 17.5 | 12.5 | 11.0 | 11.2 - - - +|Branch | SHA | date | time_1p (s) | time_4p (s) | time_8p (s) | time_12p (s) | time_16p (s) +|master | 2212c29fcea13a3aeb2daae9dbc190ca83ed9340 | 18/10/2022 | 31.3 | 17.6 | 12.5 | 11.5 | 11.2 +|fix_DEM_memory_growth | d01db9579180768ccf51d5af5cbdb73749961127 | 27/02/2024 | 32.3 | 17.5 | 12.5 | 11.0 | 11.2 +|master | a091baf2a75c0128e09edbdc07a174278b85f9bc | 30/07/2024 | 31.2 | 14.9 | 9.74 | 9.4 | 8.8 +|master | 22ef4d7a08ac50f033e479adca462598f335440c | 30/07/2024 | 31.2 | 14.9 | 9.8 | 9.4 | 8.6 diff --git a/performance_analyses/dem/short_packing_10k_particles/packing_10kparticles_dynamic.prm b/performance_analyses/dem/short_packing_10k_particles/packing_10kparticles_dynamic.prm index a67fe7cf7a..2621a310e6 100644 --- a/performance_analyses/dem/short_packing_10k_particles/packing_10kparticles_dynamic.prm +++ b/performance_analyses/dem/short_packing_10k_particles/packing_10kparticles_dynamic.prm @@ -52,9 +52,7 @@ end #--------------------------------------------------- subsection lagrangian physical properties - set gx = 0.0 - set gy = 0.0 - set gz = -9.81 + set g = 0.0, 0, -9.81 set number of particle types = 1 subsection particle type 0 set size distribution type = uniform diff --git a/performance_analyses/dem/short_packing_10k_particles/packing_10kparticles_dynamic_load_balancing.prm b/performance_analyses/dem/short_packing_10k_particles/packing_10kparticles_dynamic_load_balancing.prm index 7149a54d4a..cb41c8405a 100644 --- a/performance_analyses/dem/short_packing_10k_particles/packing_10kparticles_dynamic_load_balancing.prm +++ b/performance_analyses/dem/short_packing_10k_particles/packing_10kparticles_dynamic_load_balancing.prm @@ -56,9 +56,7 @@ end #--------------------------------------------------- subsection lagrangian physical properties - set gx = 0.0 - set gy = 0.0 - set gz = -9.81 + set g = 0.0, 0, -9.81 set number of particle types = 1 subsection particle type 0 set size distribution type = uniform diff --git a/source/dem/CMakeLists.txt b/source/dem/CMakeLists.txt index ad55e0caf9..7ceacee583 100644 --- a/source/dem/CMakeLists.txt +++ b/source/dem/CMakeLists.txt @@ -3,6 +3,7 @@ add_library(lethe-dem adaptive_sparse_contacts.cc data_containers.cc dem.cc + dem_action_manager.cc dem_contact_manager.cc dem_solver_parameters.cc distributions.cc @@ -53,6 +54,8 @@ add_library(lethe-dem ../../include/dem/contact_type.h ../../include/dem/data_containers.h ../../include/dem/dem.h + ../../include/dem/dem_action_manager.h + ../../include/dem/dem_contact_manager.h ../../include/dem/dem_solver_parameters.h ../../include/dem/distributions.h ../../include/dem/explicit_euler_integrator.h diff --git a/source/dem/adaptive_sparse_contacts.cc b/source/dem/adaptive_sparse_contacts.cc index 569b7f170a..805629926e 100644 --- a/source/dem/adaptive_sparse_contacts.cc +++ b/source/dem/adaptive_sparse_contacts.cc @@ -17,6 +17,7 @@ */ #include +#include #include @@ -24,6 +25,8 @@ template AdaptiveSparseContacts::AdaptiveSparseContacts() + : sparse_contacts_enabled(false) + , advect_particles_enabled(false) {} template @@ -31,13 +34,14 @@ void AdaptiveSparseContacts::update_local_and_ghost_cell_set( const DoFHandler &background_dh) { + if (!sparse_contacts_enabled) + return; + local_and_ghost_cells.clear(); for (const auto &cell : background_dh.active_cell_iterators()) { if (cell->is_locally_owned() || cell->is_ghost()) - { - local_and_ghost_cells.insert(cell); - } + local_and_ghost_cells.insert(cell); } } @@ -141,9 +145,25 @@ AdaptiveSparseContacts::identify_mobility_status( const unsigned int n_active_cells, MPI_Comm mpi_communicator) { + // If the sparse contacts are not enabled, exit the function + if (!sparse_contacts_enabled) + return; + // Reset cell status containers cell_mobility_status.clear(); + // Reset all status to mobile if mobility status needs to be reset + if (DEMActionManager::get_action_manager()->check_mobility_status_reset()) + { + for (auto cell : local_and_ghost_cells) + { + assign_mobility_status(cell->active_cell_index(), + mobility_status::mobile); + } + + return; + } + // Get a copy of the active & ghost cells set to iterate over and remove cell // of the set when the mobility status is known. It avoids unnecessary // iteration for next loops. We don't want to modify the original set since @@ -358,6 +378,13 @@ AdaptiveSparseContacts::update_average_velocities_acceleration( const std::vector> &force, const double dt) { + // If the sparse contacts is disabled, exit the function. + // Also, if the mobility status do not need to be reset, which only happens + // at first dem iteration of the cfd iteration + if (!sparse_contacts_enabled || + !DEMActionManager::get_action_manager()->check_mobility_status_reset()) + return; + cell_velocities_accelerations.clear(); // Tensors for velocity and acceleration * dt computation diff --git a/source/dem/dem.cc b/source/dem/dem.cc index 3f2e89b6d0..8b6f392068 100644 --- a/source/dem/dem.cc +++ b/source/dem/dem.cc @@ -6,7 +6,6 @@ #include #include #include -#include #include #include #include @@ -40,32 +39,28 @@ DEMSolver::DEMSolver(DEMSolverParameters dem_parameters) , parameters(dem_parameters) , triangulation(this->mpi_communicator) , mapping(1) - , particles_insertion_step(0) - , contact_build_number(0) + , particle_handler(triangulation, mapping, DEM::get_number_properties()) , computing_timer(this->mpi_communicator, this->pcout, TimerOutput::summary, TimerOutput::wall_times) - , particle_handler(triangulation, mapping, DEM::get_number_properties()) - , contact_detection_step(true) - , load_balance_step(true) - , checkpoint_step(true) , contact_detection_frequency( parameters.model_parameters.contact_detection_frequency) , insertion_frequency(parameters.insertion_info.insertion_frequency) - , has_periodic_boundaries(false) + , contact_build_number(0) , background_dh(triangulation) - , has_solid_objects(false) , size_distribution_object_container( parameters.lagrangian_physical_properties.particle_type_number) - , has_sparse_contacts(false) -{ - // Print simulation starting information +{} + +template +void +DEMSolver::setup_parameters() +{ // Print simulation starting information + pcout << std::endl; std::stringstream ss; - ss << "Running on " << n_mpi_processes << " rank(s)"; - announce_string(pcout, ss.str(), '*'); // Check if the output directory exists @@ -81,75 +76,109 @@ DEMSolver::DEMSolver(DEMSolverParameters dem_parameters) } } + // Get the pointer of the only instance of the action manager + action_manager = DEMActionManager::get_action_manager(); + // Change the behavior of the timer for situations when you don't want outputs if (parameters.timer.type == Parameters::Timer::Type::none) computing_timer.disable_output(); + + // Set the simulation control as transient DEM simulation_control = std::make_shared( parameters.simulation_control); - // Setup load balancing parameters + // Setup load balancing parameters and attach the correct functions to the + // signals inside the triangulation load_balancing.set_parameters(parameters.model_parameters); load_balancing.copy_references(simulation_control, triangulation, particle_handler, sparse_contacts_object); - - // Attach the correct functions to the signals inside the triangulation load_balancing.connect_weight_signals(); - + // Set the adaptive sparse contacts parameters if (parameters.model_parameters.sparse_particle_contacts) { - has_sparse_contacts = true; sparse_contacts_object.set_parameters( parameters.model_parameters.granular_temperature_threshold, parameters.model_parameters.solid_fraction_threshold, parameters.model_parameters.advect_particles); } - maximum_particle_diameter = 0; - for (unsigned int particle_type = 0; - particle_type < - parameters.lagrangian_physical_properties.particle_type_number; - particle_type++) + // Set the distribution type and initialize the neighborhood threshold + setup_distribution_type(); + + if (this_mpi_process == 0) + input_parameter_inspection(parameters, + pcout, + size_distribution_object_container); + + // Set the grid motion type + grid_motion_object = + std::make_shared>(parameters.grid_motion, + simulation_control->get_time_step()); + + // Set up the solid objects + setup_solid_objects(); + + // Check if there are periodic boundaries + for (unsigned int i_bc = 0; + i_bc < parameters.boundary_conditions.bc_types.size(); + ++i_bc) { - if (parameters.lagrangian_physical_properties.distribution_type.at( - particle_type) == - Parameters::Lagrangian::SizeDistributionType::uniform) - { - size_distribution_object_container[particle_type] = - std::make_shared( - parameters.lagrangian_physical_properties - .particle_average_diameter.at(particle_type)); - } - else if (parameters.lagrangian_physical_properties.distribution_type.at( - particle_type) == - Parameters::Lagrangian::SizeDistributionType::normal) + if (parameters.boundary_conditions.bc_types[i_bc] == + Parameters::Lagrangian::BCDEM::BoundaryType::periodic) { - size_distribution_object_container[particle_type] = - std::make_shared( - parameters.lagrangian_physical_properties - .particle_average_diameter.at(particle_type), - parameters.lagrangian_physical_properties.particle_size_std.at( - particle_type), - parameters.lagrangian_physical_properties - .seed_for_distributions[particle_type] + - this_mpi_process); + periodic_boundaries_object.set_periodic_boundaries_information( + parameters.boundary_conditions.periodic_boundary_0, + parameters.boundary_conditions.periodic_direction); + break; } - else if (parameters.lagrangian_physical_properties.distribution_type.at( - particle_type) == - Parameters::Lagrangian::SizeDistributionType::custom) + } + + // Assign gravity/acceleration + g = parameters.lagrangian_physical_properties.g; + + // If this is a start simulation + if (parameters.restart.restart) + action_manager->restart_simulation(); +} + +template +void +DEMSolver::setup_distribution_type() +{ + // Use namespace and alias to make the code more readable + using namespace Parameters::Lagrangian; + LagrangianPhysicalProperties &lpp = parameters.lagrangian_physical_properties; + + maximum_particle_diameter = 0; + for (unsigned int particle_type = 0; particle_type < lpp.particle_type_number; + particle_type++) + { + switch (lpp.distribution_type.at(particle_type)) { - size_distribution_object_container[particle_type] = - std::make_shared( - parameters.lagrangian_physical_properties.particle_custom_diameter - .at(particle_type), - parameters.lagrangian_physical_properties - .particle_custom_probability.at(particle_type), - parameters.lagrangian_physical_properties - .seed_for_distributions[particle_type] + - this_mpi_process); + case SizeDistributionType::uniform: + size_distribution_object_container[particle_type] = + std::make_shared( + lpp.particle_average_diameter.at(particle_type)); + break; + case SizeDistributionType::normal: + size_distribution_object_container[particle_type] = + std::make_shared( + lpp.particle_average_diameter.at(particle_type), + lpp.particle_size_std.at(particle_type), + lpp.seed_for_distributions[particle_type] + this_mpi_process); + break; + case SizeDistributionType::custom: + size_distribution_object_container[particle_type] = + std::make_shared( + lpp.particle_custom_diameter.at(particle_type), + lpp.particle_custom_probability.at(particle_type), + lpp.seed_for_distributions[particle_type] + this_mpi_process); + break; } + maximum_particle_diameter = std::max( maximum_particle_diameter, size_distribution_object_container[particle_type]->find_max_diameter()); @@ -159,16 +188,13 @@ DEMSolver::DEMSolver(DEMSolverParameters dem_parameters) std::pow(parameters.model_parameters.neighborhood_threshold * maximum_particle_diameter, 2); +} - if (this_mpi_process == 0) - input_parameter_inspection(parameters, - pcout, - size_distribution_object_container); - - grid_motion_object = - std::make_shared>(parameters.grid_motion, - simulation_control->get_time_step()); - +template +void +DEMSolver::setup_solid_objects() +{ + // Set up solid objects and carry them in vectors for (unsigned int i_solid = 0; i_solid < parameters.solid_objects->number_solid_surfaces; ++i_solid) @@ -185,33 +211,153 @@ DEMSolver::DEMSolver(DEMSolverParameters dem_parameters) this->parameters.solid_objects->solid_volumes[i_solid], i_solid)); } - // Generate solid objects + // Resize the mesh info containers solid_surfaces_mesh_info.resize(solid_surfaces.size()); solid_volumes_mesh_info.resize(solid_volumes.size()); - // Resize particle_floating_mesh_in_contact + // Simulation has solid objects and resize the container if ((solid_surfaces.size() + solid_volumes.size()) > 0) { - has_solid_objects = true; + action_manager->set_solid_objects_enabled(); contact_manager.particle_floating_mesh_in_contact.resize( solid_surfaces.size() + solid_volumes.size()); } +} - // Check if there are periodic boundaries - for (unsigned int i_bc = 0; - i_bc < parameters.boundary_conditions.bc_types.size(); - ++i_bc) +template +void +DEMSolver::setup_functions_and_pointers() +{ + contact_detection_iteration_check_function = + set_contact_search_iteration_function(); + + // Set insertion object type before the restart because the restart only + // rebuilds the member of the insertion object + insertion_object = set_insertion_type(); + + // Setting chosen contact force, insertion and integration methods + integrator_object = set_integrator_type(); + particle_particle_contact_force_object = + set_particle_particle_contact_force_model(parameters); + particle_wall_contact_force_object = + set_particle_wall_contact_force_model(parameters, triangulation); +} + +template +std::function +DEMSolver::set_contact_search_iteration_function() +{ + using namespace Parameters::Lagrangian; + ModelParameters::ContactDetectionMethod &contact_detection_method = + parameters.model_parameters.contact_detection_method; + + switch (contact_detection_method) { - if (parameters.boundary_conditions.bc_types[i_bc] == - Parameters::Lagrangian::BCDEM::BoundaryType::periodic) + case ModelParameters::ContactDetectionMethod::constant: + return [&] { check_contact_search_iteration_constant(); }; + case ModelParameters::ContactDetectionMethod::dynamic: + return [&] { check_contact_search_iteration_dynamic(); }; + default: + throw(std::runtime_error("Invalid contact detection method.")); + } +} + +template +std::shared_ptr> +DEMSolver::set_insertion_type() +{ + using namespace Parameters::Lagrangian; + InsertionInfo::InsertionMethod insertion_method = + parameters.insertion_info.insertion_method; + + switch (insertion_method) + { + case InsertionInfo::InsertionMethod::file: { - has_periodic_boundaries = true; - break; + return std::make_shared>( + size_distribution_object_container, triangulation, parameters); } + case InsertionInfo::InsertionMethod::list: + { + return std::make_shared>( + size_distribution_object_container, triangulation, parameters); + } + case InsertionInfo::InsertionMethod::plane: + { + return std::make_shared>( + size_distribution_object_container, triangulation, parameters); + } + case InsertionInfo::InsertionMethod::volume: + { + return std::make_shared>( + size_distribution_object_container, + triangulation, + parameters, + maximum_particle_diameter); + } + default: + throw(std::runtime_error("Invalid insertion method.")); } +} - // Assign gravity/acceleration - g = parameters.lagrangian_physical_properties.g; +template +std::shared_ptr> +DEMSolver::set_integrator_type() +{ + using namespace Parameters::Lagrangian; + ModelParameters::IntegrationMethod integration_method = + parameters.model_parameters.integration_method; + + switch (integration_method) + { + case ModelParameters::IntegrationMethod::velocity_verlet: + return std::make_shared>(); + case ModelParameters::IntegrationMethod::explicit_euler: + return std::make_shared>(); + case ModelParameters::IntegrationMethod::gear3: + return std::make_shared>(); + default: + throw(std::runtime_error("Invalid integration method.")); + } +} + +template +void +DEMSolver::setup_triangulation_dependent_parameters() +{ + // Find the smallest contact search frequency criterion between (smallest + // cell size - largest particle radius) and (security factor * (blob + // diameter - 1) * largest particle radius). This value is used in + // find_contact_detection_frequency function + smallest_contact_search_criterion = + std::min((GridTools::minimal_cell_diameter(triangulation) - + maximum_particle_diameter * 0.5), + (parameters.model_parameters.dynamic_contact_search_factor * + (parameters.model_parameters.neighborhood_threshold - 1) * + maximum_particle_diameter * 0.5)); + + // Find the smallest cell size and use this as the floating mesh mapping + // criterion. The edge case comes when the cell are completely square/cubic. + // In that case, every sides of a cell are 2^-0.5 or 3^-0.5 times the + // cell_diameter. We want to refresh the mapping each time the solid-objet + // pass through a cell or there will be late contact detection. Thus, we use + // this value. + smallest_solid_object_mapping_criterion = [&] { + if constexpr (dim == 2) // 2^-0.5 * D_c,min + return 0.70710678118 * GridTools::minimal_cell_diameter(triangulation); + if constexpr (dim == 3) // 3^-0.5 * D_c,min + return 0.57735026919 * GridTools::minimal_cell_diameter(triangulation); + }(); + + // Setup background dof + setup_background_dofs(); + + // Set up the periodic boundaries (if PBC enabled) + periodic_boundaries_object.map_periodic_cells( + triangulation, periodic_boundaries_cells_information); + + // Set up the local and ghost cells (if ASC enabled) + sparse_contacts_object.update_local_and_ghost_cell_set(background_dh); } template @@ -227,7 +373,8 @@ DEMSolver::setup_background_dofs() // Those constraints are not used for any matrix assembly in DEM solver, this // approach comes from CFD-DEM coupling where void fraction constraints are // used to achieve the periodic node mapping. - if (has_sparse_contacts && has_periodic_boundaries) + if (action_manager->check_periodic_boundaries_enabled() && + action_manager->check_sparse_contacts_enabled()) { IndexSet locally_relevant_dofs; DoFTools::extract_locally_relevant_dofs(background_dh, @@ -249,14 +396,13 @@ DEMSolver::setup_background_dofs() } } - template void DEMSolver::load_balance() { - load_balance_step = load_balancing.check_load_balance_iteration(); + load_balancing.check_load_balance_iteration(); - if (!load_balance_step) + if (!action_manager->check_load_balance()) return; TimerOutput::Scope t(this->computing_timer, "Load balancing"); @@ -270,36 +416,16 @@ DEMSolver::load_balance() // Unpack the particle handler after the mesh has been repartitioned particle_handler.unpack_after_coarsening_and_refinement(); - // Update the container with all the combinations of background and - // solid cells - for (unsigned int i_solid = 0; i_solid < solid_surfaces.size(); ++i_solid) - { - solid_surfaces_mesh_info[i_solid] = - solid_surfaces[i_solid]->map_solid_in_background_triangulation( - triangulation); - } - - for (unsigned int i_solid = 0; i_solid < solid_volumes.size(); ++i_solid) - { - solid_volumes_mesh_info[i_solid] = - solid_volumes[i_solid]->map_solid_in_background_triangulation( - triangulation); - } + // If PBC are enabled, update the periodic cells + periodic_boundaries_object.map_periodic_cells( + triangulation, periodic_boundaries_cells_information); - if (has_periodic_boundaries) - { - periodic_boundaries_object.map_periodic_cells( - triangulation, periodic_boundaries_cells_information); - - periodic_offset = - periodic_boundaries_object.get_periodic_offset_distance(); - } + // If ASC is enabled, update the local and ghost cell set + sparse_contacts_object.update_local_and_ghost_cell_set(background_dh); // Update neighbors of cells after load balance contact_manager.update_cell_neighbors(triangulation, - periodic_boundaries_cells_information, - has_periodic_boundaries, - has_solid_objects); + periodic_boundaries_cells_information); boundary_cell_object.build( triangulation, @@ -309,10 +435,9 @@ DEMSolver::load_balance() parameters.mesh.expand_particle_wall_contact_search, pcout); - if (parameters.grid_motion.motion_type != - Parameters::Lagrangian::GridMotion::MotionType::none) - boundary_cell_object.update_boundary_info_after_grid_motion( - updated_boundary_points_and_normal_vectors); + // Update the boundary information (if grid motion) + boundary_cell_object.update_boundary_info_after_grid_motion( + updated_boundary_points_and_normal_vectors); const auto average_minimum_maximum_cells = Utilities::MPI::min_max_avg(triangulation.n_active_cells(), @@ -335,82 +460,37 @@ DEMSolver::load_balance() setup_background_dofs(); } - - -template -inline std::function -DEMSolver::set_contact_search_iteration_function() -{ - if (parameters.model_parameters.contact_detection_method == - Parameters::Lagrangian::ModelParameters::ContactDetectionMethod::constant) - { - return [&] { return check_contact_search_iteration_constant(); }; - } - else if (parameters.model_parameters.contact_detection_method == - Parameters::Lagrangian::ModelParameters::ContactDetectionMethod:: - dynamic) - { - return [&] { return check_contact_search_iteration_dynamic(); }; - } - else - { - throw std::runtime_error( - "Specified contact detection method is not valid"); - } -} - template -inline bool +inline void DEMSolver::check_contact_search_iteration_dynamic() { - bool sorting_in_subdomains_step = - (particles_insertion_step || load_balance_step || contact_detection_step); - - return find_particle_contact_detection_step( + find_particle_contact_detection_step( particle_handler, simulation_control->get_time_step(), smallest_contact_search_criterion, mpi_communicator, - sorting_in_subdomains_step, displacement, (simulation_control->get_step_number() % contact_detection_frequency) == 0); } template -inline bool +inline void DEMSolver::check_contact_search_iteration_constant() { - return ( - (simulation_control->get_step_number() % contact_detection_frequency) == 0); + if ((simulation_control->get_step_number() % contact_detection_frequency) == + 0) + action_manager->contact_detection_step(); } template -bool +void DEMSolver::insert_particles() { if ((simulation_control->get_step_number() % insertion_frequency) == 1) { insertion_object->insert(particle_handler, triangulation, parameters); - return true; - } - return false; -} -template -void -DEMSolver::update_moment_of_inertia( - dealii::Particles::ParticleHandler &particle_handler, - std::vector &MOI) -{ - MOI.resize(torque.size()); - - for (auto &particle : particle_handler) - { - auto particle_properties = particle.get_properties(); - MOI[particle.get_local_index()] = - 0.1 * particle_properties[DEM::PropertiesIndex::mass] * - particle_properties[DEM::PropertiesIndex::dp] * - particle_properties[DEM::PropertiesIndex::dp]; + action_manager->particle_insertion_step(); } } @@ -443,8 +523,8 @@ DEMSolver::particle_wall_contact_force() force); } - // Particle - floating mesh contact force - if (has_solid_objects) + // Particle-solid objects contact force + if (action_manager->check_solid_objects_enabled()) // until refactor { particle_wall_contact_force_object ->calculate_particle_floating_wall_contact_force( @@ -461,7 +541,7 @@ DEMSolver::particle_wall_contact_force() parameters.lagrangian_physical_properties, force); - if (dim == 3) + if constexpr (dim == 3) { particle_point_line_contact_force_object .calculate_particle_line_contact_force( @@ -471,15 +551,33 @@ DEMSolver::particle_wall_contact_force() } } +template +void +DEMSolver::move_solid_objects() +{ + if (!action_manager->check_solid_objects_enabled()) + return; + + // Move the solid triangulations, previous time must be used here + // instead of current time. + for (auto &solid_object : solid_surfaces) + solid_object->move_solid_triangulation( + simulation_control->get_time_step(), + simulation_control->get_previous_time()); + + for (auto &solid_object : solid_volumes) + solid_object->move_solid_triangulation( + simulation_control->get_time_step(), + simulation_control->get_previous_time()); +} + template void DEMSolver::finish_simulation() { // Timer output if (parameters.timer.type == Parameters::Timer::Type::end) - { - this->computing_timer.print_summary(); - } + this->computing_timer.print_summary(); // Testing if (parameters.test.enabled) @@ -535,76 +633,6 @@ DEMSolver::finish_simulation() } } -template -std::shared_ptr> -DEMSolver::set_insertion_type(const DEMSolverParameters ¶meters) -{ - if (parameters.insertion_info.insertion_method == - Parameters::Lagrangian::InsertionInfo::InsertionMethod::file) - { - insertion_object = - std::make_shared>(size_distribution_object_container, - triangulation, - parameters); - } - else if (parameters.insertion_info.insertion_method == - Parameters::Lagrangian::InsertionInfo::InsertionMethod::list) - { - insertion_object = - std::make_shared>(size_distribution_object_container, - triangulation, - parameters); - } - else if (parameters.insertion_info.insertion_method == - Parameters::Lagrangian::InsertionInfo::InsertionMethod::plane) - { - insertion_object = std::make_shared>( - size_distribution_object_container, triangulation, parameters); - } - else if (parameters.insertion_info.insertion_method == - Parameters::Lagrangian::InsertionInfo::InsertionMethod::volume) - { - insertion_object = std::make_shared>( - size_distribution_object_container, - triangulation, - parameters, - maximum_particle_diameter); - } - else - { - throw "The chosen insertion method is invalid"; - } - return insertion_object; -} - -template -std::shared_ptr> -DEMSolver::set_integrator_type(const DEMSolverParameters ¶meters) -{ - if (parameters.model_parameters.integration_method == - Parameters::Lagrangian::ModelParameters::IntegrationMethod:: - velocity_verlet) - { - integrator_object = std::make_shared>(); - } - else if (parameters.model_parameters.integration_method == - Parameters::Lagrangian::ModelParameters::IntegrationMethod:: - explicit_euler) - { - integrator_object = std::make_shared>(); - } - else if (parameters.model_parameters.integration_method == - Parameters::Lagrangian::ModelParameters::IntegrationMethod::gear3) - { - integrator_object = std::make_shared>(); - } - else - { - throw "The chosen integration method is invalid"; - } - return integrator_object; -} - template void DEMSolver::write_output_results() @@ -675,16 +703,20 @@ template void DEMSolver::post_process_results() { - post_processing_object.write_post_processing_results( - triangulation, - grid_pvdhandler, - background_dh, - particle_handler, - parameters, - simulation_control->get_current_time(), - simulation_control->get_step_number(), - mpi_communicator, - sparse_contacts_object); + if (parameters.post_processing.Lagrangian_post_processing && + simulation_control->is_output_iteration()) + { + post_processing_object.write_post_processing_results( + triangulation, + grid_pvdhandler, + background_dh, + particle_handler, + parameters, + simulation_control->get_current_time(), + simulation_control->get_step_number(), + mpi_communicator, + sparse_contacts_object); + } } template @@ -760,129 +792,77 @@ DEMSolver::report_statistics() } template -void -DEMSolver::solve() +inline void +DEMSolver::sort_particles_into_subdomains_and_cells() { - // Reading mesh - read_mesh(parameters.mesh, - parameters.restart.restart, - pcout, - triangulation, - parameters.boundary_conditions); + particle_handler.sort_particles_into_subdomains_and_cells(); - // Store information about floating mesh/background mesh intersection - for (unsigned int i_solid = 0; i_solid < solid_surfaces.size(); ++i_solid) + // Resize the displacement, force and torque containers only if the particles + // have changed subdomains + if (action_manager->check_resize_containers()) { - solid_surfaces_mesh_info[i_solid] = - solid_surfaces[i_solid]->map_solid_in_background_triangulation( - triangulation); - } - - for (unsigned int i_solid = 0; i_solid < solid_volumes.size(); ++i_solid) - { - solid_volumes_mesh_info[i_solid] = - solid_volumes[i_solid]->map_solid_in_background_triangulation( - triangulation); - } - - // Set insertion object type before the restart because the restart only - // rebuilds the member of the insertion object - insertion_object = set_insertion_type(parameters); - - contact_detection_iteration_check_function = - set_contact_search_iteration_function(); - - if (parameters.restart.restart == true) - { - read_checkpoint(computing_timer, - parameters, - simulation_control, - particles_pvdhandler, - grid_pvdhandler, - triangulation, - particle_handler, - insertion_object, - solid_surfaces); - + // Resize and reinitialize displacement container displacement.resize(particle_handler.get_max_local_particle_index()); + // Resize and reinitialize displacement container force.resize(displacement.size()); torque.resize(displacement.size()); + MOI.resize(displacement.size()); - update_moment_of_inertia(particle_handler, MOI); - - checkpoint_step = true; - } - - // Store information about floating mesh/background mesh intersection - for (unsigned int i_solid = 0; i_solid < solid_surfaces.size(); ++i_solid) - { - solid_surfaces_mesh_info[i_solid] = - solid_surfaces[i_solid]->map_solid_in_background_triangulation( - triangulation); + // Updating moment of inertia container + for (auto &particle : particle_handler) + { + auto particle_properties = particle.get_properties(); + MOI[particle.get_local_index()] = + 0.1 * particle_properties[DEM::PropertiesIndex::mass] * + particle_properties[DEM::PropertiesIndex::dp] * + particle_properties[DEM::PropertiesIndex::dp]; + } } - for (unsigned int i_solid = 0; i_solid < solid_volumes.size(); ++i_solid) - { - solid_volumes_mesh_info[i_solid] = - solid_volumes[i_solid]->map_solid_in_background_triangulation( - triangulation); - } + // Always reset the displacement values since we are doing a search detection + std::fill(displacement.begin(), displacement.end(), 0.); - // Find the smallest contact search frequency criterion between (smallest - // cell size - largest particle radius) and (security factor * (blob - // diameter - 1) * largest particle radius). This value is used in - // find_contact_detection_frequency function - smallest_contact_search_criterion = - std::min((GridTools::minimal_cell_diameter(triangulation) - - maximum_particle_diameter * 0.5), - (parameters.model_parameters.dynamic_contact_search_factor * - (parameters.model_parameters.neighborhood_threshold - 1) * - maximum_particle_diameter * 0.5)); - - // Find the smallest cell size and use this as the floating mesh mapping - // criterion - - double mapping_criterion_constant; - if constexpr (dim == 2) - mapping_criterion_constant = 0.70710678118; // 2^-0.5 - - if constexpr (dim == 3) - mapping_criterion_constant = 0.57735026919; // 3^-0.5 + // Exchange ghost particles + particle_handler.exchange_ghost_particles(true); +} - smallest_floating_mesh_mapping_criterion = - mapping_criterion_constant * - GridTools::minimal_cell_diameter(triangulation); - // The edge case comes when the cell are completely square/cubic. In that - // case, every sides of a cell are 2^-0.5 or 3^-0.5 times the cell_diameter. - // We want to refresh the mapping each time the solid-objet pass through a - // cell or there will be late contact detection. Thus, we use this value. +template +void +DEMSolver::solve() +{ + // Set up the parameters + setup_parameters(); + // Reading mesh + read_mesh(parameters.mesh, + action_manager->check_restart_simulation(), + pcout, + triangulation, + parameters.boundary_conditions); - if (has_periodic_boundaries) - { - periodic_boundaries_object.set_periodic_boundaries_information( - parameters.boundary_conditions.periodic_boundary_0, - parameters.boundary_conditions.periodic_boundary_1, - parameters.boundary_conditions.periodic_direction); + // Set up functions and pointers according to parameters + setup_functions_and_pointers(); - periodic_boundaries_object.map_periodic_cells( - triangulation, periodic_boundaries_cells_information); + // Read checkpoint if needed + read_checkpoint(computing_timer, + parameters, + simulation_control, + particles_pvdhandler, + grid_pvdhandler, + triangulation, + particle_handler, + insertion_object, + solid_surfaces); - // Temporary offset calculation : works only for one set of periodic - // boundary on an axis. - periodic_offset = - periodic_boundaries_object.get_periodic_offset_distance(); - } + // Set up the various parameters that need the triangulation + setup_triangulation_dependent_parameters(); - // Find cell neighbors + // Build the mapping of the cell neighbors contact_manager.execute_cell_neighbors_search( - triangulation, - periodic_boundaries_cells_information, - has_periodic_boundaries, - has_solid_objects); + triangulation, periodic_boundaries_cells_information); - // Finding boundary cells with faces + // Find boundary cells with faces boundary_cell_object.build( triangulation, parameters.floating_walls, @@ -891,175 +871,109 @@ DEMSolver::solve() parameters.mesh.expand_particle_wall_contact_search, pcout); - // Setting chosen contact force, insertion and integration methods - integrator_object = set_integrator_type(parameters); - particle_particle_contact_force_object = - set_particle_particle_contact_force_model(parameters); - particle_wall_contact_force_object = - set_particle_wall_contact_force_model(parameters, triangulation); - - // Setup background dof - setup_background_dofs(); - - // DEM engine iterator: + // DEM engine iterator while (simulation_control->integrate()) { simulation_control->print_progression(pcout); if (simulation_control->is_verbose_iteration()) report_statistics(); - // Grid motion - if (parameters.grid_motion.motion_type != - Parameters::Lagrangian::GridMotion::MotionType::none) - { - grid_motion_object->move_grid(triangulation); - boundary_cell_object.update_boundary_info_after_grid_motion( - updated_boundary_points_and_normal_vectors); - } + // Move grid and update the boundary information (if grid motion) + grid_motion_object->move_grid(triangulation); + boundary_cell_object.update_boundary_info_after_grid_motion( + updated_boundary_points_and_normal_vectors); - // Keep track if particles were inserted this step - particles_insertion_step = insert_particles(); + // Insert particle if needed + insert_particles(); - if (particles_insertion_step) - displacement.resize(particle_handler.get_max_local_particle_index()); - - // Load balancing + // Load balancing (if load balancing enabled and if needed) load_balance(); - if (load_balance_step || checkpoint_step) - { - displacement.resize(particle_handler.get_max_local_particle_index()); - - if (has_sparse_contacts) - { - sparse_contacts_object.update_local_and_ghost_cell_set( - background_dh); - } - } + // Check for contact search according to the contact detection method + contact_detection_iteration_check_function(); - // Check to see if it is contact search step - contact_detection_step = contact_detection_iteration_check_function(); + // Check if solid object needs to be mapped with background mesh + // (if solid object) + find_floating_mesh_mapping_step(smallest_solid_object_mapping_criterion, + this->solid_surfaces); - bool solid_object_map_step = false; - // Check to see if floating meshes need to be mapped in background mesh - if (has_solid_objects) + // Map solid objects if the action was triggered (if solid object) + if (action_manager->check_solid_object_search()) { - solid_object_map_step = find_floating_mesh_mapping_step( - smallest_floating_mesh_mapping_criterion, this->solid_surfaces); + // Store information about floating mesh/background mesh + // intersection + for (unsigned int i_solid = 0; i_solid < solid_surfaces.size(); + ++i_solid) + { + solid_surfaces_mesh_info[i_solid] = + solid_surfaces[i_solid]->map_solid_in_background_triangulation( + triangulation); + } - if (solid_object_map_step) + for (unsigned int i_solid = 0; i_solid < solid_volumes.size(); + ++i_solid) { - // Update floating mesh information in the container manager - for (unsigned int i_solid = 0; i_solid < solid_surfaces.size(); - ++i_solid) - { - solid_surfaces_mesh_info[i_solid] = - solid_surfaces[i_solid] - ->map_solid_in_background_triangulation(triangulation); - } - - for (unsigned int i_solid = 0; i_solid < solid_volumes.size(); - ++i_solid) - { - } + solid_volumes_mesh_info[i_solid] = + solid_volumes[i_solid]->map_solid_in_background_triangulation( + triangulation); } } - contact_detection_step = contact_detection_step || solid_object_map_step; - - // Sort particles in cells - if (particles_insertion_step || load_balance_step || - contact_detection_step || checkpoint_step) + // Execute contact search if the action was triggered + if (action_manager->check_contact_search()) { // Particles displacement if passing through a periodic boundary + // (if PBC enabled) periodic_boundaries_object.execute_particles_displacement( particle_handler, periodic_boundaries_cells_information); - particle_handler.sort_particles_into_subdomains_and_cells(); - - if (has_sparse_contacts && !simulation_control->is_at_start()) - { - // Compute cell mobility for all cells after the sort particles - // into subdomains and cells and exchange ghost particles - sparse_contacts_object.identify_mobility_status( - background_dh, - particle_handler, - triangulation.n_active_cells(), - mpi_communicator); - } - - displacement.resize(particle_handler.get_max_local_particle_index()); - force.resize(displacement.size()); - torque.resize(displacement.size()); + sort_particles_into_subdomains_and_cells(); - // Updating moment of inertia container - update_moment_of_inertia(particle_handler, MOI); - - particle_handler.exchange_ghost_particles(true); - - // Reset checkpoint step - checkpoint_step = false; + // Compute cell mobility (if ASC enabled) + sparse_contacts_object.identify_mobility_status( + background_dh, + particle_handler, + triangulation.n_active_cells(), + mpi_communicator); // Execute broad search by filling containers of particle-particle // contact pair candidates and containers of particle-wall // contact pair candidates - if (!(has_sparse_contacts && contact_build_number > 1)) - { - contact_manager.execute_particle_particle_broad_search( - particle_handler, has_periodic_boundaries); - - contact_manager.execute_particle_wall_broad_search( - particle_handler, - boundary_cell_object, - solid_surfaces_mesh_info, - parameters.floating_walls, - simulation_control->get_current_time(), - has_solid_objects); - } - else - { - contact_manager.execute_particle_particle_broad_search( - particle_handler, - sparse_contacts_object, - has_periodic_boundaries); - - contact_manager.execute_particle_wall_broad_search( - particle_handler, - boundary_cell_object, - solid_surfaces_mesh_info, - parameters.floating_walls, - simulation_control->get_current_time(), - sparse_contacts_object, - has_solid_objects); - } + contact_manager.execute_particle_particle_broad_search( + particle_handler, sparse_contacts_object); - // Updating number of contact builds - contact_build_number++; + contact_manager.execute_particle_wall_broad_search( + particle_handler, + boundary_cell_object, + solid_surfaces_mesh_info, + parameters.floating_walls, + simulation_control->get_current_time(), + sparse_contacts_object); // Update contacts, remove replicates and add new contact pairs // to the contact containers when particles are exchanged between // processors - contact_manager.update_contacts(has_periodic_boundaries); + contact_manager.update_contacts(); // Updates the iterators to particles in local-local contact // containers - contact_manager.update_local_particles_in_cells( - particle_handler, load_balance_step, has_periodic_boundaries); + contact_manager.update_local_particles_in_cells(particle_handler); // Execute fine search by updating particle-particle contact // containers according to the neighborhood threshold contact_manager.execute_particle_particle_fine_search( neighborhood_threshold_squared, - has_periodic_boundaries, - periodic_offset); + periodic_boundaries_object.get_periodic_offset_distance()); - // Execute fine search by updating particle-wall contact containers - // according to the neighborhood threshold + // Execute fine search by updating particle-wall contact + // containers according to the neighborhood threshold contact_manager.execute_particle_wall_fine_search( parameters.floating_walls, simulation_control->get_current_time(), - neighborhood_threshold_squared, - has_solid_objects); + neighborhood_threshold_squared); + + // Updating number of contact builds + contact_build_number++; } else { @@ -1073,39 +987,28 @@ DEMSolver::solve() simulation_control->get_time_step(), torque, force, - periodic_offset); + periodic_boundaries_object.get_periodic_offset_distance()); + // Update the boundary points and vectors (if grid motion) // We have to update the positions of the points on boundary faces and // their normal vectors here. The update_contacts deletes the - // particle-wall contact candidate if it exists in the contact list. As - // a result, when we update the points on boundary faces and their + // particle-wall contact candidate if it exists in the contact list. + // As a result, when we update the points on boundary faces and their // normal vectors, update_contacts deletes it from the output of broad // search and they are not updated in the contact force calculations. - if (parameters.grid_motion.motion_type != - Parameters::Lagrangian::GridMotion::MotionType::none) - grid_motion_object - ->update_boundary_points_and_normal_vectors_in_contact_list( - contact_manager.particle_wall_in_contact, - updated_boundary_points_and_normal_vectors); - - // Move the solid triangulations, previous time must be used here - // instead of current time. - for (auto &solid_object : solid_surfaces) - solid_object->move_solid_triangulation( - simulation_control->get_time_step(), - simulation_control->get_previous_time()); + grid_motion_object + ->update_boundary_points_and_normal_vectors_in_contact_list( + contact_manager.particle_wall_in_contact, + updated_boundary_points_and_normal_vectors); - for (auto &solid_object : solid_volumes) - solid_object->move_solid_triangulation( - simulation_control->get_time_step(), - simulation_control->get_previous_time()); + // Move solid objects (if solid object) + move_solid_objects(); - // Particles-walls contact force: + // Particle-wall contact force particle_wall_contact_force(); - // Integration correction step (after force calculation) - // In the first step, we have to obtain location of particles at - // half-step time + // Integration of force and velocity for new location of particles + // The half step is calculated at the first iteration if (simulation_control->get_step_number() == 0) { integrator_object->integrate_half_step_location( @@ -1116,54 +1019,44 @@ DEMSolver::solve() force, MOI); } - else // Step number > 0 + else { - if (!(has_sparse_contacts && contact_build_number > 1)) - { - integrator_object->integrate(particle_handler, - g, - simulation_control->get_time_step(), - torque, - force, - MOI); - } - else - { - integrator_object->integrate(particle_handler, - g, - simulation_control->get_time_step(), - torque, - force, - MOI, - triangulation, - sparse_contacts_object); - } + integrator_object->integrate(particle_handler, + g, + simulation_control->get_time_step(), + torque, + force, + MOI, + triangulation, + sparse_contacts_object); } // Visualization if (simulation_control->is_output_iteration()) - { - write_output_results(); - } + write_output_results(); - if (parameters.forces_torques.calculate_force_torque && - (this_mpi_process == 0) && - (simulation_control->get_step_number() % - parameters.forces_torques.output_frequency == - 0) && - (parameters.forces_torques.force_torque_verbosity == - Parameters::Verbosity::verbose)) - write_forces_torques_output_locally( - forces_boundary_information[simulation_control->get_step_number()], - torques_boundary_information[simulation_control->get_step_number()]); - - // Post-processing - if (parameters.post_processing.Lagrangian_post_processing && - simulation_control->is_output_iteration()) + // Calculation of forces and torques if needed + if (parameters.forces_torques.calculate_force_torque) { - post_process_results(); + if ((this_mpi_process == 0) && + (simulation_control->get_step_number() % + parameters.forces_torques.output_frequency == + 0) && + (parameters.forces_torques.force_torque_verbosity == + Parameters::Verbosity::verbose)) + { + write_forces_torques_output_locally( + forces_boundary_information[simulation_control + ->get_step_number()], + torques_boundary_information[simulation_control + ->get_step_number()]); + } } + // Post-processing if needed + post_process_results(); + + // Write checkpoint if needed if (parameters.restart.checkpoint && simulation_control->get_step_number() % parameters.restart.frequency == @@ -1181,6 +1074,9 @@ DEMSolver::solve() pcout, mpi_communicator); } + + // Reset all trigger flags + action_manager->reset_triggers(); } finish_simulation(); diff --git a/source/dem/dem_action_manager.cc b/source/dem/dem_action_manager.cc new file mode 100644 index 0000000000..7da00eb6e3 --- /dev/null +++ b/source/dem/dem_action_manager.cc @@ -0,0 +1,12 @@ +#include + +DEMActionManager *DEMActionManager::instance = nullptr; + +DEMActionManager * +DEMActionManager::get_action_manager() +{ + if (instance == nullptr) + instance = new DEMActionManager(); + + return instance; +} diff --git a/source/dem/dem_contact_manager.cc b/source/dem/dem_contact_manager.cc index d98f1540fd..1ed19303c6 100644 --- a/source/dem/dem_contact_manager.cc +++ b/source/dem/dem_contact_manager.cc @@ -1,3 +1,4 @@ +#include #include template @@ -5,17 +6,18 @@ void DEMContactManager::execute_cell_neighbors_search( const parallel::distributed::Triangulation &triangulation, const typename DEM::dem_data_structures::periodic_boundaries_cells_info - periodic_boundaries_cells_information, - const bool has_periodic_boundaries, - const bool has_solid_objects) + periodic_boundaries_cells_information) { + // Get action manager + auto action_manager = DEMActionManager::get_action_manager(); + // Find cell neighbors cell_neighbors_object.find_cell_neighbors(triangulation, cells_local_neighbor_list, cells_ghost_neighbor_list); // Find cell periodic neighbors - if (has_periodic_boundaries) + if (action_manager->check_periodic_boundaries_enabled()) { cell_neighbors_object.find_cell_periodic_neighbors( triangulation, @@ -26,7 +28,7 @@ DEMContactManager::execute_cell_neighbors_search( } // Get total (with repetition) neighbors list for floating mesh. - if (has_solid_objects) + if (action_manager->check_solid_objects_enabled()) { cell_neighbors_object.find_full_cell_neighbors(triangulation, total_neighbor_list); @@ -35,7 +37,7 @@ DEMContactManager::execute_cell_neighbors_search( template void -DEMContactManager::update_contacts(const bool has_periodic_boundaries) +DEMContactManager::update_contacts() { // Update particle-particle contacts in local_adjacent_particles of fine // search step with local_contact_pair_candidates @@ -55,7 +57,8 @@ DEMContactManager::update_contacts(const bool has_periodic_boundaries) ContactType::ghost_particle_particle>(ghost_adjacent_particles, ghost_contact_pair_candidates); - if (has_periodic_boundaries) + if (DEMActionManager::get_action_manager() + ->check_periodic_boundaries_enabled()) { // Update periodic particle-particle contacts in // local_periodic_adjacent_particles of fine search step with @@ -130,9 +133,7 @@ DEMContactManager::update_contacts(const bool has_periodic_boundaries) template void DEMContactManager::update_local_particles_in_cells( - const Particles::ParticleHandler &particle_handler, - const bool clear_contact_structures, - const bool has_periodic_boundaries) + const Particles::ParticleHandler &particle_handler) { // Update the iterators to local particles in a map of particles update_particle_container(particle_container, &particle_handler); @@ -142,18 +143,17 @@ DEMContactManager::update_local_particles_in_cells( dim, typename DEM::dem_data_structures::adjacent_particle_pairs, ContactType::local_particle_particle>(local_adjacent_particles, - particle_container, - clear_contact_structures); + particle_container); // Update contact containers for ghost particle-particle pairs in contact update_contact_container_iterators< dim, typename DEM::dem_data_structures::adjacent_particle_pairs, ContactType::ghost_particle_particle>(ghost_adjacent_particles, - particle_container, - clear_contact_structures); + particle_container); - if (has_periodic_boundaries) + if (DEMActionManager::get_action_manager() + ->check_periodic_boundaries_enabled()) { // Update contact containers for local-local periodic particle-particle // pairs in contact @@ -161,9 +161,7 @@ DEMContactManager::update_local_particles_in_cells( dim, typename DEM::dem_data_structures::adjacent_particle_pairs, ContactType::local_periodic_particle_particle>( - local_periodic_adjacent_particles, - particle_container, - clear_contact_structures); + local_periodic_adjacent_particles, particle_container); // Update contact containers for local-ghost periodic particle-particle // pairs in contact @@ -171,9 +169,7 @@ DEMContactManager::update_local_particles_in_cells( dim, typename DEM::dem_data_structures::adjacent_particle_pairs, ContactType::ghost_periodic_particle_particle>( - ghost_periodic_adjacent_particles, - particle_container, - clear_contact_structures); + ghost_periodic_adjacent_particles, particle_container); // Update contact containers for ghost-local periodic particle-particle // pairs in contact @@ -181,9 +177,7 @@ DEMContactManager::update_local_particles_in_cells( dim, typename DEM::dem_data_structures::adjacent_particle_pairs, ContactType::ghost_local_periodic_particle_particle>( - ghost_local_periodic_adjacent_particles, - particle_container, - clear_contact_structures); + ghost_local_periodic_adjacent_particles, particle_container); } // Update contact containers for particle-wall pairs in contact @@ -231,34 +225,37 @@ template void DEMContactManager::execute_particle_particle_broad_search( dealii::Particles::ParticleHandler &particle_handler, - const bool has_periodic_boundaries) + const AdaptiveSparseContacts &sparse_contacts_object) { - particle_particle_broad_search_object.find_particle_particle_contact_pairs( - particle_handler, *this); + auto action_manager = DEMActionManager::get_action_manager(); - if (has_periodic_boundaries) + // Check if sparse contacts are enabled to use proper broad search functions + // The first broad search is the default one for sparse contacts + if (action_manager->use_default_broad_search_functions()) { particle_particle_broad_search_object - .find_particle_particle_periodic_contact_pairs(particle_handler, *this); + .find_particle_particle_contact_pairs(particle_handler, *this); + + if (action_manager->check_periodic_boundaries_enabled()) + { + particle_particle_broad_search_object + .find_particle_particle_periodic_contact_pairs(particle_handler, + *this); + } } -} - -template -void -DEMContactManager::execute_particle_particle_broad_search( - dealii::Particles::ParticleHandler &particle_handler, - const AdaptiveSparseContacts &sparse_contacts_object, - const bool has_periodic_boundaries) -{ - particle_particle_broad_search_object.find_particle_particle_contact_pairs( - particle_handler, *this, sparse_contacts_object); - - if (has_periodic_boundaries) + else { particle_particle_broad_search_object - .find_particle_particle_periodic_contact_pairs(particle_handler, - *this, - sparse_contacts_object); + .find_particle_particle_contact_pairs(particle_handler, + *this, + sparse_contacts_object); + + if (action_manager->check_periodic_boundaries_enabled()) + { + particle_particle_broad_search_object + .find_particle_particle_periodic_contact_pairs( + particle_handler, *this, sparse_contacts_object); + } } } @@ -271,107 +268,106 @@ DEMContactManager::execute_particle_wall_broad_search( solid_surfaces_mesh_info, const Parameters::Lagrangian::FloatingWalls &floating_walls, const double simulation_time, - const bool has_solid_objects) + const AdaptiveSparseContacts &sparse_contacts_object) { - // Particle-wall contact candidates - particle_wall_broad_search_object.find_particle_wall_contact_pairs( - boundary_cell_object.get_boundary_cells_information(), - particle_handler, - particle_wall_candidates); + auto action_manager = DEMActionManager::get_action_manager(); - // Particle-floating wall contact pairs - if (floating_walls.floating_walls_number > 0) + if (action_manager->use_default_broad_search_functions()) { - particle_wall_broad_search_object - .find_particle_floating_wall_contact_pairs( - boundary_cell_object.get_boundary_cells_with_floating_walls(), - particle_handler, - floating_walls, - simulation_time, - particle_floating_wall_candidates); - } - - // Particle-floating mesh broad search - if (has_solid_objects) - { - particle_wall_broad_search_object.particle_solid_surfaces_contact_search( - solid_surfaces_mesh_info, + // Particle-wall contact candidates + particle_wall_broad_search_object.find_particle_wall_contact_pairs( + boundary_cell_object.get_boundary_cells_information(), particle_handler, - particle_floating_mesh_candidates, - total_neighbor_list); - } - - particle_point_candidates = - particle_point_line_broad_search_object.find_particle_point_contact_pairs( - particle_handler, boundary_cell_object.get_boundary_cells_with_points()); - - if constexpr (dim == 3) - { - particle_line_candidates = + particle_wall_candidates); + + // Particle-floating wall contact pairs + if (floating_walls.floating_walls_number > 0) + { + particle_wall_broad_search_object + .find_particle_floating_wall_contact_pairs( + boundary_cell_object.get_boundary_cells_with_floating_walls(), + particle_handler, + floating_walls, + simulation_time, + particle_floating_wall_candidates); + } + + // Particle-floating mesh broad search + if (action_manager->check_solid_objects_enabled()) + { + particle_wall_broad_search_object + .particle_solid_surfaces_contact_search( + solid_surfaces_mesh_info, + particle_handler, + particle_floating_mesh_candidates, + total_neighbor_list); + } + + particle_point_candidates = particle_point_line_broad_search_object - .find_particle_line_contact_pairs( + .find_particle_point_contact_pairs( particle_handler, - boundary_cell_object.get_boundary_cells_with_lines()); + boundary_cell_object.get_boundary_cells_with_points()); + + if constexpr (dim == 3) + { + particle_line_candidates = + particle_point_line_broad_search_object + .find_particle_line_contact_pairs( + particle_handler, + boundary_cell_object.get_boundary_cells_with_lines()); + } } -} - -template -void -DEMContactManager::execute_particle_wall_broad_search( - const Particles::ParticleHandler &particle_handler, - BoundaryCellsInformation &boundary_cell_object, - const typename dem_data_structures::solid_surfaces_mesh_information - solid_surfaces_mesh_info, - const Parameters::Lagrangian::FloatingWalls &floating_walls, - const double simulation_time, - const AdaptiveSparseContacts &sparse_contacts_object, - const bool has_solid_objects) -{ - // Particle-wall contact candidates - particle_wall_broad_search_object.find_particle_wall_contact_pairs( - boundary_cell_object.get_boundary_cells_information(), - particle_handler, - particle_wall_candidates, - sparse_contacts_object); - - // Particle-floating wall contact pairs - if (floating_walls.floating_walls_number > 0) - { - particle_wall_broad_search_object - .find_particle_floating_wall_contact_pairs( - boundary_cell_object.get_boundary_cells_with_floating_walls(), - particle_handler, - floating_walls, - simulation_time, - particle_floating_wall_candidates, - sparse_contacts_object); - } - - // Particle-floating mesh broad search - if (has_solid_objects) + else { - particle_wall_broad_search_object.particle_solid_surfaces_contact_search( - solid_surfaces_mesh_info, + // Particle-wall contact candidates + particle_wall_broad_search_object.find_particle_wall_contact_pairs( + boundary_cell_object.get_boundary_cells_information(), particle_handler, - particle_floating_mesh_candidates, - total_neighbor_list, + particle_wall_candidates, sparse_contacts_object); - } - - particle_point_candidates = - particle_point_line_broad_search_object.find_particle_point_contact_pairs( - particle_handler, - boundary_cell_object.get_boundary_cells_with_points(), - sparse_contacts_object); - if constexpr (dim == 3) - { - particle_line_candidates = + // Particle-floating wall contact pairs + if (floating_walls.floating_walls_number > 0) + { + particle_wall_broad_search_object + .find_particle_floating_wall_contact_pairs( + boundary_cell_object.get_boundary_cells_with_floating_walls(), + particle_handler, + floating_walls, + simulation_time, + particle_floating_wall_candidates, + sparse_contacts_object); + } + + // Particle-floating mesh broad search + if (action_manager->check_solid_objects_enabled()) + { + particle_wall_broad_search_object + .particle_solid_surfaces_contact_search( + solid_surfaces_mesh_info, + particle_handler, + particle_floating_mesh_candidates, + total_neighbor_list, + sparse_contacts_object); + } + + particle_point_candidates = particle_point_line_broad_search_object - .find_particle_line_contact_pairs( + .find_particle_point_contact_pairs( particle_handler, - boundary_cell_object.get_boundary_cells_with_lines(), + boundary_cell_object.get_boundary_cells_with_points(), sparse_contacts_object); + + if constexpr (dim == 3) + { + particle_line_candidates = + particle_point_line_broad_search_object + .find_particle_line_contact_pairs( + particle_handler, + boundary_cell_object.get_boundary_cells_with_lines(), + sparse_contacts_object); + } } } @@ -379,7 +375,6 @@ template void DEMContactManager::execute_particle_particle_fine_search( const double neighborhood_threshold, - const bool has_periodic_boundaries, const Tensor<1, dim> periodic_offset) { // Fine search for local particle-particle @@ -396,7 +391,8 @@ DEMContactManager::execute_particle_particle_fine_search( ghost_contact_pair_candidates, neighborhood_threshold); - if (has_periodic_boundaries) + if (DEMActionManager::get_action_manager() + ->check_periodic_boundaries_enabled()) { // Fine search for local-local periodic particle-particle particle_particle_fine_search_object.particle_particle_fine_search( @@ -429,8 +425,7 @@ void DEMContactManager::execute_particle_wall_fine_search( const Parameters::Lagrangian::FloatingWalls &floating_walls, const double simulation_time, - const double neighborhood_threshold, - const bool has_solid_objects) + const double neighborhood_threshold) { // Particle - wall fine search particle_wall_fine_search_object.particle_wall_fine_search( @@ -447,7 +442,7 @@ DEMContactManager::execute_particle_wall_fine_search( } // Particle - floating mesh fine search - if (has_solid_objects) + if (DEMActionManager::get_action_manager()->check_solid_objects_enabled()) { particle_wall_fine_search_object.particle_floating_mesh_fine_search( particle_floating_mesh_candidates, particle_floating_mesh_in_contact); diff --git a/source/dem/explicit_euler_integrator.cc b/source/dem/explicit_euler_integrator.cc index cd3b1eabe7..1218dbb497 100644 --- a/source/dem/explicit_euler_integrator.cc +++ b/source/dem/explicit_euler_integrator.cc @@ -91,15 +91,27 @@ ExplicitEulerIntegrator::integrate( template void ExplicitEulerIntegrator::integrate( - Particles::ParticleHandler & /* particle_handler */, - const Tensor<1, 3> & /* g */, - const double /* dt */, - std::vector> & /* torque */, - std::vector> & /* force */, - const std::vector & /* MOI */, + Particles::ParticleHandler &particle_handler, + const Tensor<1, 3> &g, + const double dt, + std::vector> &torque, + std::vector> &force, + const std::vector &MOI, const parallel::distributed::Triangulation & /* triangulation */, AdaptiveSparseContacts & /* sparse_contacts_object */) { + auto action_manager = DEMActionManager::get_action_manager(); + + bool use_default_function = + !action_manager->check_sparse_contacts_enabled() || + action_manager->check_mobility_status_reset(); + + if (use_default_function) + { + integrate(particle_handler, g, dt, torque, force, MOI); + return; + } + throw std::runtime_error( "Adaptive sparse contacts are not supported with explicit Euler integrator, use Velocity Verlet integrator."); } diff --git a/source/dem/find_boundary_cells_information.cc b/source/dem/find_boundary_cells_information.cc index f7f94fcc23..4ef588aee6 100644 --- a/source/dem/find_boundary_cells_information.cc +++ b/source/dem/find_boundary_cells_information.cc @@ -1,6 +1,7 @@ #include #include +#include #include #include @@ -229,6 +230,10 @@ BoundaryCellsInformation::update_boundary_info_after_grid_motion( typename DEM::dem_data_structures::boundary_points_and_normal_vectors &updated_boundary_points_and_normal_vectors) { + // If there is no grid motion, exit the function + if (!DEMActionManager::get_action_manager()->check_grid_motion_enabled()) + return; + // Initialize a simple quadrature for on the system. This will be used to // obtain a single sample point on the boundary faces const FE_Q fe(1); diff --git a/source/dem/find_contact_detection_step.cc b/source/dem/find_contact_detection_step.cc index 3d0130fc3b..ad7f110306 100644 --- a/source/dem/find_contact_detection_step.cc +++ b/source/dem/find_contact_detection_step.cc @@ -1,21 +1,24 @@ +#include #include using namespace dealii; template -bool +void find_particle_contact_detection_step( Particles::ParticleHandler &particle_handler, const double dt, const double smallest_contact_search_criterion, MPI_Comm &mpi_communicator, - const bool sorting_in_subdomains_step, std::vector &displacement, const bool parallel_update) { - if (sorting_in_subdomains_step) - for (auto &d : displacement) - d = 0.; + // Get the action manager + auto action_manager = DEMActionManager::get_action_manager(); + // If something else has already triggered contact search, + // no need to do it again + if (action_manager->check_contact_search()) + return; double max_displacement = 0.; bool contact_detection_step = false; @@ -41,48 +44,50 @@ find_particle_contact_detection_step( // If the maximum displacement of particles exceeds criterion, this step // is a contact detection step - if (max_displacement > smallest_contact_search_criterion) - { - contact_detection_step = true; - } + contact_detection_step = max_displacement > smallest_contact_search_criterion; // Broadcasting contact detection step value to other processors if (parallel_update) { - return Utilities::MPI::logical_or(contact_detection_step, - mpi_communicator); + contact_detection_step = + Utilities::MPI::logical_or(contact_detection_step, mpi_communicator); + if (contact_detection_step) + action_manager->contact_detection_step(); } - else - return false; } -template bool +template void find_particle_contact_detection_step( Particles::ParticleHandler<2> &particle_handler, const double dt, const double smallest_contact_search_criterion, MPI_Comm &mpi_communicator, - bool sorting_in_subdomains_step, std::vector &displacement, const bool parallel_update); -template bool +template void find_particle_contact_detection_step( Particles::ParticleHandler<3> &particle_handler, const double dt, const double smallest_contact_search_criterion, MPI_Comm &mpi_communicator, - const bool sorting_in_subdomains_step, std::vector &displacement, const bool parallel_update); template -bool +void find_floating_mesh_mapping_step( const double smallest_contact_search_criterion, std::vector>> solids) { + // Get the action manager + auto action_manager = DEMActionManager::get_action_manager(); + + // If there is no solid object, no need to do anything + if (!action_manager->check_solid_objects_enabled()) + return; + bool floating_mesh_requires_map = false; for (unsigned int i_solid = 0; i_solid < solids.size(); ++i_solid) @@ -94,15 +99,16 @@ find_floating_mesh_mapping_step( displacement > smallest_contact_search_criterion; } - return floating_mesh_requires_map; + if (floating_mesh_requires_map) + action_manager->solid_objects_search_step(); } -template bool +template void find_floating_mesh_mapping_step( const double smallest_contact_search_criterion, std::vector>> solids); -template bool +template void find_floating_mesh_mapping_step( const double smallest_contact_search_criterion, std::vector>> solids); diff --git a/source/dem/gear3_integrator.cc b/source/dem/gear3_integrator.cc index 6ac7bb0e99..a578d2083c 100644 --- a/source/dem/gear3_integrator.cc +++ b/source/dem/gear3_integrator.cc @@ -81,15 +81,27 @@ for (auto particle = particle_handler.begin(); template void Gear3Integrator::integrate( - Particles::ParticleHandler & /* particle_handler */, - const Tensor<1, 3> & /* g */, - const double /* dt */, - std::vector> & /* torque */, - std::vector> & /* force */, - const std::vector & /* MOI */, + Particles::ParticleHandler &particle_handler, + const Tensor<1, 3> &g, + const double dt, + std::vector> &torque, + std::vector> &force, + const std::vector &MOI, const parallel::distributed::Triangulation & /* triangulation */, - AdaptiveSparseContacts & /* disable_contacts_object */) + AdaptiveSparseContacts & /* sparse_contacts_object */) { + auto action_manager = DEMActionManager::get_action_manager(); + + bool use_default_function = + !action_manager->check_sparse_contacts_enabled() || + action_manager->check_mobility_status_reset(); + + if (use_default_function) + { + integrate(particle_handler, g, dt, torque, force, MOI); + return; + } + throw std::runtime_error( "Adaptive sparse contacts not supported with explicit Gear 3 integrator, use Verlet integrator."); } diff --git a/source/dem/grid_motion.cc b/source/dem/grid_motion.cc index aeacaa416f..3d6789946b 100644 --- a/source/dem/grid_motion.cc +++ b/source/dem/grid_motion.cc @@ -1,3 +1,4 @@ +#include #include #include @@ -12,22 +13,40 @@ GridMotion::GridMotion( const Parameters::Lagrangian::GridMotion &grid_motion_parameters, const double dem_time_step) { - // Setting grid motion type - if (grid_motion_parameters.motion_type == - Parameters::Lagrangian::GridMotion::MotionType::rotational) - { - grid_motion = &GridMotion::move_grid_rotational; - rotation_angle = - grid_motion_parameters.grid_rotational_speed * dem_time_step; - rotation_axis = grid_motion_parameters.grid_rotational_axis; - } - else if (grid_motion_parameters.motion_type == - Parameters::Lagrangian::GridMotion< - spacedim>::MotionType::translational) + auto action_manager = DEMActionManager::get_action_manager(); + + switch (grid_motion_parameters.motion_type) { - grid_motion = &GridMotion::move_grid_translational; - shift_vector = - grid_motion_parameters.grid_translational_velocity * dem_time_step; + case Parameters::Lagrangian::GridMotion::MotionType::none: + { + grid_motion = &GridMotion::no_motion; + break; + } + case Parameters::Lagrangian::GridMotion::MotionType::rotational: + { + // Communicate to the action manager that there is a grid motion + action_manager->set_grid_motion_enabled(); + + grid_motion = &GridMotion::move_grid_rotational; + rotation_angle = + grid_motion_parameters.grid_rotational_speed * dem_time_step; + rotation_axis = grid_motion_parameters.grid_rotational_axis; + break; + } + case Parameters::Lagrangian::GridMotion< + spacedim>::MotionType::translational: + { + // Flag and communicate to the action manager that there is a + // grid motion + action_manager->set_grid_motion_enabled(); + + grid_motion = &GridMotion::move_grid_translational; + shift_vector = + grid_motion_parameters.grid_translational_velocity * dem_time_step; + break; + } + default: + AssertThrow(false, ExcNotImplemented()); } } @@ -83,6 +102,10 @@ GridMotion:: spacedim>::boundary_points_and_normal_vectors &updated_boundary_points_and_normal_vectors) { + // If there is no grid motion, exit the function + if (!DEMActionManager::get_action_manager()->check_grid_motion_enabled()) + return; + for (auto &[particle_id, pairs_in_contact_content] : particle_wall_pairs_in_contact) { diff --git a/source/dem/input_parameter_inspection.cc b/source/dem/input_parameter_inspection.cc index f5f4d89915..59fe7b5d93 100644 --- a/source/dem/input_parameter_inspection.cc +++ b/source/dem/input_parameter_inspection.cc @@ -12,7 +12,7 @@ input_parameter_inspection(const DEMSolverParameters &dem_parameters, // Getting the input parameters as local variable auto parameters = dem_parameters; auto physical_properties = dem_parameters.lagrangian_physical_properties; - double rayleigh_time_step = 1. / DBL_MIN; + double rayleigh_time_step = DBL_MAX; for (unsigned int i = 0; i < physical_properties.particle_type_number; ++i) { double shear_modulus = diff --git a/source/dem/load_balancing.cc b/source/dem/load_balancing.cc index 6a061a794b..3c4b4acc0d 100644 --- a/source/dem/load_balancing.cc +++ b/source/dem/load_balancing.cc @@ -1,4 +1,5 @@ #include +#include #include using namespace dealii; @@ -11,31 +12,23 @@ LagrangianLoadBalancing::LagrangianLoadBalancing() {} template -inline bool +inline void LagrangianLoadBalancing::check_load_balance_once() { if (simulation_control->get_step_number() == load_balance_step) - { - return true; - } - - return false; + DEMActionManager::get_action_manager()->load_balance_step(); } template -inline bool +inline void LagrangianLoadBalancing::check_load_balance_frequent() { if ((simulation_control->get_step_number() % load_balance_frequency) == 0) - { - return true; - } - - return false; + DEMActionManager::get_action_manager()->load_balance_step(); } template -inline bool +inline void LagrangianLoadBalancing::check_load_balance_dynamic() { if (simulation_control->get_step_number() % dynamic_check_frequency == 0) @@ -46,22 +39,19 @@ LagrangianLoadBalancing::check_load_balance_dynamic() unsigned int minimum_particle_number_on_proc = Utilities::MPI::min(particle_handler->n_locally_owned_particles(), mpi_communicator); + unsigned int average_particle_number_on_proc = + particle_handler->n_locally_owned_particles() / n_mpi_processes; // Execute load balancing if difference of load between processors is // larger than threshold of the load per processor if ((maximum_particle_number_on_proc - minimum_particle_number_on_proc) > - load_threshold * - (particle_handler->n_global_particles() / n_mpi_processes)) - { - return true; - } + load_threshold * average_particle_number_on_proc) + DEMActionManager::get_action_manager()->load_balance_step(); } - - return false; } template -inline bool +inline void LagrangianLoadBalancing::check_load_balance_with_sparse_contacts() { if (simulation_control->get_step_number() % dynamic_check_frequency == 0) @@ -126,14 +116,12 @@ LagrangianLoadBalancing::check_load_balance_with_sparse_contacts() // Clear and connect a new cell weight function connect_mobility_status_weight_signals(); - return true; + DEMActionManager::get_action_manager()->load_balance_step(); } } - // Clear and connect a new cell weight function + // Clear and connect a new cell weight function with new mobility status connect_mobility_status_weight_signals(); - - return false; } #if (DEAL_II_VERSION_MAJOR < 10 && DEAL_II_VERSION_MINOR < 6) diff --git a/source/dem/particle_wall_contact_force.cc b/source/dem/particle_wall_contact_force.cc index 88fa23c812..85f4e1180b 100644 --- a/source/dem/particle_wall_contact_force.cc +++ b/source/dem/particle_wall_contact_force.cc @@ -157,7 +157,7 @@ ParticleWallContactForce::calculate_force_and_torque_on_boundary( Tensor<1, 3> add_force, const Point<3> point_contact) { - if (calculate_force_torque_on_boundary == true) + if (calculate_force_torque_on_boundary) { mpi_correction_over_calculation_of_forces_and_torques(); diff --git a/source/dem/periodic_boundaries_manipulator.cc b/source/dem/periodic_boundaries_manipulator.cc index 9412391bfa..fcfb940eaa 100644 --- a/source/dem/periodic_boundaries_manipulator.cc +++ b/source/dem/periodic_boundaries_manipulator.cc @@ -9,12 +9,11 @@ #include #include -#include - using namespace dealii; template PeriodicBoundariesManipulator::PeriodicBoundariesManipulator() + : periodic_boundaries_enabled(false) {} template @@ -64,6 +63,9 @@ PeriodicBoundariesManipulator::map_periodic_cells( typename DEM::dem_data_structures::periodic_boundaries_cells_info &periodic_boundaries_cells_information) { + if (!periodic_boundaries_enabled) + return; + periodic_boundaries_cells_information.clear(); // Iterating over the active cells in the triangulation @@ -122,9 +124,10 @@ template void PeriodicBoundariesManipulator::check_and_move_particles( const periodic_boundaries_cells_info_struct &boundaries_cells_content, + const bool &particles_in_pb0_cell, typename Particles::ParticleHandler::particle_iterator_range &particles_in_cell, - bool &particles_in_pb0_cell) + bool &particle_has_been_moved) { for (auto particle = particles_in_cell.begin(); particle != particles_in_cell.end(); @@ -137,6 +140,7 @@ PeriodicBoundariesManipulator::check_and_move_particles( // the particle Point point_on_face, point_on_periodic_face; Tensor<1, dim> normal_vector, distance_between_faces; + if (particles_in_pb0_cell) { point_on_face = boundaries_cells_content.point_on_face; @@ -163,17 +167,25 @@ PeriodicBoundariesManipulator::check_and_move_particles( // Move particle outside the current cell to the periodic cell. particle_position += distance_between_faces; particle->set_location(particle_position); + + // Update flag to indicate that particle has been moved + particle_has_been_moved = true; } } } template -void +bool PeriodicBoundariesManipulator::execute_particles_displacement( const Particles::ParticleHandler &particle_handler, const typename DEM::dem_data_structures::periodic_boundaries_cells_info &periodic_boundaries_cells_information) { + if (!periodic_boundaries_enabled) + return false; + + bool particle_has_been_moved = false; + if (!periodic_boundaries_cells_information.empty()) { for (auto boundaries_cells_information_iterator = @@ -202,8 +214,9 @@ PeriodicBoundariesManipulator::execute_particles_displacement( { bool particles_in_pb0_cell = true; check_and_move_particles(boundaries_cells_content, + particles_in_pb0_cell, particles_in_cell, - particles_in_pb0_cell); + particle_has_been_moved); } } @@ -220,12 +233,15 @@ PeriodicBoundariesManipulator::execute_particles_displacement( { bool particles_in_pb0_cell = false; check_and_move_particles(boundaries_cells_content, + particles_in_pb0_cell, particles_in_periodic_cell, - particles_in_pb0_cell); + particle_has_been_moved); } } } } + + return particle_has_been_moved; } template class PeriodicBoundariesManipulator<2>; diff --git a/source/dem/read_checkpoint.cc b/source/dem/read_checkpoint.cc index 6b9feca94a..46041e5bc9 100644 --- a/source/dem/read_checkpoint.cc +++ b/source/dem/read_checkpoint.cc @@ -1,3 +1,4 @@ +#include #include using namespace dealii; @@ -15,6 +16,9 @@ read_checkpoint( std::shared_ptr> &insertion_object, std::vector>> &solid_surfaces) { + if (!DEMActionManager::get_action_manager()->check_restart_simulation()) + return; + TimerOutput::Scope timer(computing_timer, "Read checkpoint"); std::string prefix = parameters.restart.filename; simulation_control->read(prefix); diff --git a/source/dem/update_local_particle_containers.cc b/source/dem/update_local_particle_containers.cc index 591db78c92..518b3f25a0 100644 --- a/source/dem/update_local_particle_containers.cc +++ b/source/dem/update_local_particle_containers.cc @@ -1,3 +1,4 @@ +#include #include using namespace dealii; @@ -36,9 +37,13 @@ void update_contact_container_iterators( pairs_structure &pairs_in_contact, const typename DEM::dem_data_structures::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures) + &particle_container) { + // Get trigger for clearing contact structures + auto action_manager = DEMActionManager::get_action_manager(); + const bool clear_tangential_overlap = + action_manager->check_clear_tangential_overlap(); + // Loop over particle-object (particle/wall/line/point) pairs in contact for (auto pairs_in_contact_iterator = pairs_in_contact.begin(); pairs_in_contact_iterator != pairs_in_contact.end();) @@ -110,7 +115,7 @@ update_contact_container_iterators( adjacent_map_iterator->second.particle_two = particle_two_container->second; - if (clear_contact_structures) + if (clear_tangential_overlap) { adjacent_map_iterator->second.tangential_overlap.clear(); } @@ -197,8 +202,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<2>::adjacent_particle_pairs &pairs_in_contact, const typename DEM::dem_data_structures<2>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); template void update_contact_container_iterators< @@ -208,8 +212,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<3>::adjacent_particle_pairs &pairs_in_contact, const typename DEM::dem_data_structures<3>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); // Ghost particle-particle contact container template void @@ -220,8 +223,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<2>::adjacent_particle_pairs &pairs_in_contact, const typename DEM::dem_data_structures<2>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); template void update_contact_container_iterators< @@ -231,8 +233,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<3>::adjacent_particle_pairs &pairs_in_contact, const typename DEM::dem_data_structures<3>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); // Local-local particle-particle periodic contact container template void @@ -243,8 +244,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<2>::adjacent_particle_pairs &pairs_in_contact, const typename DEM::dem_data_structures<2>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); template void update_contact_container_iterators< @@ -254,8 +254,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<3>::adjacent_particle_pairs &pairs_in_contact, const typename DEM::dem_data_structures<3>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); // Local-ghost particle-particle periodic contact container template void @@ -266,8 +265,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<2>::adjacent_particle_pairs &pairs_in_contact, const typename DEM::dem_data_structures<2>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); template void update_contact_container_iterators< @@ -277,8 +275,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<3>::adjacent_particle_pairs &pairs_in_contact, const typename DEM::dem_data_structures<3>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); // Ghost-local particle-particle periodic contact container template void @@ -289,8 +286,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<2>::adjacent_particle_pairs &pairs_in_contact, const typename DEM::dem_data_structures<2>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); template void update_contact_container_iterators< @@ -300,8 +296,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<3>::adjacent_particle_pairs &pairs_in_contact, const typename DEM::dem_data_structures<3>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); // Particle-wall contact container template void @@ -312,8 +307,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<2>::particle_wall_in_contact &pairs_in_contact, const typename DEM::dem_data_structures<2>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); template void update_contact_container_iterators< @@ -323,8 +317,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<3>::particle_wall_in_contact &pairs_in_contact, const typename DEM::dem_data_structures<3>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); // Particle-floating wall contact container template void @@ -335,8 +328,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<2>::particle_wall_in_contact &pairs_in_contact, const typename DEM::dem_data_structures<2>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); template void update_contact_container_iterators< @@ -346,8 +338,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<3>::particle_wall_in_contact &pairs_in_contact, const typename DEM::dem_data_structures<3>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); // Particle-floating mesh contact container template void @@ -359,8 +350,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures< 2>::particle_floating_wall_from_mesh_in_contact &pairs_in_contact, const typename DEM::dem_data_structures<2>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); template void update_contact_container_iterators< @@ -371,8 +361,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures< 3>::particle_floating_wall_from_mesh_in_contact &pairs_in_contact, const typename DEM::dem_data_structures<3>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); // Particle-point contact container template void @@ -383,8 +372,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<2>::particle_point_line_contact_info &pairs_in_contact, const typename DEM::dem_data_structures<2>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); template void update_contact_container_iterators< @@ -394,8 +382,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<3>::particle_point_line_contact_info &pairs_in_contact, const typename DEM::dem_data_structures<3>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); // Particle-line contact container template void @@ -406,8 +393,7 @@ update_contact_container_iterators< typename DEM::dem_data_structures<2>::particle_point_line_contact_info &pairs_in_contact, const typename DEM::dem_data_structures<2>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); template void update_contact_container_iterators< @@ -417,5 +403,4 @@ update_contact_container_iterators< typename DEM::dem_data_structures<3>::particle_point_line_contact_info &pairs_in_contact, const typename DEM::dem_data_structures<3>::particle_index_iterator_map - &particle_container, - const bool clear_contact_structures); + &particle_container); diff --git a/source/dem/velocity_verlet_integrator.cc b/source/dem/velocity_verlet_integrator.cc index d3f49c39f8..24a7abdf12 100644 --- a/source/dem/velocity_verlet_integrator.cc +++ b/source/dem/velocity_verlet_integrator.cc @@ -143,6 +143,20 @@ VelocityVerletIntegrator::integrate( const parallel::distributed::Triangulation &triangulation, AdaptiveSparseContacts &sparse_contacts_object) { + auto action_manager = DEMActionManager::get_action_manager(); + + bool use_default_function = + !action_manager->check_sparse_contacts_enabled() || + action_manager->check_mobility_status_reset(); + + // Regular integration process if sparse contacts are not enabled or if this + // is the first iteration. + if (use_default_function) + { + integrate(particle_handler, g, dt, torque, force, MOI); + return; + } + // If there are advected particles, we use another function since the average // velocity and acceleration of cells are computed for mobile cells if (sparse_contacts_object.has_advected_particles()) diff --git a/source/fem-dem/cfd_dem_coupling.cc b/source/fem-dem/cfd_dem_coupling.cc index 53e75eea65..42fe80a5d9 100644 --- a/source/fem-dem/cfd_dem_coupling.cc +++ b/source/fem-dem/cfd_dem_coupling.cc @@ -12,73 +12,263 @@ #include #include +// Constructor for the class CFD-DEM class template -bool -check_contact_detection_method( - unsigned int counter, - CFDDEMSimulationParameters ¶m, - std::vector &displacement, - Particles::ParticleHandler &particle_handler, - MPI_Comm mpi_communicator, - std::shared_ptr simulation_control, - bool contact_detection_step, - bool checkpoint_step, - bool load_balance_step, - double smallest_contact_search_criterion) +CFDDEMSolver::CFDDEMSolver(CFDDEMSimulationParameters &nsparam) + : GLSVANSSolver(nsparam) + , this_mpi_process(Utilities::MPI::this_mpi_process(this->mpi_communicator)) + , n_mpi_processes(Utilities::MPI::n_mpi_processes(this->mpi_communicator)) +{} + +template +CFDDEMSolver::~CFDDEMSolver() +{} + +template +void +CFDDEMSolver::dem_setup_parameters() { - if (param.dem_parameters.model_parameters.contact_detection_method == - Parameters::Lagrangian::ModelParameters::ContactDetectionMethod::constant) + coupling_frequency = + this->cfd_dem_simulation_parameters.cfd_dem.coupling_frequency; + + dem_action_manager = DEMActionManager::get_action_manager(); + + // Initialize DEM Parameters + dem_parameters.lagrangian_physical_properties = + this->cfd_dem_simulation_parameters.dem_parameters + .lagrangian_physical_properties; + g = dem_parameters.lagrangian_physical_properties.g; + dem_parameters.boundary_conditions = + this->cfd_dem_simulation_parameters.dem_parameters.boundary_conditions; + dem_parameters.floating_walls = + this->cfd_dem_simulation_parameters.dem_parameters.floating_walls; + dem_parameters.model_parameters = + this->cfd_dem_simulation_parameters.dem_parameters.model_parameters; + dem_parameters.simulation_control = + this->cfd_dem_simulation_parameters.dem_parameters.simulation_control; + dem_parameters.post_processing = + this->cfd_dem_simulation_parameters.dem_parameters.post_processing; + dem_parameters.mesh = this->cfd_dem_simulation_parameters.dem_parameters.mesh; + dem_parameters.restart = + this->cfd_dem_simulation_parameters.dem_parameters.restart; + size_distribution_object_container.resize( + dem_parameters.lagrangian_physical_properties.particle_type_number); + + const auto parallel_triangulation = + dynamic_cast *>( + &*this->triangulation); + + // Setup load balancing parameters + load_balancing.set_parameters(dem_parameters.model_parameters); + load_balancing.copy_references(this->simulation_control, + *parallel_triangulation, + this->particle_handler, + sparse_contacts_object); + + // Attach the correct functions to the signals inside the triangulation + load_balancing.connect_weight_signals(); + + if (dem_parameters.model_parameters.sparse_particle_contacts) + { + dem_action_manager->set_sparse_contacts_enabled(); + sparse_contacts_object.set_parameters( + dem_parameters.model_parameters.granular_temperature_threshold, + dem_parameters.model_parameters.solid_fraction_threshold, + dem_parameters.model_parameters.advect_particles); + } + + setup_distribution_type(); + + dem_time_step = + this->simulation_control->get_time_step() / coupling_frequency; + + // Calculate the Rayleigh critical time step ratio + double rayleigh_time_step = DBL_MAX; + + for (unsigned int i = 0; + i < dem_parameters.lagrangian_physical_properties.particle_type_number; + ++i) { - return ((counter % param.dem_parameters.model_parameters - .contact_detection_frequency) == 0); + double youngs_modulus = dem_parameters.lagrangian_physical_properties + .youngs_modulus_particle[i]; + double poisson_ratio = + dem_parameters.lagrangian_physical_properties.poisson_ratio_particle[i]; + double density = + dem_parameters.lagrangian_physical_properties.density_particle[i]; + + double shear_modulus = youngs_modulus / (2.0 * (1.0 + poisson_ratio)); + + double min_diameter = + size_distribution_object_container.at(i)->find_min_diameter(); + + rayleigh_time_step = + std::min(M_PI_2 * min_diameter * sqrt(density / shear_modulus) / + (0.1631 * poisson_ratio + 0.8766), + rayleigh_time_step); } - else if (param.dem_parameters.model_parameters.contact_detection_method == - Parameters::Lagrangian::ModelParameters::ContactDetectionMethod:: - dynamic) + + const double time_step_rayleigh_ratio = dem_time_step / rayleigh_time_step; + this->pcout << "DEM time-step is " << time_step_rayleigh_ratio * 100 + << "% of Rayleigh time step" << std::endl; + + // Check if there are periodic boundaries + for (unsigned int i_bc = 0; + i_bc < dem_parameters.boundary_conditions.bc_types.size(); + ++i_bc) { - // The sorting into subdomain step checks whether or not the current time - // step is a step that requires sorting particles into subdomains and - // cells. - // This is applicable if any of the following three conditions apply:if - // its a load balancing step, a restart simulation step, or a contact - // detection tsep. - bool sorting_in_subdomains_step = - (checkpoint_step || load_balance_step || contact_detection_step); - - if (sorting_in_subdomains_step) - displacement.resize(particle_handler.get_max_local_particle_index()); - - contact_detection_step = find_particle_contact_detection_step( - particle_handler, - simulation_control->get_time_step() / param.cfd_dem.coupling_frequency, - smallest_contact_search_criterion, - mpi_communicator, - sorting_in_subdomains_step, - displacement); - - return contact_detection_step; + if (dem_parameters.boundary_conditions.bc_types[i_bc] == + Parameters::Lagrangian::BCDEM::BoundaryType::periodic) + { + dem_action_manager->set_periodic_boundaries_enabled(); + + periodic_boundaries_object.set_periodic_boundaries_information( + dem_parameters.boundary_conditions.periodic_boundary_0, + dem_parameters.boundary_conditions.periodic_direction); + + break; + } } - else + + // Initialize the total contact list counter + integrator_object = set_integrator_type(); + particle_particle_contact_force_object = + set_particle_particle_contact_force_model( + this->cfd_dem_simulation_parameters.dem_parameters); + + // Initialize the contact search counter + contact_search_total_number = 0; +} + +template +void +CFDDEMSolver::setup_distribution_type() +{ + // Use namespace and alias to make the code more readable + using namespace Parameters::Lagrangian; + LagrangianPhysicalProperties &lpp = + dem_parameters.lagrangian_physical_properties; + + maximum_particle_diameter = 0; + for (unsigned int particle_type = 0; particle_type < lpp.particle_type_number; + particle_type++) { - throw std::runtime_error( - "Specified contact detection method is not valid"); + switch (lpp.distribution_type.at(particle_type)) + { + case SizeDistributionType::uniform: + size_distribution_object_container[particle_type] = + std::make_shared( + lpp.particle_average_diameter.at(particle_type)); + break; + case SizeDistributionType::normal: + size_distribution_object_container[particle_type] = + std::make_shared( + lpp.particle_average_diameter.at(particle_type), + lpp.particle_size_std.at(particle_type), + lpp.seed_for_distributions[particle_type] + this_mpi_process); + break; + case SizeDistributionType::custom: + size_distribution_object_container[particle_type] = + std::make_shared( + lpp.particle_custom_diameter.at(particle_type), + lpp.particle_custom_probability.at(particle_type), + lpp.seed_for_distributions[particle_type] + this_mpi_process); + break; + } + + maximum_particle_diameter = std::max( + maximum_particle_diameter, + size_distribution_object_container[particle_type]->find_max_diameter()); } + + neighborhood_threshold_squared = + std::pow(dem_parameters.model_parameters.neighborhood_threshold * + maximum_particle_diameter, + 2); } -// Constructor for class CFD-DEM template -CFDDEMSolver::CFDDEMSolver(CFDDEMSimulationParameters &nsparam) - : GLSVANSSolver(nsparam) - , has_periodic_boundaries(false) - , has_sparse_contacts(false) - , this_mpi_process(Utilities::MPI::this_mpi_process(this->mpi_communicator)) - , n_mpi_processes(Utilities::MPI::n_mpi_processes(this->mpi_communicator)) -{} +std::shared_ptr> +CFDDEMSolver::set_integrator_type() +{ + using namespace Parameters::Lagrangian; + ModelParameters::IntegrationMethod integration_method = + dem_parameters.model_parameters.integration_method; + + switch (integration_method) + { + case ModelParameters::IntegrationMethod::velocity_verlet: + return std::make_shared>(); + case ModelParameters::IntegrationMethod::explicit_euler: + return std::make_shared>(); + case ModelParameters::IntegrationMethod::gear3: + return std::make_shared>(); + default: + throw(std::runtime_error("Invalid integration method.")); + } +} template -CFDDEMSolver::~CFDDEMSolver() -{} +void +CFDDEMSolver::initialize_dem_parameters() +{ + this->pcout << "Initializing DEM parameters" << std::endl; + + const auto parallel_triangulation = + dynamic_cast *>( + &*this->triangulation); + + // Set up the local and ghost cells (if ASC enabled) + sparse_contacts_object.update_local_and_ghost_cell_set( + this->void_fraction_dof_handler); + + particle_wall_contact_force_object = set_particle_wall_contact_force_model( + this->cfd_dem_simulation_parameters.dem_parameters, + *parallel_triangulation); + + // Finding the smallest contact search frequency criterion between (smallest + // cell size - largest particle radius) and (security factor * (blob diameter + // - 1) * largest particle radius). This value is used in + // find_contact_detection_frequency function + smallest_contact_search_criterion = + std::min((GridTools::minimal_cell_diameter(*this->triangulation) - + maximum_particle_diameter * 0.5), + (dem_parameters.model_parameters.dynamic_contact_search_factor * + (dem_parameters.model_parameters.neighborhood_threshold - 1) * + maximum_particle_diameter * 0.5)); + + // Remap periodic cells (if PBC enabled) + periodic_boundaries_object.map_periodic_cells( + *parallel_triangulation, periodic_boundaries_cells_information); + periodic_offset = periodic_boundaries_object.get_periodic_offset_distance(); + + // Find cell neighbors + contact_manager.execute_cell_neighbors_search( + *parallel_triangulation, periodic_boundaries_cells_information); + + // Find boundary cells with faces + boundary_cell_object.build( + *parallel_triangulation, + dem_parameters.floating_walls, + dem_parameters.boundary_conditions.outlet_boundaries, + this->cfd_dem_simulation_parameters.cfd_parameters.mesh + .check_for_diamond_cells, + this->cfd_dem_simulation_parameters.cfd_parameters.mesh + .expand_particle_wall_contact_search, + this->pcout); + + sort_particles_into_subdomains_and_cells(); + + // Remap periodic nodes after setup of dofs (If ASC and PBC) + if (dem_action_manager->check_periodic_boundaries_enabled() && + dem_action_manager->check_sparse_contacts_enabled()) + { + sparse_contacts_object.map_periodic_nodes( + this->void_fraction_constraints); + } + this->pcout << "Finished initializing DEM parameters" << std::endl + << "DEM time-step is " << dem_time_step << " s" << std::endl; +} template void CFDDEMSolver::read_dem() @@ -138,7 +328,7 @@ CFDDEMSolver::read_dem() << this->particle_handler.n_global_particles() << " particles are in the simulation" << std::endl; - write_DEM_output_results(); + write_dem_output_results(); } template @@ -283,7 +473,8 @@ CFDDEMSolver::read_checkpoint() this->setup_dofs(); // Remap periodic nodes after setup of dofs - if (has_periodic_boundaries && has_sparse_contacts) + if (dem_action_manager->check_periodic_boundaries_enabled() && + dem_action_manager->check_sparse_contacts_enabled()) { sparse_contacts_object.map_periodic_nodes( this->void_fraction_constraints); @@ -371,10 +562,52 @@ CFDDEMSolver::read_checkpoint() this->particle_handler.deserialize(); } +template +void +CFDDEMSolver::check_contact_detection_method(unsigned int counter) +{ + // Use namespace and alias to make the code more readable + using namespace Parameters::Lagrangian; + Parameters::Lagrangian::ModelParameters &model_parameters = + dem_parameters.model_parameters; + + switch (model_parameters.contact_detection_method) + { + case ModelParameters::ContactDetectionMethod::constant: + { + if ((counter % model_parameters.contact_detection_frequency) == 0) + dem_action_manager->contact_detection_step(); + break; + } + case ModelParameters::ContactDetectionMethod::dynamic: + { + double dt = + this->simulation_control->get_time_step() / + this->cfd_dem_simulation_parameters.cfd_dem.coupling_frequency; + + find_particle_contact_detection_step( + this->particle_handler, + dt, + smallest_contact_search_criterion, + this->mpi_communicator, + displacement); + break; + } + default: + break; + } +} + template void CFDDEMSolver::load_balance() { + load_balancing.check_load_balance_iteration(); + + // If not a load balance iteration, exit the function + if (!dem_action_manager->check_load_balance()) + return; + std::vector sol_set_transfer; sol_set_transfer.push_back(&this->present_solution); for (unsigned int i = 0; i < this->previous_solutions.size(); ++i) @@ -387,7 +620,6 @@ CFDDEMSolver::load_balance() system_trans_vectors(this->dof_handler); system_trans_vectors.prepare_for_coarsening_and_refinement(sol_set_transfer); - // Void Fraction std::vector vf_set_transfer; vf_set_transfer.push_back(&this->nodal_void_fraction_relevant); @@ -413,20 +645,13 @@ CFDDEMSolver::load_balance() parallel_triangulation->repartition(); - // Update cell neighbors - if (has_periodic_boundaries) - { - periodic_boundaries_object.map_periodic_cells( - *parallel_triangulation, periodic_boundaries_cells_information); - - periodic_offset = - periodic_boundaries_object.get_periodic_offset_distance(); - } + // If PBC are enabled remap periodic cells + periodic_boundaries_object.map_periodic_cells( + *parallel_triangulation, periodic_boundaries_cells_information); + // Update cell neighbors contact_manager.update_cell_neighbors(*parallel_triangulation, - periodic_boundaries_cells_information, - has_periodic_boundaries); - + periodic_boundaries_cells_information); boundary_cell_object.build( *parallel_triangulation, @@ -460,12 +685,17 @@ CFDDEMSolver::load_balance() this->setup_dofs(); // Remap periodic nodes after setup of dofs - if (has_periodic_boundaries && has_sparse_contacts) + if (dem_action_manager->check_periodic_boundaries_enabled() && + dem_action_manager->check_sparse_contacts_enabled()) { sparse_contacts_object.map_periodic_nodes( this->void_fraction_constraints); } + // Update the local and ghost cells (if ASC enabled) + sparse_contacts_object.update_local_and_ghost_cell_set( + this->void_fraction_dof_handler); + // Velocity Vectors std::vector x_system(1 + this->previous_solutions.size()); @@ -536,412 +766,45 @@ CFDDEMSolver::load_balance() template void -CFDDEMSolver::initialize_dem_parameters() +CFDDEMSolver::add_fluid_particle_interaction_force() { - this->pcout << "Initializing DEM parameters" << std::endl; - - const auto parallel_triangulation = - dynamic_cast *>( - &*this->triangulation); - - if (has_periodic_boundaries) + for (auto particle = this->particle_handler.begin(); + particle != this->particle_handler.end(); + ++particle) { - periodic_boundaries_object.set_periodic_boundaries_information( - dem_parameters.boundary_conditions.periodic_boundary_0, - dem_parameters.boundary_conditions.periodic_boundary_1, - dem_parameters.boundary_conditions.periodic_direction); - - periodic_boundaries_object.map_periodic_cells( - *parallel_triangulation, periodic_boundaries_cells_information); + auto particle_properties = particle->get_properties(); - // Temporary offset calculation : works only for one set of periodic - // boundary on an axis. - periodic_offset = - periodic_boundaries_object.get_periodic_offset_distance(); + types::particle_index particle_id = particle->get_local_index(); - // Initialize the flag for particle displacement in PBC - particle_displaced_in_pbc = false; + force[particle_id][0] += + particle_properties[DEM::PropertiesIndex::fem_force_x]; + force[particle_id][1] += + particle_properties[DEM::PropertiesIndex::fem_force_y]; + force[particle_id][2] += + particle_properties[DEM::PropertiesIndex::fem_force_z]; } +} - // Setup load balancing parameters - load_balancing.set_parameters(dem_parameters.model_parameters); - load_balancing.copy_references(this->simulation_control, - *parallel_triangulation, - this->particle_handler, - sparse_contacts_object); +template +void +CFDDEMSolver::add_fluid_particle_interaction_torque() +{ + for (auto particle = this->particle_handler.begin(); + particle != this->particle_handler.end(); + ++particle) + { + auto particle_properties = particle->get_properties(); - // Attach the correct functions to the signals inside the triangulation - load_balancing.connect_weight_signals(); + types::particle_index particle_id = particle->get_local_index(); - if (dem_parameters.model_parameters.sparse_particle_contacts) - { - has_sparse_contacts = true; - sparse_contacts_object.set_parameters( - dem_parameters.model_parameters.granular_temperature_threshold, - dem_parameters.model_parameters.solid_fraction_threshold, - dem_parameters.model_parameters.advect_particles); + torque[particle_id][0] += + particle_properties[DEM::PropertiesIndex::fem_torque_x]; + torque[particle_id][1] += + particle_properties[DEM::PropertiesIndex::fem_torque_y]; + torque[particle_id][2] += + particle_properties[DEM::PropertiesIndex::fem_torque_z]; } - - // Finding cell neighbors - contact_manager.execute_cell_neighbors_search( - *parallel_triangulation, - periodic_boundaries_cells_information, - has_periodic_boundaries); - - // Finding boundary cells with faces - boundary_cell_object.build( - *parallel_triangulation, - dem_parameters.floating_walls, - dem_parameters.boundary_conditions.outlet_boundaries, - this->cfd_dem_simulation_parameters.cfd_parameters.mesh - .check_for_diamond_cells, - this->cfd_dem_simulation_parameters.cfd_parameters.mesh - .expand_particle_wall_contact_search, - this->pcout); - - // Setting chosen contact force, insertion and integration methods - integrator_object = set_integrator_type(); - particle_particle_contact_force_object = - set_particle_particle_contact_force_model( - this->cfd_dem_simulation_parameters.dem_parameters); - particle_wall_contact_force_object = set_particle_wall_contact_force_model( - this->cfd_dem_simulation_parameters.dem_parameters, - *parallel_triangulation); - - this->particle_handler.sort_particles_into_subdomains_and_cells(); - - displacement.resize(this->particle_handler.get_max_local_particle_index()); - - force.resize(displacement.size()); - torque.resize(displacement.size()); - - this->particle_handler.exchange_ghost_particles(true); - - // Updating moment of inertia container - update_moment_of_inertia(this->particle_handler, MOI); - - this->pcout << "Finished initializing DEM parameters" << std::endl - << "DEM time-step is " << dem_time_step << " s" << std::endl; -} - -template -void -CFDDEMSolver::update_moment_of_inertia( - dealii::Particles::ParticleHandler &particle_handler, - std::vector &MOI) -{ - MOI.resize(torque.size()); - - for (auto &particle : particle_handler) - { - auto particle_properties = particle.get_properties(); - MOI[particle.get_local_index()] = - 0.1 * particle_properties[DEM::PropertiesIndex::mass] * - particle_properties[DEM::PropertiesIndex::dp] * - particle_properties[DEM::PropertiesIndex::dp]; - } -} - -template -std::shared_ptr> -CFDDEMSolver::set_integrator_type() -{ - if (dem_parameters.model_parameters.integration_method == - Parameters::Lagrangian::ModelParameters::IntegrationMethod:: - velocity_verlet) - { - integrator_object = std::make_shared>(); - } - else if (dem_parameters.model_parameters.integration_method == - Parameters::Lagrangian::ModelParameters::IntegrationMethod:: - explicit_euler) - { - integrator_object = std::make_shared>(); - } - else if (dem_parameters.model_parameters.integration_method == - Parameters::Lagrangian::ModelParameters::IntegrationMethod::gear3) - { - integrator_object = std::make_shared>(); - } - else - { - throw "The chosen integration method is invalid"; - } - return integrator_object; -} - -template -void -CFDDEMSolver::add_fluid_particle_interaction_force() -{ - for (auto particle = this->particle_handler.begin(); - particle != this->particle_handler.end(); - ++particle) - { - auto particle_properties = particle->get_properties(); - - types::particle_index particle_id = particle->get_local_index(); - - force[particle_id][0] += - particle_properties[DEM::PropertiesIndex::fem_force_x]; - force[particle_id][1] += - particle_properties[DEM::PropertiesIndex::fem_force_y]; - force[particle_id][2] += - particle_properties[DEM::PropertiesIndex::fem_force_z]; - } -} - -template -void -CFDDEMSolver::add_fluid_particle_interaction_torque() -{ - for (auto particle = this->particle_handler.begin(); - particle != this->particle_handler.end(); - ++particle) - { - auto particle_properties = particle->get_properties(); - - types::particle_index particle_id = particle->get_local_index(); - - torque[particle_id][0] += - particle_properties[DEM::PropertiesIndex::fem_torque_x]; - torque[particle_id][1] += - particle_properties[DEM::PropertiesIndex::fem_torque_y]; - torque[particle_id][2] += - particle_properties[DEM::PropertiesIndex::fem_torque_z]; - } -} - -template -void -CFDDEMSolver::dem_iterator(unsigned int counter) -{ - // dem_contact_build carries out the particle-particle and particle-wall - // broad and fine searches, sort_particles_into_subdomains_and_cells, and - // exchange_ghost - dem_contact_build(counter); - - // Particle-particle contact force - particle_particle_contact_force_object - ->calculate_particle_particle_contact_force( - contact_manager, dem_time_step, torque, force, periodic_offset); - - // Particles-walls contact force: - particle_wall_contact_force(); - - // Add fluid-particle interaction force to the force container - add_fluid_particle_interaction_force(); - - // Add fluid-particle interaction torque to the torque container - if (this->cfd_dem_simulation_parameters.cfd_dem.rotational_viscous_torque || - this->cfd_dem_simulation_parameters.cfd_dem.vortical_viscous_torque) - add_fluid_particle_interaction_torque(); - - // Integration correction step (after force calculation) - // In the first step, we have to obtain location of particles at half-step - // time - if (this->simulation_control->get_step_number() == 0) - { - integrator_object->integrate_half_step_location( - this->particle_handler, g, dem_time_step, torque, force, MOI); - } - else - { - if (!has_sparse_contacts) - { - integrator_object->integrate( - this->particle_handler, g, dem_time_step, torque, force, MOI); - } - else - { - if (counter > 0) - { - const auto parallel_triangulation = - dynamic_cast *>( - &*this->triangulation); - integrator_object->integrate(this->particle_handler, - g, - dem_time_step, - torque, - force, - MOI, - *parallel_triangulation, - sparse_contacts_object); - } - else // counter == 0 - { - // The cell average velocities and accelerations are updated - // from the fully computed forces at step 0, but we do not use the - // mobility status to disabled contacts at this step. The reason - // is that the current mobility status are computed for the last - // DEM time step (from the last CFD time step) are the force of - // the fluid may have significantly changed the particulate - // agitation. - - // Update the cell average velocities and accelerations - sparse_contacts_object.update_average_velocities_acceleration( - this->particle_handler, g, force, dem_time_step); - - integrator_object->integrate( - this->particle_handler, g, dem_time_step, torque, force, MOI); - } - } - } -} - -template -void -CFDDEMSolver::dem_contact_build(unsigned int counter) -{ - // Check to see if it is contact detection step - contact_detection_step = - check_contact_detection_method(counter, - this->cfd_dem_simulation_parameters, - displacement, - this->particle_handler, - this->mpi_communicator, - this->simulation_control, - contact_detection_step, - checkpoint_step, - load_balance_step, - smallest_contact_search_criterion); - - // Sort particles in cells - if (contact_search_step(counter)) - { - this->pcout << "DEM contact search at dem step " << counter << std::endl; - contact_search_counter++; - - periodic_boundaries_object.execute_particles_displacement( - this->particle_handler, periodic_boundaries_cells_information); - - this->particle_handler.sort_particles_into_subdomains_and_cells(); - - displacement.resize( - this->particle_handler.get_max_local_particle_index()); - force.resize(displacement.size()); - torque.resize(displacement.size()); - - this->particle_handler.exchange_ghost_particles(true); - - if (has_sparse_contacts) - { - if (load_balance_step || checkpoint_step || - (this->simulation_control->is_at_start() && (counter == 0))) - sparse_contacts_object.update_local_and_ghost_cell_set( - this->void_fraction_dof_handler); - - sparse_contacts_object.identify_mobility_status( - this->void_fraction_dof_handler, - this->particle_handler, - (*this->triangulation).n_active_cells(), - this->mpi_communicator, - counter); - } - - // Updating moment of inertia container - update_moment_of_inertia(this->particle_handler, MOI); - - // Execute broad search by filling containers of particle-particle - // contact pair candidates and containers of particle-wall - // contact pair candidates - if (!(has_sparse_contacts && counter > 0)) - { - contact_manager.execute_particle_particle_broad_search( - this->particle_handler, has_periodic_boundaries); - - contact_manager.execute_particle_wall_broad_search( - this->particle_handler, - boundary_cell_object, - solid_surfaces_mesh_info, - dem_parameters.floating_walls, - this->simulation_control->get_current_time()); - } - else - { - contact_manager.execute_particle_particle_broad_search( - this->particle_handler, - sparse_contacts_object, - has_periodic_boundaries); - - contact_manager.execute_particle_wall_broad_search( - this->particle_handler, - boundary_cell_object, - solid_surfaces_mesh_info, - dem_parameters.floating_walls, - this->simulation_control->get_current_time(), - sparse_contacts_object); - } - - // Update contacts, remove replicates and add new contact pairs - // to the contact containers when particles are exchanged between - // processors - contact_manager.update_contacts(has_periodic_boundaries); - - // Updates the iterators to particles in local-local contact - // containers - contact_manager.update_local_particles_in_cells(this->particle_handler, - load_balance_step, - has_periodic_boundaries); - - // Execute fine search by updating particle-particle contact - // containers regards the neighborhood threshold - contact_manager.execute_particle_particle_fine_search( - neighborhood_threshold_squared, - has_periodic_boundaries, - periodic_offset); - - // Execute fine search by updating particle-wall contact containers - // regards the neighborhood threshold - contact_manager.execute_particle_wall_fine_search( - dem_parameters.floating_walls, - this->simulation_control->get_current_time(), - neighborhood_threshold_squared); - - // Reset different steps. The contact build should be performed everytime - // we restart the simulation or everytime load balancing is performed. At - // the end of the restart step or the load balance step, and after all - // necessary contact build, vector resizing, and solution transfers have - // been performed, this functions are set to false. The checkpoint_step - // remains false for the duration of the simulation while the load - // balancing is reset everytime load balancing is called. - checkpoint_step = false; - load_balance_step = false; - particle_displaced_in_pbc = false; - } - else - { - this->particle_handler.update_ghost_particles(); - } - // TODO add DEM post-processing -} - -template -void -CFDDEMSolver::write_DEM_output_results() -{ - const std::string folder = dem_parameters.simulation_control.output_folder; - const std::string particles_solution_name = - dem_parameters.simulation_control.output_name + "_particles"; - const unsigned int iter = this->simulation_control->get_step_number(); - const double time = this->simulation_control->get_current_time(); - const unsigned int group_files = - dem_parameters.simulation_control.group_files; - - // Write particles - Visualization particle_data_out; - particle_data_out.build_patches(this->particle_handler, - properties_class.get_properties_name()); - - write_vtu_and_pvd<0, dim>(particles_pvdhandler, - particle_data_out, - folder, - particles_solution_name, - time, - iter, - group_files, - this->mpi_communicator); -} +} template void @@ -987,10 +850,36 @@ CFDDEMSolver::particle_wall_contact_force() } } +template +void +CFDDEMSolver::write_dem_output_results() +{ + const std::string folder = dem_parameters.simulation_control.output_folder; + const std::string particles_solution_name = + dem_parameters.simulation_control.output_name + "_particles"; + const unsigned int iter = this->simulation_control->get_step_number(); + const double time = this->simulation_control->get_current_time(); + const unsigned int group_files = + dem_parameters.simulation_control.group_files; + + // Write particles + Visualization particle_data_out; + particle_data_out.build_patches(this->particle_handler, + properties_class.get_properties_name()); + + write_vtu_and_pvd<0, dim>(particles_pvdhandler, + particle_data_out, + folder, + particles_solution_name, + time, + iter, + group_files, + this->mpi_communicator); +} template void -CFDDEMSolver::dem_post_process_results() +CFDDEMSolver::report_particle_statistics() { // Update statistics on contact list double number_of_list_built_since_last_log = @@ -1050,28 +939,82 @@ CFDDEMSolver::dem_post_process_results() report.set_scientific("Average", true); report.set_scientific("Total", true); - announce_string(this->pcout, "Particle statistics"); - report.write_text(std::cout, dealii::TableHandler::org_mode_table); + announce_string(this->pcout, "Particle statistics"); + report.write_text(std::cout, dealii::TableHandler::org_mode_table); + } + + if (dem_parameters.post_processing.Lagrangian_post_processing && + this->simulation_control->is_output_iteration()) + { + const auto parallel_triangulation = + dynamic_cast *>( + &*this->triangulation); + + + dem_post_processing_object.write_post_processing_results( + *parallel_triangulation, + grid_pvdhandler, + this->dof_handler, + this->particle_handler, + dem_parameters, + this->simulation_control->get_current_time(), + this->simulation_control->get_step_number(), + this->mpi_communicator, + sparse_contacts_object); + } +} + +template +void +CFDDEMSolver::print_particles_summary() +{ + this->pcout << "Particle Summary" << std::endl; + this->pcout << "id, x, y, z, v_x, v_y, v_z" << std::endl; + // Agressively force synchronization of the header line + usleep(500); + MPI_Barrier(this->mpi_communicator); + usleep(500); + MPI_Barrier(this->mpi_communicator); + + std::map> global_particles; + unsigned int current_id, current_id_max = 0; + + // Mapping of all particles & find the max id on current processor + for (auto particle = this->particle_handler.begin(); + particle != this->particle_handler.end(); + ++particle) + { + current_id = particle->get_id(); + current_id_max = std::max(current_id, current_id_max); + + global_particles.insert({current_id, particle}); } - if (dem_parameters.post_processing.Lagrangian_post_processing && - this->simulation_control->is_output_iteration()) - { - const auto parallel_triangulation = - dynamic_cast *>( - &*this->triangulation); + // Find global max particle index + unsigned int id_max = + Utilities::MPI::max(current_id_max, this->mpi_communicator); + for (unsigned int i = 0; i <= id_max; i++) + { + for (auto &iterator : global_particles) + { + unsigned int id = iterator.first; + if (id == i) + { + auto particle = iterator.second; + auto particle_properties = particle->get_properties(); + auto particle_location = particle->get_location(); - dem_post_processing_object.write_post_processing_results( - *parallel_triangulation, - grid_pvdhandler, - this->dof_handler, - this->particle_handler, - dem_parameters, - this->simulation_control->get_current_time(), - this->simulation_control->get_step_number(), - this->mpi_communicator, - sparse_contacts_object); + std::cout << std::setprecision(6) << id << " " + << particle_location << " " + << particle_properties[DEM::PropertiesIndex::v_x] << " " + << particle_properties[DEM::PropertiesIndex::v_y] << " " + << particle_properties[DEM::PropertiesIndex::v_z] + << std::endl; + } + } + usleep(500); + MPI_Barrier(this->mpi_communicator); } } @@ -1088,7 +1031,7 @@ CFDDEMSolver::postprocess_fd(bool first_iteration) // Visualization if (this->simulation_control->is_output_iteration()) { - write_DEM_output_results(); + write_dem_output_results(); } } @@ -1217,198 +1160,174 @@ CFDDEMSolver::dynamic_flow_control() } template -void -CFDDEMSolver::print_particles_summary() +inline void +CFDDEMSolver::sort_particles_into_subdomains_and_cells() { - this->pcout << "Particle Summary" << std::endl; - this->pcout << "id, x, y, z, v_x, v_y, v_z" << std::endl; - // Agressively force synchronization of the header line - usleep(500); - MPI_Barrier(this->mpi_communicator); - usleep(500); - MPI_Barrier(this->mpi_communicator); - - std::map> global_particles; - unsigned int current_id, current_id_max = 0; + this->particle_handler.sort_particles_into_subdomains_and_cells(); - // Mapping of all particles & find the max id on current processor - for (auto particle = this->particle_handler.begin(); - particle != this->particle_handler.end(); - ++particle) + // Resize the displacement, force and torque containers only if the particles + // have changed subdomains + if (dem_action_manager->check_resize_containers()) { - current_id = particle->get_id(); - current_id_max = std::max(current_id, current_id_max); - - global_particles.insert({current_id, particle}); - } - - // Find global max particle index - unsigned int id_max = - Utilities::MPI::max(current_id_max, this->mpi_communicator); + // Resize and reinitialize displacement container + displacement.resize( + this->particle_handler.get_max_local_particle_index()); + // Resize and reinitialize displacement container + force.resize(displacement.size()); + torque.resize(displacement.size()); + MOI.resize(displacement.size()); - for (unsigned int i = 0; i <= id_max; i++) - { - for (auto &iterator : global_particles) + // Updating moment of inertia container + for (auto &particle : this->particle_handler) { - unsigned int id = iterator.first; - if (id == i) - { - auto particle = iterator.second; - auto particle_properties = particle->get_properties(); - auto particle_location = particle->get_location(); - - std::cout << std::setprecision(6) << id << " " - << particle_location << " " - << particle_properties[DEM::PropertiesIndex::v_x] << " " - << particle_properties[DEM::PropertiesIndex::v_y] << " " - << particle_properties[DEM::PropertiesIndex::v_z] - << std::endl; - } + auto particle_properties = particle.get_properties(); + MOI[particle.get_local_index()] = + 0.1 * particle_properties[DEM::PropertiesIndex::mass] * + particle_properties[DEM::PropertiesIndex::dp] * + particle_properties[DEM::PropertiesIndex::dp]; } - usleep(500); - MPI_Barrier(this->mpi_communicator); } + + // Always reset the displacement values since we are doing a search detection + std::fill(displacement.begin(), displacement.end(), 0.); + + this->particle_handler.exchange_ghost_particles(true); } template void -CFDDEMSolver::dem_setup_contact_parameters() +CFDDEMSolver::dem_iterator(unsigned int counter) { - coupling_frequency = - this->cfd_dem_simulation_parameters.cfd_dem.coupling_frequency; + // dem_contact_build carries out the particle-particle and particle-wall + // broad and fine searches, sort_particles_into_subdomains_and_cells, and + // exchange_ghost + dem_contact_build(counter); - // Initialize DEM Parameters - dem_parameters.lagrangian_physical_properties = - this->cfd_dem_simulation_parameters.dem_parameters - .lagrangian_physical_properties; - g = dem_parameters.lagrangian_physical_properties.g; - dem_parameters.boundary_conditions = - this->cfd_dem_simulation_parameters.dem_parameters.boundary_conditions; - dem_parameters.floating_walls = - this->cfd_dem_simulation_parameters.dem_parameters.floating_walls; - dem_parameters.model_parameters = - this->cfd_dem_simulation_parameters.dem_parameters.model_parameters; - dem_parameters.simulation_control = - this->cfd_dem_simulation_parameters.dem_parameters.simulation_control; - dem_parameters.post_processing = - this->cfd_dem_simulation_parameters.dem_parameters.post_processing; - dem_parameters.mesh = this->cfd_dem_simulation_parameters.dem_parameters.mesh; - dem_parameters.restart = - this->cfd_dem_simulation_parameters.dem_parameters.restart; - size_distribution_object_container.reserve( - dem_parameters.lagrangian_physical_properties.particle_type_number); + // Particle-particle contact force + particle_particle_contact_force_object + ->calculate_particle_particle_contact_force( + contact_manager, dem_time_step, torque, force, periodic_offset); - maximum_particle_diameter = 0.; - for (unsigned int particle_type = 0; - particle_type < - dem_parameters.lagrangian_physical_properties.particle_type_number; - particle_type++) - { - if (dem_parameters.lagrangian_physical_properties.distribution_type.at( - particle_type) == - Parameters::Lagrangian::SizeDistributionType::uniform) - { - size_distribution_object_container.push_back( - std::make_shared( - dem_parameters.lagrangian_physical_properties - .particle_average_diameter.at(particle_type))); - } - else if (dem_parameters.lagrangian_physical_properties.distribution_type - .at(particle_type) == - Parameters::Lagrangian::SizeDistributionType::normal) - { - size_distribution_object_container.push_back( - std::make_shared( - dem_parameters.lagrangian_physical_properties - .particle_average_diameter.at(particle_type), - dem_parameters.lagrangian_physical_properties.particle_size_std - .at(particle_type), - dem_parameters.lagrangian_physical_properties - .seed_for_distributions[particle_type])); - } - else if (dem_parameters.lagrangian_physical_properties.distribution_type - .at(particle_type) == - Parameters::Lagrangian::SizeDistributionType::custom) - { - size_distribution_object_container.push_back( - std::make_shared( - dem_parameters.lagrangian_physical_properties - .particle_custom_diameter.at(particle_type), - dem_parameters.lagrangian_physical_properties - .particle_custom_probability.at(particle_type), - dem_parameters.lagrangian_physical_properties - .seed_for_distributions[particle_type])); - } + // Particles-walls contact force: + particle_wall_contact_force(); - maximum_particle_diameter = std::max( - maximum_particle_diameter, - size_distribution_object_container[particle_type]->find_max_diameter()); - } + // Add fluid-particle interaction force to the force container + add_fluid_particle_interaction_force(); - neighborhood_threshold_squared = - std::pow(dem_parameters.model_parameters.neighborhood_threshold * - maximum_particle_diameter, - 2); + // Add fluid-particle interaction torque to the torque container + if (this->cfd_dem_simulation_parameters.cfd_dem.rotational_viscous_torque || + this->cfd_dem_simulation_parameters.cfd_dem.vortical_viscous_torque) + add_fluid_particle_interaction_torque(); + // The cell average velocities and accelerations are updated + // from the fully computed forces at step 0, but we do not use the + // mobility status to disabled contacts at this step. The reason + // is that the current mobility status are computed for the last + // DEM time step (from the last CFD time step) are the force of + // the fluid may have significantly changed the particulate + // agitation. + // Update the cell average velocities and accelerations + sparse_contacts_object.update_average_velocities_acceleration( + this->particle_handler, g, force, dem_time_step); - // Finding the smallest contact search frequency criterion between (smallest - // cell size - largest particle radius) and (security factor * (blob diameter - // - 1) * largest particle radius). This value is used in - // find_contact_detection_frequency function - smallest_contact_search_criterion = - std::min((GridTools::minimal_cell_diameter(*this->triangulation) - - maximum_particle_diameter * 0.5), - (dem_parameters.model_parameters.dynamic_contact_search_factor * - (dem_parameters.model_parameters.neighborhood_threshold - 1) * - maximum_particle_diameter * 0.5)); + // Integration correction step (after force calculation) + // In the first step, we have to obtain location of particles at half-step + // time + // TODO do all DEM time step at first CFD time step are half step? + if (this->simulation_control->get_step_number() == 0) + { + integrator_object->integrate_half_step_location( + this->particle_handler, g, dem_time_step, torque, force, MOI); + } + else + { + const auto parallel_triangulation = + dynamic_cast *>( + &*this->triangulation); + integrator_object->integrate(this->particle_handler, + g, + dem_time_step, + torque, + force, + MOI, + *parallel_triangulation, + sparse_contacts_object); + } - dem_time_step = - this->simulation_control->get_time_step() / coupling_frequency; + dem_action_manager->reset_triggers(); +} - double rayleigh_time_step = 1. / DBL_MIN; +template +void +CFDDEMSolver::dem_contact_build(unsigned int counter) +{ + // For the contact search at the last CFD iteration + if (counter == (coupling_frequency - 1)) + dem_action_manager->last_dem_of_cfddem_iteration_step(); - Parameters::Lagrangian::LagrangianPhysicalProperties &physical_properties = - dem_parameters.lagrangian_physical_properties; + // Check contact detection step by detection method + check_contact_detection_method(counter); - for (unsigned int i = 0; i < physical_properties.particle_type_number; ++i) + // Sort particles in cells + if (dem_action_manager->check_contact_search()) { - double shear_modulus = - physical_properties.youngs_modulus_particle[i] / - (2.0 * (1.0 + physical_properties.poisson_ratio_particle[i])); + this->pcout << "DEM contact search at dem step " << counter << std::endl; + contact_search_counter++; - double min_diameter = - size_distribution_object_container.at(i)->find_min_diameter(); + // Execute periodic boundaries (if PBC enabled) + periodic_boundaries_object.execute_particles_displacement( + this->particle_handler, periodic_boundaries_cells_information); - rayleigh_time_step = std::min( - M_PI_2 * min_diameter * - sqrt(physical_properties.density_particle[i] / shear_modulus) / - (0.1631 * physical_properties.poisson_ratio_particle[i] + 0.8766), - rayleigh_time_step); - } + // Sort particles into subdomains and cells and reset the vectors that + // the local particle if is their index + sort_particles_into_subdomains_and_cells(); - const double time_step_rayleigh_ratio = dem_time_step / rayleigh_time_step; - this->pcout << "DEM time-step is " << time_step_rayleigh_ratio * 100 - << "% of Rayleigh time step" << std::endl; + // Identify the mobility status of particles (if ASC enabled) + sparse_contacts_object.identify_mobility_status( + this->void_fraction_dof_handler, + this->particle_handler, + (*this->triangulation).n_active_cells(), + this->mpi_communicator); - // Initialize contact detection step - contact_detection_step = false; - load_balance_step = false; + // Execute broad search by filling containers of particle-particle + // contact pair candidates and containers of particle-wall + // contact pair candidates + contact_manager.execute_particle_particle_broad_search( + this->particle_handler, sparse_contacts_object); - // Check if there are periodic boundaries - for (unsigned int i_bc = 0; - i_bc < dem_parameters.boundary_conditions.bc_types.size(); - ++i_bc) + contact_manager.execute_particle_wall_broad_search( + this->particle_handler, + boundary_cell_object, + solid_surfaces_mesh_info, + dem_parameters.floating_walls, + this->simulation_control->get_current_time(), + sparse_contacts_object); + + // Update contacts, remove replicates and add new contact pairs + // to the contact containers when particles are exchanged between + // processors + contact_manager.update_contacts(); + + // Updates the iterators to particles in local-local contact + // containers + contact_manager.update_local_particles_in_cells(this->particle_handler); + + // Execute fine search by updating particle-particle contact + // containers regards the neighborhood threshold + contact_manager.execute_particle_particle_fine_search( + neighborhood_threshold_squared, periodic_offset); + + // Execute fine search by updating particle-wall contact containers + // regards the neighborhood threshold + contact_manager.execute_particle_wall_fine_search( + dem_parameters.floating_walls, + this->simulation_control->get_current_time(), + neighborhood_threshold_squared); + } + else { - if (dem_parameters.boundary_conditions.bc_types[i_bc] == - Parameters::Lagrangian::BCDEM::BoundaryType::periodic) - { - has_periodic_boundaries = true; - break; - } + this->particle_handler.update_ghost_particles(); } - - // Initialize the total contact list counter - contact_search_total_number = 0; } template @@ -1422,12 +1341,12 @@ CFDDEMSolver::solve() true, this->cfd_dem_simulation_parameters.cfd_parameters.boundary_conditions); - dem_setup_contact_parameters(); + dem_setup_parameters(); // Reading DEM start file information if (this->cfd_dem_simulation_parameters.void_fraction->read_dem == true && - this->cfd_dem_simulation_parameters.cfd_parameters.restart_parameters - .restart == false) + !this->cfd_dem_simulation_parameters.cfd_parameters.restart_parameters + .restart) read_dem(); this->setup_dofs(); @@ -1438,32 +1357,20 @@ CFDDEMSolver::solve() .restart); // In the case the simulation is being restarted from a checkpoint file, the - // checkpoint_step parameter is set to true. This allows to perform all + // restart_simulation parameter is set to true. This allows to perform all // operations related to restarting a simulation. Once all operations have - // been performed, this checkpoint_step is reset to false. It is only set once - // and reset once since restarting only occurs once. + // been performed, this restart_simulation is reset to false. It is only set + // once and reset once since restarting only occurs once. if (this->cfd_dem_simulation_parameters.cfd_parameters.restart_parameters - .restart == true) - checkpoint_step = true; - else - checkpoint_step = false; + .restart) + dem_action_manager->restart_simulation(); // Initialize the DEM parameters and generate the required ghost particles initialize_dem_parameters(); - this->particle_handler.exchange_ghost_particles(true); - - - // Remap periodic nodes after setup of dofs - if (has_periodic_boundaries && has_sparse_contacts) - { - sparse_contacts_object.map_periodic_nodes( - this->void_fraction_constraints); - } - // Calculate first instance of void fraction once particles are set-up this->vertices_cell_mapping(); - if (!checkpoint_step) + if (!dem_action_manager->check_restart_simulation()) this->initialize_void_fraction(); while (this->simulation_control->integrate()) @@ -1497,7 +1404,8 @@ CFDDEMSolver::solve() } this->calculate_void_fraction( - this->simulation_control->get_current_time(), load_balance_step); + this->simulation_control->get_current_time(), + dem_action_manager->check_load_balance()); this->iterate(); if (this->cfd_dem_simulation_parameters.cfd_parameters.test.enabled) @@ -1538,16 +1446,11 @@ CFDDEMSolver::solve() announce_string(this->pcout, "DEM"); TimerOutput::Scope t(this->computing_timer, "DEM_Iterator"); - // Load balancing - // The input argument to this function is set to zero as this integer is - // not used for the check_load_balance_step function and is only - // important for the check_contact_search_step function. - load_balance_step = load_balancing.check_load_balance_iteration(); + // First DEM iteration of the CFD iteration + dem_action_manager->first_dem_of_cfddem_iteration_step(); - if (load_balance_step) - { - load_balance(); - } + // Load balancing if needed + load_balance(); contact_search_counter = 0; for (unsigned int dem_counter = 0; dem_counter < coupling_frequency; @@ -1559,6 +1462,7 @@ CFDDEMSolver::solve() dem_iterator(dem_counter); } } + contact_search_total_number += contact_search_counter; this->pcout << "Finished " << coupling_frequency << " DEM iterations" @@ -1571,7 +1475,7 @@ CFDDEMSolver::solve() this->GLSVANSSolver::monitor_mass_conservation(); if (this->cfd_dem_simulation_parameters.cfd_dem.particle_statistics) - dem_post_process_results(); + report_particle_statistics(); } this->finish_simulation(); diff --git a/tests/dem/find_contact_pairs.cc b/tests/dem/find_contact_pairs.cc index 3bd3b48206..0b74e954f7 100644 --- a/tests/dem/find_contact_pairs.cc +++ b/tests/dem/find_contact_pairs.cc @@ -103,8 +103,12 @@ test() Particles::ParticleIterator pit3 = particle_handler.insert_particle(particle3, pt3_info.first); + // Dummy Adaptive sparse contacts object for next call + AdaptiveSparseContacts dummy_adaptive_sparse_contacts; + // Calling broad search function - container_manager.execute_particle_particle_broad_search(particle_handler); + container_manager.execute_particle_particle_broad_search( + particle_handler, dummy_adaptive_sparse_contacts); // Output for (auto pairs_iterator = diff --git a/tests/dem/particle_particle_contact_force_linear.cc b/tests/dem/particle_particle_contact_force_linear.cc index 3e5fdeaa8f..d05265d8c3 100644 --- a/tests/dem/particle_particle_contact_force_linear.cc +++ b/tests/dem/particle_particle_contact_force_linear.cc @@ -157,7 +157,10 @@ test() particle_iterator; } - container_manager.execute_particle_particle_broad_search(particle_handler); + // Dummy Adaptive sparse contacts object and particle-particle broad search + AdaptiveSparseContacts dummy_adaptive_sparse_contacts; + container_manager.execute_particle_particle_broad_search( + particle_handler, dummy_adaptive_sparse_contacts); // Calling fine search container_manager.execute_particle_particle_fine_search( diff --git a/tests/dem/particle_particle_contact_force_nonlinear.cc b/tests/dem/particle_particle_contact_force_nonlinear.cc index f665da3931..6131b607d4 100644 --- a/tests/dem/particle_particle_contact_force_nonlinear.cc +++ b/tests/dem/particle_particle_contact_force_nonlinear.cc @@ -15,7 +15,7 @@ */ /** - * @ brief In this test, the performance of non-linear (Hertzian) + * @brief In this test, the performance of non-linear (Hertzian) * particle-particle contact force is checked. */ @@ -156,7 +156,10 @@ test() particle_iterator; } - container_manager.execute_particle_particle_broad_search(particle_handler); + // Dummy Adaptive sparse contacts object and particle-particle broad search + AdaptiveSparseContacts dummy_adaptive_sparse_contacts; + container_manager.execute_particle_particle_broad_search( + particle_handler, dummy_adaptive_sparse_contacts); // Calling fine search container_manager.execute_particle_particle_fine_search( diff --git a/tests/dem/particle_particle_contact_on_two_processors.cc b/tests/dem/particle_particle_contact_on_two_processors.cc index ffb579e9ec..59d688baa7 100644 --- a/tests/dem/particle_particle_contact_on_two_processors.cc +++ b/tests/dem/particle_particle_contact_on_two_processors.cc @@ -208,12 +208,13 @@ test() particle_handler.exchange_ghost_particles(); - container_manager.update_local_particles_in_cells(particle_handler, - false); + container_manager.update_local_particles_in_cells(particle_handler); - // Calling broad search + // Dummy Adaptive sparse contacts object and particle-particle broad + // search + AdaptiveSparseContacts dummy_adaptive_sparse_contacts; container_manager.execute_particle_particle_broad_search( - particle_handler); + particle_handler, dummy_adaptive_sparse_contacts); // Calling fine search container_manager.execute_particle_particle_fine_search( diff --git a/tests/dem/particle_particle_fine_search.cc b/tests/dem/particle_particle_fine_search.cc index f1822ae545..2a46f07823 100644 --- a/tests/dem/particle_particle_fine_search.cc +++ b/tests/dem/particle_particle_fine_search.cc @@ -123,7 +123,10 @@ test() particle_iterator; } - container_manager.execute_particle_particle_broad_search(particle_handler); + // Dummy Adaptive sparse contacts object and particle-particle broad search + AdaptiveSparseContacts dummy_adaptive_sparse_contacts; + container_manager.execute_particle_particle_broad_search( + particle_handler, dummy_adaptive_sparse_contacts); // Calling fine search container_manager.execute_particle_particle_fine_search( diff --git a/tests/dem/two_particles_multiple_contacts_parallel.cc b/tests/dem/two_particles_multiple_contacts_parallel.cc index b29b40addd..a657091dff 100644 --- a/tests/dem/two_particles_multiple_contacts_parallel.cc +++ b/tests/dem/two_particles_multiple_contacts_parallel.cc @@ -1,5 +1,3 @@ -#include - #include #include @@ -186,12 +184,13 @@ test() particle_handler.exchange_ghost_particles(); - container_manager.update_local_particles_in_cells(particle_handler, - false); + container_manager.update_local_particles_in_cells(particle_handler); - // Calling broad search + // Dummy Adaptive sparse contacts object and particle-particle broad + // search + AdaptiveSparseContacts dummy_adaptive_sparse_contacts; container_manager.execute_particle_particle_broad_search( - particle_handler); + particle_handler, dummy_adaptive_sparse_contacts); // Calling fine search container_manager.execute_particle_particle_fine_search(