-
Notifications
You must be signed in to change notification settings - Fork 1
/
heun.go
56 lines (46 loc) · 1.07 KB
/
heun.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
package engine
import (
"github.com/kuchkin/mumax3-gneb/cuda"
"github.com/kuchkin/mumax3-gneb/util"
"math"
)
// Adaptive Heun solver.
type Heun struct{}
// Adaptive Heun method, can be used as solver.Step
func (_ *Heun) Step() {
y := M.Buffer()
dy0 := cuda.Buffer(VECTOR, y.Size())
defer cuda.Recycle(dy0)
if FixDt != 0 {
Dt_si = FixDt
}
dt := float32(Dt_si * GammaLL)
util.Assert(dt > 0)
// stage 1
torqueFn(dy0)
cuda.Madd2(y, y, dy0, 1, dt) // y = y + dt * dy
// stage 2
dy := cuda.Buffer(3, y.Size())
defer cuda.Recycle(dy)
Time += Dt_si
torqueFn(dy)
err := cuda.MaxVecDiff(dy0, dy) * float64(dt)
// adjust next time step
if err < MaxErr || Dt_si <= MinDt || FixDt != 0 { // mindt check to avoid infinite loop
// step OK
cuda.Madd3(y, y, dy, dy0, 1, 0.5*dt, -0.5*dt)
M.normalize()
NSteps++
adaptDt(math.Pow(MaxErr/err, 1./2.))
setLastErr(err)
setMaxTorque(dy)
} else {
// undo bad step
util.Assert(FixDt == 0)
Time -= Dt_si
cuda.Madd2(y, y, dy0, 1, -dt)
NUndone++
adaptDt(math.Pow(MaxErr/err, 1./3.))
}
}
func (_ *Heun) Free() {}