Miguel escobar

Project.toml

@@ -2,3 +2,11 @@ name = "MyPkg"
 uuid = "b59f1ed4-4c8e-4fc7-a035-90d7cbd7aaf3"
 authors = ["Carlos Castillo Passi <[email protected]>"]
 version = "0.1.0"
+
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
+[compat]
+LinearAlgebra = "1.10.0"
+StaticArrays = "1.9.10"

benchmarks/benchmark.jl

@@ -11,13 +11,13 @@ tmax = 3.0
 
 # Define the benchmark
 suite = BenchmarkGroup()
-suite["function1"][github_username] = @benchmarkable solve(m0, Δt, tmax, ForwardEuler())
-# suite["function2"][github_username] = @benchmarkable solve(m0, Δt, t_max, params, ForwardEuler())
+#suite["function1"][github_username] = @benchmarkable solve(m0, Δt, tmax, ForwardEuler())
+suite["function2"][github_username] = @benchmarkable solve(m0, Δt, tmax, ForwardEuler())
 
 # Tune and run the benchmark
 tune!(suite)
 results = run(suite, verbose=true)
 median_results = median(results)
-median_results["function1"][github_username].gctime = std(results["function1"][github_username]).time
-# median_results["function2"][github_username].gctime = std(results["function2"][github_username]).time
+#median_results["function1"][github_username].gctime = std(results["function1"][github_username]).time
+median_results["function2"][github_username].gctime = std(results["function2"][github_username]).time
 BenchmarkTools.save("benchmarks/benchmark_results.json", median_results)

benchmarks/benchmark_results.json

@@ -1 +0,0 @@
-

src/MyPkg.jl

@@ -1,5 +1,69 @@
 module MyPkg
+using LinearAlgebra: cross
+using StaticArrays
 
-# Your code here
+export solve, step, Theoretical, ForwardEuler
 
-end  # module MyPkg
+const gamma::Float32 = 2pi * 42.58 * 1e6
+const M0::Float32 = 1.0
+const T1::Float32 = 1.0
+const T2::Float32 = 0.5
+const Bz::Float32 = 1e-7
+const gammaBz::Float32 = gamma * Bz
+
+const M0T1::Float32 = -M0/T1
+
+const Ti = SA{Float32}[1/T2, 1/T2, 1/T1]
+const m0t1 = SA{Float32}[0.0,0.0, M0/T1]
+const B = SA{Float32}[0.0, 0.0, Bz]
+const tmax::Float32 = 3.0
+const m0 = SA{Float32}[M0, 0.0, 0.0]
+
+struct ForwardEuler
+
+end
+
+struct Theoretical
+
+end
+
+function step(dt, m, gammaBzdt, M0T1dt, Tidt, method::ForwardEuler)
+    return m .+ bloch(m, gammaBzdt, M0T1dt, Tidt)
+end
+
+function aux(a, gammaBzdt, M0T1dt)
+    return @SVector [a[2] * gammaBzdt, -a[1] * gammaBzdt, M0T1dt]
+  end
+
+function bloch(m, gammaBzdt, M0T1dt, Tidt)
+    return aux(m, gammaBzdt, M0T1dt) .- (Tidt .* m) 
+end
+
+function zeros_via_calloc(::Type{T}, dims::Integer...) where T
+    ptr = Ptr{T}(Libc.calloc(prod(dims), sizeof(T)))
+    return unsafe_wrap(Array{T}, ptr, dims; own=true)
+ end
+
+
+function solve(m0, dt, tmax, method)
+    Nsteps = Int(ceil(tmax/dt))
+    m = SVector{3, Float32}(m0)
+    mt = zeros_via_calloc(Float32, ceil(Int64,tmax/dt) + 1, 3)#zeros(Float64, (ceil(Int64,tmax/dt) + 1, 3))
+    gammaBzdt, M0T1dt, Tidt = gammaBz * dt, M0T1 * dt, Ti .* dt
+    @inbounds for i in 1:Nsteps
+        m = step(dt, m,gammaBzdt, M0T1dt, Tidt, method)
+        mt[i, :] .= m
+    end
+    return mt
+end
+
+function solve(M0, dt, tmax, method::Theoretical)
+    mt = zeros(3, Int(ceil(tmax/dt)) + 1)
+    t = 0:dt:tmax
+    mt[1, :] .= M0 .* cos.(gamma .* Bz .* t) .* exp.(-t ./ T2)
+    mt[2, :] .= -M0 .* sin.(gamma .* Bz .* t) .* exp.(-t ./ T2)
+    mt[3, :] .= M0 .* (1 .- exp.(-t ./ T1))
+    return mt
+end
+
+end

test/Project.toml

@@ -1,2 +1,3 @@
 [deps]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/runtests.jl

@@ -1,12 +1,19 @@
 using Test
 using MyPkg
+using StaticArrays
 
-# Inputs
-m0   = [1.0, 0.0, 0.0]
-Δt   = 0.001
-tmax = 3.0
+const gamma = 2pi * 42.58 * 1e6
+const M0 = 1.0
+const T1 = 1.0
+const T2 = 0.5
+const Bz = 1e-7
 
-expected_result    = solve(m0, Δt, tmax, Theoretical())
-numerical_solution = solve(m0, Δt, tmax, ForwardEuler())
+const Ti = SA[1/T2, 1/T2, 1/T1]
+const m0t1 = SA[0.0,0.0, M0/T1]
+const B = SA[0.0, 0.0, Bz]
+const tmax = 3.0
+const m0 = SA[M0, 0.0, 0.0]
 
-@test numerical_solution ≈ expected_result
+dt = 1e-3
+
+@test abs(solve(m0, dt, 3.0, ForwardEuler())[1, end] - solve(M0, dt, 3.0, Theoretical())[end, 1]) <= 2e-3