Added new micromagnetic Anisotropy and DMI energy classes

fidimag/common/c_clib.pyx

@@ -266,10 +266,17 @@ cdef extern from "c_energy.h":
 
     cdef cppclass ExchangeEnergy(Energy):
         ExchangeEnergy() except +
-
         void setup(double * A, MicroSim * sim)
         # double *A
 
+    cdef cppclass AnisotropyEnergy(Energy):
+        AnisotropyEnergy() except +
+        void setup(double * Ku, double * axis, MicroSim * sim)
+
+    cdef cppclass DMIEnergy(Energy):
+        ExchangeEnergy() except +
+        void setup(double * D, double * dmi_vector, int n_dmis, MicroSim * sim)
+
 # cdef extern from "c_energy.cpp":
 #     pass
 
@@ -311,7 +318,6 @@ cdef class PyExchangeEnergy(PyEnergy):
     cdef ExchangeEnergy *_thisptr
     # Try cinit:
     def __cinit__(self, double [:] A, PyMicroSim sim):
-        print("In Python B")
 
         self._thisptr = self.thisptr = new ExchangeEnergy()
         self._thisptr.setup(&A[0], sim.thisptr)
@@ -323,6 +329,37 @@ cdef class PyExchangeEnergy(PyEnergy):
         del self._thisptr
 
 
+cdef class PyAnisotropyEnergy(PyEnergy):
+    cdef AnisotropyEnergy *_thisptr
+    # Try cinit:
+    def __cinit__(self, double [:] Ku, double [:] axis, PyMicroSim sim):
+
+        self._thisptr = self.thisptr = new AnisotropyEnergy()
+        self._thisptr.setup(&Ku[0], &axis[0], sim.thisptr)
+
+    def setup(self, double [:] Ku, double [:] axis, PyMicroSim sim):
+        self._thisptr.setup(&Ku[0], &axis[0], sim.thisptr)
+
+    def __dealloc__(self):
+        del self._thisptr
+
+
+cdef class PyDMIEnergy(PyEnergy):
+    cdef DMIEnergy *_thisptr
+    # Try cinit:
+    def __cinit__(self, double [:] D, double [:] dmi_vector, n_dmis,
+                  PyMicroSim sim):
+
+        self._thisptr = self.thisptr = new DMIEnergy()
+        self._thisptr.setup(&D[0], &dmi_vector[0], n_dmis, sim.thisptr)
+
+    def setup(self, double [:] D, double [:] dmi_vector, n_dmis, PyMicroSim sim):
+        self._thisptr.setup(&D[0], &dmi_vector[0], n_dmis, sim.thisptr)
+
+    def __dealloc__(self):
+        del self._thisptr
+
+
 # Simulation class ------------------------------------------------------------
 
 cdef extern from "c_micro_sim.h":

native/include/c_energy.h

@@ -33,26 +33,48 @@ class ExchangeEnergy : public Energy {
         this->interaction_id = 1;
     };
     ~ExchangeEnergy() {std::cout << "Killing Exchange\n";};
+    double *A;
     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);
-// };
+class AnisotropyEnergy : public Energy {
+public:
+    AnisotropyEnergy() {
+        std::cout << "Instatiating AnisotropyEnergy class; at " << this << "\n";
+        this->interaction_id = 1;
+    };
+    ~AnisotropyEnergy() {std::cout << "Killing Anisotropy\n";};
+    double *Ku;
+    double *axis;
+    void setup(double *Ku, double *axis, MicroSim * sim) {
+        _setup(sim);
+        this->Ku = Ku;
+        this->axis = axis;
+    }
+    void compute_field(double t);
+};
+
+
+class DMIEnergy : public Energy {
+public:
+    DMIEnergy() {
+        std::cout << "Instatiating DMIEnergy class; at " << this << "\n";
+        this->interaction_id = 1;
+    };
+    ~DMIEnergy() {std::cout << "Killing DMI\n";};
+    double *D;
+    double *dmi_vector;
+    int n_dmis;
+    void setup(double * D, double *dmi_vector, int n_dmis, MicroSim * sim) {
+        _setup(sim);
+        this->D = D;
+        this->dmi_vector = dmi_vector;
+        this->n_dmis = n_dmis;
+    }
+    void compute_field(double t);
+};

native/src/m_anis.cpp

@@ -1,4 +1,8 @@
 #include "m_clib.h"
+#include "c_energy.h"
+#include "c_constants.h"
+
+// TO BE DEPRECATED::
 
 void compute_uniaxial_anis(double * m, double * field, double * energy, double * Ms_inv, 
 			   double * Ku, double * axis, int nx, int ny, int nz) {
@@ -27,3 +31,28 @@ void compute_uniaxial_anis(double * m, double * field, double * energy, double *
   }
 }
 
+// ----------------------------------------------------------------------------
+
+void AnisotropyEnergy::compute_field(double t) {
+
+    for (int i = 0; i < n; i++) {
+
+        int j = 3 * i;
+
+        if (Ms_inv[i] == 0.0){
+          field[j] = 0;
+          field[j + 1] = 0;
+          field[j + 2] = 0;
+          energy[i] = 0;
+          continue;
+        }
+
+        double m_u = spin[j] * axis[j] + spin[j + 1] * axis[j + 1] + spin[j + 2] * axis[j + 2];
+
+        field[j]     = 2 * Ku[i] * m_u * Ms_inv[i] * MU0_INV * axis[j];
+        field[j + 1] = 2 * Ku[i] * m_u * Ms_inv[i] * MU0_INV * axis[j + 1];
+        field[j + 2] = 2 * Ku[i] * m_u * Ms_inv[i] * MU0_INV * axis[j + 2];
+
+        energy[i] = Ku[i] * (1 - m_u * m_u); 
+  }
+}

native/src/m_dmi.cpp

@@ -3,6 +3,8 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <math.h>
+#include "c_energy.h"
+#include "c_constants.h"
 
 /* Functions to Compute the micromagnetic Dzyaloshinskii Moriya interaction
  *  field and energy using the
@@ -364,3 +366,100 @@ void dmi_field(double * m, double * field,
         field[3 * i + 2] = fz * Ms_inv[i] * MU0_INV;
     }
 }
+
+// ----------------------------------------------------------------------------
+// Using the new Energy class::
+
+
+void DMIEnergy::compute_field(double t) {
+
+    /* These are for the DMI prefactor or coefficient */
+    double dxs[6] = {dx, dx, dy, dy, dz, dz};
+
+    /* Here we iterate through every mesh node */
+    for (int i = 0; i < n; i++) {
+        double DMIc;
+        double fx = 0, fy = 0, fz = 0;
+        int idnm;     // 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;
+            energy[i] = 0;
+            continue;
+        }
+
+        /* Here we iterate through the neighbours. Remember:
+         * j = 0, 1, 2, 3, 4, 5  --> -x, +x, -y, +y, -z, +z
+         Previously we had n_dmi_neighbours as a parameter,
+         but now we just skip if the vector entries are zero
+         instead, because we can have cases such as in D_n
+         where the D1 component has -x,+x,-y,+y,0,0
+         and the D2 component has 0,0,0,0,-z,+z.
+        */
+        for (int j = 0; j < 6; j++) {
+            /* Add the DMI field x times for every DMI constant */
+            for (int k = 0; k < n_dmis; k++) {
+
+                // starting index of the DMI vector for this neighbour (j)
+                // (remember we have 18 comps of dmi_vector per DMI constant)
+                int ngbr_idx_D = k * 18 + 3 * j;
+
+                /* We skip the neighbour if
+                   (a) it doesn't exist (ngbs[idn + j] = -1)
+                   (b) there is no material there
+                   (c) DMI value is zero there
+                */
+                if ((ngbs[idn + j] != -1) && (Ms_inv[ngbs[idn + j]] !=  0)) {
+                /* We do here:  -(D / dx_i) * ( r_{ij} X M_{j} )
+                 * The cross_i function gives the i component of the
+                 * cross product. The coefficient is computed according
+                 * to the DMI strength of the current lattice site.
+                 * For the denominator, for example, if j=2 or 3, then
+                 * dxs[j] = dy
+                 */
+
+                    /* check the DMI vector exists */
+                    if (dmi_vector[ngbr_idx_D    ] != 0 ||
+                        dmi_vector[ngbr_idx_D + 1] != 0 ||
+                        dmi_vector[ngbr_idx_D + 2] != 0   ) {
+                        /* Notice the x component of the cross product of +-x
+                         * times anything is zero (similar for the other comps) */
+
+                        // Get DMI constant from neighbour and scale
+                        // DMI components are in consecutive order per neighbour,
+                        // i.e. if 2 DMI consts: [D0 D1 D0 D1 .... ]
+                        DMIc = -D[n_dmis * ngbs[idn + j] + k] / dxs[j];
+
+                        idnm = 3 * ngbs[idn + j]; // index for magnetisation
+
+                        fx += DMIc * cross_x(dmi_vector[ngbr_idx_D],
+                                             dmi_vector[ngbr_idx_D + 1],
+                                             dmi_vector[ngbr_idx_D + 2],
+                                             spin[idnm], spin[idnm + 1], spin[idnm + 2]);
+                        fy += DMIc * cross_y(dmi_vector[ngbr_idx_D],
+                                             dmi_vector[ngbr_idx_D + 1],
+                                             dmi_vector[ngbr_idx_D + 2],
+                                             spin[idnm], spin[idnm + 1], spin[idnm + 2]);
+                        fz += DMIc * cross_z(dmi_vector[ngbr_idx_D],
+                                             dmi_vector[ngbr_idx_D + 1],
+                                             dmi_vector[ngbr_idx_D + 2],
+                                             spin[idnm], spin[idnm + 1], spin[idnm + 2]);
+
+                    }
+                }
+            }  // Close for loop through n of DMI constants
+        }  // Close for loop through neighbours per mesh site
+
+        /* Energy as: (-mu0 * Ms / 2) * [ H_dmi * m ]   */
+        energy[i] = -0.5 * (fx * spin[3 * i] + fy * spin[3 * i + 1]
+        		+ fz * spin[3 * i + 2]);
+
+        /* Update the field H_dmi which has the same structure than *m */
+        field[3 * i]     = fx * Ms_inv[i] * MU0_INV;
+        field[3 * i + 1] = fy * Ms_inv[i] * MU0_INV;
+        field[3 * i + 2] = fz * Ms_inv[i] * MU0_INV;
+    }
+}