-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmd_lib_SoA.c
151 lines (135 loc) · 5.13 KB
/
md_lib_SoA.c
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
/*
Name: Zefan Wu
Student Number: u6474528
Course: COMP3320
Assignment Number: 1
Name of this file: md_lib_SoA.c
I declare that the material I am submitting in this file is entirely my own
work. I have not collaborated with anyone to produce it, nor have I copied it,
in part or in full, from work produced by someone else. I have not given access
to this material to any other student in this course.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>
#include <sys/types.h>
#include <sys/times.h>
#include <unistd.h>
#include <time.h>
#ifdef __SSE2__
#include <immintrin.h>
#endif
// Struct containting x, y, z components of a vector using SoA format
typedef struct _vector_array {
double *restrict x;
double *restrict y;
double *restrict z;
} VectorArray;
// Struct containting the position, velocity, acceleration along with the number
// of atoms n for a set of atoms
typedef struct _atoms {
VectorArray position;
VectorArray velocity;
VectorArray acceleration;
VectorArray old_acceleration;
int n;
} Atoms;
void update_position(double dt, Atoms *restrict atoms) {
for (int r1 = 0; r1 < atoms->n; r1++)
{
atoms->position.x[r1] = atoms->position.x[r1] + dt * (atoms->velocity.x[r1] + dt * atoms->acceleration.x[r1] / 2);
atoms->position.y[r1] = atoms->position.y[r1] + dt * (atoms->velocity.y[r1] + dt * atoms->acceleration.y[r1] / 2);
atoms->position.z[r1] = atoms->position.z[r1] + dt * (atoms->velocity.z[r1] + dt * atoms->acceleration.z[r1] / 2);
}
return;
}
void calculate_pe(Atoms *restrict atoms, double *potential_energy) {
// calculate new PE
for (int r1 = 0; r1 < atoms->n - 1; r1++)
{
double r1_x = atoms->position.x[r1];
double r1_y = atoms->position.y[r1];
double r1_z = atoms->position.z[r1];
for (int r2 = r1 + 1; r2 < atoms->n; r2++)
{
double r2_x = atoms->position.x[r2];
double r2_y = atoms->position.y[r2];
double r2_z = atoms->position.z[r2];
double dis_r1r2 = sqrt(pow(r1_x - r2_x, 2) + pow(r1_y - r2_y, 2) + pow(r1_z - r2_z, 2));
*potential_energy += 1 / pow(dis_r1r2, 12) - 2 / pow(dis_r1r2, 6);
}
}
return;
}
void update_force(Atoms *restrict atoms) {
for (int r1 = 0; r1 < atoms->n; r1++)
{
double r1_x = atoms->position.x[r1];
double r1_y = atoms->position.y[r1];
double r1_z = atoms->position.z[r1];
// record old acceleration
atoms->old_acceleration.x[r1] = atoms->acceleration.x[r1];
atoms->old_acceleration.y[r1] = atoms->acceleration.y[r1];
atoms->old_acceleration.z[r1] = atoms->acceleration.z[r1];
atoms->acceleration.x[r1] = 0;
atoms->acceleration.y[r1] = 0;
atoms->acceleration.z[r1] = 0;
for (int r2 = 0; r2 < atoms->n; r2++)
{
if(r1 != r2)
{
double r2_x = atoms->position.x[r2];
double r2_y = atoms->position.y[r2];
double r2_z = atoms->position.z[r2];
// distance between r1 and r2
double dis_r1r2 = sqrt(pow(r1_x - r2_x, 2) + pow(r1_y - r2_y, 2) + pow(r1_z - r2_z, 2));
// acceleration(t + dt) of r1 in x, y, z directions
atoms->acceleration.x[r1] += (12 / pow(dis_r1r2, 14) - 12 / pow(dis_r1r2, 8)) * (r1_x - r2_x);
atoms->acceleration.y[r1] += (12 / pow(dis_r1r2, 14) - 12 / pow(dis_r1r2, 8)) * (r1_y - r2_y);
atoms->acceleration.z[r1] += (12 / pow(dis_r1r2, 14) - 12 / pow(dis_r1r2, 8)) * (r1_z - r2_z);
}
}
}
return;
}
void update_volocities(double dt, Atoms *restrict atoms, double *kinetic_energy) {
for (int r1 = 0; r1 < atoms->n; r1++){
// update r1 velocity for x, y, z directions
atoms->velocity.x[r1] = atoms->velocity.x[r1] + dt * (atoms->acceleration.x[r1] + atoms->old_acceleration.x[r1]) / 2;
atoms->velocity.y[r1] = atoms->velocity.y[r1] + dt * (atoms->acceleration.y[r1] + atoms->old_acceleration.y[r1]) / 2;
atoms->velocity.z[r1] = atoms->velocity.z[r1] + dt * (atoms->acceleration.z[r1] + atoms->old_acceleration.z[r1]) / 2;
// calculate the KE
*kinetic_energy += (pow(atoms->velocity.x[r1], 2) + pow(atoms->velocity.y[r1], 2) + pow(atoms->velocity.z[r1], 2)) / 2;
}
return;
}
// Update all atoms by a single timestep of length dt and return the kinetic and
// potential energy of the system
void update_timestep(double dt, Atoms *restrict atoms, double *kinetic_energy,
double *potential_energy) {
*potential_energy = 0;
*kinetic_energy = 0;
/*** YOUR CODE HERE FOR STEP 2 ***/
update_position(dt, atoms);
// calculate new PE
calculate_pe(atoms, potential_energy);
// update new acceleration and record old acceleration
update_force(atoms);
// update velocities
update_volocities(dt, atoms, kinetic_energy);
return;
}
// Put any one time initialization code here - this will be called once at the
// start of your program. You do not need to add code here if you have no
// initialization requirements.
void initial_config(double dt, Atoms *restrict atoms,
double *potential_energy) {
*potential_energy = 0;
/*** YOUR CODE HERE FOR STEP 2 ***/
// calculate intial PE
calculate_pe(atoms, potential_energy);
// update new acceleration and record old acceleration
update_force(atoms);
return;
}