Hi Paul,
many thanks for your answer. I have attached the input file.
It also contains a fix that I have written myself, 'grid_analysis' and that I am using the calculate macroscopic quantities for a grid of cells. I do not know, if you want to make the effort of running it, but in that case I can send the source code of that fix as well.
cheers,
Marco
Paul Crozier wrote:
Channel_r0.1_T2_210x10x10.in (1.38 KB)
Marco,
The main problem was with your temp/rescale command. You switched N and Tstop. See: http://lammps.sandia.gov/doc/fix_temp_rescale.html
Instead of:
fix temp all temp/rescale 2.0 2.0 10 0.0 0.5
You should have put:
fix temp all temp/rescale 10 2.0 2.0 0.0 0.5
The order of the variables may have been changed on you in one of the LAMMPS upgrades.
So since your simulations appear to have been running from T=2 to T=10 rather than at the constant correct T=2 setpoint, the default timestep was too small, and the dynamics were not right. I saw strange periodic behavior when I ran with your settings, but the strange periodic behavior seems to go away when running at T=2. I’ll paste my simple script below.
Also, you had to use “neighbor 5.0 bin” to avoid dangerous reneighborings, but the default setting of “neighbor 0.3 bin” should be fine when running at T=2. And your simulations will run much faster.
What is “atom_style energy”? I think that should be “atom_style atomic”.
Paul
dimension 3
units lj
atom_style atomic
boundary f p p
lattice sc 0.1
region SimBox block -5.0 205.0 0 10 0 10 units box
create_box 1 SimBox
mass 1 1.0
pair_style lj/cut 2.5
pair_modify shift yes
pair_coeff 1 1 1.0 1.0
#neighbor 5.0 bin
neigh_modify delay 4 every 4 check yes
create_atoms 1 SimBox
velocity all create 2.0 34556 mom yes
fix int all nve
fix walls all wall/reflect xlo xhi
fix temp all temp/rescale 10 2.0 2.0 0.0 0.5
#fix temp all temp/rescale 2.0 2.0 10 0.0 0.5
variable b equal vcm(all,x)
variable a loop 10000
label loop
thermo 100
run 1000
print The average x-directional velocity is now: $b
next a
jump wall.in loop
Paul,
This is a bit more lengthy version of the first email.
For the neighbor list "skin" (neighbor <skin> bin), is there somewhere a list with recommondations, for different densities / temperatures?
Regarding the temp/rescale problem: I did not use the most recent version. In my version the temp/rescale was still like
fix ID group-ID temp/rescale Tstart Tstop N window fraction
I should have changed this in the version that I had sent to you. I have no upgraded to the newest version.
The "atom_style energy" is something that I have put in. It just defines an extra field for the potential energy of each atom. This field is updated during the normal force
calculation. I need it, because I work on hybrid MD-Continuum methods. For this I have to calculate cell averages of density, momentum and total energy for the continuum solver. To get the total energy of the cell, I obviously need Epot of each atom.
Naturally, the cell averages are subject fluctuation in time. In order to obtain an averages that are within a limit (like 5%) from the "true" value, one has to choose two parameters:
- Number of atoms within one cell
- the sample frequency (what is the optimal number time step between each measurment)
- time average over how many measurements
The number of atoms for my cases is quite low (100 -1000). This is mainly because I test a lot and want to the system to run fast.
As far as I know the sample frequency should be based on the correlation time of the variable of interest.
To find the optimal sample frequency I have used the block method, that is described in Allen&Tildesley (from a paper by Flyberg et. al . I think).
From the block method one should get a curve that converges to a plateau, which shows a useful sample frequency.
So I applied this method to cell averaged variables(density, momentum), but the result was not clear.
When I tried to find out why, I realised that there were quite low frequent flucutations.
Calculating the the time correlation function directly I obtain a correlation time of around 40,000 time steps. This seemed very long and is not practicle to use, because I would have to run for a very long time in order to calculate reliable averages.
So I performed a discrete fourier transformation to obtain the power spectrum.
I may be wrong, and please correct in that case: For random fluctuations the power spectra should look random as well, and energies should be distributed over the entire frequency range and not only at low frequencies. But here, all the energy is in low frequent waves and not "uniformly" distributed over the spectrum.
I have put together the graphs for the some example runs. My basic idea was to use small free flow cases to test the method. So all the setups are periodic in y and z direction and with "continuum boundary conditions" in x directions. For fluid at rest, I just used walls to confine the fluid, which are later replaced by the continuum boundary conditions. However, the long time fluctuations are also present for peridicity in all directions (so the boundary is not the problem).
gas (r=0.01) at T=2, simulation box based on a lattice of 200x5x5 simulated over 10,000,000 time steps:
1a) velocity of the entire channel: www.cranfield.ac.uk/~ei2989/lammps/vx/Powerapectra_r0.01_T2_all_460x23x23.jpg
1b) velocity of the entire channel: www.cranfield.ac.uk/~ei2989/lammps/vx/Powerapectra_r0.01_T2_cell_460x23x23.jpg
A 2D case of the same system (simulation box based on a lattice of 200x200x5), 1,000,000 time steps:
2) velocity of a cell: www.cranfield.ac.uk/~ei2989/lammps/vx/PowerSpectra_r0.01_T2_vx_cell_465x465x46.jpg
A smaller case with r=0.1 and T=2 , simulation based on a lattice of 100x5x5 time steps:
3a) velocitity of the entire channel: www.cranfield.ac.uk/~ei2989/lammps/vx/PowerSpectra_r0.1_T2_vx_all_215x23x23.jpg
3b) velocitity of a cell in the middle: www.cranfield.ac.uk/~ei2989/lammps/vx/PowerSpectra_r0.1_T2_vx_cell_215x23x23.jpg
The most extrem is 3a), as here almost the entire energy is contained in one wave length. But for all other cases the fluctuations are low frequent.
Marco,
Using the default skin distance is a good first approximation. If you get “dangerous rebuilds”, then increasing the skin distance is one way to get rid of them. I don’t know of a list of recommended skin distances for various T / density values.
I’d certainly recommend updating to the most recent version of LAMMPS since there have been many bug fixes that could possibly affect your results. We have a much harder time helping users who don’t have the current version and have made substantial changes to LAMMPS as you have.
The only thing, that I actually need know is why the fluctuations are
mainly low frequent and how big a system or/and a subcell has to be to
obtain avoid mainly low frequent fluctuations.
Since the walls are artificial and introduce discontinuities into the net momentum of the system, you very well may get strange net velocity behavior. This would be minimized somewhat by using a bigger system with more particles.
Paul
