-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmain.cpp
167 lines (139 loc) · 6.88 KB
/
main.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
//
// Created by egor on 1/25/24.
//
#include <iostream>
#include <vector>
#include <chrono>
#include <Eigen/Eigen>
#include <libgran/contact_force/contact_force.h>
#include <libgran/hamaker_force/hamaker_force.h>
#include <libgran/granular_system/granular_system.h>
#include "writer.h"
#include "test/compute_energy.h"
// The driver program is responsible for:
// 1) Instantiating the force models
// 2) Instantiating the step handler
// 3) Initializing the initial positions and velocities
// 4) Initializing sinter bridges and sinter-contact map
// 4) Instantiating granular system
// 5) Performing the time steps
int main() {
// 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 n_thermo_dumps = 10000;
const size_t dump_period = n_steps / n_dumps;
const size_t thermo_dump_period = n_steps / n_thermo_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;
// A sphere
const double r_sphere = 8.0 * r_part;
const Eigen::Vector3d sphere_1_center = Eigen::Vector3d::Zero();
const Eigen::Vector3d sphere_2_center = Eigen::Vector3d::Zero() + 2.5 * Eigen::Vector3d::UnitX() * r_sphere;
const size_t max_num_parts_per_dim = size_t(2.0 * r_sphere / r_part);
// Create sphere 1
for (size_t i = 0; i < max_num_parts_per_dim; i ++) {
double x = sphere_1_center[0] - r_sphere + r_part + 2.0 * r_part * double(i);
for (size_t j = 0; j < max_num_parts_per_dim; j ++) {
double y = sphere_1_center[1] - r_sphere + r_part + 2.0 * r_part * double(j);
for (size_t m = 0; m < max_num_parts_per_dim; m ++) {
double z = sphere_1_center[2] - r_sphere + r_part + 2.0 * r_part * double(m);
if ((Eigen::Vector3d{x, y, z} - sphere_1_center).norm() < r_sphere) {
x0.emplace_back(x, y, z);
v0.emplace_back(0.5, 0, 0);
}
}
}
}
// Create sphere 2
for (size_t i = 0; i < max_num_parts_per_dim; i ++) {
double x = sphere_2_center[0] - r_sphere + r_part + 2.0 * r_part * double(i);
for (size_t j = 0; j < max_num_parts_per_dim; j ++) {
double y = sphere_2_center[1] - r_sphere + r_part + 2.0 * r_part * double(j);
for (size_t m = 0; m < max_num_parts_per_dim; m ++) {
double z = sphere_2_center[2] - r_sphere + r_part + 2.0 * r_part * double(m);
if ((Eigen::Vector3d{x, y, z} - sphere_2_center).norm() < r_sphere) {
x0.emplace_back(x, y, z);
v0.emplace_back(-0.5, 0, 0);
}
}
}
}
std::cout << "Number of particles: " << x0.size() << std::endl;
// Initialize the remaining buffers
theta0.resize(x0.size());
omega0.resize(x0.size());
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
// 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<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(x0,
v0, theta0, omega0, 0.0, Eigen::Vector3d::Zero(), 0.0,
step_handler_instance, binary_force_functors, unary_force_functors);
// Output stream for data
std::ofstream ofs("energies.dat");
ofs << "t\tE_trs\tE_rot\tE_tot\tP\tL" << std::endl;
auto start_time = std::chrono::high_resolution_clock::now();
for (size_t n = 0; n < n_steps; n ++) {
if (n % dump_period == 0) {
std::cout << "Dump #" << n / dump_period << std::endl;
write_particles("run", system.get_x(), system.get_theta(), r_part);
}
if (n % thermo_dump_period == 0) {
ofs << double(n) * dt << "\t"
<< compute_translational_kinetic_energy(system.get_v(), mass) << "\t"
<< compute_rotational_kinetic_energy(system.get_omega(), inertia) << "\t"
<< compute_total_kinetic_energy(system.get_v(), system.get_omega(), mass, inertia) << "\t"
<< compute_linear_momentum(system.get_v(), mass) << "\t"
<< compute_angular_momentum(system.get_x(), system.get_v(), mass, system.get_omega(), inertia) << std::endl;
}
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;
return 0;
}