Skip to content

Commit 00887c9

Browse files
davidcortesortunodavidcortesortuno
committed
Documented C code for the calculation of the right hand side of the LLG equation. Since this code is both used by atomistic and micromagnetic evolvers, we should move it to the common folder in the future
1 parent 3d6b29d commit 00887c9

File tree

1 file changed

+162
-65
lines changed
  • fidimag/atomistic/lib

1 file changed

+162
-65
lines changed

‎fidimag/atomistic/lib/llg.c‎

Lines changed: 162 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -1,70 +1,167 @@
11
#include "clib.h"
22

3-
void llg_rhs(double *dm_dt, double *m, double *h, double *alpha, int *pins,
4-
double gamma, int n, int do_precession, double default_c) {
3+
/* The right hand side of the LLG equation for the CVOde solver. This can be
4+
* used both for the micromagnetic and atomistic codes since m or S are unitless
5+
* and the prefactors keep the same structure.
6+
*
7+
* The LLG equation has the structure:
8+
* ( * is for dot or scalar product)
9+
*
10+
* dm -gamma
11+
* ---- = -------- ( m X H_eff + a * m X ( m x H_eff ) )
12+
* dt 2
13+
* ( 1 + a )
14+
*
15+
* where a is the Gilbert damping constant, gamma is the gyromagnetic
16+
* ratio ( gamma = 1.76e11 for a free electron; for the micromagnetic
17+
* case, we use gamma_0 = mu_0 * gamma = 2.21e5, for a free electron),
18+
* m is the magnetisation vector and H_eff is the effective field.
19+
*
20+
* In our calculation, we usually compute:
21+
22+
* m x (m x H_eff) = ( m * H_eff) m - (m * m) H_eff
23+
24+
* Then, if we define the perpendicular component of the effective
25+
* field as:
26+
27+
* H_perp = H_eff - ( m * H_eff) m
28+
29+
* we can write
30+
31+
* dm -gamma
32+
* ---- = -------- ( m X H_perp - a * H_perp )
33+
* dt 2
34+
* ( 1 + a )
35+
36+
* since m X m = 0 and, theoretically, m * m = 1 (at zero temperature m has
37+
* fixed length). However, for the second term to the right hand side,
38+
* proportional to a, it is better to keep the m X ( m x H_eff ) term with (m *
39+
* m) for stability purposes, thus we use
40+
*
41+
42+
* H_perp = (m * m) H_eff - ( m * H_eff) m
43+
44+
* for both terms.
45+
*
46+
* Additionally, to preserve the magnetisation length, we need to correct the
47+
* dm / dt term for every time step, adding the following term to the right
48+
* hand size expression of the LLG:
49+
50+
* dm dm 2
51+
* ---- ---> ---- + c * ( 1 - m ) m
52+
* dt dt
53+
54+
* with
55+
_________
56+
* / 2
57+
* c = 6 * / ( dm )
58+
* / ( ---- )
59+
* \/ ( dt )
60+
61+
* The correction must be introduced since numerically, m can change length
62+
* during a step of the integration (which would not occur if the integration
63+
* step is infinitely small), deviating from the real answer. If we just
64+
* rescaled m, we would have to recompute the effective field (m changes) and
65+
* also the solution would present jumps due to the rescaling. With the term
66+
* specified above, m stays close to 1, giving a more continuous solution, and
67+
* also the term stays as zero if m is exactly 1 (notice that if m increases,
68+
* the correction is negative and it decreases m length; similarly if m
69+
* decreases). The prefactor c must be specified because the correction term
70+
* must be sufficiently strong to affect the solution. Accordingly, we can
71+
* think of dm/dt as a kind of velocity that is proportional to the change of
72+
* rate of m, hence using its magnitude, the correction is stronger for large
73+
* deviations and weak for small deviations. The factor 6 is added ad-hoc,
74+
* which seems to work well when computing the solutions, but its specification
75+
* stills requires a more strict proof. It is worth mentioning that the norm of
76+
* dm/dt changes the time scaling by a factor proportional to 1/t, therefore in
77+
* the future we could try to estimate its
78+
* influence with more mathematical/numerical rigour and analyse an optimal
79+
* value for the prefactor (6).
80+
*/
81+
82+
void llg_rhs(double * dm_dt, double * m, double * h, double * alpha, int * pins,
83+
double gamma, int n, int do_precession, double default_c) {
584

6-
int i, j, k;
85+
int i, j, k;
786

887
double coeff, mm, mh, c;
9-
double hpi,hpj,hpk;
10-
11-
#pragma omp parallel for private(i,j,k,coeff,mm, mh, c, hpi,hpj,hpk)
12-
for (int id = 0; id < n; id++) {
13-
i = 3*id;
14-
j = i+1;
15-
k = i+2;
16-
17-
if (pins[id]>0){
18-
dm_dt[i] = 0;
19-
dm_dt[j] = 0;
20-
dm_dt[k] = 0;
21-
continue;
22-
}
23-
24-
coeff = -gamma/(1.0+alpha[id]*alpha[id]);
25-
26-
mm = m[i]*m[i] + m[j]*m[j] + m[k]*m[k];
27-
mh = m[i]*h[i] + m[j]*h[j] + m[k]*h[k];
28-
29-
//suppose m is normalised, i.e., mm=1; hp=mm.h-mh.m=-mx(mxh)
30-
hpi = mm*h[i] - mh*m[i];
31-
hpj = mm*h[j] - mh*m[j];
32-
hpk = mm*h[k] - mh*m[k];
33-
//IMPORTANT: never ignore mm!!!
34-
//what we found is that if we igonre mm, i.e. using
35-
// hpi = h[i] - mh*m[i]
36-
// hpj = h[j] - mh*m[j];
37-
// hpk = h[k] - mh*m[k];
38-
//then the standard problem 4 failed to converge ?!!
39-
double mth0=0, mth1=0, mth2=0;
88+
double hpi, hpj, hpk;
89+
90+
#pragma omp parallel for private(i,j,k,coeff,mm, mh, c, hpi,hpj,hpk)
91+
for (int id = 0; id < n; id++) {
92+
// Indexes for the 3 components of the spin (magnetic moment)
93+
// at the i-th lattice (mesh) site --> x, y, z
94+
i = 3 * id;
95+
j = i + 1;
96+
k = i + 2;
97+
98+
// Pinned spins do not follow the dynamical equation
99+
if (pins[id] > 0) {
100+
dm_dt[i] = 0;
101+
dm_dt[j] = 0;
102+
dm_dt[k] = 0;
103+
continue;
104+
}
40105

106+
coeff = -gamma / (1.0 + alpha[id] * alpha[id]);
107+
108+
// Dot products
109+
mm = m[i] * m[i] + m[j] * m[j] + m[k] * m[k];
110+
mh = m[i] * h[i] + m[j] * h[j] + m[k] * h[k];
111+
112+
// Usually, m is normalised, i.e., mm=1;
113+
// so hp = mm.h - mh.m = -m x (m x h)
114+
// We set here the perpendicular componenet of the field
115+
// but using the (m * m) product
116+
hpi = mm * h[i] - mh * m[i];
117+
hpj = mm * h[j] - mh * m[j];
118+
hpk = mm * h[k] - mh * m[k];
119+
120+
// IMPORTANT: do not ignore mm !!
121+
// What we've found is that if we igonre mm, i.e. using
122+
// hpi = h[i] - mh * m[i];
123+
// hpj = h[j] - mh * m[j];
124+
// hpk = h[k] - mh * m[k];
125+
// the micromagnetic standard problem 4 failed to converge (?)
126+
//
127+
// NOTE (Fri 08 Jul 2016 13:58): In fact, the problem converges but with 2 less
128+
// decimals of accuracy, compared with the OOMMF calculation
129+
double mth0 = 0, mth1 = 0, mth2 = 0;
130+
131+
// The first term: m x H_eff = m x H_perp
41132
if (do_precession){
42-
mth0 = cross_x(m[i],m[j],m[k],hpi,hpj,hpk);
43-
mth1 = cross_y(m[i],m[j],m[k],hpi,hpj,hpk);
44-
mth2 = cross_z(m[i],m[j],m[k],hpi,hpj,hpk);
133+
mth0 = cross_x(m[i], m[j], m[k], hpi, hpj, hpk);
134+
mth1 = cross_y(m[i], m[j], m[k], hpi, hpj, hpk);
135+
mth2 = cross_z(m[i], m[j], m[k], hpi, hpj, hpk);
45136
}
46137

47-
dm_dt[i] = coeff*(mth0 - hpi * alpha[id]);
48-
dm_dt[j] = coeff*(mth1 - hpj * alpha[id]);
49-
dm_dt[k] = coeff*(mth2 - hpk * alpha[id]);
50-
51-
// in future, we will try the new method to integrate the LLG equation,
52-
// A mixed mid-point Runge-Kutta like scheme for the integration of Landau-Lifshitz equation
53-
// Journal of Applied Physics 115, 17D101 (2014)
54-
// if possible, we can combine it with adaptive step size, don't know how to do but it's worth a try.
55-
56-
if (default_c<0){
57-
c = 6*sqrt(dm_dt[i]*dm_dt[i]+dm_dt[j]*dm_dt[j]+dm_dt[k]*dm_dt[k]);
58-
}else{
59-
c = default_c;
138+
// The RHS term of the LLG equation
139+
dm_dt[i] = coeff * (mth0 - hpi * alpha[id]);
140+
dm_dt[j] = coeff * (mth1 - hpj * alpha[id]);
141+
dm_dt[k] = coeff * (mth2 - hpk * alpha[id]);
142+
143+
// In future, we will try the new method to integrate the LLG equation,
144+
// A mixed mid-point Runge-Kutta like scheme for the integration of
145+
// Landau-Lifshitz equation Journal of Applied Physics 115, 17D101
146+
// (2014) if possible, we can combine it with adaptive step size, don't
147+
// know how to do but it's worth a try.
148+
149+
if (default_c < 0){
150+
c = 6 * sqrt(dm_dt[i] * dm_dt[i] +
151+
dm_dt[j] * dm_dt[j] +
152+
dm_dt[k] * dm_dt[k]
153+
);
154+
} else {
155+
c = default_c;
60156
}
61157
//printf("%0.15g %0.15g\n", c, default_c);
62158

63-
dm_dt[i] += c*(1-mm)*m[i];
64-
dm_dt[j] += c*(1-mm)*m[j];
65-
dm_dt[k] += c*(1-mm)*m[k];
159+
// Correct the RHS term to keep m normalised
160+
dm_dt[i] += c * (1 - mm) * m[i];
161+
dm_dt[j] += c * (1 - mm) * m[j];
162+
dm_dt[k] += c * (1 - mm) * m[k];
66163

67-
}
164+
}
68165

69166
}
70167

@@ -106,9 +203,9 @@ void llg_rhs_jtimes(double *jtn, double *m, double *h, double *mp, double *hp, d
106203

107204
/*
108205
if (default_c>0){
109-
jtn[i] += default_c *((1-mm)*mp[i]-2*mmp*m[i]);
110-
jtn[j] += default_c *((1-mm)*mp[j]-2*mmp*m[j]);
111-
jtn[k] += default_c *((1-mm)*mp[k]-2*mmp*m[k]);
206+
jtn[i] += default_c *((1-mm)*mp[i]-2*mmp*m[i]);
207+
jtn[j] += default_c *((1-mm)*mp[j]-2*mmp*m[j]);
208+
jtn[k] += default_c *((1-mm)*mp[k]-2*mmp*m[k]);
112209
}*/
113210

114211
}
@@ -118,16 +215,16 @@ void llg_rhs_jtimes(double *jtn, double *m, double *h, double *mp, double *hp, d
118215

119216
void llg_s_rhs(double *dm_dt, double *m, double *h, double *alpha, double *chi, double gamma, int n) {
120217

121-
int i, j, k;
218+
int i, j, k;
122219

123-
double mth0, mth1, mth2;
124-
double coeff = - gamma;
220+
double mth0, mth1, mth2;
221+
double coeff = - gamma;
125222
double mm, c;
126223
double hpi,hpj,hpk,mh;
127224

128-
for (i = 0; i < n; i++) {
129-
j = i + n;
130-
k = j + n;
225+
for (i = 0; i < n; i++) {
226+
j = i + n;
227+
k = j + n;
131228

132229
if (chi[i] > 0.0){
133230
mth0 = coeff * (m[j] * h[k] - m[k] * h[j]);
@@ -147,7 +244,7 @@ void llg_s_rhs(double *dm_dt, double *m, double *h, double *alpha, double *chi,
147244

148245
}else{ //do LL equation, similiar to the function of llg_rhs
149246

150-
mm = m[i] * m[i] + m[j] * m[j] + m[k] * m[k];
247+
mm = m[i] * m[i] + m[j] * m[j] + m[k] * m[k];
151248
mh = m[i]*h[i] + m[j]*h[j] + m[k]*h[k];
152249

153250
hpi = mm*h[i] - mh*m[i];
@@ -163,6 +260,6 @@ void llg_s_rhs(double *dm_dt, double *m, double *h, double *alpha, double *chi,
163260
dm_dt[k] = coeff*(mth2 - hpk * alpha[i]);
164261
}
165262

166-
}
263+
}
167264

168265
}

0 commit comments

Comments
 (0)