I am running into issues while trying to write a new compute style.
The compute outputs a global vector that is supposed to contain the
sum of pairwise forces between a set of atoms the line joining each
pair of which, intersects a plane. (I am trying to discretize the method
of planes (MOP) approach to find the pressure distribution on a plane).
I am using fix ave/time to time average this global vector that my compute
returns. What I see is nearly the expected pressure distribution together
with strange, unexpected spikes. I suspect that this is possibly associated
with memory allocation for this global vector done by fix ave/time.
I think this because I see that the values of the vector obtained using the two
different averaging windows, on post processing to correspond to averages
over the same set of timesteps, give different values (different enough to not
appear like round off error).
The first set of values in file2 should be obtained on averaging
the first 10 sets of values in file1, but that is not the case.
And in any case, the spurious spikes produced by the ave/time
happens at random locations based on the window in ave/time,
which is making me think that it is not related to the physics of
the problem, but has to do with the memory/implementation.
I can send more details on the compute (and will send a pull request
with it once I have it working correctly) if that would be helpful
in diagnosing the problem. I am allocating memory
for the global vector (which could have a size upto 100*8 bytes)
inside the compute. Is that what is causing the problem?
I would really, really appreciate any input you may have.
Thank you very much and happy new year!
I am running into issues while trying to write a new compute style.
The compute outputs a global vector that is supposed to contain the
sum of pairwise forces between a set of atoms the line joining each
pair of which, intersects a plane. (I am trying to discretize the method
of planes (MOP) approach to find the pressure distribution on a plane).
I am using fix ave/time to time average this global vector that my compute
returns. What I see is nearly the expected pressure distribution together
with strange, unexpected spikes. I suspect that this is possibly associated
with memory allocation for this global vector done by fix ave/time.
I think this because I see that the values of the vector obtained using the
two
different averaging windows, on post processing to correspond to averages
over the same set of timesteps, give different values (different enough to
not
appear like round off error).
The first set of values in file2 should be obtained on averaging
the first 10 sets of values in file1, but that is not the case.
And in any case, the spurious spikes produced by the ave/time
happens at random locations based on the window in ave/time,
which is making me think that it is not related to the physics of
the problem, but has to do with the memory/implementation.
I can send more details on the compute (and will send a pull request
with it once I have it working correctly) if that would be helpful
in diagnosing the problem. I am allocating memory
for the global vector (which could have a size upto 100*8 bytes)
inside the compute. Is that what is causing the problem?
my crystal is completely worn out, so i cannot answer this kind of
question specifically.
in fact, you yourself are in the best position to answer your own question.
just have your compute output the raw global vector data to a file and
check it rather than using
if you suspect, memory allocation is a problem, you can try running
under valgrind's memory checker or use gcc's address sanitizer.
You don’t need to worry about fix ave/time having a memory issue,
since it has been in LAMMPS for a long time. Your new compute
has to allocate the global vector it produces and declare its
length so that fix ave/time can use it. It should normally be a static
length (defined in the constructor), and thus only be allocated once
(also in the constructor), and only freed in the destructor when
the compute is deleted (e.g. end of run).
I made a compute to do just that with stresses and heat fluxes for a paper several months ago. If your plane is spatially defined and your material at a high enough temperature spikes in the data can happen from the kinetic value of fluxes as two atoms that used to intersect in a timestep no longer intersect the plane in the next.
I responded with this message but I think it wasn’t sent on the same thread, so I am sending it again:
Thank you very much for your responses!
@Adrian: I did suspect a method-related instability causing the spikes
but I didn’t think of the reasoning you suggested.
It sounds quite plausible and I would be extremely grateful if you could
point me to the right literature that suggests this possibility. I read
some papers that suggested that adding up the kinetic term nullifies
the spikes but that doesn’t apply in my case since I checked that there
are no atom crossings.
@Steve: I am dynamically allocating memory for the vector (of a
static length) in the constructor of the compute style and it gets
deallocated in the destructor. There is no other vector/array that is dynamically allocated inside the compute.
@Axel: Thanks for your suggestions! I will get back to you after trying them