diff --git a/.travis.yml b/.travis.yml index b27503df..c08c9cab 100644 --- a/.travis.yml +++ b/.travis.yml @@ -13,8 +13,8 @@ matrix: - julia: nightly after_success: - - julia --project -e 'using Pkg; Pkg.add("Coverage"); Pkg.dir("GeophysicalFlows"); using Coverage; Coveralls.submit(Coveralls.process_folder())' - - julia --project -e 'using Pkg; Pkg.add("Coverage"); Pkg.dir("GeophysicalFlows"); using Coverage; Codecov.submit(Codecov.process_folder())' + - julia --project -e 'using Pkg; Pkg.add("Coverage"); import GeophysicalFlows; joinpath(dirname(pathof(GeophysicalFlows)), ".."); using Coverage; Coveralls.submit(Coveralls.process_folder());' + - julia --project -e 'using Pkg; Pkg.add("Coverage"); import GeophysicalFlows; joinpath(dirname(pathof(GeophysicalFlows)), ".."); using Coverage; Codecov.submit(Codecov.process_folder())' jobs: include: diff --git a/Manifest.toml b/Manifest.toml index f3cfe932..a4a68be2 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -6,6 +6,12 @@ git-tree-sha1 = "380e36c66edfa099cd90116b24c1ce8cafccac40" uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c" version = "0.4.1" +[[Adapt]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "82dab828020b872fa9efd3abec1152b075bc7cbf" +uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +version = "1.0.0" + [[AxisAlgorithms]] deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"] git-tree-sha1 = "a4d07a1c313392a77042855df46c5f534076fab9" @@ -27,12 +33,35 @@ git-tree-sha1 = "c7361ce8a2129f20b0e05a89f7070820cfed6648" uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232" version = "0.5.6" +[[CEnum]] +git-tree-sha1 = "62847acab40e6855a9b5905ccb99c2b5cf6b3ebb" +uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" +version = "0.2.0" + [[CSTParser]] deps = ["Tokenize"] git-tree-sha1 = "c69698c3d4a7255bc1b4bc2afc09f59db910243b" uuid = "00ebfdb7-1f24-5e51-bd34-a7502290713f" version = "0.6.2" +[[CUDAapi]] +deps = ["Libdl", "Logging"] +git-tree-sha1 = "9b2b4b71d6b7f946c9689bb4dea03ff92e3c7091" +uuid = "3895d2a7-ec45-59b8-82bb-cfc6a382f9b3" +version = "1.1.0" + +[[CUDAdrv]] +deps = ["CUDAapi", "Libdl", "Printf"] +git-tree-sha1 = "9ce99b5732c70e06ed97c042187baed876fb1698" +uuid = "c5f51814-7f29-56b8-a69c-e4d8f6be1fde" +version = "3.1.0" + +[[CUDAnative]] +deps = ["Adapt", "CUDAapi", "CUDAdrv", "DataStructures", "InteractiveUtils", "LLVM", "Libdl", "Logging", "Printf", "TimerOutputs"] +git-tree-sha1 = "52ae1ce10ebfa686e227655c47b19add89308623" +uuid = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17" +version = "2.3.1" + [[CodecZlib]] deps = ["BinaryProvider", "Libdl", "TranscodingStreams"] git-tree-sha1 = "05916673a2627dd91b4969ff8ba6941bc85a960e" @@ -69,6 +98,18 @@ git-tree-sha1 = "f9780daf3fea51ad2d7b7ed4f480ac34ab36587f" uuid = "a2441757-f6aa-5fb2-8edb-039e3f45d037" version = "0.9.2" +[[Crayons]] +deps = ["Test"] +git-tree-sha1 = "f621b8ef51fd2004c7cf157ea47f027fdeac5523" +uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" +version = "4.0.0" + +[[CuArrays]] +deps = ["AbstractFFTs", "Adapt", "CUDAapi", "CUDAdrv", "CUDAnative", "GPUArrays", "LinearAlgebra", "MacroTools", "NNlib", "Printf", "Random", "Requires", "SparseArrays", "TimerOutputs"] +git-tree-sha1 = "46b48742a84bb839e74215b7e468a4a1c6ba30f9" +uuid = "3a865a2d-5b23-5a0f-bc46-62713ec82fae" +version = "1.2.1" + [[DataStructures]] deps = ["InteractiveUtils", "OrderedCollections"] git-tree-sha1 = "0809951a1774dc724da22d26e4289bbaab77809a" @@ -88,10 +129,10 @@ deps = ["Random", "Serialization", "Sockets"] uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" [[FFTW]] -deps = ["AbstractFFTs", "BinaryProvider", "Compat", "Conda", "Libdl", "LinearAlgebra", "Reexport", "Test"] -git-tree-sha1 = "29cda58afbf62f35b1a094882ad6c745a47b2eaa" +deps = ["AbstractFFTs", "BinaryProvider", "Conda", "Libdl", "LinearAlgebra", "Reexport", "Test"] +git-tree-sha1 = "e1a479d3c972f20c9a70563eec740bbfc786f515" uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" -version = "0.2.4" +version = "0.3.0" [[FileIO]] deps = ["Pkg"] @@ -99,22 +140,34 @@ git-tree-sha1 = "351f001a78aa1b7ad2696e386e110b5abd071c71" uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" version = "1.0.7" +[[FillArrays]] +deps = ["LinearAlgebra", "Random", "SparseArrays"] +git-tree-sha1 = "8fba6ddaf66b45dec830233cea0aae43eb1261ad" +uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" +version = "0.6.4" + [[FixedPointNumbers]] git-tree-sha1 = "d14a6fa5890ea3a7e5dcab6811114f132fec2b4b" uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" version = "0.6.1" [[FourierFlows]] -deps = ["Coverage", "FFTW", "Interpolations", "JLD2", "LinearAlgebra", "Printf", "Random", "Reexport", "Requires", "Statistics"] -git-tree-sha1 = "e430df633d39e02e772856e6cb745b99b069de19" +deps = ["CUDAapi", "Coverage", "CuArrays", "FFTW", "Interpolations", "JLD2", "LinearAlgebra", "Printf", "Random", "Reexport", "Requires", "Statistics"] +git-tree-sha1 = "153b295514c016e30b300b9e65b0615061976166" uuid = "2aec4490-903f-5c70-9b11-9bed06a700e1" -version = "0.3.1" +version = "0.3.2" + +[[GPUArrays]] +deps = ["Adapt", "FFTW", "FillArrays", "LinearAlgebra", "Printf", "Random", "Serialization", "StaticArrays", "Test"] +git-tree-sha1 = "b5009ac44b141ded5e6f04c4db83807970f56e91" +uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7" +version = "1.0.2" [[HTTP]] deps = ["Base64", "Dates", "IniFile", "MbedTLS", "Sockets"] -git-tree-sha1 = "c4a527dba1d26add0e85946e1a53f42a1b343acc" +git-tree-sha1 = "f1e1b417a34cf73a70cbed919915bf8f8bad1806" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" -version = "0.8.5" +version = "0.8.6" [[IniFile]] deps = ["Test"] @@ -144,6 +197,12 @@ git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e" uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" version = "0.21.0" +[[LLVM]] +deps = ["CEnum", "Libdl", "Printf", "Unicode"] +git-tree-sha1 = "4a05f742837779a00bd8c9a18da6817367c4245d" +uuid = "929cbde3-209d-540e-8aea-75f648917ca0" +version = "1.3.0" + [[LaTeXStrings]] deps = ["Compat"] git-tree-sha1 = "7ab9b8788cfab2bdde22adf9004bda7ad9954b6c" @@ -182,6 +241,12 @@ version = "0.7.0" [[Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" +[[NNlib]] +deps = ["Libdl", "LinearAlgebra", "Requires", "Statistics", "TimerOutputs"] +git-tree-sha1 = "0c667371391fc6bb31f7f12f96a56a17098b3de8" +uuid = "872c559c-99b0-510c-b3b7-b6c96a88d5cd" +version = "0.6.0" + [[OffsetArrays]] git-tree-sha1 = "1af2f79c7eaac3e019a0de41ef63335ff26a0a57" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" @@ -263,10 +328,10 @@ deps = ["LinearAlgebra", "Random"] uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [[SpecialFunctions]] -deps = ["BinDeps", "BinaryProvider", "Libdl", "Test"] -git-tree-sha1 = "0b45dc2e45ed77f445617b99ff2adf0f5b0f23ea" +deps = ["BinDeps", "BinaryProvider", "Libdl"] +git-tree-sha1 = "3bdd374b6fd78faf0119b8c5d538788dbf910c6e" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "0.7.2" +version = "0.8.0" [[StaticArrays]] deps = ["LinearAlgebra", "Random", "Statistics"] @@ -282,6 +347,12 @@ uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" deps = ["Distributed", "InteractiveUtils", "Logging", "Random"] uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +[[TimerOutputs]] +deps = ["Crayons", "Printf", "Test", "Unicode"] +git-tree-sha1 = "b80671c06f8f8bae08c55d67b5ce292c5ae2660c" +uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +version = "0.5.0" + [[Tokenize]] git-tree-sha1 = "dfcdbbfb2d0370716c815cbd6f8a364efb6f42cf" uuid = "0796e94c-ce3b-5d07-9a54-7f471281c624" diff --git a/Project.toml b/Project.toml index 316f3478..29deccfe 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.3.0" [deps] Coverage = "a2441757-f6aa-5fb2-8edb-039e3f45d037" +CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" FourierFlows = "2aec4490-903f-5c70-9b11-9bed06a700e1" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" @@ -17,12 +18,14 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +[compat] +FFTW = "≥ 0.2.4" +FourierFlows = "≥ 0.3.2" +JLD2 = "≥ 0.1.2" +julia = "≥ 1.0.0" + [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["Test"] - -[compat] -FourierFlows = "≥ 0.3.0" -julia = "≥ 1.0.0" diff --git a/examples/twodturb_mcwilliams1984.jl b/examples/twodturb_mcwilliams1984.jl index 3a4e3816..71be770c 100644 --- a/examples/twodturb_mcwilliams1984.jl +++ b/examples/twodturb_mcwilliams1984.jl @@ -6,14 +6,16 @@ import GeophysicalFlows.TwoDTurb import GeophysicalFlows.TwoDTurb: energy, enstrophy import GeophysicalFlows: peakedisotropicspectrum +dev = CPU() # Device (CPU/GPU) + # Parameters n = 256 L = 2π -nnu = 2 - nu = 0.0 + nν = 2 + ν = 0.0 dt = 1e-2 nsteps = 5000 -nsubs = 200 + nsubs = 200 # Files filepath = "." @@ -26,7 +28,7 @@ if isfile(filename); rm(filename); end if !isdir(plotpath); mkdir(plotpath); end # Initialize problem -prob = TwoDTurb.Problem(; nx=n, Lx=L, ny=n, Ly=L, nu=nu, nnu=nnu, dt=dt, stepper="FilteredRK4") +prob = TwoDTurb.Problem(; nx=n, Lx=L, ny=n, Ly=L, ν=ν, nν=nν, dt=dt, stepper="FilteredRK4", dev=dev) sol, cl, vs, gr, filter = prob.sol, prob.clock, prob.vars, prob.grid, prob.timestepper.filter x, y = gridpoints(gr) @@ -35,7 +37,7 @@ x, y = gridpoints(gr) # that reproduces the results of the paper by McWilliams (1984) seed!(1234) k0, E0 = 6, 0.5 -zetai = peakedisotropicspectrum(gr, k0, E0, mask=filter) +zetai = peakedisotropicspectrum(gr, k0, E0, mask=filter) TwoDTurb.set_zeta!(prob, zetai) # Create Diagnostic -- energy and enstrophy are functions imported at the top. @@ -45,8 +47,8 @@ diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will # be updated every timestep. # Create Output -get_sol(prob) = prob.sol # extracts the Fourier-transformed solution -get_u(prob) = irfft(im*gr.l.*gr.invKrsq.*sol, gr.nx) +get_sol(prob) = Array(prob.sol) # extracts the Fourier-transformed solution +get_u(prob) = Array(irfft(im*gr.l.*gr.invKrsq.*sol, gr.nx)) out = Output(prob, filename, (:sol, get_sol), (:u, get_u)) saveproblem(out) @@ -89,7 +91,7 @@ while cl.step < nsteps println(log) plot_output(prob, fig, axs; drawcolorbar=false) end - +println("finished") plot_output(prob, fig, axs; drawcolorbar=false) savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), cl.step) diff --git a/examples/twodturb_randomdecay.jl b/examples/twodturb_randomdecay.jl index f6f1aa03..c7c696e1 100644 --- a/examples/twodturb_randomdecay.jl +++ b/examples/twodturb_randomdecay.jl @@ -7,18 +7,19 @@ using FFTW: rfft import GeophysicalFlows.TwoDTurb - -nx = 256 # Resolution -Lx = 2π # Domain size -nu = 1e-6 # Viscosity -nnu = 1 # Order of (hyper-)viscosity. nnu=1 means Laplacian -dt = 0.1 # Timestep -nint = 200 # Number of steps between plots -ntot = 10nint # Number of total timesteps - + dev = CPU() # Device (CPU/GPU) + + nx = 256 # Resolution + Lx = 2π # Domain size + ν = 1e-6 # Viscosity + nν = 1 # Order of (hyper-)viscosity. nν=1 means Laplacian + dt = 0.1 # Timestep +nint = 200 # Number of steps between plots +ntot = 10nint # Number of total timesteps + # Define problem -prob = TwoDTurb.Problem(nx=nx, Lx=Lx, nu=nu, nnu=nnu, dt=dt, stepper="FilteredRK4") -TwoDTurb.set_zeta!(prob, rand(nx, nx)) +prob = TwoDTurb.Problem(nx=nx, Lx=Lx, ν=ν, nν=nν, dt=dt, stepper="FilteredRK4", dev=dev) +TwoDTurb.set_zeta!(prob, ArrayType(dev)(rand(nx, nx))) cl, vs, gr = prob.clock, prob.vars, prob.grid x, y = gridpoints(gr) @@ -49,7 +50,7 @@ while cl.step < ntot end # Plot the radial energy spectrum -E = @. 0.5*(vs.u^2 + vs.v^2) # energy density +E = @. 0.5*(vs.u^2 + vs.v^2) # energy density Eh = rfft(E) kr, Ehr = FourierFlows.radialspectrum(Eh, gr, refinement=1) diff --git a/examples/twodturb_stochasticforcing.jl b/examples/twodturb_stochasticforcing.jl index 687f7bdb..a3ca7566 100644 --- a/examples/twodturb_stochasticforcing.jl +++ b/examples/twodturb_stochasticforcing.jl @@ -6,21 +6,23 @@ using Printf: @printf import GeophysicalFlows.TwoDTurb import GeophysicalFlows.TwoDTurb: energy, enstrophy, dissipation, work, drag - n, L = 256, 2π -nu, nnu = 1e-7, 2 -mu, nmu = 1e-1, 0 -dt, tf = 0.005, 0.2/mu -nt = round(Int, tf/dt) -ns = 4 + dev = CPU() # Device (CPU/GPU) + + n, L = 256, 2π + ν, nν = 1e-7, 2 + μ, nμ = 1e-1, 0 +dt, tf = 0.005, 0.2/μ + nt = round(Int, tf/dt) + ns = 4 # Forcing kf, dkf = 12.0, 2.0 # forcing central wavenumber, wavenumber width ε = 0.1 # energy injection rate -gr = TwoDGrid(n, L) +gr = TwoDGrid(dev, n, L) x, y = gridpoints(gr) -Kr = [ gr.kr[i] for i=1:gr.nkr, j=1:gr.nl] +Kr = ArrayType(dev)([ gr.kr[i] for i=1:gr.nkr, j=1:gr.nl]) force2k = @. exp(-(sqrt(gr.Krsq)-kf)^2/(2*dkf^2)) force2k[gr.Krsq .< 2.0^2 ] .= 0 @@ -32,14 +34,14 @@ force2k .= ε/ε0 * force2k seed!(1234) function calcF!(Fh, sol, t, cl, v, p, g) - eta = exp.(2π*im*rand(Float64, size(sol)))/sqrt(cl.dt) + eta = ArrayType(dev)(exp.(2π*im*rand(typeof(gr.Lx), size(sol)))/sqrt(cl.dt)) eta[1, 1] = 0 @. Fh = eta*sqrt(force2k) nothing end -prob = TwoDTurb.Problem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, nmu=nmu, dt=dt, stepper="RK4", - calcF=calcF!, stochastic=true) +prob = TwoDTurb.Problem(nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt, stepper="RK4", + calcF=calcF!, stochastic=true, dev=dev) sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid @@ -55,7 +57,7 @@ function makeplot(prob, diags) TwoDTurb.updatevars!(prob) E, D, W, R = diags - t = round(mu*cl.t, digits=2) + t = round(μ*cl.t, digits=2) sca(axs[1]); cla() pcolormesh(x, y, v.zeta) xlabel(L"$x$") @@ -80,25 +82,25 @@ function makeplot(prob, diags) # If the Ito interpretation was used for the work # then we need to add the drift term: I[ii2] + ε - plot(mu*E.t[ii], W[ii2], label=L"work ($W$)") # Ito - # plot(mu*E.t[ii], W[ii2] , label=L"work ($W$)") # Stratonovich - plot(mu*E.t[ii], ε .+ 0*E.t[ii], "--", label=L"ensemble mean work ($\langle W\rangle $)") - # plot(mu*E.t[ii], -D[ii], label="dissipation (\$D\$)") - plot(mu*E.t[ii], -R[ii], label=L"drag ($D=2\mu E$)") - plot(mu*E.t[ii], 0*E.t[ii], "k:", linewidth=0.5) + plot(μ*E.t[ii], W[ii2], label=L"work ($W$)") # Ito + # plot(μ*E.t[ii], W[ii2] , label=L"work ($W$)") # Stratonovich + plot(μ*E.t[ii], ε .+ 0*E.t[ii], "--", label=L"ensemble mean work ($\langle W\rangle $)") + # plot(μ*E.t[ii], -D[ii], label="dissipation (\$D\$)") + plot(μ*E.t[ii], -R[ii], label=L"drag ($D=2\mu E$)") + plot(μ*E.t[ii], 0*E.t[ii], "k:", linewidth=0.5) ylabel("Energy sources and sinks") xlabel(L"$\mu t$") legend(fontsize=10) sca(axs[2]); cla() - plot(mu*E.t[ii], total[ii], label=L"computed $W-D$") - plot(mu*E.t[ii], dEdt, "--k", label=L"numerical $dE/dt$") + plot(μ*E.t[ii], total[ii], label=L"computed $W-D$") + plot(μ*E.t[ii], dEdt, "--k", label=L"numerical $dE/dt$") ylabel(L"$dE/dt$") xlabel(L"$\mu t$") legend(fontsize=10) sca(axs[4]); cla() - plot(mu*E.t[ii], residual, "c-", label=L"residual $dE/dt$ = computed $-$ numerical") + plot(μ*E.t[ii], residual, "c-", label=L"residual $dE/dt$ = computed $-$ numerical") xlabel(L"$\mu t$") legend(fontsize=10) diff --git a/src/twodturb.jl b/src/twodturb.jl index cc39fbdc..44063aa4 100644 --- a/src/twodturb.jl +++ b/src/twodturb.jl @@ -22,11 +22,6 @@ using FourierFlows: getfieldspecs, structvarsexpr, parsevalsum, parsevalsum2 abstract type TwoDTurbVars <: AbstractVars end -const physicalvars = [:zeta, :u, :v] -const transformvars = [ Symbol(var, :h) for var in physicalvars ] -const forcedvars = [:Fh] -const stochforcedvars = [:prevsol] - nothingfunction(args...) = nothing """ @@ -42,22 +37,23 @@ function Problem(; Ly = Lx, dt = 0.01, # Drag and/or hyper-/hypo-viscosity - nu = 0, - nnu = 1, - mu = 0, - nmu = 0, + ν = 0, + nν = 1, + μ = 0, + nμ = 0, # Timestepper and eqn options stepper = "RK4", calcF = nothingfunction, stochastic = false, - T = Float64 + T = Float64, + dev = CPU() ) - gr = TwoDGrid(nx, Lx, ny, Ly; T=T) - pr = Params{T}(nu, nnu, mu, nmu, calcF) - vs = calcF == nothingfunction ? Vars(gr) : (stochastic ? StochasticForcedVars(gr) : ForcedVars(gr)) + gr = TwoDGrid(dev, nx, Lx, ny, Ly; T=T) + pr = Params{T}(ν, nν, μ, nμ, calcF) + vs = calcF == nothingfunction ? Vars(dev, gr) : (stochastic ? StochasticForcedVars(dev, gr) : ForcedVars(dev, gr)) eq = Equation(pr, gr) - FourierFlows.Problem(eq, stepper, dt, gr, vs, pr) + FourierFlows.Problem(eq, stepper, dt, gr, vs, pr, dev) end @@ -66,18 +62,18 @@ end # ---------- """ - Params(nu, nnu, mu, nmu, calcF!) + Params(ν, nν, μ, nμ, calcF!) Returns the params for two-dimensional turbulence. """ struct Params{T} <: AbstractParams - nu::T # Vorticity viscosity - nnu::Int # Vorticity hyperviscous order - mu::T # Bottom drag or hypoviscosity - nmu::Int # Order of hypodrag - calcF!::Function # Function that calculates the forcing F + ν :: T # Vorticity viscosity + nν :: Int # Vorticity hyperviscous order + μ :: T # Bottom drag or hypoviscosity + nμ :: Int # Order of hypodrag + calcF! :: Function # Function that calculates the forcing F end -Params(nu, nnu) = Params(nu, nnu, typeof(nu)(0), 0, nothingfunction) +Params(ν, nν) = Params(ν, nν, typeof(ν)(0), 0, nothingfunction) # --------- @@ -90,7 +86,7 @@ Params(nu, nnu) = Params(nu, nnu, typeof(nu)(0), 0, nothingfunction) Returns the equation for two-dimensional turbulence with params p and grid g. """ function Equation(p::Params, g::AbstractGrid{T}) where T - L = @. - p.nu*g.Krsq^p.nnu - p.mu*g.Krsq^p.nmu + L = @. - p.ν*g.Krsq^p.nν - p.μ*g.Krsq^p.nμ L[1, 1] = 0 FourierFlows.Equation(L, calcN!, g) end @@ -100,51 +96,54 @@ end # Vars # ---- -varspecs = cat( - getfieldspecs(physicalvars, :(Array{T,2})), - getfieldspecs(transformvars, :(Array{Complex{T},2})), - dims=1) - -forcedvarspecs = cat(varspecs, getfieldspecs(forcedvars, :(Array{Complex{T},2})), dims=1) -stochforcedvarspecs = cat(forcedvarspecs, getfieldspecs(stochforcedvars, :(Array{Complex{T},2})), dims=1) +struct Vars{Aphys, Atrans, F, P} <: TwoDTurbVars + zeta :: Aphys + u :: Aphys + v :: Aphys + zetah :: Atrans + uh :: Atrans + vh :: Atrans + Fh :: F + prevsol :: P +end -# Construct Vars types -eval(structvarsexpr(:Vars, varspecs; parent=:TwoDTurbVars)) -eval(structvarsexpr(:ForcedVars, forcedvarspecs; parent=:TwoDTurbVars)) -eval(structvarsexpr(:StochasticForcedVars, stochforcedvarspecs; parent=:TwoDTurbVars)) +const ForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, Nothing} +const StochasticForcedVars = Vars{<:AbstractArray, <:AbstractArray, <:AbstractArray, <:AbstractArray} """ - Vars(g) + Vars(dev, g) -Returns the vars for unforced two-dimensional turbulence with grid g. +Returns the vars for unforced two-dimensional turbulence on device dev and with +grid g. """ -function Vars(g; T=typeof(g.Lx)) - @createarrays T (g.nx, g.ny) zeta u v - @createarrays Complex{T} (g.nkr, g.nl) sol zetah uh vh - Vars(zeta, u, v, zetah, uh, vh) +function Vars(::Dev, g::AbstractGrid{T}) where {Dev, T} + @devzeros Dev T (g.nx, g.ny) zeta u v + @devzeros Dev Complex{T} (g.nkr, g.nl) zetah uh vh + Vars(zeta, u, v, zetah, uh, vh, nothing, nothing) end """ - ForcedVars(g) + ForcedVars(dev, g) -Returns the vars for forced two-dimensional turbulence with grid g. +Returns the vars for forced two-dimensional turbulence on device dev and with +grid g. """ -function ForcedVars(g; T=typeof(g.Lx)) - v = Vars(g; T=T) - Fh = zeros(Complex{T}, (g.nkr, g.nl)) - ForcedVars(getfield.(Ref(v), fieldnames(typeof(v)))..., Fh) +function ForcedVars(dev::Dev, g::AbstractGrid{T}) where {Dev, T} + @devzeros Dev T (g.nx, g.ny) zeta u v + @devzeros Dev Complex{T} (g.nkr, g.nl) zetah uh vh Fh + Vars(zeta, u, v, zetah, uh, vh, Fh, nothing) end """ - StochasticForcedVars(g; T) + StochasticForcedVars(dev, g) -Returns the vars for stochastically forced two-dimensional turbulence with grid -g. +Returns the vars for stochastically forced two-dimensional turbulence on device +dev and with grid g. """ -function StochasticForcedVars(g; T=typeof(g.Lx)) - v = ForcedVars(g; T=T) - prevsol = zeros(Complex{T}, (g.nkr, g.nl)) - StochasticForcedVars(getfield.(Ref(v), fieldnames(typeof(v)))..., prevsol) +function StochasticForcedVars(dev::Dev, g::AbstractGrid{T}) where {Dev, T} + @devzeros Dev T (g.nx, g.ny) zeta u v + @devzeros Dev Complex{T} (g.nkr, g.nl) zetah uh vh Fh prevsol + Vars(zeta, u, v, zetah, uh, vh, Fh, prevsol) end @@ -158,8 +157,8 @@ end Calculates the advection term. """ function calcN_advection!(N, sol, t, cl, v, p, g) - @. v.uh = im * g.l * g.invKrsq * sol - @. v.vh = -im * g.kr * g.invKrsq * sol + @. v.uh = im * g.l * g.invKrsq * sol + @. v.vh = - im * g.kr * g.invKrsq * sol @. v.zetah = sol ldiv!(v.u, g.rfftplan, v.uh) @@ -172,7 +171,7 @@ function calcN_advection!(N, sol, t, cl, v, p, g) mul!(v.uh, g.rfftplan, v.u) # \hat{u*zeta} mul!(v.vh, g.rfftplan, v.v) # \hat{v*zeta} - @. N = -im*g.kr*v.uh - im*g.l*v.vh + @. N = - im*g.kr*v.uh - im*g.l*v.vh nothing end @@ -211,9 +210,9 @@ Update the vars in v on the grid g with the solution in sol. """ function updatevars!(prob) v, g, sol = prob.vars, prob.grid, prob.sol - v.zetah .= sol - @. v.uh = im * g.l * g.invKrsq * sol - @. v.vh = -im * g.kr * g.invKrsq * sol + @. v.zetah = sol + @. v.uh = im * g.l * g.invKrsq * sol + @. v.vh = - im * g.kr * g.invKrsq * sol ldiv!(v.zeta, g.rfftplan, deepcopy(v.zetah)) ldiv!(v.u, g.rfftplan, deepcopy(v.uh)) ldiv!(v.v, g.rfftplan, deepcopy(v.vh)) @@ -260,13 +259,13 @@ end """ dissipation(prob) -Returns the domain-averaged dissipation rate. nnu must be >= 1. +Returns the domain-averaged dissipation rate. nν must be >= 1. """ @inline function dissipation(prob) sol, v, p, g = prob.sol, prob.vars, prob.params, prob.grid - @. v.uh = g.Krsq^(p.nnu-1) * abs2(sol) + @. v.uh = g.Krsq^(p.nν-1) * abs2(sol) v.uh[1, 1] = 0 - p.nu/(g.Lx*g.Ly)*parsevalsum(v.uh, g) + p.ν/(g.Lx*g.Ly)*parsevalsum(v.uh, g) end """ @@ -281,7 +280,7 @@ Returns the domain-averaged rate of work of energy by the forcing Fh. end @inline function work(sol, v::StochasticForcedVars, g) - @. v.uh = g.invKrsq * (v.prevsol + sol)/2.0 * conj(v.Fh) # Stratonovich + @. v.uh = g.invKrsq * (v.prevsol + sol)/2 * conj(v.Fh) # Stratonovich # @. v.uh = g.invKrsq * v.prevsol * conj(v.Fh) # Ito 1/(g.Lx*g.Ly)*parsevalsum(v.uh, g) end @@ -291,13 +290,13 @@ end """ drag(prob) -Returns the extraction of domain-averaged energy by drag/hypodrag mu. +Returns the extraction of domain-averaged energy by drag/hypodrag μ. """ @inline function drag(prob) sol, v, p, g = prob.sol, prob.vars, prob.params, prob.grid - @. v.uh = g.Krsq^(p.nmu-1) * abs2(sol) + @. v.uh = g.Krsq^(p.nμ-1) * abs2(sol) v.uh[1, 1] = 0 - p.mu/(g.Lx*g.Ly)*parsevalsum(v.uh, g) + p.μ/(g.Lx*g.Ly)*parsevalsum(v.uh, g) end end # module diff --git a/src/utils.jl b/src/utils.jl index b975afe6..70fe83b7 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -4,14 +4,15 @@ Return the 2D vorticity field of the Lamb dipole with strength `U` and radius `R`, centered on `center=(xc, yc)` and on the grid `g`. The default value of `center` is the middle of the grid. """ -function lambdipole(U, R, g; center=(mean(g.x), mean(g.y))) +function lambdipole(U, R, g::AbstractGrid{T, A}; center=(mean(g.x), mean(g.y))) where {T, A} firstzero = 3.8317059702075123156 k = firstzero/R # dipole wavenumber for radius R in terms of first zero of besselj q0 = -2U*k/besselj(0, k*R) # dipole amplitude for strength U and radius R xc, yc = center r = @. sqrt( (g.x-xc)^2 + (g.y-yc)^2 ) - q = @. q0 * besselj(1, k*r) * (g.y-yc)/r + besselj1 = A([besselj(1, k*r[i, j]) for i=1:g.nx, j=1:g.ny]) + q = @. q0 * besselj1 * (g.y-yc)/r @. q[r >= R] = 0 q end @@ -22,20 +23,23 @@ Generate a real and random two-dimensional vorticity field q(x, y) with a Fourier spectrum peaked around a central non-dimensional wavenumber kpeak and normalized so that its total energy is E0. """ -function peakedisotropicspectrum(g::TwoDGrid, kpeak::Real, E0::Real; mask=ones(size(g.Krsq)), allones=false) +function peakedisotropicspectrum(g::AbstractGrid{T, A}, kpeak::Real, E0::Real; mask=ones(size(g.Krsq)), allones=false) where {T, A} if g.Lx !== g.Ly - error("the domain is not square") + error("the domain is not square") else k0 = kpeak*2π/g.Lx modk = sqrt.(g.Krsq) - psik = zeros(g.nk, g.nl) + psik = A(zeros(T, (g.nk, g.nl))) psik = @. (modk^2 * (1 + (modk/k0)^4))^(-0.5) psik[1, 1] = 0.0 - psih = @. randn(Complex{typeof(g.Lx)}, size(g.Krsq))*psik + phases = randn(Complex{T}, size(g.Krsq)) + phases_real, phases_imag = real.(phases), imag.(phases) + phases = A(phases_real) + im*A(phases_imag) + psih = @. phases*psik if allones; psih = psik; end - psih = psih.*mask + psih = psih.*A(mask) Ein = real(sum(g.Krsq.*abs2.(psih)/(g.nx*g.ny)^2)) psih = psih*sqrt(E0/Ein) - q = -irfft(g.Krsq.*psih, g.nx) + q = A(-irfft(g.Krsq.*psih, g.nx)) end end diff --git a/test/runtests.jl b/test/runtests.jl index e08a817b..8754aa85 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,11 +1,9 @@ -#!/usr/bin/env julia - using FourierFlows, - Test, Statistics, Random, - FFTW + FFTW, + Test import # use 'import' rather than 'using' for submodules to keep namespace clean GeophysicalFlows.TwoDTurb, @@ -16,6 +14,11 @@ import # use 'import' rather than 'using' for submodules to keep namespace clean using FourierFlows: parsevalsum, xmoment, ymoment using GeophysicalFlows: lambdipole, peakedisotropicspectrum +# the devices on which tests will run +devices = (CPU(),) +@has_cuda devices = (CPU(), GPU()) +@has_cuda using CuArrays + const rtol_lambdipole = 1e-2 # tolerance for lamb dipole tests const rtol_multilayerqg = 1e-13 # tolerance for multilayerqg forcing tests const rtol_twodturb = 1e-13 # tolerance for twodturb forcing tests @@ -28,24 +31,33 @@ cfl(prob) = cfl(prob.vars.u, prob.vars.v, prob.clock.dt, prob.grid.dx) # Run tests testtime = @elapsed begin -@testset "Utils" begin - include("test_utils.jl") +for dev in devices + + println("testing on "*string(typeof(dev))) + + @testset "Utils" begin + include("test_utils.jl") + + @test testpeakedisotropicspectrum(dev) + @test_throws ErrorException("the domain is not square") testpeakedisotropicspectrum_rectangledomain() + end + + @testset "TwoDTurb" begin + include("test_twodturb.jl") + + @test test_twodturb_advection(0.0005, "ForwardEuler", dev) + @test test_twodturb_lambdipole(256, 1e-3, dev) + @test test_twodturb_stochasticforcingbudgets(dev) + @test test_twodturb_deterministicforcingbudgets(dev) + @test test_twodturb_energyenstrophy(dev) + @test test_twodturb_problemtype(Float32) + @test TwoDTurb.nothingfunction() == nothing + end - @test testpeakedisotropicspectrum() - @test_throws ErrorException("the domain is not square") testpeakedisotropicspectrum_rectangledomain() end -@testset "TwoDTurb" begin - include("test_twodturb.jl") - @test test_twodturb_advection(0.0005, "ForwardEuler") - @test test_twodturb_lambdipole(256, 1e-3) - @test test_twodturb_stochasticforcingbudgets() - @test test_twodturb_deterministicforcingbudgets() - @test test_twodturb_energyenstrophy() - @test test_twodturb_problemtype(Float32) - @test TwoDTurb.nothingfunction() == nothing -end +println("rest of tests only on CPU") @testset "BarotropicQG" begin include("test_barotropicqg.jl") diff --git a/test/test_twodturb.jl b/test/test_twodturb.jl index b7e963f2..5be8e293 100644 --- a/test/test_twodturb.jl +++ b/test/test_twodturb.jl @@ -1,11 +1,11 @@ -function test_twodturb_lambdipole(n, dt; L=2π, Ue=1, Re=L/20, nu=0.0, nnu=1, ti=L/Ue*0.01, nm=3) +function test_twodturb_lambdipole(n, dt, dev::Device=CPU(); L=2π, Ue=1, Re=L/20, ν=0.0, nν=1, ti=L/Ue*0.01, nm=3) nt = round(Int, ti/dt) - prob = TwoDTurb.Problem(nx=n, Lx=L, nu=nu, nnu=nnu, dt=dt, stepper="FilteredRK4") + prob = TwoDTurb.Problem(nx=n, Lx=L, ν=ν, nν=nν, dt=dt, stepper="FilteredRK4", dev=dev) zeta₀ = lambdipole(Ue, Re, prob.grid) TwoDTurb.set_zeta!(prob, zeta₀) - xzeta = zeros(nm) # centroid of abs(zeta) - Ue_m = zeros(nm) # measured dipole speed + xzeta = zeros(nm) # centroid of abs(zeta) + Ue_m = zeros(nm) # measured dipole speed x, y = gridpoints(prob.grid) zeta = prob.vars.zeta @@ -20,19 +20,19 @@ function test_twodturb_lambdipole(n, dt; L=2π, Ue=1, Re=L/20, nu=0.0, nnu=1, ti isapprox(Ue, mean(Ue_m[2:end]), rtol=rtol_lambdipole) end -function test_twodturb_stochasticforcingbudgets(; n=256, L=2π, dt=0.005, nu=1e-7, nnu=2, mu=1e-1, nmu=0, tf=0.1/mu) +function test_twodturb_stochasticforcingbudgets(dev::Device=CPU(); n=256, L=2π, dt=0.005, ν=1e-7, nν=2, μ=1e-1, nμ=0, tf=0.1/μ) nt = round(Int, tf/dt) - # Forcing + # Forcing parameters kf, dkf = 12.0, 2.0 ε = 0.1 - gr = TwoDGrid(n, L) + gr = TwoDGrid(dev, n, L) x, y = gridpoints(gr) - Kr = [ gr.kr[i] for i=1:gr.nkr, j=1:gr.nl] + Kr = ArrayType(dev)([ gr.kr[i] for i=1:gr.nkr, j=1:gr.nl]) - force2k = zero(gr.Krsq) + force2k = ArrayType(dev)(zero(gr.Krsq)) @. force2k = exp(-(sqrt(gr.Krsq)-kf)^2/(2*dkf^2)) @. force2k[gr.Krsq .< 2.0^2 ] = 0 @. force2k[gr.Krsq .> 20.0^2 ] = 0 @@ -43,14 +43,14 @@ function test_twodturb_stochasticforcingbudgets(; n=256, L=2π, dt=0.005, nu=1e- Random.seed!(1234) function calcF!(Fh, sol, t, cl, v, p, g) - eta = exp.(2π*im*rand(Float64, size(sol)))/sqrt(cl.dt) + eta = ArrayType(dev)(exp.(2π*im*rand(Float64, size(sol)))/sqrt(cl.dt)) eta[1, 1] = 0.0 @. Fh = eta * sqrt(force2k) nothing end - prob = TwoDTurb.Problem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, nmu=nmu, dt=dt, - stepper="RK4", calcF=calcF!, stochastic=true) + prob = TwoDTurb.Problem(nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt, + stepper="RK4", calcF=calcF!, stochastic=true, dev=dev) sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid; @@ -65,7 +65,7 @@ function test_twodturb_stochasticforcingbudgets(; n=256, L=2π, dt=0.005, nu=1e- TwoDTurb.updatevars!(prob) E, D, W, R = diags - t = round(mu*cl.t, digits=2) + t = round(μ*cl.t, digits=2) i₀ = 1 dEdt = (E[(i₀+1):E.i] - E[i₀:E.i-1])/cl.dt @@ -83,28 +83,26 @@ function test_twodturb_stochasticforcingbudgets(; n=256, L=2π, dt=0.005, nu=1e- end -function test_twodturb_deterministicforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7, nnu=2, mu=1e-1, nmu=0) +function test_twodturb_deterministicforcingbudgets(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, μ=1e-1, nμ=0) n, L = 256, 2π - nu, nnu = 1e-7, 2 - mu, nmu = 1e-1, 0 - dt, tf = 0.005, 0.1/mu + ν, nν = 1e-7, 2 + μ, nμ = 1e-1, 0 + dt, tf = 0.005, 0.1/μ nt = round(Int, tf/dt) - g = TwoDGrid(n, L) - - gr = TwoDGrid(n, L) + + gr = TwoDGrid(dev, n, L) x, y = gridpoints(gr) # Forcing = 0.01cos(4x)cos(5y)cos(2t) - - f = @. 0.01*cos(4*x)*cos(5*y) + f = @. 0.01cos(4x)*cos(5y) fh = rfft(f) - function calcF!(Fh, sol, t, cl, v, p, g) - @. Fh = fh*cos(2*t) + function calcF!(Fh, sol, t, cl, v, p, g::AbstractGrid{T, A}) where {T, A} + Fh = fh*cos(2t) nothing end - prob = TwoDTurb.Problem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, nmu=nmu, dt=dt, - stepper="RK4", calcF=calcF!, stochastic=false) + prob = TwoDTurb.Problem(nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt, + stepper="RK4", calcF=calcF!, stochastic=false, dev=dev) sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid @@ -121,7 +119,7 @@ function test_twodturb_deterministicforcingbudgets(; n=256, dt=0.01, L=2π, nu=1 TwoDTurb.updatevars!(prob) E, D, W, R = diags - t = round(mu*cl.t, digits=2) + t = round(μ*cl.t, digits=2) i₀ = 1 dEdt = (E[(i₀+1):E.i] - E[i₀:E.i-1])/cl.dt @@ -142,25 +140,25 @@ Tests the advection term in the twodturb module by timestepping a test problem with timestep dt and timestepper identified by the string stepper. The test problem is derived by picking a solution ζf (with associated streamfunction ψf) for which the advection term J(ψf, ζf) is non-zero. Next, a -forcing Ff is derived according to Ff = ∂ζf/∂t + J(ψf, ζf) - nuΔζf. One solution +forcing Ff is derived according to Ff = ∂ζf/∂t + J(ψf, ζf) - νΔζf. One solution to the vorticity equation forced by this Ff is then ζf. (This solution may not be realized, at least at long times, if it is unstable.) """ -function test_twodturb_advection(dt, stepper; n=128, L=2π, nu=1e-2, nnu=1, mu=0.0, nmu=0) +function test_twodturb_advection(dt, stepper, dev::Device=CPU(); n=128, L=2π, ν=1e-2, nν=1, μ=0.0, nμ=0) n, L = 128, 2π - nu, nnu = 1e-2, 1 - mu, nmu = 0.0, 0 + ν, nν = 1e-2, 1 + μ, nμ = 0.0, 0 tf = 1.0 nt = round(Int, tf/dt) - gr = TwoDGrid(n, L) + gr = TwoDGrid(dev, n, L) x, y = gridpoints(gr) psif = @. sin(2x)*cos(2y) + 2sin(x)*cos(3y) zetaf = @. -8sin(2x)*cos(2y) - 20sin(x)*cos(3y) Ff = @. -( - nu*( 64sin(2x)*cos(2y) + 200sin(x)*cos(3y) ) + ν*( 64sin(2x)*cos(2y) + 200sin(x)*cos(3y) ) + 8*( cos(x)*cos(3y)*sin(2x)*sin(2y) - 3cos(2x)*cos(2y)*sin(x)*sin(3y) ) ) @@ -172,7 +170,7 @@ function test_twodturb_advection(dt, stepper; n=128, L=2π, nu=1e-2, nnu=1, mu=0 nothing end - prob = TwoDTurb.Problem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, nmu=nmu, dt=dt, stepper=stepper, calcF=calcF!, stochastic=false) + prob = TwoDTurb.Problem(nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt, stepper=stepper, calcF=calcF!, stochastic=false, dev=dev) sol, cl, p, v, g = prob.sol, prob.clock, prob.params, prob.vars, prob.grid TwoDTurb.set_zeta!(prob, zetaf) @@ -182,10 +180,10 @@ function test_twodturb_advection(dt, stepper; n=128, L=2π, nu=1e-2, nnu=1, mu=0 isapprox(prob.vars.zeta, zetaf, rtol=rtol_twodturb) end -function test_twodturb_energyenstrophy() +function test_twodturb_energyenstrophy(dev::Device=CPU()) nx, Lx = 128, 2π ny, Ly = 128, 3π - gr = TwoDGrid(nx, Lx, ny, Ly) + gr = TwoDGrid(dev, nx, Lx, ny, Ly) x, y = gridpoints(gr) k0, l0 = gr.k[2], gr.l[2] # fundamental wavenumbers @@ -195,7 +193,7 @@ function test_twodturb_energyenstrophy() energy_calc = 29/9 enstrophy_calc = 2701/162 - prob = TwoDTurb.Problem(nx=nx, Lx=Lx, ny=ny, Ly=Ly, stepper="ForwardEuler") + prob = TwoDTurb.Problem(nx=nx, Lx=Lx, ny=ny, Ly=Ly, stepper="ForwardEuler", dev=dev) sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid; @@ -205,15 +203,15 @@ function test_twodturb_energyenstrophy() energyzeta0 = TwoDTurb.energy(prob) enstrophyzeta0 = TwoDTurb.enstrophy(prob) - params = TwoDTurb.Params(p.nu, p.nnu) + params = TwoDTurb.Params(p.ν, p.nν) (isapprox(energyzeta0, energy_calc, rtol=rtol_twodturb) && isapprox(enstrophyzeta0, enstrophy_calc, rtol=rtol_twodturb) && TwoDTurb.addforcing!(prob.timestepper.N, sol, cl.t, cl, v, p, g)==nothing && p == params) end -function test_twodturb_problemtype(T=Float32) - prob = TwoDTurb.Problem(T=T) +function test_twodturb_problemtype(T=Float32, dev::Device=CPU()) + prob = TwoDTurb.Problem(T=T, dev=dev) (typeof(prob.sol)==Array{Complex{T},2} && typeof(prob.grid.Lx)==T && typeof(prob.grid.x)==Array{T,2} && typeof(prob.vars.u)==Array{T,2}) end diff --git a/test/test_utils.jl b/test/test_utils.jl index a4c366c5..db1bc898 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -1,9 +1,9 @@ """ Test the peakedisotropicspectrum function. """ -function testpeakedisotropicspectrum() +function testpeakedisotropicspectrum(dev::Device=CPU()) n, L = 128, 2π - gr = TwoDGrid(n, L) + gr = TwoDGrid(dev, n, L) k0, E0 = 6, 0.5 qi = peakedisotropicspectrum(gr, k0, E0; allones=true) ρ, qhρ = FourierFlows.radialspectrum(rfft(qi).*gr.invKrsq, gr)