Skip to content

Commit 2f19fcf

Browse files
Made energy class setup to depend on the micro sim class. Updated test. Moved exchange class functions to micro exch file
1 parent e8b15b0 commit 2f19fcf

File tree

6 files changed

+231
-220
lines changed

6 files changed

+231
-220
lines changed

fidimag/common/c_clib.pyx

Lines changed: 9 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -249,11 +249,6 @@ cdef extern from "c_energy.h":
249249

250250
void compute_field(double t)
251251
double compute_energy()
252-
void setup(int nx, int ny, int nz, double dx, double dy, double dz,
253-
double unit_length, double *spin, double *Ms, double *Ms_inv,
254-
double *coordinates, int *ngbs,
255-
double *energy, double *field
256-
)
257252

258253
# # Not using these variables from Cython:
259254
# bool set_up
@@ -272,7 +267,7 @@ cdef extern from "c_energy.h":
272267
cdef cppclass ExchangeEnergy(Energy):
273268
ExchangeEnergy() except +
274269

275-
void init(double *A)
270+
void setup(double * A, MicroSim * sim)
276271
# double *A
277272

278273
# cdef extern from "c_energy.cpp":
@@ -304,17 +299,6 @@ cdef class PyEnergy:
304299
def compute_energy(self, time):
305300
return self.thisptr.compute_energy()
306301

307-
def setup(self, nx, ny, nz, dx, dy, dz, unit_length,
308-
double [:] spin, double [:] Ms, double [:] Ms_inv,
309-
double [:, :] coordinates, int [:, :] neighbours,
310-
double [:] energy, double [:] field):
311-
312-
return self.thisptr.setup(nx, ny, nz, dx, dy, dz, unit_length,
313-
&spin[0], &Ms[0], &Ms_inv[0],
314-
&coordinates[0, 0], &neighbours[0, 0],
315-
&energy[0], &field[0]
316-
)
317-
318302
def add_interaction_to_sim(self, PyMicroSim sim):
319303
sim.thisptr.add_interaction(<void *> self.thisptr,
320304
self.thisptr.interaction_id)
@@ -326,11 +310,14 @@ cdef class PyEnergy:
326310
cdef class PyExchangeEnergy(PyEnergy):
327311
cdef ExchangeEnergy *_thisptr
328312
# Try cinit:
329-
def __cinit__(self, double [:] A):
313+
def __cinit__(self, double [:] A, PyMicroSim sim):
330314
print("In Python B")
331315

332316
self._thisptr = self.thisptr = new ExchangeEnergy()
333-
self._thisptr.init(&A[0])
317+
self._thisptr.setup(&A[0], sim.thisptr)
318+
319+
def setup(self, double [:] A, PyMicroSim sim):
320+
self._thisptr.setup(&A[0], sim.thisptr)
334321

335322
def __dealloc__(self):
336323
del self._thisptr
@@ -356,6 +343,9 @@ cdef extern from "c_micro_sim.h":
356343

357344

358345
cdef class PyMicroSim(object):
346+
"""
347+
Wrapper for the C++ MicroSim simulation class
348+
"""
359349
cdef MicroSim *thisptr
360350
# Try cinit:
361351
def __cinit__(self):

native/include/c_energy.h

Lines changed: 24 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,7 @@
1+
#pragma once
12
#include<cmath>
23
#include<iostream>
4+
#include "c_micro_sim.h"
35

46
class Energy {
57
public:
@@ -18,25 +20,39 @@ class Energy {
1820
double *coordinates;
1921
int *ngbs;
2022
double compute_energy();
21-
void setup(int nx, int ny, int nz, double dx, double dy, double dz,
22-
double unit_length, double *spin, double *Ms, double *Ms_inv,
23-
double *coordinates, int *ngbs,
24-
double *energy, double *field
25-
);
23+
// Will get the parameters from a simulation class
24+
void _setup(MicroSim * sim);
2625
virtual void compute_field(double t) {};
2726
};
2827

28+
2929
class ExchangeEnergy : public Energy {
3030
public:
3131
ExchangeEnergy() {
3232
std::cout << "Instatiating ExchangeEnergy class; at " << this << "\n";
3333
this->interaction_id = 1;
3434
};
35-
~ExchangeEnergy() {std::cout << "Killing B\n";};
36-
void init(double *A) {
37-
this->set_up = false;
35+
~ExchangeEnergy() {std::cout << "Killing Exchange\n";};
36+
void setup(double *A, MicroSim * sim) {
37+
_setup(sim);
3838
this->A = A;
3939
}
4040
double *A;
4141
void compute_field(double t);
4242
};
43+
44+
45+
// class AnisotropyEnergy : public Energy {
46+
// public:
47+
// AnisotropyEnergy() {
48+
// std::cout << "Instatiating AnisotropyEnergy class; at " << this << "\n";
49+
// this->interaction_id = 1;
50+
// };
51+
// ~AnisotropyEnergy() {std::cout << "Killing Anisotropy\n";};
52+
// void init(double *Ku) {
53+
// this->set_up = false;
54+
// this->Ku = Ku;
55+
// }
56+
// double *A;
57+
// void compute_field(double t);
58+
// };

native/include/c_micro_sim.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
#pragma once
12
#include<cmath>
23
#include<iostream>
34
#include<vector>

native/src/c_energy.cpp

Lines changed: 17 additions & 180 deletions
Original file line numberDiff line numberDiff line change
@@ -9,189 +9,26 @@ double Energy::compute_energy() {
99
return sum * (dx * dy * dz * std::pow(unit_length, 3));
1010
}
1111

12+
void Energy::_setup(MicroSim * sim) {
1213

13-
void Energy::setup(int nx, int ny, int nz, double dx, double dy, double dz,
14-
double unit_length, double *spin, double *Ms, double *Ms_inv,
15-
double *coordinates, int *ngbs,
16-
double *energy, double *field
17-
) {
18-
19-
this->nx = nx;
20-
this->ny = ny;
21-
this->nz = nz;
22-
this->n = nx * ny * nz;
23-
this->dx = dx;
24-
this->dy = dy;
25-
this->dz = dz;
26-
this->unit_length = unit_length;
14+
this->nx = sim->nx;
15+
this->ny = sim->ny;
16+
this->nz = sim->nz;
17+
this->n = sim->nx * ny * nz;
18+
this->dx = sim->dx;
19+
this->dy = sim->dy;
20+
this->dz = sim->dz;
21+
this->unit_length = sim->unit_length;
2722

2823
// Arrays
29-
this->spin = spin;
30-
this->Ms = Ms;
31-
this->Ms_inv = Ms_inv;
32-
this->coordinates = coordinates;
33-
this->ngbs = ngbs;
34-
this->energy = energy;
35-
this->field = field;
24+
this->spin = sim->spin;
25+
this->Ms = sim->Ms;
26+
this->Ms_inv = sim->Ms_inv;
27+
this->coordinates = sim->coordinates;
28+
this->ngbs = sim->ngbs;
29+
this->energy = sim->energy;
30+
this->field = sim->field;
31+
// this->pins = sim->pins;
3632

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

0 commit comments

Comments
 (0)