From 651891f95a771c0ddc2d3daaffbbd0bad3a60a4a Mon Sep 17 00:00:00 2001 From: Alberto Ramos Date: Wed, 13 Dec 2023 14:45:45 +0100 Subject: [PATCH] Added documentation for most modules Only Spinors and Dirac are missing. --- docs/make.jl | 6 ++++- docs/src/fields.md | 9 +++---- docs/src/flow.md | 47 ++++++++++++++++++++++++++++++++++++ docs/src/groups.md | 4 ++++ docs/src/io.md | 14 +++++++++++ docs/src/sf.md | 9 +++++++ docs/src/space.md | 5 +++- docs/src/ym.md | 59 ++++++++++++++++++++++++++++++++++++++++++++++ src/LatticeGPU.jl | 1 + src/MD/MD.jl | 20 ++++++++++++++++ src/YM/YM.jl | 36 ++++++++++++++++++++++++++-- src/YM/YMact.jl | 13 +++++++--- src/YM/YMfields.jl | 5 ++++ src/YM/YMflow.jl | 57 ++++++++++++++++++++++++++++++++++++++++---- src/YM/YMhmc.jl | 24 +++++++++++++++++-- src/YM/YMsf.jl | 8 +++++-- 16 files changed, 298 insertions(+), 19 deletions(-) create mode 100644 docs/src/flow.md create mode 100644 docs/src/io.md create mode 100644 docs/src/sf.md create mode 100644 docs/src/ym.md diff --git a/docs/make.jl b/docs/make.jl index 60b7af4..f52fe3e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -9,6 +9,10 @@ makedocs(sitename="LatticeGPU", modules=[LatticeGPU], doctest=true, "LatticeGPU.jl" => "index.md", "Space-time" => "space.md", "Groups and algebras" => "groups.md", - "Fields" => "fields.md" + "Fields" => "fields.md", + "Yang-Mills" => "ym.md", + "Gradient flow" => "flow.md", + "Schrödinger Functional" => "sf.md", + "Input Output" => "io.md" ], repo = "https://igit.ific.uv.es/alramos/latticegpu.jl") diff --git a/docs/src/fields.md b/docs/src/fields.md index b68e110..1dc42ab 100644 --- a/docs/src/fields.md +++ b/docs/src/fields.md @@ -10,10 +10,11 @@ gauge fields `SU3`, for scalar fields `Float64`). We have: direction. - `N` scalar fields: `N` elemental types at each spacetime point. -For all these fields the spacetime point are ordered in memory -according to the point-in-block and block indices (see -[`SpaceParm`](@ref)). An execption is the [`scalar_field_point`](@ref) -fields. +Fields can have **naturaL indexing**, where the memory layout follows +the point-in-block and block indices (see +[`SpaceParm`](@ref)). Fields can also have **lexicographic indexing**, +where points are labelled by a D-dimensional index (see [`scalar_field_point`](@ref)). + ## Initialization diff --git a/docs/src/flow.md b/docs/src/flow.md new file mode 100644 index 0000000..519afe7 --- /dev/null +++ b/docs/src/flow.md @@ -0,0 +1,47 @@ +# Gradient flow + +The gradient flow equations can be integrated in two different ways: +1. Using a fixed step-size integrator. In this approach one fixes the + step size $\epsilon$ and the links are evolved from + $V_\mu(t)$ to $V_\mu(t +\epsilon)$ using some integration + scheme. +1. Using an adaptive step-size integrator. In this approach one fixes + the tolerance $h$ and the links are evolved for a time $t_{\rm + end}$ (i.e. from $V_\mu(t)$ to $V_\mu(t +t_{\rm end})$) + with the condition that the maximum error while advancing is not + larger than $h$. + +In general adaptive step size integrators are much more efficient, but +one loses the possibility to measure flow quantities at the +intermediate times $\epsilon, 2\epsilon, 3\epsilon,...$. Adaptive +step size integrators are ideal for finite size scaling studies, while +a mix of both integrators is the most efficient approach in scale +setting applications. + +## Integration schemes + +```@docs +FlowIntr +wfl_euler +zfl_euler +wfl_rk2 +zfl_rk2 +wfl_rk3 +zfl_rk3 +``` + +## Integrating the flow equations + +```@docs +flw +flw_adapt +``` + +## Observables + +```@docs +Eoft_plaq +Eoft_clover +Qtop +``` + diff --git a/docs/src/groups.md b/docs/src/groups.md index 8426a1a..508ff32 100644 --- a/docs/src/groups.md +++ b/docs/src/groups.md @@ -153,6 +153,10 @@ projalg ## Generic `Algebra` methods ```@docs +dot +norm +norm2 +normalize exp expm alg2mat diff --git a/docs/src/io.md b/docs/src/io.md new file mode 100644 index 0000000..581e543 --- /dev/null +++ b/docs/src/io.md @@ -0,0 +1,14 @@ + +# Input/Output + +## Configurations + +Routines to read/write and import gauge configurations. +```@docs +read_cnfg +save_cnfg +import_bsfqcd +import_lex64 +import_cern64 +``` + diff --git a/docs/src/sf.md b/docs/src/sf.md new file mode 100644 index 0000000..9b9daad --- /dev/null +++ b/docs/src/sf.md @@ -0,0 +1,9 @@ +# Schödinger Functional + +Specific SF observables and routines + +```@docs +setbndfield +sfcoupling +``` + diff --git a/docs/src/space.md b/docs/src/space.md index 4db86c7..de3a478 100644 --- a/docs/src/space.md +++ b/docs/src/space.md @@ -3,7 +3,10 @@ D-dimensional lattice points are labeled by two ordered integer numbers: the point-in-block index ($$b$$ in the figure below) and the -block index ($$r$$ in the figure below). The routines [`up`](@ref) and +block index ($$r$$ in the figure below). This is called **natural +indexing**, in contrast with the **lexicographic indexing** where points on +the lattice are represented by a D-dimensional `CartesianIndex`. +The routines [`up`](@ref) and [`dw`](@ref) allow you to displace to the neighboring points of the lattice. ![D dimensional lattice points are labeled by its diff --git a/docs/src/ym.md b/docs/src/ym.md new file mode 100644 index 0000000..5886f32 --- /dev/null +++ b/docs/src/ym.md @@ -0,0 +1,59 @@ + +# Simulating Yang-Mills on the lattice + +```@docs +GaugeParm +YMworkspace +ztwist +``` + +## Gauge actions and forces + +Routines to compute the gauge action. +```@docs +gauge_action +``` + +Routines to compute the force derived from gauge actions. + +```@docs +force_gauge +``` + +### Force field refresh + +Algebra fields with **natural indexing** can be randomized. +```@docs +randomize! +``` + + +## Basic observables + +Some basic observable. +```@docs +plaquette +``` + +## HMC simulations + +### Integrating the EOM + +```@docs +IntrScheme +leapfrog +omf2 +omf4 +MD! +``` + +### HMC algorithm + +```@docs +hamiltonian +HMC! +``` + + + + diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl index 99306dd..0224489 100644 --- a/src/LatticeGPU.jl +++ b/src/LatticeGPU.jl @@ -40,6 +40,7 @@ include("YM/YM.jl") using .YM export ztwist export YMworkspace, GaugeParm, force0_wilson!, field, field_pln, randomize!, zero!, norm2 +export force_gauge, MD! export gauge_action, hamiltonian, plaquette, HMC!, OMF4! export Eoft_clover, Eoft_plaq, Qtop export FlowIntr, wfl_euler, zfl_euler, wfl_rk2, zfl_rk2, wfl_rk3, zfl_rk3 diff --git a/src/MD/MD.jl b/src/MD/MD.jl index 1b62775..f7560f2 100644 --- a/src/MD/MD.jl +++ b/src/MD/MD.jl @@ -24,6 +24,11 @@ const r1omf2 = 0.1931833275037836 const r2omf2 = 0.5 const r3omf2 = 1 - 2*r1omf2 +""" + struct IntrScheme{N, T} + +Integrator for the molecular dynamics. +""" struct IntrScheme{N, T} r::NTuple{N, T} eps::T @@ -31,8 +36,23 @@ struct IntrScheme{N, T} end +""" + omf2(::Type{T}, eps, ns) + +Second order Omelyan integrator with `eps` stepsize and `ns` steps. +""" omf2(::Type{T}, eps, ns) where T = IntrScheme{3,T}((r1omf2,r2omf2,r3omf2), eps, ns) +""" + omf4(::Type{T}, eps, ns) + +Fourth order Omelyan integrator with `eps` stepsize and `ns` steps. +""" omf4(::Type{T}, eps, ns) where T = IntrScheme{6,T}((r1omf4,r2omf4,r3omf4,r4omf4,r5omf4,r6omf4), eps, ns) +""" + leapfrog(::Type{T}, eps, ns) + +Leapfrog integrator with `eps` stepsize and `ns` steps. +""" leapfrog(::Type{T}, eps, ns) where T = IntrScheme{2,T}((0.5,1.0), eps, ns) diff --git a/src/YM/YM.jl b/src/YM/YM.jl index 490a434..bdfc098 100644 --- a/src/YM/YM.jl +++ b/src/YM/YM.jl @@ -20,6 +20,19 @@ using ..MD import Base.show +""" + struct GaugeParm{T,G,N} + +Structure containning the parameters of a pure gauge simulation. These are: +- beta: Type `T`. The bare coupling of the simulation +- c0: Type `T`. LatticeGPU supports the simulation of gauge actions made of 1x1 Wilson Loops and 2x1 Wilson loops. The parameter c0 defines the coefficient on the simulation of the 1x1 loops. Some common choices are: + - c0=1: Wilson plaquette action + - c0=5/3: Tree-level improved Lüscher-Weisz action. + - c0=3.648: Iwasaki gauge action +- cG: Tuple (`T`, `T`). Boundary improvement parameters. +- ng: `Int64`. Rank of the gauge group. +- Ubnd: Boundary field for SF boundary conditions +""" struct GaugeParm{T,G,N} beta::T c0::T @@ -63,6 +76,21 @@ function Base.show(io::IO, gp::GaugeParm{T, G, N}) where {T,G,N} return nothing end +""" + struct YMworkspace{T} + +Structure containing memory workspace that is resused by different routines in order to avoid allocating/deallocating time. +The parameter `T` represents the precision of the simulation (i.e. single/double). The structure contains the following components +- GRP: Group being simulated +- ALG: Corresponding Algebra +- PRC: Precision (i.e. `T`) +- frc1: Algebra field with natural indexing. +- frc2: Algebra field with natural indexing. +- mom: Algebra field with natural indexing. +- U1: Group field with natural indexing. +- cm: Complex field with lexicographic indexing. +- rm: Real field with lexicographic indexing. +""" struct YMworkspace{T} GRP ALG @@ -110,7 +138,11 @@ function Base.show(io::IO, ymws::YMworkspace) return nothing end +""" + function ztwist(gp::GaugeParm{T,G}, lp::SpaceParm{N,M,B,D}[, ipl]) +Returns the twist factor. If a plane index is passed, returns the twist factor as a complex{T}. If this is not provided, returns a tuple, containing the factor of each plane. +""" function ztwist(gp::GaugeParm{T,G}, lp::SpaceParm{N,M,B,D}) where {T,G,N,M,B,D} function plnf(ipl) @@ -133,10 +165,10 @@ include("YMfields.jl") export randomize!, zero!, norm2 include("YMact.jl") -export krnl_plaq!, force0_wilson! +export krnl_plaq!, force_gauge, force_wilson include("YMhmc.jl") -export gauge_action, hamiltonian, plaquette, HMC!, OMF4! +export gauge_action, hamiltonian, plaquette, HMC!, MD! include("YMflow.jl") export FlowIntr, flw, flw_adapt diff --git a/src/YM/YMact.jl b/src/YM/YMact.jl index c9bec0f..35e965b 100644 --- a/src/YM/YMact.jl +++ b/src/YM/YMact.jl @@ -335,9 +335,9 @@ function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, Ubnd, cG, end """ - function force_wilson(ymws::YMworkspace, U, lp::SpaceParm) + function force_gauge(ymws::YMworkspace, U, gp::GaugeParm, lp::SpaceParm) -Computes the force deriving from the Wilson plaquette action, without +Computes the force deriving from an improved action with parameter `c0`, without the prefactor 1/g0^2, and assign it to the workspace force `ymws.frc1` """ function force_gauge(ymws::YMworkspace, U, c0, cG, gp::GaugeParm, lp::SpaceParm) @@ -354,8 +354,15 @@ function force_gauge(ymws::YMworkspace, U, c0, cG, gp::GaugeParm, lp::SpaceParm) end return nothing end - force_gauge(ymws::YMworkspace, U, c0, gp, lp) = force_gauge(ymws, U, c0, gp.cG[1], gp, lp) +force_gauge(ymws::YMworkspace, U, gp, lp) = force_gauge(ymws, U, gp.c0, gp.cG[1], gp, lp) + +""" + function force_wilson(ymws::YMworkspace, U, gp::GaugeParm, lp::SpaceParm) + +Computes the force deriving from the Wilson plaquette action, without +the prefactor 1/g0^2, and assign it to the workspace force `ymws.frc1` +""" force_wilson(ymws::YMworkspace, U, gp::GaugeParm, lp::SpaceParm) = force_gauge(ymws, U, 1, gp, lp) force_wilson(ymws::YMworkspace, U, cG, gp::GaugeParm, lp::SpaceParm) = force_gauge(ymws, U, 1, gp.cG[1], gp, lp) diff --git a/src/YM/YMfields.jl b/src/YM/YMfields.jl index c6e5433..e22625a 100644 --- a/src/YM/YMfields.jl +++ b/src/YM/YMfields.jl @@ -9,6 +9,11 @@ ### created: Thu Jul 15 15:16:47 2021 ### +""" + function randomize!(f, lp::SpaceParm, ymws::YMworkspace) + +Given an algebra field with natural indexing, this routine sets the components to random Gaussian distributed values. If SF boundary conditions are used, the force at the boundaries is set to zero. +""" function randomize!(f, lp::SpaceParm, ymws::YMworkspace) if ymws.ALG == SU2alg diff --git a/src/YM/YMflow.jl b/src/YM/YMflow.jl index 8da9b35..0ab2926 100644 --- a/src/YM/YMflow.jl +++ b/src/YM/YMflow.jl @@ -10,6 +10,11 @@ ### +""" + struct FlowIntr{N,T} + +Structure containing info about a particular flow integrator +""" struct FlowIntr{N,T} r::T e0::NTuple{N,T} @@ -26,11 +31,46 @@ struct FlowIntr{N,T} end # pre-defined integrators +""" + wfl_euler(::Type{T}, eps::T, tol::T) + +Euler scheme integrator for the Wilson Flow. The fixed step size is given by `eps` and the tolerance for the adaptive integrators by `tol`. +""" wfl_euler(::Type{T}, eps::T, tol::T) where T = FlowIntr{0,T}(one(T),(),(),false,one(T),eps,tol,one(T)/200,one(T)/10,9/10) + +""" + zfl_euler(::Type{T}, eps::T, tol::T) + +Euler scheme integrator for the Zeuthen flow. The fixed step size is given by `eps` and the tolerance for the adaptive integrators by `tol`. +""" zfl_euler(::Type{T}, eps::T, tol::T) where T = FlowIntr{0,T}(one(T),(),(),true, (one(T)*5)/3,eps,tol,one(T)/200,one(T)/10,9/10) + +""" + wfl_rk2(::Type{T}, eps::T, tol::T) + +Second order Runge-Kutta integrator for the Wilson flow. The fixed step size is given by `eps` and the tolerance for the adaptive integrators by `tol`. +""" wfl_rk2(::Type{T}, eps::T, tol::T) where T = FlowIntr{1,T}(one(T)/2,(-one(T)/2,),(one(T),),false,one(T),eps,tol,one(T)/200,one(T)/10,9/10) + +""" + zfl_rk2(::Type{T}, eps::T, tol::T) + +Second order Runge-Kutta integrator for the Zeuthen flow. The fixed step size is given by `eps` and the tolerance for the adaptive integrators by `tol`. +""" zfl_rk2(::Type{T}, eps::T, tol::T) where T = FlowIntr{1,T}(one(T)/2,(-one(T)/2,),(one(T),),true, (one(T)*5)/3,eps,tol,one(T)/200,one(T)/10,9/10) + +""" + wfl_rk3(::Type{T}, eps::T, tol::T) + +Third order Runge-Kutta integrator for the Wilson flow. The fixed step size is given by `eps` and the tolerance for the adaptive integrators by `tol`. +""" wfl_rk3(::Type{T}, eps::T, tol::T) where T = FlowIntr{2,T}(one(T)/4,(-17/36,-one(T)),(8/9,3/4),false,one(T),eps,tol,one(T)/200,one(T)/10,9/10) + +""" + Zfl_rk3(::Type{T}, eps::T, tol::T) + +Third order Runge-Kutta integrator for the Zeuthen flow. The fixed step size is given by `eps` and the tolerance for the adaptive integrators by `tol`. +""" zfl_rk3(::Type{T}, eps::T, tol::T) where T = FlowIntr{2,T}(one(T)/4,(-17/36,-one(T)),(8/9,3/4),true, (one(T)*5)/3,eps,tol,one(T)/200,one(T)/10,9/10) function Base.show(io::IO, int::FlowIntr{N,T}) where {N,T} @@ -122,6 +162,11 @@ function krnl_add_zth!(frc, frc2::AbstractArray{TA}, U::AbstractArray{TG}, lp::S return nothing end +""" + function flw(U, int::FlowIntr{NI,T}, ns::Int64, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) + +Integrates the flow equations with the integration scheme defined by `int` performing `ns` steps with fixed step size. The configuration `U` is overwritten. +""" function flw(U, int::FlowIntr{NI,T}, ns::Int64, eps, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) where {NI,T} @timeit "Integrating flow equations" begin for i in 1:ns @@ -152,6 +197,11 @@ flw(U, int::FlowIntr{NI,T}, ns::Int64, gp::GaugeParm, lp::SpaceParm, ymws::YMwor # Adaptive step size integrators ## +""" + function flw_adapt(U, int::FlowIntr{NI,T}, tend::T, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) + +Integrates the flow equations with the integration scheme defined by `int` using the adaptive step size integrator up to `tend` with the tolerance defined in `int`. The configuration `U` is overwritten. +""" function flw_adapt(U, int::FlowIntr{NI,T}, tend::T, epsini::T, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) where {NI,T} eps = int.eps_ini @@ -201,7 +251,7 @@ flw_adapt(U, int::FlowIntr{NI,T}, tend::T, gp::GaugeParm, lp::SpaceParm, ymws::Y """ function Eoft_plaq([Eslc,] U, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) -Measure the action density `E(t)` using the plaquette discretization. If the argument `Eslc` +Measure the action density `E(t)` using the plaquette discretization. If the argument `Eslc` is given the contribution for each Euclidean time slice and plane are returned. """ function Eoft_plaq(Eslc, U, gp::GaugeParm{T,G,NN}, lp::SpaceParm{N,M,B,D}, ymws::YMworkspace) where {T,G,NN,N,M,B,D} @@ -277,10 +327,9 @@ function krnl_plaq_pln!(plx, U::AbstractArray{T}, Ubnd, ztw, ipl, lp::SpaceParm{ end """ - Qtop([Qslc,] U, lp, ymws) + Qtop([Qslc,] U, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace) -Measure the topological charge `Q` of the configuration `U`. If the argument `Qslc` is present -the contribution for each Euclidean time slice are returned. +Measure the topological charge `Q` of the configuration `U` using the clover definition of the field strength tensor. If the argument `Qslc` is present the contribution for each Euclidean time slice are returned. Only wors in 4D. """ function Qtop(Qslc, U, gp::GaugeParm, lp::SpaceParm{4,M,B,D}, ymws::YMworkspace) where {M,B,D} diff --git a/src/YM/YMhmc.jl b/src/YM/YMhmc.jl index c973875..57d2d22 100644 --- a/src/YM/YMhmc.jl +++ b/src/YM/YMhmc.jl @@ -13,7 +13,7 @@ function gauge_action(U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace) -Returns the value of the gauge plaquette action for the configuration U. The parameters `\beta` and `c0` are taken from the `gp` structure. +Returns the value of the gauge action for the configuration U. The parameters `\beta` and `c0` are taken from the `gp` structure. """ function gauge_action(U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}) where T <: AbstractFloat @@ -37,6 +37,11 @@ function gauge_action(U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}) whe return S end +""" + function plaquette(U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace) + +Computes the average plaquette for the configuration `U`. +""" function plaquette(U, lp::SpaceParm{N,M,B,D}, gp::GaugeParm, ymws::YMworkspace) where {N,M,B,D} ztw = ztwist(gp, lp) @@ -48,7 +53,12 @@ function plaquette(U, lp::SpaceParm{N,M,B,D}, gp::GaugeParm, ymws::YMworkspace) return CUDA.mapreduce(real, +, ymws.cm)/(prod(lp.iL)*lp.npls) end - + +""" + function hamiltonian(mom, U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace) + +Returns the Energy ``H = \\frac{p^2}{2}+S[U]``, where the momenta field is given by `mom` and the configuration by `U`. +""" function hamiltonian(mom, U, lp, gp, ymws) @timeit "Computing Hamiltonian" begin K = CUDA.mapreduce(norm2, +, mom)/2 @@ -58,6 +68,11 @@ function hamiltonian(mom, U, lp, gp, ymws) return K+V end +""" + HMC!(U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace; noacc=false) + +Performs a HMC step (molecular dynamics integration and accept/reject step). The configuration `U` is updated ans function returns the energy violation and if the configuration was accepted in a tuple. +""" function HMC!(U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}; noacc=false) where T @timeit "HMC trayectory" begin @@ -92,6 +107,11 @@ function HMC!(U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspac end HMC!(U, eps, ns, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}; noacc=false) where T = HMC!(U, omf4(T, eps, ns), lp, gp, ymws; noacc=noacc) +""" + function MD!(mom, U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace) + +Performs the integration of a molecular dynamics trajectory starting from the momentum field `mom` and the configuration `U` according to the integrator described by `int`. +""" function MD!(mom, U, int::IntrScheme{NI, T}, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}) where {NI, T <: AbstractFloat} @timeit "MD evolution" begin diff --git a/src/YM/YMsf.jl b/src/YM/YMsf.jl index 6ecc916..abd93d9 100644 --- a/src/YM/YMsf.jl +++ b/src/YM/YMsf.jl @@ -10,9 +10,9 @@ ### """ - sfcoupling(U, lp::SpaceParm{N,M,B,D}, gp::GaugeParm, ymws::YMworkspace) where {N,M,B,D} + sfcoupling(U, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace) -Measures the Schrodinger Functional coupling `ds/d\eta` and `d^2S/d\eta d\nu`. +Measures the Schrodinger Functional coupling ``{\\rm d}S/{\\rm d}\\eta`` and ``{\\rm d}^2S/{\\rm d}\\eta d\nu``. """ function sfcoupling(U, lp::SpaceParm{N,M,B,D}, gp::GaugeParm, ymws::YMworkspace) where {N,M,B,D} @@ -89,7 +89,11 @@ end return exp(X) end +""" + function setbndfield(U, phi, lp::SpaceParm) +Sets abelian boundary fields with phases `phi[1]` and `phi[2]` to the configuration `U` at time salice ``x_0=0``. +""" function setbndfield(U, phi, lp::SpaceParm{N,M,B,D}) where {N,M,B,D} CUDA.@sync begin