fix a bug for stt (i1*i2 may overflow)

src/llg/llg.jl

@@ -70,8 +70,9 @@ compute tau = (u.nabla) m.
 
     #x-direction
     i1::Int32 = ngbs[1, I] #we assume that i1<0 for the area with Ms=0
-    i2::Int32 = ngbs[2, I]
-    factor::T = i1 * i2 > 0 ? 1 / (2 * dx) : 1 / dx
+    i2::Int32 = ngbs[2, I] 
+    # i1 * i2 may overflow
+    factor::T = (i1 > 0 && i2 > 0) ? 1 / (2 * dx) : 1 / dx
     i1 < 0 && (i1 = I)
     i2 < 0 && (i2 = I)
     j1 = 3 * i1 - 2
@@ -84,7 +85,7 @@ compute tau = (u.nabla) m.
     #y-direction
     i1 = ngbs[3, I]
     i2 = ngbs[4, I]
-    factor = i1 * i2 > 0 ? 1 / (2 * dy) : 1 / dy
+    factor = (i1 > 0 && i2 > 0) ? 1 / (2 * dy) : 1 / dy
     i1 < 0 && (i1 = I)
     i2 < 0 && (i2 = I)
     j1 = 3 * i1 - 2
@@ -97,7 +98,7 @@ compute tau = (u.nabla) m.
     #z-direction
     i1 = ngbs[5, I]
     i2 = ngbs[6, I]
-    factor = i1 * i2 > 0 ? 1 / (2 * dz) : 1 / dz
+    factor = (i1 > 0 && i2 > 0) ? 1 / (2 * dz) : 1 / dz
     i1 < 0 && (i1 = I)
     i2 < 0 && (i2 = I)
     j1 = 3 * i1 - 2

src/micro/add_field.jl

@@ -249,7 +249,7 @@ function add_dmi(sim::MicroSim, D::NumberOrTupleOrArrayOrFunction; name="dmi", t
             error("D2d only support uniform DMI!")
         end
     else
-        error("Supported DMI type:", "interfacial", "bulk", "D2d")
+        error("Supported DMI type:", " interfacial", " bulk", " D2d")
     end
 
     push!(sim.interactions, dmi)

src/vtk.jl

@@ -38,10 +38,11 @@ function save_vtk(sim::AbstractSim, fname::String; fields::Array{String,1}=Strin
     nx, ny, nz = mesh.nx, mesh.ny, mesh.nz
     xyz = zeros(Float32, 3, nx + 1, ny + 1, nz + 1)
     dx, dy, dz = mesh.dx, mesh.dy, mesh.dz
+    scale_factor = 10^floor(log10(dx))
     for k in 1:(nz + 1), j in 1:(ny + 1), i in 1:(nx + 1)
-        xyz[1, i, j, k] = (i - 0.5 - nx / 2) * dx
-        xyz[2, i, j, k] = (j - 0.5 - ny / 2) * dy
-        xyz[3, i, j, k] = (k - 0.5 - nz / 2) * dz
+        xyz[1, i, j, k] = (i - 0.5 - nx / 2) * dx / scale_factor
+        xyz[2, i, j, k] = (j - 0.5 - ny / 2) * dy / scale_factor
+        xyz[3, i, j, k] = (k - 0.5 - nz / 2) * dz / scale_factor
     end
     vtk = vtk_grid(fname, xyz)
     spin = zeros(eltype(sim.spin), length(sim.spin))
@@ -61,6 +62,9 @@ function save_vtk(sim::AbstractSim, fname::String; fields::Array{String,1}=Strin
             end
         end
     end
+    if scale_factor != 1.0
+        vtk["scale_factor", VTKFieldData()] = string(scale_factor)
+    end
     return vtk_save(vtk)
 end
 
@@ -69,10 +73,11 @@ function save_vtk_points(sim::AbstractSim, fname::String; fields::Array{String,1
     nx, ny, nz = mesh.nx, mesh.ny, mesh.nz
     xyz = zeros(Float32, 3, nx, ny, nz)
     dx, dy, dz = mesh.dx, mesh.dy, mesh.dz
+    scale_factor = 10^floor(log10(dx))
     for k in 1:nz, j in 1:ny, i in 1:nx
-        xyz[1, i, j, k] = (i - 0.5 - nx / 2) * dx
-        xyz[2, i, j, k] = (j - 0.5 - ny / 2) * dy
-        xyz[3, i, j, k] = (k - 0.5 - nz / 2) * dz
+        xyz[1, i, j, k] = (i - 0.5 - nx / 2) * dx / scale_factor
+        xyz[2, i, j, k] = (j - 0.5 - ny / 2) * dy / scale_factor
+        xyz[3, i, j, k] = (k - 0.5 - nz / 2) * dz / scale_factor
     end
     vtk = vtk_grid(fname, xyz)
 
@@ -92,6 +97,9 @@ function save_vtk_points(sim::AbstractSim, fname::String; fields::Array{String,1
             end
         end
     end
+    if scale_factor != 1.0
+        vtk["scale_factor", VTKFieldData()] = string(scale_factor)
+    end
     return vtk_save(vtk)
 end
 

test/test_interlayer.jl

@@ -15,6 +15,22 @@ function m0_fun(i, j, k, dx, dy, dz)
     return (1, 1, 1)
 end
 
+function init_single_skx(i, j, k, dx, dy, dz)
+    r = 20
+    if k == 3
+        if (i - 10)^2 + (j - 30)^2 < r^2
+            return (0, -0.01, -1)
+        end
+        return (0, 0.01, 1)
+    end
+    if k == 1
+        if (i - 10)^2 + (j - 30)^2 < r^2
+            return (0, 0.01, 1)
+        end
+        return (0, -0.01, -1)
+    end
+end
+
 function spatial_Ms(i, j, k, dx, dy, dz)
     if k == 1 || k == 3
         return 8e5
@@ -38,7 +54,8 @@ function test_interlayer()
     MicroMagnetic.effective_field(sim, sim.spin)
 
     expected_exch = [4, 5, 6, 0, 0, 0, 1, 2, 3, 0, 0, 0] .* J / (mu_0 * Ms * dz)
-    expected_exch = vcat([repeat(expected_exch[i:i+2], 15) for i in 1:3:length(expected_exch)]...)
+    expected_exch = vcat([repeat(expected_exch[i:(i + 2)], 15)
+                          for i in 1:3:length(expected_exch)]...)
 
     fx1 = 1 / (mu_0 * Ms * dz) * MicroMagnetic.cross_x(Dx, Dy, Dz, 4.0, 5.0, 6.0)
     fy1 = 1 / (mu_0 * Ms * dz) * MicroMagnetic.cross_y(Dx, Dy, Dz, 4.0, 5.0, 6.0)
@@ -48,11 +65,39 @@ function test_interlayer()
     fz2 = 1 / (mu_0 * Ms * dz) * MicroMagnetic.cross_z(1.0, 2.0, 3.0, Dx, Dy, Dz)
 
     expected_dmi = [fx1, fy1, fz1, 0, 0, 0, fx2, fy2, fz2, 0, 0, 0]
-    expected_dmi = vcat([repeat(expected_dmi[i:i+2], 15) for i in 1:3:length(expected_dmi)]...)
+    expected_dmi = vcat([repeat(expected_dmi[i:(i + 2)], 15)
+                         for i in 1:3:length(expected_dmi)]...)
 
     @test isapprox(expected_exch, Array(exch.field))
     @test isapprox(expected_dmi, Array(dmi.field))
 end
 
+function test_interlayer_exch()
+    J = 1e-5
+    Ms = 8e5
+    dz = 2e-9
+
+    nx, ny, nz = 30, 80, 3
+    mesh = FDMesh(; dx=2e-9, dy=2e-9, dz=2e-9, nx=nx, ny=ny, nz=nz, pbc="xy")
+    sim = Sim(mesh)
+    set_Ms(sim, spatial_Ms)
+    init_m0(sim, init_single_skx)
+
+    exch = add_exch_int(sim, J; k1=1, k2=3)
+
+    MicroMagnetic.effective_field(sim, sim.spin)
+
+    f = Array(exch.field)
+    f = reshape(f, (3, nx, ny, nz))
+    f1 = f[:, :, :, 1]
+    f3 = f[:, :, :, 3]
+
+    s = f1 .+ f3
+
+    @test maximum(s) == 0
+    @test minimum(s) == 0
+end
+
 @using_gpu()
-test_functions("Interlayer", test_interlayer)
+
+test_functions("Interlayer", test_interlayer, test_interlayer_exch)