(draft) Feature: Time Quantity

engine/ext_timequantity.go

@@ -0,0 +1,77 @@
+package engine
+
+import (
+	"github.com/mumax/3/cuda"
+	"github.com/mumax/3/data"
+	"github.com/mumax/3/script"
+)
+
+func init() {
+	DeclROnly("timeQ_s", TimeQ_s, "Current simulation time as a scalar Quantity (s)")
+	DeclROnly("timeQ_v", TimeQ_v, "Current simulation time as a vector Quantity (s)")
+	DeclROnly("timeQ_s2", TimeQ_s2, "Current simulation time as a scalar Quantity (s)")
+	DeclROnly("timeQ_v2", TimeQ_v2, "Current simulation time as a vector Quantity (s)")
+	DeclFunc("ScalarFuncQ", ScalarFuncQ, "Wraps a scalar function of time (e.g. sin(t)) as a Quantity usable in field expressions")
+	DeclFunc("VectorFuncQ", VectorFuncQ, "Wraps three scalar functions of time as a vector Quantity (e.g. VectorFuncQ(sin(t), 0, cos(t)))")
+}
+
+var TimeQ_s = NewScalarField("TimeQ", "s", "Simulation time", func(dst *data.Slice) {
+	cuda.Memset(dst, float32(Time))
+})
+
+var TimeQ_v = NewVectorField("TimeQVec", "s", "Simulation time (vector)", func(dst *data.Slice) {
+	t := float32(Time)
+	for c := 0; c < 3; c++ {
+		cuda.Memset(dst.Comp(c), t)
+	}
+})
+
+// A second potential implementation:
+
+type timeQuantity struct {
+	nComp int
+}
+
+var TimeQ_s2 Quantity = &timeQuantity{1}
+var TimeQ_v2 Quantity = &timeQuantity{3}
+
+func (t *timeQuantity) NComp() int { return t.nComp }
+
+func (t *timeQuantity) EvalTo(dst *data.Slice) {
+	v := float32(Time)
+	for c := 0; c < t.nComp; c++ {
+		cuda.Memset(dst.Comp(c), v)
+	}
+}
+
+// functions of time as Quantities:
+
+type scalarFuncQuantity struct {
+	f script.ScalarFunction
+}
+
+func (q *scalarFuncQuantity) NComp() int { return 1 }
+
+func (q *scalarFuncQuantity) EvalTo(dst *data.Slice) {
+	cuda.Memset(dst.Comp(0), float32(q.f.Float()))
+}
+
+func ScalarFuncQ(f script.ScalarFunction) Quantity {
+	return &scalarFuncQuantity{f}
+}
+
+type vectorFuncQuantity struct {
+	fx, fy, fz script.ScalarFunction
+}
+
+func (q *vectorFuncQuantity) NComp() int { return 3 }
+
+func (q *vectorFuncQuantity) EvalTo(dst *data.Slice) {
+	cuda.Memset(dst.Comp(0), float32(q.fx.Float()))
+	cuda.Memset(dst.Comp(1), float32(q.fy.Float()))
+	cuda.Memset(dst.Comp(2), float32(q.fz.Float()))
+}
+
+func VectorFuncQ(fx, fy, fz script.ScalarFunction) Quantity {
+	return &vectorFuncQuantity{fx, fy, fz}
+}

test/timequantity.mx3

@@ -0,0 +1,121 @@
+/*
+	Compare the TimeQuantity implementation to the existing B_ext.Add() functionality
+*/
+
+Nx := 128
+Ny := 64
+Nz := 2
+c := 5e-9
+setgridsize(Nx, Ny, Nz)
+setcellsize(c, c, c)
+Aex = 10e-12
+Msat = 600e3
+m = uniform(1, 0, 0)
+
+ff := 5e9
+omega := 2*pi*ff
+lambdacus := 137e-9
+kk := (2*pi) / lambdacus
+
+// ===================================================
+// variant 1: timeQ_v
+// ===================================================
+mask := newVectorMask(Nx, Ny, Nz)
+for i := 0; i < Nx; i++ {
+    for j := 0; j < Ny; j++ {
+        for k := 0; k < Nz; k++ {
+            r := index2coord(i, j, k)
+            x := r.X()
+            mask.setVector(i, j, k, vector(0, 0, x))
+        }
+    }
+}
+xmask := CustomQuantity(mask)
+
+innertterm := Mul(Const(-2*pi*ff), timeQ_v)
+innterxterm := Mul(Const(kk), xmask)
+innertotal := Add(innterxterm , innertimeterm)
+sinwave := Qsin(innertotal)
+v1term := Mul(Constvector(0.00, 0.00, 0.01), sinwave)
+AddFieldTerm(v1term)
+
+run(1e-12)
+ref := B_custom.Comp(2).Average()
+print("variant 1 (timeQ_v) B_custom.z average:", ref)
+
+RemoveCustomFields()
+t = 0
+m = uniform(1, 0, 0)
+
+// ===================================================
+// variant 2: VectorFuncQ
+// ===================================================
+mask_sin := newVectorMask(Nx, Ny, Nz)
+for i := 0; i < Nx; i++ {
+    for j := 0; j < Ny; j++ {
+        for k := 0; k < Nz; k++ {
+            r := index2coord(i, j, k)
+            x := r.X()
+            mask_sin.setVector(i, j, k, vector(0, 0, sin(kk*x)))
+        }
+    }
+}
+xmask_sin := CustomQuantity(mask_sin)
+
+mask_cos := newVectorMask(Nx, Ny, Nz)
+for i := 0; i < Nx; i++ {
+    for j := 0; j < Ny; j++ {
+        for k := 0; k < Nz; k++ {
+            r := index2coord(i, j, k)
+            x := r.X()
+            mask_cos.setVector(i, j, k, vector(0, 0, cos(kk*x)))
+        }
+    }
+}
+xmask_cos := CustomQuantity(mask_cos)
+
+sinTerm := Mul(xmask_sin, VectorFuncQ(0, 0, cos(omega*t)))
+cosTerm := Mul(xmask_cos, VectorFuncQ(0, 0, -sin(omega*t)))
+AddFieldTerm(Mul(Const(0.01), Add(sinTerm, cosTerm)))
+
+run(1e-12)
+v2 := B_custom.Comp(2).Average()
+print("variant 2 (VectorFuncQ) B_custom.z average:", v2)
+Expect("VectorFuncQ vs timeQ_v z", v2, ref, 1e-7)
+
+RemoveCustomFields()
+t = 0
+m = uniform(1, 0, 0)
+
+// ===================================================
+// variant 3: classical B_ext.add
+// ===================================================
+mask2 := newVectorMask(Nx, Ny, Nz)
+for i := 0; i < Nx; i++ {
+    for j := 0; j < Ny; j++ {
+        for k := 0; k < Nz; k++ {
+            r := index2coord(i, j, k)
+            x := r.X()
+            mask2.setVector(i, j, k, vector(0, 0, sin(kk*x)))
+        }
+    }
+}
+B_ext.add(mask2, 0.01*cos(omega*t))
+
+mask3 := newVectorMask(Nx, Ny, Nz)
+for i := 0; i < Nx; i++ {
+    for j := 0; j < Ny; j++ {
+        for k := 0; k < Nz; k++ {
+            r := index2coord(i, j, k)
+            x := r.X()
+            mask3.setVector(i, j, k, vector(0, 0, cos(kk*x)))
+        }
+    }
+}
+B_ext.add(mask3, -0.01*sin(omega*t))
+
+run(1e-12)
+v3 := B_ext.Comp(2).Average()
+print("variant 3 (B_ext.add) B_ext.z average:", v3)
+Expect("B_ext.add vs timeQ_v z", v3, ref, 1e-7)
+