Nbody structs #31
Nbody structs #31
Conversation
Codecov Report
@@ Coverage Diff @@
## master #31 +/- ##
===========================================
+ Coverage 70.37% 82.44% +12.07%
===========================================
Files 4 10 +6
Lines 108 336 +228
===========================================
+ Hits 76 277 +201
- Misses 32 59 +27
Continue to review full report at Codecov.
|
src/nbody_base.jl
Outdated
| using StaticArrays, DiffEqBase | ||
| #Structures for bodies and systems under symulations | ||
|
|
||
| # Fields for position, velocity and mass of a particle should be required for every descendant |
There was a problem hiding this comment.
This might be a good use of https://github.com/tbreloff/ConcreteAbstractions.jl . We can fork and register if it makes sense. Or @def
There was a problem hiding this comment.
That is actually cool! Good, I'll check what I can do with @def first.
src/nbody_base.jl
Outdated
| function ode_system!(du, u, p, t) | ||
| du[:, 1:n] = @view u[:, n+1:2n]; | ||
| @inbounds for i=1:n | ||
| du[:, n+i] .= acceleration(u[:, 1:n], i, n, simulation.sysmem, simulation.external_electric_field); |
There was a problem hiding this comment.
Make acceleration non-allocating.
src/nbody_base.jl
Outdated
| external_gravitational_field) | ||
| end | ||
|
|
||
| # Function, scpecifically designed for a system of ChargedParticles |
There was a problem hiding this comment.
I don't get this part. What's going on here? Are you using different acceleration functions based on the existence of some fields and not other?
There was a problem hiding this comment.
This is unclear, I agree. No, not based on the presence or absence of external fields.
Actually, the second implementation is wrong - I want to pass NBodySystem into acceleration function and there to decide what forces act on i-th particle. In addition to the forces from other particles the external fields also affect the i-th particle. But then I realized, that I need the positions of bodies. So here is a mess, I agree.
The better definition of this function would be something like this:
function acceleration(
u,
i,
system :: NBodySystem,
external_electric_field,
external_magnetic_field,
external_gravitational_field)
+endThere was a problem hiding this comment.
Looks like passing the NBodySimulation could make sense?
There was a problem hiding this comment.
I already pass almost all its body, yes. Right now passing NBodySimulation makes sence, indeed. I'll think about it: there might be more properties and fields for simulation. By "simulation" I understand conditions of an experiment, but not properties of the system under test/experiment/simulation.
src/nbody_base.jl
Outdated
| end | ||
|
|
||
| # Function as this one can be used in a force-oriented approach of calculations | ||
| function acceleration( |
There was a problem hiding this comment.
Is there a way at the top level to have the potential be fully user-definable?
There was a problem hiding this comment.
The potential of what? Just to make it clear)
As I understand, during an arbitrary n-body simulation, a user will have a set of bodies/particles with some physical properties, which are significant in his/her situation. Then, for i-th body in the system user can define a potential function which can be used for calculation of the force acting on this i-th body (as the gradient of the potential).
Right now the suggestion is to write a custom acceleration function and pass it into one of the NBodySimulation constructors. For handling this, in the code where converting NBodySimulation to DEProblem takes place, it is necessary to choose a standard acceleration or a custom one, in case it was defined by a user.
Though it doesn't make big difference, whether we work with forces or potentials, I understand that potentials might be more convenient in some situation. I'll prepare interfaces for them too.
|
A general comment before I go through the code. Fundamentally, any N-body system is driven by a set of potentials (which determine the forces). In the general case (which we should support) potentials can depend on (in order of increasing rarity) coordinates, time, velocities and accelerations. What the potentials depend on fundamentally restricts our choice of solvers. For this package to be an integral part of DiffEq, we need the API to communicate these restrictions. Do the existing DiffEq solvers have traits for this sort of thing? |
|
There aren't traits, but traits could be added to the solvers so that way they error when some solver assumption is violated. Right now they pass through. For example, many RKN methods make a documented velocity-independence assumption in some part, but if you violate this it just runs with order loss. Expanding the trait system we have is never bad, just takes man hours. |
|
Most of the existing RKNs will still work, though, in the sense of producing results that are correct, just at a lower order of accuracy, right? |
yes. Actually, all of them will, likely with order 1. |
|
Ok, then any traits we introduce for this should probably only produce warnings, not outright error out. |
|
Okay, having gone through the code I think the general idea of it is okay, although many of the examples probably bear little resemblance to how they're going to look in the end ( |
+ Periodic boundary conditions + Lennard-Jones potential
|
@ChrisRackauckas, I haven't corrected all the places afrter your review yet, but I surely keep them in mind. |
|
No worries |
- Prepared functions for temperature as an average kinetic energy of particles. - Wrote a plot recipe for simulation results (right now it is a scatter plot). - Next step of refactoring: some files rearrangement. - Added a tests.
src/nbody_base.jl
Outdated
| @@ -207,7 +185,41 @@ function pairwise_lennard_jones_acceleration(rs, | |||
| return -24*system.lj_parameters.ϵ/system.lj_parameters.σ*accel | |||
| end | |||
|
|
|||
| function temperature(result :: SimulationResult, time :: Float64) | |||
| function get_velocity(result :: SimulationResult, time, i = 0) | |||
There was a problem hiding this comment.
this doesn't look right
There was a problem hiding this comment.
oh, probably not done here yet.
There was a problem hiding this comment.
I am confused, what did you mean? I have two variants:
timeandiinget_velocityare without type assertioniinget_velocityrepresents particle for which I obtain the velocity. And if i=0, then get velocity values for all the particles in a system - probably it is better to have two specific functions.
There was a problem hiding this comment.
Also, the name result for a function argument is a bad choice.
- Codded an iterator for animation simulation results - Added representation of some objects via Base.show() - Changed implementation of PBC according to Allen
- created a file structure for the interface types - refactored code - changed the gravitational simulation types for the new interface - added example of a custom potential usage
| end | ||
|
|
||
| # this should be a method for integrators designed for the SecondOrderODEProblem (It is worth somehow to sort them from other algorithms) | ||
| function run_simulation(s::NBodySimulation, alg_type::Union{VelocityVerlet,DPRKN6,Yoshida6}, args...; kwargs...) |
There was a problem hiding this comment.
There are too many algorithms in JuliaDiffEq) It is cool, but I am not sure which of them would be good for N-Body simulations. VelocityVerlet is the most mentioned one. I am just not acquainted with others yet, that is why I've chosen these three for a while.
There was a problem hiding this comment.
Don't choose any. Just leave it open without a type parameter. One of the interesting things will be to test different algorithms. This has no effect on performance anyways.
There was a problem hiding this comment.
In fact, I put some algorithms here in order to distinguish when to transform an N-body-system into ODEProblem or into SecondOrderODEProblem. Via the multiple dispatch mechanism. I could group all the appropriate algorithms for SecondOrderODEProblem under a particular parent, though I am not sure if it is ok for the current structure of JuliaDiffEq. Another approach is to add an argument like transformaton_order::Integer = 2.
src/nbody_basic_potentials.jl
Outdated
| end | ||
|
|
||
| struct LennardJonesParameters <: PotentialParameters | ||
| ϵ::AbstractFloat |
There was a problem hiding this comment.
All variables in a struct must have a concrete type. You can do this by parámetrising the type.
There was a problem hiding this comment.
You are right, it would be better for performance (I have just read the Performance tips again).
Ok.
- corrected types of fields in structures - small refactoring
|
I faced the following behavior of systems under simulations. Sometimes calculations run fast, but sometimes they slow down dramatically. In the animation particles are shown not with their real size: markers show just centres of the particles. |
|
Can you do that last plot in log scale? |
|
Ok, but for another round of simulations: I lost those previous data due to a sudden crash. |
|
what you're seeing is instability. The cause of this instability is dependent on the speed of how growth which would be seen in log scale. |
- reworked PBC handling - tried to add ManifoldProjection constraints
|
After having found what slowed down calculations, I still struggle with instabilities. Here I tried to simulate the behavior of 4 particles with zero initial velocities, but at such distance between them, that corresponds to the average distance of that liquid argon molecules: As you can see, at the beginning particles start to repeal from each other, but then the upper particle goes out of the cell borders and its image appears at the bottom of the cell. This image pushed away the particle laying above and then a chain reaction takes place. I tried to use @ChrisRackauckas, @dextorious can you advise me what to try next? My variants are thermostating and/or controlling velocities. |
|
How many steps is it taking before diverging like that? |
|
There are 224 steps in the solution. I tried to calculate it in real and reduced units: all the same. |
|
Try it on this small time span with arbitrary precision units. |
|
If I understood correctly: here is total energy at time span (2.0e-12, 3.0e-12) beginning from the 200th time step: |
|
@ChrisRackauckas, probably I didn't understand you correctly. Did you mean simulate the system on the time span beginning when diverging happens? If I run simulation from, say, the 200th time step, it will break after 224 steps again. So, it seems, that there is no difference in what is the starting time of simulations |
|
What happens if you increase the initial inter-particle distance slightly (say, by 20% or so) and increase the size of the box by 2-3x? Does the instability only appear after the periodic boundary conditions are applied? Thermostatting is not a prerequisite of stable LJ fluid dynamics. An intrusive thermostat might well hide the instability, but that's just masking whatever issue is causing this - we need to find the root cause. |
|
Do the simulation from the beginning but using bigfloats for state and time. There is a chance that this is really slow (since it's arbitrary precision) and doesn't diverge, so make sure you solve it on a short timespan (but long enough to find issues). That would rule out catestrophic cancellation as the issue. Also reduce the tolerance a lot, to like |
|
Are you sure you are calculating the interaction between particle i and the closest version of particle j accounting for the PBCs? From a quick scan of the code I am not sure line 99 of src/nbody_simulation_result.jl changes anything of note, so particles only start to effect particles on the other side of the box when they jump over. Hope this helps. |
|
@jgreener64 , you are right! Right now It works incorrectly. Thanks. |
|
Ok, the problem was in PBC. I think I repaired it. Here is the plot for the temperature of the system of 1000 particles (T = 120.0 K, σ = 3.4e-10 m, τ = 1e-14 s) for 400 steps: And here is the energy graphic: @ChrisRackauckas , @dextorious , I would like to finish and close this pull request as an accomplished step. How do you think, what else should be done for that? |
|
LGTM. The only reservation I have is that the old NBodyProblem was deleted but the replacement isn't feature-complete yet (i.e. it cannot take an arbitrary potential function) so there's currently no deprecation path. It should be high on the priorities to get some replacement in this newer API. I'll give some time for @dextorious to review and then we can merge. |
|
@ChrisRackauckas , actually |
| system::sType | ||
| tspan::Tuple{Float64,Float64} | ||
| boundary_conditions::bcType | ||
| external_electric_field |
There was a problem hiding this comment.
the way that there are applied doesn't cause performance issues? If it's just a single modifying call then it's fine. From what I can see it looks fine?
There was a problem hiding this comment.
Right now it isn't implemented: I left it here just to show where I want to set these fields up. In the future, I think, these fields will depend on time and coordinates
|
Okay, I've gone through the code and profiled the argon test a bit. I didn't find any major bugs, but it's pretty clear that further optimization will be in order before we can run this on real systems of appreciable size. The periodic boundaries seem to work fine now, good job on solving that. On what comes next, my personal view is that it probably makes sense to implement one of the simpler water models (SPC/E or TIP3, for instance) before investing serious effort in optimization because: 1) it involves considerably more functionality and may expose some API issues that would require refactoring and 2) it's much easier to benchmark, since every major MD package out there has published performance numbers on standard water boxes. Still, this is debatable. If you find that performance issues inhibit your productivity too much, it may make sense to optimize sooner. Thermostats would then be the next piece on the agenda, but the simpler thermostats don't really have a significant impact on the performance, so those can be added at any point. As for the merge decision, that's up to whatever standard @ChrisRackauckas wants to maintain for the master branches of DiffEq packages. It's pretty clear that this API shouldn't be considered stable and some elements will probably change substantially as we work towards more complex problems, but it does work and, in my view, represents a great start. |
|
I think it's always best to do smaller PRs and iterate over time. This is fine on master. Lots of unfinished stuff can be in master and it's only "released" when documented. |
CustomPotentialParametersand in test/nbody_custom_potential_test.jl