Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

v0.1.1 Correction #3

Merged
merged 3 commits into from
Sep 7, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "MHDFlows"
uuid = "1d939cba-ab73-4bc0-975c-87d4c856e1f9"
authors = ["Ka Wai HO <woodyho13@gmail.com> "]
version = "0.1.0"
version = "0.1.1"

[deps]
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
Expand Down
139 changes: 139 additions & 0 deletions src/utils/IC.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,139 @@
# ----------
# General Function Module, providing function for setup IC of the problem
# ----------
include("utils.jl")

"""
Construct a Cylindrical Mask Function χ for VP method
Keyword arguments
=================
- `grid`: MHDFlows problem's grid
- `R₂` : Outwards radius boundary
- `R₁` : Inwards radius boundary
$(TYPEDFIELDS)
"""
function Cylindrical_Mask_Function(grid;R₂=0.82π,R₁=0.0π)
nx,ny,nz = grid.nx,grid.ny,grid.nz;
x,y,z = grid.x,grid.y,grid.z;
S = BitArray(undef, nx::Int,ny::Int,nz::Int);

for k ∈ 1:nz::Int, j ∈ 1:ny::Int,i ∈ 1:nx::Int
xᵢ,yᵢ,zᵢ = x[i]::AbstractFloat,y[j]::AbstractFloat,z[k]::AbstractFloat;
Rᵢ = √(xᵢ^2+yᵢ^2);
# S = 0 if inside fluid domain while S = 1 in the solid domain
S[i,j,k] = (R₂ >= Rᵢ >= R₁) ? 0 : 1;
end
return S
end
"""
function of setting up the initial condition of the problem
Keyword arguments
=================
- `prob`: MHDFlows problem
- `ux/uy/uz` : velocity in real space
- `bx/by/bz` : B-field in real space
- `U₀x/U₀y/U₀z/B₀x/B₀y/B₀z` : VP method parameter
$(TYPEDFIELDS)
"""
function SetUpProblemIC!(prob; ux = [], uy = [], uz =[],
bx = [], by = [], bz =[],
U₀x= [], U₀y= [], U₀z=[],
B₀x= [], B₀y= [], B₀z=[])
sol = prob.sol;
vars = prob.vars;
grid = prob.grid;
params = prob.params;
# Copy the data to both output and solution array
for (uᵢ,prob_uᵢ,uᵢind) in zip([ux,uy,uz],[vars.ux,vars.uy,vars.uz],
[params.ux_ind,params.uy_ind,params.uz_ind])
if uᵢ != []
@views sol₀ = sol[:, :, :, uᵢind];
copyto!(prob_uᵢ,uᵢ);
mul!(sol₀ , grid.rfftplan, prob_uᵢ);
end
end
if prob.flag.b
for (bᵢ,prob_bᵢ,bᵢind) in zip([bx,by,bz],[vars.bx,vars.by,vars.bz],
[params.bx_ind,params.by_ind,params.bz_ind])
if bᵢ != []
@views sol₀ = sol[:, :, :, bᵢind];
copyto!(prob_bᵢ,bᵢ);
mul!(sol₀ , grid.rfftplan, prob_bᵢ);
end
end
end
if prob.flag.vp
for (Uᵢ,prob_Uᵢ) in zip([U₀x,U₀y,U₀z],
[params.U₀x,params.U₀y,params.U₀z])
Uᵢ == [] ? nothing : copyto!(prob_Uᵢ,Uᵢ);
end

if prob.flag.b
for (Bᵢ,prob_Bᵢ) in zip([B₀x,B₀y,B₀z],
[params.B₀x,params.B₀y,params.B₀z])
Bᵢ == [] ? nothing : copyto!(prob_Bᵢ,Bᵢ);
end
end
end
return nothing;

end


"""
Construct a Div Free Spectra Vector Map
Keyword arguments
=================
- `Nx/Ny/Nz`: size of the Vector Map
- `k0` : Slope of the Map
- `b` : Anisotropy of the Map
$(TYPEDFIELDS)
"""
function DivFreeSpectraMap( Nx::Int, Ny::Int, Nz::Int;
P = 1, k0 = -5/3/2, b = 1, T = Float64)
grid = GetSimpleThreeDGrid(Nx,Ny,Nz, T = T)
return DivFreeSpectraMap( grid; P = P, k0 = k0, b = b);
end

function DivFreeSpectraMap( grid;
P = 1, k0 = -5/3/2, b = 1)

T = eltype(grid);
@devzeros typeof(CPU()) Complex{T} (grid.nkr,grid.nl,grid.nm) eⁱᶿ Fk Fxh Fyh Fzh
@devzeros typeof(CPU()) T (grid.nx ,grid.ny,grid.nz) Fx Fy Fz

kx,ky,kz = grid.kr,grid.l,grid.m;
Lx,Ly,Lz = grid.Lx,grid.Ly,grid.Lz;
dx,dy,dz = grid.dx,grid.dy,grid.dz;
k⁻¹ = @. √(grid.invKrsq);
k = @. √(grid.Krsq);
k⊥ = @. √(kx^2 + ky^2);
dk⁻² = @. 1/(k+1)^2;
Fk = @. k.^(k0);
Fk[1,1,1] = 0.0;
# Fk[k.>5] .= 0;
∫Fkdk = sum(@. Fk*dk⁻²);
A = sqrt(P*3*(Lx/dx)*(Ly/dy)*(Lz/dz)/∫Fkdk*(1/dx/dy/dz));
Fk*=A;

e1x = @. ky/k⊥;
e1y = @. -kx/k⊥;
e2x = @. kx*kz/k⊥*k⁻¹;
e2y = @. ky*kz/k⊥*k⁻¹;
e2z = @. -k⊥*k⁻¹;
e1x[isnan.(e1x)] .= 0;
e1y[isnan.(e1y)] .= 0;
e2x[isnan.(e2x)] .= 0;
e2y[isnan.(e2y)] .= 0;

# Work out the random conponement
eⁱᶿ .= exp.(im.*rand(T,grid.nkr,grid.nl,grid.nm)*2π);
@. Fxh += Fk*eⁱᶿ*e2x;
@. Fyh += Fk*eⁱᶿ*e2y;
@. Fzh += Fk*eⁱᶿ*e2z;
ldiv!(Fx, grid.rfftplan, deepcopy(Fxh));
ldiv!(Fy, grid.rfftplan, deepcopy(Fyh));
ldiv!(Fz, grid.rfftplan, deepcopy(Fzh));
return Fx,Fy,Fz;

end
12 changes: 6 additions & 6 deletions src/utils/MHDAnalysis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,10 @@
# ----------
include("utils.jl")


# Scale Decomposition FUnction
function ScaleDecomposition(B1::Array;kf=[1,5],Lx = 2π,T=Float32)
nx,ny,nz = size(B1);
grid = ThreeDGrid(nx, Lx, T = T);
grid = GetSimpleThreeDGrid(nx, Lx, T = T);
cB1 = ScaleDecomposition(B1,grid;kf=kf)
return cB1;
end
Expand Down Expand Up @@ -38,7 +37,7 @@ end

function ScaleDecomposition(B1::Array,B2::Array,B3::Array;kf=[1,5],Lx = 2π,T=Float32)
nx,ny,nz = size(B1);
grid = ThreeDGrid(nx, Lx, T = T);
grid = GetSimpleThreeDGrid(nx, Lx, T = T);

cB1,cB2,cB3 = ScaleDecomposition(B1,B2,B3,grid;kf=kf)
return cB1,cB2,cB3;
Expand Down Expand Up @@ -94,7 +93,7 @@ end
function VectorPotential(B1::Array{T,3},B2::Array{T,3},B3::Array{T,3};L=2π) where T
# Wrapper of actual Vector Potential function
nx,ny,nz = size(B1);
grid = ThreeDGrid(nx, L, T = T);
grid = GetSimpleGetSimpleThreeDGrid(nx, L; T = T);
A1,A2,A3 = VectorPotential(B1,B2,B3,grid);
return A1,A2,A3;
end
Expand Down Expand Up @@ -140,7 +139,7 @@ end
#Checking for anagular momentum using center point as a reference (r = 0 at the center)
function getL(iv::Array{T,3},jv::Array{T,3},kv::Array{T,3};L=2π) where T
nx,ny,nz = size(iv);
grid = ThreeDGrid(nx, L, T = T);
grid = GetSimpleThreeDGrid(nx, L, T = T);
Lᵢ,Lⱼ,Lₖ = getL(iv,jv,kv,grid)
return Lᵢ,Lⱼ,Lₖ
end
Expand Down Expand Up @@ -181,7 +180,8 @@ end
function spectralline(A::Array{T,3};Lx=2π) where T
nx,ny,nz = size(A);
Ak = zeros(Complex{T},div(nx,2)+1,ny,nz);
k²,rfftplan = Getk²_and_fftplan(nx,Lx;T=T);
grid = GetSimpleGetSimpleThreeDGrid(nx,Lx;T=T);
k²,rfftplan = grid.Krsq,grid.rfftplan;
mul!(Ak,rfftplan,A);
kk = @. √(k²::Array{T,3});
krmax = round(Int,maximum(kk)+1);
Expand Down
8 changes: 4 additions & 4 deletions src/utils/VectorCalculus.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ function Curl(B1::Array,B2::Array,B3::Array;
Lx = 2π, Ly = Lx, Lz = Lx,T = Float32)
# Wrapper for Curl Function
nx,ny,nz = size(B1);
grid = ThreeDGrid(nx, Lx, ny, Ly, nz, Lz, T = T);
grid = GetSimpleThreeDGrid(nx, Lx, ny, Ly, nz, Lz, T = T);
cB1,cB2,cB3 = Curl(B1,B2,B3,grid)
return cB1,cB2,cB3;
end
Expand Down Expand Up @@ -43,7 +43,7 @@ end
function Div(B1::Array,B2::Array,B3::Array;
Lx = 2π, Ly = Lx, Lz = Lx,T = Float32)
nx,ny,nz = size(B1);
grid = ThreeDGrid(nx, Lx, ny, Ly, nz, Lz, T = T);
grid = GetSimpleGrid(nx, Lx, ny, Ly, nz, Lz, T = T);
cB1 = Div(B1,B2,B3,grid);

return cB1
Expand Down Expand Up @@ -76,7 +76,7 @@ end
function ∂i(B1::Array, direction; Lx = 2π, Ly = Lx, Lz = Lx,T = eltype(B1))
nx,ny,nz = size(B1);
dev = typeof(B1) <: Array ? CPU() : GPU();
grid = ThreeDGrid(dev, nx, Lx, ny, Ly, nz, Lz, T = T);
grid = GetSimpleThreeDGrid(nx, Lx, ny, Ly, nz, Lz, T = T);
return ∂i(B1,grid, direction);
end

Expand Down Expand Up @@ -114,7 +114,7 @@ end

function LaplaceSolver(B; Lx=2π, Ly = Lx, Lz = Lx, T = Float32)
nx,ny,nz = size(B);
grid = ThreeDGrid(nx, Lx, ny, Ly, nz, Lz, T = T);
grid = GetSimpleThreeDGrid(nx, Lx, ny, Ly, nz, Lz, T = T);
Φ = LaplaceSolver(B,grid);
return Φ
end
Expand Down
59 changes: 39 additions & 20 deletions src/utils/utils.jl
Original file line number Diff line number Diff line change
@@ -1,20 +1,39 @@
function Getk²_and_fftplan(nx = 64, Lx = 2π, ny = nx, Ly = Lx, nz = nx, Lz = Lx;
nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE,
T=Float64, ArrayType=Array)

nk = nx
nl = ny
nm = nz
nkr = Int(nx/2 + 1)
# Wavenubmer
k = ArrayType{T}(reshape( fftfreq(nx, 2π/Lx*nx), (nk, 1, 1)))
l = ArrayType{T}(reshape( fftfreq(ny, 2π/Ly*ny), (1, nl, 1)))
m = ArrayType{T}(reshape( fftfreq(nz, 2π/Lz*nz), (1, 1, nm)))
kr = ArrayType{T}(reshape(rfftfreq(nx, 2π/Lx*nx), (nkr, 1, 1)))
k² = @. k^2 + l^2 + m^2;

FFTW.set_num_threads(nthreads);
rfftplan = plan_rfft(ArrayType{T, 3}(undef, nx, ny, nz))

return k², rfftplan;
end
mutable struct SimpleGrid{i,i²,plan}
k :: i
l :: i
m :: i
kr :: i
Ksq :: i²
invKsq :: i²
Krsq :: i²
invKrsq :: i²
rfftplan :: plan
end
Base.eltype(grid::SimpleGrid) = eltype(grid.k);

function GetSimpleThreeDGrid( nx = 64, Lx = 2π, ny = nx, Ly = Lx, nz = nx, Lz = Lx;
nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE,
T=Float64, ArrayType=Array)
nk = nx
nl = ny
nm = nz
nkr = Int(nx/2 + 1)
# Wavenubmer
k = ArrayType{T}(reshape( fftfreq(nx, 2π/Lx*nx), (nk, 1, 1)))
l = ArrayType{T}(reshape( fftfreq(ny, 2π/Ly*ny), (1, nl, 1)))
m = ArrayType{T}(reshape( fftfreq(nz, 2π/Lz*nz), (1, 1, nm)))
kr = ArrayType{T}(reshape(rfftfreq(nx, 2π/Lx*nx), (nkr, 1, 1)))

Ksq = @. k^2 + l^2 + m^2
invKsq = @. 1 / Ksq
invKsq[1, 1, 1] = 0;

Krsq = @. kr^2 + l^2 + m^2
invKrsq = @. 1 / Krsq
invKrsq[1, 1, 1] = 0;

FFTW.set_num_threads(nthreads);
rfftplan = plan_rfft(ArrayType{T, 3}(undef, nx, ny, nz))

return SimpleGrid(k,l,m,kr,Ksq, invKsq, Krsq, invKrsq, rfftplan);
end