Rework extract functions.

[Joroks]

I've reworked my previous PR with a more focused scope. It's still quite large so I've tried my best to make it as polished as possible to make reviewing it easier / reduce the amount of revisions needed for it to be merged.

I do feel a bit sorry for comming out of nowhere and suddendly putting this huge amount of extra workload on your table 😅. I'm really gratefull for the work you've put into this thus far! Please take as much time as you need to review this.

[Joroks]

I would almost consider not exporting them.

And then name then NONE instead of LAMMPS_NONE and STYLE_GLOBAL instead of LMP_STYLE_GLOBAL.
Then it's consitent in user code as LAMMPS.NONE and LAMMPS.STYLE_GLOBAL

@vchuravy

I've quickly implemented it in the PR to get a feel for how it would look like and I'm not quite sure if I like this change.

It works well enough for _LMP_DATATYPE but that's about it. The constants for the other enums get really long to type out. I think the better solution is to do this:

const LAMMPS_NONE = _LMP_DATATYPE{API.LAMMPS_NONE}()
...
const TYPE_SCALAR = _LMP_TYPE{API.LMP_TYPE_SCALAR}()
...
const VAR_EQUAL = _LMP_VARIABLE{API.LMP_VAR_EQUAL}()
...
const STYLE_GLOBAL = API.LMP_STYLE_GLOBAL

and export all of them. I originaly just kept the naming from API.jl for _LMP_STYLE_CONST because I wasn't wrapping them in TypeEnum. But removing the LMP_ at the begining would give us a nice and consistent naming scheme.

[vchuravy]

the better solution is to do this

I think that's fine. Better than exporting the STRING

[vchuravy]

Besides the two small nits, this looks great to me.

I would also like to invite you as a collaborator to this repository. The expectation would still be to use PRs and to wait/ask for review before merging large changes, but by now you rewritten most of the code here :)

[Joroks]

I would also like to invite you as a collaborator to this repository. The expectation would still be to use PRs and to wait/ask for review before merging large changes, but by now you rewritten most of the code here :)

I'll gladly accept this invitation!

[vchuravy]

copy should imply own? We can almost never use own=true since we need to use lammps_free for the deallocation

[Joroks]

As far as I can tell, we don't need to use lammps_free for deallocation. It just calls free internally without doing anything else. Thats the same thing that would happen If we take ownership of the memory if I'm not mistaken. I've taken this directly from the lammps source code:

void lammps_free(void *ptr)
{
  free(ptr);
}

The own keyword is used for the extract_variable method, specifically when it's a Atom style variable:

  if lmp_variable == VAR_ATOM
      ndata = extract_setting(lmp, "nlocal")
      ptr = _reinterpret(LAMMPS_DOUBLE, void_ptr)
      # lammps expects us to take ownership of the data
      return _extract(ptr, ndata; own=true)
  end

because the documentation states:

For atom-style variables, the return value is a pointer to an allocated block of storage of double of the length
atom->nlocal. Since the data returned are a copy, the location will persist, but its content will not be updated in case the variable is re-evaluated. To avoid a memory leak, this pointer needs to be freed after use in the calling program.

That means that the value is already a copy and copying it again would be superfluous.

I've also (crudely) verified, that this does indeed not cause a memory leak by running the following script:

lmp = LMP()

command(lmp, """
         atom_modify map yes
         region cell block 0 3 0 3 0 3
         create_box 1 cell
         lattice sc 1
         create_atoms 1 region cell
         mass 1 1
        
         variable var atom id
    """)

for _ in 1:10^7
    extract_variable(lmp, "var", VAR_ATOM)
end

If I set own to false, my memory usage rises by about 2GB. If I set own to true the memory consumption stays stable.
So I'd say its safe to let the GC handle freeing the memory if we're supposed to take ownership.

I have also verified that the data stays allocated even if we close the lammps instance. I've even written this as a test:

lmp = LMP(["-screen", "none"])
command(lmp, """
    atom_modify map yes
    region cell block 0 3 0 3 0 3
    create_box 1 cell
    lattice sc 1
    create_atoms 1 region cell
    mass 1 1

    variable var atom id
""")

var = extract_variable(lmp, "var", VAR_ATOM)
var_copy = copy(var)
LAMMPS.close!(lmp)

@test var == var_copy

In summary, we use copy=true, if the ownership stays with LAMMPS but we want a copy and own=true if the ownership is handed to us, which is usually because LAMMPS has made a copy for us

[vchuravy]

t just calls free internally without doing anything else. Thats the same thing that would happen If we take ownership of the memory if I'm not mistaken.

This is sadly not true. In shared libraries the malloc and free maybe be resolved to different implementations. So you can never rely that the functions are the same in two shared libraries.

if we're supposed to take ownership.

So that's actually a good question. When are we supposed to take ownership?

As far as I understood extract we should never take ownershi, but if the underlying API is allocating memory for us, then we will need to do a copy + lammps_free. Which sucks.

[vchuravy]

Okay according to the docs https://docs.lammps.org/Library_objects.html
atom style requires us to call lammps_free.

[Joroks]

I've thought about registering lammps_free as a finalizer on our array. I think that could work?

[Joroks]

I've done:

finalizer(result) do r
    API.lammps_free(r)
    @async println("finializer has run!")
end

and then the following test script:

lmp = LMP()

command(lmp, """
   atom_modify map yes
   region cell block 0 3 0 3 0 3
   create_box 1 cell
   lattice sc 1
   create_atoms 1 region cell
   mass 1 1

   variable var atom id
""")

var1 = extract_variable(lmp, "var", VAR_ATOM)
var2 = reshape(var1, 1, 27)

var1 = nothing
GC.gc() # finallizer is not run

var2 = nothing
GC.gc() # finallizer is run

[Joroks]

I think the finalizer is the best solution. Should we merge this PR?
Then as soon as 1.11 gets released, we should revisit this. It might become necessary to force 1.11 or higher if we need to register the finalizer on the Memory?

[vchuravy]

So on 1.11

julia> var2 = reshape(var1, 1, 27)
1×27 Matrix{Float64}:
 1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0  9.0  10.0  11.0  12.0  13.0  14.0  …  17.0  18.0  19.0  20.0  21.0  22.0  23.0  24.0  25.0  26.0  27.0

julia> var1 = nothing

julia> GC.gc() # finallizer is run
Finalizer run

julia> var2
1×27 Matrix{Float64}:
 1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0  9.0  10.0  11.0  12.0  13.0  14.0  …  18.0  19.0  20.0  21.0  22.0  23.0  24.0  25.0  26.0  1.265e-321

Note the data flipping.

[vchuravy]

I opened an upstream issue, but I don't expect this behavior to material change before 1.11 gets released

JuliaLang/julia#54985

So if we want to use this finalizer version, we would need to use @static if VERSION >= v1.11-dev and attach the finalizer to result.ref.mem, or unsafe_wrap it as a Memory and then create an Array from that, that is probably clearer.

[Joroks]

this now puts a lot of stress on the GC (in 1.11)

after running my test loop (with 10^8 iterations) i have to manually trigger GC twice for the memory to be freed completely.

this doesn't happen if I attach the finalizer to the array instead of the memory