RMPCDMD
poiseuille_flow.f90
Go to the documentation of this file.
1 ! This file is part of RMPCDMD
2 ! Copyright (c) 2015-2016 Pierre de Buyl and contributors
3 ! License: BSD 3-clause (see file LICENSE)
4 
21 
23  use rmpcdmd_module
24  use hdf5
25  use h5md_module
26  use threefry_module
27  use parsetext
28  use iso_c_binding
29  use omp_lib
30  implicit none
31 
32  type(threefry_rng_t), allocatable :: state(:)
33  type(pto) :: config
34 
35  type(cell_system_t) :: solvent_cells
36  type(particle_system_t) :: solvent
37 
38  type(profile_t) :: tz
39  type(histogram_t) :: rhoz
40  type(profile_t) :: vx
41 
42  type(h5md_file_t) :: datafile
43  type(h5md_element_t) :: elem
44  type(h5md_element_t) :: elem_tz, elem_tz_count, elem_vx_count
45  type(h5md_element_t) :: elem_rhoz
46  type(h5md_element_t) :: elem_vx
47  type(h5md_element_t) :: elem_T
48  type(h5md_element_t) :: elem_v_com
49  integer(HID_T) :: box_group
50  integer(HID_T) :: fields_group
51  type(particle_system_io_t) :: solvent_io
52 
53  integer :: i, L(3), error, N, n_threads
54  integer :: rho
55  integer :: N_loop, N_therm
56 
57  double precision :: v_com(3), wall_v(3,2), wall_t(2)
58  double precision :: gravity_field(3)
59  double precision :: T, set_temperature, tau
60  double precision :: alpha
61  logical :: thermostat
62  logical :: do_hydro
63  type(args_t) :: args
64 
65  args = get_input_args()
66  call ptparse(config, args%input_file, 11)
67 
68  n_threads = omp_get_max_threads()
69  allocate(state(n_threads))
70  call threefry_rng_init(state, args%seed)
71 
72  call h5open_f(error)
73 
74  l = ptread_ivec(config, 'L', 3)
75  rho = ptread_i(config, 'rho')
76  n = rho *l(1)*l(2)*l(3)
77 
78  set_temperature = ptread_d(config, 'T')
79  thermostat = ptread_l(config, 'thermostat')
80  do_hydro = ptread_l(config, 'do_hydro')
81 
82  tau = ptread_d(config, 'tau')
83  alpha = ptread_d(config, 'alpha')
84  n_therm = ptread_i(config, 'N_therm')
85  n_loop = ptread_i(config, 'N_loop')
86  gravity_field = 0
87  gravity_field(1) = ptread_d(config, 'g')
88 
89  call solvent% init(n)
90 
91  do i=1, solvent% Nmax
92  solvent% vel(1,i) = threefry_normal(state(1))
93  solvent% vel(2,i) = threefry_normal(state(1))
94  solvent% vel(3,i) = threefry_normal(state(1))
95  end do
96  solvent%vel = solvent%vel*sqrt(set_temperature)
97  v_com = sum(solvent% vel, dim=2) / size(solvent% vel, dim=2)
98  solvent% vel = solvent% vel - spread(v_com, dim=2, ncopies=size(solvent% vel, dim=2))
99 
100  solvent% force = 0
101  solvent% species = 1
102  call solvent% random_placement(l*1.d0, state=state(1))
103 
104  call solvent_cells%init(l, 1.d0, has_walls=.true.)
105  solvent_cells% origin(3) = -0.5d0
106 
107  call solvent_cells%count_particles(solvent% pos)
108 
109  call datafile% create(args%output_file, 'RMPCDMD:poiseuille_flow', &
110  rmpcdmd_revision, 'Pierre de Buyl')
111 
112  call ptkill(config)
113 
114  call tz% init(0.d0, solvent_cells% edges(3), l(3))
115  call rhoz% init(0.d0, solvent_cells% edges(3), l(3))
116  call vx% init(0.d0, solvent_cells% edges(3), l(3))
117 
118  solvent_io%force_info%store = .false.
119  solvent_io%species_info%store = .false.
120  solvent_io%position_info%store = .true.
121  solvent_io%position_info%mode = ior(h5md_linear,h5md_store_time)
122  solvent_io%position_info%step = 1
123  solvent_io%position_info%time = tau
124  solvent_io%image_info = solvent_io%position_info
125  solvent_io%velocity_info = solvent_io%position_info
126  solvent_io%id_info = solvent_io%position_info
127  call solvent_io%init(datafile, 'solvent', solvent)
128 
129  call h5gcreate_f(solvent_io%group, 'box', box_group, error)
130  call h5md_write_attribute(box_group, 'dimension', 3)
131  call h5md_write_attribute(box_group, 'boundary', ['periodic', 'periodic', 'periodic'])
132  call elem% create_fixed(box_group, 'edges', l*1.d0)
133  call h5gclose_f(box_group, error)
134 
135  call h5gcreate_f(datafile%id, 'fields', fields_group, error)
136  call elem_tz% create_time(fields_group, 'tz', tz% data, ior(h5md_time, h5md_store_time))
137  call elem_tz_count% create_time(fields_group, 'tz_count', tz% count, ior(h5md_time, h5md_store_time))
138  call elem_vx% create_time(fields_group, 'vx', tz% data, ior(h5md_time, h5md_store_time))
139  call elem_vx_count% create_time(fields_group, 'vx_count', vx% count, ior(h5md_time, h5md_store_time))
140  call elem_rhoz% create_time(fields_group, 'rhoz', rhoz% data, ior(h5md_time, h5md_store_time))
141  call elem_t% create_time(datafile% observables, 'temperature', t, ior(h5md_time, h5md_store_time))
142  call elem_v_com% create_time(datafile% observables, 'center_of_mass_velocity', v_com, ior(h5md_time, h5md_store_time))
143  call h5gclose_f(fields_group, error)
144 
145  call solvent% sort(solvent_cells)
146 
147  solvent%force(1,:) = gravity_field(1)
148  solvent%force_old(1,:) = gravity_field(1)
149 
150  wall_v = 0
151  wall_t = set_temperature
152  solvent_cells%bc = [ periodic_bc, periodic_bc, bounce_back_bc ]
153  do i = 1, n_therm
154  call mpcd_stream_xforce_yzwall(solvent, solvent_cells, tau, gravity_field(1))
155  call md_vel(solvent, tau)
156  call apply_pbc(solvent, solvent_cells%edges)
157  call solvent_cells%random_shift(state(1))
158  call solvent% sort(solvent_cells)
159  call wall_mpcd_step(solvent, solvent_cells, state, &
160  wall_temperature=wall_t, wall_v=wall_v, wall_n=[rho, rho], thermostat=thermostat, &
161  bulk_temperature=set_temperature, alpha=alpha, keep_cell_v=do_hydro)
162  end do
163 
164  do i = 1, n_loop
165  if (mod(i,100)==0) then
166  write(*,*) i
167  end if
168  call mpcd_stream_xforce_yzwall(solvent, solvent_cells, tau, gravity_field(1))
169  call md_vel(solvent, tau)
170  call apply_pbc(solvent, solvent_cells%edges)
171  call solvent_cells%random_shift(state(1))
172  call solvent% sort(solvent_cells)
173  call wall_mpcd_step(solvent, solvent_cells, state, &
174  wall_temperature=wall_t, wall_v=wall_v, wall_n=[rho, rho], thermostat=thermostat, &
175  bulk_temperature=set_temperature, alpha=alpha, keep_cell_v=do_hydro)
176  v_com = sum(solvent% vel, dim=2) / size(solvent% vel, dim=2)
177 
178 
179  t = compute_temperature(solvent, solvent_cells, tz)
180  call elem_v_com%append(v_com, i, i*tau)
181  call elem_t% append(t, i, i*tau)
182 
183  call compute_rho(solvent, rhoz)
184  call compute_vx(solvent, vx)
185 
186  if (modulo(i, 50) == 0) then
187  call tz% norm()
188  call elem_tz% append(tz% data, i, i*tau)
189  call elem_tz_count% append(tz% count, i, i*tau)
190  call tz% reset()
191  call vx% norm()
192  call elem_vx% append(vx% data, i, i*tau)
193  call elem_vx_count% append(vx% count, i, i*tau)
194  call vx% reset()
195  rhoz% data = rhoz% data / (50.d0 * rhoz% dx)
196  call elem_rhoz% append(rhoz% data, i, i*tau)
197  rhoz% data = 0
198  end if
199 
200  if (modulo(i,100)==0) then
201  call solvent_io%position%append(solvent% pos)
202  call solvent_io%image%append(solvent% image)
203  call solvent_io%id%append(solvent% id)
204  call solvent_io%velocity%append(solvent% vel)
205  call h5fflush_f(datafile%id, h5f_scope_global_f, error)
206  end if
207 
208  end do
209 
210  call elem_tz% close()
211  call elem_tz_count% close()
212  call elem_rhoz% close()
213  call elem_vx% close()
214  call elem_vx_count% close()
215  call elem_t% close()
216  call elem_v_com% close()
217 
218  call datafile% close()
219 
220  call h5close_f(error)
221 
222  write(*,'(a16,f8.3)') solvent%time_stream%name, solvent%time_stream%total
223  write(*,'(a16,f8.3)') solvent%time_step%name, solvent%time_step%total
224  write(*,'(a16,f8.3)') solvent%time_count%name, solvent%time_count%total
225  write(*,'(a16,f8.3)') solvent%time_sort%name, solvent%time_sort%total
226  write(*,'(a16,f8.3)') solvent%time_ct%name, solvent%time_ct%total
227  write(*,'(a16,f8.3)') 'total ', &
228  solvent%time_stream%total + solvent%time_step%total + solvent%time_count%total +&
229  solvent%time_sort%total + solvent%time_ct%total
230 
231 end program poiseuille_flow
program poiseuille_flow
Simulate a forced flow between two plates.