Skip to content

Commit a1256a6

Browse files
rpeprpep
committed
Add initial Energy class implementation
1 parent e273782 commit a1256a6

File tree

3 files changed

+191
-0
lines changed

3 files changed

+191
-0
lines changed

‎fidimag/common/vectormath.c

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
#error Do not use this file, it is the result of a failed Cython compilation.

‎fidimag/micro/energ.pyx

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
# distutils: language = c++
2+
3+
cdef extern from "c_energy.h":
4+

‎native/src/c_energy.cpp

Lines changed: 186 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,186 @@
1+
#include "c_energy.h"
2+
3+
double Energy::compute_energy() {
4+
double sum = 0;
5+
for(int i = 0; i < n; i++) {
6+
sum += energy[i];
7+
}
8+
return sum * (dx * dy * dz * std::pow(unit_length, 3));
9+
}
10+
11+
12+
void Energy::setup(int nx, int ny, int nz,
13+
double dx, double dy, double dz,
14+
double unit_length,
15+
double *spin, double *Ms,
16+
double *Ms_inv) {
17+
this->nx = nx;
18+
this->ny = ny;
19+
this->nz = nz;
20+
this->dx = dx;
21+
this->dy = dy;
22+
this->dz = dz;
23+
this->unit_length = unit_length;
24+
this->spin = spin;
25+
this->Ms = Ms;
26+
this->Ms_inv = Ms_inv;
27+
set_up = true;
28+
}
29+
30+
void Exchange::compute_field(double t) {
31+
/* Compute the micromagnetic exchange field and energy using the
32+
* matrix of neighbouring spins and a second order approximation
33+
* for the derivative
34+
*
35+
* Ms_inv :: Array with the (1 / Ms) values for every mesh node.
36+
* The values are zero for points with Ms = 0 (no material)
37+
*
38+
* A :: Exchange constant
39+
*
40+
* dx, dy, dz :: Mesh spacings in the corresponding directions
41+
*
42+
* n :: Number of mesh nodes
43+
*
44+
* ngbs :: The array of neighbouring spins, which has (6 * n)
45+
* entries. Specifically, it contains the indexes of
46+
* the neighbours of every mesh node, in the following order:
47+
* -x, +x, -y, +y, -z, +z
48+
*
49+
* Thus, the array is like:
50+
* | 0-x, 0+x, 0-y, 0+y, 0-z, 0+z, 1-x, 1+x, 1-y, ... |
51+
* i=0 i=1 ...
52+
*
53+
* where 0-y is the index of the neighbour of the 0th spin,
54+
* in the -y direction, for example. The index value for a
55+
* neighbour where Ms = 0, is evaluated as -1. The array
56+
* automatically gives periodic boundaries.
57+
*
58+
* A basic example is a 3 x 3 two dimensional mesh with PBCs
59+
* in the X and Y direction:
60+
*
61+
* +-----------+
62+
* | 6 | 7 | 8 |
63+
* +-----------+
64+
* | 3 | 4 | 5 |
65+
* +-----------+
66+
* | 0 | 1 | 2 |
67+
* +-----------+
68+
*
69+
* so, the first 6 entries (neighbours of the 0th mesh node)
70+
* of the array would be: [ 2 1 6 3 -1 -1 ... ]
71+
* (-1 since there is no material in +-z, and a '2' first,
72+
* since it is the left neighbour which is the PBC in x, etc..)
73+
*
74+
* For the exchange computation, the field is defined as:
75+
* H_ex = (2 * A / (mu0 * Ms)) * nabla^2 (mx, my, mz)
76+
*
77+
* Therefore, for the i-th mesh node (spin), we approximate the
78+
* derivatives as:
79+
* nabla^2 mx = (1 / dx^2) * ( spin[i-x] - 2 * spin[i] + spin[i+x] ) +
80+
* (1 / dy^2) * ( spin[i-y] - 2 * spin[i] + spin[i+y] ) +
81+
* (1 / dz^2) * ( spin[i-z] - 2 * spin[i] + spin[i+z] )
82+
*
83+
* Where i-x is the neighbour in the -x direction. This is similar
84+
* for my and mz.
85+
* We can notice that the sum is the same if we do:
86+
* ( spin[i-x] - spin[i] ) + ( spin[i+x] - spin[i] )
87+
* so we can iterate through the neighbours and perform the sum with the
88+
* corresponding coefficient 1 /dx, 1/dy or 1/dz
89+
*
90+
* The *m array contains the spins as:
91+
* [mx0, my0, mz0, mx1, my1, mz1, mx2, ...]
92+
* so if we want the starting position of the magnetisation for the
93+
* i-th spin, we only have to do (3 * i) for mx, (3 * i + 1) for my
94+
* and (3 * i + 2) for mz
95+
*
96+
*
97+
* IMPORTANT: The ex field usually has the structure:
98+
* 2 * A / (mu0 Ms ) * (Second derivative of M)
99+
* When discretising the derivative, it carries a "2" in the
100+
* denominator which we "cancel" with the "2" in the prefactor,
101+
* hence we do not put it explicitly in the calculations
102+
*
103+
* So, when computing the energy: (-1/2) * mu * Ms * H_ex
104+
* we only put the 0.5 factor and don't worry about the "2"s in the
105+
* field
106+
*
107+
*/
108+
109+
/* Here we iterate through every mesh node */
110+
for (int i = 0; i < n; i++) {
111+
/* Define the coefficients */
112+
double ax = 2 * A[i] / (dx * dx);
113+
double ay = 2 * A[i] / (dy * dy);
114+
double az = 2 * A[i] / (dz * dz);
115+
116+
double fx = 0, fy = 0, fz = 0;
117+
int idnm = 0; // Index for the magnetisation matrix
118+
int idn = 6 * i; // index for the neighbours
119+
120+
/* Set a zero field for sites without magnetic material */
121+
if (Ms_inv[i] == 0.0){
122+
field[3 * i] = 0;
123+
field[3 * i + 1] = 0;
124+
field[3 * i + 2] = 0;
125+
continue;
126+
}
127+
128+
/* Here we iterate through the neighbours */
129+
for (int j = 0; j < 6; j++) {
130+
/* Remember that index=-1 is for sites without material */
131+
if (ngbs[idn + j] >= 0) {
132+
/* Magnetisation of the neighbouring spin since ngbs gives
133+
* the neighbour's index */
134+
idnm = 3 * ngbs[idn + j];
135+
136+
/* Check that the magnetisation of the neighbouring spin
137+
* is larger than zero */
138+
if (Ms_inv[ngbs[idn + j]] > 0){
139+
140+
/* Neighbours in the -x and +x directions
141+
* giving: ( spin[i-x] - spin[i] ) + ( spin[i+x] - spin[i] )
142+
* when ngbs[idn + j] > 0 for j = 0 and j=1
143+
* If, for example, there is no
144+
* neighbour at -x (j=0) in the 0th node (no PBCs),
145+
* the second derivative would only be avaluated as:
146+
* (1 / dx * dx) * ( spin[i+x] - spin[i] )
147+
* which, according to
148+
* [M.J. Donahue and D.G. Porter; Physica B, 343, 177-183 (2004)]
149+
* when performing the integration of the energy, we still
150+
* have error of the order O(dx^2)
151+
* This same applies for the other directions
152+
*/
153+
if (j == 0 || j == 1) {
154+
fx += ax * (spin[idnm] - spin[3 * i]);
155+
fy += ax * (spin[idnm + 1] - spin[3 * i + 1]);
156+
fz += ax * (spin[idnm + 2] - spin[3 * i + 2]);
157+
}
158+
/* Neighbours in the -y and +y directions */
159+
else if (j == 2 || j == 3) {
160+
fx += ay * (spin[idnm] - spin[3 * i]);
161+
fy += ay * (spin[idnm + 1] - spin[3 * i + 1]);
162+
fz += ay * (spin[idnm + 2] - spin[3 * i + 2]);
163+
}
164+
/* Neighbours in the -z and +z directions */
165+
else if (j == 4 || j == 5) {
166+
fx += az * (spin[idnm] - spin[3 * i]);
167+
fy += az * (spin[idnm + 1] - spin[3 * i + 1]);
168+
fz += az * (spin[idnm + 2] - spin[3 * i + 2]);
169+
}
170+
else {
171+
continue;
172+
}
173+
}
174+
}
175+
}
176+
177+
/* Energy as: (-mu0 * Ms / 2) * [ H_ex * m ] */
178+
energy[i] = -0.5 * (fx * spin[3 * i] + fy * spin[3 * i + 1]
179+
+ fz * spin[3 * i + 2]);
180+
181+
/* Update the field H_ex which has the same structure than *m */
182+
field[3 * i] = fx * Ms_inv[i] * MU0_INV;
183+
field[3 * i + 1] = fy * Ms_inv[i] * MU0_INV;
184+
field[3 * i + 2] = fz * Ms_inv[i] * MU0_INV;
185+
}
186+
}

0 commit comments

Comments
 (0)