|
1 | 1 | #include "common_clib.h"
|
2 | 2 | #include "math.h"
|
3 | 3 |
|
4 |
| -void sd_update_spin (double *spin, double *spin_last, |
| 4 | +void sd_update_spin (double *spin, double *spin_last, double *magnetisation, |
5 | 5 | double *mxH, double *mxmxH, double *mxmxH_last, double *tau,
|
6 | 6 | int* pins, int n) {
|
7 | 7 |
|
8 | 8 | // Update the field -------------------------------------------------------
|
9 | 9 | #pragma omp parallel for
|
10 | 10 | for (int i = 0; i < n; i++) {
|
11 | 11 |
|
12 |
| - int spin_idx; |
13 |
| - double mxH_sq; |
14 |
| - double factor_plus; |
15 |
| - double factor_minus; |
16 |
| - double new_spin[3]; |
| 12 | + if (magnetisation[i] > 0.0) { |
17 | 13 |
|
18 |
| - spin_idx = 3 * i; |
| 14 | + int spin_idx; |
| 15 | + double mxH_sq; |
| 16 | + double factor_plus; |
| 17 | + double factor_minus; |
| 18 | + double new_spin[3]; |
19 | 19 |
|
20 |
| - // Copy spin from previous (counter) step |
21 |
| - for (int j = 0; j < 3; j++) { |
22 |
| - spin_last[spin_idx + j] = spin[spin_idx + j]; |
23 |
| - } |
| 20 | + spin_idx = 3 * i; |
24 | 21 |
|
25 |
| - mxH_sq = 0; |
26 |
| - for (int j = 0; j < 3; j++) { |
27 |
| - mxH_sq += mxH[spin_idx + j] * mxH[spin_idx + j]; |
28 |
| - } |
| 22 | + // Copy spin from previous (counter) step |
| 23 | + for (int j = 0; j < 3; j++) { |
| 24 | + spin_last[spin_idx + j] = spin[spin_idx + j]; |
| 25 | + } |
| 26 | + |
| 27 | + mxH_sq = 0; |
| 28 | + for (int j = 0; j < 3; j++) { |
| 29 | + mxH_sq += mxH[spin_idx + j] * mxH[spin_idx + j]; |
| 30 | + } |
29 | 31 |
|
30 |
| - factor_plus = 4 + tau[i] * tau[i] * mxH_sq; |
31 |
| - factor_minus = 4 - tau[i] * tau[i] * mxH_sq; |
| 32 | + factor_plus = 4 + tau[i] * tau[i] * mxH_sq; |
| 33 | + factor_minus = 4 - tau[i] * tau[i] * mxH_sq; |
32 | 34 |
|
33 |
| - for (int j = 0; j < 3; j++) { |
34 |
| - new_spin[j] = factor_minus * spin[spin_idx + j] |
35 |
| - - 4 * tau[i] * mxmxH[spin_idx + j]; |
36 |
| - new_spin[j] = new_spin[j] / factor_plus; |
| 35 | + for (int j = 0; j < 3; j++) { |
| 36 | + new_spin[j] = factor_minus * spin[spin_idx + j] |
| 37 | + - 4 * tau[i] * mxmxH[spin_idx + j]; |
| 38 | + new_spin[j] = new_spin[j] / factor_plus; |
37 | 39 |
|
38 |
| - spin[spin_idx + j] = new_spin[j]; |
39 |
| - } |
40 |
| - } |
| 40 | + spin[spin_idx + j] = new_spin[j]; |
| 41 | + } |
| 42 | + } // close if Ms or mu_s > 0 |
| 43 | + } // close for |
41 | 44 | normalise(spin, pins, n);
|
42 | 45 | }
|
43 | 46 |
|
44 |
| -void sd_compute_step (double *spin, double *spin_last, double *field, |
| 47 | +void sd_compute_step (double *spin, double *spin_last, double *magnetisation, double *field, |
45 | 48 | double *mxH, double *mxmxH, double *mxmxH_last, double *tau,
|
46 | 49 | int *pins, int n, int counter, double tmin, double tmax) {
|
47 | 50 |
|
48 | 51 | #pragma omp parallel for
|
49 | 52 | for (int i = 0; i < n; i++) {
|
50 | 53 |
|
51 |
| - int spin_idx; |
52 |
| - double ds[3], dy[3]; |
53 |
| - double num, den, res, sign; |
54 |
| - |
55 |
| - spin_idx = 3 * i; |
56 |
| - |
57 |
| - // Copy mxmxH from previous (counter) step |
58 |
| - for (int j = 0; j < 3; j++) { |
59 |
| - mxmxH_last[spin_idx + j] = mxmxH[spin_idx + j]; |
60 |
| - } |
61 |
| - |
62 |
| - // Compute the torques |
63 |
| - |
64 |
| - mxH[spin_idx] = cross_x(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], |
65 |
| - field[spin_idx], |
66 |
| - field[spin_idx + 1], |
67 |
| - field[spin_idx + 2]); |
68 |
| - |
69 |
| - mxH[spin_idx + 1] = cross_y(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], |
70 |
| - field[spin_idx], |
71 |
| - field[spin_idx + 1], |
72 |
| - field[spin_idx + 2]); |
73 |
| - |
74 |
| - mxH[spin_idx + 2] = cross_z(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], |
75 |
| - field[spin_idx], |
76 |
| - field[spin_idx + 1], |
77 |
| - field[spin_idx + 2]); |
78 |
| - |
79 |
| - mxmxH[spin_idx] = cross_x(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], |
80 |
| - mxH[spin_idx], mxH[spin_idx + 1], mxH[spin_idx + 2]); |
81 |
| - mxmxH[spin_idx + 1] = cross_y(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], |
82 |
| - mxH[spin_idx], mxH[spin_idx + 1], mxH[spin_idx + 2]); |
83 |
| - mxmxH[spin_idx + 2] = cross_z(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], |
84 |
| - mxH[spin_idx], mxH[spin_idx + 1], mxH[spin_idx + 2]); |
85 |
| - |
86 |
| - // Compute terms for the calculation of the time step tau |
87 |
| - for (int j = 0; j < 3; j++) { |
88 |
| - ds[j] = spin[spin_idx + j] - spin_last[spin_idx + j]; |
89 |
| - dy[j] = mxmxH[spin_idx + j] - mxmxH_last[spin_idx + j]; |
90 |
| - } |
91 |
| - |
92 |
| - num = 0; |
93 |
| - den = 0; |
94 |
| - if (counter % 2 == 0) { |
| 54 | + if (magnetisation[i] > 0.0) { |
| 55 | + int spin_idx; |
| 56 | + double ds[3], dy[3]; |
| 57 | + double num, den, res, sign; |
| 58 | + |
| 59 | + spin_idx = 3 * i; |
| 60 | + |
| 61 | + // Copy mxmxH from previous (counter) step |
95 | 62 | for (int j = 0; j < 3; j++) {
|
96 |
| - num += ds[j] * ds[j]; |
97 |
| - den += ds[j] * dy[j]; |
| 63 | + mxmxH_last[spin_idx + j] = mxmxH[spin_idx + j]; |
98 | 64 | }
|
99 |
| - } |
100 |
| - else { |
| 65 | + |
| 66 | + // Compute the torques |
| 67 | + |
| 68 | + mxH[spin_idx] = cross_x(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], |
| 69 | + field[spin_idx], |
| 70 | + field[spin_idx + 1], |
| 71 | + field[spin_idx + 2]); |
| 72 | + |
| 73 | + mxH[spin_idx + 1] = cross_y(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], |
| 74 | + field[spin_idx], |
| 75 | + field[spin_idx + 1], |
| 76 | + field[spin_idx + 2]); |
| 77 | + |
| 78 | + mxH[spin_idx + 2] = cross_z(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], |
| 79 | + field[spin_idx], |
| 80 | + field[spin_idx + 1], |
| 81 | + field[spin_idx + 2]); |
| 82 | + |
| 83 | + mxmxH[spin_idx] = cross_x(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], |
| 84 | + mxH[spin_idx], mxH[spin_idx + 1], mxH[spin_idx + 2]); |
| 85 | + mxmxH[spin_idx + 1] = cross_y(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], |
| 86 | + mxH[spin_idx], mxH[spin_idx + 1], mxH[spin_idx + 2]); |
| 87 | + mxmxH[spin_idx + 2] = cross_z(spin[spin_idx], spin[spin_idx + 1], spin[spin_idx + 2], |
| 88 | + mxH[spin_idx], mxH[spin_idx + 1], mxH[spin_idx + 2]); |
| 89 | + |
| 90 | + // Compute terms for the calculation of the time step tau |
101 | 91 | for (int j = 0; j < 3; j++) {
|
102 |
| - num += ds[j] * dy[j]; |
103 |
| - den += dy[j] * dy[j]; |
| 92 | + ds[j] = spin[spin_idx + j] - spin_last[spin_idx + j]; |
| 93 | + dy[j] = mxmxH[spin_idx + j] - mxmxH_last[spin_idx + j]; |
| 94 | + } |
| 95 | + |
| 96 | + num = 0; |
| 97 | + den = 0; |
| 98 | + if (counter % 2 == 0) { |
| 99 | + for (int j = 0; j < 3; j++) { |
| 100 | + num += ds[j] * ds[j]; |
| 101 | + den += ds[j] * dy[j]; |
| 102 | + } |
| 103 | + } |
| 104 | + else { |
| 105 | + for (int j = 0; j < 3; j++) { |
| 106 | + num += ds[j] * dy[j]; |
| 107 | + den += dy[j] * dy[j]; |
| 108 | + } |
104 | 109 | }
|
105 |
| - } |
106 |
| - |
107 |
| - // Criteria for the evaluation of tau is in line 96 of: |
108 |
| - //https://github.com/MicroMagnum/MicroMagnum/blob/minimizer/src/magnum/micromagnetics/micro_magnetics_solver.py |
109 |
| - if (den == 0.0) { |
110 |
| - res = tmax; |
111 |
| - } |
112 |
| - else { |
113 |
| - res = num / den; |
114 |
| - } |
115 |
| - |
116 |
| - sign = (res > 0) ? 1 : ((res < 0) ? -1 : 0); |
117 |
| - tau[i] = fmax(fmin(fabs(res), tmax), tmin) * sign; |
118 |
| - |
119 |
| - // In MuMax3 they only define a specific value: |
120 |
| - // https://github.com/mumax/3/blob/master/engine/minimizer.go |
121 |
| - // tau[i] = res; |
122 |
| - } |
| 110 | + |
| 111 | + // Criteria for the evaluation of tau is in line 96 of: |
| 112 | + //https://github.com/MicroMagnum/MicroMagnum/blob/minimizer/src/magnum/micromagnetics/micro_magnetics_solver.py |
| 113 | + if (den == 0.0) { |
| 114 | + res = tmax; |
| 115 | + } |
| 116 | + else { |
| 117 | + res = num / den; |
| 118 | + } |
| 119 | + |
| 120 | + sign = (res > 0) ? 1 : ((res < 0) ? -1 : 0); |
| 121 | + tau[i] = fmax(fmin(fabs(res), tmax), tmin) * sign; |
| 122 | + |
| 123 | + // In MuMax3 they only define a specific value: |
| 124 | + // https://github.com/mumax/3/blob/master/engine/minimizer.go |
| 125 | + // tau[i] = res; |
| 126 | + |
| 127 | + } // Close if Ms or mu_s > 0 |
| 128 | + } // close for |
123 | 129 | }
|
0 commit comments