From 35072f6a317c6837229f579d7ca56fb1c89aebb4 Mon Sep 17 00:00:00 2001 From: ElricleNecro Date: Mon, 5 Aug 2013 13:04:09 +0200 Subject: [PATCH 1/8] Move types declaration to types.go and function attach to ByDist and Search in types_func.go --- octree/octree.go | 42 ------------------------------------------ octree/types.go | 7 +++++++ octree/types_func.go | 36 ++++++++++++++++++++++++++++++++++++ 3 files changed, 43 insertions(+), 42 deletions(-) create mode 100644 octree/types_func.go diff --git a/octree/octree.go b/octree/octree.go index d5fea4e..5cc7a5b 100644 --- a/octree/octree.go +++ b/octree/octree.go @@ -89,48 +89,6 @@ func (e *Node) Create(NbMin int32) { } } -type Search struct { - Radius float64 - Part rg.Particule -} - -type ByDist []Search - -func (e ByDist) Len() int { - return len(e) -} - -func (e ByDist) Swap(i, j int) { - e[i], e[j] = e[j], e[i] -} - -func (e ByDist) Less(i, j int) bool { - return e[i].Radius < e[j].Radius -} - -func Insert(e []Search, to Search) { - var tmp, bak Search - var i int - - // We search where to place the particule : - for i, _ = range e { - if e[i].Radius > to.Radius { - tmp = e[i] - e[i] = to - i++ - break - } - } - - // We move all particule from where we place 'to' to the end of the array - // (we don't keep those who get out from the array) : - for ; i < len(e); i++ { - bak = e[i] - e[i] = tmp - tmp = bak - } -} - func (e *Node) Nearest(part rg.Particule, NbVois int) []Search { var res []Search = make([]Search, NbVois) diff --git a/octree/types.go b/octree/types.go index ca65217..c4ca288 100644 --- a/octree/types.go +++ b/octree/types.go @@ -14,3 +14,10 @@ type Node struct { type WriteStringer interface { WriteString(s string) (ret int, err error) } + +type Search struct { + Radius float64 + Part rg.Particule +} + +type ByDist []Search diff --git a/octree/types_func.go b/octree/types_func.go new file mode 100644 index 0000000..ffe4748 --- /dev/null +++ b/octree/types_func.go @@ -0,0 +1,36 @@ +package octree + +func (e ByDist) Len() int { + return len(e) +} + +func (e ByDist) Swap(i, j int) { + e[i], e[j] = e[j], e[i] +} + +func (e ByDist) Less(i, j int) bool { + return e[i].Radius < e[j].Radius +} + +func Insert(e []Search, to Search) { + var tmp, bak Search + var i int + + // We search where to place the particule : + for i, _ = range e { + if e[i].Radius > to.Radius { + tmp = e[i] + e[i] = to + i++ + break + } + } + + // We move all particule from where we place 'to' to the end of the array + // (we don't keep those who get out from the array) : + for ; i < len(e); i++ { + bak = e[i] + e[i] = tmp + tmp = bak + } +} From de35309cb1939e5cd319ed18485e77e64cb5ff2e Mon Sep 17 00:00:00 2001 From: ElricleNecro Date: Mon, 5 Aug 2013 14:26:38 +0200 Subject: [PATCH 2/8] Add function and new fields to the Node struct for processing the gravitational potential. --- octree/octree.go | 71 ++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 60 insertions(+), 11 deletions(-) diff --git a/octree/octree.go b/octree/octree.go index 5cc7a5b..c763700 100644 --- a/octree/octree.go +++ b/octree/octree.go @@ -1,14 +1,17 @@ //-tags: main package octree -import rg "github.com/ElricleNecro/WisiGo/ReadGadget" -import "fmt" -import m "math" +import ( + "fmt" + rg "github.com/ElricleNecro/WisiGo/ReadGadget" + m "math" +) var ( UseGoRoutine = false //Use of Go Routine for the tree creation? // If you set this to true, don't forget to set runtime.GOMAXPROCS // to the correct value. + Ggrav = 6.67e-11 //Gravitational constant express in International system of units (L^3 M^{-1}T^{-2}). ) // This function is used to create a tree using a list of the particle in the node, the center of the cube and his size. Her essential case of use will be for create the root of the tree. @@ -68,20 +71,24 @@ func (e *Node) Create(NbMin int32) { //We are calculating the particle contain in the Node : NbUtil := t1.setPart() - println((*t1).level, NbUtil, NbUse) - print((*t1).level, " ") - for z := 0; z < 3; z++ { - print("] ", (*t1).Center[z]-(*t1).Size/2., "; ", (*t1).Center[z]+(*t1).Size/2., "] ") - } - println("") + //println((*t1).level, NbUtil, NbUse) + //print((*t1).level, " ") + //for z := 0; z < 3; z++ { + //print("] ", (*t1).Center[z]-(*t1).Size/2., "; ", (*t1).Center[z]+(*t1).Size/2., "] ") + //} + //println("") //If there is particle in the node, we are actualizing the Particle slice, launching calculation of the son and preparing the next brother. if NbUtil != 0 { (*t1).Part = e.Part[NbUse:(NbUse + NbUtil)] + (*t1).Mass = 0.0 + for _, v := range (*t1).Part { + (*t1).Mass += float64(v.Mass) + } - fmt.Println("Going down!") + //fmt.Println("Going down!") (*t1).Create(NbMin) - fmt.Println("Going up!") + //fmt.Println("Going up!") NbUse += NbUtil t1 = &(*t1).Frere @@ -127,6 +134,48 @@ func (e *Node) SearchNeighbor(part rg.Particule, searchy []Search) { } } +func (e *Node) GetBrother() *Node { + return e.Frere +} + +func (e *Node) GetSon() *Node { + return e.Fils +} + +func (e *Node) TotalPotential(opening float64) float64 { + var pot float64 = 0. + for _, v := range e.Part { + pot += e.Potential(v, opening) + } + return pot +} + +func (e *Node) Potential(part rg.Particule, opening float64) float64 { + //We can process the Node and it has sons: + if (float64)(e.Size)/e.Dist(part) >= opening && e.Fils != nil { + var pot float64 = 0.0 + for t1 := e.Fils; t1.Frere != nil; t1 = t1.Frere { + pot += t1.Potential(part, opening) + } + return pot + } + + //The previous condition is wrong, so we have two possibilities: + if e.Fils == nil { + //we are in a leaf: + var pot float64 = 0. + for _, v := range e.Part { + if v.Dist(part) != 0.0 { + pot += Ggrav * (float64)(v.Mass/v.Dist(part)) + } + } + return pot + } else { + //the node is too far or too small to be processed, we do some approximation: + return Ggrav * e.Mass / m.Sqrt((float64)((e.Center[0]-part.Pos[0])*(e.Center[0]-part.Pos[0])+(e.Center[1]-part.Pos[1])*(e.Center[1]-part.Pos[1])+(e.Center[2]-part.Pos[2])*(e.Center[2]-part.Pos[2]))) //e.Dist(part) + } +} + func (e *Node) setPart() (NbUse int32) { NbUse = 0 for i, _ := range e.Part { From 9ed4840a5a486d178455b5167b875003c6ed2174 Mon Sep 17 00:00:00 2001 From: ElricleNecro Date: Mon, 5 Aug 2013 14:27:03 +0200 Subject: [PATCH 3/8] Add test function for the potential calculation. --- octree/octree_test.go | 41 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 40 insertions(+), 1 deletion(-) diff --git a/octree/octree_test.go b/octree/octree_test.go index acad83b..9034652 100644 --- a/octree/octree_test.go +++ b/octree/octree_test.go @@ -1,7 +1,6 @@ package octree import ( - //"fmt" rg "github.com/ElricleNecro/WisiGo/ReadGadget" "math/rand" "os" @@ -245,3 +244,43 @@ func TestSetPart(t *testing.T) { root.Create(1) root.setPart() } + +func TestPotential(t *testing.T) { + rdm := rand.New(rand.NewSource(int64(33))) + nb_part := 10 + part := make([]rg.Particule, nb_part) + var testy rg.Particule + + for i := 0; i < nb_part; i++ { + for j, _ := range part[i].Pos { + part[i].Pos[j] = 2.0*rdm.Float32() - 1.0 + part[i].Id = int64(i + 1) + part[i].Mass = 1.0 + } + } + + center := [...]float32{0., 0., 0.} + root := New(part, center, 2.0) + Ggrav = 1.0 + + testy.Mass = 1.0 + testy.Id = 0 + for i, _ := range testy.Pos { + testy.Pos[i] = 0. + } + + pot := root.Potential(testy, 0.0) + pot2 := bruteforcePotential(part, testy) + + if pot != pot2 { + t.Fatal("Tree Code is different from brute force:", pot, " != ", pot2) + } +} + +func bruteforcePotential(part []rg.Particule, testy rg.Particule) float64 { + var pot float64 = 0.0 + for _, v := range part { + pot += Ggrav * (float64)(v.Mass/v.Dist(testy)) + } + return pot +} From e6aa93b9a27c4c7d3e99e42305f378a1639d44a0 Mon Sep 17 00:00:00 2001 From: ElricleNecro Date: Mon, 5 Aug 2013 14:27:58 +0200 Subject: [PATCH 4/8] Add Mass fields to the Node struct for having the gravitational potential going a little bit faster. --- octree/types.go | 1 + 1 file changed, 1 insertion(+) diff --git a/octree/types.go b/octree/types.go index c4ca288..03778a5 100644 --- a/octree/types.go +++ b/octree/types.go @@ -7,6 +7,7 @@ type Node struct { Parent, Fils, Frere *Node //These is the Main Node to Now : the parent, the first son and the next brother. Center [3]float32 //Center of the cube. Size float32 //Side size of the cube. + Mass float64 //Total Mass contain into the Cube. Part []rg.Particule //List of the particle contain in the node. level int64 //Level of the tree. } From 2ab04ad202053d59182a051d1d3ebd18d78aeea6 Mon Sep 17 00:00:00 2001 From: ElricleNecro Date: Mon, 5 Aug 2013 14:32:55 +0200 Subject: [PATCH 5/8] Add TestTotalPotential function for test purpose --- octree/octree_test.go | 39 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 38 insertions(+), 1 deletion(-) diff --git a/octree/octree_test.go b/octree/octree_test.go index 9034652..10e3048 100644 --- a/octree/octree_test.go +++ b/octree/octree_test.go @@ -277,10 +277,47 @@ func TestPotential(t *testing.T) { } } +func TestTotalPotential(t *testing.T) { + rdm := rand.New(rand.NewSource(int64(33))) + nb_part := 10 + part := make([]rg.Particule, nb_part) + var testy rg.Particule + + for i := 0; i < nb_part; i++ { + for j, _ := range part[i].Pos { + part[i].Pos[j] = 2.0*rdm.Float32() - 1.0 + part[i].Id = int64(i + 1) + part[i].Mass = 1.0 + } + } + + center := [...]float32{0., 0., 0.} + root := New(part, center, 2.0) + Ggrav = 1.0 + + testy.Mass = 1.0 + testy.Id = 0 + for i, _ := range testy.Pos { + testy.Pos[i] = 0. + } + + pot := root.TotalPotential(0.0) + var pot2 float64 = 0.0 + for _, v := range part { + pot2 += bruteforcePotential(part, v) + } + + if pot != pot2 { + t.Fatal("Tree Code is different from brute force:", pot, " != ", pot2) + } +} + func bruteforcePotential(part []rg.Particule, testy rg.Particule) float64 { var pot float64 = 0.0 for _, v := range part { - pot += Ggrav * (float64)(v.Mass/v.Dist(testy)) + if v.Dist(testy) != 0.0 { + pot += Ggrav * (float64)(v.Mass/v.Dist(testy)) + } } return pot } From d07a571b332fb570ccde2cb7e0472d1b0ab4f534 Mon Sep 17 00:00:00 2001 From: ElricleNecro Date: Mon, 5 Aug 2013 14:38:24 +0200 Subject: [PATCH 6/8] Modified the way we walk in the tree in function SearchNeighbor and Potential. Now all Necessary Node are walked in. --- octree/octree.go | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/octree/octree.go b/octree/octree.go index c763700..803ac29 100644 --- a/octree/octree.go +++ b/octree/octree.go @@ -128,7 +128,10 @@ func (e *Node) SearchNeighbor(part rg.Particule, searchy []Search) { return } if e.Fils != nil { - e.Fils.SearchNeighbor(part, searchy[:]) + //e.Fils.SearchNeighbor(part, searchy[:]) + for t1 := e.Fils; t1 != nil; t1 = t1.Frere { + t1.SearchNeighbor(part, searchy[:]) + } } else { e.fill_neighboorhood(part, searchy[:]) } @@ -154,7 +157,7 @@ func (e *Node) Potential(part rg.Particule, opening float64) float64 { //We can process the Node and it has sons: if (float64)(e.Size)/e.Dist(part) >= opening && e.Fils != nil { var pot float64 = 0.0 - for t1 := e.Fils; t1.Frere != nil; t1 = t1.Frere { + for t1 := e.Fils; t1 != nil; t1 = t1.Frere { pot += t1.Potential(part, opening) } return pot From 1d82e00cc8103e6148457b7b80d198185fb460aa Mon Sep 17 00:00:00 2001 From: ElricleNecro Date: Mon, 3 Nov 2014 16:38:26 +0100 Subject: [PATCH 7/8] Forget a 0.5 factor in the potential. --- octree/octree.go | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/octree/octree.go b/octree/octree.go index 803ac29..f0e8f12 100644 --- a/octree/octree.go +++ b/octree/octree.go @@ -150,7 +150,7 @@ func (e *Node) TotalPotential(opening float64) float64 { for _, v := range e.Part { pot += e.Potential(v, opening) } - return pot + return pot * 0.5 } func (e *Node) Potential(part rg.Particule, opening float64) float64 { From 5e5e2bf0001d13f622834dd6f9eb47a624a5e091 Mon Sep 17 00:00:00 2001 From: ElricleNecro Date: Mon, 3 Nov 2014 16:39:44 +0100 Subject: [PATCH 8/8] First attempt to bind my IOGadget library. --- ReadGadget/IOGadget.go | 81 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) create mode 100644 ReadGadget/IOGadget.go diff --git a/ReadGadget/IOGadget.go b/ReadGadget/IOGadget.go new file mode 100644 index 0000000..6b615ff --- /dev/null +++ b/ReadGadget/IOGadget.go @@ -0,0 +1,81 @@ +package ReadGadget + +import ( + "binary" + "fmt" + "log" + "os" +) + +type Gadget struct { + name string + format int + Part []Particule + header Header + log *log.Logger +} + +func New(name string, format int) *Gadget { + res := new(Gadget) + res.name = name + res.format = format + return res +} + +func (c *Gadget) SetLogger(val *log.Logger) { + c.log = val +} + +func (c *Gadget) GetLogger() *log.Logger { + return c.log +} + +type UnsupportedFormat int + +func (c UnsupportedFormat) Error() string { + return fmt.Sprint("Unsupported format: ", c) +} + +type NotImplemented string + +func (c NotImplemented) Error() string { + return fmt.Sprint("You are trying to use a not implemented method: ", c) +} + +func (c *Gadget) Read() error { + switch c.format { + case 1: + return c.read_format1() + case 2: + return c.read_format2() + case 3: + return c.read_format3() + default: + return UnsupportedFormat(c.format) + } + +} + +func (c *Gadget) read_format1() error { + var size, nb_files, pc, pc_new int = 0, 1, 0, 0 + var err error = nil + var file *os.File + + for i := 0; i < nb_files; i++ { + if file, err = os.Open(c.name); err != nil { + return err + } + + pc = pc_new + } + + return NotImplemented("Read()->read_format1()") +} + +func (c *Gadget) read_format2() error { + return NotImplemented("Read()->read_format2()") +} + +func (c *Gadget) read_format3() error { + return NotImplemented("Read()->read_format3()") +}