diff --git a/CHANGELOG.md b/CHANGELOG.md index 1b5d7ee350..69739e30e7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,12 @@ All notable changes to the Lethe project will be documented in this file. The format is based on [Keep a Changelog](http://keepachangelog.com/). +## [Master] - 2023-12-11 + +### Fixed + +- MINOR Solid objects can now be restarted adequately in DEM. They will resume at the position they had at the end of the simulation [#959](https://github.com/lethe-cfd/lethe/pull/959) + ## [Master] - 2023-11-27 ### Removed diff --git a/applications_tests/lethe-particles/CMakeLists.txt b/applications_tests/lethe-particles/CMakeLists.txt index 7af2122804..1b3eb83be6 100644 --- a/applications_tests/lethe-particles/CMakeLists.txt +++ b/applications_tests/lethe-particles/CMakeLists.txt @@ -22,6 +22,18 @@ file(COPY periodic_gmsh_files/pipe.msh DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/ file(COPY load_balancing_solid_object_files/square.msh DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/load_balancing_solid_object.${_build_type}") file(COPY moving_float_files/square.msh DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/moving_float.${_build_type}") +file(COPY moving_receptacle_files/restart.particles DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/restart_moving_receptacle.${_build_type}") +file(COPY moving_receptacle_files/restart.pvdhandler DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/restart_moving_receptacle.${_build_type}") +file(COPY moving_receptacle_files/restart.simulationcontrol DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/restart_moving_receptacle.${_build_type}") +file(COPY moving_receptacle_files/restart.triangulation DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/restart_moving_receptacle.${_build_type}") +file(COPY moving_receptacle_files/restart.triangulation.info DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/restart_moving_receptacle.${_build_type}") +file(COPY moving_receptacle_files/restart.triangulation_fixed.data DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/restart_moving_receptacle.${_build_type}") +file(COPY moving_receptacle_files/restart.triangulation_variable.data DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/restart_moving_receptacle.${_build_type}") +file(COPY moving_receptacle_files/square.msh DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/restart_moving_receptacle.${_build_type}") +file(COPY moving_receptacle_files/restart.solid_object.00.displacement DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/restart_moving_receptacle.${_build_type}") +file(COPY moving_receptacle_files/restart.solid_object.00.dof DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/restart_moving_receptacle.${_build_type}") +file(COPY moving_receptacle_files/restart.solid_object.00.pvdhandler DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/restart_moving_receptacle.${_build_type}") + deal_ii_pickup_tests() if(CMAKE_BUILD_TYPE STREQUAL "Debug") diff --git a/applications_tests/lethe-particles/moving_receptacle_files/moving_receptacle.prm b/applications_tests/lethe-particles/moving_receptacle_files/moving_receptacle.prm new file mode 100644 index 0000000000..2bfdf37690 --- /dev/null +++ b/applications_tests/lethe-particles/moving_receptacle_files/moving_receptacle.prm @@ -0,0 +1,134 @@ +# Listing of Parameters +#---------------------- + +set dimension = 3 + +#--------------------------------------------------- +# Simulation Control +#--------------------------------------------------- + +subsection simulation control + set time step = 5e-5 + set time end = .7 + set log frequency = 500 + set output frequency = 0 + set output path = ./output/ + set output boundaries = false +end + +#--------------------------------------------------- +# Model parameters +#--------------------------------------------------- + +subsection model parameters + subsection contact detection + set contact detection method = dynamic + set dynamic contact search size coefficient = 0.9 + set neighborhood threshold = 1.3 + end + subsection load balancing + set load balance method = none + set frequency = 10000 + end + set particle particle contact force method = hertz_mindlin_limit_overlap + set rolling resistance torque method = constant_resistance + set particle wall contact force method = nonlinear + set integration method = velocity_verlet +end + +#--------------------------------------------------- +# Physical Properties +#--------------------------------------------------- + +subsection lagrangian physical properties + set gx = 9.81 + set gy = 0 + set gz = 0.0 + set number of particle types = 1 + subsection particle type 0 + set size distribution type = uniform + set diameter = 0.01 + set number of particles = 100 + set density particles = 2560 + set young modulus particles = 1e6 + set poisson ratio particles = 0.2 + set restitution coefficient particles = 0.5 + set friction coefficient particles = 0.5 + set rolling friction particles = 0.3 + end + set young modulus wall = 1e6 + set poisson ratio wall = 0.2 + set restitution coefficient wall = 0.9 + set friction coefficient wall = 0.5 + set rolling friction wall = 0.3 +end + +#--------------------------------------------------- +# Insertion Info +#--------------------------------------------------- + +subsection insertion info + set insertion method = volume + set inserted number of particles at each time step = 600 + set insertion frequency = 10000 + set insertion box minimum x = 0.0025 + set insertion box minimum y = 0.0025 + set insertion box minimum z = 0.0025 + set insertion box maximum x = 0.597 + set insertion box maximum y = 0.340 + set insertion box maximum z = 0.397 + set insertion distance threshold = 1.075 + set insertion random number range = 0.025 + set insertion random number seed = 19 +end + +#--------------------------------------------------- +# Insertion Info +#--------------------------------------------------- + +subsection restart + # Checkpointing parameters + set checkpoint = true + set frequency = 14000 + + # Output/input files + set filename = restart + + # Restarting parameter + set restart = false +end + + +#--------------------------------------------------- +# Mesh +#--------------------------------------------------- + +subsection mesh + set type = dealii + set grid type = subdivided_hyper_rectangle + set grid arguments = 2,1,1 : 0,0,0 : 0.8,0.4,0.4 : false + set initial refinement = 0 +end + +#--------------------------------------------------- +# Solid Objects +#--------------------------------------------------- + +subsection solid objects + set number of solids = 1 + subsection solid object 0 + subsection mesh + set type = gmsh + set file name = square.msh + set simplex = true + set initial refinement = 0 + end + + subsection translational velocity + set Function expression = if(t>0.4,if(t<0.6,0.1,0),0) ; 0 ; 0 + end + subsection angular velocity + set Function expression = 0 ; 0 ; 0 + end + end +end diff --git a/applications_tests/lethe-particles/moving_receptacle_files/restart.particles b/applications_tests/lethe-particles/moving_receptacle_files/restart.particles new file mode 100644 index 0000000000..261058b2e0 --- /dev/null +++ b/applications_tests/lethe-particles/moving_receptacle_files/restart.particles @@ -0,0 +1 @@ +0 0 100 100 100 diff --git a/applications_tests/lethe-particles/moving_receptacle_files/restart.pvdhandler b/applications_tests/lethe-particles/moving_receptacle_files/restart.pvdhandler new file mode 100644 index 0000000000..61fef595f1 --- /dev/null +++ b/applications_tests/lethe-particles/moving_receptacle_files/restart.pvdhandler @@ -0,0 +1,2 @@ +0 +Time File diff --git a/applications_tests/lethe-particles/moving_receptacle_files/restart.simulationcontrol b/applications_tests/lethe-particles/moving_receptacle_files/restart.simulationcontrol new file mode 100644 index 0000000000..c95c3d7d6a --- /dev/null +++ b/applications_tests/lethe-particles/moving_receptacle_files/restart.simulationcontrol @@ -0,0 +1,8 @@ +Simulation control +dt_0 5e-05 +dt_1 5e-05 +dt_2 5e-05 +dt_3 5e-05 +CFL 0 +Time 0.7 +Iter 14000 diff --git a/applications_tests/lethe-particles/moving_receptacle_files/restart.solid_object.00.displacement b/applications_tests/lethe-particles/moving_receptacle_files/restart.solid_object.00.displacement new file mode 100644 index 0000000000..f913907d39 --- /dev/null +++ b/applications_tests/lethe-particles/moving_receptacle_files/restart.solid_object.00.displacement @@ -0,0 +1 @@ +22 serialization::archive 17 0 0 0 0 15 2.00000000000007880e-02 0.00000000000000000e+00 0.00000000000000000e+00 2.00000000000007880e-02 0.00000000000000000e+00 0.00000000000000000e+00 2.00000000000007880e-02 0.00000000000000000e+00 0.00000000000000000e+00 2.00000000000007880e-02 0.00000000000000000e+00 0.00000000000000000e+00 2.00000000000007880e-02 0.00000000000000000e+00 0.00000000000000000e+00 diff --git a/applications_tests/lethe-particles/moving_receptacle_files/restart.solid_object.00.dof b/applications_tests/lethe-particles/moving_receptacle_files/restart.solid_object.00.dof new file mode 100644 index 0000000000..fb34e48dbd --- /dev/null +++ b/applications_tests/lethe-particles/moving_receptacle_files/restart.solid_object.00.dof @@ -0,0 +1 @@ +22 serialization::archive 17 0 0 0 0 0 0 0 0 0 0 0 0 0 0 15 15 0 0 0 0 1 0 0 0 0 15 0 1 15 0 1 0 15 0 0 1 0 1 0 0 15 0 1 15 0 0 0 1 0 0 0 3 15 0 12 13 14 0 1 2 3 4 5 9 10 11 6 7 8 0 0 0 0 1 0 3 6 0 0 3 6 9 12 15 9 0 0 0 0 0 0 0 0 0 0 5 0 0 0 0 0 0 4 34 FESystem<2,3>[FE_SimplexP<2>(1)^3] 23 Policy::Sequential<2,3> diff --git a/applications_tests/lethe-particles/moving_receptacle_files/restart.solid_object.00.pvdhandler b/applications_tests/lethe-particles/moving_receptacle_files/restart.solid_object.00.pvdhandler new file mode 100644 index 0000000000..61fef595f1 --- /dev/null +++ b/applications_tests/lethe-particles/moving_receptacle_files/restart.solid_object.00.pvdhandler @@ -0,0 +1,2 @@ +0 +Time File diff --git a/applications_tests/lethe-particles/moving_receptacle_files/restart.triangulation b/applications_tests/lethe-particles/moving_receptacle_files/restart.triangulation new file mode 100644 index 0000000000..2d579810c0 Binary files /dev/null and b/applications_tests/lethe-particles/moving_receptacle_files/restart.triangulation differ diff --git a/applications_tests/lethe-particles/moving_receptacle_files/restart.triangulation.info b/applications_tests/lethe-particles/moving_receptacle_files/restart.triangulation.info new file mode 100644 index 0000000000..ef7f73fbd5 --- /dev/null +++ b/applications_tests/lethe-particles/moving_receptacle_files/restart.triangulation.info @@ -0,0 +1,2 @@ +version nproc n_attached_fixed_size_objs n_attached_variable_size_objs n_coarse_cells +5 1 0 1 2 diff --git a/applications_tests/lethe-particles/moving_receptacle_files/restart.triangulation_fixed.data b/applications_tests/lethe-particles/moving_receptacle_files/restart.triangulation_fixed.data new file mode 100644 index 0000000000..d423448477 Binary files /dev/null and b/applications_tests/lethe-particles/moving_receptacle_files/restart.triangulation_fixed.data differ diff --git a/applications_tests/lethe-particles/moving_receptacle_files/restart.triangulation_variable.data b/applications_tests/lethe-particles/moving_receptacle_files/restart.triangulation_variable.data new file mode 100644 index 0000000000..d934303259 Binary files /dev/null and b/applications_tests/lethe-particles/moving_receptacle_files/restart.triangulation_variable.data differ diff --git a/applications_tests/lethe-particles/moving_receptacle_files/square.geo b/applications_tests/lethe-particles/moving_receptacle_files/square.geo new file mode 100644 index 0000000000..b6a4eda4fc --- /dev/null +++ b/applications_tests/lethe-particles/moving_receptacle_files/square.geo @@ -0,0 +1,24 @@ +// Define a variable + +lc = 1; +H=0.5; +W=0.4; +x0=0.6; + + +Point(0) = {x0, 0, 0, lc}; +Point(1) = {x0, 0, W, lc}; +Point(2) = {x0, H, W, lc}; +Point(3) = {x0, H, 0, lc}; + +Line(0)={0,1}; +Line(1)={1,2}; +Line(2)={2,3}; +Line(3)={3,0}; + +Line Loop(1) = {0,1,2,3}; + +Plane Surface(1) = {1} ; + +Physical Surface(0) = {1}; + diff --git a/applications_tests/lethe-particles/moving_receptacle_files/square.msh b/applications_tests/lethe-particles/moving_receptacle_files/square.msh new file mode 100644 index 0000000000..ce25a606d4 --- /dev/null +++ b/applications_tests/lethe-particles/moving_receptacle_files/square.msh @@ -0,0 +1,41 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$Entities +4 4 1 0 +0 0.6 0 0 0 +1 0.6 0 0.4 0 +2 0.6 0.5 0.4 0 +3 0.6 0.5 0 0 +0 0.6 0 0 0.6 0 0.4 0 2 0 -1 +1 0.6 0 0.4 0.6 0.5 0.4 0 2 1 -2 +2 0.6 0.5 0 0.6 0.5 0.4 0 2 2 -3 +3 0.6 0 0 0.6 0.5 0 0 2 3 0 +1 0.6 0 0 0.6 0.5 0.4 1 0 4 0 1 2 3 +$EndEntities +$Nodes +5 5 1 5 +0 0 0 1 +1 +0.6 0 0 +0 1 0 1 +2 +0.6 0 0.4 +0 2 0 1 +3 +0.6 0.5 0.4 +0 3 0 1 +4 +0.6 0.5 0 +2 1 0 1 +5 +0.6 0.25 0.2 +$EndNodes +$Elements +1 4 1 4 +2 1 2 4 +1 2 3 5 +2 4 1 5 +3 1 2 5 +4 3 4 5 +$EndElements diff --git a/applications_tests/lethe-particles/restart_moving_receptacle.output b/applications_tests/lethe-particles/restart_moving_receptacle.output new file mode 100644 index 0000000000..c7d8c2f54a --- /dev/null +++ b/applications_tests/lethe-particles/restart_moving_receptacle.output @@ -0,0 +1,130 @@ + +********************* +Running on 1 rank(s) +********************* +DEM time-step is 3.69227% of Rayleigh time step +Reading triangulation + +Finished reading triangulation +Warning: expansion of particle-wall contact list is disabled. +This feature is useful in geometries with concave boundaries. + +***************************************************************** +Transient iteration: 14500 Time: 0.725 Time step: 5e-05 +***************************************************************** +| Variable | Min | Max | Average | Total | +| Contact list generation | 0.0000e+00 | 8.0000e+00 | 2.7586e-01 | 8.0000e+00 | +| Velocity magnitude | 2.0385e-05 | 4.6928e-01 | 1.7556e-02 | 1.7556e+00 | +| Angular velocity magnitude | 1.4224e-03 | 1.3591e+02 | 3.0792e+00 | 3.0792e+02 | +| Translational kinetic energy | 2.7849e-13 | 1.4760e-04 | 2.8498e-06 | 2.8498e-04 | +| Rotational kinetic energy | 2.7121e-14 | 2.4758e-04 | 3.2443e-06 | 3.2443e-04 | + +***************************************************************** +Transient iteration: 15000 Time: 0.75 Time step: 5e-05 +***************************************************************** +| Variable | Min | Max | Average | Total | +| Contact list generation | 0.0000e+00 | 1.0000e+01 | 6.0000e-01 | 1.8000e+01 | +| Velocity magnitude | 3.7519e-06 | 6.7722e-01 | 1.8740e-02 | 1.8740e+00 | +| Angular velocity magnitude | 3.5347e-03 | 1.3591e+02 | 2.9847e+00 | 2.9847e+02 | +| Translational kinetic energy | 9.4344e-15 | 3.0737e-04 | 4.8233e-06 | 4.8233e-04 | +| Rotational kinetic energy | 1.6747e-13 | 2.4758e-04 | 3.2513e-06 | 3.2513e-04 | + +***************************************************************** +Transient iteration: 15500 Time: 0.775 Time step: 5e-05 +***************************************************************** +| Variable | Min | Max | Average | Total | +| Contact list generation | 0.0000e+00 | 1.1000e+01 | 9.3548e-01 | 2.9000e+01 | +| Velocity magnitude | 3.1242e-06 | 6.0538e-01 | 1.4793e-02 | 1.4793e+00 | +| Angular velocity magnitude | 5.2241e-03 | 4.3128e+01 | 1.6292e+00 | 1.6292e+02 | +| Translational kinetic energy | 6.5418e-15 | 2.4562e-04 | 3.8443e-06 | 3.8443e-04 | +| Rotational kinetic energy | 3.6581e-13 | 2.4932e-05 | 7.0570e-07 | 7.0570e-05 | + +***************************************************************** +Transient iteration: 16000 Time: 0.8 Time step: 5e-05 +***************************************************************** +| Variable | Min | Max | Average | Total | +| Contact list generation | 0.0000e+00 | 1.1000e+01 | 1.1875e+00 | 3.8000e+01 | +| Velocity magnitude | 5.1992e-06 | 3.2067e-01 | 8.5319e-03 | 8.5319e-01 | +| Angular velocity magnitude | 2.5888e-03 | 4.0267e+01 | 1.0967e+00 | 1.0967e+02 | +| Translational kinetic energy | 1.8117e-14 | 6.8918e-05 | 1.2100e-06 | 1.2100e-04 | +| Rotational kinetic energy | 8.9831e-14 | 2.1734e-05 | 4.0814e-07 | 4.0814e-05 | + +***************************************************************** +Transient iteration: 16500 Time: 0.825 Time step: 5e-05 +***************************************************************** +| Variable | Min | Max | Average | Total | +| Contact list generation | 0.0000e+00 | 1.1000e+01 | 1.3030e+00 | 4.3000e+01 | +| Velocity magnitude | 6.1075e-06 | 2.6943e-01 | 4.7210e-03 | 4.7210e-01 | +| Angular velocity magnitude | 1.3688e-03 | 4.0267e+01 | 6.1168e-01 | 6.1168e+01 | +| Translational kinetic energy | 2.5000e-14 | 4.8653e-05 | 5.6392e-07 | 5.6392e-05 | +| Rotational kinetic energy | 2.5114e-14 | 2.1734e-05 | 2.3115e-07 | 2.3115e-05 | + +***************************************************************** +Transient iteration: 17000 Time: 0.85 Time step: 5e-05 +***************************************************************** +| Variable | Min | Max | Average | Total | +| Contact list generation | 0.0000e+00 | 1.1000e+01 | 1.4412e+00 | 4.9000e+01 | +| Velocity magnitude | 2.7954e-06 | 4.5913e-01 | 7.0493e-03 | 7.0493e-01 | +| Angular velocity magnitude | 3.1862e-03 | 4.0267e+01 | 4.9273e-01 | 4.9273e+01 | +| Translational kinetic energy | 5.2373e-15 | 1.4128e-04 | 1.6556e-06 | 1.6556e-04 | +| Rotational kinetic energy | 1.3608e-13 | 2.1734e-05 | 2.1890e-07 | 2.1890e-05 | + +***************************************************************** +Transient iteration: 17500 Time: 0.875 Time step: 5e-05 +***************************************************************** +| Variable | Min | Max | Average | Total | +| Contact list generation | 0.0000e+00 | 1.1000e+01 | 1.6000e+00 | 5.6000e+01 | +| Velocity magnitude | 2.0135e-06 | 4.2505e-01 | 5.9059e-03 | 5.9059e-01 | +| Angular velocity magnitude | 1.1980e-03 | 2.8965e+01 | 3.6601e-01 | 3.6601e+01 | +| Translational kinetic energy | 2.7170e-15 | 1.2109e-04 | 1.3096e-06 | 1.3096e-04 | +| Rotational kinetic energy | 1.9237e-14 | 1.1246e-05 | 1.1363e-07 | 1.1363e-05 | + +***************************************************************** +Transient iteration: 18000 Time: 0.9 Time step: 5e-05 +***************************************************************** +| Variable | Min | Max | Average | Total | +| Contact list generation | 0.0000e+00 | 1.1000e+01 | 1.6667e+00 | 6.0000e+01 | +| Velocity magnitude | 1.7462e-06 | 1.3054e-01 | 2.4122e-03 | 2.4122e-01 | +| Angular velocity magnitude | 2.9066e-03 | 2.8965e+01 | 3.3427e-01 | 3.3427e+01 | +| Translational kinetic energy | 2.0435e-15 | 1.1421e-05 | 1.8021e-07 | 1.8021e-05 | +| Rotational kinetic energy | 1.1324e-13 | 1.1246e-05 | 1.1249e-07 | 1.1249e-05 | + +***************************************************************** +Transient iteration: 18500 Time: 0.925 Time step: 5e-05 +***************************************************************** +| Variable | Min | Max | Average | Total | +| Contact list generation | 0.0000e+00 | 1.1000e+01 | 1.7027e+00 | 6.3000e+01 | +| Velocity magnitude | 9.9271e-08 | 1.6718e-01 | 3.0282e-03 | 3.0282e-01 | +| Angular velocity magnitude | 3.3061e-03 | 4.0013e+00 | 8.4518e-02 | 8.4518e+00 | +| Translational kinetic energy | 6.6047e-18 | 1.8731e-05 | 2.9386e-07 | 2.9386e-05 | +| Rotational kinetic energy | 1.4652e-13 | 2.1460e-07 | 2.1813e-09 | 2.1813e-07 | + +***************************************************************** +Transient iteration: 19000 Time: 0.95 Time step: 5e-05 +***************************************************************** +| Variable | Min | Max | Average | Total | +| Contact list generation | 0.0000e+00 | 1.1000e+01 | 1.7368e+00 | 6.6000e+01 | +| Velocity magnitude | 4.0188e-06 | 1.3331e-01 | 2.3277e-03 | 2.3277e-01 | +| Angular velocity magnitude | 2.1452e-03 | 4.0013e+00 | 8.4088e-02 | 8.4088e+00 | +| Translational kinetic energy | 1.0824e-14 | 1.1911e-05 | 1.7588e-07 | 1.7588e-05 | +| Rotational kinetic energy | 6.1685e-14 | 2.1460e-07 | 2.1818e-09 | 2.1818e-07 | + +***************************************************************** +Transient iteration: 19500 Time: 0.975 Time step: 5e-05 +***************************************************************** +| Variable | Min | Max | Average | Total | +| Contact list generation | 0.0000e+00 | 1.1000e+01 | 1.7692e+00 | 6.9000e+01 | +| Velocity magnitude | 1.2660e-06 | 1.6892e-01 | 2.3707e-03 | 2.3707e-01 | +| Angular velocity magnitude | 1.1465e-03 | 1.5682e-01 | 4.3282e-02 | 4.3282e+00 | +| Translational kinetic energy | 1.0742e-15 | 1.9124e-05 | 2.1616e-07 | 2.1616e-05 | +| Rotational kinetic energy | 1.7619e-14 | 3.2966e-10 | 3.4652e-11 | 3.4652e-09 | + +***************************************************************** +Transient iteration: 20000 Time: 1 Time step: 5e-05 +***************************************************************** +| Variable | Min | Max | Average | Total | +| Contact list generation | 0.0000e+00 | 1.1000e+01 | 1.7500e+00 | 7.0000e+01 | +| Velocity magnitude | 2.1238e-06 | 9.7668e-02 | 1.1679e-03 | 1.1679e-01 | +| Angular velocity magnitude | 3.6932e-04 | 7.2734e-01 | 5.3802e-02 | 5.3802e+00 | +| Translational kinetic energy | 3.0231e-15 | 6.3932e-06 | 6.5134e-08 | 6.5134e-06 | +| Rotational kinetic energy | 1.8283e-15 | 7.0911e-09 | 1.1543e-10 | 1.1543e-08 | diff --git a/applications_tests/lethe-particles/restart_moving_receptacle.prm b/applications_tests/lethe-particles/restart_moving_receptacle.prm new file mode 100644 index 0000000000..9f890ab95d --- /dev/null +++ b/applications_tests/lethe-particles/restart_moving_receptacle.prm @@ -0,0 +1,132 @@ +# Listing of Parameters +#---------------------- + +set dimension = 3 + +#--------------------------------------------------- +# Simulation Control +#--------------------------------------------------- + +subsection simulation control + set time step = 5e-5 + set time end = 1 + set log frequency = 500 + set output frequency = 0 +end + +#--------------------------------------------------- +# Model parameters +#--------------------------------------------------- + +subsection model parameters + subsection contact detection + set contact detection method = dynamic + set dynamic contact search size coefficient = 0.9 + set neighborhood threshold = 1.3 + end + subsection load balancing + set load balance method = none + set frequency = 10000 + end + set particle particle contact force method = hertz_mindlin_limit_overlap + set rolling resistance torque method = constant_resistance + set particle wall contact force method = nonlinear + set integration method = velocity_verlet +end + +#--------------------------------------------------- +# Physical Properties +#--------------------------------------------------- + +subsection lagrangian physical properties + set gx = 9.81 + set gy = 0 + set gz = 0.0 + set number of particle types = 1 + subsection particle type 0 + set size distribution type = uniform + set diameter = 0.01 + set number of particles = 100 + set density particles = 2560 + set young modulus particles = 1e6 + set poisson ratio particles = 0.2 + set restitution coefficient particles = 0.5 + set friction coefficient particles = 0.5 + set rolling friction particles = 0.3 + end + set young modulus wall = 1e6 + set poisson ratio wall = 0.2 + set restitution coefficient wall = 0.9 + set friction coefficient wall = 0.5 + set rolling friction wall = 0.3 +end + +#--------------------------------------------------- +# Insertion Info +#--------------------------------------------------- + +subsection insertion info + set insertion method = volume + set inserted number of particles at each time step = 0 + set insertion frequency = 10000 + set insertion box minimum x = 0.0025 + set insertion box minimum y = 0.0025 + set insertion box minimum z = 0.0025 + set insertion box maximum x = 0.597 + set insertion box maximum y = 0.340 + set insertion box maximum z = 0.397 + set insertion distance threshold = 1.075 + set insertion random number range = 0.025 + set insertion random number seed = 19 +end + +#--------------------------------------------------- +# Insertion Info +#--------------------------------------------------- + +subsection restart + # Checkpointing parameters + set checkpoint = false + set frequency = 10000 + + # Output/input files + set filename = restart + + # Restarting parameter + set restart = true +end + + +#--------------------------------------------------- +# Mesh +#--------------------------------------------------- + +subsection mesh + set type = dealii + set grid type = subdivided_hyper_rectangle + set grid arguments = 2,1,1 : 0,0,0 : 0.8,0.4,0.4 : false + set initial refinement = 0 +end + +#--------------------------------------------------- +# Solid Objects +#--------------------------------------------------- + +subsection solid objects + set number of solids = 1 + subsection solid object 0 + subsection mesh + set type = gmsh + set file name = square.msh + set simplex = true + set initial refinement = 0 + end + + subsection translational velocity + set Function expression = if(t>0.4,if(t<0.6,0.1,0),0) ; 0 ; 0 + end + subsection angular velocity + set Function expression = 0 ; 0 ; 0 + end + end +end diff --git a/include/core/serial_solid.h b/include/core/serial_solid.h index 6118b64e27..55eb2f2669 100644 --- a/include/core/serial_solid.h +++ b/include/core/serial_solid.h @@ -147,7 +147,7 @@ class SerialSolid * * @return The angular velocity of the solid object */ - Tensor<1, 3> + inline Tensor<1, 3> get_angular_velocity() const { return this->current_angular_velocity; @@ -158,7 +158,7 @@ class SerialSolid * * @return The center of rotation of the solid object */ - Point<3> + inline Point<3> get_center_of_rotation() const { if constexpr (spacedim == 3) @@ -175,7 +175,7 @@ class SerialSolid * * @return id of the solid */ - unsigned int + inline unsigned int get_solid_id() const { return id; @@ -204,11 +204,11 @@ class SerialSolid write_output_results(std::shared_ptr simulation_control); /** - * @brief read solid base triangulation checkpoint + * @brief read solid base triangulation checkpoint and replaces the + * triangulation at the location where it was when it was checkpointed * * @param prefix_name The prefix of the checkpoint of the simulation * - * @warning This is currently not implemented */ void read_checkpoint(std::string prefix_name); @@ -217,8 +217,6 @@ class SerialSolid * @brief write solid base triangulation checkpoint * * @param prefix_name The prefix of the checkpoint of the simulation - * - * @warning This is currently not implemented */ void write_checkpoint(std::string prefix_name); @@ -242,6 +240,15 @@ class SerialSolid displacement_since_mapped = 0; } + /** + * @brief Moves the vertices of the solid triangulation using the displacement + * vector. This is mostly used when restarting simulations which contain + * a serial solid in order to replace it at the position it was when the + * simulation finished. + */ + void + move_solid_triangulation_with_displacement(); + /** * @brief Rotates the grid by a given angle around an axis * diff --git a/include/dem/read_checkpoint.h b/include/dem/read_checkpoint.h index 46f0d20712..e9b122dec8 100644 --- a/include/dem/read_checkpoint.h +++ b/include/dem/read_checkpoint.h @@ -12,12 +12,10 @@ * the top level of the Lethe distribution. * * --------------------------------------------------------------------- - - * - * Author: Shahab Golshan, Polytechnique Montreal, 2019 */ #include +#include #include @@ -49,15 +47,18 @@ using namespace std; * @param grid_pvdhandler PVD handler for post-processing * @param triangulation Triangulation * @param particle_handler Particle handler + * @param solid_objects Vector of solids objects used in DEM simulations */ template void -read_checkpoint(TimerOutput &computing_timer, - const DEMSolverParameters &dem_parameters, - std::shared_ptr &simulation_control, - PVDHandler &particles_pvdhandler, - PVDHandler &grid_pvdhandler, - parallel::distributed::Triangulation &triangulation, - Particles::ParticleHandler &particle_handler); +read_checkpoint( + TimerOutput &computing_timer, + const DEMSolverParameters &dem_parameters, + std::shared_ptr &simulation_control, + PVDHandler &particles_pvdhandler, + PVDHandler &grid_pvdhandler, + parallel::distributed::Triangulation &triangulation, + Particles::ParticleHandler &particle_handler, + std::vector>> &solid_objects); #endif /* read_checkpoint_h */ diff --git a/include/dem/write_checkpoint.h b/include/dem/write_checkpoint.h index 7bc8fec0a7..7d7d35a0ee 100644 --- a/include/dem/write_checkpoint.h +++ b/include/dem/write_checkpoint.h @@ -18,6 +18,7 @@ */ #include +#include #include @@ -49,20 +50,23 @@ using namespace std; * @param grid_pvdhandler PVD handler for post-processing * @param triangulation Triangulation * @param particle_handler Particle handler + * @param solid_objects Vector of solids objects used in DEM simulations * @param pcout Printing in parallel * @param mpi_communicator */ template void -write_checkpoint(TimerOutput &computing_timer, - const DEMSolverParameters &dem_parameters, - std::shared_ptr &simulation_control, - PVDHandler &particles_pvdhandler, - PVDHandler &grid_pvdhandler, - parallel::distributed::Triangulation &triangulation, - Particles::ParticleHandler &particle_handler, - const ConditionalOStream &pcout, - MPI_Comm &mpi_communicator); +write_checkpoint( + TimerOutput &computing_timer, + const DEMSolverParameters &dem_parameters, + std::shared_ptr &simulation_control, + PVDHandler &particles_pvdhandler, + PVDHandler &grid_pvdhandler, + parallel::distributed::Triangulation &triangulation, + Particles::ParticleHandler &particle_handler, + std::vector>> &solid_objects, + const ConditionalOStream &pcout, + MPI_Comm &mpi_communicator); #endif /* write_checkpoint_h */ diff --git a/source/core/serial_solid.cc b/source/core/serial_solid.cc index d7f99435cc..f0c017ace5 100644 --- a/source/core/serial_solid.cc +++ b/source/core/serial_solid.cc @@ -429,6 +429,40 @@ SerialSolid::displace_solid_triangulation() } } +template +void +SerialSolid::move_solid_triangulation_with_displacement() +{ + const unsigned int dofs_per_cell = displacement_fe->dofs_per_cell; + + // Set of vertices which were displaced already. A set is used + // for efficiency when searching over it. + std::set dof_vertex_displaced; + std::vector local_dof_indices(dofs_per_cell); + + // Loop over all cells + for (const auto &cell : displacement_dh.active_cell_iterators()) + { + for (unsigned int vertex = 0; + vertex < cell->reference_cell().n_vertices(); + ++vertex) + { + if (!dof_vertex_displaced.count(cell->vertex_index(vertex))) + { + const auto dof_index = cell->vertex_dof_index(vertex, 0); + Point &vertex_position = cell->vertex(vertex); + + for (unsigned d = 0; d < spacedim; ++d) + { + vertex_position[d] += displacement[dof_index + d]; + } + + dof_vertex_displaced.insert(cell->vertex_index(vertex)); + } + } + } +} + template void @@ -467,7 +501,8 @@ SerialSolid::write_output_results( const std::string folder = simulation_control->get_output_path(); const std::string solution_name = simulation_control->get_output_name() + - "_solid_" + Utilities::int_to_string(id, 2); + ".solid_object." + + Utilities::int_to_string(id, 2); const unsigned int iter = simulation_control->get_step_number(); const double time = simulation_control->get_current_time(); const unsigned int group_files = simulation_control->get_group_files(); @@ -484,77 +519,57 @@ SerialSolid::write_output_results( template void -SerialSolid::write_checkpoint(std::string /*prefix*/) +SerialSolid::write_checkpoint(std::string prefix) { - // SolutionTransfer, spacedim> system_trans_vectors( - // this->displacement_dh); - - // std::vector *> sol_set_transfer; - // sol_set_transfer.push_back(&displacement); - - // system_trans_vectors.prepare_for_serialization(sol_set_transfer); + // Checkpoint the DOF Handler + { + std::string file_name = + prefix + ".solid_object." + Utilities::int_to_string(id, 2) + ".dof"; + std::ofstream ofs(file_name); + boost::archive::text_oarchive oa(ofs); + displacement_dh.save(oa, 0); + } + // Checkpoint the displacement vector which we will use to shift the + // triangulation + { + std::string file_name = prefix + ".solid_object." + + Utilities::int_to_string(id, 2) + ".displacement"; + std::ofstream ofs(file_name); + boost::archive::text_oarchive oa(ofs); + displacement.save(oa, 0); + } - // if (auto tria = dynamic_cast *>( - // this->solid_tria.get())) - // { - // std::string triangulationName = prefix + ".triangulation"; - // tria->save(prefix + ".triangulation"); - // } + // Re-read pvd handler from output files + pvdhandler.save(prefix + "_solid_object_" + Utilities::int_to_string(id, 2)); } template void -SerialSolid::read_checkpoint(std::string /*prefix*/) +SerialSolid::read_checkpoint(std::string prefix) { - throw(std::runtime_error( - "Restarting DEM simulations with floating meshes is currently not possible.")); - // Setup an un-refined triangulation before loading - // setup_triangulation(true); - - // Read the triangulation from the checkpoint - // const std::string filename = prefix + ".triangulation"; - // std::ifstream in(filename.c_str()); - // if (!in) - // AssertThrow(false, - // ExcMessage( - // std::string("You are trying to read a solid triangulation, " - // "but the restart file <") + - // filename + "> does not appear to exist!")); - - // try - // { - // if (auto tria = dynamic_cast - // *>( - // this->solid_tria.get())) - // tria->load(filename.c_str()); - // } - // catch (...) - // { - // AssertThrow(false, - // ExcMessage("Cannot open snapshot mesh file or read the " - // "triangulation stored there.")); - // } - - // Setup dof-handler for solid and displacement - // solid_dh.distribute_dofs(*fe); - // setup_displacement(); - - - //// Read displacement vector - // std::vector *> x_system(1); - // x_system[0] = &(displacement); - - // SolutionTransfer, spacedim> system_trans_vectors( - // this->displacement_dh); - - // system_trans_vectors.deserialize(x_system); - // displacement; - - //// Reset triangulation position using displacement vector - // move_solid_triangulation_with_displacement(); - - - // We did not checkpoint particles, we re-create them from scratch + // Import dof handler + { + // displacement_dh.distribute_dofs(fe); + std::string file_name = + prefix + ".solid_object." + Utilities::int_to_string(id, 2) + ".dof"; + std::ifstream ifs(file_name); + boost::archive::text_iarchive ia(ifs); + displacement_dh.load(ia, 0); + } + // Import nodal counts + { + std::string file_name = prefix + ".solid_object." + + Utilities::int_to_string(id, 2) + ".displacement"; + std::ifstream ifs(file_name); + boost::archive::text_iarchive ia(ifs); + displacement.load(ia, 0); + } + + // Re-read pvd handler from output files + pvdhandler.read(prefix + ".solid_object." + Utilities::int_to_string(id, 2)); + + // Reset triangulation position using displacement vector + move_solid_triangulation_with_displacement(); } diff --git a/source/dem/dem.cc b/source/dem/dem.cc index e26fccdb4a..fb09532782 100644 --- a/source/dem/dem.cc +++ b/source/dem/dem.cc @@ -1115,13 +1115,6 @@ DEMSolver::solve() triangulation, parameters.boundary_conditions); - // Store information about floating mesh/background mesh intersection - for (unsigned int i_solid = 0; i_solid < solids.size(); ++i_solid) - { - floating_mesh_info[i_solid] = - solids[i_solid]->map_solid_in_background_triangulation(triangulation); - } - if (parameters.restart.restart == true) { read_checkpoint(computing_timer, @@ -1130,7 +1123,8 @@ DEMSolver::solve() particles_pvdhandler, grid_pvdhandler, triangulation, - particle_handler); + particle_handler, + solids); displacement.resize(particle_handler.get_max_local_particle_index()); force.resize(displacement.size()); @@ -1141,6 +1135,13 @@ DEMSolver::solve() checkpoint_step = true; } + // Store information about floating mesh/background mesh intersection + for (unsigned int i_solid = 0; i_solid < solids.size(); ++i_solid) + { + floating_mesh_info[i_solid] = + solids[i_solid]->map_solid_in_background_triangulation(triangulation); + } + // 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 @@ -1482,6 +1483,7 @@ DEMSolver::solve() grid_pvdhandler, triangulation, particle_handler, + solids, pcout, mpi_communicator); } diff --git a/source/dem/read_checkpoint.cc b/source/dem/read_checkpoint.cc index 50b5e56d11..75d3591379 100644 --- a/source/dem/read_checkpoint.cc +++ b/source/dem/read_checkpoint.cc @@ -4,13 +4,15 @@ using namespace dealii; template void -read_checkpoint(TimerOutput &computing_timer, - const DEMSolverParameters ¶meters, - std::shared_ptr &simulation_control, - PVDHandler &particles_pvdhandler, - PVDHandler &grid_pvdhandler, - parallel::distributed::Triangulation &triangulation, - Particles::ParticleHandler &particle_handler) +read_checkpoint( + TimerOutput &computing_timer, + const DEMSolverParameters ¶meters, + std::shared_ptr &simulation_control, + PVDHandler &particles_pvdhandler, + PVDHandler &grid_pvdhandler, + parallel::distributed::Triangulation &triangulation, + Particles::ParticleHandler &particle_handler, + std::vector>> &solid_objects) { TimerOutput::Scope timer(computing_timer, "read_checkpoint"); std::string prefix = parameters.restart.filename; @@ -58,6 +60,12 @@ read_checkpoint(TimerOutput &computing_timer, // Unpack the information in the particle handler particle_handler.deserialize(); + + // Load solid objects + for (unsigned int i = 0; i < solid_objects.size(); ++i) + { + solid_objects[i]->read_checkpoint(prefix); + } } template void @@ -67,7 +75,8 @@ read_checkpoint(TimerOutput &computing_timer, PVDHandler &particles_pvdhandler, PVDHandler &grid_pvdhandler, parallel::distributed::Triangulation<2> &triangulation, - Particles::ParticleHandler<2> &particle_handler); + Particles::ParticleHandler<2> &particle_handler, + std::vector>> &solid_objects); template void read_checkpoint(TimerOutput &computing_timer, @@ -76,4 +85,5 @@ read_checkpoint(TimerOutput &computing_timer, PVDHandler &particles_pvdhandler, PVDHandler &grid_pvdhandler, parallel::distributed::Triangulation<3> &triangulation, - Particles::ParticleHandler<3> &particle_handler); + Particles::ParticleHandler<3> &particle_handler, + std::vector>> &solid_objects); diff --git a/source/dem/write_checkpoint.cc b/source/dem/write_checkpoint.cc index 32ccb86b04..6efc80da06 100644 --- a/source/dem/write_checkpoint.cc +++ b/source/dem/write_checkpoint.cc @@ -4,15 +4,17 @@ using namespace dealii; template void -write_checkpoint(TimerOutput &computing_timer, - const DEMSolverParameters ¶meters, - std::shared_ptr &simulation_control, - PVDHandler &particles_pvdhandler, - PVDHandler &grid_pvdhandler, - parallel::distributed::Triangulation &triangulation, - Particles::ParticleHandler &particle_handler, - const ConditionalOStream &pcout, - MPI_Comm &mpi_communicator) +write_checkpoint( + TimerOutput &computing_timer, + const DEMSolverParameters ¶meters, + std::shared_ptr &simulation_control, + PVDHandler &particles_pvdhandler, + PVDHandler &grid_pvdhandler, + parallel::distributed::Triangulation &triangulation, + Particles::ParticleHandler &particle_handler, + std::vector>> &solid_objects, + const ConditionalOStream &pcout, + MPI_Comm &mpi_communicator) { TimerOutput::Scope timer(computing_timer, "write_checkpoint"); @@ -43,6 +45,12 @@ write_checkpoint(TimerOutput &computing_timer, std::string particle_filename = prefix + ".particles"; std::ofstream output(particle_filename.c_str()); output << oss.str() << std::endl; + + // Checkpoint the serial solid objects one by one + for (unsigned int i = 0; i < solid_objects.size(); ++i) + { + solid_objects[i]->write_checkpoint(prefix); + } } template void @@ -53,8 +61,9 @@ write_checkpoint(TimerOutput &computing_timer, PVDHandler &grid_pvdhandler, parallel::distributed::Triangulation<2> &triangulation, Particles::ParticleHandler<2> &particle_handler, - const ConditionalOStream &pcout, - MPI_Comm &mpi_communicator); + std::vector>> &solid_objects, + const ConditionalOStream &pcout, + MPI_Comm &mpi_communicator); template void write_checkpoint(TimerOutput &computing_timer, @@ -64,5 +73,6 @@ write_checkpoint(TimerOutput &computing_timer, PVDHandler &grid_pvdhandler, parallel::distributed::Triangulation<3> &triangulation, Particles::ParticleHandler<3> &particle_handler, - const ConditionalOStream &pcout, - MPI_Comm &mpi_communicator); + std::vector>> &solid_objects, + const ConditionalOStream &pcout, + MPI_Comm &mpi_communicator);