Browse Source

Measure mean time per step of single-dimer simulation

Peter Colberg 1 year ago
parent
commit
bef77e2274

+ 58
- 0
examples/single_dimer/benchmark/config.lua View File

@@ -0,0 +1,58 @@
1
+-- Edge lengths of simulation domain.
2
+L = {50, 50, 50}
3
+-- Periodic boundary conditions.
4
+periodic = {true, true, true}
5
+-- Number of dimer spheres.
6
+Nd = 2
7
+-- Diameters of dimer spheres.
8
+diameter = {C = 4.0, N = 8.0}
9
+-- Solvent number density.
10
+rho = 9.0
11
+-- Number of solvent particles.
12
+Ns = math.floor(rho*L[1]*L[2]*L[3])
13
+-- Cutoff radius of potential in units of σ.
14
+r_cut = 2^(1/6)
15
+-- Minimum number of steps between solvent particle sorts.
16
+sort = {interval = 500}
17
+-- Solvent temperature.
18
+temp = 1/6
19
+-- Masses of neutrally buoyant dimer spheres.
20
+mass = {C = rho * (math.pi/6) * diameter.C^3, N = rho * (math.pi/6) * diameter.N^3}
21
+-- Distance of dimer spheres.
22
+bond = (diameter.C + diameter.N)/2 * r_cut
23
+-- Potential well depths of dimer spheres with solvent particles.
24
+epsilon = {CA = 1.0, CB = 1.0, NA = 1.0, NB = 0.1}
25
+-- Potential diameters.
26
+sigma = {C = diameter.C/2, N = diameter.N/2}
27
+-- Forward and backward reaction of solvent particles.
28
+reaction = {radius = math.max(L[1], L[2], L[3])/2}
29
+-- Verlet neighbour list skin.
30
+skin = 1.0
31
+-- Number of placeholders per dimer-sphere neighbour list.
32
+nlist_size = {ds = 5000}
33
+-- Number of bins per dimension for dimer spheres.
34
+nbin = {math.floor(L[1]/2), math.floor(L[2]/2), math.floor(L[3]/2)}
35
+-- Number of placeholders per bin.
36
+bin_size = 2
37
+-- Multi-particle collision dynamics.
38
+mpcd = {interval = 50, bin_size = 50}
39
+-- Integration time-step.
40
+timestep = 0.01
41
+-- Measurement interval.
42
+interval = 10000
43
+-- Number of steps.
44
+nsteps = 20*interval
45
+-- Initial measurement step.
46
+initial = 10*interval
47
+-- OpenCL platform.
48
+platform = 1
49
+-- OpenCL device per platform.
50
+device = 1
51
+-- Random number generator seed.
52
+seed = "0x0788baa97a44f354"
53
+-- JSON output filename.
54
+output = "single_dimer.json"
55
+-- JSON file author.
56
+author = {name = "Berni Julian Alder", email = "bjalder@example.org"}
57
+-- JSON file creator.
58
+creator = {name = "single_dimer.lua", version = "1.0"}

+ 220
- 0
examples/single_dimer/benchmark/single_dimer.lua View File

@@ -0,0 +1,220 @@
1
+#!/usr/bin/env luajit
2
+------------------------------------------------------------------------------
3
+-- Measure mean time per step of many-dimer simulation.
4
+-- Copyright © 2013–2017 Peter Colberg.
5
+-- Distributed under the MIT license. (See accompanying file LICENSE.)
6
+------------------------------------------------------------------------------
7
+
8
+-- Check uses of undeclared global variables.
9
+require("strict")
10
+
11
+local nm = {
12
+  config = require("nanomotor.config"),
13
+  box = require("nanomotor.box"),
14
+  domain = require("nanomotor.domain"),
15
+  position = require("nanomotor.position"),
16
+  velocity = require("nanomotor.velocity"),
17
+  properties = require("nanomotor.properties"),
18
+  sort = require("nanomotor.sort"),
19
+  neigh = require("nanomotor.neigh"),
20
+  lj = require("nanomotor.lj"),
21
+  rattle = require("nanomotor.rattle"),
22
+  verlet = require("nanomotor.verlet"),
23
+  mpcd = require("nanomotor.mpcd"),
24
+  observe = require("nanomotor.observe"),
25
+  random = require("nanomotor.random"),
26
+  stats = require("nanomotor.statistics"),
27
+}
28
+
29
+local cl = require("opencl")
30
+local ffi = require("ffi")
31
+local syscall = require("syscall")
32
+local json = require("cjson")
33
+
34
+-- Read configuration file.
35
+local args = nm.config({}, arg[1])
36
+
37
+-- Seed random number generator.
38
+math.randomseed(args.seed)
39
+
40
+-- Create OpenCL context on chosen OpenCL device.
41
+local platform = cl.get_platforms()[args.platform]
42
+print(platform:get_info("version"))
43
+print(platform:get_info("name"))
44
+print(platform:get_info("vendor"))
45
+print(platform:get_info("extensions"))
46
+local device = platform:get_devices()[args.device]
47
+print(device:get_info("name"))
48
+print(device:get_info("vendor"))
49
+print(device:get_info("extensions"))
50
+local context = cl.create_context({device})
51
+local queue = context:create_command_queue(device)
52
+
53
+-- Create simulation box with given boundary conditions.
54
+local box = nm.box(args)
55
+
56
+-- Allocate zero-filled particle buffers.
57
+local dom = nm.domain(queue, box, args)
58
+
59
+-- Place sphere-dimer at random position and orientation.
60
+do
61
+  local rd, spd, idd = dom.rd, dom.spd, dom.idd
62
+  local L, bond = box.L_global, args.bond
63
+  for i = 0, dom.Nd-1, 2 do
64
+    rd[i].x = math.random() * L[1]
65
+    rd[i].y = math.random() * L[2]
66
+    rd[i].z = math.random() * L[3]
67
+    local x, y, z = nm.random.sphere()
68
+    x = rd[i].x + bond*x
69
+    y = rd[i].y + bond*y
70
+    z = rd[i].z + bond*z
71
+    rd[i+1].x, rd[i+1].y, rd[i+1].z = box.minimage(x, y, z)
72
+    spd[i].s0, spd[i+1].s0 = 0, 1
73
+    idd[i], idd[i+1] = i, i+1
74
+  end
75
+end
76
+
77
+-- Randomly place solvent particles excluding dimer-sphere volumes.
78
+nm.position.random(dom, args)()
79
+-- Pick solvent velocities from Maxwell-Boltzmann distribution.
80
+nm.velocity.gaussian(dom, args)()
81
+
82
+-- Computes thermodynamic properties.
83
+local properties = nm.properties(dom)
84
+
85
+-- Create observables.
86
+local observables = {}
87
+
88
+-- Equilibrate system and measure mean time per step.
89
+local acc = nm.stats.accum2()
90
+do
91
+  local initial, nsteps, interval = args.initial, args.nsteps, args.interval
92
+  local self = {step = 0}
93
+
94
+  function self.observe()
95
+    while self.step <= initial do
96
+      coroutine.yield()
97
+      local p = properties()
98
+      io.write("#", self.step, " ", p.en_tot, "\n")
99
+      io.flush()
100
+      self.step = self.step + interval
101
+    end
102
+    while self.step <= nsteps do
103
+      queue:finish()
104
+      local t1 = syscall.clock_gettime("monotonic")
105
+      coroutine.yield()
106
+      queue:finish()
107
+      local t2 = syscall.clock_gettime("monotonic")
108
+      acc((t2.time - t1.time) / interval)
109
+      local p = properties()
110
+      io.write("#", self.step, " ", p.en_tot, "\n")
111
+      io.flush()
112
+      self.step = self.step + interval
113
+    end
114
+  end
115
+
116
+  table.insert(observables, self)
117
+end
118
+
119
+-- Initialise observables.
120
+local observe = nm.observe(observables)
121
+
122
+-- Run simulation.
123
+local sort = nm.sort(dom)
124
+local neigh = nm.neigh(dom, args)
125
+local lj = nm.lj(dom, neigh, args)
126
+local rattle = nm.rattle(dom, args)
127
+local verlet = nm.verlet(dom, args)
128
+local mpcd = nm.mpcd(dom, args.mpcd)
129
+
130
+neigh.pre_update()
131
+sort.update()
132
+neigh.update()
133
+lj.update()
134
+neigh.post_update()
135
+verlet.post_neigh()
136
+rattle.post_neigh()
137
+lj.post_update()
138
+
139
+local mpcd_interval = args.mpcd.interval
140
+local sort_interval = args.sort.interval
141
+local sort_step = dom.step
142
+
143
+while true do
144
+  local step = observe(dom.step)
145
+  if step == nil then break end
146
+  dom.step = dom.step + 1
147
+  verlet.integrate()
148
+  rattle.integrate()
149
+  while true do
150
+    lj.post_integrate()
151
+    if dom.step % mpcd_interval == 0 then
152
+      verlet.stream()
153
+    end
154
+    if verlet.check_neigh() or rattle.check_neigh() then
155
+      if verlet.check_stream() then
156
+        verlet.stream()
157
+      end
158
+      neigh.pre_update()
159
+      if dom.step >= sort_step + sort_interval then
160
+        sort.update()
161
+        sort_step = dom.step
162
+      end
163
+      neigh.update()
164
+      neigh.post_update()
165
+      verlet.post_neigh()
166
+      rattle.post_neigh()
167
+    end
168
+    lj.update()
169
+    if dom.step % mpcd_interval == 0 then
170
+      verlet.finalize()
171
+      mpcd.collide()
172
+      verlet.post_collide()
173
+    end
174
+    if dom.step == step then break end
175
+    dom.step = dom.step + 1
176
+    verlet.integrate()
177
+    lj.post_update()
178
+    rattle.finalize()
179
+    rattle.integrate()
180
+  end
181
+  if verlet.check_stream() then
182
+    verlet.stream()
183
+  end
184
+  if verlet.check_finalize() then
185
+    verlet.finalize()
186
+  end
187
+  lj.post_update()
188
+  rattle.finalize()
189
+end
190
+
191
+local result = {
192
+  parameters = args,
193
+  platform = {
194
+    version = platform:get_info("version"),
195
+    name = platform:get_info("name"),
196
+    vendor = platform:get_info("vendor"),
197
+    extensions = platform:get_info("extensions"),
198
+  },
199
+  device = {
200
+    version = device:get_info("version"),
201
+    name = device:get_info("name"),
202
+    vendor = device:get_info("vendor"),
203
+    driver_version = device:get_info("driver_version"),
204
+    extensions = device:get_info("extensions"),
205
+    address_bits = device:get_info("address_bits"),
206
+    global_mem_size = device:get_info("global_mem_size"),
207
+  },
208
+  time = {
209
+    count = nm.stats.count(acc),
210
+    mean = nm.stats.mean(acc),
211
+    error = nm.stats.sem(acc),
212
+  },
213
+}
214
+
215
+-- Write JSON output file.
216
+do
217
+  local f = assert(io.open(args.output, "w"))
218
+  assert(f:write(json.encode(result)))
219
+  assert(f:close())
220
+end

Loading…
Cancel
Save