Add 4th Order UniaxialAnisotropy

Author: rpep

examples/micromagnetic/skyrmion/Ku_H/main.py

@@ -3,10 +3,10 @@
 import matplotlib.pyplot as plt
 
 import numpy as np
-from micro import Sim
-from common import CuboidMesh
-from micro import UniformExchange, Demag, DMI, UniaxialAnisotropy
-from micro import Zeeman, TimeZeeman
+from fidimag.micro import Sim
+from fidimag.common import CuboidMesh
+from fidimag.micro import UniformExchange, Demag, DMI, UniaxialAnisotropy
+from fidimag.micro import Zeeman, TimeZeeman
 from fidimag.common.fileio import DataReader
 
 mu0 = 4 * np.pi * 1e-7

fidimag/micro/__init__.py

@@ -3,7 +3,7 @@
 from .exchange import UniformExchange
 from .exchange_rkky import ExchangeRKKY
 from .zeeman import Zeeman, TimeZeeman
-from .anisotropy import UniaxialAnisotropy
+from .anisotropy import UniaxialAnisotropy, UniaxialAnisotropy4
 from .dmi import DMI
 from .baryakhtar import LLBar, LLBarFull
 from .simple_demag import SimpleDemag

fidimag/micro/anisotropy.py

@@ -38,3 +38,47 @@ def compute_field(self, t=0, spin=None):
                                             self.ny,
                                             self.nz)
         return self.field
+
+
+class UniaxialAnisotropy4(Energy):
+
+    """
+    4th Order Uniaxial Anisotropy
+    """
+
+    def __init__(self, K1, K2, axis=(1, 0, 0), name='UAnisotropy4'):
+        self.K1 = K1
+        self.K2 = K2
+        self.name = name
+        self.jac = True
+        self.axis = axis
+
+    def setup(self, mesh, spin, Ms, Ms_inv):
+        super(UniaxialAnisotropy4, self).setup(mesh, spin, Ms, Ms_inv)
+
+        self._K1 = helper.init_scalar(self.K1, self.mesh)
+        self._K2 = helper.init_scalar(self.K2, self.mesh)
+        self._axis = helper.init_vector(self.axis, self.mesh, 3, norm=True)
+
+    def compute_field(self, t=0, spin=None):
+        if spin is not None:
+            m = spin
+        else:
+            m = self.spin
+        micro_clib.compute_anisotropy4_micro(m,
+                                             self.field,
+                                             self.energy,
+                                             self.Ms_inv,
+                                             self._K1,
+                                             self._K2,
+                                             self._axis,
+                                             self.nx,
+                                             self.ny,
+                                             self.nz)
+
+
+
+
+
+        return self.field
+

fidimag/micro/lib/anis.c

@@ -30,3 +30,52 @@ void compute_uniaxial_anis(double *restrict m, double *restrict field, double *r
 
 }
 
+
+void compute_uniaxial4_anis(double *restrict m, double *restrict field, double *restrict energy, double *restrict Ms_inv, 
+    double *restrict K1, double *restrict K2, double *restrict axis, int nx, int ny, int nz) {
+    
+    int n = nx * ny * nz;
+
+    // Follows calculation of OOMMF extension by Hans and Richard Boardman
+    // http://www.soton.ac.uk/~fangohr/software/oxs_uniaxial4/download/uniaxialanisotropy4.cc
+
+    #pragma omp parallel for
+    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 k1 = K1[i];
+        double k2 = K2[i];
+
+        double field_mult1 = MU0_INV * 2.0 * k1 * Ms_inv[i];
+        double field_mult2 = MU0_INV * 4.0 * k2 * Ms_inv[i];
+        double m_dot_u = m[j] * axis[j] + m[j + 1] * axis[j + 1] + m[j + 2] * axis[j + 2];
+
+        if (k1 <= 0) {
+            field[j + 0] = (field_mult1*m_dot_u) * axis[j + 0] + (field_mult2 * m_dot_u*m_dot_u*m_dot_u) * axis[j + 0];
+            field[j + 1] = (field_mult1*m_dot_u) * axis[j + 1] + (field_mult2 * m_dot_u*m_dot_u*m_dot_u) * axis[j + 1];
+            field[j + 2] = (field_mult1*m_dot_u) * axis[j + 2] + (field_mult2 * m_dot_u*m_dot_u*m_dot_u) * axis[j + 2];
+            energy[i] = -k1*m_dot_u*m_dot_u - k2*m_dot_u*m_dot_u*m_dot_u*m_dot_u;
+        }
+
+        else {
+            double u_x_m[3];
+            u_x_m[0] = cross_x(axis[j], axis[j+1], axis[j+2], m[j], m[j+1], m[j+2]);
+            u_x_m[1] = cross_y(axis[j], axis[j+1], axis[j+2], m[j], m[j+1], m[j+2]);
+            u_x_m[2] = cross_z(axis[j], axis[j+1], axis[j+2], m[j], m[j+1], m[j+2]);
+            double u_x_m_mag2 = u_x_m[1]*u_x_m[1] + u_x_m[1]*u_x_m[1] + u_x_m[2]*u_x_m[2];
+            field[j + 0] = (field_mult1*m_dot_u) * axis[j + 0] + (field_mult2*m_dot_u*m_dot_u*m_dot_u) * axis[j + 0];
+            field[j + 1] = (field_mult1*m_dot_u) * axis[j + 1] + (field_mult2*m_dot_u*m_dot_u*m_dot_u) * axis[j + 1];
+            field[j + 2] = (field_mult1*m_dot_u) * axis[j + 2] + (field_mult2*m_dot_u*m_dot_u*m_dot_u) * axis[j + 2];
+            energy[i] = (k1 + 2*k2)*u_x_m_mag2 - k2*u_x_m_mag2*u_x_m_mag2;
+        }
+    }
+
+}

fidimag/micro/lib/micro_clib.h

@@ -28,5 +28,8 @@ void compute_exch_field_rkky_micro(double *m, double *field, double *energy, dou
 void compute_uniaxial_anis(double *restrict m, double *restrict field, double *restrict energy, double *restrict Ms_inv,
 	double *restrict Ku, double *restrict axis, int nx, int ny, int nz);
 
+void compute_uniaxial4_anis(double *restrict m, double *restrict field, double *restrict energy, double *restrict Ms_inv, 
+    double *restrict K1, double *restrict K2, double *restrict axis, int nx, int ny, int nz);
+
 double skyrmion_number(double *restrict spin, double *restrict charge,
                        int nx, int ny, int nz, int *restrict ngbs);

fidimag/micro/lib/micro_clib.pyx

@@ -21,10 +21,20 @@ cdef extern from "micro_clib.h":
                                double *Ku, double *axis,
                                int nx, int ny, int nz)
 
+
+    void compute_uniaxial4_anis(double *m, double *field,
+                               double *energy, double *Ms_inv,
+                               double *K1, double *K2,
+                               double *axis,
+                               int nx, int ny, int nz)
+
+
     double skyrmion_number(double *m, double *charge,
                            int nx, int ny, int nz, int *ngbs)
 
 
+
+
 def compute_exchange_field_micro(double [:] m,
                                  double [:] field,
                                  double [:] energy,
@@ -74,6 +84,20 @@ def compute_anisotropy_micro(double [:] m,
     compute_uniaxial_anis(&m[0], &field[0], &energy[0], &Ms_inv[0],
                           &Ku[0], &axis[0], nx, ny, nz)
 
+
+def compute_anisotropy4_micro(double [:] m,
+                             double [:] field,
+                             double [:] energy,
+                             double [:] Ms_inv,
+                             double [:] K1,
+                             double [:] K2,
+                             double [:] axis,
+                             nx, ny, nz):
+
+    compute_uniaxial4_anis(&m[0], &field[0], &energy[0], &Ms_inv[0],
+                          &K1[0], &K2[0], &axis[0], nx, ny, nz)
+
+
 def compute_skyrmion_number(double [:] m,
                             double [:] charge,
                             nx, ny, nz,