Skip to content

Commit 52a0598

Browse files
rpeprpep
authored
Merge pull request #119 from computationalmodelling/refactor_magnetisation
Refactor magnetisation
2 parents 446e98a + 6683963 commit 52a0598

31 files changed

+543
-201
lines changed

‎.gitignore

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,12 +10,14 @@
1010

1111
# ignore automatically generated cython files
1212
fidimag/atomistic/lib/clib.c
13+
fidimag/common/lib/common_clib.c
1314
fidimag/common/dipolar/dipolar.c
1415
fidimag/common/neb/neb_clib.c
1516
fidimag/common/neb_method/*clib.c
1617
fidimag/common/sundials/cvode.c
1718
fidimag/micro/lib/baryakhtar/baryakhtar_clib.c
1819
fidimag/micro/lib/micro_clib.c
20+
fidimag/user/example/example.c
1921

2022
# ignore .cache from pytest
2123
.cache

‎fidimag/atomistic/anisotropy.py

Lines changed: 12 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -50,8 +50,8 @@ def __init__(self, Ku, axis=(1, 0, 0), name='Anisotropy'):
5050
self.axis = axis
5151
self.jac = True
5252

53-
def setup(self, mesh, spin, mu_s):
54-
super(Anisotropy, self).setup(mesh, spin, mu_s)
53+
def setup(self, mesh, spin, mu_s, mu_s_inv):
54+
super(Anisotropy, self).setup(mesh, spin, mu_s, mu_s_inv)
5555

5656
self._Ku = helper.init_scalar(self.Ku, self.mesh)
5757
self._axis = helper.init_vector(self.axis, self.mesh, True)
@@ -65,13 +65,14 @@ def compute_field(self, t=0, spin=None):
6565

6666
clib.compute_anisotropy(m,
6767
self.field,
68+
self.mu_s_inv,
6869
self.energy,
6970
self._Ku,
7071
self._axis,
7172
self.n
7273
)
7374

74-
return self.field * self.mu_s_inv
75+
return self.field
7576

7677

7778
class CubicAnisotropy(Energy):
@@ -84,9 +85,8 @@ def __init__(self, Kc, name='CubicAnisotropy'):
8485
self.name = name
8586
self.jac = True
8687

87-
88-
def setup(self, mesh, spin, mu_s):
89-
super(CubicAnisotropy, self).setup(mesh, spin, mu_s)
88+
def setup(self, mesh, spin, mu_s, mu_s_inv):
89+
super(CubicAnisotropy, self).setup(mesh, spin, mu_s, mu_s_inv)
9090
self._Kc = helper.init_scalar(self.Kc, self.mesh)
9191

9292
def compute_field(self, t=0, spin=None):
@@ -96,9 +96,10 @@ def compute_field(self, t=0, spin=None):
9696
m = self.spin
9797

9898
clib.compute_anisotropy_cubic(m,
99-
self.field,
100-
self.energy,
101-
self._Kc,
102-
self.n)
99+
self.field,
100+
self.mu_s_inv,
101+
self.energy,
102+
self._Kc,
103+
self.n)
103104

104-
return self.field * self.mu_s_inv
105+
return self.field

‎fidimag/atomistic/atomistic_driver.py

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -133,7 +133,11 @@ def set_default_options(self, gamma=1, mu_s=1, alpha=0.1):
133133
# When we create the simulation, mu_s is set to the default value. This
134134
# is overriden when calling the set_mu_s method from the Simulation
135135
# class or when setting mu_s directly (property)
136+
# David Tue 19 Jun 2018: Not a very clear thing to do, we must set a
137+
# WARNING
136138
self._mu_s[:] = mu_s
139+
self._mu_s_inv[:] = 1. / mu_s
140+
137141
self.mu_s_const = mu_s
138142
self.gamma = gamma
139143
self.do_precession = True

‎fidimag/atomistic/demag.py

Lines changed: 7 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,9 @@
11
import fidimag.extensions.dipolar as clib
22
import numpy as np
3+
from .energy import Energy
34

45

5-
class Demag(object):
6+
class Demag(Energy):
67
"""
78
89
Energy class for the demagnetising field (a.k.a. dipolar interactions,
@@ -37,24 +38,13 @@ def __init__(self, name='Demag'):
3738
self.name = name
3839
self.jac = True
3940

40-
def setup(self, mesh, spin, mu_s):
41-
self.mesh = mesh
42-
self.dx = mesh.dx
43-
self.dy = mesh.dy
44-
self.dz = mesh.dz
45-
self.nx = mesh.nx
46-
self.ny = mesh.ny
47-
self.nz = mesh.nz
48-
self.spin = spin
49-
self.n = mesh.n
50-
self.field = np.zeros(3 * self.n, dtype=np.float)
51-
unit_length = mesh.unit_length
52-
self.mu_s_scale = np.zeros(mesh.n, dtype=np.float)
53-
54-
# note that the 1e-7 comes from \frac{\mu_0}{4\pi}
55-
self.scale = 1e-7 / unit_length**3
41+
def setup(self, mesh, spin, mu_s, mu_s_inv):
42+
super(Demag, self).setup(mesh, spin, mu_s, mu_s_inv)
43+
self.scale = 1e-7 / mesh.unit_length**3
5644

5745
# could be wrong, needs carefully tests!!!
46+
# David Tue 19 Jun 2018: This variable is updated in the SIM class in
47+
# case mu_s changes
5848
self.mu_s_scale = mu_s * self.scale
5949

6050
self.demag = clib.FFTDemag(self.dx, self.dy, self.dz,

‎fidimag/atomistic/demag_full.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -18,8 +18,8 @@ def __init__(self, name='DemagFull'):
1818
self.name = name
1919
self.jac = True
2020

21-
def setup(self, mesh, spin, mu_s):
22-
super(DemagFull, self).setup(mesh, spin, mu_s)
21+
def setup(self, mesh, spin, mu_s, mu_s_inv):
22+
super(DemagFull, self).setup(mesh, spin, mu_s, mu_s_inv)
2323

2424
unit_length = mesh.unit_length
2525
self.mu_s_scale = np.zeros(mesh.n, dtype=np.float)

‎fidimag/atomistic/demag_hexagonal.py

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,7 @@ def __init__(self, name='DemagHexagonal'):
5151
self.name = name
5252
self.jac = True
5353

54-
def setup(self, mesh, spin, mu_s):
54+
def setup(self, mesh, spin, mu_s, mu_s_inv):
5555

5656
if mesh.mesh_type != 'hexagonal':
5757
raise Exception('This interaction is only defined'
@@ -89,6 +89,8 @@ def setup(self, mesh, spin, mu_s):
8989
self.scale = 1e-7 / unit_length ** 3
9090

9191
# could be wrong, needs carefully tests!!!
92+
# David Tue 19 Jun 2018: this variable is upated in the sim class
93+
# in case mu_s changes
9294
self.mu_s_scale = mu_s * self.scale
9395

9496
# This seems not to be necessary

‎fidimag/atomistic/dmi.py

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -81,8 +81,8 @@ def __init__(self, D, name='DMI', dmi_type='bulk'):
8181

8282
self.jac = True
8383

84-
def setup(self, mesh, spin, mu_s):
85-
super(DMI, self).setup(mesh, spin, mu_s)
84+
def setup(self, mesh, spin, mu_s, mu_s_inv):
85+
super(DMI, self).setup(mesh, spin, mu_s, mu_s_inv)
8686

8787
if self.mesh_type == 'hexagonal':
8888
self.n_ngbs_dmi = 6
@@ -120,6 +120,7 @@ def compute_field(self, t=0, spin=None):
120120
if self.dmi_type == 'bulk':
121121
clib.compute_dmi_field(m,
122122
self.field,
123+
self.mu_s_inv,
123124
self.energy,
124125
self._D,
125126
self.neighbours,
@@ -131,6 +132,7 @@ def compute_field(self, t=0, spin=None):
131132

132133
clib.compute_dmi_field_interfacial(m,
133134
self.field,
135+
self.mu_s_inv,
134136
self.energy,
135137
self.D,
136138
self.neighbours,
@@ -140,7 +142,7 @@ def compute_field(self, t=0, spin=None):
140142
self.DMI_vector
141143
)
142144

143-
return self.field * self.mu_s_inv
145+
return self.field
144146

145147
def compute_energy_direct(self):
146148
"""

‎fidimag/atomistic/energy.py

Lines changed: 2 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ class Energy(object):
1212
1313
"""
1414

15-
def setup(self, mesh, spin, mu_s):
15+
def setup(self, mesh, spin, mu_s, mu_s_inv):
1616
self.mesh = mesh
1717
self.dx = mesh.dx * mesh.unit_length
1818
self.dy = mesh.dy * mesh.unit_length
@@ -30,16 +30,8 @@ def setup(self, mesh, spin, mu_s):
3030

3131
self.field = np.zeros(3 * self.n, dtype=np.float)
3232
self.energy = np.zeros(mesh.n, dtype=np.float)
33-
self.mu_s_inv = np.zeros(3 * self.n, dtype=np.float)
3433

35-
self.mu_s_inv.shape = (-1, 3)
36-
for i in range(mesh.n):
37-
if self.mu_s[i] == 0.0:
38-
self.mu_s_inv[i, :] = 0
39-
else:
40-
self.mu_s_inv[i, :] = 1.0 / self.mu_s[i]
41-
42-
self.mu_s_inv.shape = (-1,)
34+
self.mu_s_inv = mu_s_inv
4335

4436
self.xperiodic, self.yperiodic = (mesh.periodicity[0],
4537
mesh.periodicity[1])

‎fidimag/atomistic/exchange.py

Lines changed: 11 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -86,8 +86,8 @@ def __init__(self, J, name='Exchange'):
8686
self.name = name
8787
self.jac = False
8888

89-
def setup(self, mesh, spin, mu_s):
90-
super(Exchange, self).setup(mesh, spin, mu_s)
89+
def setup(self, mesh, spin, mu_s, mu_s_inv):
90+
super(Exchange, self).setup(mesh, spin, mu_s, mu_s_inv)
9191

9292
# Uniform exchange ----------------------------------------------------
9393
if isinstance(self.J, (int, float)):
@@ -128,21 +128,23 @@ def compute_field_spatial(self, t=0, spin=None):
128128

129129
clib.compute_exchange_field_spatial(m,
130130
self.field,
131+
self.mu_s_inv,
131132
self.energy,
132133
self._J,
133134
self.neighbours,
134135
self.n,
135136
self.n_ngbs
136137
)
137138

138-
return self.field * self.mu_s_inv
139+
return self.field
139140

140141
def compute_field_uniform(self, t=0, spin=None):
141142

142143
m = spin if spin is not None else self.spin
143144

144145
clib.compute_exchange_field(m,
145146
self.field,
147+
self.mu_s_inv,
146148
self.energy,
147149
self.Jx,
148150
self.Jy,
@@ -152,13 +154,16 @@ def compute_field_uniform(self, t=0, spin=None):
152154
self.n_ngbs
153155
)
154156

155-
return self.field * self.mu_s_inv
157+
return self.field
156158

157159
def compute_field_full(self, t=0, spin=None):
158160

159161
m = spin if spin is not None else self.spin
160162

161-
clib.compute_full_exchange_field(m, self.field, self.energy,
163+
clib.compute_full_exchange_field(m,
164+
self.field,
165+
self.mu_s_inv,
166+
self.energy,
162167
self._J,
163168
self.neighbours,
164169
self.n, self.mesh.n_ngbs,
@@ -167,7 +172,7 @@ def compute_field_full(self, t=0, spin=None):
167172
self.mesh._sum_ngbs_shell
168173
)
169174

170-
return self.field * self.mu_s_inv
175+
return self.field
171176

172177

173178
class UniformExchange(Exchange):

‎fidimag/atomistic/lib/anis.c

Lines changed: 29 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,10 @@
11
#include "clib.h"
22

33

4-
void compute_anis(double *restrict spin, double *restrict field, double *restrict energy,
5-
double *restrict Ku, double *restrict axis, int n) {
4+
void compute_anis(double *restrict spin, double *restrict field,
5+
double *restrict mu_s_inv,
6+
double *restrict energy,
7+
double *restrict Ku, double *restrict axis, int n) {
68

79
/* Remember that the magnetisation order is
810
* mx1, my1, mz1, mx2, my2, mz2, mx3,...
@@ -28,28 +30,43 @@ void compute_anis(double *restrict spin, double *restrict field, double *restric
2830

2931
energy[i] = -Ku[i] * (m_u * m_u);
3032

33+
// Scale field by 1/mu_s
34+
field[3 * i] *= mu_s_inv[i];
35+
field[3 * i + 1] *= mu_s_inv[i];
36+
field[3 * i + 2] *= mu_s_inv[i];
37+
3138
}
3239

3340
}
3441

3542

36-
void compute_anis_cubic(double *restrict spin, double *restrict field, double *restrict energy,
37-
double *restrict Kc, int n) {
43+
void compute_anis_cubic(double *restrict spin, double *restrict field,
44+
double *restrict mu_s_inv,
45+
double *restrict energy,
46+
double *restrict Kc, int n) {
3847

3948
/* Remember that the magnetisation order is
4049
* mx1, my1, mz1, mx2, my2, mz2, mx3,...
4150
* so we get the corresponding components multiplying
4251
* by 3 in every iteration.
4352
*
4453
*/
45-
#pragma omp parallel for
46-
for (int i = 0; i < n; i++) {
47-
int j = 3*i;
48-
field[j] = - 4*Kc[i]*spin[j]*spin[j]*spin[j];
49-
field[j+1] = - 4*Kc[i]*spin[j+1]*spin[j+1]*spin[j+1];
50-
field[j+2] = - 4*Kc[i]*spin[j+2]*spin[j+2]*spin[j+2];
51-
52-
energy[i] = -0.25*(field[j]*spin[j]+field[j+1]*spin[j+1]+field[j+2]*spin[j+2]) ;
54+
#pragma omp parallel for
55+
for (int i = 0; i < n; i++) {
56+
int j = 3 * i;
57+
field[j] = - 4 * Kc[i] * spin[j] * spin[j] * spin[j];
58+
field[j+1] = - 4 * Kc[i] * spin[j+1] * spin[j+1] * spin[j+1];
59+
field[j+2] = - 4 * Kc[i] * spin[j+2] * spin[j+2] * spin[j+2];
60+
61+
energy[i] = -0.25 * (field[j] * spin[j] +
62+
field[j+1] * spin[j+1] +
63+
field[j+2] * spin[j+2]
64+
);
65+
66+
// Scale field by 1/mu_s
67+
field[3 * i] *= mu_s_inv[i];
68+
field[3 * i + 1] *= mu_s_inv[i];
69+
field[3 * i + 2] *= mu_s_inv[i];
5370

5471
}
5572

0 commit comments

Comments
 (0)