From ba65b4778745ffe6b75c9dbeca3ff0fc2501884d Mon Sep 17 00:00:00 2001 From: Alberto Ramos Date: Wed, 6 Oct 2021 18:27:43 +0200 Subject: [PATCH 1/3] Bug corrected --- src/LatticeGPU.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl index af2aa78..337e5bb 100644 --- a/src/LatticeGPU.jl +++ b/src/LatticeGPU.jl @@ -25,10 +25,10 @@ using .Space export SpaceParm export up, dw, updw, global_point -include("Fields/Fields.kl") +include("Fields/Fields.jl") +using .Fields export vector_field, scalar_field, nscalar_field - include("YM/YM.jl") using .YM From 18bb584fdfb110549558b0a7377c22c4fa014e93 Mon Sep 17 00:00:00 2001 From: Alberto Ramos Date: Wed, 6 Oct 2021 18:30:43 +0200 Subject: [PATCH 2/3] Bugs corrected --- src/Fields/Fields.jl | 12 +++++++++++- src/YM/YM.jl | 3 ++- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/src/Fields/Fields.jl b/src/Fields/Fields.jl index ae665f3..f5a5f4a 100644 --- a/src/Fields/Fields.jl +++ b/src/Fields/Fields.jl @@ -9,7 +9,17 @@ ### created: Wed Oct 6 17:37:03 2021 ### +module Fields + +using CUDA +using ..Space + vector_field(::Type{T}, lp::SpaceParm) where {T} = CuArray{T, 3}(undef, lp.bsz, lp.ndim, lp.rsz) -scalar_field(::Type{T}, lp::SpaceParm) where {T} = CuArray{T, 3}(undef, lp.bsz, lp.rsz) +scalar_field(::Type{T}, lp::SpaceParm) where {T} = CuArray{T, 2}(undef, lp.bsz, lp.rsz) nscalar_field(::Type{T}, n, lp::SpaceParm) where {T} = CuArray{T, 3}(undef, lp.bsz, n, lp.rsz) + +export vector_field, scalar_field, nscalar_field + +end + diff --git a/src/YM/YM.jl b/src/YM/YM.jl index ebdf41d..f257b70 100644 --- a/src/YM/YM.jl +++ b/src/YM/YM.jl @@ -15,6 +15,7 @@ module YM using CUDA, Random, StructArrays using ..Space using ..Groups +using ..Fields import Base.show @@ -85,7 +86,7 @@ end include("YMfields.jl") -export field, field_pln, randomize!, zero!, norm2 +export randomize!, zero!, norm2 include("YMact.jl") export krnl_plaq!, force0_wilson! From 0ae90b253437a7bf427f482eed9698cbc358986d Mon Sep 17 00:00:00 2001 From: Alberto Ramos Date: Sat, 9 Oct 2021 00:40:04 +0200 Subject: [PATCH 3/3] Added module for integrators --- src/LatticeGPU.jl | 5 ++++ src/MD/MD.jl | 35 ++++++++++++++++++++++++++++ src/YM/YM.jl | 1 + src/YM/YMhmc.jl | 59 ++++++++++++++++++++--------------------------- 4 files changed, 66 insertions(+), 34 deletions(-) create mode 100644 src/MD/MD.jl diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl index 337e5bb..c920257 100644 --- a/src/LatticeGPU.jl +++ b/src/LatticeGPU.jl @@ -29,6 +29,11 @@ include("Fields/Fields.jl") using .Fields export vector_field, scalar_field, nscalar_field +include("MD/MD.jl") +using .MD +export IntrScheme +export omf4, leapfrog + include("YM/YM.jl") using .YM diff --git a/src/MD/MD.jl b/src/MD/MD.jl new file mode 100644 index 0000000..a67d43b --- /dev/null +++ b/src/MD/MD.jl @@ -0,0 +1,35 @@ +### +### "THE BEER-WARE LICENSE": +### Alberto Ramos wrote this file. As long as you retain this +### notice you can do whatever you want with this stuff. If we meet some +### day, and you think this stuff is worth it, you can buy me a beer in +### return. +### +### file: MD.jl +### created: Fri Oct 8 21:53:14 2021 +### + +module MD + +# Dalla Brida / Luscher coefficients of +# fourth order integrator +const r1 = 0.08398315262876693 +const r2 = 0.25397851084105950 +const r3 = 0.68223653357190910 +const r4 = -0.03230286765269967 +const r5 = 0.5-r1-r3 +const r6 = 1.0-2.0*(r2+r4) + +struct IntrScheme{N, T} + r::NTuple{N, T} + eps::T + ns::Int64 +end + +omf4(::Type{T}, eps, ns) where T = IntrScheme{6,T}((r1,r2,r3,r4,r5,r6), eps, ns) +leapfrog(::Type{T}, eps, ns) where T = IntrScheme{2,T}((0.5,1.0,0.5), eps, ns) + +export IntrScheme, omf4, leapfrog + + +end diff --git a/src/YM/YM.jl b/src/YM/YM.jl index f257b70..4fecea5 100644 --- a/src/YM/YM.jl +++ b/src/YM/YM.jl @@ -16,6 +16,7 @@ using CUDA, Random, StructArrays using ..Space using ..Groups using ..Fields +using ..MD import Base.show diff --git a/src/YM/YMhmc.jl b/src/YM/YMhmc.jl index 55bfbbc..07579e2 100644 --- a/src/YM/YMhmc.jl +++ b/src/YM/YMhmc.jl @@ -48,14 +48,15 @@ function hamiltonian(mom, U, lp, gp, ymws) return K+V end -function HMC!(U, eps, ns, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace; noacc=false) +function HMC!(U, eps, ns, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}; noacc=false) where T + int = omf4(T, eps, ns) ymws.U1 .= U randomize!(ymws.mom, lp, ymws) hini = hamiltonian(ymws.mom, U, lp, gp, ymws) - OMF4!(ymws.mom, U, eps, ns, lp, gp, ymws) + MD!(ymws.mom, U, int, lp, gp, ymws) dh = hamiltonian(ymws.mom, U, lp, gp, ymws) - hini pacc = exp(-dh) @@ -76,40 +77,30 @@ function HMC!(U, eps, ns, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace; noacc return dh, acc end -function OMF4!(mom, U, eps, ns, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T}) where T <: AbstractFloat +function MD!(mom, U, int::IntrScheme{NI, T}, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T}) where {NI, T <: AbstractFloat} - r1::T = 0.08398315262876693 - r2::T = 0.2539785108410595 - r3::T = 0.6822365335719091 - r4::T = -0.03230286765269967 - r5::T = 0.5-r1-r3 - r6::T = 1.0-2.0*(r2+r4) - - ee = eps*gp.beta/gp.ng + ee = int.eps*gp.beta/gp.ng force_gauge(ymws, U, gp.c0, lp) - for i in 1:ns - mom .= mom .+ (r1*ee) .* ymws.frc1 - U .= expm.(U, mom, eps*r2) - - force_gauge(ymws, U, gp.c0, lp) - mom .= mom .+ (r3*ee) .* ymws.frc1 - U .= expm.(U, mom, eps*r4) - - force_gauge(ymws, U, gp.c0, lp) - mom .= mom .+ (r5*ee) .* ymws.frc1 - U .= expm.(U, mom, eps*r6) - - force_gauge(ymws, U, gp.c0, lp) - mom .= mom .+ (r5*ee) .* ymws.frc1 - U .= expm.(U, mom, eps*r4) - - force_gauge(ymws, U, gp.c0, lp) - mom .= mom .+ (r3*ee) .* ymws.frc1 - U .= expm.(U, mom, eps*r2) - - force_gauge(ymws, U, gp.c0, lp) - mom .= mom .+ (r1*ee) .* ymws.frc1 + mom .= mom .+ (int.r[1]*ee) .* ymws.frc1 + for i in 1:int.ns + k = 2 + off = 1 + for j in 1:NI-1 + U .= expm.(U, mom, int.eps*int.r[k]) + if k == NI + off = -1 + end + k += off + + force_gauge(ymws, U, gp.c0, lp) + if (i < int.ns) && (k == 1) + mom .= mom .+ (2*int.r[k]*ee) .* ymws.frc1 + else + mom .= mom .+ (int.r[k]*ee) .* ymws.frc1 + end + k += off + end end - + return nothing end