Backup

CMakeLists.txt

@@ -60,7 +60,7 @@ add_executable(alt_sintered_test test/alt_sintered.cpp writer.cpp test/mass_dist
 # Particle is collided with a plane
 add_executable(surface_test test/surface.cpp writer.cpp test/compute_energy.cpp)
 # Same as hamaker_test, but with truncation and neighbor tracking
-add_executable(hamaker_neighbor_test test/hamaker_neighbor.cpp)
+add_executable(hamaker_neighbor_test test/hamaker_neighbor.cpp test/compute_energy.cpp test/mass_distribution.cpp)
 
 add_executable(sintered_2_test test/sintered_2.cpp writer.cpp)
 

include/libgran/granular_system/granular_system_neighbor_list.h

@@ -12,10 +12,10 @@
 #ifndef LIBGRAN_USE_OMP
 #define binary_system_implementation rotational_binary_system
 #include <libtimestep/rotational_system/rotational_binary_system.h>
-#pragma message("libgran: using C++ 17 parallel algorithms")
+#error "Neighbor lists are not supported with C++ 17 algorithms yet"
 #else
-#define binary_system_implementation rotational_binary_system_omp
-#include <libtimestep/rotational_system/rotational_binary_system_omp.h>
+#define binary_system_implementation rotational_binary_system_neighbors_omp
+#include <libtimestep/rotational_system/rotational_binary_system_neighbors_omp.h>
 #pragma message("libgran: using OpenMP parallelization")
 #endif //LIBGRAN_USE_OMP
 
@@ -42,23 +42,23 @@ template <
         typename step_handler_t,
         typename binary_force_functor_container_t,
         typename unary_force_functor_container_t>
-class granular_system_neignbor_list : public binary_system_implementation<field_value_t, real_t, integrator_t, step_handler_t,
-        granular_system_neignbor_list<field_value_t, real_t, integrator_t, step_handler_t,
+class granular_system_neighbor_list : public binary_system_implementation<field_value_t, real_t, integrator_t, step_handler_t,
+        granular_system_neighbor_list<field_value_t, real_t, integrator_t, step_handler_t,
                 binary_force_functor_container_t, unary_force_functor_container_t>,
         (std::tuple_size<decltype(unary_force_functor_container_t::unary_force_functors)>::value > 0)> {
 public:
     typedef std::vector<field_value_t> field_container_t;
 
-    granular_system_neignbor_list(granular_system_neignbor_list const &) = delete;
+    granular_system_neighbor_list(granular_system_neighbor_list const &) = delete;
 
-    granular_system_neignbor_list(field_container_t x0, field_container_t v0,
+    granular_system_neighbor_list(real_t r_verlet, field_container_t x0, field_container_t v0,
                     field_container_t theta0, field_container_t omega0,
                     real_t t0, field_value_t field_zero, real_t real_zero, step_handler_t<field_container_t, field_value_t> & step_handler,
                     binary_force_functor_container_t binary_force_functors, unary_force_functor_container_t unary_force_functors) :
             binary_system_implementation<field_value_t, real_t, integrator_t, step_handler_t,
-                    granular_system_neignbor_list<field_value_t, real_t, integrator_t, step_handler_t, binary_force_functor_container_t, unary_force_functor_container_t>,
+                    granular_system_neighbor_list<field_value_t, real_t, integrator_t, step_handler_t, binary_force_functor_container_t, unary_force_functor_container_t>,
                     (std::tuple_size<decltype(unary_force_functor_container_t::unary_force_functors)>::value > 0)>
-                    (std::move(x0), std::move(v0), std::move(theta0),
+                    (x0.size(), r_verlet, std::move(x0), std::move(v0), std::move(theta0),
                      std::move(omega0), t0, field_zero, real_zero, *this, step_handler),
             binary_force_functors(std::move(binary_force_functors)),
             unary_force_functors(std::move(unary_force_functors)) {}
@@ -82,8 +82,6 @@ class granular_system_neignbor_list : public binary_system_implementation<field_
     }
 
 private:
-    const real_t r_verlet;  // Verlet radius for neighbor list computation
-
     binary_force_functor_container_t binary_force_functors;
     unary_force_functor_container_t unary_force_functors;
 };

test/hamaker_neighbor.cpp

@@ -22,3 +22,128 @@
 #include <libgran/contact_force/contact_force.h>
 #include <libgran/hamaker_force/hamaker_force.h>
 #include <libgran/granular_system/granular_system_neighbor_list.h>
+
+#include "compute_energy.h"
+#include "mass_distribution.h"
+
+int main() {
+    enable_fp_exceptions();
+
+    // General simulation parameters
+    const double dt = 1e-13;
+    const double t_tot = 3.0e-7;
+    const auto n_steps = size_t(t_tot / dt);
+    const size_t n_dumps = 300;
+    const size_t dump_period = n_steps / n_dumps;
+
+    // General parameters
+    const double rho = 1700.0;
+    const double r_part = 1.4e-8;
+    const double mass = 4.0 / 3.0 * M_PI * pow(r_part, 3.0) * rho;
+    const double inertia = 2.0 / 5.0 * mass * pow(r_part, 2.0);
+
+    // Parameters for the contact model
+    const double k = 10000.0;
+    const double gamma_n = 5.0e-9;
+    const double mu = 1.0;
+    const double phi = 1.0;
+    const double mu_o = 0.1;
+    const double gamma_t = 0.2 * gamma_n;
+    const double gamma_r = 0.05 * gamma_n;
+    const double gamma_o = 0.05 * gamma_n;
+
+    // Parameters for the Van der Waals model
+    const double A = 1.0e-20;
+    const double h0 = 1.0e-9;
+
+    // Initialize the particles
+    std::vector<Eigen::Vector3d> x0, v0, theta0, omega0;
+    // 18 particles in the x-y plane
+    for (size_t i = 0; i < 3; i ++) {
+        auto y = -3.0 * r_part + double(i) * 2.0 * r_part;
+        for (size_t j = 0; j < 3; j ++) {
+            auto x = -3.0 * r_part + double(j) * 2.0 * r_part;
+            x0.emplace_back(x, y, 0.0);
+            x0.emplace_back(x, y, -2.0*r_part);
+        }
+    }
+
+    // 1 particle above the plane
+    x0.emplace_back(0, 0, 3.0*r_part);
+
+    // Initialize the remaining buffers
+    v0.resize(x0.size());
+    theta0.resize(x0.size());
+    omega0.resize(x0.size());
+    std::fill(v0.begin(), v0.end(), Eigen::Vector3d::Zero());
+    std::fill(theta0.begin(), theta0.end(), Eigen::Vector3d::Zero());
+    std::fill(omega0.begin(), omega0.end(), Eigen::Vector3d::Zero());
+
+    // Set the initial velocity of the above-plane particle
+    v0[v0.size() - 1] = {1.0, 0.0, -10.0};
+
+    // Create an instance of step_handler
+    // Using field type Eigen::Vector3d with container std::vector
+    rotational_step_handler<std::vector<Eigen::Vector3d>, Eigen::Vector3d> step_handler_instance;
+
+    // Create an instance of contact force model
+    // Using field type Eigen::Vector3d with real type double
+    contact_force_functor<Eigen::Vector3d, double> contact_force_model(x0.size(),
+                                                                       k, gamma_n, k, gamma_t, mu, phi, k, gamma_r, mu_o, phi, k, gamma_o, mu_o, phi,
+                                                                       r_part, mass, inertia, dt, Eigen::Vector3d::Zero(), 0.0);
+
+    // Create an instance of Hamaker model
+    hamaker_functor<Eigen::Vector3d, double> hamaker_model(A, h0,
+                                                           r_part, mass, Eigen::Vector3d::Zero(), 0.0);
+
+    binary_force_functor_container<Eigen::Vector3d, double,
+            contact_force_functor<Eigen::Vector3d, double>,
+            hamaker_functor<Eigen::Vector3d, double>>
+            binary_force_functors{contact_force_model, hamaker_model};
+
+    unary_force_functor_container<Eigen::Vector3d, double>
+            unary_force_functors;
+
+    // Create an instance of granular_system using the contact force model
+    // Using velocity Verlet integrator for rotational systems and a default
+    // step handler for rotational systems
+    granular_system_neighbor_list<Eigen::Vector3d, double, rotational_velocity_verlet_half,
+            rotational_step_handler,
+            binary_force_functor_container<Eigen::Vector3d, double,
+                    contact_force_functor<Eigen::Vector3d, double>,
+                    hamaker_functor<Eigen::Vector3d, double>>,
+            unary_force_functor_container<Eigen::Vector3d, double>> system(5.0 * r_part, x0,
+                                                                           v0, theta0, omega0, 0.0, Eigen::Vector3d::Zero(), 0.0,
+                                                                           step_handler_instance, binary_force_functors, unary_force_functors);
+
+    auto start_time = std::chrono::high_resolution_clock::now();
+    for (size_t n = 0; n < n_steps; n ++) {
+        if (n % 20 == 0)
+            system.update_neighbor_list();
+        system.do_step(dt);
+    }
+    auto end_time = std::chrono::high_resolution_clock::now();
+    auto duration = std::chrono::duration_cast<std::chrono::seconds>(end_time - start_time).count();
+    std::cout << "Elapsed time: " << duration << " sec" << std::endl;
+
+    double target_linear_momentum = 1.96373e-19;
+    double target_ke = 8.33959e-20;
+    double target_rg = 3.56253e-08;
+
+    double linear_momentum_tolerance = 1.0; // Percent
+    double ke_tolerance = 1.0; // Percent
+    double rg_tolerance = 1.0; // Percent
+
+    double linear_momentum = compute_linear_momentum(system.get_v(), mass);
+    double ke = compute_total_kinetic_energy(system.get_v(), system.get_omega(), mass, inertia);
+    double rg = radius_of_gyration(system.get_x());
+
+    if (abs(ke - target_ke) / target_ke > ke_tolerance / 100.0 ||
+        abs(linear_momentum - target_linear_momentum) / target_linear_momentum > linear_momentum_tolerance / 100.0 ||
+        abs(rg - target_rg) / target_rg > rg_tolerance / 100.0)
+        return EXIT_FAILURE;
+
+    return 0;
+
+    return 0;
+}