Merge pull request #135 from computationalmodelling/demag-2d

Demag 2d

Author: davidcortesortuno

fidimag/atomistic/hexagonal_mesh.py

@@ -258,12 +258,14 @@ def init_grid(self):
         vertex_counter = 0
         vertices = []  # list of tuples of coordinates
         hexagons = []  # list of tuples of vertices
+        self.corners = []
         for j in range(self.ny):
             for i in range(self.nx):
                 index = self._index(i, j)
                 x, y = self.coordinates[index][0], self.coordinates[index][1]
                 # self.radius is the inradius while self.h/2  is the circumradius
                 corners = self.hexagon_corners(x, y, self.h * 0.5)
+                self.corners.append(corners)
                 hexagon = []
                 # We'll go through the corners in a counter-clockwise direction.
                 # For each corner, we think about if it's a "new" vertex, or
@@ -328,6 +330,8 @@ def init_grid(self):
                     hexagon.append(vertex_counter)
                     vertex_counter += 1
                 hexagons.append(hexagon)
+
+        self.corners = np.array(self.corners)
         return np.array(vertices), np.array(hexagons)
 
     def index(self, i, j):

fidimag/micro/demag.py

@@ -50,26 +50,25 @@ def setup(self, mesh, spin, Ms, Ms_inv):
                 asymptotic_radius = options['asymptotic_radius']
             if 'dipolar_radius' in options:
                 dipolar_radius = options['dipolar_radius']
-            if 'tensor_file_name' in options:
-                tensor_file_name = options['tensor_file_name']
-
-            if tensor_file_name:
-                if not (os.path.exists(tensor_file_name+'.npz')):
-                    self.demag.compute_tensors_2dpbc(tensors, pbc_2d_error, sample_repeat_nx, sample_repeat_ny, dipolar_radius)
-                else:
-                    npzfile = np.load(tensor_file_name+'.npz')
-                    geo_info = npzfile['geo']
-                    geo_arr = np.array([self.nx,self.ny, self.nz, self.dx, self.dy, self.dz], dtype=np.float)
-                    #print 'info', geo_info
-                    if not np.allclose(geo_arr, geo_info):
-                        self.demag.compute_tensors_2dpbc(tensors, pbc_2d_error, sample_repeat_nx, sample_repeat_ny, dipolar_radius)
-                    else:
-                        tensors = npzfile['tensors']
-                geo_arr = np.array([self.nx, self.ny, self.nz, self.dx, self.dy, self.dz], dtype=np.float)
-                np.savez(tensor_file_name+'.npz', geo=geo_arr, tensors=tensors)
-
-            else:
-                self.demag.compute_tensors_2dpbc(tensors, pbc_2d_error, sample_repeat_nx, sample_repeat_ny, dipolar_radius)
+            #if 'tensor_file_name' in options:
+            #    tensor_file_name = options['tensor_file_name']
+
+            #if tensor_file_name:
+                #if not (os.path.exists(tensor_file_name+'.npz')):
+                #    self.demag.compute_tensors_2dpbc(tensors, pbc_2d_error, sample_repeat_nx, sample_repeat_ny, dipolar_radius)
+                #else:
+                #    npzfile = np.load(tensor_file_name+'.npz')
+                #    geo_info = npzfile['geo']
+                #    geo_arr = np.array([self.nx,self.ny, self.nz, self.dx, self.dy, self.dz], dtype=np.float)
+                #    #print 'info', geo_info
+                #    if not np.allclose(geo_arr, geo_info):
+                #        self.demag.compute_tensors_2dpbc(tensors, pbc_2d_error, sample_repeat_nx, sample_repeat_ny, dipolar_radius)
+                #    else:
+                #        tensors = npzfile['tensors']
+                #geo_arr = np.array([self.nx, self.ny, self.nz, self.dx, self.dy, self.dz], dtype=np.float)
+                #np.savez(tensor_file_name+'.npz', geo=geo_arr, tensors=tensors)
+
+            self.demag.compute_tensors_2dpbc(tensors, pbc_2d_error, sample_repeat_nx, sample_repeat_ny, dipolar_radius)
 
             #print tensors
             self.demag.fill_demag_tensors(tensors)

tests/test_demag_libraries.py

@@ -4,7 +4,7 @@
 we only use cuboid meshes, so we only have the FFT approach
 for the dipolar field)
 """
-
+import fidimag
 from fidimag.atomistic import Sim
 from fidimag.common import CuboidMesh
 from fidimag.atomistic.hexagonal_mesh import HexagonalMesh
@@ -222,6 +222,35 @@ def test_hexagonal_demags_2D():
     assert (DemagFull_energy - demag_2fft_energy) < 1e-10
 
 
+
+def test_demag_2d_pbc():
+    """
+    Attempt to check that demag with 2d_pbc option does
+    not give a nonsensical answer.
+    """
+    A=1.3e-11
+    Ms=8.6e5
+    n = 40
+    d = 2.5
+    
+    mesh = fidimag.common.CuboidMesh(nx=n, ny=n, nz=1, dx=d, dy=d, dz=d, unit_length=1e-9, periodicity=(True, True, False))
+    sim = fidimag.micro.Sim(mesh, name="pbc_2d_bug")
+    sim.set_Ms(Ms)
+    sim.set_m((0, 0, 1.0), normalise=True)
+    demag = fidimag.micro.Demag(pbc_2d=True)
+
+    sim.add(demag)
+    sim.compute_effective_field(0)
+
+    assert not np.isnan(demag.demag.tensor_xx).any()
+    assert not np.isnan(demag.demag.tensor_xy).any()
+    assert not np.isnan(demag.demag.tensor_xz).any()
+    assert not np.isnan(demag.demag.tensor_yy).any()
+    assert not np.isnan(demag.demag.tensor_yz).any()
+    assert not np.isnan(demag.demag.tensor_zz).any()
+    assert not np.isnan(sim.field).any(), "NaN in demag array"
+    assert 2 == 1
+
 if __name__ == '__main__':
     test_hexagonal_demags_1Dchain()
     test_cuboid_demags_1Dchain()