Fix Gadget2 read and units #8

elehcim · 2021-10-27T07:58:36Z

I added a set of units (the default for Gadget2) which is used when reading/writing files.
I corrected the units of Potential (linked to JuliaAstroSim/PhysicalParticles.jl#18 ).

Also I fixed the usage of StructArray and added some tests.

I'm available this period to contribute on this.

I have a couple of ideas, don't know if you have already thought about that.
For example I'd like to define a getindex method for retrieving only a certain collection from a StructArray of Stars.
I'm not super expert of Julia, but I really want to learn more about it.
Also, there is the need to add Temperature and Energy and star-particle related fields, e.g.: age, chemical composition...
Before going forward, let me know if you have any thought about this. I'll probably open some issues to discuss.

when trying to read past the file to identify new blocks

Do not convert data

src/Gadget.jl

islent · 2021-10-28T14:27:59Z

src/Gadget.jl

@@ -1,3 +1,6 @@
+@unit Gadget2Mass "1e10M⊙" Gadget2MassUnit 1e10*u"Msun" false
+const uGadget2 = [u"kpc",u"kpc/km*s",u"A",u"K",u"cd",Gadget2Mass,u"mol"]


Great design!!!

islent · 2021-10-28T14:33:06Z

src/Gadget.jl

+        vx = read(f, Float32) * 1.0 * uVel
+        vy = read(f, Float32) * 1.0 * uVel
+        vz = read(f, Float32) * 1.0 * uVel


We are assuming the snapshot velocity unit to be km/s (which is not true in some cases), and there need a unit conversion if uVel is not km/s. So I suggest this to be unchanged:

Suggested change

vx = read(f, Float32) * 1.0 * uVel

vy = read(f, Float32) * 1.0 * uVel

vz = read(f, Float32) * 1.0 * uVel

vx = uconvert(uVel, read(f, Float32) * 1.0 * u"km/s")

vy = uconvert(uVel, read(f, Float32) * 1.0 * u"km/s")

vz = uconvert(uVel, read(f, Float32) * 1.0 * u"km/s") #TODO: user defined snapshot unit

Even better, km/s should be assignable (TODO)

islent · 2021-10-28T14:39:41Z

src/Gadget.jl

+                Mass[i] = read(f, Float32) * uMass
            end
        else # set using header
            for i in start:tail
-                Mass[i] = header.mass[type] * 1.0e10 * uMass
+                Mass[i] = header.mass[type] * uMass


The new version only works for uGadget2. I suggest a more general version:

Suggested change

Mass[i] = read(f, Float32) * uMass

end

else # set using header

for i in start:tail

Mass[i] = header.mass[type] * 1.0e10 * uMass

Mass[i] = header.mass[type] * uMass

Mass[i] = uconvert(uMass, read(f, Float32) * 1.0e10 * u"Msun")

end

else # set using header

for i in start:tail

Mass[i] = uconvert(uMass, header.mass[type] * 1.0e10 * u"Msun") #TODO: user defined snapshot unit

islent · 2021-10-28T14:41:12Z

src/Gadget.jl

@@ -210,7 +234,7 @@ function read_Density!(f::Union{IOStream,Stream{format"Gadget2"}}, data::StructA
    temp1 = read(f, Int32)
    Density = data.Density
    for i in 1:NumGas
-        Density[i] = read(f, Float32) * 10e10uDensity
+        Density[i] = read(f, Float32) * 1.0 * uDensity


Suggested change

Density[i] = read(f, Float32) * 1.0 * uDensity

Density[i] = uconvert(uDensity, read(f, Float32) * 1.0 * u"Msun/kpc^3") #TODO: user defined snapshot unit

islent · 2021-10-28T14:41:46Z

src/Gadget.jl

@@ -234,7 +258,7 @@ function read_POT!(f::Union{IOStream,Stream{format"Gadget2"}}, data::StructArray
    temp1 = read(f, Int32)
    Potential = data.Potential
    for i in 1:length(data)
-        Potential[i] = uconvert(uPot, read(f, Float32) * u"km^2/s^2" * data.Mass[i])
+        Potential[i] = read(f, Float32) * 1.0 * uPot


Suggested change

Potential[i] = read(f, Float32) * 1.0 * uPot

Potential[i] = uconvert(uPot, read(f, Float32) * u"km^2/s^2" * data.Mass[i]) #TODO: user defined snapshot unit

islent · 2021-10-28T14:42:31Z

src/Gadget.jl

+        accx = read(f, Float32) * 1.0 * uAcc
+        accy = read(f, Float32) * 1.0 * uAcc
+        accz = read(f, Float32) * 1.0 * uAcc


Suggested change

accx = read(f, Float32) * 1.0 * uAcc

accy = read(f, Float32) * 1.0 * uAcc

accz = read(f, Float32) * 1.0 * uAcc

accx = uconvert(uAcc, read(f, Float32) * 1.0 * u"km^2/kpc/s^2")

accy = uconvert(uAcc, read(f, Float32) * 1.0 * u"km^2/kpc/s^2")

accz = uconvert(uAcc, read(f, Float32) * 1.0 * u"km^2/kpc/s^2") #TODO: user defined snapshot unit

islent · 2021-10-28T14:44:04Z

src/Gadget.jl

@@ -257,7 +281,7 @@ function read_ACCE!(f::Union{IOStream,Stream{format"Gadget2"}}, data::StructArra
    end
 end

-function read_gadget2_particle(f::Union{IOStream,Stream{format"Gadget2"}}, header::HeaderGadget2, units = uAstro;
+function read_gadget2_particle(f::Union{IOStream,Stream{format"Gadget2"}}, header::HeaderGadget2, units = uGadget2;


Changing the default unit may cause errors, so I have to add a TODO mark:

Suggested change

function read_gadget2_particle(f::Union{IOStream,Stream{format"Gadget2"}}, header::HeaderGadget2, units = uGadget2;

function read_gadget2_particle(f::Union{IOStream,Stream{format"Gadget2"}}, header::HeaderGadget2, units = uGadget2;

#TODO check uGadget2 unit

Related to: JuliaAstroSim/PhysicalParticles.jl#22

islent · 2021-10-29T11:17:28Z

The error comes from julia 1.3.0 not supporting keyword abbreviation. Ready to merge.

elehcim · 2021-10-29T11:40:32Z

ok, thanks. I didn't finish to address the units comments. I see the point of them being configurable, I support that idea in general.
In in the manual, in fact, the internal code units are the ones defined in uGadget2. So I'd expect the files to have those units.
Maybe for the moment we can leave them like this.

As a side note: I see these days more and more use of hdf5 as output format for Gadget2 (which is indeed more flexible). This file has units embedded in its header, so configurable units are needed in a while. This is for the (far) future (for now I don't even have example files).

I'm planning to create a Gadget2Particle type inside AstroIO to contain the fields needed by Gadget2 files.
In this way, PhysicalParticles.jl and AstroIO.jl are a bit less interdependent (e.g. a fact that is causing the CI to fail at the moment because the package PhysicalParticles.jl used now in the tests is a bit outdated).

islent · 2021-10-29T11:52:45Z

I'm planning to create a Gadget2Particle type inside AstroIO to contain the fields needed by Gadget2 files.

This is the reason that I'm using Star or AbstractParticle type while implementing gravity solvers. Other users may prefer particle types like Electron, Photon, etc. Let abstract types to guarantee the most possible reusability.

a fact that is causing the CI to fail at the moment because the package PhysicalParticles.jl used now in the tests is a bit outdated

Yes, local tests succeeded. A little wierd to me though.

islent · 2021-10-29T11:55:25Z

For reference, these are critical fields for Gadget-like gravity solver:

    Pos::PVector2D{P}
    Vel::PVector2D{V}
    Acc::PVector2D{A}
    Mass::M
    ID::I
    Collection::Collection

    Ti_endstep::I  # Next integer step on the timeline.
    Ti_begstep::I  # Present integer step on the timeline.
    GravCost::I    # For each two-particle interaction, GravCost += 1
    
    Potential::E   # Particle potential in the force field
    OldAcc::A      # Save the normalization of acceleration of last step. Useful in Tree n-body method.

elehcim added 9 commits October 26, 2021 15:22

Typos

0a85339

Use StructArray properly

de15f14

Convenience function to read only a file header

4b89d1b

Fix EOF error

a1a5537

when trying to read past the file to identify new blocks

Use uGadget2 set of units

5a3c3eb

Do not convert data

Potential in energy unit per unit mass

9737484

Fix mass and density reading

fc323a9

add tests

2ff135b

register Gadget unit and add test

18e23f6

islent self-requested a review October 28, 2021 13:52

islent requested changes Oct 28, 2021

View reviewed changes

islent reviewed Oct 28, 2021

View reviewed changes

islent requested changes Oct 28, 2021

View reviewed changes

islent mentioned this pull request Oct 28, 2021

Define a getindex method for retrieving only a certain collection from a StructArray of Stars? JuliaAstroSim/PhysicalParticles.jl#23

Closed

elehcim added 2 commits October 29, 2021 12:51

no need to filter on GAS particles

c14ed9c

Add test on filter on particle Collection using getingex

ab380b9

Related to: JuliaAstroSim/PhysicalParticles.jl#22

islent merged commit c5e784e into JuliaAstroSim:master Oct 29, 2021

This pull request was closed.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix Gadget2 read and units #8

Fix Gadget2 read and units #8

elehcim commented Oct 27, 2021

islent Oct 28, 2021

islent Oct 28, 2021

islent Oct 28, 2021

islent Oct 28, 2021

islent Oct 28, 2021

islent Oct 28, 2021

islent Oct 28, 2021

islent commented Oct 29, 2021

elehcim commented Oct 29, 2021

islent commented Oct 29, 2021 •

edited

Loading

islent commented Oct 29, 2021

		@@ -1,3 +1,6 @@
		@unit Gadget2Mass "1e10M⊙" Gadget2MassUnit 1e10*u"Msun" false
		const uGadget2 = [u"kpc",u"kpc/km*s",u"A",u"K",u"cd",Gadget2Mass,u"mol"]

	Density[i] = read(f, Float32) * 1.0 * uDensity
	Density[i] = uconvert(uDensity, read(f, Float32) * 1.0 * u"Msun/kpc^3") #TODO: user defined snapshot unit

	Potential[i] = read(f, Float32) * 1.0 * uPot
	Potential[i] = uconvert(uPot, read(f, Float32) * u"km^2/s^2" * data.Mass[i]) #TODO: user defined snapshot unit

	function read_gadget2_particle(f::Union{IOStream,Stream{format"Gadget2"}}, header::HeaderGadget2, units = uGadget2;
	function read_gadget2_particle(f::Union{IOStream,Stream{format"Gadget2"}}, header::HeaderGadget2, units = uGadget2;
	#TODO check uGadget2 unit

Fix Gadget2 read and units #8

Fix Gadget2 read and units #8

Conversation

elehcim commented Oct 27, 2021

islent Oct 28, 2021

Choose a reason for hiding this comment

islent Oct 28, 2021

Choose a reason for hiding this comment

islent Oct 28, 2021

Choose a reason for hiding this comment

islent Oct 28, 2021

Choose a reason for hiding this comment

islent Oct 28, 2021

Choose a reason for hiding this comment

islent Oct 28, 2021

Choose a reason for hiding this comment

islent Oct 28, 2021

Choose a reason for hiding this comment

islent commented Oct 29, 2021

elehcim commented Oct 29, 2021

islent commented Oct 29, 2021 • edited Loading

islent commented Oct 29, 2021

islent commented Oct 29, 2021 •

edited

Loading