(Almost) Done with distributed fields. EO-Precon'ed Wilson-Clover sti…

…ll elusive

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)

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

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

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)

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)

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,
             """

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)

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)`
 

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)`