Skip to content

Commit

Permalink
Made energy class setup to depend on the micro sim class. Updated tes…
Browse files Browse the repository at this point in the history 
…t. Moved exchange class functions to micro exch file
  • Loading branch information
@davidcortesortuno
davidcortesortuno committed Jun 12, 2019
1 parent e8b15b0 commit 2f19fcf
Show file tree
Hide file tree
Showing 6 changed files with 231 additions and 220 deletions.
28 changes: 9 additions & 19 deletions fidimag/common/c_clib.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -249,11 +249,6 @@ cdef extern from "c_energy.h":

void compute_field(double t)
double compute_energy()
void setup(int nx, int ny, int nz, double dx, double dy, double dz,
double unit_length, double *spin, double *Ms, double *Ms_inv,
double *coordinates, int *ngbs,
double *energy, double *field
)

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

void init(double *A)
void setup(double * A, MicroSim * sim)
# double *A

# cdef extern from "c_energy.cpp":
Expand Down Expand Up @@ -304,17 +299,6 @@ cdef class PyEnergy:
def compute_energy(self, time):
return self.thisptr.compute_energy()

def setup(self, nx, ny, nz, dx, dy, dz, unit_length,
double [:] spin, double [:] Ms, double [:] Ms_inv,
double [:, :] coordinates, int [:, :] neighbours,
double [:] energy, double [:] field):

return self.thisptr.setup(nx, ny, nz, dx, dy, dz, unit_length,
&spin[0], &Ms[0], &Ms_inv[0],
&coordinates[0, 0], &neighbours[0, 0],
&energy[0], &field[0]
)

def add_interaction_to_sim(self, PyMicroSim sim):
sim.thisptr.add_interaction(<void *> self.thisptr,
self.thisptr.interaction_id)
Expand All @@ -326,11 +310,14 @@ cdef class PyEnergy:
cdef class PyExchangeEnergy(PyEnergy):
cdef ExchangeEnergy *_thisptr
# Try cinit:
def __cinit__(self, double [:] A):
def __cinit__(self, double [:] A, PyMicroSim sim):
print("In Python B")

self._thisptr = self.thisptr = new ExchangeEnergy()
self._thisptr.init(&A[0])
self._thisptr.setup(&A[0], sim.thisptr)

def setup(self, double [:] A, PyMicroSim sim):
self._thisptr.setup(&A[0], sim.thisptr)

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


cdef class PyMicroSim(object):
"""
Wrapper for the C++ MicroSim simulation class
"""
cdef MicroSim *thisptr
# Try cinit:
def __cinit__(self):
Expand Down
32 changes: 24 additions & 8 deletions native/include/c_energy.h
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
#pragma once
#include<cmath>
#include<iostream>
#include "c_micro_sim.h"

class Energy {
public:
Expand All @@ -18,25 +20,39 @@ class Energy {
double *coordinates;
int *ngbs;
double compute_energy();
void setup(int nx, int ny, int nz, double dx, double dy, double dz,
double unit_length, double *spin, double *Ms, double *Ms_inv,
double *coordinates, int *ngbs,
double *energy, double *field
);
// Will get the parameters from a simulation class
void _setup(MicroSim * sim);
virtual void compute_field(double t) {};
};


class ExchangeEnergy : public Energy {
public:
ExchangeEnergy() {
std::cout << "Instatiating ExchangeEnergy class; at " << this << "\n";
this->interaction_id = 1;
};
~ExchangeEnergy() {std::cout << "Killing B\n";};
void init(double *A) {
this->set_up = false;
~ExchangeEnergy() {std::cout << "Killing Exchange\n";};
void setup(double *A, MicroSim * sim) {
_setup(sim);
this->A = A;
}
double *A;
void compute_field(double t);
};


// class AnisotropyEnergy : public Energy {
// public:
// AnisotropyEnergy() {
// std::cout << "Instatiating AnisotropyEnergy class; at " << this << "\n";
// this->interaction_id = 1;
// };
// ~AnisotropyEnergy() {std::cout << "Killing Anisotropy\n";};
// void init(double *Ku) {
// this->set_up = false;
// this->Ku = Ku;
// }
// double *A;
// void compute_field(double t);
// };
1 change: 1 addition & 0 deletions native/include/c_micro_sim.h
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
#pragma once
#include<cmath>
#include<iostream>
#include<vector>
Expand Down
197 changes: 17 additions & 180 deletions native/src/c_energy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,189 +9,26 @@ double Energy::compute_energy() {
return sum * (dx * dy * dz * std::pow(unit_length, 3));
}

void Energy::_setup(MicroSim * sim) {

void Energy::setup(int nx, int ny, int nz, double dx, double dy, double dz,
double unit_length, double *spin, double *Ms, double *Ms_inv,
double *coordinates, int *ngbs,
double *energy, double *field
) {

this->nx = nx;
this->ny = ny;
this->nz = nz;
this->n = nx * ny * nz;
this->dx = dx;
this->dy = dy;
this->dz = dz;
this->unit_length = unit_length;
this->nx = sim->nx;
this->ny = sim->ny;
this->nz = sim->nz;
this->n = sim->nx * ny * nz;
this->dx = sim->dx;
this->dy = sim->dy;
this->dz = sim->dz;
this->unit_length = sim->unit_length;

// Arrays
this->spin = spin;
this->Ms = Ms;
this->Ms_inv = Ms_inv;
this->coordinates = coordinates;
this->ngbs = ngbs;
this->energy = energy;
this->field = field;
this->spin = sim->spin;
this->Ms = sim->Ms;
this->Ms_inv = sim->Ms_inv;
this->coordinates = sim->coordinates;
this->ngbs = sim->ngbs;
this->energy = sim->energy;
this->field = sim->field;
// this->pins = sim->pins;

set_up = true;
}

void ExchangeEnergy::compute_field(double t) {
/* Compute the micromagnetic exchange field and energy using the
* matrix of neighbouring spins and a second order approximation
* for the derivative
*
* Ms_inv :: Array with the (1 / Ms) values for every mesh node.
* The values are zero for points with Ms = 0 (no material)
*
* A :: Exchange constant
*
* dx, dy, dz :: Mesh spacings in the corresponding directions
*
* n :: Number of mesh nodes
*
* ngbs :: The array of neighbouring spins, which has (6 * n)
* entries. Specifically, it contains the indexes of
* the neighbours of every mesh node, in the following order:
* -x, +x, -y, +y, -z, +z
*
* Thus, the array is like:
* | 0-x, 0+x, 0-y, 0+y, 0-z, 0+z, 1-x, 1+x, 1-y, ... |
* i=0 i=1 ...
*
* where 0-y is the index of the neighbour of the 0th spin,
* in the -y direction, for example. The index value for a
* neighbour where Ms = 0, is evaluated as -1. The array
* automatically gives periodic boundaries.
*
* A basic example is a 3 x 3 two dimensional mesh with PBCs
* in the X and Y direction:
*
* +-----------+
* | 6 | 7 | 8 |
* +-----------+
* | 3 | 4 | 5 |
* +-----------+
* | 0 | 1 | 2 |
* +-----------+
*
* so, the first 6 entries (neighbours of the 0th mesh node)
* of the array would be: [ 2 1 6 3 -1 -1 ... ]
* (-1 since there is no material in +-z, and a '2' first,
* since it is the left neighbour which is the PBC in x, etc..)
*
* For the exchange computation, the field is defined as:
* H_ex = (2 * A / (mu0 * Ms)) * nabla^2 (mx, my, mz)
*
* Therefore, for the i-th mesh node (spin), we approximate the
* derivatives as:
* nabla^2 mx = (1 / dx^2) * ( spin[i-x] - 2 * spin[i] + spin[i+x] ) +
* (1 / dy^2) * ( spin[i-y] - 2 * spin[i] + spin[i+y] ) +
* (1 / dz^2) * ( spin[i-z] - 2 * spin[i] + spin[i+z] )
*
* Where i-x is the neighbour in the -x direction. This is similar
* for my and mz.
* We can notice that the sum is the same if we do:
* ( spin[i-x] - spin[i] ) + ( spin[i+x] - spin[i] )
* so we can iterate through the neighbours and perform the sum with the
* corresponding coefficient 1 /dx, 1/dy or 1/dz
*
* The *m array contains the spins as:
* [mx0, my0, mz0, mx1, my1, mz1, mx2, ...]
* so if we want the starting position of the magnetisation for the
* i-th spin, we only have to do (3 * i) for mx, (3 * i + 1) for my
* and (3 * i + 2) for mz
*
*
* IMPORTANT: The ex field usually has the structure:
* 2 * A / (mu0 Ms ) * (Second derivative of M)
* When discretising the derivative, it carries a "2" in the
* denominator which we "cancel" with the "2" in the prefactor,
* hence we do not put it explicitly in the calculations
*
* So, when computing the energy: (-1/2) * mu * Ms * H_ex
* we only put the 0.5 factor and don't worry about the "2"s in the
* field
*
*/

/* Here we iterate through every mesh node */
for (int i = 0; i < n; i++) {
/* Define the coefficients */
double ax = 2 * A[i] / (dx * dx);
double ay = 2 * A[i] / (dy * dy);
double az = 2 * A[i] / (dz * dz);

double fx = 0, fy = 0, fz = 0;
int idnm = 0; // Index for the magnetisation matrix
int idn = 6 * i; // index for the neighbours

/* Set a zero field for sites without magnetic material */
if (Ms_inv[i] == 0.0){
field[3 * i] = 0;
field[3 * i + 1] = 0;
field[3 * i + 2] = 0;
continue;
}

/* Here we iterate through the neighbours */
for (int j = 0; j < 6; j++) {
/* Remember that index=-1 is for sites without material */
if (ngbs[idn + j] >= 0) {
/* Magnetisation of the neighbouring spin since ngbs gives
* the neighbour's index */
idnm = 3 * ngbs[idn + j];

/* Check that the magnetisation of the neighbouring spin
* is larger than zero */
if (Ms_inv[ngbs[idn + j]] > 0){

/* Neighbours in the -x and +x directions
* giving: ( spin[i-x] - spin[i] ) + ( spin[i+x] - spin[i] )
* when ngbs[idn + j] > 0 for j = 0 and j=1
* If, for example, there is no
* neighbour at -x (j=0) in the 0th node (no PBCs),
* the second derivative would only be avaluated as:
* (1 / dx * dx) * ( spin[i+x] - spin[i] )
* which, according to
* [M.J. Donahue and D.G. Porter; Physica B, 343, 177-183 (2004)]
* when performing the integration of the energy, we still
* have error of the order O(dx^2)
* This same applies for the other directions
*/
if (j == 0 || j == 1) {
fx += ax * (spin[idnm] - spin[3 * i]);
fy += ax * (spin[idnm + 1] - spin[3 * i + 1]);
fz += ax * (spin[idnm + 2] - spin[3 * i + 2]);
}
/* Neighbours in the -y and +y directions */
else if (j == 2 || j == 3) {
fx += ay * (spin[idnm] - spin[3 * i]);
fy += ay * (spin[idnm + 1] - spin[3 * i + 1]);
fz += ay * (spin[idnm + 2] - spin[3 * i + 2]);
}
/* Neighbours in the -z and +z directions */
else if (j == 4 || j == 5) {
fx += az * (spin[idnm] - spin[3 * i]);
fy += az * (spin[idnm + 1] - spin[3 * i + 1]);
fz += az * (spin[idnm + 2] - spin[3 * i + 2]);
}
else {
continue;
}
}
}
}

/* Energy as: (-mu0 * Ms / 2) * [ H_ex * m ] */
energy[i] = -0.5 * (fx * spin[3 * i] + fy * spin[3 * i + 1]
+ fz * spin[3 * i + 2]);

/* Update the field H_ex which has the same structure than *m */
field[3 * i] = fx * Ms_inv[i] * MU_0_INV;
field[3 * i + 1] = fy * Ms_inv[i] * MU_0_INV;
field[3 * i + 2] = fz * Ms_inv[i] * MU_0_INV;

}
}
Loading

0 comments on commit 2f19fcf

Please sign in to comment.