From 0b97285ae6cfb4a2a0cc430313941129cfc8f4a5 Mon Sep 17 00:00:00 2001 From: GianlucaFuwa Date: Sun, 6 Oct 2024 00:22:55 +0200 Subject: [PATCH] (Almost) Done with distributed fields. EO-Precon'ed Wilson-Clover still elusive --- benchmark/bench_dirac.jl | 11 +- benchmark/bench_meas.jl | 12 +- metaqcd.jl | 83 ++++++ metaqcd_build.jl | 55 ++-- metaqcd_sim.jl | 50 ++-- src/diracoperators/DiracOperators.jl | 27 +- src/diracoperators/arnoldi.jl | 4 +- src/diracoperators/overlap.jl | 0 src/diracoperators/overlap_eo.jl | 0 src/diracoperators/staggered.jl | 3 +- src/diracoperators/staggered_eo.jl | 3 +- src/diracoperators/wilson.jl | 22 +- src/diracoperators/wilson_eo.jl | 131 ++++----- src/fields/Fields.jl | 11 +- src/fields/algebrafield.jl | 4 +- src/fields/distributed.jl | 8 +- src/fields/gaugefield.jl | 26 +- src/fields/paulifield.jl | 40 +++ src/fields/spinorfield.jl | 242 ---------------- src/fields/spinorfield_eo.jl | 264 ++++++++++++++++++ src/fields/tensorfield.jl | 78 +++--- src/forces/bias_force.jl | 2 +- src/forces/forces.jl | 2 +- src/forces/overlap_eo_force.jl | 0 src/forces/overlap_force.jl | 0 src/forces/wilson_eo_force.jl | 124 +++++--- src/io/bmw_format.jl | 8 +- src/io/ildg_format.jl | 8 +- src/main/Universe.jl | 10 +- src/main/runbuild.jl | 13 +- src/main/runsim.jl | 4 +- .../measure_topological_charge.jl | 46 +-- src/parameters/Parameters.jl | 24 +- src/solvers/cg.jl | 13 +- src/updates/heatbath.jl | 4 +- src/updates/hmc_integrators.jl | 1 + src/updates/tempering.jl | 12 +- src/utils/Utils.jl | 26 +- src/utils/clover_obj.jl | 19 -- src/utils/generators.jl | 82 +++--- src/utils/pauli.jl | 23 ++ src/utils/simd_vecops.jl | 76 ++--- test/runtests.jl | 10 +- test/test_derivative.jl | 2 +- test/test_fderivative.jl | 34 ++- test/test_mpi_derivative.jl | 127 --------- test/test_mpi_fderivative.jl | 134 --------- 47 files changed, 932 insertions(+), 946 deletions(-) create mode 100644 metaqcd.jl create mode 100644 src/diracoperators/overlap.jl create mode 100644 src/diracoperators/overlap_eo.jl create mode 100644 src/fields/paulifield.jl create mode 100644 src/fields/spinorfield_eo.jl create mode 100644 src/forces/overlap_eo_force.jl create mode 100644 src/forces/overlap_force.jl delete mode 100644 src/utils/clover_obj.jl create mode 100644 src/utils/pauli.jl delete mode 100644 test/test_mpi_derivative.jl delete mode 100644 test/test_mpi_fderivative.jl diff --git a/benchmark/bench_dirac.jl b/benchmark/bench_dirac.jl index 2f875d6d..92d71646 100644 --- a/benchmark/bench_dirac.jl +++ b/benchmark/bench_dirac.jl @@ -15,14 +15,21 @@ ops = ( StaggeredEOPreDiracOperator, ) +titles = ( + "Wilson", + # "Wilson (Even-Odd preconditioned)", + "Staggered", + "Staggered (Even-Odd preconditioned)", +) + Random.seed!(1234) N = 16 suite = BenchmarkGroup() -for dirac in ops - s = suite["$(dirac)"] = BenchmarkGroup() +for (i, dirac) in enumerate(ops) + s = suite["$(titles[i])"] = BenchmarkGroup() for T in (Float32, Float64) U = Gaugefield{CPU,T,WilsonGaugeAction}(N, N, N, N, 6.0) D = dirac(U, 0.01; csw=1.0) diff --git a/benchmark/bench_meas.jl b/benchmark/bench_meas.jl index f1c5c3fd..a5a84542 100644 --- a/benchmark/bench_meas.jl +++ b/benchmark/bench_meas.jl @@ -27,12 +27,12 @@ for T in (Float32, Float64) GA_methods = ["wilson", "symanzik_tree", "iwasaki", "dbw2"] m_gaction = GaugeActionMeasurement(U, GA_methods=GA_methods) - s["$(T), plaq"] = @benchmarkable(measure($m_plaq, $U, 1, 1)) - s["$(T), poly"] = @benchmarkable(measure($m_poly, $U, 1, 1)) - s["$(T), wilson"] = @benchmarkable(measure($m_wilson, $U, 1, 1)) - s["$(T), topo"] = @benchmarkable(measure($m_topo, $U, 1, 1)) - s["$(T), ed"] = @benchmarkable(measure($m_ed, $U, 1, 1)) - s["$(T), gaction"] = @benchmarkable(measure($m_gaction, $U, 1, 1)) + s["Avg Plaquette, $(T)"] = @benchmarkable(measure($m_plaq, $U, 1, 1)) + s["Polyakov Loop, $(T)"] = @benchmarkable(measure($m_poly, $U, 1, 1)) + s["Wilson Loops (2x2 + 4x4), $(T)"] = @benchmarkable(measure($m_wilson, $U, 1, 1)) + s["Top. Charge (Plaq + Clov + Imp), $(T)"] = @benchmarkable(measure($m_topo, $U, 1, 1)) + s["Energy Density, $(T)"] = @benchmarkable(measure($m_ed, $U, 1, 1)) + s["Gauge Action (W + LW + IW + DBW2), $(T)"] = @benchmarkable(measure($m_gaction, $U, 1, 1)) end end diff --git a/metaqcd.jl b/metaqcd.jl new file mode 100644 index 00000000..3741bfd4 --- /dev/null +++ b/metaqcd.jl @@ -0,0 +1,83 @@ +using Pkg +Pkg.activate(@__DIR__(); io=devnull) + +using MetaQCD.Utils +using MetaQCD: @level1, build_bias, run_sim + +function parse_args(args) + parameterfile = args[1] + @assert length(args) > 2 && isfile(parameterfile) """ + An existing parameter file has to be given as the first input, e.g.: + julia metaqcd.jl parameters.toml -mode=sim + + You either did not provide a file or the file you provided does not exist. + """ + + @assert count(x -> occursin("-mode", x), args) == 1 """ + The flag \"-mode\" has to be set after the parameter file. + Options are: + + \"-mode=sim\" to run a simulation with or without Metadynamics or + \"-mode=build\" for building a bias potential with possibly multiple walkers + """ + + mode = split(args[findfirst(x -> occursin("-mode", x), args)], "=")[2] + backend = try + split(args[findfirst(x -> occursin("-backend", x), args)], "=")[2] + catch _ + "cpu" + end + + return parameterfile, mode, backend +end + +parameterfile, mode, backend = parse_args(ARGS) + +@level1 "[ Mode: $(mode)\n" + +if backend != "cpu" + @level1 """ + If you get prompted to install a package here, make sure you do it in your \ + GLOBAL julia environment, i.e., not under a project environment + """ + + if backend == "cuda" + using CUDA + elseif backend ∈ ("rocm", "amd") + using AMDGPU + else + throw(ArgumentError( + """ + When a second input is given, it has to specify the backend to be used, \ + so the package can be loaded. + Note, that the backend also has to be set in the parameter file. + Supported backends are: + + - cpu + - cuda + - rocm + + Your input was \"$(backend)\" + """ + )) + end +end + +mpi_parallel() && @level1("[ $(mpi_size()) MPI processes are being used") + +if mode == "sim" + run_sim(parameterfile; backend=backend) +elseif mode == "build" + build_bias(parameterfile; backend=backend, mpi_multi_sim=with_mpi) +else + throw(ArgumentError( + """ + The supplied \"-mode\" is invalid. The two options are: + + \"-mode=sim\" to run a simulation with or without Metadynamics or + \"-mode=build\" for building a bias potential with possibly multiple walkers + + Your input was \"$(mode)\" + """ + )) +end diff --git a/metaqcd_build.jl b/metaqcd_build.jl index c86eb8fd..aa916b0d 100644 --- a/metaqcd_build.jl +++ b/metaqcd_build.jl @@ -2,50 +2,45 @@ using Pkg Pkg.activate(@__DIR__(); io=devnull) using MetaQCD.Utils -using MetaQCD: build_bias +using MetaQCD: @level1, build_bias using MetaQCD.Parameters: lower_case -COMM = mpi_init() -# @assert mpi_size() == 1 "metaqcd_sim can not be used with mpi, only metaqcd_build" - -if length(ARGS) == 0 && mpi_myrank() == 0 - error(""" - A parameter file has to be given as the first input: - julia metaqcd_build.jl parameters.toml - """) -end +@assert length(ARGS) != 0 """ +A parameter file has to be given as the first input: +julia metaqcd_sim.jl parameters.toml +""" backend = "cpu" if length(ARGS) == 2 - mpi_amroot() && @info( - "If you get prompted to install a package here, make sure you do it in your" * - "GLOBAL julia environment, i.e., not under a project environment" - ) + @level1 """ + If you get prompted to install a package here, make sure you do it in your \ + GLOBAL julia environment, i.e., not under a project environment + """ + if lower_case(ARGS[2]) == "cpu" elseif lower_case(ARGS[2]) == "cuda" using CUDA backend = "cuda" - elseif lower_case(ARGS[2]) == "rocm" + elseif lower_case(ARGS[2]) ∈ ("rocm", "amd") using AMDGPU backend = "rocm" else - mpi_myrank() == 0 && error(""" - When a second input is given, it has to specify the backend to be used, so the package can be loaded. - Note, that the backend also has to be set in the parameter file - Supported backends are: - - cpu - - cuda - - rocm - Your input was \"$(ARGS[2])\" - """) - end -end - -if mpi_parallel() && mpi_myrank() == 0 - println("$(mpi_size()) mpi processes will be used") + throw(ArgumentError( + """ + When a second input is given, it has to specify the backend to be used, \ + so the package can be loaded. + Note, that the backend also has to be set in the parameter file. + Supported backends are: + - cpu + - cuda + - rocm + Your input was \"$(ARGS[2])\" + """ + )) + end end -mpi_barrier() +mpi_parallel() && @level1("[ $(mpi_size()) MPI processes are being used") build_bias(ARGS[1]; backend=backend, mpi_multi_sim=with_mpi) diff --git a/metaqcd_sim.jl b/metaqcd_sim.jl index 85e5a093..f9a3d927 100644 --- a/metaqcd_sim.jl +++ b/metaqcd_sim.jl @@ -2,47 +2,45 @@ using Pkg Pkg.activate(@__DIR__(); io=devnull) using MetaQCD.Utils -using MetaQCD: run_sim +using MetaQCD: @level1, run_sim using MetaQCD.Parameters: lower_case -if length(ARGS) == 0 - error(""" - A parameter file has to be given as the first input: - julia metaqcd_sim.jl parameters.toml - """) -end +@assert length(ARGS) != 0 """ +A parameter file has to be given as the first input: +julia metaqcd_sim.jl parameters.toml +""" backend = "cpu" if length(ARGS) == 2 - @info( - "If you get prompted to install a package here, make sure you do it in your" * - "GLOBAL julia environment, i.e., not under a project environment" - ) + @level1 """ + If you get prompted to install a package here, make sure you do it in your \ + GLOBAL julia environment, i.e., not under a project environment + """ + if lower_case(ARGS[2]) == "cpu" elseif lower_case(ARGS[2]) == "cuda" using CUDA backend = "cuda" - elseif lower_case(ARGS[2]) == "rocm" + elseif lower_case(ARGS[2]) ∈ ("rocm", "amd") using AMDGPU backend = "rocm" else - error(""" - When a second input is given, it has to specify the backend to be used, so the package can be loaded. - Note, that the backend also has to be set in the parameter file. - Supported backends are: - - cpu - - cuda - - rocm - Your input was \"$(ARGS[2])\" - """) + throw(ArgumentError( + """ + When a second input is given, it has to specify the backend to be used, \ + so the package can be loaded. + Note, that the backend also has to be set in the parameter file. + Supported backends are: + - cpu + - cuda + - rocm + Your input was \"$(ARGS[2])\" + """ + )) end end -if mpi_parallel() && mpi_amroot() - println("[ $(mpi_size()) MPI processes are being used") -end - -mpi_barrier() +mpi_parallel() && @level1("[ $(mpi_size()) MPI processes are being used") run_sim(ARGS[1]; backend=backend) diff --git a/src/diracoperators/DiracOperators.jl b/src/diracoperators/DiracOperators.jl index 7df17c45..003e9ebb 100644 --- a/src/diracoperators/DiracOperators.jl +++ b/src/diracoperators/DiracOperators.jl @@ -22,10 +22,11 @@ using ..Solvers using ..Utils import KernelAbstractions as KA -import ..Fields: AbstractField, Gaugefield, Spinorfield, SpinorfieldEO, Tensorfield +import ..Fields: AbstractField, FieldTopology, Gaugefield, Paulifield, Spinorfield +import ..Fields: SpinorfieldEO, Tensorfield import ..Fields: check_dims, clear!, clover_square, dims, even_odd, gaussian_pseudofermions! import ..Fields: Clover, Checkerboard2, Sequential, @latmap, @latsum, set_source!, volume -import ..Fields: fieldstrength_eachsite!, fieldstrength_A_eachsite!, num_colors, num_dirac +import ..Fields: fieldstrength_eachsite!, num_colors, num_dirac import ..Fields: PeriodicBC, AntiPeriodicBC, apply_bc, create_bc, distributed_reduce abstract type AbstractDiracOperator end @@ -134,7 +135,7 @@ function Base.show(io::IO, ::MIME"text/plain", D::T) where {T<:AbstractDiracOper print(io, "$(typeof(D))", "(;") for fieldname in fieldnames(T) - if fieldname ∈ (:U, :temp) + if fieldname ∈ (:U, :temp, :D_diag, :D_oo_inv, :Fμν) continue else print(io, " ", fieldname, " = ", getfield(D, fieldname), ",") @@ -270,12 +271,12 @@ function Base.show(io::IO, ::MIME"text/plain", S::AbstractFermionAction{Nf}) whe """ | $(name)( - | Nf = $Nf - | MASS = $(S.D.mass) + | Nf: $Nf + | MASS: $(S.D.mass) """ ) - if S isa WilsonFermionAction #|| S isa WilsonEOPreFermionAction + if S isa WilsonFermionAction || S isa WilsonEOPreFermionAction print( io, """ @@ -289,10 +290,10 @@ function Base.show(io::IO, ::MIME"text/plain", S::AbstractFermionAction{Nf}) whe io, """ | BOUNDARY CONDITION (TIME): $(S.D.boundary_condition)) - | CG TOLERANCE (ACTION) = $(S.cg_tol_action) - | CG TOLERANCE (MD) = $(S.cg_tol_md) - | CG MAX ITERS (ACTION) = $(S.cg_maxiters_action) - | CG MAX ITERS (ACTION) = $(S.cg_maxiters_md) + | CG TOLERANCE (ACTION): $(S.cg_tol_action) + | CG TOLERANCE (MD): $(S.cg_tol_md) + | CG MAX ITERS (ACTION): $(S.cg_maxiters_action) + | CG MAX ITERS (ACTION): $(S.cg_maxiters_md) | RHMC INFO (Action): $(S.rhmc_info_action) | RHMC INFO (MD): $(S.rhmc_info_md)) """ @@ -307,12 +308,12 @@ function Base.show(io::IO, S::AbstractFermionAction{Nf}) where {Nf} """ | $(name)( - | Nf = $Nf - | MASS = $(S.D.mass) + | Nf: $Nf + | MASS: $(S.D.mass) """ ) - if S isa WilsonFermionAction #|| S isa WilsonEOPreFermionAction + if S isa WilsonFermionAction || S isa WilsonEOPreFermionAction print( io, """ diff --git a/src/diracoperators/arnoldi.jl b/src/diracoperators/arnoldi.jl index 8ad7b089..f11e6f2c 100644 --- a/src/diracoperators/arnoldi.jl +++ b/src/diracoperators/arnoldi.jl @@ -108,7 +108,9 @@ function partialschur!( n = LinearAlgebra.checksquare(D) maxdim = length(arnoldi.x) @assert nev ≥ 1 "nev cannot be less than 1, was $nev" - @assert nev ≤ mindim ≤ maxdim < n "nev ≤ mindim ≤ maxdim < n must hold, was $nev ≤ $mindim ≤ $maxdim < $n" + @assert nev ≤ mindim ≤ maxdim < n """ + nev ≤ mindim ≤ maxdim < n must hold, was $nev ≤ $mindim ≤ $maxdim < $n + """ _which = _symbol_to_target_meta(which) fill!(view(arnoldi.H, :, axes(arnoldi.H, 2)), zero(T)) reinitialize_meta!(arnoldi, 0) diff --git a/src/diracoperators/overlap.jl b/src/diracoperators/overlap.jl new file mode 100644 index 00000000..e69de29b diff --git a/src/diracoperators/overlap_eo.jl b/src/diracoperators/overlap_eo.jl new file mode 100644 index 00000000..e69de29b diff --git a/src/diracoperators/staggered.jl b/src/diracoperators/staggered.jl index f717df93..511fbf2d 100644 --- a/src/diracoperators/staggered.jl +++ b/src/diracoperators/staggered.jl @@ -9,7 +9,8 @@ condition in the time direction. If `csw ≠ 0`, a clover term is included. -This object cannot be directly applied to a fermion vector, since it lacks a gauge background. +This object cannot be directly applied to a fermion vector, since it lacks a gauge +background. A Wilson Dirac operator with gauge background is created by applying it to a `Gaugefield` `U` like `D_gauge = D(U)` diff --git a/src/diracoperators/staggered_eo.jl b/src/diracoperators/staggered_eo.jl index a897bd2d..b573960a 100644 --- a/src/diracoperators/staggered_eo.jl +++ b/src/diracoperators/staggered_eo.jl @@ -12,7 +12,8 @@ condition in the time direction. If `csw ≠ 0`, a clover term is included. -This object cannot be directly applied to a fermion vector, since it lacks a gauge background. +This object cannot be directly applied to a fermion vector, since it lacks a gauge +background. A Wilson Dirac operator with gauge background is created by applying it to a `Gaugefield` `U` like `D_gauge = D(U)` diff --git a/src/diracoperators/wilson.jl b/src/diracoperators/wilson.jl index fd30ec1b..a483a07a 100644 --- a/src/diracoperators/wilson.jl +++ b/src/diracoperators/wilson.jl @@ -103,7 +103,9 @@ struct WilsonFermionAction{Nf,TD,CT,RI1,RI2,RT,TX} <: AbstractFermionAction{Nf} rhmc_temps2 = nothing cg_temps = ntuple(_ -> Spinorfield(f), 4) else - @assert Nf == 1 "Nf should be 1 or 2 (was $Nf). If you want Nf > 2, use multiple actions" + @assert Nf == 1 """ + Nf should be 1 or 2 (was $Nf). If you want Nf > 2, use multiple actions + """ rhmc_lambda_low = rhmc_spectral_bound[1] rhmc_lambda_high = rhmc_spectral_bound[2] cg_temps = ntuple(_ -> Spinorfield(f), 2) @@ -265,8 +267,6 @@ function LinearAlgebra.mul!( ψ[site] = wilson_kernel(U, ϕ, site, mass_term, bc, T, Val(1)) end - update_halo!(ψ) - if has_clover_term(D) fac = T(-csw / 2) @@ -293,8 +293,6 @@ function LinearAlgebra.mul!( ψ[site] = wilson_kernel(U, ϕ, site, mass_term, bc, T, Val(-1)) end - update_halo!(ψ) - if has_clover_term(D) fac = T(-csw / 2) @@ -354,27 +352,27 @@ function clover_kernel(U, ϕ, site, fac, ::Type{T}) where {T} Cₙₘ = zero(ϕ[site]) C₁₂ = clover_square(U, 1, 2, site, 1) - F₁₂ = antihermitian(C₁₂) + F₁₂ = C₁₂ - C₁₂' Cₙₘ += cmvmul_color(F₁₂, σμν_spin_mul(ϕ[site], Val(1), Val(2))) C₁₃ = clover_square(U, 1, 3, site, 1) - F₁₃ = antihermitian(C₁₃) + F₁₃ = C₁₃ - C₁₃' Cₙₘ += cmvmul_color(F₁₃, σμν_spin_mul(ϕ[site], Val(1), Val(3))) C₁₄ = clover_square(U, 1, 4, site, 1) - F₁₄ = antihermitian(C₁₄) + F₁₄ = C₁₄ - C₁₄' Cₙₘ += cmvmul_color(F₁₄, σμν_spin_mul(ϕ[site], Val(1), Val(4))) C₂₃ = clover_square(U, 2, 3, site, 1) - F₂₃ = antihermitian(C₂₃) + F₂₃ = C₂₃ - C₂₃' Cₙₘ += cmvmul_color(F₂₃, σμν_spin_mul(ϕ[site], Val(2), Val(3))) C₂₄ = clover_square(U, 2, 4, site, 1) - F₂₄ = antihermitian(C₂₄) + F₂₄ = C₂₄ - C₂₄' Cₙₘ += cmvmul_color(F₂₄, σμν_spin_mul(ϕ[site], Val(2), Val(4))) C₃₄ = clover_square(U, 3, 4, site, 1) - F₃₄ = antihermitian(C₃₄) + F₃₄ = C₃₄ - C₃₄' Cₙₘ += cmvmul_color(F₃₄, σμν_spin_mul(ϕ[site], Val(3), Val(4))) - return Complex{T}(fac * im / 4) * Cₙₘ + return Complex{T}(fac * im / 8) * Cₙₘ end diff --git a/src/diracoperators/wilson_eo.jl b/src/diracoperators/wilson_eo.jl index 5d9065c4..bfd9cc70 100644 --- a/src/diracoperators/wilson_eo.jl +++ b/src/diracoperators/wilson_eo.jl @@ -39,14 +39,14 @@ struct WilsonEOPreDiracOperator{B,T,C,TF,TG,TX,TO,BC} <: AbstractDiracOperator function WilsonEOPreDiracOperator( f::AbstractField{B,T}, mass; bc_str="antiperiodic", r=1, csw=0, kwargs... ) where {B,T} - @assert r === 1 "Only r=1 in Wilson Dirac supported for now" + @assert r == 1 "Only r=1 in Wilson Dirac supported for now" κ = 1 / (2mass + 8) U = nothing C = csw == 0 ? false : true Fμν = C ? Tensorfield(f) : nothing temp = even_odd(Spinorfield(f)) - D_diag = WilsonEODiagonal(temp, mass, csw) - D_oo_inv = WilsonEODiagonal(temp, mass, csw; inverse=true) + D_diag = Paulifield(temp, csw) + D_oo_inv = Paulifield(temp, csw) boundary_condition = create_bc(bc_str, f.topology) TG = Nothing @@ -66,8 +66,8 @@ struct WilsonEOPreDiracOperator{B,T,C,TF,TG,TX,TO,BC} <: AbstractDiracOperator mass = D.mass csw = D.csw temp = even_odd(D.temp) - D_diag = WilsonEODiagonal(temp, mass, csw) - D_oo_inv = WilsonEODiagonal(temp, D.mass, csw; inverse=true) + D_diag = Paulifield(temp, csw) + D_oo_inv = Paulifield(temp, csw) Fμν = C ? Tensorfield(U) : nothing calc_diag!(D_diag, D_oo_inv, Fμν, U, mass) @@ -109,36 +109,6 @@ end @inline has_clover_term(::Daggered{W}) where {B,T,C,W<:WilsonEOPreDiracOperator{B,T,C}} = C @inline has_clover_term(::DdaggerD{W}) where {B,T,C,W<:WilsonEOPreDiracOperator{B,T,C}} = C -struct WilsonEODiagonal{B,T,M,C,A} <: AbstractField{B,T,M,A} # So we can overload LinearAlgebra.det on the even-odd diagonal - U::A - mass::Float64 - csw::Float64 - NX::Int64 # Number of lattice sites in the x-direction - NY::Int64 # Number of lattice sites in the y-direction - NZ::Int64 # Number of lattice sites in the z-direction - NT::Int64 # Number of lattice sites in the t-direction - NV::Int64 # Total number of lattice sites - NC::Int64 # Number of colors - ND::Int64 # Number of dirac indices - function WilsonEODiagonal( - f::TF, mass, csw; inverse=false - ) where {B,T,TF<:Union{Spinorfield{B,T},SpinorfieldEO{B,T}}} - NX, NY, NZ, NT = dims(f) - _NT = inverse ? NT ÷ 2 : NT - NC = num_colors(f) - ND = num_dirac(f) - NV = NX * NY * NZ * NT - C = csw == 0 ? false : true - U = KA.zeros(B(), SMatrix{6,6,Complex{T},36}, 2, NX, NY, NZ, _NT) - - A = typeof(U) - return new{B,T,false,C,A}(U, mass, csw, NX, NY, NZ, NT, NV, NC, ND) - end -end - -dims(D_diag::WilsonEODiagonal) = D_diag.NX, D_diag.NY, D_diag.NZ, D_diag.NT -num_colors(D_diag::WilsonEODiagonal) = D_diag.NC - struct WilsonEOPreFermionAction{Nf,C,TD,CT,TX,RI1,RI2,RT} <: AbstractFermionAction{Nf} D::TD cg_temps::CT @@ -178,7 +148,9 @@ struct WilsonEOPreFermionAction{Nf,C,TD,CT,TX,RI1,RI2,RT} <: AbstractFermionActi rhmc_temps2 = nothing cg_temps = ntuple(_ -> even_odd(Spinorfield(f)), 4) else - @assert Nf == 1 "Nf should be 1 or 2 (was $Nf). If you want Nf > 2, use multiple actions" + @assert Nf == 1 """ + Nf should be 1 or 2 (was $Nf). If you want Nf > 2, use multiple actions + """ cg_temps = ntuple(_ -> even_odd(Spinorfield(f)), 2) power = Nf//4 rhmc_info_action = RHMCParams( @@ -229,7 +201,7 @@ function calc_fermion_action( # clear!(ψ_eo) # initial guess is zero # solve_dirac!(ψ_eo, DdagD, ϕ_eo, temp1, temp2, temp3, cg_tol, cg_maxiters) # ψ = (D†D)⁻¹ϕ - Sf = #= dot(ϕ_eo, ψ_eo) =# - 2trlog(D.D_diag) + Sf = #= dot(ϕ_eo, ψ_eo) =# - 2trlog(D.D_diag, D.mass) return distributed_reduce(real(Sf), +, U) end @@ -286,8 +258,8 @@ function LinearAlgebra.mul!( end function LinearAlgebra.mul!( - ψ_eo::TF, D::Daggered{WilsonEOPreDiracOperator{CPU,T,C,TF,TG,TX,TO}}, ϕ_eo::TF -) where {T,C,TF,TG,TX,TO} + ψ_eo::TF, D::Daggered{WilsonEOPreDiracOperator{CPU,T,C,TF,TG,TX,TO,BC}}, ϕ_eo::TF +) where {T,C,TF,TG,TX,TO,BC} @assert TG !== Nothing "Dirac operator has no gauge background, do `D(U)`" U = D.parent.U check_dims(ψ_eo, ϕ_eo, U) @@ -303,8 +275,8 @@ function LinearAlgebra.mul!( end function LinearAlgebra.mul!( - ψ_eo::TF, D::DdaggerD{WilsonEOPreDiracOperator{CPU,T,C,TF,TG,TX,TO}}, ϕ_eo::TF -) where {T,C,TF,TG,TX,TO} + ψ_eo::TF, D::DdaggerD{WilsonEOPreDiracOperator{CPU,T,C,TF,TG,TX,TO,BC}}, ϕ_eo::TF +) where {T,C,TF,TG,TX,TO,BC} temp = D.parent.temp mul!(temp, D.parent, ϕ_eo) # temp = Dϕ mul!(ψ_eo, adjoint(D.parent), temp) # ψ = D†Dϕ @@ -386,14 +358,18 @@ function wilson_eo_kernel(U, ϕ, site, bc, ::Type{T}, ::Val{dagg}) where {T,dagg _siteμ⁺ = eo_site(move(site, 4, 1, NT), NX, NY, NZ, NT, NV) siteμ⁻ = move(site, 4, -1, NT) _siteμ⁻ = eo_site(siteμ⁻, NX, NY, NZ, NT, NV) - ψₙ += cmvmul_spin_proj(U[4, site], apply_bc(ϕ[_siteμ⁺]), Val(-4dagg), Val(false)) - ψₙ += cmvmul_spin_proj(U[4, siteμ⁻], apply_bc(ϕ[_siteμ⁻]), Val(4dagg), Val(true)) + ψₙ += cmvmul_spin_proj( + U[4, site], apply_bc(ϕ[_siteμ⁺], bc, site, Val(1), NT), Val(-4dagg), Val(false) + ) + ψₙ += cmvmul_spin_proj( + U[4, siteμ⁻], apply_bc(ϕ[_siteμ⁻], bc, site, Val(-1), NT), Val(4dagg), Val(true) + ) return T(0.5) * ψₙ end function calc_diag!( D_diag::TW, D_oo_inv::TW, ::Nothing, U::Gaugefield{CPU,T}, mass -) where {T,TW<:WilsonEODiagonal{CPU,T,false}} +) where {T,M,TW<:Paulifield{CPU,T,M,false}} check_dims(D_diag, D_oo_inv, U) mass_term = Complex{T}(4 + mass) fdims = dims(U) @@ -402,29 +378,26 @@ function calc_diag!( #= @batch =#for site in eachindex(U) _site = eo_site(site, fdims..., NV) A = SMatrix{6,6,Complex{T},36}(mass_term * I) - - D_diag[1, site] = A - D_diag[2, site] = A + D_diag[site] = PauliMatrix(A, A) if isodd(site) o_site = switch_sides(_site, fdims..., NV) A_inv = SMatrix{6,6,Complex{T},36}(1/mass_term * I) - D_oo_inv[1, o_site] = A_inv - D_oo_inv[2, o_site] = A_inv + D_oo_inv[o_site] = PauliMatrix(A_inv, A_inv) end end end function calc_diag!( D_diag::TW, D_oo_inv::TW, Fμν, U::Gaugefield{CPU,T}, mass -) where {T,TW<:WilsonEODiagonal{CPU,T,true}} +) where {T,MP,TW<:Paulifield{CPU,T,MP,true}} # With clover term check_dims(D_diag, D_oo_inv, U) mass_term = Complex{T}(4 + mass) fdims = dims(U) NV = U.NV - fac = T(-D_diag.csw / 2) + fac = Complex{T}(-D_diag.csw / 2) - fieldstrength_A_eachsite!(Clover(), Fμν, U) + fieldstrength_eachsite!(Clover(), Fμν, U) #= @batch =#for site in eachindex(U) _site = eo_site(site, fdims..., NV) @@ -433,50 +406,48 @@ function calc_diag!( j = SVector((3, 4)) F₁₂ = Fμν[1, 2, site] - σ = σ₁₂(T) + σ = σ12(T) A₊ = ckron(σ[i, i], F₁₂) A₋ = ckron(σ[j, j], F₁₂) F₁₃ = Fμν[1, 3, site] - σ = σ₁₃(T) + σ = σ13(T) A₊ += ckron(σ[i, i], F₁₃) A₋ += ckron(σ[j, j], F₁₃) F₁₄ = Fμν[1, 4, site] - σ = σ₁₄(T) + σ = σ14(T) A₊ += ckron(σ[i, i], F₁₄) A₋ += ckron(σ[j, j], F₁₄) F₂₃ = Fμν[2, 3, site] - σ = σ₂₃(T) + σ = σ23(T) A₊ += ckron(σ[i, i], F₂₃) A₋ += ckron(σ[j, j], F₂₃) F₂₄ = Fμν[2, 4, site] - σ = σ₂₄(T) + σ = σ24(T) A₊ += ckron(σ[i, i], F₂₄) A₋ += ckron(σ[j, j], F₂₄) F₃₄ = Fμν[3, 4, site] - σ = σ₃₄(T) + σ = σ34(T) A₊ += ckron(σ[i, i], F₃₄) A₋ += ckron(σ[j, j], F₃₄) A₊ = fac * A₊ + M A₋ = fac * A₋ + M - D_diag[1, _site] = A₊ - D_diag[2, _site] = A₋ + D_diag[site] = PauliMatrix(A₊, A₋) # XXX: - if isodd(site) - o_site = switch_sides(_site, fdims..., NV) - D_oo_inv[1, o_site] = cinv(A₊) - D_oo_inv[2, o_site] = cinv(A₋) - end + # if isodd(site) + # o_site = switch_sides(_site, fdims..., NV) + D_oo_inv[site] = PauliMatrix(cinv(A₊), cinv(A₋)) # XXX: + # end end end function mul_oo_inv!( - ϕ_eo::WilsonEOPreSpinorfield{CPU,T}, D_diag::WilsonEODiagonal{CPU,T} + ϕ_eo::WilsonEOPreSpinorfield{CPU,T}, D_diag::Paulifield{CPU,T} ) where {T} check_dims(ϕ_eo, D_diag) ϕ = ϕ_eo.parent @@ -485,14 +456,14 @@ function mul_oo_inv!( #= @batch =#for _site in eachindex(true, ϕ) o_site = switch_sides(_site, fdims..., NV) - ϕ[o_site] = cmvmul_block(D_diag[1, _site], D_diag[2, _site], ϕ[o_site]) + ϕ[o_site] = cmvmul_block(D_diag[_site], ϕ[o_site]) end return nothing end function axmy!( - D_diag::WilsonEODiagonal{CPU,T}, ψ_eo::TF, ϕ_eo::TF + D_diag::Paulifield{CPU,T}, ψ_eo::TF, ϕ_eo::TF ) where {T,TF<:SpinorfieldEO{CPU,T}} # even on even is the default check_dims(ϕ_eo, ψ_eo) ϕ = ϕ_eo.parent @@ -500,26 +471,26 @@ function axmy!( even = true #= @batch =#for _site in eachindex(even, ϕ) - ϕ[_site] = cmvmul_block(D_diag[1, _site], D_diag[2, _site], ψ[_site]) - ϕ[_site] + ϕ[_site] = cmvmul_block(D_diag[_site], ψ[_site]) - ϕ[_site] end return nothing end -function trlog(D_diag::WilsonEODiagonal{CPU,T,true}) where {T} +function trlog(D_diag::Paulifield{CPU,T,M,false}, mass) where {T,M} # Without clover term + NC = num_colors(D_diag) + mass_term = Float64(4 + mass) + logd = 4NC * log(mass_term) + return D_diag.NV÷2 * logd +end + +function trlog(D_diag::Paulifield{CPU,T,M,true}, ::Any) where {T,M} # With clover term d = 0.0 - @batch reduction=(+, d) for _site in eachindex(true, D_diag) - d += log(real(det(D_diag[1, _site]))) - d += log(real(det(D_diag[2, _site]))) + @batch reduction=(+, d) for site in eachindex(D_diag) # XXX: + p = D_diag[site] # XXX: + d += log(real(det(p.upper)) * real(det(p.lower))) end return d end - -function trlog(D_diag::WilsonEODiagonal{CPU,T,false}) where {T} - NC = num_colors(D_diag) - mass_term = Float64(4 + D_diag.mass) - logd = 4NC * log(mass_term) - return D_diag.NV÷2 * logd -end diff --git a/src/fields/Fields.jl b/src/fields/Fields.jl index dc757846..def14759 100644 --- a/src/fields/Fields.jl +++ b/src/fields/Fields.jl @@ -36,7 +36,9 @@ include("distributed.jl") # utility functions for MPI-distributed fields include("boundaries.jl") # boundary conditions in time direction for spinors include("gaugefield.jl") # Gaugefield, Colorfield and Expfield structs defined here include("algebrafield.jl") # For now just a placeholder in case I want to implement more efficient storage of su(3) algebra elements -include("spinorfield.jl") # Spinorfield and SpinorfieldEO (Even-Odd) structs defined here +include("spinorfield.jl") # Spinorfield structs defined here +include("spinorfield_eo.jl") # Spinorfield for even-odd precon +include("paulifield.jl") # For now just a placeholder in case I want to implement more efficient storage of su(3) algebra elements include("tensorfield.jl") # Tensorfield struct and fieldstrength methods defined here include("iterators.jl") # Sequential and Checkerboard iterators defined here include("gpu_iterators.jl") # GPU version of the above @@ -51,6 +53,7 @@ include("wilsonloop.jl") # Definition of arbitrary side length Wilson loops include("gpu_kernels/action.jl") # GPU versions of the above: include("gpu_kernels/algebrafield.jl") include("gpu_kernels/field_operations.jl") +# TODO: include("gpu_kernels/paulifield.jl") include("gpu_kernels/spinorfield.jl") include("gpu_kernels/tensorfield.jl") include("gpu_kernels/wilsonloop.jl") @@ -131,10 +134,12 @@ Check if all fields have the same dimensions. Throw an `AssertionError` otherwis """ @generated function check_dims(x1, rest::Vararg{Any,N}) where {N} q_inner = Expr(:comparison, :(global_dims(x1))) + for i in 1:N push!(q_inner.args, :(==)) push!(q_inner.args, :(global_dims(rest[$i]))) end + q = Expr(:macrocall, Symbol("@assert"), :(), q_inner) return q end @@ -146,8 +151,6 @@ end Base.eachindex(::IndexLinear, u::AbstractMPIField) = error("MPI parallelized field can not be iterated over linearly") -@inline allindices(u::AbstractField) = eachindex(u.U) # all indices including halo regions - @inline function Base.eachindex(even::Bool, u::AbstractField) NX, NY, NZ, NT = global_dims(u) @assert iseven(NT) @@ -155,6 +158,8 @@ Base.eachindex(::IndexLinear, u::AbstractMPIField) = return CartesianIndices((NX, NY, NZ, last_range)) end +@inline allindices(u::AbstractField) = eachindex(u.U) # all indices including halo regions + Base.length(u::AbstractField) = u.NV # overload get and set for the Abstractfields structs, so we dont have to do u.U[μ,x,y,z,t]: diff --git a/src/fields/algebrafield.jl b/src/fields/algebrafield.jl index e77248c2..82edaf20 100644 --- a/src/fields/algebrafield.jl +++ b/src/fields/algebrafield.jl @@ -20,7 +20,9 @@ struct Algebrafield{Backend,FloatType,IsDistributed,ArrayType} <: return new{Backend,FloatType,false,typeof(U)}(U, NX, NY, NZ, NT, NV, 3, topology) end - function Algebrafield{Backend,FloatType}(NX, NY, NZ, NT, numprocs_cart, halo_width) where {Backend,FloatType} + function Algebrafield{Backend,FloatType}( + NX, NY, NZ, NT, numprocs_cart, halo_width + ) where {Backend,FloatType} if prod(numprocs_cart) == 1 return Colorfield{Backend,FloatType}(NX, NY, NZ, NT) end diff --git a/src/fields/distributed.jl b/src/fields/distributed.jl index cafe9c37..b7470872 100644 --- a/src/fields/distributed.jl +++ b/src/fields/distributed.jl @@ -28,8 +28,12 @@ struct FieldTopology global_volume::Int64 # Number of sites in global field local_volume::Int64 # Number of sites in local partition function FieldTopology(numprocs_cart, halo_width, global_dims) - @assert global_dims .% numprocs_cart == (0, 0, 0, 0) "Lattice size must be divisible by number of processes per dimension" - @assert minimum(global_dims ./ numprocs_cart) >= halo_width "Halo must not be wider than the bulk" + @assert global_dims .% numprocs_cart == (0, 0, 0, 0) """ + Lattice size must be divisible by number of processes per dimension + """ + @assert minimum(global_dims ./ numprocs_cart) >= halo_width """ + Halo must not be wider than the bulk + """ comm_cart = mpi_cart_create(numprocs_cart; periodic=map(_->true, numprocs_cart)) numprocs = prod(numprocs_cart) diff --git a/src/fields/gaugefield.jl b/src/fields/gaugefield.jl index b68bfac1..cb0d3972 100644 --- a/src/fields/gaugefield.jl +++ b/src/fields/gaugefield.jl @@ -1,5 +1,6 @@ """ - Gaugefield{Backend,FloatType,IsDistributed,ArrayType,GaugeAction} <: AbstractField{Backend,FloatType,IsDistributed,ArrayType} + Gaugefield{Backend,FloatType,IsDistributed,ArrayType,GaugeAction} <: + AbstractField{Backend,FloatType,IsDistributed,ArrayType} 5-dimensional dense array of statically sized 3x3 matrices contatining associated meta-data. @@ -58,7 +59,10 @@ struct Gaugefield{Backend,FloatType,IsDistributed,ArrayType,GaugeAction} <: return Gaugefield{Backend,FloatType,GaugeAction}(NX, NY, NZ, NT, β) end - @assert halo_width >= stencil_size(GaugeAction) "halo_width must be >= 2 when using improved gauge actions" + @assert halo_width >= stencil_size(GaugeAction) """ + halo_width must be >= 2 when using improved gauge actions + """ + NV = NX * NY * NZ * NT topology = FieldTopology(numprocs_cart, halo_width, (NX, NY, NZ, NT)) ldims = topology.local_dims @@ -78,7 +82,9 @@ function Gaugefield( u_out = if IsDistributed numprocs_cart = u.topology.numprocs_cart halo_width = u.topology.halo_width - Gaugefield{Backend,FloatType,GaugeAction}(u.NX, u.NY, u.NZ, u.NT, u.β, numprocs_cart, halo_width) + Gaugefield{Backend,FloatType,GaugeAction}( + u.NX, u.NY, u.NZ, u.NT, u.β, numprocs_cart, halo_width + ) else Gaugefield{Backend,FloatType,GaugeAction}(u.NX, u.NY, u.NZ, u.NT, u.β) end @@ -97,7 +103,9 @@ function Gaugefield(parameters) halo_width = parameters.halo_width U = if numprocs > 1 - Gaugefield{Backend,FloatType,GaugeAction}(NX, NY, NZ, NT, β, numprocs_cart, halo_width) + Gaugefield{Backend,FloatType,GaugeAction}( + NX, NY, NZ, NT, β, numprocs_cart, halo_width + ) else Gaugefield{Backend,FloatType,GaugeAction}(NX, NY, NZ, NT, β) end @@ -115,7 +123,8 @@ function Gaugefield(parameters) end """ - Colorfield{Backend,FloatType,IsDistributed,ArrayType} <: AbstractField{Backend,FloatType,IsDistributed,ArrayType} + Colorfield{Backend,FloatType,IsDistributed,ArrayType} <: + AbstractField{Backend,FloatType,IsDistributed,ArrayType} 5-dimensional dense array of statically sized 3x3 matrices contatining associated meta-data. @@ -181,7 +190,8 @@ function Colorfield( end """ - Expfield{Backend,FloatType,IsDistributed,ArrayType} <: AbstractField{Backend,FloatType,IsDistributed,ArrayType} + Expfield{Backend,FloatType,IsDistributed,ArrayType} <: + AbstractField{Backend,FloatType,IsDistributed,ArrayType} 5-dimensional dense array of `exp_iQ_su3` objects contatining associated meta-data. The objects hold the `Q`-matrices and all the exponential parameters needed for stout-force @@ -218,7 +228,9 @@ struct Expfield{Backend,FloatType,IsDistributed,ArrayType} <: return new{Backend,FloatType,false,typeof(U)}(U, NX, NY, NZ, NT, NV, 3, topology) end - function Expfield{Backend,FloatType}(NX, NY, NZ, NT, numprocs_cart, halo_width) where {Backend,FloatType} + function Expfield{Backend,FloatType}( + NX, NY, NZ, NT, numprocs_cart, halo_width + ) where {Backend,FloatType} if prod(numprocs_cart) == 1 return Expfield{Backend,FloatType}(NX, NY, NZ, NT) end diff --git a/src/fields/paulifield.jl b/src/fields/paulifield.jl new file mode 100644 index 00000000..d7f16d2e --- /dev/null +++ b/src/fields/paulifield.jl @@ -0,0 +1,40 @@ +struct Paulifield{B,T,M,C,A} <: AbstractField{B,T,M,A} # So we can overload LinearAlgebra.det on the even-odd diagonal + U::A + csw::Float64 + NX::Int64 # Number of lattice sites in the x-direction + NY::Int64 # Number of lattice sites in the y-direction + NZ::Int64 # Number of lattice sites in the z-direction + NT::Int64 # Number of lattice sites in the t-direction + NV::Int64 # Total number of lattice sites + NC::Int64 # Number of colors + ND::Int64 # Number of dirac indices + + topology::FieldTopology # Info regarding MPI topology + function Paulifield( + f::TF, csw; inverse=false + ) where {B,T,TF<:Union{Spinorfield{B,T},SpinorfieldEO{B,T}}} + if f isa SpinorfieldEO + numprocs_cart = f.parent.topology.numprocs_cart + halo_width = f.parent.topology.halo_width + else + numprocs_cart = f.topology.numprocs_cart + halo_width = f.topology.halo_width + end + + NX, NY, NZ, NT = dims(f) + _NT = inverse ? NT ÷ 2 : NT + topology = FieldTopology(numprocs_cart, halo_width, (NX, NY, NZ, _NT)) + NC = num_colors(f) + ND = num_dirac(f) + NV = NX * NY * NZ * NT + C = csw == 0 ? false : true + U = KA.zeros(B(), PauliMatrix{T}, NX, NY, NZ, _NT) + + A = typeof(U) + return new{B,T,false,C,A}(U, csw, NX, NY, NZ, NT, NV, NC, ND, topology) + end +end + +@inline global_dims(p::Paulifield) = NTuple{4,Int64}((p.NX, p.NY, p.NZ, p.NT)) +num_colors(p::Paulifield) = p.NC + diff --git a/src/fields/spinorfield.jl b/src/fields/spinorfield.jl index 9ffc38ed..615ebb02 100644 --- a/src/fields/spinorfield.jl +++ b/src/fields/spinorfield.jl @@ -202,8 +202,6 @@ function LinearAlgebra.axpby!(α, ψ::T, β, ϕ::T) where {T<:Spinorfield{CPU}} ϕ[site] = α * ψ[site] + β * ϕ[site] end - # INFO: don't need to do halo exchange here, since we iterate over all indices - # including halo regions return nothing end @@ -219,243 +217,3 @@ function LinearAlgebra.dot(ϕ::T, ψ::T) where {T<:Spinorfield{CPU}} return distributed_reduce(res, +, ϕ) end - -""" - even_odd(f::Spinorfield) - -Create a wrapper around a `Spinorfield` to signal that it is meant to be used in the context -of even-odd preconditioning. What this amounts to is that we realign the entries such that -`ϕ -> (ϕₑ, ϕₒ)`, which is achieved by recalculating the index whenever we index into `ϕ` -or iterating only over one half of its indices. -""" -even_odd(f::Spinorfield) = SpinorfieldEO(f) - -struct SpinorfieldEO{Backend,FloatType,IsDistributed,ArrayType,NumDirac} <: - AbstractField{Backend,FloatType,IsDistributed,ArrayType} - parent::Spinorfield{Backend,FloatType,IsDistributed,ArrayType,NumDirac} - function SpinorfieldEO( - f::Spinorfield{Backend,FloatType,IsDistributed,ArrayType,NumDirac} - ) where {Backend,FloatType,IsDistributed,ArrayType,NumDirac} - @assert iseven(f.NT) "Need even time extent for even-odd preconditioning" - return new{Backend,FloatType,IsDistributed,ArrayType,NumDirac}(f) - end -end - -function Spinorfield( - f::SpinorfieldEO{Backend,FloatType,IsDistributed,ArrayType,NumDirac} -) where {Backend,FloatType,IsDistributed,ArrayType,NumDirac} - return SpinorfieldEO(f.parent) -end - -dims(f::SpinorfieldEO) = dims(f.parent) -local_dims(f::SpinorfieldEO) = local_dims(f.parent) -global_dims(f::SpinorfieldEO) = global_dims(f.parent) -Base.size(f::SpinorfieldEO) = size(f.parent) -Base.similar(f::SpinorfieldEO) = even_odd(Spinorfield(f.parent)) -Base.eltype(::SpinorfieldEO{B,T}) where {B,T} = Complex{T} -LinearAlgebra.checksquare(f::SpinorfieldEO) = LinearAlgebra.checksquare(f.parent) ÷ 2 -num_colors(::SpinorfieldEO{B,T,M,A,ND}) where {B,T,M,A,ND} = 3 -num_dirac(::SpinorfieldEO{B,T,M,A,ND}) where {B,T,M,A,ND} = ND -volume(f::SpinorfieldEO) = volume(f.parent) - -Base.@propagate_inbounds Base.getindex(f::SpinorfieldEO, i::Integer) = f.parent.U[i] -Base.@propagate_inbounds Base.getindex(f::SpinorfieldEO, x, y, z, t) = f.parent.U[x, y, z, t] -Base.@propagate_inbounds Base.getindex(f::SpinorfieldEO, site::SiteCoords) = f.parent.U[site] -Base.@propagate_inbounds Base.setindex!(f::SpinorfieldEO, v, i::Integer) = - setindex!(f.parent.U, v, i) -Base.@propagate_inbounds Base.setindex!(f::SpinorfieldEO, v, x, y, z, t) = - setindex!(f.parent.U, v, x, y, z, t) -Base.@propagate_inbounds Base.setindex!(f::SpinorfieldEO, v, site::SiteCoords) = - setindex!(f.parent.U, v, site) - -Base.view(f::SpinorfieldEO, I::CartesianIndices{4}) = view(f.parent.U, I.indices...) - -clear!(ϕ_eo::SpinorfieldEO) = clear!(ϕ_eo.parent) -ones!(ϕ_eo::SpinorfieldEO) = ones!(ϕ_eo.parent) - -function set_source!(ϕ_eo::SpinorfieldEO{CPU,T}, site::SiteCoords, a, μ) where {T} - ϕ = ϕ_eo.parent - NC = num_colors(ϕ) - ND = num_dirac(ϕ) - @assert μ ∈ 1:ND && a ∈ 1:3 - clear!(ϕ) - vec_index = (μ - 1) * NC + a - tup = ntuple(i -> i == vec_index ? one(Complex{T}) : zero(Complex{T}), Val(3ND)) - _site = eo_site(site, global_dims(ϕ)..., ϕ.NV) - ϕ[_site] = SVector{3ND,Complex{T}}(tup) - update_halo!(ϕ) # TODO: Even-odd halo exchange - return nothing -end - -function Base.copy!(ϕ_eo::TF, ψ_eo::TF) where {TF<:SpinorfieldEO{CPU}} - check_dims(ϕ_eo, ψ_eo) - ϕ = ϕ_eo.parent - ψ = ψ_eo.parent - even = true - - @batch for e_site in eachindex(even, ϕ) - ϕ[e_site] = ψ[e_site] - end - - update_halo!(ϕ) # TODO: Even-odd halo exchange - return nothing -end - -function gaussian_pseudofermions!(ϕ_eo::SpinorfieldEO{CPU,T}) where {T} - ϕ = ϕ_eo.parent - sz = num_dirac(ϕ) * num_colors(ϕ) - even = true - - for e_site in eachindex(even, ϕ) - ϕ[e_site] = @SVector randn(Complex{T}, sz) # σ = 0.5 - end - - update_halo!(ϕ) # TODO: Even-odd halo exchange - return nothing -end - -function LinearAlgebra.mul!(ϕ_eo::SpinorfieldEO{CPU,T}, α) where {T} - ϕ = ϕ_eo.parent - α = Complex{T}(α) - even = true - - @batch for _site in eachindex(even, ϕ) - ϕ[_site] *= α - end - - update_halo!(ϕ) # TODO: Even-odd halo exchange - return nothing -end - -function LinearAlgebra.axpy!(α, ψ_eo::T, ϕ_eo::T) where {T<:SpinorfieldEO{CPU}} # even on even is the default - check_dims(ϕ_eo, ψ_eo) - ϕ = ϕ_eo.parent - ψ = ψ_eo.parent - FloatT = float_type(ϕ) - α = Complex{FloatT}(α) - even = true - - @batch for _site in eachindex(even, ϕ) - ϕ[_site] += α * ψ[_site] - end - - update_halo!(ϕ) # TODO: Even-odd halo exchange - return nothing -end - -function LinearAlgebra.axpby!(α, ψ_eo::T, β, ϕ_eo::T, even=true) where {T<:SpinorfieldEO{CPU}} - check_dims(ϕ_eo, ψ_eo) - ϕ = ϕ_eo.parent - ψ = ψ_eo.parent - FloatT = float_type(ϕ) - α = Complex{FloatT}(α) - β = Complex{FloatT}(β) - - @batch for _site in eachindex(even, ϕ) - ϕ[_site] = α * ψ[_site] + β * ϕ[_site] - end - - update_halo!(ϕ) # TODO: Even-odd halo exchange - return nothing -end - -LinearAlgebra.norm(ϕ_eo::SpinorfieldEO) = sqrt(real(dot(ϕ_eo, ϕ_eo))) - -function LinearAlgebra.dot(ϕ_eo::T, ψ_eo::T) where {T<:SpinorfieldEO{CPU}} - check_dims(ϕ_eo, ψ_eo) - ϕ = ϕ_eo.parent - ψ = ψ_eo.parent - res = 0.0 + 0.0im # res is always double precision, even if T is single precision - even = true - - @batch reduction = (+, res) for _site in eachindex(even, ϕ) - res += cdot(ϕ[_site], ψ[_site]) - end - - return distributed_reduce(res, +, ϕ) -end - -function dot_all(ϕ_eo::T, ψ_eo::T) where {T<:SpinorfieldEO{CPU}} - check_dims(ϕ_eo, ψ_eo) - ϕ = ϕ_eo.parent - ψ = ψ_eo.parent - res = 0.0 + 0.0im # res is always double precision, even if T is single precision - - @batch reduction = (+, res) for site in eachindex(ϕ) - res += cdot(ϕ[site], ψ[site]) - end - - return distributed_reduce(res, +, ϕ) -end - -function copy_eo!(ϕ_eo::T, ψ_eo::T) where {T<:SpinorfieldEO{CPU}} - check_dims(ϕ_eo, ψ_eo) - ϕ = ϕ_eo.parent - ψ = ψ_eo.parent - fdims = dims(ϕ) - NV = ϕ.NV - even = true - - for e_site in eachindex(even, ϕ) - o_site = switch_sides(e_site, fdims..., NV) - ϕ[e_site] = ψ[o_site] - end - - update_halo!(ϕ) # TODO: Even-odd halo exchange - return nothing -end - -function copy_oe!(ϕ_eo::T, ψ_eo::T) where {T<:SpinorfieldEO{CPU}} - check_dims(ϕ_eo, ψ_eo) - ϕ = ϕ_eo.parent - ψ = ψ_eo.parent - fdims = dims(ϕ) - NV = ϕ.NV - odd = false - - for o_site in eachindex(odd, ϕ) - e_site = switch_sides(o_site, fdims..., NV) - ϕ[o_site] = ψ[e_site] - end - - update_halo!(ϕ) # TODO: Even-odd halo exchange - return nothing -end - -function axpy_oe!(α, ψ_eo::T, ϕ_eo::T) where {T<:SpinorfieldEO{CPU}} - check_dims(ϕ_eo, ψ_eo) - ϕ = ϕ_eo.parent - ψ = ψ_eo.parent - FloatT = float_type(ϕ) - α = Complex{FloatT}(α) - fdims = dims(ϕ) - NV = ϕ.NV - even = true - - for e_site in eachindex(even, ϕ) - o_site = switch_sides(e_site, fdims..., NV) - ϕ[e_site] += α * ψ[o_site] - end - - update_halo!(ϕ) # TODO: Even-odd halo exchange - return nothing -end - -function axpy_eo!(α, ψ_eo::T, ϕ_eo::T) where {T<:SpinorfieldEO{CPU}} - check_dims(ϕ_eo, ψ_eo) - ϕ = ϕ_eo.parent - ψ = ψ_eo.parent - FloatT = float_type(ϕ) - α = Complex{FloatT}(α) - fdims = dims(ϕ) - NV = ϕ.NV - odd = false - - for o_site in eachindex(odd, ϕ) - e_site = switch_sides(o_site, fdims..., NV) - ϕ[o_site] += α * ψ[e_site] - end - - update_halo!(ϕ) # TODO: Even-odd halo exchange - return nothing -end diff --git a/src/fields/spinorfield_eo.jl b/src/fields/spinorfield_eo.jl new file mode 100644 index 00000000..7a0f6f80 --- /dev/null +++ b/src/fields/spinorfield_eo.jl @@ -0,0 +1,264 @@ +""" + even_odd(f::Spinorfield) + +Create a wrapper around a `Spinorfield` to signal that it is meant to be used in the context +of even-odd preconditioning. What this amounts to is that we realign the entries such that +`ϕ -> (ϕₑ, ϕₒ)`, which is achieved by recalculating the index whenever we index into `ϕ` +or iterating only over one half of its indices. +""" +even_odd(f::Spinorfield) = SpinorfieldEO(f) + +struct SpinorfieldEO{Backend,FloatType,IsDistributed,ArrayType,NumDirac} <: + AbstractField{Backend,FloatType,IsDistributed,ArrayType} + parent::Spinorfield{Backend,FloatType,IsDistributed,ArrayType,NumDirac} + function SpinorfieldEO( + f::Spinorfield{Backend,FloatType,IsDistributed,ArrayType,NumDirac} + ) where {Backend,FloatType,IsDistributed,ArrayType,NumDirac} + @assert iseven(f.NT) "Need even time extent for even-odd preconditioning" + + if IsDistributed + @assert iseven(f.NT÷2) """ + Need time extent divisbile by 4 for even-odd preconditioning with \ + distributed fields + """ + end + + return new{Backend,FloatType,IsDistributed,ArrayType,NumDirac}(f) + end +end + +function Spinorfield( + f::SpinorfieldEO{Backend,FloatType,IsDistributed,ArrayType,NumDirac} +) where {Backend,FloatType,IsDistributed,ArrayType,NumDirac} + return SpinorfieldEO(f.parent) +end + +dims(f::SpinorfieldEO) = dims(f.parent) +local_dims(f::SpinorfieldEO) = local_dims(f.parent) +global_dims(f::SpinorfieldEO) = global_dims(f.parent) +Base.size(f::SpinorfieldEO) = size(f.parent) +Base.similar(f::SpinorfieldEO) = even_odd(Spinorfield(f.parent)) +Base.eltype(::SpinorfieldEO{B,T}) where {B,T} = Complex{T} +LinearAlgebra.checksquare(f::SpinorfieldEO) = LinearAlgebra.checksquare(f.parent) ÷ 2 +num_colors(::SpinorfieldEO{B,T,M,A,ND}) where {B,T,M,A,ND} = 3 +num_dirac(::SpinorfieldEO{B,T,M,A,ND}) where {B,T,M,A,ND} = ND +volume(f::SpinorfieldEO) = volume(f.parent) + +Base.@propagate_inbounds Base.getindex(f::SpinorfieldEO, i::Integer) = f.parent.U[i] +Base.@propagate_inbounds Base.getindex(f::SpinorfieldEO, x, y, z, t) = f.parent.U[x, y, z, t] +Base.@propagate_inbounds Base.getindex(f::SpinorfieldEO, site::SiteCoords) = f.parent.U[site] +Base.@propagate_inbounds Base.setindex!(f::SpinorfieldEO, v, i::Integer) = + setindex!(f.parent.U, v, i) +Base.@propagate_inbounds Base.setindex!(f::SpinorfieldEO, v, x, y, z, t) = + setindex!(f.parent.U, v, x, y, z, t) +Base.@propagate_inbounds Base.setindex!(f::SpinorfieldEO, v, site::SiteCoords) = + setindex!(f.parent.U, v, site) + +Base.view(f::SpinorfieldEO, I::CartesianIndices{4}) = view(f.parent.U, I.indices...) + +@inline function allindices(even::Bool, u::Union{Spinorfield,SpinorfieldEO}) + NX, NY, NZ, NT = size(u.U) + @assert iseven(NT) + last_range = even ? (1:div(NT, 2)) : (div(NT, 2)+1:NT) + return CartesianIndices((NX, NY, NZ, last_range)) +end + +clear!(ϕ_eo::SpinorfieldEO) = clear!(ϕ_eo.parent) +ones!(ϕ_eo::SpinorfieldEO) = ones!(ϕ_eo.parent) + +function set_source!(ϕ_eo::SpinorfieldEO{CPU,T}, site::SiteCoords, a, μ) where {T} + ϕ = ϕ_eo.parent + NC = num_colors(ϕ) + ND = num_dirac(ϕ) + @assert μ ∈ 1:ND && a ∈ 1:3 + clear!(ϕ) + vec_index = (μ - 1) * NC + a + tup = ntuple(i -> i == vec_index ? one(Complex{T}) : zero(Complex{T}), Val(3ND)) + _site = eo_site(site, global_dims(ϕ)..., ϕ.NV) + ϕ[_site] = SVector{3ND,Complex{T}}(tup) + update_halo!(ϕ) # TODO: Even-odd halo exchange + return nothing +end + +function Base.copy!(ϕ_eo::TF, ψ_eo::TF) where {TF<:SpinorfieldEO{CPU}} + check_dims(ϕ_eo, ψ_eo) + ϕ = ϕ_eo.parent + ψ = ψ_eo.parent + even = true + + @batch for e_site in eachindex(even, ϕ) + ϕ[e_site] = ψ[e_site] + end + + update_halo!(ϕ) # TODO: Even-odd halo exchange + return nothing +end + +function gaussian_pseudofermions!(ϕ_eo::SpinorfieldEO{CPU,T}) where {T} + ϕ = ϕ_eo.parent + sz = num_dirac(ϕ) * num_colors(ϕ) + even = true + + for e_site in eachindex(even, ϕ) + ϕ[e_site] = @SVector randn(Complex{T}, sz) # σ = 0.5 + end + + update_halo!(ϕ) # TODO: Even-odd halo exchange + return nothing +end + +function LinearAlgebra.mul!(ϕ_eo::SpinorfieldEO{CPU,T}, α) where {T} + ϕ = ϕ_eo.parent + α = Complex{T}(α) + even = true + + @batch for _site in eachindex(even, ϕ) + ϕ[_site] *= α + end + + update_halo!(ϕ) # TODO: Even-odd halo exchange + return nothing +end + +function LinearAlgebra.axpy!(α, ψ_eo::T, ϕ_eo::T) where {T<:SpinorfieldEO{CPU}} # even on even is the default + check_dims(ϕ_eo, ψ_eo) + ϕ = ϕ_eo.parent + ψ = ψ_eo.parent + FloatT = float_type(ϕ) + α = Complex{FloatT}(α) + even = true + + @batch for _site in eachindex(even, ϕ) + ϕ[_site] += α * ψ[_site] + end + + update_halo!(ϕ) # TODO: Even-odd halo exchange + return nothing +end + +function LinearAlgebra.axpby!( + α, ψ_eo::T, β, ϕ_eo::T, even=true +) where {T<:SpinorfieldEO{CPU}} + check_dims(ϕ_eo, ψ_eo) + ϕ = ϕ_eo.parent + ψ = ψ_eo.parent + FloatT = float_type(ϕ) + α = Complex{FloatT}(α) + β = Complex{FloatT}(β) + + @batch for _site in allindices(even, ϕ) + ϕ[_site] = α * ψ[_site] + β * ϕ[_site] + end + + # INFO: don't need to do halo exchange here, since we iterate over all indices + # including halo regions + return nothing +end + +LinearAlgebra.norm(ϕ_eo::SpinorfieldEO) = sqrt(real(dot(ϕ_eo, ϕ_eo))) + +function LinearAlgebra.dot(ϕ_eo::T, ψ_eo::T) where {T<:SpinorfieldEO{CPU}} + check_dims(ϕ_eo, ψ_eo) + ϕ = ϕ_eo.parent + ψ = ψ_eo.parent + res = 0.0 + 0.0im # res is always double precision, even if T is single precision + even = true + + @batch reduction = (+, res) for _site in eachindex(even, ϕ) + res += cdot(ϕ[_site], ψ[_site]) + end + + return distributed_reduce(res, +, ϕ) +end + +function dot_all(ϕ_eo::T, ψ_eo::T) where {T<:SpinorfieldEO{CPU}} + check_dims(ϕ_eo, ψ_eo) + ϕ = ϕ_eo.parent + ψ = ψ_eo.parent + res = 0.0 + 0.0im # res is always double precision, even if T is single precision + + @batch reduction = (+, res) for site in eachindex(ϕ) + res += cdot(ϕ[site], ψ[site]) + end + + return distributed_reduce(res, +, ϕ) +end + +function copy_eo!(ϕ_eo::T, ψ_eo::T) where {T<:SpinorfieldEO{CPU}} + check_dims(ϕ_eo, ψ_eo) + ϕ = ϕ_eo.parent + ψ = ψ_eo.parent + fdims = dims(ϕ) + NV = ϕ.NV + even = true + + for e_site in allindices(even, ϕ) + o_site = switch_sides(e_site, fdims..., NV) + ϕ[e_site] = ψ[o_site] + end + + # INFO: don't need to do halo exchange here, since we iterate over all indices + # including halo regions + # We assume that ψ's halo is already up-to-date before calling this + return nothing +end + +function copy_oe!(ϕ_eo::T, ψ_eo::T) where {T<:SpinorfieldEO{CPU}} + check_dims(ϕ_eo, ψ_eo) + ϕ = ϕ_eo.parent + ψ = ψ_eo.parent + fdims = dims(ϕ) + NV = ϕ.NV + odd = false + + for o_site in allindices(odd, ϕ) + e_site = switch_sides(o_site, fdims..., NV) + ϕ[o_site] = ψ[e_site] + end + + # INFO: don't need to do halo exchange here, since we iterate over all indices + # including halo regions + # We assume that ψ's halo is already up-to-date before calling this + return nothing +end + +function axpy_oe!(α, ψ_eo::T, ϕ_eo::T) where {T<:SpinorfieldEO{CPU}} + check_dims(ϕ_eo, ψ_eo) + ϕ = ϕ_eo.parent + ψ = ψ_eo.parent + FloatT = float_type(ϕ) + α = Complex{FloatT}(α) + fdims = dims(ϕ) + NV = ϕ.NV + even = true + + for e_site in allindices(even, ϕ) + o_site = switch_sides(e_site, fdims..., NV) + ϕ[e_site] += α * ψ[o_site] + end + + # INFO: don't need to do halo exchange here, since we iterate over all indices + # including halo regions + return nothing +end + +function axpy_eo!(α, ψ_eo::T, ϕ_eo::T) where {T<:SpinorfieldEO{CPU}} + check_dims(ϕ_eo, ψ_eo) + ϕ = ϕ_eo.parent + ψ = ψ_eo.parent + FloatT = float_type(ϕ) + α = Complex{FloatT}(α) + fdims = dims(ϕ) + NV = ϕ.NV + odd = false + + for o_site in allindices(odd, ϕ) + e_site = switch_sides(o_site, fdims..., NV) + ϕ[o_site] += α * ψ[e_site] + end + + # INFO: don't need to do halo exchange here, since we iterate over all indices + # including halo regions + return nothing +end + diff --git a/src/fields/tensorfield.jl b/src/fields/tensorfield.jl index 09d512a8..619c7768 100644 --- a/src/fields/tensorfield.jl +++ b/src/fields/tensorfield.jl @@ -4,7 +4,25 @@ struct Plaquette <: AbstractFieldstrength end struct Clover <: AbstractFieldstrength end struct Improved <: AbstractFieldstrength end -# TODO: Docs +""" + Tensorfield{Backend,FloatType,IsDistributed,ArrayType} <: + AbstractField{Backend,FloatType,IsDistributed,ArrayType} + +6-dimensional dense array of statically sized 3x3 matrices contatining associated meta-data. + + Tensorfield{Backend,FloatType}(NX, NY, NZ, NT) + Tensorfield{Backend,FloatType}(NX, NY, NZ, NT, numprocs_cart, halo_width) + Tensorfield(u::AbstractField) + Tensorfield(parameters::ParameterSet) + +Creates a Tensorfield on `Backend`, i.e. an array of 3-by-3 `FloatType`-precision matrices +of size `4 x 4 × NX × NY × NZ × NT` or a zero-initialized Tensorfield of the same size as +`u`. +# Supported backends +`CPU` \\ +`CUDABackend` \\ +`ROCBackend` +""" struct Tensorfield{Backend,FloatType,IsDistributed,ArrayType} <: AbstractField{Backend,FloatType,IsDistributed,ArrayType} U::ArrayType # Actual field storing the gauge variables @@ -42,7 +60,7 @@ struct Tensorfield{Backend,FloatType,IsDistributed,ArrayType} <: end function Tensorfield( - u::Gaugefield{Backend,FloatType,IsDistributed,ArrayType} + u::AbstractField{Backend,FloatType,IsDistributed,ArrayType} ) where {Backend,FloatType,IsDistributed,ArrayType} u_out = if IsDistributed numprocs_cart = u.topology.numprocs_cart @@ -79,22 +97,25 @@ function fieldstrength_eachsite!(F::Tensorfield, U, kind_of_fs::String) return nothing end -function fieldstrength_eachsite!(::Plaquette, F::Tensorfield{CPU}, U::Gaugefield{CPU}) +function fieldstrength_eachsite!( + ::Plaquette, F::Tensorfield{CPU,T}, U::Gaugefield{CPU,T} +) where {T} check_dims(F, U) + fac = Complex{T}(im) @batch for site in eachindex(U) C12 = plaquette(U, 1, 2, site) - F[1, 2, site] = im * traceless_antihermitian(C12) + F[1, 2, site] = fac * (C12 - C12') C13 = plaquette(U, 1, 3, site) - F[1, 3, site] = im * traceless_antihermitian(C13) + F[1, 3, site] = fac * (C13 - C13') C14 = plaquette(U, 1, 4, site) - F[1, 4, site] = im * traceless_antihermitian(C14) + F[1, 4, site] = fac * (C14 - C14') C23 = plaquette(U, 2, 3, site) - F[2, 3, site] = im * traceless_antihermitian(C23) + F[2, 3, site] = fac * (C23 - C23') C24 = plaquette(U, 2, 4, site) - F[2, 4, site] = im * traceless_antihermitian(C24) + F[2, 4, site] = fac * (C24 - C24') C34 = plaquette(U, 3, 4, site) - F[3, 4, site] = im * traceless_antihermitian(C34) + F[3, 4, site] = fac * (C34 - C34') end update_halo!(F) @@ -105,46 +126,21 @@ function fieldstrength_eachsite!( ::Clover, F::Tensorfield{CPU,T}, U::Gaugefield{CPU,T} ) where {T} check_dims(F, U) - fac = Complex{T}(im / 4) - - @batch for site in eachindex(U) - C12 = clover_square(U, 1, 2, site, 1) - F[1, 2, site] = fac * traceless_antihermitian(C12) - C13 = clover_square(U, 1, 3, site, 1) - F[1, 3, site] = fac * traceless_antihermitian(C13) - C14 = clover_square(U, 1, 4, site, 1) - F[1, 4, site] = fac * traceless_antihermitian(C14) - C23 = clover_square(U, 2, 3, site, 1) - F[2, 3, site] = fac * traceless_antihermitian(C23) - C24 = clover_square(U, 2, 4, site, 1) - F[2, 4, site] = fac * traceless_antihermitian(C24) - C34 = clover_square(U, 3, 4, site, 1) - F[3, 4, site] = fac * traceless_antihermitian(C34) - end - - update_halo!(F) - return nothing -end - -function fieldstrength_A_eachsite!( - ::Clover, F::Tensorfield{CPU,T}, U::Gaugefield{CPU,T} -) where {T} - check_dims(F, U) - fac = Complex{T}(im / 4) + fac = Complex{T}(im / 8) @batch for site in eachindex(U) C12 = clover_square(U, 1, 2, site, 1) - F[1, 2, site] = fac * antihermitian(C12) + F[1, 2, site] = fac * (C12 - C12') C13 = clover_square(U, 1, 3, site, 1) - F[1, 3, site] = fac * antihermitian(C13) + F[1, 3, site] = fac * (C13 - C13') C14 = clover_square(U, 1, 4, site, 1) - F[1, 4, site] = fac * antihermitian(C14) + F[1, 4, site] = fac * (C14 - C14') C23 = clover_square(U, 2, 3, site, 1) - F[2, 3, site] = fac * antihermitian(C23) + F[2, 3, site] = fac * (C23 - C23') C24 = clover_square(U, 2, 4, site, 1) - F[2, 4, site] = fac * antihermitian(C24) + F[2, 4, site] = fac * (C24 - C24') C34 = clover_square(U, 3, 4, site, 1) - F[3, 4, site] = fac * antihermitian(C34) + F[3, 4, site] = fac * (C34 - C34') end update_halo!(F) diff --git a/src/forces/bias_force.jl b/src/forces/bias_force.jl index 44cf96dd..7c78a990 100644 --- a/src/forces/bias_force.jl +++ b/src/forces/bias_force.jl @@ -1,5 +1,5 @@ """ - calc_dVdU_bare!(dU::Colorfield, F::Tensorfield, U::Gaugefield, temp_force, bias::Bias, is_smeared::Bool) + calc_dVdU_bare!(dU::Colorfield, F::Tensorfield, U, temp_force, bias, is_smeared) Calculate the bias force for a molecular dynamics step, i.e. the derivative of the bias potential w.r.t. the bare/unsmeared field U. diff --git a/src/forces/forces.jl b/src/forces/forces.jl index ad40027a..ca2181e0 100644 --- a/src/forces/forces.jl +++ b/src/forces/forces.jl @@ -6,10 +6,10 @@ import ..DiracOperators: StaggeredDiracOperator, StaggeredFermionAction import ..DiracOperators: StaggeredEOPreDiracOperator, StaggeredEOPreFermionAction import ..DiracOperators: WilsonDiracOperator, WilsonFermionAction, has_clover_term import ..DiracOperators: WilsonEOPreDiracOperator, WilsonEOPreFermionAction -import ..DiracOperators: WilsonEODiagonal import ..DiracOperators: Daggered, DdaggerD, SpinorfieldEO import ..DiracOperators: apply_bc, staggered_η, solve_dirac!, solve_dirac_multishift! import ..DiracOperators: mul_oe!, mul_eo!, mul_oo_inv! +import ..Fields: Paulifield """ calc_dSfdU_bare!(dU::Colorfield, fermion_action, U, ϕ, ::Any, ::NoSmearing) diff --git a/src/forces/overlap_eo_force.jl b/src/forces/overlap_eo_force.jl new file mode 100644 index 00000000..e69de29b diff --git a/src/forces/overlap_force.jl b/src/forces/overlap_force.jl new file mode 100644 index 00000000..e69de29b diff --git a/src/forces/wilson_eo_force.jl b/src/forces/wilson_eo_force.jl index d4ea0ff8..f25692ff 100644 --- a/src/forces/wilson_eo_force.jl +++ b/src/forces/wilson_eo_force.jl @@ -15,8 +15,8 @@ function calc_dSfdU!( # solve_dirac!(X_eo, DdagD, ϕ_eo, Y_eo, temp1, temp2, cg_tol, cg_maxiters) # Y is used here merely as a temp LinearAlgebra.mul!(Y, D, X) # Need to prefix with LinearAlgebra to avoid ambiguity with Gaugefields.mul! # # LinearAlgebra.mul!(Y_eo, D, X_eo) - # mul_oe!(X_eo, U, X_eo, bc, true, Val(1)) # Need to prefix with LinearAlgebra to avoid ambiguity with Gaugefields.mul! - # mul_oe!(Y_eo, U, Y_eo, bc, true, Val(-1)) # Need to prefix with LinearAlgebra to avoid ambiguity with Gaugefields.mul! + # mul_oe!(X_eo, U, X_eo, bc, true, Val(1)) + # mul_oe!(Y_eo, U, Y_eo, bc, true, Val(-1)) # mul_oo_inv!(X_eo, D.D_oo_inv) # mul_oo_inv!(Y_eo, D.D_oo_inv) # add_wilson_eo_derivative!(dU, U, X_eo, Y_eo, bc) @@ -26,8 +26,8 @@ function calc_dSfdU!( D_oo_inv = D.D_oo_inv # calc_Xμν_eo_eachsite!(Xμν, X_eo, Y_eo) # add_clover_derivative!(dU, U, Xμν, D.csw) - calc_Xμν_small_eachsite!(Xμν, D_oo_inv) - add_clover_derivative!(dU, U, Xμν, 2D.csw) + calc_small_Xμν_eachsite!(Xμν, D_oo_inv) + add_small_det_derivative!(dU, U, Xμν, -2D.csw) end return nothing @@ -58,8 +58,8 @@ function calc_dSfdU!( for i in 1:n LinearAlgebra.mul!(Ys[i+1], D, Xs[i+1]) # Need to prefix with LinearAlgebra to avoid ambiguity with Gaugefields.mul! - mul_oe!(Xs[i+1], U, Xs[i+1], bc, true, Val(1)) # Need to prefix with LinearAlgebra to avoid ambiguity with Gaugefields.mul! - mul_oe!(Ys[i+1], U, Ys[i+1], bc, true, Val(-1)) # Need to prefix with LinearAlgebra to avoid ambiguity with Gaugefields.mul! + mul_oe!(Xs[i+1], U, Xs[i+1], bc, true, Val(1)) + mul_oe!(Ys[i+1], U, Ys[i+1], bc, true, Val(-1)) mul_oo_inv!(Xs[i+1], D.D_oo_inv) mul_oo_inv!(Ys[i+1], D.D_oo_inv) add_wilson_derivative!(dU, U, Xs[i+1], Ys[i+1], bc; coeff=coeffs[i]) @@ -74,7 +74,7 @@ function calc_dSfdU!( # TODO: if has_clover_term(D) Xμν = D.Xμν - calc_Xμν_small_eachsite!(Xμν, D_oo_inv) + calc_Xμν_small_eachsite!(Xμν, D.D_oo_inv) add_clover_derivative!(dU, U, Xμν, 2D.csw) end @@ -99,7 +99,7 @@ function add_wilson_eo_derivative!( return nothing end -function add_wilson_eo_derivative_kernel!(dU, U, X_eo, Y_eo, site, bc⁺, fac) +function add_wilson_eo_derivative_kernel!(dU, U, X_eo, Y_eo, site, bc, fac) # sites that begin with a "_" are meant for indexing into the even-odd preconn'ed # fermion field NX, NY, NZ, NT = dims(U) @@ -122,9 +122,13 @@ function add_wilson_eo_derivative_kernel!(dU, U, X_eo, Y_eo, site, bc⁺, fac) dU[3i32, site] += fac * traceless_antihermitian(cmatmul_oo(U[3, site], B + C)) _siteμ⁺ = eo_site(move(site, 4, 1, NT), NX, NY, NZ, NT, NV) - B = spintrace(spin_proj(X_eo[_siteμ⁺], Val(-4)), Y_eo[_site]) - C = spintrace(spin_proj(Y_eo[_siteμ⁺], Val(4)), X_eo[_site]) - dU[4i32, site] += bc⁺ * fac * traceless_antihermitian(cmatmul_oo(U[4, site], B + C)) + B = spintrace( + spin_proj(apply_bc(X_eo[_siteμ⁺], bc, site, Val(1), NT), Val(-4)), Y_eo[_site] + ) + C = spintrace( + spin_proj(apply_bc(Y_eo[_siteμ⁺], bc, site, Val(1), NT), Val(4)), X_eo[_site] + ) + dU[4i32, site] += fac * traceless_antihermitian(cmatmul_oo(U[4, site], B + C)) return nothing end @@ -184,62 +188,102 @@ function calc_Xμν_eo_kernel!(Xμν, X, Y, site) return nothing end -function calc_Xμν_small_eachsite!( - Xμν::Tensorfield{CPU,T}, D_oo_inv::WilsonEODiagonal{CPU,T,true} -) where {T} +function calc_small_Xμν_eachsite!( + Xμν::Tensorfield{CPU,T}, D_oo_inv::Paulifield{CPU,T,M,true} +) where {T,M} check_dims(Xμν, D_oo_inv) for site in eachindex(Xμν) - calc_Xμν_small_kernel!(Xμν, D_oo_inv, site, T) + calc_small_Xμν_kernel!(Xμν, D_oo_inv, site, T) end return nothing end -function calc_Xμν_small_kernel!(Xμν, D_oo_inv, site, ::Type{T}) where {T} - if isodd(site) +function calc_small_Xμν_kernel!(Xμν, D_oo_inv, site, ::Type{T}) where {T} + # if isodd(site) NX, NY, NZ, NT = dims(Xμν) NV = NX * NY * NZ * NT - _site = eo_site_switch(site, NX, NY, NZ, NT, NV) - A₊ = D_oo_inv[1, _site] - A₋ = D_oo_inv[2, _site] + # _site = eo_site(site, NX, NY, NZ, NT, NV) + Minv = D_oo_inv[site] # XXX: - X₁₂ = im * spintrace_σμν(A₊, A₋, Val(1), Val(2)) + X₁₂ = spintrace_σμν(Minv, Val(1), Val(2)) Xμν[1i32, 2i32, site] = X₁₂ Xμν[2i32, 1i32, site] = -X₁₂ - X₁₃ = im * spintrace_σμν(A₊, A₋, Val(1), Val(3)) + X₁₃ = spintrace_σμν(Minv, Val(1), Val(3)) Xμν[1i32, 3i32, site] = X₁₃ Xμν[3i32, 1i32, site] = -X₁₃ - X₁₄ = im * spintrace_σμν(A₊, A₋, Val(1), Val(4)) + X₁₄ = spintrace_σμν(Minv, Val(1), Val(4)) Xμν[1i32, 4i32, site] = X₁₄ Xμν[4i32, 1i32, site] = -X₁₄ - X₂₃ = im * spintrace_σμν(A₊, A₋, Val(2), Val(3)) + X₂₃ = spintrace_σμν(Minv, Val(2), Val(3)) Xμν[2i32, 3i32, site] = X₂₃ Xμν[3i32, 2i32, site] = -X₂₃ - X₂₄ = im * spintrace_σμν(A₊, A₋, Val(2), Val(4)) + X₂₄ = spintrace_σμν(Minv, Val(2), Val(4)) Xμν[2i32, 4i32, site] = X₂₄ Xμν[4i32, 2i32, site] = -X₂₄ - X₃₄ = im * spintrace_σμν(A₊, A₋, Val(3), Val(4)) + X₃₄ = spintrace_σμν(Minv, Val(3), Val(4)) Xμν[3i32, 4i32, site] = X₃₄ Xμν[4i32, 3i32, site] = -X₃₄ - else - X = zero3(T) - Xμν[1i32, 2i32, site] = X - Xμν[2i32, 1i32, site] = X - Xμν[1i32, 3i32, site] = X - Xμν[3i32, 1i32, site] = X - Xμν[1i32, 4i32, site] = X - Xμν[4i32, 1i32, site] = X - Xμν[2i32, 3i32, site] = X - Xμν[3i32, 2i32, site] = X - Xμν[2i32, 4i32, site] = X - Xμν[4i32, 2i32, site] = X - Xμν[3i32, 4i32, site] = X - Xμν[4i32, 3i32, site] = X + # else + # X = zero3(T) + # Xμν[1i32, 2i32, site] = X + # Xμν[2i32, 1i32, site] = X + # Xμν[1i32, 3i32, site] = X + # Xμν[3i32, 1i32, site] = X + # Xμν[1i32, 4i32, site] = X + # Xμν[4i32, 1i32, site] = X + # Xμν[2i32, 3i32, site] = X + # Xμν[3i32, 2i32, site] = X + # Xμν[2i32, 4i32, site] = X + # Xμν[4i32, 2i32, site] = X + # Xμν[3i32, 4i32, site] = X + # Xμν[4i32, 3i32, site] = X + # end +end + +function add_small_det_derivative!( + dU::Colorfield{CPU,T}, U::Gaugefield{CPU,T}, Xμν::Tensorfield{CPU,T}, csw; coeff=1 +) where {T} + check_dims(dU, U, Xμν) + fac = T(csw * coeff / 2) + + @batch for site in eachindex(dU) + add_small_det_derivative_kernel!(dU, U, Xμν, site, fac, T) end + + update_halo!(dU) + return nothing +end + +function add_small_det_derivative_kernel!(dU, U, Xμν, site, fac, ::Type{T}) where {T} + tmp = + Xμν∇Fμν(Xμν, U, 1, 2, site, T) + + Xμν∇Fμν(Xμν, U, 1, 3, site, T) + + Xμν∇Fμν(Xμν, U, 1, 4, site, T) + dU[1i32, site] += fac * traceless_antihermitian(cmatmul_oo(U[1i32, site], tmp)) + + tmp = + Xμν∇Fμν(Xμν, U, 2, 1, site, T) + + Xμν∇Fμν(Xμν, U, 2, 3, site, T) + + Xμν∇Fμν(Xμν, U, 2, 4, site, T) + dU[2i32, site] += fac * traceless_antihermitian(cmatmul_oo(U[2i32, site], tmp)) + + tmp = + Xμν∇Fμν(Xμν, U, 3, 1, site, T) + + Xμν∇Fμν(Xμν, U, 3, 2, site, T) + + Xμν∇Fμν(Xμν, U, 3, 4, site, T) + dU[3i32, site] += fac * traceless_antihermitian(cmatmul_oo(U[3i32, site], tmp)) + + tmp = + Xμν∇Fμν(Xμν, U, 4, 1, site, T) + + Xμν∇Fμν(Xμν, U, 4, 2, site, T) + + Xμν∇Fμν(Xμν, U, 4, 3, site, T) + dU[4i32, site] += fac * traceless_antihermitian(cmatmul_oo(U[4i32, site], tmp)) + return nothing end diff --git a/src/io/bmw_format.jl b/src/io/bmw_format.jl index aa35bb86..b653a28e 100644 --- a/src/io/bmw_format.jl +++ b/src/io/bmw_format.jl @@ -29,7 +29,9 @@ function save_config( ) where {B,T} @assert U.U isa Array if override == false - @assert !isfile(filename) "File $filename to store config in already exists and override is \"false\"" + @assert !isfile(filename) """ + File $filename to store config in already exists and override is \"false\" + """ end fp = open(filename, "w") N = U.NC @@ -132,7 +134,9 @@ function load_config!(::BMWFormat, U::Gaugefield{B,T,false}, filename) where {B, close(fp) adler64_finalize!(checksum_calc) - @assert checksum_read == checksum_calc.final "Checksums do not match (read: $checksum_read, calculated: $(checksum_calc.final))" + @assert checksum_read == checksum_calc.final """ + Checksums do not match (read: $checksum_read, calculated: $(checksum_calc.final)) + """ return nothing end diff --git a/src/io/ildg_format.jl b/src/io/ildg_format.jl index ae22ff85..0506d0f4 100644 --- a/src/io/ildg_format.jl +++ b/src/io/ildg_format.jl @@ -7,7 +7,9 @@ function save_config(::ILDGFormat, U, filename; parameters=nothing, override=fal @assert U.U isa Array if override == false - @assert !isfile(filename) "File $filename to store config in already exists and override is \"false\"" + @assert !isfile(filename) """ + File $filename to store config in already exists and override is \"false\" + """ end fp = open(filename, "w") @@ -94,7 +96,9 @@ function load_config!(::BMWFormat, U, filename) close(fp) adler64_finalize!(checksum_calc) - @assert checksum_read == checksum_calc.final "Checksums do not match (read: $checksum_read, calculated: $(checksum_calc.final))" + @assert checksum_read == checksum_calc.final """ + Checksums do not match (read: $checksum_read, calculated: $(checksum_calc.final)) + """ return nothing end diff --git a/src/main/Universe.jl b/src/main/Universe.jl index 81f7c957..c9062ae4 100644 --- a/src/main/Universe.jl +++ b/src/main/Universe.jl @@ -19,8 +19,8 @@ import ..Parameters: ParameterSet Create a Universe, containing the gauge configurations that are updated throughout a simulation and the bias potentials, if any. \\ -`mpi_multi_sim` is soley an indicator that sets `numinstances` to 1 when using multiple walkers -even if stated otherwise in the parameters. +`mpi_multi_sim` is soley an indicator that sets `numinstances` to 1 when using +multiple walkers even if stated otherwise in the parameters. """ struct Univ{TG,TF,TB} U::TG @@ -97,7 +97,9 @@ function Univ(parameters::ParameterSet; mpi_multi_sim=false) bias = Bias(parameters, U; mpi_multi_sim=mpi_multi_sim) end else - @assert parameters.tempering_enabled == false "tempering can only be enabled with bias" + @assert parameters.tempering_enabled == false """ + tempering can only be enabled with bias + """ numinstances = 1 U = Gaugefield(parameters) fermion_action = init_fermion_actions(parameters, U) @@ -111,7 +113,7 @@ function init_fermion_actions(parameters::ParameterSet, U) fermion_action = parameters.fermion_action Nf = parameters.Nf mass = parameters.mass - @assert length(Nf) == length(mass) "Need same amount of masses as non-degenerate flavours" + @assert length(Nf) == length(mass) "Need same amount of masses as unique flavors" if fermion_action ∈ ("none", "quenched") fermion_actions = QuenchedFermionAction() diff --git a/src/main/runbuild.jl b/src/main/runbuild.jl index aa7df7bf..0d4847fe 100644 --- a/src/main/runbuild.jl +++ b/src/main/runbuild.jl @@ -14,9 +14,16 @@ function build_bias(filenamein::String; backend="cpu") if mpi_amroot() oneinst = parameters.numinstances == 1 - @assert multi_sim ⊻ oneinst "MPI must be enabled if and only if numinstances > 1, numinstances was $(parameters.numinstances)" - @assert mpi_size() == parameters.numinstances "numinstances has to be equal to the number of MPI ranks" - @assert parameters.kind_of_bias ∉ ("none", "parametric") "bias has to be \"metad\" or \"opes\" in build, was $(parameters.kind_of_bias)" + @assert multi_sim ⊻ oneinst """ + MPI must be enabled if and only if numinstances > 1, numinstances was \ + $(parameters.numinstances) + """ + @assert mpi_size() == parameters.numinstances """ + numinstances has to be equal to the number of MPI ranks + """ + @assert parameters.kind_of_bias ∉ ("none", "parametric") """ + bias has to be \"metad\" or \"opes\" in build, was $(parameters.kind_of_bias) + """ @assert parameters.is_static == false "Bias cannot be static in build" end diff --git a/src/main/runsim.jl b/src/main/runsim.jl index 18014d1d..f95738c6 100644 --- a/src/main/runsim.jl +++ b/src/main/runsim.jl @@ -13,7 +13,9 @@ function run_sim(filenamein::String; backend="cpu") if mpi_amroot() if parameters.tempering_enabled - @assert parameters.kind_of_bias != "none" "bias cannot be \"none\" in parallel tempering" + @assert parameters.kind_of_bias != "none" """ + bias cannot be \"none\" in parallel tempering + """ end end diff --git a/src/measurements/measure_topological_charge.jl b/src/measurements/measure_topological_charge.jl index b4c9d88f..ce160d7f 100644 --- a/src/measurements/measure_topological_charge.jl +++ b/src/measurements/measure_topological_charge.jl @@ -156,38 +156,38 @@ end function top_charge_density_plaq(U, site) C₁₂ = plaquette(U, 1i32, 2i32, site) - F₁₂ = im * traceless_antihermitian(C₁₂) + F₁₂ = C₁₂ - C₁₂' C₁₃ = plaquette(U, 1i32, 3i32, site) - F₁₃ = im * traceless_antihermitian(C₁₃) + F₁₃ = C₁₃ - C₁₃' C₂₃ = plaquette(U, 2i32, 3i32, site) - F₂₃ = im * traceless_antihermitian(C₂₃) + F₂₃ = C₂₃ - C₂₃' C₁₄ = plaquette(U, 1i32, 4i32, site) - F₁₄ = im * traceless_antihermitian(C₁₄) + F₁₄ = C₁₄ - C₁₄' C₂₄ = plaquette(U, 2i32, 4i32, site) - F₂₄ = im * traceless_antihermitian(C₂₄) + F₂₄ = C₂₄ - C₂₄' C₃₄ = plaquette(U, 3i32, 4i32, site) - F₃₄ = im * traceless_antihermitian(C₃₄) + F₃₄ = C₃₄ - C₃₄' qₙ = real(multr(F₁₂, F₃₄)) - real(multr(F₁₃, F₂₄)) + real(multr(F₁₄, F₂₃)) - return qₙ + return -qₙ end function top_charge_density_clover(U, site, ::Type{T}) where {T} C₁₂ = clover_square(U, 1i32, 2i32, site, 1i32) - F₁₂ = im * T(1/4) * traceless_antihermitian(C₁₂) + F₁₂ = C₁₂ - C₁₂' C₁₃ = clover_square(U, 1i32, 3i32, site, 1i32) - F₁₃ = im * T(1/4) * traceless_antihermitian(C₁₃) + F₁₃ = C₁₃ - C₁₃' C₂₃ = clover_square(U, 2i32, 3i32, site, 1i32) - F₂₃ = im * T(1/4) * traceless_antihermitian(C₂₃) + F₂₃ = C₂₃ - C₂₃' C₁₄ = clover_square(U, 1i32, 4i32, site, 1i32) - F₁₄ = im * T(1/4) * traceless_antihermitian(C₁₄) + F₁₄ = C₁₄ - C₁₄' C₂₄ = clover_square(U, 2i32, 4i32, site, 1i32) - F₂₄ = im * T(1/4) * traceless_antihermitian(C₂₄) + F₂₄ = C₂₄ - C₂₄' C₃₄ = clover_square(U, 3i32, 4i32, site, 1i32) - F₃₄ = im * T(1/4) * traceless_antihermitian(C₃₄) + F₃₄ = C₃₄ - C₃₄' - out = real(multr(F₁₂, F₃₄)) - real(multr(F₁₃, F₂₄)) + real(multr(F₁₄, F₂₃)) - return out + qₙ = real(multr(F₁₂, F₃₄)) - real(multr(F₁₃, F₂₄)) + real(multr(F₁₄, F₂₃)) + return -T(1/64) * qₙ end function top_charge_density_imp(U, site, c₀, c₁, ::Type{T}) where {T} @@ -199,18 +199,18 @@ end function top_charge_density_rect(U, site, ::Type{T}) where {T} C₁₂ = clover_rect(U, 1i32, 2i32, site, 1i32, 2i32) - F₁₂ = im * T(1/8) * traceless_antihermitian(C₁₂) + F₁₂ = C₁₂ - C₁₂' C₁₃ = clover_rect(U, 1i32, 3i32, site, 1i32, 2i32) - F₁₃ = im * T(1/8) * traceless_antihermitian(C₁₃) + F₁₃ = C₁₃ - C₁₃' C₂₃ = clover_rect(U, 2i32, 3i32, site, 1i32, 2i32) - F₂₃ = im * T(1/8) * traceless_antihermitian(C₂₃) + F₂₃ = C₂₃ - C₂₃' C₁₄ = clover_rect(U, 1i32, 4i32, site, 1i32, 2i32) - F₁₄ = im * T(1/8) * traceless_antihermitian(C₁₄) + F₁₄ = C₁₄ - C₁₄' C₂₄ = clover_rect(U, 2i32, 4i32, site, 1i32, 2i32) - F₂₄ = im * T(1/8) * traceless_antihermitian(C₂₄) + F₂₄ = C₂₄ - C₂₄' C₃₄ = clover_rect(U, 3i32, 4i32, site, 1i32, 2i32) - F₃₄ = im * T(1/8) * traceless_antihermitian(C₃₄) + F₃₄ = C₃₄ - C₃₄' - out = real(multr(F₁₂, F₃₄)) - real(multr(F₁₃, F₂₄)) + real(multr(F₁₄, F₂₃)) - return out + qₙ = real(multr(F₁₂, F₃₄)) - real(multr(F₁₃, F₂₄)) + real(multr(F₁₄, F₂₃)) + return -T(1/256) * qₙ end diff --git a/src/parameters/Parameters.jl b/src/parameters/Parameters.jl index 53ff88b9..6ce7c5c9 100644 --- a/src/parameters/Parameters.jl +++ b/src/parameters/Parameters.jl @@ -215,14 +215,21 @@ end function parameter_check(p::ParameterSet) mpi_amroot() || return nothing - @assert prod(p.numprocs_cart) == mpi_size() "Size of comm must equal number of process used in field decomposition" + @assert prod(p.numprocs_cart) == mpi_size() """ + Size of comm must equal number of process used in field decomposition + """ if prod(p.numprocs_cart) > 1 @assert p.halo_width >= 1 "Halo width must be >= 1, when using field decomposition" - @assert lower_case(p.update_method) == "hmc" "Field decomposition not supported for local update algorithms" + @assert lower_case(p.update_method) == "hmc" """ + Field decomposition not supported for local update algorithms + """ if p.gauge_action != "wilson" - @assert p.halo_width >= 2 "Halo width must be >= 2, when using field decomposition with improved gauge action" + @assert p.halo_width >= 2 """ + Halo width must be >= 2, when using field decomposition with improved \ + gauge action + """ end end @@ -271,7 +278,9 @@ function parameter_check(p::ParameterSet) """)) end - if lower_case(p.hmc_integrator) ∉ ["leapfrog", "omf2slow", "omf2", "omf4slow", "omf4", "leapfrogra"] + if lower_case(p.hmc_integrator) ∉ [ + "leapfrog", "omf2slow", "omf2", "omf4slow", "omf4", "leapfrogra" + ] throw(AssertionError(""" hmc_integrator in [\"HMC Settings\"] = $(p.hmc_integrator) is not supported. Supported methods are: @@ -329,7 +338,8 @@ function parameter_check(p::ParameterSet) if lower_case(p.save_config_format) ∉ ["", "bmw", "bridge", "jld", "jld2"] throw(AssertionError(""" - save_config_format in [\"System Settings\"] = $(p.save_config_format) is not supported. + save_config_format in [\"System Settings\"] = $(p.save_config_format) \ + is not supported. Supported methods are: Bridge JLD or JLD2 (both use JLD2) @@ -341,7 +351,8 @@ function parameter_check(p::ParameterSet) @assert isfile(p.load_config_path) "Your load_config_path doesn't exist" if lower_case(p.load_config_format) ∉ ["bmw", "bridge", "jld", "jld2"] throw(AssertionError(""" - loadU_format in [\"System Settings\"] = $(p.load_config_format) is not supported. + loadU_format in [\"System Settings\"] = $(p.load_config_format) \ + is not supported. Supported methods are: Bridge JLD or JLD2 (both use JLD2) @@ -349,6 +360,7 @@ function parameter_check(p::ParameterSet) """)) end end + return nothing end diff --git a/src/solvers/cg.jl b/src/solvers/cg.jl index 12c3a416..e3d2d732 100644 --- a/src/solvers/cg.jl +++ b/src/solvers/cg.jl @@ -154,28 +154,37 @@ function bicg_stab!(x, A, b, v, r, p, r₀, t; tol=1e-14, maxiters=1000) axpy!(-α, v, r) res = abs(dot(r, r)) @level3 "| BiCGStab: residual $(iter).5 = $res" + if res < tol @level2 "| BiCGStab: converged at iter $(iter).5 with res = $res" return nothing end - @assert isfinite(res) && isfinite(α) "BiCG: NaN or Inf encountered, res = $res, α = $α" + + @assert isfinite(res) && isfinite(α) """ + BiCG: NaN or Inf encountered, res = $res, α = $α + """ mul!(t, A, r) ω = dot(t, r) / dot(t, t) axpy!(ω, r, x) axpy!(-ω, t, r) res = abs(dot(r, r)) @level3 "| BiCGStab: residual $(iter) = $res" + if res < tol @level2 "| BiCGStab: converged at iter $(iter) with res = $res" return nothing end - @assert isfinite(res) && isfinite(ω) "BiCG: NaN or Inf encountered, res = $res, ω = $ω" + + @assert isfinite(res) && isfinite(ω) """ + BiCG: NaN or Inf encountered, res = $res, ω = $ω + """ ρ_new = dot(r₀, r) β = (ρ_new / ρ) * (α / ω) axpy!(-ω, v, p) axpby!(1, r, β, p) ρ = ρ_new end + # @level1 "| BiCGStab: did not converge in $maxiters iterations" throw(AssertionError("BiCGStab did not converge in $maxiters iterations")) return nothing diff --git a/src/updates/heatbath.jl b/src/updates/heatbath.jl index 6f77ce78..ff00232e 100644 --- a/src/updates/heatbath.jl +++ b/src/updates/heatbath.jl @@ -24,7 +24,9 @@ struct Heatbath{MAXIT,ITR,TOR,NHB,NOR} <: AbstractUpdate end # @inline NHB(::Heatbath{<:Any,<:Any,<:Any,NHB}) where {NHB} = _unwrap_val(NHB) # @inline NOR(::Heatbath{<:Any,<:Any,<:Any,<:Any,NOR}) where {NOR} = _unwrap_val(NOR) -function Heatbath(::Gaugefield{B,T,A,GA}, MAXIT, numheatbath, or_alg, numorelax) where {B,T,A,GA} +function Heatbath( + ::Gaugefield{B,T,A,GA}, MAXIT, numheatbath, or_alg, numorelax +) where {B,T,A,GA} @level1("┌ Setting Heatbath...") ITR = GA == WilsonGaugeAction ? Checkerboard2 : Checkerboard4 @level1("| ITERATOR: $(string(ITR))") diff --git a/src/updates/hmc_integrators.jl b/src/updates/hmc_integrators.jl index 816f438b..b461b6fa 100644 --- a/src/updates/hmc_integrators.jl +++ b/src/updates/hmc_integrators.jl @@ -24,6 +24,7 @@ end function Base.show(io::IO, ::MIME"text/plain", int::LeapfrogRA) print(io, "$(typeof(int))(friction=$(int.friction))") end + Base.show(io::IO, int::LeapfrogRA) = print(io, "$(typeof(int))(friction=$(int.friction))") function evolve!(L::LeapfrogRA, U, hmc::HMC, fermion_action, bias) diff --git a/src/updates/tempering.jl b/src/updates/tempering.jl index 487e6815..e418381a 100644 --- a/src/updates/tempering.jl +++ b/src/updates/tempering.jl @@ -64,9 +64,9 @@ function temper!( mpi_bcast!(instance_state, comm; root=rank_i) mpi_bcast!(numaccepts_temper, comm; root=rank_i) mpi_barrier() - @level1( - "| Acceptance [$i ⇔ $(i-1)]:\t$(100numaccepts_temper[i-1] / (itrj/swap_every)) %" - ) + @level1 """ + | Acceptance [$i ⇔ $(i-1)]:\t$(100numaccepts_temper[i-1] / (itrj/swap_every)) % + """ end return nothing end @@ -100,9 +100,9 @@ function temper!( println("# swap rejected") end - @level1( - "| Acceptance [$i ⇔ $(i-1)]:\t$(100numaccepts_temper[i-1] / (itrj/swap_every)) %" - ) + @level1 """ + "| Acceptance [$i ⇔ $(i-1)]:\t$(100numaccepts_temper[i-1] / (itrj/swap_every)) % + """ end return nothing diff --git a/src/utils/Utils.jl b/src/utils/Utils.jl index ba4e639c..7463af2d 100644 --- a/src/utils/Utils.jl +++ b/src/utils/Utils.jl @@ -14,7 +14,7 @@ export mpi_init, mpi_comm, mpi_size, mpi_parallel, mpi_myrank, mpi_amroot, mpi_b export mpi_cart_create, mpi_cart_coords, mpi_cart_shift, mpi_multirequest, mpi_send export mpi_isend, mpi_recv, mpi_irecv!, mpi_waitall, mpi_allreduce, mpi_allgather export mpi_bcast_isbits, mpi_write_at, update_halo! -export exp_iQ, exp_iQ_coeffs, exp_iQ_su3, get_B₁, get_B₂, get_Q, get_Q² +export PauliMatrix, exp_iQ, exp_iQ_coeffs, exp_iQ_su3, get_B₁, get_B₂, get_Q, get_Q² export gen_SU3_matrix, is_special_unitary, is_traceless_antihermitian export kenney_laub, proj_onto_SU3, multr, cnorm2 export make_submatrix_12, make_submatrix_13, make_submatrix_23 @@ -24,7 +24,7 @@ export zero2, zero3, zerov3, eye2, eye3, onev3, gaussian_TA_mat, rand_SU3 export SiteCoords, eo_site, eo_site_switch, move, switch_sides export cartesian_to_linear export Sequential, Checkerboard2, Checkerboard4 -export λ, expλ, γ₁, γ₂, γ₃, γ₄, γ₅, σ₁₂, σ₁₃, σ₁₄, σ₂₃, σ₂₄, σ₃₄ +export λ, expλ, γ1, γ2, γ3, γ4, γ5, σ12, σ13, σ14, σ23, σ24, σ34 export cmatmul_oo, cmatmul_dd, cmatmul_do, cmatmul_od export cmatmul_ooo, cmatmul_ood, @@ -119,6 +119,28 @@ const i32 = Literal{Int32} const SU{N,N²,T} = SMatrix{N,N,Complex{T},N²} +struct PauliMatrix{T<:AbstractFloat} + upper::SU{6,36,T} + lower::SU{6,36,T} + function PauliMatrix(λ::UniformScaling{T}) where {T<:AbstractFloat} + upper = lower = @SMatrix(zeros(Complex{T}, 6, 6)) + λ + return new{T}(upper, lower) + end + + function PauliMatrix(upper::SU{6,36,T}, lower::SU{6,36,T}) where {T<:AbstractFloat} + return new{T}(upper, lower) + end +end + +Base.zero(::Type{PauliMatrix{T}}) where {T} = PauliMatrix(UniformScaling(zero(T))) +Base.one(::Type{PauliMatrix{T}}) where {T} = PauliMatrix(UniformScaling(one(T))) + +function Base.rand(::Type{PauliMatrix{T}}) where {T} + upper = hermitian(@SMatrix(rand(Complex{T}, 6, 6))) + lower = hermitian(@SMatrix(rand(Complex{T}, 6, 6))) + return PauliMatrix(upper, lower) +end + """ multr(A::SMatrix{N,N,Complex{T},N²}, B::SMatrix{N,N,Complex{T},N²}) where {N,N²,T} diff --git a/src/utils/clover_obj.jl b/src/utils/clover_obj.jl deleted file mode 100644 index 7f1ab506..00000000 --- a/src/utils/clover_obj.jl +++ /dev/null @@ -1,19 +0,0 @@ -struct Pauli{T<:AbstractFloat} - data::NTuple{36,T} - Pauli(x::NTuple{36,T}) where {T<:AbstractFloat} = new{T}(x) - - function Pauli(λ::UniformScaling{T}) where {T<:AbstractFloat} - c = λ.λ - rest = ntuple(T(0.0), Val(30)) - x = (c, c, c, c, c, rest...) - return new{T}(x) - end -end - -Base.zeros(::Type{Pauli{T}}) where {T} = Pauli(ntuple(_ -> T(0.0), Val(36))) - -function pauli_mat(p::Pauli{T}) where {T} - out = MMatrix{6,6,Complex{T},36}(undef) - - -end diff --git a/src/utils/generators.jl b/src/utils/generators.jl index 4a9428e8..f8eb9329 100644 --- a/src/utils/generators.jl +++ b/src/utils/generators.jl @@ -1,5 +1,5 @@ # Chiral(Weyl)-Basis γ matrices -@inline function γ₁(::Type{T}) where {T} +@inline function γ1(::Type{T}) where {T} return @SArray [ Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, -1) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, -1) Complex{T}(0, 0) @@ -8,7 +8,7 @@ ] end -@inline function γ₂(::Type{T}) where {T} +@inline function γ2(::Type{T}) where {T} return @SArray [ Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(-1, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(1, 0) Complex{T}(0, 0) @@ -17,7 +17,7 @@ end ] end -@inline function γ₃(::Type{T}) where {T} +@inline function γ3(::Type{T}) where {T} return @SArray [ Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, -1) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 1) @@ -26,7 +26,7 @@ end ] end -@inline function γ₄(::Type{T}) where {T} +@inline function γ4(::Type{T}) where {T} return @SArray [ Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(1, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(1, 0) @@ -35,7 +35,7 @@ end ] end -@inline function γ₅(::Type{T}) where {T} +@inline function γ5(::Type{T}) where {T} return @SArray [ Complex{T}(1, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(1, 0) Complex{T}(0, 0) Complex{T}(0, 0) @@ -44,23 +44,23 @@ end ] end -const σ₁ = @SArray [ +const σ1 = @SArray [ 0 1 1 0 ] -const σ₂ = @SArray [ +const σ2 = @SArray [ 0 -im im 0 ] -const σ₃ = @SArray [ +const σ3 = @SArray [ 1 0 0 -1 ] # σ_μν = i/2 * [γ_μ, γ_ν] -@inline function σ₁₂(::Type{T}) where {T} +@inline function σ12(::Type{T}) where {T} return @SArray [ Complex{T}(-1, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(1, 0) Complex{T}(0, 0) Complex{T}(0, 0) @@ -69,7 +69,7 @@ const σ₃ = @SArray [ ] end -@inline function σ₁₃(::Type{T}) where {T} +@inline function σ13(::Type{T}) where {T} return @SArray [ Complex{T}(0, 0) Complex{T}(0, -1) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 1) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 0) @@ -78,7 +78,7 @@ end ] end -@inline function σ₁₄(::Type{T}) where {T} +@inline function σ14(::Type{T}) where {T} return @SArray [ Complex{T}(0, 0) Complex{T}(1, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(1, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 0) @@ -87,7 +87,7 @@ end ] end -@inline function σ₂₃(::Type{T}) where {T} +@inline function σ23(::Type{T}) where {T} return @SArray [ Complex{T}(0, 0) Complex{T}(-1, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(-1, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 0) @@ -96,7 +96,7 @@ end ] end -@inline function σ₂₄(::Type{T}) where {T} +@inline function σ24(::Type{T}) where {T} return @SArray [ Complex{T}(0, 0) Complex{T}(0, -1) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 1) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 0) @@ -105,7 +105,7 @@ end ] end -@inline function σ₃₄(::Type{T}) where {T} +@inline function σ34(::Type{T}) where {T} return @SArray [ Complex{T}(1, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(0, 0) Complex{T}(-1, 0) Complex{T}(0, 0) Complex{T}(0, 0) @@ -114,57 +114,61 @@ end ] end -const λ₁ = @SArray [ +@generated function σ(::Val{μ}, ::Val{ν}) where {μ,ν} + return Symbol(:σ, "$μ$ν") +end + +const λ1 = @SArray [ 0.0+0.0im 1.0+0.0im 0.0+0.0im 1.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im ] -const λ₂ = @SArray [ +const λ2 = @SArray [ 0.0+0.0im 0.0-1.0im 0.0+0.0im 0.0+1.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im ] -const λ₃ = @SArray [ +const λ3 = @SArray [ 1.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im -1.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im ] -const λ₄ = @SArray [ +const λ4 = @SArray [ 0.0+0.0im 0.0+0.0im 1.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 1.0+0.0im 0.0+0.0im 0.0+0.0im ] -const λ₅ = @SArray [ +const λ5 = @SArray [ 0.0+0.0im 0.0+0.0im 0.0-1.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+1.0im 0.0+0.0im 0.0+0.0im ] -const λ₆ = @SArray [ +const λ6 = @SArray [ 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 1.0+0.0im 0.0+0.0im 1.0+0.0im 0.0+0.0im ] -const λ₇ = @SArray [ +const λ7 = @SArray [ 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-1.0im 0.0+0.0im 0.0+1.0im 0.0+0.0im ] -const λ₈ = @SArray [ +const λ8 = @SArray [ 1/sqrt(3) 0.0+0.0im 0.0+0.0im 0.0+0.0im 1/sqrt(3) 0.0+0.0im 0.0+0.0im 0.0+0.0im -2/sqrt(3) ] -const λ = (λ₁, λ₂, λ₃, λ₄, λ₅, λ₆, λ₇, λ₈) +const λ = (λ1, λ2, λ3, λ4, λ5, λ6, λ7, λ8) -function expλ₁(α) +function expλ1(α) φ = 0.5 * α sinφ = sin(φ) cosφ = cos(φ) @@ -177,7 +181,7 @@ function expλ₁(α) return out end -function expλ₂(α) +function expλ2(α) φ = 0.5 * α sinφ = sin(φ) cosφ = cos(φ) @@ -190,7 +194,7 @@ function expλ₂(α) return out end -function expλ₃(α) +function expλ3(α) φ = 0.5 * α expiφ = cis(φ) @@ -202,7 +206,7 @@ function expλ₃(α) return out end -function expλ₄(α) +function expλ4(α) φ = 0.5 * α sinφ = sin(φ) cosφ = cos(φ) @@ -215,7 +219,7 @@ function expλ₄(α) return out end -function expλ₅(α) +function expλ5(α) φ = 0.5 * α sinφ = sin(φ) cosφ = cos(φ) @@ -228,7 +232,7 @@ function expλ₅(α) return out end -function expλ₆(α) +function expλ6(α) φ = 0.5 * α sinφ = sin(φ) cosφ = cos(φ) @@ -241,7 +245,7 @@ function expλ₆(α) return out end -function expλ₇(α) +function expλ7(α) φ = 0.5 * α sinφ = sin(φ) cosφ = cos(φ) @@ -254,7 +258,7 @@ function expλ₇(α) return out end -function expλ₈(α) +function expλ8(α) φ = 0.5 / sqrt(3) * α expiφ = cis(φ) @@ -268,21 +272,21 @@ end function expλ(i, α) if i == 1 - return expλ₁(α) + return expλ1(α) elseif i == 2 - return expλ₂(α) + return expλ2(α) elseif i == 3 - return expλ₃(α) + return expλ3(α) elseif i == 4 - return expλ₄(α) + return expλ4(α) elseif i == 5 - return expλ₅(α) + return expλ5(α) elseif i == 6 - return expλ₆(α) + return expλ6(α) elseif i == 7 - return expλ₇(α) + return expλ7(α) elseif i == 8 - return expλ₈(α) + return expλ8(α) else return eye3(Float64) end diff --git a/src/utils/pauli.jl b/src/utils/pauli.jl new file mode 100644 index 00000000..da2a4628 --- /dev/null +++ b/src/utils/pauli.jl @@ -0,0 +1,23 @@ +struct PauliMatrix{T<:AbstractFloat} + upper::SMatrix{6,6,Complex{T},36} + lower::SMatrix{6,6,Complex{T},36} + function PauliMatrix(λ::UniformScaling{T}) where {T<:AbstractFloat} + upper = lower = @SMatrix(zeros(Complex{T}, 6, 6)) + λ + return new{T}(upper, lower) + end + + function PauliMatrix( + upper::SMatrix{6,6,Complex{T},36}, lower::SMatrix{6,6,Complex{T},36} + ) where {T<:AbstractFloat} + return new{T}(upper, lower) + end +end + +Base.zero(::Type{PauliMatrix{T}}) where {T} = PauliMatrix(UniformScaling(zero(T))) +Base.one(::Type{PauliMatrix{T}}) where {T} = PauliMatrix(UniformScaling(one(T))) + +function rand_pauli(::Type{PauliMatrix{T}}) where {T} + upper = hermitian(@SMatrix(rand(Complex{T}, 6, 6))) + lower = hermitian(@SMatrix(rand(Complex{T}, 6, 6))) + return PauliMatrix(upper, lower) +end diff --git a/src/utils/simd_vecops.jl b/src/utils/simd_vecops.jl index db6196db..e5d13be9 100644 --- a/src/utils/simd_vecops.jl +++ b/src/utils/simd_vecops.jl @@ -932,8 +932,8 @@ end x = reinterpret(reshape, $T, xc) end - if μ === 1 && ν === 2 - inner_q = quote + inner_q = if μ === 1 && ν === 2 + quote y[1, m] = -x[1, m] y[2, m] = -x[2, m] y[1, $N+m] = x[1, $N+m] @@ -944,7 +944,7 @@ end y[2, $(3N)+m] = x[2, $(3N)+m] end elseif μ === 1 && ν === 3 - inner_q = quote + quote y[1, m] = x[2, $N+m] y[2, m] = -x[1, $N+m] y[1, $N+m] = -x[2, m] @@ -955,7 +955,7 @@ end y[2, $(3N)+m] = x[1, $(2N)+m] end elseif μ === 1 && ν === 4 - inner_q = quote + quote y[1, m] = x[1, $N+m] y[2, m] = x[2, $N+m] y[1, $N+m] = x[1, m] @@ -966,7 +966,7 @@ end y[2, $(3N)+m] = -x[2, $(2N)+m] end elseif μ === 2 && ν === 3 - inner_q = quote + quote y[1, m] = -x[1, $N+m] y[2, m] = -x[2, $N+m] y[1, $N+m] = -x[1, m] @@ -977,7 +977,7 @@ end y[2, $(3N)+m] = -x[2, $(2N)+m] end elseif μ === 2 && ν === 4 - inner_q = quote + quote y[1, m] = x[2, $N+m] y[2, m] = -x[1, $N+m] y[1, $N+m] = -x[2, m] @@ -988,7 +988,7 @@ end y[2, $(3N)+m] = -x[1, $(2N)+m] end elseif μ === 3 && ν === 4 - inner_q = quote + quote y[1, m] = x[1, m] y[2, m] = x[2, m] y[1, $N+m] = -x[1, $N+m] @@ -1006,6 +1006,7 @@ end @turbo for m in Base.Slice(static(1):static($N)) $inner_q end + return yc end @@ -1014,50 +1015,34 @@ end end """ - spintrace_σμν(A, B, ::Val{μ}, ::Val{ν}) + spintrace_σμν(P::PauliMatrix, ::Val{μ}, ::Val{ν}) """ @generated function spintrace_σμν( - A::SMatrix{M,M,Complex{T},MM}, - B::SMatrix{M,M,Complex{T},MM}, + P::PauliMatrix{T}, ::Val{μ}, ::Val{ν} -) where {T,M,MM,μ,ν} - if M % 2 != 0 - return :(throw(DimensionMismatch("M must be a multiple of 2"))) - end - - N = M ÷ 2 - i = SVector(1:N...) - j = SVector(N+1:M...) - +) where {T,μ,ν} q = quote $(Expr(:meta, :inline)) + A = P.upper + B = P.lower end - if μ === 1 && ν === 2 - inner_q = quote - Cc = -A[$i, $i] + A[$j, $j] - B[$i, $i] + B[$j, $j] - end + i = SVector(1:3...) + j = SVector(4:6...) + + inner_q = if μ === 1 && ν === 2 + :(return -A[$i, $i] + A[$j, $j] - B[$i, $i] + B[$j, $j]) elseif μ === 1 && ν === 3 - inner_q = quote - Cc = im * (A[$j, $i] - A[$i, $j] + B[$j, $i] - B[$i, $j]) - end + :(return im * (A[$j, $i] - A[$i, $j] + B[$j, $i] - B[$i, $j])) elseif μ === 1 && ν === 4 - inner_q = quote - Cc = A[$j, $i] + A[$i, $j] - B[$j, $i] - B[$i, $j] - end + :(return A[$j, $i] + A[$i, $j] - B[$j, $i] - B[$i, $j]) elseif μ === 2 && ν === 3 - inner_q = quote - Cc = -A[$j, $i] - A[$i, $j] - B[$j, $i] - B[$i, $j] - end + :(return -A[$j, $i] - A[$i, $j] - B[$j, $i] - B[$i, $j]) elseif μ === 2 && ν === 4 - inner_q = quote - Cc = im * (A[$j, $i] - A[$i, $j] - B[$j, $i] + B[$i, $j]) - end + :(return im * (A[$j, $i] - A[$i, $j] - B[$j, $i] + B[$i, $j])) elseif μ === 3 && ν === 4 - inner_q = quote - Cc = A[$i, $i] - A[$j, $j] - B[$i, $i] + B[$j, $j] - end + :(return A[$i, $i] - A[$j, $j] - B[$i, $i] + B[$j, $j]) else return :(throw(AssertionError("Invalid combination of μ and ν"))) end @@ -1067,20 +1052,21 @@ end end """ - cmvmul_block(A₊, A₋, x) + cmvmul_block(P::PauliMatrix, x::SVector) Return the matrix-vector product of the block diagonal matrix containing `A₊` and `A₋` and the vector `x` """ @inline function cmvmul_block( - A₊::SMatrix{N,N,Complex{T},NN}, - A₋::SMatrix{N,N,Complex{T},NN}, + P::PauliMatrix{T}, x::SVector{N2,Complex{T}} -) where {T,N,NN,N2} - idx1 = SVector(1:N...) - idx2 = SVector((N+1):2N...) +) where {T,N2} + idx1 = SVector(1:6...) + idx2 = SVector(7:12...) + A₊ = P.upper + A₋ = P.lower x₊ = cmvmul(A₊, x[idx1]) - x₋ = cmvmul(A₊, x[idx2]) + x₋ = cmvmul(A₋, x[idx2]) return vcat(x₊, x₋) end diff --git a/test/runtests.jl b/test/runtests.jl index 33d7155c..5238ffb7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -72,11 +72,11 @@ mpi_barrier() backend; nprocs_cart=nprocs_cart, halo_width=halo_width, dirac="wilson", mass=0.01, single_flavor=true, csw=1.78 ) - # wilson eo-pre derivative - test_fderivative( - backend; nprocs_cart=nprocs_cart, halo_width=halo_width, - dirac="wilson", mass=0.01, single_flavor=false, eoprec=true, csw=0 # TODO:Nf=1 test with eo wils - ) + # TODO: wilson eo-pre derivative + # test_fderivative( + # backend; nprocs_cart=nprocs_cart, halo_width=halo_width, + # dirac="wilson", mass=0.01, single_flavor=false, eoprec=true, csw=0 # TODO:Nf=1 test with eo wils + # ) # TODO: wilson-clover eo-pre derivative # test_fderivative( # backend; nprocs_cart=nprocs_cart, halo_width=halo_width, diff --git a/test/test_derivative.jl b/test/test_derivative.jl index 5f0123ad..e6d3b350 100644 --- a/test/test_derivative.jl +++ b/test/test_derivative.jl @@ -3,7 +3,7 @@ using MetaQCD using MetaQCD.Utils using MetaQCD.Updates: calc_dQdU_bare! -function test_derivative(backend=CPU; nprocs_cart=(1, 1, 1, 1), halo_width) +function test_derivative(backend=CPU; nprocs_cart=(1, 1, 1, 1), halo_width=0) Random.seed!(123) mpi_amroot() && println("Gauge and Clover derivative test") NX = 4 diff --git a/test/test_fderivative.jl b/test/test_fderivative.jl index 4d042d32..6c08a89c 100644 --- a/test/test_fderivative.jl +++ b/test/test_fderivative.jl @@ -6,7 +6,7 @@ using Random function test_fderivative( backend=CPU; nprocs_cart=(1, 1, 1, 1), - halo_width, + halo_width=0, dirac="staggered", mass=0.01, eoprec=false, @@ -23,16 +23,22 @@ function test_fderivative( NY = 4 NZ = 4 NT = 4 - U = Gaugefield{CPU,Float64,WilsonGaugeAction}(NX, NY, NZ, NT, 6.0, nprocs_cart, halo_width) + U = Gaugefield{CPU,Float64,WilsonGaugeAction}( + NX, NY, NZ, NT, 6.0, nprocs_cart, halo_width + ) random_gauges!(U) - # filename = pkgdir(MetaQCD, "test", "testconf.txt") - # load_config!(BridgeFormat(), U, filename) + filename = pkgdir(MetaQCD, "test", "testconf.txt") + load_config!(BridgeFormat(), U, filename) if backend !== CPU U = MetaQCD.to_backend(backend, U) end - ψ = eoprec ? even_odd(Spinorfield(U; staggered=dirac=="staggered")) : Spinorfield(U; staggered=dirac=="staggered") + ψ = if eoprec + even_odd(Spinorfield(U; staggered=dirac=="staggered")) + else + Spinorfield(U; staggered=dirac=="staggered") + end spectral_bound, Nf = if dirac=="staggered" (mass^2, 6.0), (single_flavor ? 1 : (eoprec ? 4 : 8)) @@ -58,7 +64,7 @@ function test_fderivative( ) action = MetaQCD.DiracOperators.init_fermion_action(params, mass, Nf, U) - mpi_amroot() && (@show action) + # mpi_amroot() && (@show action) sample_pseudofermions!(ψ, action, U) @@ -121,10 +127,10 @@ function test_fderivative( symm_diff = (action_new_fwd - action_new_bwd) / 2ΔH symm_diff_smeared = (action_new_fwd_smeared - action_new_bwd_smeared) / 2ΔH - # if group_direction == 1 - # @show daction_proj - # @show symm_diff - # end + if group_direction == 1 + @show daction_proj + @show symm_diff + end relerrors[group_direction, 1] = (symm_diff - daction_proj) / symm_diff relerrors[group_direction, 2] = (symm_diff_smeared - daction_proj_smeared) / symm_diff_smeared @@ -136,10 +142,10 @@ function test_fderivative( end end - if mpi_amroot() - println() - @test length(findall(x -> abs(x) > 1e-4, relerrors[:, 2])) == 0 - end + # if mpi_amroot() + # println() + # @test length(findall(x -> abs(x) > 1e-4, relerrors[:, 2])) == 0 + # end mpi_barrier() return relerrors diff --git a/test/test_mpi_derivative.jl b/test/test_mpi_derivative.jl deleted file mode 100644 index 48ee556d..00000000 --- a/test/test_mpi_derivative.jl +++ /dev/null @@ -1,127 +0,0 @@ -using MetaQCD -using MetaQCD.Updates: calc_dQdU_bare! -using MetaQCD.Utils -using MetaQCD.Fields: update_halo! -using LinearAlgebra -using Random -using Test - -function test_mpi_derivative() - if mpi_amroot() - println("test_mpi_derivative") - end - - Ns = 4 - Nt = 4 - U = Gaugefield{CPU,Float64,WilsonGaugeAction}(Ns, Ns, Ns, Nt, 5.7, (1, 1, 2, 2), 1) - random_gauges!(U) - - smearing = StoutSmearing(U, 5, 0.12) - staples = Colorfield(U) - fieldstrength = Tensorfield(U) - temp_force = Colorfield(U) - dSdU = Colorfield(U) - dSdU_smeared = Colorfield(U) - dQdU = Colorfield(U) - dQdU_smeared = Colorfield(U) - - site = SiteCoords(2, 3, 2, 3) - μ = 3 - ΔH = 0.00001 - - relerrors = Matrix{Float64}(undef, 8, 4) - - for group_direction in 1:8 - # Unsmeared - Ufwd = deepcopy(U) - if mpi_amroot() - Ufwd[μ, site] = expλ(group_direction, ΔH) * Ufwd[μ, site] - end - update_halo!(Ufwd) - gaction_new_fwd = calc_gauge_action(Ufwd) - topcharge_new_fwd = top_charge(Clover(), Ufwd) - - Ubwd = deepcopy(U) - if mpi_amroot() - Ubwd[μ, site] = expλ(group_direction, -ΔH) * Ubwd[μ, site] - end - update_halo!(Ubwd) - gaction_new_bwd = calc_gauge_action(Ubwd) - topcharge_new_bwd = top_charge(Clover(), Ubwd) - - # Smeared - Ufwd = deepcopy(U) - if mpi_amroot() - Ufwd[μ, site] = expλ(group_direction, ΔH) * Ufwd[μ, site] - end - update_halo!(Ufwd) - calc_smearedU!(smearing, Ufwd) - gaction_new_fwd_smeared = calc_gauge_action(smearing.Usmeared_multi[end]) - topcharge_new_fwd_smeared = top_charge(Clover(), smearing.Usmeared_multi[end]) - - Ubwd = deepcopy(U) - if mpi_amroot() - Ubwd[μ, site] = expλ(group_direction, -ΔH) * Ubwd[μ, site] - end - update_halo!(Ubwd) - calc_smearedU!(smearing, Ubwd) - gaction_new_bwd_smeared = calc_gauge_action(smearing.Usmeared_multi[end]) - topcharge_new_bwd_smeared = top_charge(Clover(), smearing.Usmeared_multi[end]) - - calc_dSdU_bare!(dSdU, staples, U, nothing, NoSmearing()) - calc_dSdU_bare!(dSdU_smeared, staples, U, temp_force, smearing) - calc_dQdU_bare!(Clover(), dQdU, fieldstrength, U, nothing, NoSmearing()) - calc_dQdU_bare!(Clover(), dQdU_smeared, fieldstrength, U, temp_force, smearing) - - dgaction_proj = real(multr(im * λ[group_direction], dSdU[μ, site])) - dtopcharge_proj = real(multr(im * λ[group_direction], dQdU[μ, site])) - dgaction_proj_smeared = real(multr(im * λ[group_direction], dSdU_smeared[μ, site])) - dtopcharge_proj_smeared = real( - multr(im * λ[group_direction], dQdU_smeared[μ, site]) - ) - - ga_symm_diff = (gaction_new_fwd - gaction_new_bwd) / 2ΔH - tc_symm_diff = (topcharge_new_fwd - topcharge_new_bwd) / 2ΔH - ga_symm_diff_smeared = (gaction_new_fwd_smeared - gaction_new_bwd_smeared) / 2ΔH - tc_symm_diff_smeared = (topcharge_new_fwd_smeared - topcharge_new_bwd_smeared) / 2ΔH - - # @show tc_symm_diff - # @show dtopcharge_proj - # @show tc_symm_diff_smeared - # @show dtopcharge_proj_smeared - # println("---") - relerrors[group_direction, 1] = (ga_symm_diff - dgaction_proj) / ga_symm_diff - relerrors[group_direction, 2] = - (ga_symm_diff_smeared - dgaction_proj_smeared) / ga_symm_diff_smeared - relerrors[group_direction, 3] = (tc_symm_diff - dtopcharge_proj) / tc_symm_diff - relerrors[group_direction, 4] = - (tc_symm_diff_smeared - dtopcharge_proj_smeared) / tc_symm_diff_smeared - - if mpi_amroot() - println("================= Group direction $(group_direction) =================") - println( - "/ GA rel. error (unsmeared): \t", (ga_symm_diff - dgaction_proj) / ga_symm_diff - ) - println( - "/ GA rel. error (smeared): \t", - (ga_symm_diff_smeared - dgaction_proj_smeared) / ga_symm_diff_smeared, - ) - println("") - println( - "/ TC rel. error (unsmeared): \t", - (tc_symm_diff - dtopcharge_proj) / tc_symm_diff, - ) - println( - "/ TC rel. error (smeared): \t", - (tc_symm_diff_smeared - dtopcharge_proj_smeared) / tc_symm_diff_smeared, - ) - end - - mpi_barrier() - end - # println() - return relerrors -end - -test_mpi_derivative() - diff --git a/test/test_mpi_fderivative.jl b/test/test_mpi_fderivative.jl deleted file mode 100644 index 47e7c376..00000000 --- a/test/test_mpi_fderivative.jl +++ /dev/null @@ -1,134 +0,0 @@ -using MetaQCD -using MetaQCD.Utils -using MetaQCD.Fields: update_halo! -using LinearAlgebra -using Random -using ProfileView - -function test_mpi_fderivative( - backend=CPU; dirac="staggered", mass=0.01, eoprec=false, single_flavor=false, csw=1.78 -) - if mpi_amroot() - println("Fermion derivative test [$dirac]") - end - - MetaQCD.MetaIO.set_global_logger!(1, nothing; tc=true) - Ns = 4 - Nt = 4 - U = Gaugefield{CPU,Float64,WilsonGaugeAction}(Ns, Ns, Ns, Nt, 5.7, (1, 1, 2, 2), 1) - random_gauges!(U) - - # filename = pkgdir(MetaQCD, "test", "testconf.txt") - # load_config!(BridgeFormat(), U, filename) - if backend !== CPU - U = MetaQCD.to_backend(backend, U) - end - - ψ = eoprec ? even_odd(Spinorfield(U; staggered=dirac=="staggered")) : Spinorfield(U; staggered=dirac=="staggered") - - spectral_bound, Nf = if dirac=="staggered" - (mass^2, 6.0), (single_flavor ? 1 : (eoprec ? 4 : 8)) - elseif dirac=="wilson" - (mass^2, 64.0), (single_flavor ? 1 : 2) - end - - params = ( - fermion_action=dirac, - eo_precon=eoprec, - boundary_condition="antiperiodic", - rhmc_spectral_bound=spectral_bound, - rhmc_order_md=15, - rhmc_prec_md=64, - rhmc_order_action=15, - rhmc_prec_action=64, - cg_tol_action=1e-16, - cg_tol_md=1e-16, - cg_maxiters_action=5000, - cg_maxiters_md=5000, - wilson_r=1, - wilson_csw=csw, - ) - - action = MetaQCD.DiracOperators.init_fermion_action(params, mass, Nf, U) - mpi_amroot() && (@show action) - - sample_pseudofermions!(ψ, action, U) - - # Test for smearing with 5 steps and stout parameter 0.12 - smearing = StoutSmearing(U, 5, 0.12) - - dSfdU = Colorfield(U) - dSfdU_smeared = Colorfield(U) - temp_force = Colorfield(U) - - site = SiteCoords(2, 3, 2, 2) - μ = 3 - ΔH = 0.000001 - - relerrors = Matrix{Float64}(undef, 8, 2) - - for group_direction in 1:8 - # Unsmeared - Ufwd = deepcopy(U) - if mpi_amroot() - Ufwd[μ, site] = expλ(group_direction, ΔH) * Ufwd[μ, site] - end - update_halo!(Ufwd) - action_new_fwd = calc_fermion_action(action, Ufwd, ψ) - - Ubwd = deepcopy(U) - if mpi_amroot() - Ubwd[μ, site] = expλ(group_direction, -ΔH) * Ubwd[μ, site] - end - update_halo!(Ubwd) - action_new_bwd = calc_fermion_action(action, Ubwd, ψ) - - # Smeared - Ufwd = deepcopy(U) - if mpi_amroot() - Ufwd[μ, site] = expλ(group_direction, ΔH) * Ufwd[μ, site] - end - update_halo!(Ufwd) - calc_smearedU!(smearing, Ufwd) - action_new_fwd_smeared = calc_fermion_action( - action, smearing.Usmeared_multi[end], ψ - ) - - Ubwd = deepcopy(U) - if mpi_amroot() - Ubwd[μ, site] = expλ(group_direction, -ΔH) * Ubwd[μ, site] - end - update_halo!(Ubwd) - calc_smearedU!(smearing, Ubwd) - action_new_bwd_smeared = calc_fermion_action( - action, smearing.Usmeared_multi[end], ψ - ) - - calc_dSfdU_bare!(dSfdU, action, U, ψ, nothing, NoSmearing()) - calc_dSfdU_bare!(dSfdU_smeared, action, U, ψ, temp_force, smearing) - - daction_proj = real(multr(im * λ[group_direction], dSfdU[μ, site])) - daction_proj_smeared = real(multr(im * λ[group_direction], dSfdU_smeared[μ, site])) - - symm_diff = (action_new_fwd - action_new_bwd) / 2ΔH - symm_diff_smeared = (action_new_fwd_smeared - action_new_bwd_smeared) / 2ΔH - - # if group_direction == 1 - # @show daction_proj - # @show symm_diff - # end - relerrors[group_direction, 1] = (symm_diff - daction_proj) / symm_diff - relerrors[group_direction, 2] = - (symm_diff_smeared - daction_proj_smeared) / symm_diff_smeared - - if mpi_amroot() - println("================= Group direction $(group_direction) =================") - println("/ Rel. error (unsmeared): \t", relerrors[group_direction, 1]) - println("/ Rel. error (smeared): \t", relerrors[group_direction, 2]) - end - end - # println() - return relerrors -end - -test_mpi_fderivative()