Merge in OpenMP fixes, etc

Author: rpep

examples/mc/skx.py

@@ -14,7 +14,7 @@ def run(mesh):
     mc = MonteCarlo(mesh, name='test1')
     mc.set_m(random_m)
     J = 50*const.k_B
-    mc.set_options(H=[0,0,0.0], J=J, D=0.27*J, T=4.0)
+    mc.set_options(H=[0,0,1.0], J=J, D=0.27*J, T=4.0)
     mc.run(steps=100000, save_m_steps=None, save_vtk_steps=50000, save_data_steps=50000)
 
 if __name__=='__main__':

examples/micromagnetic/PRB_88_184422/fidimag/prb.py

@@ -32,7 +32,7 @@ def test_prb88_184422():
 
     sim.set_m((0, 0, 1))
 
-    sim.add(UniformExchange(A))
+    sim.add(UniformExchange(A=A))
     sim.add(DMI(-D, type='interfacial'))
     sim.add(UniaxialAnisotropy(K, axis=(0, 0, 1)))
 

examples/micromagnetic/dw_stt/main.py

@@ -3,12 +3,12 @@
 import matplotlib.pyplot as plt
 
 import numpy as np
-from micro import Sim
-from common import CuboidMesh
+from fidimag.micro import Sim
+from fidimag.common import CuboidMesh
 
 # The energies, we can use DMI in a future simulation
-from micro import UniformExchange
-from micro import UniaxialAnisotropy
+from fidimag.micro import UniformExchange
+from fidimag.micro import UniaxialAnisotropy
 # from micro import DMI
 
 mu0 = 4 * np.pi * 1e-7
@@ -21,11 +21,11 @@ def init_m(pos):
     x = pos[0]
 
     if x < 400:
-        return (1, 0, 0)
+        return (-1, 0, 0)
     elif 400 <= x < 500:
         return (0, 1, 1)
     else:
-        return (-1, 0, 0)
+        return (1, 0, 0)
 
 
 def relax_system(mesh):
@@ -38,7 +38,7 @@ def relax_system(mesh):
     sim.driver.alpha = 0.5
     sim.driver.gamma = 2.211e5
     sim.Ms = 8.6e5
-    sim.do_precession = False
+    sim.driver.do_precession = False
 
     # The initial state passed as a function
     sim.set_m(init_m)
@@ -59,8 +59,8 @@ def relax_system(mesh):
     # ONE_DEGREE_PER_NS = 17453292.52
 
     # Start relaxation and save the state in m0.npy
-    sim.relax(dt=1e-14, stopping_dmdt=0.00001, max_steps=5000,
-              save_m_steps=None, save_vtk_steps=None)
+    sim.relax(dt=1e-14, stopping_dmdt=0.01, max_steps=5000,
+              save_m_steps=None, save_vtk_steps=50)
 
     np.save('m0.npy', sim.spin)
 
@@ -106,8 +106,7 @@ def deal_plot(_list, output_file):
 
 
 # THIS NEEDS REVISION
-def deal_plot_dynamics(spin, m_component='mz',
-                       output_file):
+def deal_plot_dynamics(spin, m_component='mz', output_file='output.pdf'):
     """
     The dynamics of the m_component of the i-th spin
     of the chain, from the magnetisation snapshot
@@ -140,10 +139,10 @@ def deal_plot_dynamics(spin, m_component='mz',
 def excite_system(mesh):
 
     # Specify the stt dynamics in the simulation
-    sim = Sim(mesh, name='dyn', driver='llg_stt')
+    sim = Sim(mesh, name='dyn2', driver='llg_stt')
 
     sim.driver.set_tols(rtol=1e-12, atol=1e-14)
-    sim.driver.alpha = 0.05
+    sim.driver.alpha = 0.2
     sim.driver.gamma = 2.211e5
     sim.Ms = 8.6e5
 
@@ -163,16 +162,16 @@ def excite_system(mesh):
 
     # Set the current in the x direction, in A / m
     # beta is the parameter in the STT torque
-    sim.jx = -1e12
-    sim.beta = 1
+    sim.driver.jx = -1e12
+    sim.driver.beta = 0.01
 
     # The simulation will run for 5 ns and save
     # 500 snapshots of the system in the process
     ts = np.linspace(0, 5e-9, 501)
 
     for t in ts:
         print 'time', t
-        sim.run_until(t)
+        sim.driver.run_until(t)
         sim.save_vtk()
         sim.save_m()
 

examples/micromagnetic/vortex/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
-from micro import Zeeman, TimeZeeman
+from fidimag.micro import Sim
+from fidimag.common import CuboidMesh
+from fidimag.micro import UniformExchange, Demag
+from fidimag.micro import Zeeman, TimeZeeman
 from fidimag.common.fileio import DataReader
 
 mu0 = 4 * np.pi * 1e-7
@@ -46,7 +46,7 @@ def relax_system(mesh):
     sim.driver.alpha = 0.5
     sim.driver.gamma = 2.211e5
     sim.Ms = spatial_Ms
-    sim.do_precession = False
+    sim.driver.do_precession = False
 
     sim.set_m(init_m)
     # sim.set_m(np.load('m0.npy'))
@@ -55,7 +55,7 @@ def relax_system(mesh):
     exch = UniformExchange(A=A)
     sim.add(exch)
 
-    demag = Demag()
+    demag = Demag(pbc_2d=True)
     sim.add(demag)
 
     mT = 795.7747154594767
@@ -84,7 +84,7 @@ def excite_system(mesh):
     exch = UniformExchange(A=A)
     sim.add(exch)
 
-    demag = Demag()
+    demag = Demag(pbc_2d=True)
     sim.add(demag)
 
     mT = 795.7747154594767
@@ -110,9 +110,9 @@ def gaussian_fun(t):
 
     mesh = CuboidMesh(nx=80, ny=80, nz=2, dx=2.5, dy=2.5, dz=5.0, unit_length=1e-9)
 
-    # relax_system(mesh)
+    relax_system(mesh)
 
-    excite_system(mesh)
+    #excite_system(mesh)
 
     # apply_field1(mesh)
     # deal_plot()

fidimag/atomistic/lib/mc.c

@@ -165,8 +165,8 @@ void run_step_mc(mt19937_state *state, double *spin, double *new_spin,
       delta_E += compute_deltaE_exchange_DMI_hexagonal(
           &spin[0], &new_spin[0], &ngbs[0], n_ngbs, J, D, i, new_i);
     } else {
-      delta_E += compute_deltaE_exchange_DMI(&spin[0], &new_spin[0], &ngbs[0], n_ngbs,
-                                             &nngbs[0], J, J1, D, D1, i, new_i);
+      delta_E += compute_deltaE_exchange_DMI(&spin[0], &new_spin[0], &ngbs[0], 
+                                             &nngbs[0], n_ngbs, J, J1, D, D1, i, new_i);
     }
 
     update = 0;

fidimag/common/dipolar/demag.c

@@ -5,7 +5,7 @@
 #include "demagcoef.h"
 
 
-double Nxxdipole(double x, double y, double z) {
+inline double Nxxdipole(double x, double y, double z) {
 	double x2 = x * x;
 	double y2 = y * y;
 	double z2 = z * z;
@@ -16,7 +16,7 @@ double Nxxdipole(double x, double y, double z) {
 	return -(2 * x2 - y2 - z2) / (R * R * r);
 }
 
-double Nxydipole(double x, double y, double z) {
+inline double Nxydipole(double x, double y, double z) {
 	double R = x * x + y * y + z * z;
 	if (R == 0)
 		return 0.0;
@@ -54,9 +54,16 @@ void compute_dipolar_tensors(fft_demag_plan *plan) {
 	int lenx = plan->lenx;
 	int leny = plan->leny;
 	int lenz = plan->lenz;
-    int lenxy = lenx * leny;
-
-	
+        int lenxy = lenx * leny;
+	// Parallelising this like this
+	// means that z should be the largest index
+	// in order to get better performance from threading.
+	// The data writing is not very clever here; we're 
+        // going to be invalidating the cache a lot. It would be better
+        // maybe to split this into six seperate loops for each component
+        // of the tensor in order that each thread is working on a smaller
+        // memory address range
+	#pragma omp parallel for private(j, i, x, y, z, id) schedule(dynamic, 32)
 	for (k = 0; k < lenz; k++) {
 		for (j = 0; j < leny; j++) {
 			for (i = 0; i < lenx; i++) {
@@ -81,7 +88,7 @@ void compute_dipolar_tensors(fft_demag_plan *plan) {
 void compute_demag_tensors(fft_demag_plan *plan) {
 
 	int i, j, k, id;
-	double x, y, z;
+	double x, y, z, radius_sq;
 
 	int nx = plan->nx;
 	int ny = plan->ny;
@@ -97,7 +104,7 @@ void compute_demag_tensors(fft_demag_plan *plan) {
 
 	double length = pow(dx*dy*dz, 1/3.0);
 	double asymptotic_radius_sq = pow(26.0*length,2.0);
-
+	#pragma omp parallel for private(j, i, id, x, y, z, radius_sq) schedule(dynamic, 32)
 	for (k = 0; k < lenz; k++) {
 		for (j = 0; j < leny; j++) {
 			for (i = 0; i < lenx; i++) {
@@ -107,7 +114,7 @@ void compute_demag_tensors(fft_demag_plan *plan) {
 				y = (j - ny + 1) * dy;
 				z = (k - nz + 1) * dz;
 
-				double radius_sq = x*x+y*y+z*z;
+				radius_sq = x*x+y*y+z*z;
 
 				if (radius_sq>asymptotic_radius_sq){
 					//printf("%g %g %g %g %g %g\n",x,y,z,dx,dy,dz);
@@ -338,7 +345,7 @@ void compute_fields(fft_demag_plan *restrict plan, double *restrict spin, double
 	//print_r("hz", plan->hz, plan->total_length);
 
 	double scale = -1.0  / plan->total_length;
-
+        #pragma omp parallel for private(j, i, id1, id2) schedule(dynamic, 32)
 	for (k = 0; k < nz; k++) {
 		for (j = 0; j < ny; j++) {
 			for (i = 0; i < nx; i++) {

fidimag/micro/__init__.py

@@ -1,8 +1,10 @@
 from .sim import Sim
 from .demag import Demag
 from .exchange import UniformExchange
+from .exchange_rkky import ExchangeRKKY
 from .zeeman import Zeeman, TimeZeeman
 from .anisotropy import UniaxialAnisotropy
 from .dmi import DMI
 from .baryakhtar import LLBar, LLBarFull
 from .simple_demag import SimpleDemag
+

fidimag/micro/dmi.py

@@ -52,7 +52,8 @@ class DMI(Energy):
 
     D       :: DMI vector norm which can be specified as an int, float, (X * n)
                array (X=6 for bulk DMI and X=4 for interfacial DMI), (n) array
-               or spatially dependent scalar field function.
+               or spatially dependent scalar field function. The units are 
+               Joules / ( meter **2 ).
 
                int, float: D will have the same magnitude for every NN of the
                spins at every mesh node, given by this magnitude

fidimag/micro/exchange.py

@@ -6,7 +6,15 @@
 class UniformExchange(Energy):
 
     """
+    UniformExchange(A, name='UniformExchange')
+
     Compute the exchange field in micromagnetics.
+    
+    Inputs:
+        A: float
+            A is the exchange stiffness constant measured in 
+            Joules / Meter (J / M)
+    
     """
 
     def __init__(self, A, name='UniformExchange'):

fidimag/micro/exchange_rkky.py

@@ -0,0 +1,61 @@
+import fidimag.extensions.micro_clib as micro_clib
+from .energy import Energy
+#from constant import mu_0
+
+
+class ExchangeRKKY(Energy):
+
+    """
+    RKKYExchange(sigma, Delta, z_bottom=0, z_top = -1, name='RKKYExchange')
+
+    Compute the RKKY-style exchange interaction defined by
+
+    E = sigma/Delta *(1-m_t * m_b)
+
+    where E is the energy density, sigma is the surface exchange coefficient between the two surfaces, 
+    Delata is the space thickness. m_t and m_b are the unit vectors of the top and bottom layers, respectively.
+    
+    Inputs:
+        sigma: float
+            sigma is the surface exchange stiffness constant.
+        Delta: float
+            Delta is the the space thickness in Meter.
+        z_bottom: int
+            z_bottom is the index of the bottom layer
+        z_top: int
+            z_top is the index of the top layer.
+    
+    """
+
+    def __init__(self, sigma, Delta=1e-9, z_bottom=0, z_top = -1, name='RKKYExchange'):
+        self.sigma = sigma/Delta
+        self.name = name
+        self.z_bottom = z_bottom
+        self.z_top = z_top
+        self.jac = True
+
+    def setup(self, mesh, spin, Ms):
+        super(ExchangeRKKY, self).setup(mesh, spin, Ms)
+        if self.z_top<0:
+            self.z_top+= self.nz
+ 
+
+    def compute_field(self, t=0, spin=None):
+        if spin is not None:
+            m = spin
+        else:
+            m = self.spin
+
+        micro_clib.compute_exchange_field_micro_rkky(m,
+                                                self.field,
+                                                self.energy,
+                                                self.Ms_inv,
+                                                self.sigma,
+                                                self.nx,
+                                                self.ny,
+                                                self.nz,
+                                                self.z_bottom,
+                                                self.z_top
+                                                )
+
+        return self.field