From 6f0ed7415ddd6d0c0ab81aa3d24420a349a70a4e Mon Sep 17 00:00:00 2001 From: Navid Constantinou Date: Mon, 22 Oct 2018 17:21:45 +1100 Subject: [PATCH 1/2] simplifies addforcing function --- src/barotropicqg.jl | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/src/barotropicqg.jl b/src/barotropicqg.jl index a67fc82b..d28cd764 100644 --- a/src/barotropicqg.jl +++ b/src/barotropicqg.jl @@ -210,7 +210,7 @@ function calcN_advection!(N, sol, t, s, v, p, g) ldiv!(v.zeta, g.rfftplan, v.zetah) ldiv!(v.u, g.rfftplan, v.uh) - v.psih .= deepcopy(v.vh) # FFTW's irfft destroys its input; v.vh is needed for N[1, 1] + v.psih .= v.vh # FFTW's irfft destroys its input; v.vh is needed for N[1, 1] ldiv!(v.v, g.rfftplan, v.psih) @. v.q = v.zeta + p.eta @@ -227,6 +227,12 @@ end function calcN!(N, sol, t, s, v, p, g) calcN_advection!(N, sol, t, s, v, p, g) addforcing!(N, t, s, v, p, g) + if p.calcFU != nothingfunction + # 'Nonlinear' term for U with topographic correlation. + # Note: < v*eta > = sum( conj(vh)*eta ) / (nx^2*ny^2) if fft is used + # while < v*eta > = 2*sum( conj(vh)*eta ) / (nx^2*ny^2) if rfft is used + N[1, 1] = p.calcFU(t) + 2*sum(conj(v.vh).*p.etah).re / (g.nx^2.0*g.ny^2.0) + end nothing end @@ -235,12 +241,6 @@ addforcing!(N, t, s, v::Vars, p, g) = nothing function addforcing!(N, t, s, v::ForcedVars, p, g) p.calcFq!(v.Fqh, t, s, v, p, g) @. N += v.Fqh - if p.calcFU != nothingfunction - # 'Nonlinear' term for U with topographic correlation. - # Note: < v*eta > = sum( conj(vh)*eta ) / (nx^2*ny^2) if fft is used - # while < v*eta > = 2*sum( conj(vh)*eta ) / (nx^2*ny^2) if rfft is used - N[1, 1] = p.calcFU(t) + 2*sum(conj(v.vh).*p.etah).re / (g.nx^2.0*g.ny^2.0) - end nothing end @@ -250,13 +250,6 @@ function addforcing!(N, t, s, v::StochasticForcedVars, p, g) p.calcFq!(v.Fqh, t, s, v, p, g) end @. N += v.Fqh - - if p.calcFU != nothingfunction - # 'Nonlinear' term for U with topographic correlation. - # Note: < v*eta > = sum( conj(vh)*eta ) / (nx^2*ny^2) if fft is used - # while < v*eta > = 2*sum( conj(vh)*eta ) / (nx^2*ny^2) if rfft is used - N[1, 1] = p.calcFU(t) + 2*sum(conj(v.vh).*p.etah).re / (g.nx^2.0*g.ny^2.0) - end nothing end From f1d90eae8b0f01307636dda29ed475621b6ebd2f Mon Sep 17 00:00:00 2001 From: Navid Constantinou Date: Mon, 22 Oct 2018 17:21:45 +1100 Subject: [PATCH 2/2] adds zonal mean psi in Vars --- src/barotropicqgql.jl | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/barotropicqgql.jl b/src/barotropicqgql.jl index 58571f06..4b632d6a 100644 --- a/src/barotropicqgql.jl +++ b/src/barotropicqgql.jl @@ -27,7 +27,7 @@ import FFTW: rfft abstract type BarotropicQGQLVars <: AbstractVars end abstract type BarotropicQGQLForcedVars <: BarotropicQGQLVars end -const physicalvars = [:zeta, :psi, :u, :v, :uzeta, :vzeta, :U, :Zeta] +const physicalvars = [:zeta, :psi, :u, :v, :uzeta, :vzeta, :U, :Zeta, :Psi] const transformvars = [ Symbol(var, :h) for var in physicalvars ] const forcedvars = [:Fh] const stochforcedvars = [:prevsol] @@ -150,6 +150,7 @@ mutable struct Vars{T} <: BarotropicQGQLVars zeta::Array{T,2} Zeta::Array{T,2} psi::Array{T,2} + Psi::Array{T,2} Nz::Array{Complex{T},2} NZ::Array{Complex{T},2} uh::Array{Complex{T},2} @@ -158,12 +159,13 @@ mutable struct Vars{T} <: BarotropicQGQLVars zetah::Array{Complex{T},2} Zetah::Array{Complex{T},2} psih::Array{Complex{T},2} + Psih::Array{Complex{T},2} end function Vars(g; T=typeof(g.Lx)) - @createarrays T (g.nx, g.ny) u v U uzeta vzeta zeta Zeta psi - @createarrays Complex{T} (g.nkr, g.nl) N NZ uh vh Uh zetah Zetah psih - Vars(u, v, U, uzeta, vzeta, zeta, Zeta, psi, N, NZ, uh, vh, Uh, zetah, Zetah, psih) + @createarrays T (g.nx, g.ny) u v U uzeta vzeta zeta Zeta psi Psi + @createarrays Complex{T} (g.nkr, g.nl) N NZ uh vh Uh zetah Zetah psih Psih + Vars(u, v, U, uzeta, vzeta, zeta, Zeta, psi, Psi, N, NZ, uh, vh, Uh, zetah, Zetah, psih, Psih) end """ @@ -180,6 +182,7 @@ mutable struct ForcedVars{T} <: BarotropicQGQLForcedVars zeta::Array{T,2} Zeta::Array{T,2} psi::Array{T,2} + Psi::Array{T,2} Nz::Array{Complex{T},2} NZ::Array{Complex{T},2} uh::Array{Complex{T},2} @@ -188,6 +191,7 @@ mutable struct ForcedVars{T} <: BarotropicQGQLForcedVars zetah::Array{Complex{T},2} Zetah::Array{Complex{T},2} psih::Array{Complex{T},2} + Psih::Array{Complex{T},2} Fh::Array{Complex{T},2} end @@ -207,6 +211,7 @@ mutable struct StochasticForcedVars{T} <: BarotropicQGQLForcedVars zeta::Array{T,2} Zeta::Array{T,2} psi::Array{T,2} + Psi::Array{T,2} Nz::Array{Complex{T},2} NZ::Array{Complex{T},2} uh::Array{Complex{T},2} @@ -215,6 +220,7 @@ mutable struct StochasticForcedVars{T} <: BarotropicQGQLForcedVars zetah::Array{Complex{T},2} Zetah::Array{Complex{T},2} psih::Array{Complex{T},2} + Psih::Array{Complex{T},2} Fh::Array{Complex{T},2} prevsol::Array{Complex{T},2} end @@ -310,6 +316,7 @@ function updatevars!(s, v, p, g) @. v.Zetah = s.sol @. v.Zetah[abs.(g.Kr).>0] = 0 + @. v.Psih = -v.Zetah * g.invKKrsq @. v.psih = -v.zetah * g.invKKrsq @. v.uh = -im * g.l * v.psih @. v.vh = im * g.kr * v.psih @@ -318,6 +325,7 @@ function updatevars!(s, v, p, g) ldiv!(v.zeta, g.rfftplan, deepcopy(v.zetah)) ldiv!(v.Zeta, g.rfftplan, deepcopy(v.Zetah)) ldiv!(v.psi, g.rfftplan, v.psih) + ldiv!(v.Psi, g.rfftplan, v.Psih) ldiv!(v.u, g.rfftplan, deepcopy(v.uh)) ldiv!(v.v, g.rfftplan, deepcopy(v.vh)) ldiv!(v.U, g.rfftplan, deepcopy(v.Uh))