-
Notifications
You must be signed in to change notification settings - Fork 0
/
heun.go
52 lines (45 loc) · 1.06 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
package engine
import (
"github.com/mumax/3/cuda"
"github.com/mumax/3/data"
"github.com/mumax/3/util"
"math"
)
// Adaptive Heun method, can be used as solver.Step
func HeunStep(s *solver, y *data.Slice) {
dy0 := cuda.Buffer(VECTOR, y.Size())
defer cuda.Recycle(dy0)
dt := float32(s.Dt_si * *s.dt_mul)
util.Assert(dt > 0)
// stage 1
s.torqueFn(dy0)
s.NEval++
cuda.Madd2(y, y, dy0, 1, dt) // y = y + dt * dy
// stage 2
dy := cuda.Buffer(3, y.Size())
defer cuda.Recycle(dy)
Time += s.Dt_si
s.torqueFn(dy)
s.NEval++
// determine error
err := 0.0
if s.FixDt == 0 { // time step not fixed
err = cuda.MaxVecDiff(dy0, dy) * float64(dt)
}
// adjust next time step
if err < s.MaxErr || s.Dt_si <= s.MinDt { // mindt check to avoid infinite loop
// step OK
cuda.Madd3(y, y, dy, dy0, 1, 0.5*dt, -0.5*dt)
s.postStep()
s.NSteps++
s.adaptDt(math.Pow(s.MaxErr/err, 1./2.))
s.LastErr = err
} else {
// undo bad step
util.Assert(s.FixDt == 0)
Time -= s.Dt_si
cuda.Madd2(y, y, dy0, 1, -dt)
s.NUndone++
s.adaptDt(math.Pow(s.MaxErr/err, 1./3.))
}
}