Merge pull request #135 from mumax/feature/bubbleshift

Feature/bubbleshift

engine/ext_bubblepos.go

@@ -12,41 +12,40 @@ var (
 
 func bubblePos() []float64 {
 	m := M.Buffer()
+	n := Mesh().Size()
+	c := Mesh().CellSize()
 	mz := m.Comp(Z).HostCopy().Scalars()[0]
 
-	posx, posy := 0, 0
+	posx, posy := 0., 0.
+	sign := magsign(M.GetCell(0, n[Y]/2, n[Z]/2)[Z]) //TODO make more robust with temperature?
 
 	{
-		max := float32(-1e32)
+		var magsum float32
+		var weightedsum float32
+
 		for iy := range mz {
-			var sum float32
-			for ix := range mz[iy] {
-				sum += mz[iy][ix]
-			}
-			if sum > max {
-				posy = iy
-				max = sum
+			for ix := range mz[0] {
+				magsum += ((mz[iy][ix]*float32(-1*sign) + 1.) / 2.)
+				weightedsum += ((mz[iy][ix]*float32(-1*sign) + 1.) / 2.) * float32(iy)
 			}
 		}
+		posy = float64(weightedsum / magsum)
 	}
 
 	{
-		max := float32(-1e32)
+		var magsum float32
+		var weightedsum float32
+
 		for ix := range mz[0] {
-			var sum float32
 			for iy := range mz {
-				sum += mz[iy][ix]
-			}
-			if sum > max {
-				posx = ix
-				max = sum
+				magsum += ((mz[iy][ix]*float32(-1*sign) + 1.) / 2.)
+				weightedsum += ((mz[iy][ix]*float32(-1*sign) + 1.) / 2.) * float32(ix)
 			}
 		}
+		posx = float64(weightedsum / magsum)
 	}
 
-	c := Mesh().CellSize()
-	n := Mesh().Size()
-	return []float64{float64(posx-n[X]/2)*c[X] + GetShiftPos(), float64(posy-n[Y]/2) * c[Y], 0}
+	return []float64{(posx-float64(n[X]/2))*c[X] + GetShiftPos(), (posy-float64(n[Y]/2))*c[Y] + GetShiftYPos(), 0}
 }
 
 var (

engine/ext_centerbubble.go

@@ -0,0 +1,71 @@
+package engine
+
+import (
+	"github.com/mumax/3/data"
+	"math"
+)
+
+func init() {
+	DeclFunc("ext_centerBubble", CenterBubble, "centerBubble shifts m after each step to keep the bubble position close to the center of the window")
+}
+
+func centerBubble() {
+	M := &M
+	n := Mesh().Size()
+
+	m := M.Buffer()
+	mz := m.Comp(Z).HostCopy().Scalars()[0]
+
+	posx, posy := 0., 0.
+	sign := magsign(M.GetCell(0, n[Y]/2, n[Z]/2)[Z]) //TODO make more robust with temperature?
+
+	{
+		var magsum float32
+		var weightedsum float32
+
+		for iy := range mz {
+			for ix := range mz[0] {
+				magsum += ((mz[iy][ix]*float32(-1*sign) + 1.) / 2.)
+				weightedsum += ((mz[iy][ix]*float32(-1*sign) + 1.) / 2.) * float32(iy)
+			}
+		}
+		posy = float64(weightedsum / magsum)
+	}
+
+	{
+		var magsum float32
+		var weightedsum float32
+
+		for ix := range mz[0] {
+			for iy := range mz {
+				magsum += ((mz[iy][ix]*float32(-1*sign) + 1.) / 2.)
+				weightedsum += ((mz[iy][ix]*float32(-1*sign) + 1.) / 2.) * float32(ix)
+			}
+		}
+		posx = float64(weightedsum / magsum)
+	}
+
+	zero := data.Vector{0, 0, 0}
+	if ShiftMagL == zero || ShiftMagR == zero || ShiftMagD == zero || ShiftMagU == zero {
+		ShiftMagL[Z] = float64(sign)
+		ShiftMagR[Z] = float64(sign)
+		ShiftMagD[Z] = float64(sign)
+		ShiftMagU[Z] = float64(sign)
+	}
+	dx := int(math.Floor(float64(n[X]/2) - posx))
+	dy := int(math.Floor(float64(n[Y]/2) - posy))
+
+	//put bubble to center
+	if dx != 0 {
+		Shift(dx)
+	}
+	if dy != 0 {
+		YShift(dy)
+	}
+
+}
+
+// This post-step function centers the simulation window on a bubble
+func CenterBubble() {
+	PostStep(func() { centerBubble() })
+}

engine/geom.go

@@ -234,6 +234,36 @@ func (g *geom) shift(dx int) {
 
 }
 
+func (g *geom) shiftY(dy int) {
+	// empty mask, nothing to do
+	if g == nil || g.buffer.IsNil() {
+		return
+	}
+
+	// allocated mask: shift
+	s := g.buffer
+	s2 := cuda.Buffer(1, g.Mesh().Size())
+	defer cuda.Recycle(s2)
+	newv := float32(1) // initially fill edges with 1's
+	cuda.ShiftY(s2, s, dy, newv, newv)
+	data.Copy(s, s2)
+
+	n := Mesh().Size()
+	y1, y2 := shiftDirtyRange(dy)
+
+	for iz := 0; iz < n[Z]; iz++ {
+		for ix := 0; ix < n[X]; ix++ {
+			for iy := y1; iy < y2; iy++ {
+				r := Index2Coord(ix, iy, iz) // includes shift
+				if !g.shape(r[X], r[Y], r[Z]) {
+					cuda.SetCell(g.buffer, 0, ix, iy, iz, 0) // a bit slowish, but hardly reached
+				}
+			}
+		}
+	}
+
+}
+
 // x range that needs to be refreshed after shift over dx
 func shiftDirtyRange(dx int) (x1, x2 int) {
 	nx := Mesh().Size()[X]

engine/regions.go

@@ -238,6 +238,31 @@ func (b *Regions) shift(dx int) {
 	}
 }
 
+func (b *Regions) shiftY(dy int) {
+	// TODO: return if no regions defined
+	r1 := b.Gpu()
+	r2 := cuda.NewBytes(b.Mesh().NCell()) // TODO: somehow recycle
+	defer r2.Free()
+	newreg := byte(0) // new region at edge
+	cuda.ShiftBytesY(r2, r1, b.Mesh(), dy, newreg)
+	r1.Copy(r2)
+
+	n := Mesh().Size()
+	y1, y2 := shiftDirtyRange(dy)
+
+	for iz := 0; iz < n[Z]; iz++ {
+		for ix := 0; ix < n[X]; ix++ {
+			for iy := y1; iy < y2; iy++ {
+				r := Index2Coord(ix, iy, iz) // includes shift
+				reg := b.get(r)
+				if reg != 0 {
+					b.SetCell(ix, iy, iz, reg) // a bit slowish, but hardly reached
+				}
+			}
+		}
+	}
+}
+
 func (r *Regions) Mesh() *data.Mesh { return Mesh() }
 
 func prod(s [3]int) int {

engine/shift.go

@@ -6,23 +6,26 @@ import (
 )
 
 var (
-	TotalShift                      float64                        // accumulated window shift (X) in meter
-	ShiftMagL, ShiftMagR            data.Vector                    // when shifting m, put these value at the left/right edge.
-	ShiftM, ShiftGeom, ShiftRegions bool        = true, true, true // should shift act on magnetization, geometry, regions?
+	TotalShift, TotalYShift                    float64                        // accumulated window shift (X and Y) in meter
+	ShiftMagL, ShiftMagR, ShiftMagU, ShiftMagD data.Vector                    // when shifting m, put these value at the left/right edge.
+	ShiftM, ShiftGeom, ShiftRegions            bool        = true, true, true // should shift act on magnetization, geometry, regions?
 )
 
 func init() {
 	DeclFunc("Shift", Shift, "Shifts the simulation by +1/-1 cells along X")
 	DeclVar("ShiftMagL", &ShiftMagL, "Upon shift, insert this magnetization from the left")
 	DeclVar("ShiftMagR", &ShiftMagR, "Upon shift, insert this magnetization from the right")
+	DeclVar("ShiftMagU", &ShiftMagU, "Upon shift, insert this magnetization from the top")
+	DeclVar("ShiftMagD", &ShiftMagD, "Upon shift, insert this magnetization from the bottom")
 	DeclVar("ShiftM", &ShiftM, "Whether Shift() acts on magnetization")
 	DeclVar("ShiftGeom", &ShiftGeom, "Whether Shift() acts on geometry")
 	DeclVar("ShiftRegions", &ShiftRegions, "Whether Shift() acts on regions")
 	DeclVar("TotalShift", &TotalShift, "Amount by which the simulation has been shifted (m).")
 }
 
 // position of the window lab frame
-func GetShiftPos() float64 { return -TotalShift }
+func GetShiftPos() float64  { return -TotalShift }
+func GetShiftYPos() float64 { return -TotalYShift }
 
 // shift the simulation window over dx cells in X direction
 func Shift(dx int) {
@@ -48,3 +51,28 @@ func shiftMag(m *data.Slice, dx int) {
 		data.Copy(comp, m2) // str0 ?
 	}
 }
+
+// shift the simulation window over dy cells in Y direction
+func YShift(dy int) {
+	TotalYShift += float64(dy) * Mesh().CellSize()[Y] // needed to re-init geom, regions
+	if ShiftM {
+		shiftMagY(M.Buffer(), dy)
+	}
+	if ShiftRegions {
+		regions.shiftY(dy)
+	}
+	if ShiftGeom {
+		geometry.shiftY(dy)
+	}
+	M.normalize()
+}
+
+func shiftMagY(m *data.Slice, dy int) {
+	m2 := cuda.Buffer(1, m.Size())
+	defer cuda.Recycle(m2)
+	for c := 0; c < m.NComp(); c++ {
+		comp := m.Comp(c)
+		cuda.ShiftY(m2, comp, dy, float32(ShiftMagU[c]), float32(ShiftMagD[c]))
+		data.Copy(comp, m2) // str0 ?
+	}
+}

engine/util.go

@@ -132,7 +132,7 @@ func Index2Coord(ix, iy, iz int) data.Vector {
 	n := m.Size()
 	c := m.CellSize()
 	x := c[X]*(float64(ix)-0.5*float64(n[X]-1)) - TotalShift
-	y := c[Y] * (float64(iy) - 0.5*float64(n[Y]-1))
+	y := c[Y]*(float64(iy)-0.5*float64(n[Y]-1)) - TotalYShift
 	z := c[Z] * (float64(iz) - 0.5*float64(n[Z]-1))
 	return data.Vector{x, y, z}
 }

test/bubbleshiftpos.mx3

@@ -0,0 +1,40 @@
+SetMesh(128, 128, 1, 1e-9,1e-9,0.4e-9, 1, 1, 0)
+
+Msat  =580e3
+Aex   = 15e-12
+enabledemag=false
+alpha = 0.1
+Ku1=0.59e6
+anisU=vector(0,0,1)
+Dind=0.0034089785
+
+shiftregions=true
+
+maxregion:=255
+seed:=17
+ext_makegrains(10e-9, maxregion, seed)
+
+for i:=0; i<maxRegion; i++{
+        for j:=i+1; j<maxRegion; j++{
+                ext_ScaleExchange(i, j, 0.5)
+        }
+}
+
+
+m =neelskyrmion(1, -1)
+run(5e-11)
+t=0.
+minimize()
+
+setsolver(5)
+maxerr=0.001
+
+ext_centerbubble()
+
+Xi=0.2
+j=vector(5e12*cos(2*pi*1e9*t),5e12*sin(2*pi*1e9*t),0)
+
+run(2e-10)
+TOL:=1e-11
+expectv("position", ext_bubblepos.average(), vector(-7.202214431762697e-08,-8.216819381713864e-08,0), TOL)
+