custom anisotropy works

cuda/madd.go

@@ -6,6 +6,7 @@ import (
 )
 
 // multiply: dst[i] = a[i] * b[i]
+// a and b must have the same number of components
 func Mul(dst, a, b *data.Slice) {
 	N := dst.Len()
 	nComp := dst.NComp()

engine/customfield.go

@@ -6,14 +6,15 @@ import (
 	"fmt"
 	"github.com/mumax/3/cuda"
 	"github.com/mumax/3/data"
+	"github.com/mumax/3/util"
 )
 
 var (
 	B_custom       = NewVectorField("B_custom", "T", "User-defined field", AddCustomField)
 	Edens_custom   = NewScalarField("Edens_custom", "J/m3", "Energy density of user-defined field.", AddCustomEnergyDensity)
 	E_custom       = NewScalarValue("E_custom", "J", "total energy of user-defined field", GetCustomEnergy)
-	customTerms    []VectorField
-	customEnergies []ScalarField
+	customTerms    []field // vector
+	customEnergies []field // scalar
 )
 
 func init() {
@@ -25,8 +26,14 @@ func init() {
 	DeclFunc("ConstVector", ConstVector, "Constant, uniform vector")
 }
 
-func AddFieldTerm(b outputField) {
-	customTerms = append(customTerms, AsVectorField(b))
+type field interface {
+	Slice() (*data.Slice, bool)
+	NComp() int
+}
+
+// AddFieldTerm adds a function to B_eff
+func AddFieldTerm(b field) {
+	customTerms = append(customTerms, b)
 }
 
 // AddCustomField evaluates the user-defined custom field terms
@@ -57,53 +64,32 @@ func GetCustomEnergy() float64 {
 
 type constValue struct {
 	value []float64
-	*info
 }
 
-func (c *constValue) Mesh() *data.Mesh {
-	return M.Mesh()
-}
+func (c *constValue) NComp() int { return len(c.value) }
 
 func (d *constValue) Slice() (*data.Slice, bool) {
-	buf := cuda.Buffer(d.NComp(), d.Mesh().Size())
+	buf := cuda.Buffer(d.NComp(), Mesh().Size())
 	for c, v := range d.value {
 		cuda.Memset(buf.Comp(c), float32(v))
 	}
 	return buf, true
 }
 
-func (d *constValue) average() []float64 {
-	return d.value
+func Const(v float64) field {
+	return &constValue{[]float64{v}}
 }
 
-func Const(v float64) outputField {
-	// TODO: unit?
-	return &constValue{[]float64{v}, &info{
-		name:  fmt.Sprint(v),
-		unit:  "?", // TODO
-		nComp: 1,
-	}}
-}
-
-func ConstVector(x, y, z float64) outputField {
-	// TODO: unit?
-	return &constValue{[]float64{x, y, z}, &info{
-		name:  fmt.Sprintf("(%v,%v,%v)", x, y, z),
-		unit:  "?", // TODO
-		nComp: 3,
-	}}
+func ConstVector(x, y, z float64) field {
+	return &constValue{[]float64{x, y, z}}
 }
 
 // fieldOp holds the abstract functionality for operations
 // (like add, multiply, ...) on space-dependend quantites
 // (like M, B_sat, ...)
 type fieldOp struct {
-	a, b outputField
-	*info
-}
-
-func (o fieldOp) Mesh() *data.Mesh {
-	return o.a.Mesh()
+	a, b  field
+	nComp int
 }
 
 func (o fieldOp) NComp() int {
@@ -117,15 +103,12 @@ type dotProduct struct {
 // DotProduct creates a new quantity that is the dot product of
 // quantities a and b. E.g.:
 // 	DotProct(&M, &B_ext)
-func Dot(a, b outputField) outputField {
-	return &dotProduct{fieldOp{a, b, &info{
-		nComp: 1,
-		name:  a.Name() + "_dot_" + b.Name(),
-		unit:  a.Unit() + "*" + b.Unit()}}}
+func Dot(a, b field) field {
+	return &dotProduct{fieldOp{a, b, 1}}
 }
 
 func (d *dotProduct) Slice() (*data.Slice, bool) {
-	slice := cuda.Buffer(d.NComp(), d.Mesh().Size())
+	slice := cuda.Buffer(d.NComp(), Mesh().Size())
 	cuda.Zero(slice)
 	A, r := d.a.Slice()
 	if r {
@@ -139,70 +122,116 @@ func (d *dotProduct) Slice() (*data.Slice, bool) {
 	return slice, true
 }
 
-func (d *dotProduct) average() []float64 {
-	return qAverageUniverse(d)
-}
-
-func (d *dotProduct) Average() float64 {
-	return d.average()[0]
-}
-
 type pointwiseMul struct {
 	fieldOp
 }
 
-func Mul(a, b outputField) outputField {
-	return &pointwiseMul{fieldOp{a, b, &info{
-		nComp: a.NComp(),
-		name:  a.Name() + "_mul_" + b.Name(),
-		unit:  a.Unit() + "*" + b.Unit()}}}
+func Mul(a, b field) field {
+	nComp := -1
+	switch {
+	case a.NComp() == b.NComp():
+		nComp = a.NComp() // vector*vector, scalar*scalar
+	case a.NComp() == 1:
+		nComp = b.NComp() // scalar*something
+	case b.NComp() == 1:
+		nComp = a.NComp() // something*scalar
+	default:
+		panic(fmt.Sprintf("Cannot point-wise multiply %v components by %v components", a.NComp(), b.NComp()))
+	}
+
+	return &pointwiseMul{fieldOp{a, b, nComp}}
 }
 
 func (d *pointwiseMul) Slice() (*data.Slice, bool) {
-	slice := cuda.Buffer(d.NComp(), d.Mesh().Size())
-	cuda.Zero(slice)
-	A, r := d.a.Slice()
+	buf := cuda.Buffer(d.NComp(), Mesh().Size())
+	cuda.Zero(buf)
+	a, r := d.a.Slice()
 	if r {
-		defer cuda.Recycle(A)
+		defer cuda.Recycle(a)
 	}
-	B, r := d.b.Slice()
+	b, r := d.b.Slice()
 	if r {
-		defer cuda.Recycle(B)
+		defer cuda.Recycle(b)
 	}
-	cuda.Mul(slice, A, B)
-	return slice, true
+
+	switch {
+	case a.NComp() == b.NComp():
+		mulNN(buf, a, b) // vector*vector, scalar*scalar
+	case a.NComp() == 1:
+		mul1N(buf, a, b)
+	case b.NComp() == 1:
+		mul1N(buf, b, a)
+	default:
+		panic(fmt.Sprintf("Cannot point-wise multiply %v components by %v components", a.NComp(), b.NComp()))
+	}
+
+	return buf, true
+}
+
+// mulNN pointwise multiplies two N-component vectors,
+// yielding an N-component vector stored in dst.
+func mulNN(dst, a, b *data.Slice) {
+	cuda.Mul(dst, a, b)
 }
 
-func (d *pointwiseMul) average() []float64 {
-	return qAverageUniverse(d)
+// mul1N pointwise multiplies a scalar (1-component) with an N-component vector,
+// yielding an N-component vector stored in dst.
+func mul1N(dst, a, b *data.Slice) {
+	util.Assert(a.NComp() == 1)
+	util.Assert(dst.NComp() == b.NComp())
+	for c := 0; c < dst.NComp(); c++ {
+		cuda.Mul(dst.Comp(c), a, b.Comp(c))
+	}
 }
 
 type pointwiseDiv struct {
 	fieldOp
 }
 
-func Div(a, b outputField) outputField {
-	return &pointwiseDiv{fieldOp{a, b, &info{
-		nComp: a.NComp(),
-		name:  a.Name() + "_mul_" + b.Name(),
-		unit:  a.Unit() + "*" + b.Unit()}}}
+func Div(a, b field) field {
+	nComp := -1
+	switch {
+	case a.NComp() == b.NComp():
+		nComp = a.NComp() // vector/vector, scalar/scalar
+	case b.NComp() == 1:
+		nComp = a.NComp() // something/scalar
+	default:
+		panic(fmt.Sprintf("Cannot point-wise divide %v components by %v components", a.NComp(), b.NComp()))
+	}
+	return &pointwiseDiv{fieldOp{a, b, nComp}}
 }
 
 func (d *pointwiseDiv) Slice() (*data.Slice, bool) {
-	slice := cuda.Buffer(d.NComp(), d.Mesh().Size())
-	cuda.Zero(slice)
-	A, r := d.a.Slice()
+	buf := cuda.Buffer(d.NComp(), Mesh().Size())
+	a, r := d.a.Slice()
 	if r {
-		defer cuda.Recycle(A)
+		defer cuda.Recycle(a)
 	}
-	B, r := d.b.Slice()
+	b, r := d.b.Slice()
 	if r {
-		defer cuda.Recycle(B)
+		defer cuda.Recycle(b)
 	}
-	cuda.Div(slice, A, B)
-	return slice, true
+
+	switch {
+	case a.NComp() == b.NComp():
+		divNN(buf, a, b) // vector*vector, scalar*scalar
+	case b.NComp() == 1:
+		divN1(buf, a, b)
+	default:
+		panic(fmt.Sprintf("Cannot point-wise divide %v components by %v components", a.NComp(), b.NComp()))
+	}
+
+	return buf, true
 }
 
-func (d *pointwiseDiv) average() []float64 {
-	return qAverageUniverse(d)
+func divNN(dst, a, b *data.Slice) {
+	cuda.Div(dst, a, b)
+}
+
+func divN1(dst, a, b *data.Slice) {
+	util.Assert(dst.NComp() == a.NComp())
+	util.Assert(b.NComp() == 1)
+	for c := 0; c < dst.NComp(); c++ {
+		cuda.Div(dst.Comp(c), a.Comp(c), b)
+	}
 }

engine/effectivefield.go

@@ -6,8 +6,9 @@ import "github.com/mumax/3/data"
 
 var B_eff = NewVectorField("B_eff", "T", "Effective field", SetEffectiveField)
 
-// Sets dst to the current effective field (T).
-// TODO: extensible slice
+// Sets dst to the current effective field, in Tesla.
+// This is the sum of all effective field terms,
+// like demag, exchange, ...
 func SetEffectiveField(dst *data.Slice) {
 	SetDemagField(dst)    // set to B_demag...
 	AddExchangeField(dst) // ...then add other terms
@@ -16,4 +17,5 @@ func SetEffectiveField(dst *data.Slice) {
 	if !relaxing {
 		B_therm.AddTo(dst)
 	}
+	AddCustomField(dst)
 }

test/custom_anisotropy.mx3

@@ -0,0 +1,41 @@
+/*
+	Test custom field implementation.
+	Like uniaxialanisotropy.mx3, but with custom anisotropy implementation.
+*/
+
+setgridsize(64, 64, 1)
+setcellsize(4e-9, 4e-9, 2e-9)
+
+Aex   = 13e-12
+alpha = 1
+M     = uniform(1, 1, 0)
+
+// Custom anisotropy, easy, in-plane
+Msat  = 1100e3
+K := 0.5e6
+u := ConstVector(1, 0, 0)
+
+prefactor := Const( (2 * K) / (Msat.Average()))
+MyAnis := Mul(prefactor, Mul( Dot(u, m), u))
+AddFieldTerm(MyAnis)
+
+B_ext = vector(0, 0.00, 0)
+relax()
+expect("my", m.average()[1], 0.000, 1e-3)
+
+B_ext = vector(0, 0.01, 0)
+relax()
+expect("my", m.average()[1], 0.011, 1e-3)
+
+B_ext = vector(0, 0.03, 0)
+relax()
+expect("my", m.average()[1], 0.033, 1e-3)
+
+B_ext = vector(0, 0.10, 0)
+relax()
+expect("my", m.average()[1], 0.110, 1e-3)
+
+B_ext = vector(0, 0.30, 0)
+relax()
+expect("my", m.average()[1], 0.331, 1e-3)
+

test/custom_std4.mx3

@@ -0,0 +1,27 @@
+/*
+	Micromagnetic standard problem 4 (a),
+	using a custom-defined zeeman field.
+*/
+
+	setgridsize(128, 32, 1)
+	setcellsize(500e-9/128, 125e-9/32, 3e-9)
+
+	Msat  = 800e3
+	Aex   = 13e-12
+	alpha = 0.02
+	m  = uniform(1, .1, 0)
+
+	minimize()
+	save(m)
+	TOL := 1e-5
+	expectv("m", m.average(), vector(0.9669684171676636,  0.1252732127904892, 0), TOL)
+
+	tableautosave(10e-12)
+	autosave(m, 100e-12)
+	autosnapshot(m, 50e-12)
+
+	myField := div( constvector(-24.6, 4.3, 0), Const(1e3) )
+	AddFieldTerm(myField)
+
+	run(1e-9)
+	expectv("m", m.average(), vector(-0.9846124053001404, 0.12604089081287384, 0.04327124357223511), TOL)