diff --git a/docs/make.jl b/docs/make.jl index f52fe3e..890b3f1 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -13,6 +13,8 @@ makedocs(sitename="LatticeGPU", modules=[LatticeGPU], doctest=true, "Yang-Mills" => "ym.md", "Gradient flow" => "flow.md", "Schrödinger Functional" => "sf.md", + "Dirac" => "dirac.md", + "Solvers" => "solvers.md", "Input Output" => "io.md" ], repo = "https://igit.ific.uv.es/alramos/latticegpu.jl") diff --git a/docs/src/dirac.md b/docs/src/dirac.md new file mode 100644 index 0000000..75a330a --- /dev/null +++ b/docs/src/dirac.md @@ -0,0 +1,112 @@ + +# Dirac operator + +The module `Dirac` has the necessary stuctures and function +to simulate non-dynamical 4-dimensional Wilson fermions. + +There are two main data structures in this module, the structure [`DiracParam`](@ref) + +```@docs +DiracParam +``` + +and the workspace [`DiracWorkspace`](@ref) + +```@docs +DiracWorkspace +``` + +The workspace stores four fermion fields, namely `.sr`, `.sp`, `.sAp` and `.st`, used +for different purposes. If the representation is either `SU2fund` of `SU3fund`, an extra +field with values in `U2alg`/`U3alg` is created to store the clover, used for the improvent. + +## Functions + +The functions [`Dw!`](@ref), [`g5Dw!`](@ref) and [`DwdagDw!`](@ref) are all related to the +Wilson-Dirac operator. + +The action of the Dirac operator `Dw!` is the following: + +```math +D_w\psi (\vec{x} = x_1,x_2,x_3,x_4) = (4 + m_0 + i \mu \gamma_5)\psi(\vec{x}) - +``` +```math + - \frac{1}{2}\sum_{\mu = 1}^4 \theta (\mu) (1-\gamma_\mu) U_\mu(\vec{x}) \psi(\vec{x} + \hat{\mu}) + \theta^* (\mu) (1 + \gamma_\mu) U^{-1}_\mu(\vec{x} - \hat{\mu}) \psi(\vec{x} - \hat{\mu}) +``` + +where $$m_0$$ and $$\theta$$ are respectively the values `.m0` and `.th` of [`DiracParam`](@ref). +Note that $$|\theta(\mu)|=1$$ is not built into the code, so it should be imposed explicitly. + +Additionally, if |`dpar.csw`| > 1.0E-10, the clover term is assumed to be stored in `ymws.csw`, which +can be done via the [`Csw!`](@ref) function. In this case we have the Sheikholeslami–Wohlert (SW) term +in `Dw!`: + +```math + \delta D_w^{sw} = \frac{i}{2}c_{sw} \sum_{\pi = 1}^6 F^{cl}_\pi \sigma_\pi \psi(\vec{x}) +``` +where the $$\sigma$$ matrices are those described in the `Spinors` module and the index $$\pi$$ runs +as specified in `lp.plidx`. + +If the boudary conditions, defined in `lp`, are either `BC_SF_ORBI,D` or `BC_SF_AFWB`, the +improvement term + +```math + \delta D_w^{SF} = (c_t -1) (\delta_{x_4,a} \psi(\vec{x}) + \delta_{x_4,T-a} \psi(\vec{x})) +``` +is added. Since the time-slice $$t=T$$ is not stored, this accounts to modifying the second +and last time-slice. + +Note that the Dirac operator for SF boundary conditions assumes that the value of the field +in the first time-slice is zero. To enforce this, we have the function + +```@docs +SF_bndfix! +``` + +Note that this is not enforced in the Dirac operators, so if the field `so` does not satisfy SF +boundary conditions, it will not (in general) satisfy them after applying [`Dw!`](@ref) +or [`g5Dw!`](@ref). This function is called for the function [`DwdagDw!`](@ref), so in this case +`so` will always be a proper SF field after calling this function. + +The function [`Csw!`](@ref) is used to store the clover in `dws.csw`. It is computed +according to the expression + +```math +F_{\mu,\nu} = \frac{1}{8} (Q_{\mu \nu} - Q_{\nu \mu}) +``` + +where +```math +Q_{\mu\nu} = U_\mu(\vec{x})U_{\nu}(x+\mu)U_{\mu}^{-1}(\vec{x}+\nu)U_{\nu}(\vec{x}) + +U_{\nu}^{-1}(\vec{x}-\nu) U_\mu (\vec{x}-\nu) U_{\nu}(\vec{x} +\mu - \nu) U^{-1}_{\mu}(\vec{x}) + +``` +```math ++ U^{-1}_{\mu}(x-\mu)U_\nu^{-1}(\vec{x} - \mu - \nu)U_\mu(\vec{x} - \mu - \nu)U_\nu^{-1}(x-\nu) + +``` +```math ++U_{\nu}(\vec{x})U_{\mu}^{-1}(\vec{x} + \nu - \mu)U^{-1}_{\nu}(\vec{x} - \mu)U_\mu(\vec{x}-\mu) + +``` + +The correspondence between the tensor field and the GPU-Array is the following: +```math +F[b,1,r] \to F_{41}(b,r) ,\quad F[b,2,r] \to F_{42}(b,r) ,\quad F[b,3,r] \to F_{43}(b,r) +``` +```math +F[b,4,r] \to F_{31}(b,r) ,\quad F[b,5,r] \to F_{32}(b,r) ,\quad F[b,6,r] \to F_{21}(b,r) +``` +where $$(b,r)$$ labels the lattice points as explained in the module `Space` + +The function [`pfrandomize!`](@ref), userfull for stochastic sources, is also present. It +randomizes a fermion field either in all the space or in a specifit time-slice. + +The generic interface of these functions reads + +```@docs +Dw! +g5Dw! +DwdagDw! +Csw! +pfrandomize! +mtwmdpar +``` diff --git a/docs/src/solvers.md b/docs/src/solvers.md new file mode 100644 index 0000000..b6fc606 --- /dev/null +++ b/docs/src/solvers.md @@ -0,0 +1,86 @@ + +# Solvers + +The module `Solvers` contains the functions to invert the Dirac +operator as well as functions to obtain specific propagators. + + +## CG.jl + +The function [`CG!`](@ref) implements the Conjugate gradient +algorith for the operator A + +```@docs +CG! +``` + +where the tolerance is normalized with respect to $$|$$``si``$$|^2$$. +The operator A must have the same input structure as all the Dirac +operators. If the maximum number of iterations `maxiter` is reached, +the function will throw an error. The estimation for $$A^{-1}x$$ is +stored in ``si``, and the number of iterations is returned. + + +Note that all the fermion field in ``dws`` are used +inside the function and will be modified. In particular, the final residue +is given by $$|$$``dws.sr``$$|^2$$. + + + +## Propagators.jl + +In this file, we define a couple of useful functions to obtain certain +propagators. + +```@docs +propagator! +``` + +Internally, this function solves the equation + +```math + D_w^\dagger D_w \psi = \gamma_5 D_w \gamma_5 \eta +``` + +where $$\eta$$ is either a point-source with the specified color and spin +or a random source in a time-slice and stores the value in ``pro``. +To solve this equation, the [`CG!`](@ref) function is used. + + +For the case of SF boundary conditions, we have the boundary-to-bulk +propagator, implemented by the function [`bndpropagator!`](@ref) + +```@docs +bndpropagator! +``` + +This propagator is defined by the equation: + +```math +D_W S(x) = \frac{c_t}{\sqrt{V}} \delta_{x_0,1} U_0^\dagger(0,\vec{x}) P_+ +``` + +The analog for the other boundary is implemented in the function [`Tbndpropagator!`](@ref) + +```@docs +Tbndpropagator! +``` + +defined by the equation: + +```math +D_W R(x) = \frac{c_t}{\sqrt{V}} \delta_{x_0,T-1} U_0(T-1,\vec{x}) P_- +``` + +Where $$P_\pm = (1 \pm \gamma_0)/2$$. The boundary-to-boundary +propagator + +```math + - \frac{c_t}{\sqrt{V}} \sum_\vec{x} U_0 ^\dagger (T-1,\vec{x}) P_+ S(T-1,\vec{x}) +``` + +is computed by the function [`bndtobnd`](@ref) + +```@docs +bndtobnd +``` diff --git a/docs/src/spinors.md b/docs/src/spinors.md new file mode 100644 index 0000000..16675bc --- /dev/null +++ b/docs/src/spinors.md @@ -0,0 +1,86 @@ +# Spinors + +The module Spinors defines the necessary functions for the structure `Spinor{NS,G}`, +which is a NS-tuple with values in G. + +The functions `norm`, `norm2`, `dot`, `*`, `/`, `/`, `+`, `-`, `imm` and `mimm`, +if defined for G, are extended to Spinor{NS,G} for general NS. + +For the 4d case where NS = 4 there are some specific functions to implement different +operations with the gamma matrices. The convention for these matrices is + + +```math +\gamma _4 = \left( + \begin{array}{cccc} + 0 & 0 & -1 & 0\\ + 0 & 0 & 0 & -1\\ + -1 & 0 & 0 & 0\\ + 0 & -1 & 0 & 0\\ + \end{array} + \right) + \quad + \gamma_1 = \left( + \begin{array}{cccc} + 0 & 0 & 0 & -i\\ + 0 & 0 & -i & 0\\ + 0 & i & 0 & 0\\ + i & 0 & 0 & 0\\ + \end{array} + \right) +``` +```math + \gamma _2 = \left( + \begin{array}{cccc} + 0 & 0 & 0 & -1\\ + 0 & 0 & 1 & 0\\ + 0 & 1 & 0 & 0\\ + -1 & 0 & 0 & 0\\ + \end{array} + \right) + \quad + \gamma_3 = \left( + \begin{array}{cccc} + 0 & 0 & -i & 0\\ + 0 & 0 & 0 & i\\ + i & 0 & 0 & 0\\ + 0 & -i & 0 & 0\\ + \end{array} + \right) +``` + + +The function [`dmul`](@ref) implements the multiplication over the $$\gamma$$ matrices + +```@docs +dmul +``` + +The function [`pmul`](@ref) implements the $$ (1 \pm \gamma_N) $$ proyectors. The functions +[`gpmul`](@ref) and [`gdagpmul`](@ref) do the same and then multiply each element by `g`and +g^-1 repectively. + +```@docs +pmul +gpmul +gdagpmul +``` + +## Some examples + +Here we just display some examples for these functions. We display it with `ComplexF64` +instead of `SU3fund` or `SU2fund` for simplicity. + + +```@setup exs +import Pkg # hide +Pkg.activate("/home/alberto/code/julia/LatticeGPU/") # hide +using LatticeGPU # hide +``` +```@repl exs +spin = Spinor{4,Complex{Float64}}((1.0,im*0.5,2.3,0.0)) +println(spin) +println(dmul(Gamma{4},spin)) +println(pmul(Pgamma{2,-1},spin)) + +``` diff --git a/src/Dirac/Dirac.jl b/src/Dirac/Dirac.jl index 2d78098..74eb63c 100644 --- a/src/Dirac/Dirac.jl +++ b/src/Dirac/Dirac.jl @@ -19,20 +19,31 @@ using ..Fields using ..YM using ..Spinors +""" + struct DiracParam{T,R} +Stores the parameters of the Dirac operator. It can be generated via the constructor `function DiracParam{T}(::Type{R},m0,csw,th,tm,ct)`. The first argument can be ommited and is taken to be `SU3fund`. +The parameters are: + +- `m0::T` : Mass of the fermion +- `csw::T` : Improvement coefficient for the Csw term +- `th{Ntuple{4,Complex{T}}}` : Phase for the fermions included in the boundary conditions, reabsorbed in the Dirac operator. +- `tm` : Twisted mass parameter +- `ct` : Boundary improvement term, only used for Schrödinger Funtional boundary conditions. +""" struct DiracParam{T,R} m0::T csw::T th::NTuple{4,Complex{T}} + tm::T ct::T - - function DiracParam{T}(::Type{R},m0,csw,th,ct) where {T,R} - return new{T,R}(m0,csw,th,ct) + function DiracParam{T}(::Type{R},m0,csw,th,tm,ct) where {T,R} + return new{T,R}(m0,csw,th,tm,ct) end - function DiracParam{T}(m0,csw,th,ct) where {T} - return new{T,SU3fund}(m0,csw,th,ct) + function DiracParam{T}(m0,csw,th,tm,ct) where {T} + return new{T,SU3fund}(m0,csw,th,tm,ct) end end function Base.show(io::IO, dpar::DiracParam{T,R}) where {T,R} @@ -40,11 +51,24 @@ function Base.show(io::IO, dpar::DiracParam{T,R}) where {T,R} println(io, "Wilson fermions in the: ", R, " representation") println(io, " - Bare mass: ", dpar.m0," // Kappa = ",0.5/(dpar.m0+4)) println(io, " - Csw : ", dpar.csw) - println(io, " - c_t: ", dpar.ct) println(io, " - Theta: ", dpar.th) + println(io, " - Twisted mass: ", dpar.tm) + println(io, " - c_t: ", dpar.ct) return nothing end + +""" + struct DiracWorkspace{T} + +Workspace needed to work with fermion fields. It contains four scalar fermion fields and, for the SU2fund and SU3fund, a U(N) field to store the clover term. + +It can be created with the constructor `DiracWorkspace(::Type{G}, ::Type{T}, lp::SpaceParm{4,6,B,D})`. For example: + +dws = DiracWorkspace(SU2fund,Float64,lp); +dws = DiracWorkspace(SU3fund,Float64,lp); + +""" struct DiracWorkspace{T} sr sp @@ -146,18 +170,26 @@ function krnl_csw!(csw::AbstractArray{T}, U, Ubnd, ipl, lp::SpaceParm{4,M,B,D}) return nothing end + +""" + function Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4,6,B,D}) + +Computes the Dirac operator (with the Wilson term) `\`\``D_w``\`\` with gauge field U and parameters `dpar` of the field `si` and stores it in `so`. +If `dpar.csw` is different from zero, the clover term should be stored in `dws.csw` via the Csw! function and is automatically included in the operator. + +""" function Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4,6,B,D}) where {B,D} if abs(dpar.csw) > 1.0E-10 @timeit "Dw" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dwimpr!(so, U, si, dws.csw, dpar.m0, dpar.th, dpar.csw, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dwimpr!(so, U, si, dws.csw, dpar.m0, dpar.tm, dpar.th, dpar.csw, lp) end end else @timeit "Dw" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw!(so, U, si, dpar.m0, dpar.th, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw!(so, U, si, dpar.m0, dpar.tm, dpar.th, lp) end end end @@ -165,7 +197,7 @@ function Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4,6 return nothing end -function krnl_Dwimpr!(so, U, si, Fcsw, m0, th, csw, lp::SpaceParm{4,6,B,D}) where {B,D} +function krnl_Dwimpr!(so, U, si, Fcsw, m0, tm, th, csw, lp::SpaceParm{4,6,B,D}) where {B,D} b = Int64(CUDA.threadIdx().x); r = Int64(CUDA.blockIdx().x) @@ -180,8 +212,8 @@ function krnl_Dwimpr!(so, U, si, Fcsw, m0, th, csw, lp::SpaceParm{4,6,B,D}) wher @inbounds begin - so[b,r] = (4+m0)*si[b,r] + 0.5*csw*im*( Fcsw[b,1,r]*dmul(Gamma{10},si[b,r]) + Fcsw[b,2,r]*dmul(Gamma{11},si[b,r]) + Fcsw[b,3,r]*dmul(Gamma{12},si[b,r]) - +Fcsw[b,4,r]*dmul(Gamma{15},si[b,r]) + Fcsw[b,5,r]*dmul(Gamma{14},si[b,r]) + Fcsw[b,6,r]*dmul(Gamma{13},si[b,r])) + so[b,r] = (4+m0)*si[b,r]+ im*tm*dmul(Gamma{5},si[b,r]) + 0.5*csw*im*( Fcsw[b,1,r]*dmul(Gamma{10},si[b,r]) + Fcsw[b,2,r]*dmul(Gamma{11},si[b,r]) + Fcsw[b,3,r]*dmul(Gamma{12},si[b,r]) + +Fcsw[b,4,r]*dmul(Gamma{15},si[b,r]) + Fcsw[b,5,r]*dmul(Gamma{14},si[b,r]) + Fcsw[b,6,r]*dmul(Gamma{13},si[b,r])) so[b,r] -= 0.5*(th[1]*gpmul(Pgamma{1,-1},U[b,1,r],si[bu1,ru1]) +conj(th[1])*gdagpmul(Pgamma{1,+1},U[bd1,1,rd1],si[bd1,rd1]) + th[2]*gpmul(Pgamma{2,-1},U[b,2,r],si[bu2,ru2]) +conj(th[2])*gdagpmul(Pgamma{2,+1},U[bd2,2,rd2],si[bd2,rd2]) + @@ -193,7 +225,7 @@ function krnl_Dwimpr!(so, U, si, Fcsw, m0, th, csw, lp::SpaceParm{4,6,B,D}) wher return nothing end -function krnl_Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D} +function krnl_Dw!(so, U, si, m0, tm, th, lp::SpaceParm{4,6,B,D}) where {B,D} b = Int64(CUDA.threadIdx().x); r = Int64(CUDA.blockIdx().x) @@ -208,7 +240,7 @@ function krnl_Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D} @inbounds begin - so[b,r] = (4+m0)*si[b,r] + so[b,r] = (4+m0)*si[b,r] + im*tm*dmul(Gamma{5},si[b,r]) so[b,r] -= 0.5*(th[1]*gpmul(Pgamma{1,-1},U[b,1,r],si[bu1,ru1]) +conj(th[1])*gdagpmul(Pgamma{1,+1},U[bd1,1,rd1],si[bd1,rd1]) + th[2]*gpmul(Pgamma{2,-1},U[b,2,r],si[bu2,ru2]) +conj(th[2])*gdagpmul(Pgamma{2,+1},U[bd2,2,rd2],si[bd2,rd2]) + @@ -225,13 +257,13 @@ function Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::Union{SpacePa if abs(dpar.csw) > 1.0E-10 @timeit "Dw" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dwimpr!(so, U, si, dws.csw, dpar.m0, dpar.th, dpar.csw, dpar.ct, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dwimpr!(so, U, si, dws.csw, dpar.m0, dpar.tm, dpar.th, dpar.csw, dpar.ct, lp) end end else @timeit "Dw" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw!(so, U, si, dpar.m0, dpar.th, dpar.ct, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_Dw!(so, U, si, dpar.m0, dpar.tm, dpar.th, dpar.ct, lp) end end end @@ -239,7 +271,7 @@ function Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::Union{SpacePa return nothing end -function krnl_Dwimpr!(so, U, si, Fcsw, m0, th, csw, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D} +function krnl_Dwimpr!(so, U, si, Fcsw, m0, tm, th, csw, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D} # The field si is assumed to be zero at t = 0 @@ -258,8 +290,8 @@ function krnl_Dwimpr!(so, U, si, Fcsw, m0, th, csw, ct, lp::Union{SpaceParm{4,6, @inbounds begin - so[b,r] = (4+m0)*si[b,r] + 0.5*csw*im*( Fcsw[b,1,r]*dmul(Gamma{10},si[b,r]) + Fcsw[b,2,r]*dmul(Gamma{11},si[b,r]) + Fcsw[b,3,r]*dmul(Gamma{12},si[b,r]) - +Fcsw[b,4,r]*dmul(Gamma{15},si[b,r]) + Fcsw[b,5,r]*dmul(Gamma{14},si[b,r]) + Fcsw[b,6,r]*dmul(Gamma{13},si[b,r])) + so[b,r] = (4+m0)*si[b,r] + im*tm*dmul(Gamma{5},si[b,r]) + 0.5*csw*im*( Fcsw[b,1,r]*dmul(Gamma{10},si[b,r]) + Fcsw[b,2,r]*dmul(Gamma{11},si[b,r]) + Fcsw[b,3,r]*dmul(Gamma{12},si[b,r]) + +Fcsw[b,4,r]*dmul(Gamma{15},si[b,r]) + Fcsw[b,5,r]*dmul(Gamma{14},si[b,r]) + Fcsw[b,6,r]*dmul(Gamma{13},si[b,r])) so[b,r] -= 0.5*(th[1]*gpmul(Pgamma{1,-1},U[b,1,r],si[bu1,ru1]) +conj(th[1])*gdagpmul(Pgamma{1,+1},U[bd1,1,rd1],si[bd1,rd1]) + @@ -276,7 +308,7 @@ function krnl_Dwimpr!(so, U, si, Fcsw, m0, th, csw, ct, lp::Union{SpaceParm{4,6, return nothing end -function krnl_Dw!(so, U, si, m0, th, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D} +function krnl_Dw!(so, U, si, m0, tm, th, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D} # The field si is assumed to be zero at t = 0 @@ -295,7 +327,7 @@ function krnl_Dw!(so, U, si, m0, th, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},S @inbounds begin - so[b,r] = (4+m0)*si[b,r] + so[b,r] = (4+m0)*si[b,r] + im*tm*dmul(Gamma{5},si[b,r]) so[b,r] -= 0.5*(th[1]*gpmul(Pgamma{1,-1},U[b,1,r],si[bu1,ru1]) +conj(th[1])*gdagpmul(Pgamma{1,+1},U[bd1,1,rd1],si[bd1,rd1]) + th[2]*gpmul(Pgamma{2,-1},U[b,2,r],si[bu2,ru2]) +conj(th[2])*gdagpmul(Pgamma{2,+1},U[bd2,2,rd2],si[bd2,rd2]) + th[3]*gpmul(Pgamma{3,-1},U[b,3,r],si[bu3,ru3]) +conj(th[3])*gdagpmul(Pgamma{3,+1},U[bd3,3,rd3],si[bd3,rd3]) + @@ -310,18 +342,24 @@ function krnl_Dw!(so, U, si, m0, th, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},S return nothing end +""" + function g5Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4,6,B,D}) + +Computes \`\` \\gamma_5 \`\` times the Dirac operator (with the Wilson term) with gauge field U and parameters `dpar` of the field `si` and stores it in `so`. +If `dpar.csw` is different from zero, the clover term should be stored in `dws.csw` via the Csw! function and is automatically included in the operator. +""" function g5Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4,6,B,D}) where {B,D} if abs(dpar.csw) > 1.0E-10 @timeit "g5Dw" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(so, U, si, dws.csw, dpar.m0, dpar.th, dpar.csw, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(so, U, si, dws.csw, dpar.m0, dpar.tm, dpar.th, dpar.csw, lp) end end else @timeit "g5Dw" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, si, dpar.m0, dpar.th, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, si, dpar.m0, dpar.tm, dpar.th, lp) end end end @@ -329,7 +367,7 @@ function g5Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4 return nothing end -function krnl_g5Dwimpr!(so, U, si, Fcsw, m0, th, csw, lp::SpaceParm{4,6,B,D}) where {B,D} +function krnl_g5Dwimpr!(so, U, si, Fcsw, m0, tm, th, csw, lp::SpaceParm{4,6,B,D}) where {B,D} b = Int64(CUDA.threadIdx().x); r = Int64(CUDA.blockIdx().x) @@ -352,13 +390,13 @@ function krnl_g5Dwimpr!(so, U, si, Fcsw, m0, th, csw, lp::SpaceParm{4,6,B,D}) wh th[3]*gpmul(Pgamma{3,-1},U[b,3,r],si[bu3,ru3]) +conj(th[3])*gdagpmul(Pgamma{3,+1},U[bd3,3,rd3],si[bd3,rd3]) + th[4]*gpmul(Pgamma{4,-1},U[b,4,r],si[bu4,ru4]) +conj(th[4])*gdagpmul(Pgamma{4,+1},U[bd4,4,rd4],si[bd4,rd4]) ) - so[b,r] = dmul(Gamma{5}, so[b,r]) + so[b,r] = dmul(Gamma{5}, so[b,r])+ im*tm*si[b,r] end return nothing end -function krnl_g5Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D} +function krnl_g5Dw!(so, U, si, m0, tm, th, lp::SpaceParm{4,6,B,D}) where {B,D} b = Int64(CUDA.threadIdx().x); r = Int64(CUDA.blockIdx().x) @@ -380,7 +418,7 @@ function krnl_g5Dw!(so, U, si, m0, th, lp::SpaceParm{4,6,B,D}) where {B,D} th[3]*gpmul(Pgamma{3,-1},U[b,3,r],si[bu3,ru3]) +conj(th[3])*gdagpmul(Pgamma{3,+1},U[bd3,3,rd3],si[bd3,rd3]) + th[4]*gpmul(Pgamma{4,-1},U[b,4,r],si[bu4,ru4]) +conj(th[4])*gdagpmul(Pgamma{4,+1},U[bd4,4,rd4],si[bd4,rd4]) ) - so[b,r] = dmul(Gamma{5}, so[b,r]) + so[b,r] = dmul(Gamma{5}, so[b,r]) + im*tm*si[b,r] end return nothing @@ -391,13 +429,13 @@ function g5Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::Union{Space if abs(dpar.csw) > 1.0E-10 @timeit "g5Dw" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(so, U, si, dws.csw, dpar.m0, dpar.th, dpar.csw, dpar.ct, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(so, U, si, dws.csw, dpar.m0, dpar.tm, dpar.th, dpar.csw, dpar.ct, lp) end end else @timeit "g5Dw" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, si, dpar.m0, dpar.th, dpar.ct, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, si, dpar.m0, dpar.tm, dpar.th, dpar.ct, lp) end end end @@ -405,7 +443,7 @@ function g5Dw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::Union{Space return nothing end -function krnl_g5Dwimpr!(so, U, si, Fcsw, m0, th, csw, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D} +function krnl_g5Dwimpr!(so, U, si, Fcsw, m0, tm, th, csw, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D} # The field si is assumed to be zero at t = 0 @@ -439,12 +477,12 @@ function krnl_g5Dwimpr!(so, U, si, Fcsw, m0, th, csw, ct, lp::Union{SpaceParm{4, end end - so[b,r] = dmul(Gamma{5}, so[b,r]) + so[b,r] = dmul(Gamma{5}, so[b,r])+ im*tm*si[b,r] return nothing end -function krnl_g5Dw!(so, U, si, m0, th, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D} +function krnl_g5Dw!(so, U, si, m0, tm, th, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D} # The field si is assumed to be zero at t = 0 @@ -475,11 +513,17 @@ function krnl_g5Dw!(so, U, si, m0, th, ct, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D} end end - so[b,r] = dmul(Gamma{5}, so[b,r]) + so[b,r] = dmul(Gamma{5}, so[b,r]) + im*tm*si[b,r] return nothing end +""" + function DwdagDw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpaceParm{4,6,B,D}) + +Applies the operator \`\` \\gamma_5 D_w \`\` twice to `si` and stores the result in `so`. This is equivalent to appling the operator \`\` D_w^\\dagger D_w \`\` +The Dirac operator is the same as in the functions `Dw!` and `g5Dw!` +""" function DwdagDw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D} if abs(dpar.csw) > 1.0E-10 @@ -487,30 +531,32 @@ function DwdagDw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::Union{Sp @timeit "g5Dw" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(dws.st, U, si, dws.csw, dpar.m0, dpar.th, dpar.csw, dpar.ct, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(dws.st, U, si, dws.csw, dpar.m0, dpar.tm, dpar.th, dpar.csw, dpar.ct, lp) end end - + SF_bndfix!(dws.st,lp) @timeit "g5Dw" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(so, U, dws.st, dws.csw, dpar.m0, dpar.th, dpar.csw, dpar.ct, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(so, U, dws.st, dws.csw, dpar.m0, -dpar.tm, dpar.th, dpar.csw, dpar.ct, lp) end end + SF_bndfix!(so,lp) end else @timeit "DwdagDw" begin @timeit "g5Dw" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(dws.st, U, si, dpar.m0, dpar.th, dpar.ct, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(dws.st, U, si, dpar.m0, dpar.tm, dpar.th, dpar.ct, lp) end end - + SF_bndfix!(dws.st,lp) @timeit "g5Dw" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, dws.st, dpar.m0, dpar.th, dpar.ct, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, dws.st, dpar.m0, -dpar.tm, dpar.th, dpar.ct, lp) end end + SF_bndfix!(so,lp) end end @@ -524,13 +570,13 @@ function DwdagDw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpacePar @timeit "g5Dw" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(dws.st, U, si, dws.csw, dpar.m0, dpar.th, dpar.csw, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(dws.st, U, si, dws.csw, dpar.m0, dpar.tm, dpar.th, dpar.csw, lp) end end @timeit "g5Dw" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(so, U, dws.st, dws.csw, dpar.m0, dpar.th, dpar.csw, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dwimpr!(so, U, dws.st, dws.csw, dpar.m0, -dpar.tm, dpar.th, dpar.csw, lp) end end end @@ -539,13 +585,13 @@ function DwdagDw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpacePar @timeit "g5Dw" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(dws.st, U, si, dpar.m0, dpar.th, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(dws.st, U, si, dpar.m0, dpar.tm, dpar.th, lp) end end @timeit "g5Dw" begin CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, dws.st, dpar.m0, dpar.th, lp) + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_g5Dw!(so, U, dws.st, dpar.m0, -dpar.tm, dpar.th, lp) end end end @@ -554,12 +600,27 @@ function DwdagDw!(so, U, si, dpar::DiracParam, dws::DiracWorkspace, lp::SpacePar return nothing end +""" + function mtwmdpar(dpar::DiracParam) +Returns `dpar` with oposite value of the twisted mass. +""" +function mtwmdpar(dpar::DiracParam{P,R}) where {P,R} + return DiracParam{P}(R,dpar.m0,dpar.csw,dpar.th,-dpar.tm,dpar.ct) +end + + +""" + SF_bndfix!(sp, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) + +Sets all the values of `sp` in the first time slice to zero. +""" function SF_bndfix!(sp, lp::Union{SpaceParm{4,6,BC_SF_ORBI,D},SpaceParm{4,6,BC_SF_AFWB,D}}) where {D} - CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_sfbndfix!(sp, lp) + @timeit "SF boundary fix" begin + CUDA.@sync begin + CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_sfbndfix!(sp, lp) + end end - return nothing end @@ -643,7 +704,7 @@ function krnl_assign_pf_su2!(f::AbstractArray, p , lp::SpaceParm, t::Int64) return nothing end -export Dw!, g5Dw!, DwdagDw!, SF_bndfix!, Csw!, pfrandomize! +export Dw!, g5Dw!, DwdagDw!, SF_bndfix!, Csw!, pfrandomize!, mtwmdpar include("DiracIO.jl") diff --git a/src/Groups/AlgebraU2.jl b/src/Groups/AlgebraU2.jl index 4e88638..48a6ec9 100644 --- a/src/Groups/AlgebraU2.jl +++ b/src/Groups/AlgebraU2.jl @@ -1,12 +1,23 @@ +""" + struct U2alg{T} <: Algebra + +Elements of the `U(2)` Algebra. The type `T <: AbstractFloat` can be used to define single or double precision elements. +""" struct U2alg{T} <: Algebra u11::T u22::T u12::Complex{T} end + +""" + antsym(a::SU2{T}) where T <: AbstractFloat + +Returns the antisymmetrization of the SU2 element `a`, that is `\`\` `a - a^{\\dagger}` `\`. This method returns al element of `U2alg{T}`. +""" function antsym(a::SU2{T}) where T <: AbstractFloat return U2alg{T}(2.0*imag(a.t1),-2.0*imag(a.t1),2.0*a.t2) end diff --git a/src/Groups/AlgebraU3.jl b/src/Groups/AlgebraU3.jl index fa03579..ba0bbd4 100644 --- a/src/Groups/AlgebraU3.jl +++ b/src/Groups/AlgebraU3.jl @@ -1,6 +1,10 @@ +""" + struct U3alg{T} <: Algebra +Elements of the `U(3)` Algebra. The type `T <: AbstractFloat` can be used to define single or double precision elements. +""" struct U3alg{T} <: Algebra u11::T u22::T @@ -10,6 +14,11 @@ struct U3alg{T} <: Algebra u23::Complex{T} end +""" + antsym(a::SU3{T}) where T <: AbstractFloat + +Returns the antisymmetrization of the SU3 element `a`, that is `\`\` `a - a^{\\dagger}` `\`. This method returns al element of `U3alg{T}`. +""" function antsym(a::SU3{T}) where T <: AbstractFloat t1 = 2.0*imag(a.u11) t2 = 2.0*imag(a.u22) diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl index 198da9a..63c408d 100644 --- a/src/LatticeGPU.jl +++ b/src/LatticeGPU.jl @@ -46,19 +46,19 @@ export Eoft_clover, Eoft_plaq, Qtop export FlowIntr, wfl_euler, zfl_euler, wfl_rk2, zfl_rk2, wfl_rk3, zfl_rk3 export flw, flw_adapt export sfcoupling, bndfield, setbndfield -export import_lex64, import_cern64, import_bsfqcd, save_cnfg, read_cnfg +export import_lex64, import_cern64, import_bsfqcd, save_cnfg, read_cnfg, read_gp include("Spinors/Spinors.jl") using .Spinors -export Spinor, Pgamma +export Spinor, Pgamma, Gamma export imm, mimm export pmul, gpmul, gdagpmul, dmul include("Dirac/Dirac.jl") using .Dirac export DiracWorkspace, DiracParam -export Dw!, g5Dw!, DwdagDw!, SF_bndfix!, Csw!, pfrandomize! +export Dw!, g5Dw!, DwdagDw!, SF_bndfix!, Csw!, pfrandomize!, mtwmdpar export read_prop, save_prop, read_dpar include("Solvers/Solvers.jl") diff --git a/src/Solvers/CG.jl b/src/Solvers/CG.jl index 8e86c7b..ffec5d4 100644 --- a/src/Solvers/CG.jl +++ b/src/Solvers/CG.jl @@ -9,11 +9,6 @@ ### created: Tue Nov 30 11:10:57 2021 ### -""" - function CG! - -Solves the linear equation `Ax = si` -""" function krnl_dot!(sum,fone,ftwo) b=Int64(CUDA.threadIdx().x) r=Int64(CUDA.blockIdx().x) @@ -23,7 +18,7 @@ function krnl_dot!(sum,fone,ftwo) return nothing end -function field_dot(fone::AbstractArray,ftwo::AbstractArray,sumf,lp) where {T} +function field_dot(fone::AbstractArray,ftwo::AbstractArray,sumf,lp) CUDA.@sync begin CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_dot!(sumf,fone,ftwo) @@ -32,6 +27,12 @@ function field_dot(fone::AbstractArray,ftwo::AbstractArray,sumf,lp) where {T} return sum(sumf) end + +""" + function CG!(si, U, A, dpar::DiracParam, lp::SpaceParm, dws::DiracWorkspace{T}, maxiter::Int64 = 10, tol=1.0) + +Solves the linear equation `Ax = si` +""" function CG!(si, U, A, dpar::DiracParam, lp::SpaceParm, dws::DiracWorkspace{T}, maxiter::Int64 = 10, tol=1.0) where {T} dws.sr .= si @@ -74,4 +75,4 @@ function CG!(si, U, A, dpar::DiracParam, lp::SpaceParm, dws::DiracWorkspace{T}, end return niter -end \ No newline at end of file +end diff --git a/src/Solvers/Propagators.jl b/src/Solvers/Propagators.jl index dc42f5b..58008b2 100644 --- a/src/Solvers/Propagators.jl +++ b/src/Solvers/Propagators.jl @@ -25,7 +25,7 @@ function propagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::Space CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(dws.sp) end - g5Dw!(pro,U,dws.sp,dpar,dws,lp) + g5Dw!(pro,U,dws.sp,mtwmdpar(dpar),dws,lp) CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol) return nothing @@ -46,11 +46,7 @@ function propagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::Space CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(dws.sp) end - g5Dw!(pro,U,dws.sp,dpar,dws,lp) - - CUDA.@sync begin - CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(dws.sp) - end + g5Dw!(pro,U,dws.sp,mtwmdpar(dpar),dws,lp) CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol) return nothing @@ -60,7 +56,7 @@ end function bndpropagator!(pro,U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) -Saves the propagator in from the t=0 boundary to the bulk for the SF boundary conditions for a source with color 'c' and spin 's'. The factor c_t is included while the factor 1/sqrt(V) is not. +Saves the propagator from the t=0 boundary to the bulk for the SF boundary conditions for a source with color 'c' and spin 's' in 'pro'. The factor c_t is included while the factor 1/sqrt(V) is not. For the propagator from T to the bulk, use the function Tbndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) """ @@ -85,6 +81,7 @@ function bndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::Sp return nothing end + SF_bndfix!(pro,lp) fill!(dws.sp,zero(eltype(scalar_field(Spinor{4,SU3fund{Float64}},lp)))) CUDA.@sync begin @@ -95,15 +92,15 @@ function bndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::Sp CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(dws.sp) end - g5Dw!(pro,U,dpar.ct*dws.sp,dpar,dws,lp) + g5Dw!(pro,U,dpar.ct*dws.sp,mtwmdpar(dpar),dws,lp) CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol) - return pro + return nothing end """ - function Tbndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) + function Tbndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) Returns the propagator from the t=T boundary to the bulk for the SF boundary conditions for a source with color 'c' and spin 's'. The factor c_t is included while the factor 1/sqrt(V) is not. For the propagator from t=0 to the bulk, use the function bndpropagator(U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::SpaceParm{4,6,1,D}, maxiter::Int64, tol::Float64, c::Int64, s::Int64) @@ -128,7 +125,8 @@ function Tbndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::S return nothing end - + + SF_bndfix!(pro,lp) fill!(dws.sp,zero(eltype(scalar_field(Spinor{4,SU3fund{Float64}},lp)))) CUDA.@sync begin @@ -139,10 +137,10 @@ function Tbndpropagator!(pro, U, dpar::DiracParam{T}, dws::DiracWorkspace, lp::S CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(dws.sp) end - g5Dw!(pro,U,dpar.ct*dws.sp,dpar,dws,lp) + g5Dw!(pro,U,dpar.ct*dws.sp,mtwmdpar(dpar),dws,lp) CG!(pro,U,DwdagDw!,dpar,lp,dws,maxiter,tol) - return pro + return nothing end diff --git a/src/Spinors/Spinors.jl b/src/Spinors/Spinors.jl index 6b8a46b..1cd29ea 100644 --- a/src/Spinors/Spinors.jl +++ b/src/Spinors/Spinors.jl @@ -169,7 +169,7 @@ end """ - gpmul(pgamma{N,S}, g::G, a::Spinor) G <: Group + gpmul(Pgamma{N,S}, g::G, a::Spinor) G <: Group Returns ``g(1+s\\gamma_N)a`` """ @@ -226,7 +226,7 @@ end end """ - gdagpmul(pgamma{N,S}, g::G, a::Spinor) G <: Group + gdagpmul(Pgamma{N,S}, g::G, a::Spinor) G <: Group Returns ``g^+ (1+s\\gamma_N)a`` """ @@ -284,33 +284,33 @@ end # dummy structs for dispatch: -# Basis of \\Gamma_n +# Basis of \\gamma_n struct Gamma{N} end """ - dmul(n::Int64, a::Spinor) + dmul(Gamma{n}, a::Spinor) -Returns ``\\Gamma_n a`` +Returns ``\\gamma_n a`` -indexing for Dirac basis ``\\Gamma_n``: +indexing for Dirac basis ``\\gamma_n``: - 1 gamma1 - 2 gamma2 - 3 gamma3 - 4 gamma0 - 5 gamma5 - 6 gamma1 gamma5 - 7 gamma2 gamma5 - 8 gamma3 gamma5 - 9 gamma0 gamma5 -10 sigma01 -11 sigma02 -12 sigma03 -13 sigma21 -14 sigma32 -15 sigma31 -16 identity + 1 gamma1; + 2 gamma2; + 3 gamma3; + 4 gamma0; + 5 gamma5; + 6 gamma1 gamma5; + 7 gamma2 gamma5; + 8 gamma3 gamma5; + 9 gamma0 gamma5; +10 sigma01; +11 sigma02; +12 sigma03; +13 sigma21; +14 sigma32; +15 sigma31; +16 identity; """ @inline dmul(::Type{Gamma{1}}, a::Spinor{NS,G}) where {NS,G} = Spinor{NS,G}((mimm(a.s[4]), mimm(a.s[3]), imm(a.s[2]), imm(a.s[1]))) diff --git a/src/YM/YM.jl b/src/YM/YM.jl index bdfc098..16165e4 100644 --- a/src/YM/YM.jl +++ b/src/YM/YM.jl @@ -179,6 +179,6 @@ include("YMsf.jl") export sfcoupling, bndfield, setbndfield include("YMio.jl") -export import_lex64, import_cern64, import_bsfqcd, save_cnfg, read_cnfg +export import_lex64, import_cern64, import_bsfqcd, save_cnfg, read_cnfg, read_gp end diff --git a/src/YM/YMio.jl b/src/YM/YMio.jl index 41e6132..f492497 100644 --- a/src/YM/YMio.jl +++ b/src/YM/YMio.jl @@ -297,3 +297,50 @@ function import_cern64(fname, ibc, lp::SpaceParm; log=true) return CuArray(Ucpu) end + + + +""" + read_gp(fname::String) + +Reads Gauge parameters from file `fname` using the native (BDIO) format. Returns GaugeParm and SpaceParm. +""" +function read_gp(fname::String) + + UID_HDR = 14 + fb = BDIO_open(fname, "r") + while BDIO_get_uinfo(fb) != UID_HDR + BDIO_seek!(fb) + end + ihdr = Vector{Int32}(undef, 2) + BDIO_read(fb, ihdr) + if (ihdr[1] != convert(Int32, 1653996111)) && (ihdr[2] != convert(Int32, 2)) + error("Wrong file format [header]") + end + + run = BDIO.BDIO_read_str(fb) + + while BDIO_get_uinfo(fb) != 1 + BDIO_seek!(fb) + end + + ifoo = Vector{Int32}(undef, 4) + BDIO_read(fb, ifoo) + ndim = convert(Int64, ifoo[1]) + npls = convert(Int64, round(ndim*(ndim-1)/2)) + ibc = convert(Int64, ifoo[2]) + nf = ifoo[4] + + ifoo = Vector{Int32}(undef, ndim+convert(Int32, npls)) + BDIO_read(fb, ifoo) + iL = ntuple(i -> convert(Int64, ifoo[i]),ndim) + ntw = ntuple(i -> convert(Int64, ifoo[i+ndim]), npls) + + dfoo = Vector{Float64}(undef, 4) + BDIO_read(fb, dfoo) + + lp = SpaceParm{ndim}(iL, (4,4,4,4), ibc, ntw) + gp = GaugeParm{Float64}(SU3{Float64}, dfoo[1], dfoo[2]) + + return gp, lp +end diff --git a/test/dirac/test_fp_fa.jl b/test/dirac/test_fp_fa.jl index 6610bb4..ad93afb 100644 --- a/test/dirac/test_fp_fa.jl +++ b/test/dirac/test_fp_fa.jl @@ -8,59 +8,84 @@ using TimerOutputs function fP_test(;theta = (0.5,0.7,1.0,0.0), m = 1.3, size = (8,8,8,16),prec = 1.0e-16) -@timeit "fP inversion (x12)" begin + @timeit "fP inversion (x12)" begin -lp = SpaceParm{4}(size,(4,4,4,4),1,(0,0,0,0,0,0)); -exptheta = exp.(im.*theta./lp.iL); + lp = SpaceParm{4}(size,(4,4,4,4),1,(0,0,0,0,0,0)); + exptheta = exp.(im.*theta./lp.iL); -dpar = DiracParam{Float64}(SU3fund,m,0.0,exptheta,1.0); -dws = DiracWorkspace(SU3fund,Float64,lp); + dpar = DiracParam{Float64}(SU3fund,m,0.0,exptheta,0.0,1.0); + dws = DiracWorkspace(SU3fund,Float64,lp); -U = fill!(vector_field(SU3{Float64},lp),one(SU3{Float64})); -psi = scalar_field(Spinor{4,SU3fund{Float64}},lp); + U = fill!(vector_field(SU3{Float64},lp),one(SU3{Float64})); + psi = scalar_field(Spinor{4,SU3fund{Float64}},lp); -res = zeros(lp.iL[4]) + res = zeros(lp.iL[4]) -for s in 1:4 for c in 1:3 - bndpropagator!(psi,U,dpar,dws,lp,1000,prec,c,s); + for s in 1:4 for c in 1:3 + bndpropagator!(psi,U,dpar,dws,lp,1000,prec,c,s); - for t in 1:lp.iL[4] - #for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3] - i=abs(rand(Int))%lp.iL[1] +1;j=abs(rand(Int))%lp.iL[2] +1;k=abs(rand(Int))%lp.iL[3] +1; - CUDA.@allowscalar (res[t] += norm2(psi[point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)...])/2) - #end end end - #res[t] = res[t]/(lp.iL[1]*lp.iL[2]*lp.iL[3]) + for t in 1:lp.iL[4] + #for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3] + i=abs(rand(Int))%lp.iL[1] +1;j=abs(rand(Int))%lp.iL[2] +1;k=abs(rand(Int))%lp.iL[3] +1; + CUDA.@allowscalar (res[t] += norm2(psi[point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)...])/2) + #end end end + #res[t] = res[t]/(lp.iL[1]*lp.iL[2]*lp.iL[3]) + + end + + end end end -end end + @timeit "fP analitical solution" begin -end + #THEORETICAL SOLUTION: hep-lat/9606016 eq (2.33) -@timeit "fP analitical solution" begin + res_th = zeros(lp.iL[4]) - #THEORETICAL SOLUTION: hep-lat/9606016 eq (2.33) + pp3 = ntuple(i -> theta[i]/lp.iL[i],3) + omega = 2 * asinh(0.5* sqrt(( sum((sin.(pp3)).^2) + (m + 2*(sum((sin.(pp3./2)).^2) ))^2) / (1+m+2*(sum((sin.(pp3./2)).^2) )) ) ) + pp = (-im*omega,pp3...) + Mpp = m + 2* sum((sin.(pp./2)).^2) + Rpp = Mpp*(1-exp(-2*omega*lp.iL[4])) + sinh(omega) * (1+exp(-2*omega*lp.iL[4])) - res_th = zeros(lp.iL[4]) + for i in 2:lp.iL[4] + res_th[i] = (2*3*sinh(omega)/(Rpp^2)) * ( (Mpp + sinh(omega))*exp(-2*omega*(i-1)) - (Mpp - sinh(omega))*exp(-2*omega*(2*lp.iL[4]- (i - 1))) ) + end - pp3 = ntuple(i -> theta[i]/lp.iL[i],3) - omega = 2 * asinh(0.5* sqrt(( sum((sin.(pp3)).^2) + (m + 2*(sum((sin.(pp3./2)).^2) ))^2) / (1+m+2*(sum((sin.(pp3./2)).^2) )) ) ) - pp = (-im*omega,pp3...) - Mpp = m + 2* sum((sin.(pp./2)).^2) - Rpp = Mpp*(1-exp(-2*omega*lp.iL[4])) + sinh(omega) * (1+exp(-2*omega*lp.iL[4])) - - for i in 2:lp.iL[4] - res_th[i] = (2*3*sinh(omega)/(Rpp^2)) * ( (Mpp + sinh(omega))*exp(-2*omega*(i-1)) - (Mpp - sinh(omega))*exp(-2*omega*(2*lp.iL[4]- (i - 1))) ) end - -end return sum(abs.(res-res_th)) end function fA_test(;theta = (0.5,0.7,1.0,0.0), m = 1.3, size = (8,8,8,16),prec = 1.0e-16) -@timeit "fA inversion (x12)" begin + @timeit "fA inversion (x12)" begin + + lp = SpaceParm{4}(size,(4,4,4,4),1,(0,0,0,0,0,0)); + exptheta = exp.(im.*theta./lp.iL); + + dpar = DiracParam{Float64}(SU3fund,m,0.0,exptheta,0.0,1.0); + dws = DiracWorkspace(SU3fund,Float64,lp); + + U = fill!(vector_field(SU3{Float64},lp),one(SU3{Float64})); + psi = scalar_field(Spinor{4,SU3fund{Float64}},lp); + + res = im*zeros(lp.iL[4]) + + for s in 1:4 for c in 1:3 + bndpropagator!(psi,U,dpar,dws,lp,1000,prec,c,s); + + for t in 1:lp.iL[4] + #for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3] + i=abs(rand(Int))%lp.iL[1] +1;j=abs(rand(Int))%lp.iL[2] +1;k=abs(rand(Int))%lp.iL[3] +1; + CUDA.@allowscalar (res[t] += -dot(psi[point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)...],dmul(Gamma{4},psi[point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)...]))/2) + #end end end + #res[t] = res[t]/(lp.iL[1]*lp.iL[2]*lp.iL[3]) + + end + + end end lp = SpaceParm{4}(size,(4,4,4,4),1,(0,0,0,0,0,0)); exptheta = exp.(im.*theta./lp.iL); @@ -103,11 +128,26 @@ end res_th[i] = (6/(Rpp^2)) * ( 2*(Mpp - sinh(omega))*(Mpp + sinh(omega))*exp(-2*omega*lp.iL[4]) - Mpp*((Mpp + sinh(omega))*exp(-2*omega*(i-1)) + (Mpp - sinh(omega))*exp(-2*omega*(2*lp.iL[4]- (i - 1))))) end - -end - + #THEORETICAL SOLUTION: hep-lat/9606016 eq (2.32) + + @timeit "fA analitical solution" begin + res_th = zeros(lp.iL[4]) + + pp3 = ntuple(i -> theta[i]/lp.iL[i],3) + omega = 2 * asinh(0.5* sqrt(( sum((sin.(pp3)).^2) + (m + 2*(sum((sin.(pp3./2)).^2) ))^2) / (1+m+2*(sum((sin.(pp3./2)).^2) )) ) ) + pp = (-im*omega,pp3...) + Mpp = m + 2* sum((sin.(pp./2)).^2) + Rpp = Mpp*(1-exp(-2*omega*lp.iL[4])) + sinh(omega) * (1+exp(-2*omega*lp.iL[4])) + + for i in 2:lp.iL[4] + res_th[i] = (6/(Rpp^2)) * ( 2*(Mpp - sinh(omega))*(Mpp + sinh(omega))*exp(-2*omega*lp.iL[4]) + - Mpp*((Mpp + sinh(omega))*exp(-2*omega*(i-1)) + (Mpp - sinh(omega))*exp(-2*omega*(2*lp.iL[4]- (i - 1))))) + end + + end + return sum(abs.(res-res_th)) - + end diff --git a/test/dirac/test_solver_plw.jl b/test/dirac/test_solver_plw.jl index 195a1a3..a4de0c7 100644 --- a/test/dirac/test_solver_plw.jl +++ b/test/dirac/test_solver_plw.jl @@ -7,7 +7,7 @@ using LatticeGPU, CUDA, TimerOutputs function Dwpw_test(;p=0,s=1,c=1) lp = SpaceParm{4}((16,16,16,16), (4,4,4,4), 0, (0,0,0,0,0,0)) gp = GaugeParm{Float64}(SU3{Float64}, 6.0, 1.0) -dpar = DiracParam{Float64}(SU3fund,1.3,0.0,(1.0,1.0,1.0,1.0),0.0) +dpar = DiracParam{Float64}(SU3fund,1.3,0.0,(1.0,1.0,1.0,1.0),0.0,0.0) dws = DiracWorkspace(SU3fund,Float64,lp); p==0 ? p = Int.(round.(lp.iL.*rand(4),RoundUp)) : nothing @@ -96,15 +96,15 @@ end begin -dif = 0.0 +diff = 0.0 for i in 1:3 for j in 1:4 - dif += Dwpw_test(c=i,s=j) + global diff += Dwpw_test(c=i,s=j) end end -if dif < 1.0e-15 - print("Dwpl test passed with average error ", dif/12,"!\n") +if diff < 1.0e-15 + print("Dwpl test passed with average error ", diff/12,"!\n") else - error("Dwpl test failed with difference: ",dif,"\n") + error("Dwpl test failed with difference: ",diff,"\n") end diff --git a/test/dirac/test_solver_rand.jl b/test/dirac/test_solver_rand.jl index c06de0a..0714d8f 100644 --- a/test/dirac/test_solver_rand.jl +++ b/test/dirac/test_solver_rand.jl @@ -9,7 +9,7 @@ using CUDA, LatticeGPU, TimerOutputs lp = SpaceParm{4}((16,16,16,16), (4,4,4,4), 0, (0,0,0,0,0,0)) gp = GaugeParm{Float64}(SU3{Float64}, 6.0, 1.0) ymws = YMworkspace(SU3, Float64, lp) - dpar = DiracParam{Float64}(SU3fund,2.3,0.0,(1.0,1.0,1.0,1.0),0.0) + dpar = DiracParam{Float64}(SU3fund,2.3,0.0,(1.0,1.0,1.0,1.0),0.0,0.0) dws = DiracWorkspace(SU3fund,Float64,lp); randomize!(ymws.mom, lp, ymws) @@ -32,7 +32,8 @@ end CUDA.@sync begin CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnlg5!(rpsi) end -g5Dw!(prop,U,rpsi,dpar,dws,lp) + +g5Dw!(prop,U,rpsi,mtwmdpar(dpar),dws,lp) CG!(prop,U,DwdagDw!,dpar,lp,dws,10000,1.0e-14) Dw!(dws.sp,U,prop,dpar,dws,lp) diff --git a/test/runtests.jl b/test/runtests.jl index fe59102..2f68f70 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,3 +1,6 @@ -#include("SAD/test_sad.jl") +include("SAD/test_sad.jl") include("flow/test_adapt.jl") +include("dirac/test_fp_fa.jl") +include("dirac/test_solver_plw.jl") +include("dirac/test_solver_rand.jl")