-
Notifications
You must be signed in to change notification settings - Fork 1
/
exchange.go
218 lines (190 loc) · 7.39 KB
/
exchange.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
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
package engine
// transforms Beff into Forces
import (
"math"
"unsafe"
"github.com/kuchkin/mumax3-gneb/cuda"
"github.com/kuchkin/mumax3-gneb/cuda/cu"
"github.com/kuchkin/mumax3-gneb/data"
"github.com/kuchkin/mumax3-gneb/util"
)
var (
Aex = NewScalarParam("Aex", "J/m", "Exchange stiffness", &lex2)
Dind = NewScalarParam("Dind", "J/m2", "Interfacial Dzyaloshinskii-Moriya strength", &din2)
Dbulk = NewScalarParam("Dbulk", "J/m2", "Bulk Dzyaloshinskii-Moriya strength", &dbulk2)
lex2 exchParam // inter-cell Aex
din2 exchParam // inter-cell Dind
dbulk2 exchParam // inter-cell Dbulk
B_exch = NewVectorField("B_exch", "T", "Exchange field", AddExchangeField)
E_exch = NewScalarValue("E_exch", "J", "Total exchange energy (including the DMI energy)", GetExchangeEnergy)
Edens_exch = NewScalarField("Edens_exch", "J/m3", "Total exchange energy density (including the DMI energy density)", AddExchangeEnergyDensity)
// Average exchange coupling with neighbors. Useful to debug inter-region exchange
ExchCoupling = NewScalarField("ExchCoupling", "arb.", "Average exchange coupling with neighbors", exchangeDecode)
DindCoupling = NewScalarField("DindCoupling", "arb.", "Average DMI coupling with neighbors", dindDecode)
OpenBC = false
//gneb-related parameters
// GNEB2D = false
// GNEB3D = false
JZ = 1.0
// NumOfImages = NewScalarParam("NumOfImages", "dimless", "Number of Images", &NumOfImages)
// NumOfImages exchParam // number of images
)
var AddExchangeEnergyDensity = makeEdensAdder(&B_exch, -0.5) // TODO: normal func
func init() {
registerEnergy(GetExchangeEnergy, AddExchangeEnergyDensity)
DeclFunc("ext_ScaleExchange", ScaleInterExchange, "Re-scales exchange coupling between two regions.")
DeclFunc("ext_InterExchange", InterExchange, "Sets exchange coupling between two regions.")
DeclFunc("ext_ScaleDind", ScaleInterDind, "Re-scales Dind coupling between two regions.")
DeclFunc("ext_InterDind", InterDind, "Sets Dind coupling between two regions.")
DeclVar("OpenBC", &OpenBC, "Use open boundary conditions (default=false)")
// DeclVar("GNEB2D", &GNEB2D, "GNEB?")
// DeclVar("GNEB3D", &GNEB3D, "GNEB?")
DeclVar("JZ", &JZ, "JZ")
lex2.init(Aex)
din2.init(Dind)
dbulk2.init(Dbulk)
}
// Adds the current exchange field to dst
func AddExchangeField(dst *data.Slice) {
inter := !Dind.isZero()
bulk := !Dbulk.isZero()
ms := Msat.MSlice()
defer ms.Recycle()
switch {
case !inter && !bulk:
cuda.AddExchange(dst, M.Buffer(), lex2.Gpu(), ms, regions.Gpu(), M.Mesh(), float32(JZ))
case inter && !bulk:
Refer("mulkers2017")
cuda.AddDMI(dst, M.Buffer(), lex2.Gpu(), din2.Gpu(), ms, regions.Gpu(), M.Mesh(), OpenBC) // dmi+exchange
case bulk && !inter:
cuda.AddDMIBulk(dst, M.Buffer(), lex2.Gpu(), dbulk2.Gpu(), ms, regions.Gpu(), M.Mesh(), OpenBC) // dmi+exchange
// TODO: add ScaleInterDbulk and InterDbulk
case inter && bulk:
util.Fatal("Cannot have interfacial-induced DMI and bulk DMI at the same time")
}
}
// Set dst to the average exchange coupling per cell (average of lex2 with all neighbors).
func exchangeDecode(dst *data.Slice) {
cuda.ExchangeDecode(dst, lex2.Gpu(), regions.Gpu(), M.Mesh())
}
// Set dst to the average dmi coupling per cell (average of din2 with all neighbors).
func dindDecode(dst *data.Slice) {
cuda.ExchangeDecode(dst, din2.Gpu(), regions.Gpu(), M.Mesh())
}
// Returns the current exchange energy in Joules.
func GetExchangeEnergy() float64 {
return -0.5 * cellVolume() * dot(&M_full, &B_exch)
}
// Scales the heisenberg exchange interaction between region1 and 2.
// Scale = 1 means the harmonic mean over the regions of Aex.
func ScaleInterExchange(region1, region2 int, scale float64) {
lex2.setScale(region1, region2, scale)
}
// Sets the exchange interaction between region 1 and 2.
func InterExchange(region1, region2 int, value float64) {
lex2.setInter(region1, region2, value)
}
// Scales the DMI interaction between region 1 and 2.
func ScaleInterDind(region1, region2 int, scale float64) {
din2.setScale(region1, region2, scale)
}
// Sets the DMI interaction between region 1 and 2.
func InterDind(region1, region2 int, value float64) {
din2.setInter(region1, region2, value)
}
// stores interregion exchange stiffness and DMI
// the interregion exchange/DMI by default is the harmonic mean (scale=1, inter=0)
type exchParam struct {
parent *RegionwiseScalar
lut [NREGION * (NREGION + 1) / 2]float32 // harmonic mean of regions (i,j)
scale [NREGION * (NREGION + 1) / 2]float32 // extra scale factor for lut[SymmIdx(i, j)]
inter [NREGION * (NREGION + 1) / 2]float32 // extra term for lut[SymmIdx(i, j)]
gpu cuda.SymmLUT // gpu copy of lut, lazily transferred when needed
gpu_ok, cpu_ok bool // gpu cache up-to date with lut source
}
// to be called after Aex or scaling changed
func (p *exchParam) invalidate() {
p.cpu_ok = false
p.gpu_ok = false
}
func (p *exchParam) init(parent *RegionwiseScalar) {
for i := range p.scale {
p.scale[i] = 1 // default scaling
p.inter[i] = 0 // default additional interexchange term
}
p.parent = parent
}
// Get a GPU mirror of the look-up table.
// Copies to GPU first only if needed.
func (p *exchParam) Gpu() cuda.SymmLUT {
p.update()
if !p.gpu_ok {
p.upload()
}
return p.gpu
}
// sets the interregion exchange/DMI using a specified value (scale = 0)
func (p *exchParam) setInter(region1, region2 int, value float64) {
p.scale[symmidx(region1, region2)] = float32(0.)
p.inter[symmidx(region1, region2)] = float32(value)
p.invalidate()
}
// sets the interregion exchange/DMI by rescaling the harmonic mean (inter = 0)
func (p *exchParam) setScale(region1, region2 int, scale float64) {
p.scale[symmidx(region1, region2)] = float32(scale)
p.inter[symmidx(region1, region2)] = float32(0.)
p.invalidate()
}
func (p *exchParam) update() {
if !p.cpu_ok {
ex := p.parent.cpuLUT()
for i := 0; i < NREGION; i++ {
exi := ex[0][i]
for j := i; j < NREGION; j++ {
exj := ex[0][j]
I := symmidx(i, j)
p.lut[I] = p.scale[I]*exchAverage(exi, exj) + p.inter[I]
}
}
p.gpu_ok = false
p.cpu_ok = true
}
}
func (p *exchParam) upload() {
// alloc if needed
if p.gpu == nil {
p.gpu = cuda.SymmLUT(cuda.MemAlloc(int64(len(p.lut)) * cu.SIZEOF_FLOAT32))
}
lut := p.lut // Copy, to work around Go 1.6 cgo pointer limitations.
cuda.MemCpyHtoD(unsafe.Pointer(p.gpu), unsafe.Pointer(&lut[0]), cu.SIZEOF_FLOAT32*int64(len(p.lut)))
p.gpu_ok = true
}
// Index in symmetric matrix where only one half is stored.
// (!) Code duplicated in exchange.cu
func symmidx(i, j int) int {
if j <= i {
return i*(i+1)/2 + j
} else {
return j*(j+1)/2 + i
}
}
// Returns the intermediate value of two exchange/dmi strengths.
// If both arguments have the same sign, the average mean is returned. If the arguments differ in sign
// (which is possible in the case of DMI), the geometric mean of the geometric and arithmetic mean is
// used. This average is continuous everywhere, monotonic increasing, and bounded by the argument values.
func exchAverage(exi, exj float32) float32 {
if exi*exj >= 0.0 {
if exi == 0 || exj == 0 {
return 0
} else {
return 2 / (1/exi + 1/exj)
}
// return 2.*exi*exj / (exi + exj)
// return 0.5*(exi + exj)
} else {
exi_, exj_ := float64(exi), float64(exj)
sign := math.Copysign(1, exi_+exj_)
magn := math.Sqrt(math.Sqrt(-exi_*exj_) * math.Abs(exi_+exj_) / 2)
return float32(sign * magn)
}
}