v 0.1.3 Update

Manifest.toml

@@ -117,7 +117,7 @@ uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 deps = ["CUDA", "DocStringExtensions", "FFTW", "Interpolations", "JLD2", "LinearAlgebra", "Random", "Reexport", "Statistics"]
 git-tree-sha1 = "4d54809841acd3a5c20c135f94b92f23e75f0307"
 uuid = "2aec4490-903f-5c70-9b11-9bed06a700e1"
-version = "0.9.2"
+version = "0.10.1"
 
 [[deps.GPUArrays]]
 deps = ["Adapt", "GPUArraysCore", "LLVM", "LinearAlgebra", "Printf", "Random", "Reexport", "Serialization", "Statistics"]

Project.toml

@@ -1,25 +1,25 @@
-name = "MHDFlows"
-uuid = "1d939cba-ab73-4bc0-975c-87d4c856e1f9"
-authors = ["Ka Wai HO <woodyho13@gmail.com> "]
-version = "0.1.2"
-
-[deps]
-CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
-DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-FourierFlows = "2aec4490-903f-5c70-9b11-9bed06a700e1"
-FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
-Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-
-[compat]
-CUDA = "^1, ^2.4.2, 3.0.0 - 3.6.4, ^3.7.1"
-DocStringExtensions = "^0.8"
-FFTW = "^1"
-FourierFlows = "^0.9"
-Reexport = "^0.2, ^1"
-SpecialFunctions = "^0.10, ^1, 2"
-julia = "^1.5.3"
-HDF5 = "^0.14.3"
+name = "MHDFlows"
+uuid = "1d939cba-ab73-4bc0-975c-87d4c856e1f9"
+authors = ["Ka Wai HO <woodyho13@gmail.com> "]
+version = "0.1.3"
+
+[deps]
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+FourierFlows = "2aec4490-903f-5c70-9b11-9bed06a700e1"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+
+[compat]
+CUDA = "^1, ^2.4.2, 3.0.0 - 3.6.4, ^3.7.1"
+DocStringExtensions = "^0.8"
+FFTW = "^1"
+FourierFlows = "^0.10.1"
+Reexport = "^0.2, ^1"
+julia = "^1.5.3"
+HDF5 = "^0.14.3"
+ProgressMeter = "^1.7.2"

README.md

@@ -1,5 +1,4 @@
 # MHDFlows.jl
-
 ![Julia flow](img/TG_Instability.gif)
 
 Three Dimensional Magnetohydrodynamic(MHD) pseudospectral solvers written in Julia language with <a href="http://github.com/FourierFlows/FourierFlows.jl">FourierFlows.jl</a>. This solver support the following features:
@@ -12,12 +11,10 @@ Three Dimensional Magnetohydrodynamic(MHD) pseudospectral solvers written in Jul
 This package leverages the [FourierFlows.jl](http://github.com/FourierFlows/FourierFlows.jl) package to set up the module. The main purpose of MHDFlows.jl aims to solve the portable 3D MHD problems on personal computer instead of cluster. Utilizing the Nvidia CUDA technology, the MHDFlows.jl could solve the front-end MHD turbulence research problems in the order of few-ten minutes by using a mid to high end gaming display card (see Memory usage & speed section). Feel free to modify yourself for your own research purpose.
 
 ## Version No.
-
-v 0.1.0
+v 0.1.3
 
 ## Installation Guide & compatibility 
-
-The current version is tested on v1.5.3/1.7.3/1.8 version.
+The current version is tested on v1.5.3/1.7.3/1.8.2 version.
 
 Currently, you have two way of installing MHDFlows.jl
 
@@ -31,9 +28,7 @@ Currently, you have two way of installing MHDFlows.jl
    ```
 
 
-
 ## Scalability 
-
 The MHD Solver could either run on CPU or GPU. The scalability is same as Fourierflows, which is restricted to either a single CPU or single GPU. This restriction may change in the future depending on the development of FourierFlows. If you are running the package using GPU, the CUDA package is needed for the computation. Check out [CUDA.jl](https://juliagpu.github.io/CUDA.jl/stable/lib/driver/#Device-Management) for more detail. 
 
 **Important note**: As this is a 3D solver, the memory usage and computational time scale as `N^3` (at least). Beware of the memory usage especially when you are using the GPU. 
@@ -44,9 +39,9 @@ The MHD Solver could either run on CPU or GPU. The scalability is same as Fourie
 
 For GPU users, here are some useful numbers of memory requirement for choosing the resolution of the simulation. You may end up getting higher resolution for the same memory.
 
-| Memory Size | Maximum Resolution ($N^3$ )   |
-| ----------- | ----------------------------- |
-| 6 GB        | $256^3$ (pure HD simulation)  |
+| Memory Size | Maximum Resolution ( $N^3$ )    |
+| ----------- | ------------------------------ |
+| 6 GB        | $256^3$ (pure HD simulation) |
 | 10 GB       | $300^3$ (pure MHD simulation) |
 
 **Speed**
@@ -72,13 +67,10 @@ Environment: WSL2 in Win11 (Ubuntu 18.04 LTS through jupyter-lab)
 | NVIDIA RTX 3080 10GB        | 0.041s | 0.060s | 0.15s   | 1.23s   |
 
 ## Example
-
 Few examples were set up to illustrate the workflow of using this package. See `example\` for more detail.  The documentation is work in progress and will be available in the future. 
 
 ## Developer
-
 MHDFlows is currently developed by [Ka Wai HO@UW-Madison Astronomy](https://scholar.google.com/citations?user=h2j8wbYAAAAJ&hl=en).
 
 ## Citing
-
 A paper can be cited elsewhere in the future :slightly_smiling_face:. Feel free to cite the GitHub page right now. 

src/MHDFlows.jl

@@ -3,11 +3,11 @@ module MHDFlows
 using 
   CUDA,
   Statistics,
-  SpecialFunctions,
   Reexport,
   DocStringExtensions,
   HDF5,
-  FFTW
+  FFTW,
+  ProgressMeter
 
 @reexport using FourierFlows
 
@@ -28,10 +28,12 @@ include("Solver/MHDSolver_VP.jl")
 include("DiagnosticWrapper.jl")
 include("integrator.jl")
 include("datastructure.jl")
+include("utils/utils.jl");
 include("utils/VectorCalculus.jl")
 include("utils/MHDAnalysis.jl")
 include("utils/GeometryFunction.jl")
 include("utils/IC.jl")
+include("utils/UserInterface.jl")
 
 
 #pgen module

src/Problems.jl

@@ -124,11 +124,12 @@ function MHDFLowsProblem(eqn::FourierFlows.Equation, stepper, dt, grid::Abstract
                  vars=EmptyVars, params=EmptyParams, dev::Device=CPU(); 
                  BFlag = false, VPFlag = false, DyeFlag = false, usr_func = [],
                  stepperkwargs...) where T
+
   clock = FourierFlows.Clock{T}(dt, 0, 0)
 
   timestepper = FourierFlows.TimeStepper(stepper, eqn, dt, dev);
 
-  sol = FourierFlows.devzeros(dev, eqn.T, eqn.dims);
+  sol = zeros(dev, eqn.T, eqn.dims);
 
   flag = Flag(BFlag, VPFlag);
 
@@ -147,7 +148,7 @@ show(io::IO, problem::MHDFlowsProblem) =
     		      "  │     └── user function: "*CheckFunction(problem.usr_func),'\n',
     		      "  │                        ",'\n',
               "  │     Features           ",'\n',  
-              "  │     ├─────────── grid: grid (on " * FourierFlows.griddevice(problem.grid) * ")", '\n',
+              "  │     ├─────────── grid: grid (on " * string(typeof(problem.grid.device)) * ")", '\n',
               "  │     ├───── parameters: params", '\n',
               "  │     ├────── variables: vars", '\n',
               "  └─────├─── state vector: sol", '\n',

src/Solver/HDSolver.jl

@@ -32,28 +32,25 @@ function UᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="x")
   	a    = 1;
   	kₐ   = grid.kr;
   	k⁻²  = grid.invKrsq;
-  	uᵢh  = vars.uxh;
   	∂uᵢh∂t = @view N[:,:,:,params.ux_ind];
 
   elseif direction == "y"
 
   	a    = 2;
   	kₐ   = grid.l;
   	k⁻²  = grid.invKrsq;
-  	uᵢh  = vars.uyh;
   	∂uᵢh∂t = @view N[:,:,:,params.uy_ind];
 
   elseif direction == "z"
 
   	a    = 3;
   	kₐ   = grid.m;
   	k⁻²  = grid.invKrsq;
-  	uᵢh  = vars.uzh;
   	∂uᵢh∂t = @view N[:,:,:,params.uz_ind];
 
   else
 
-  	@warn "Warning : Unknown direction is declerad"
+  	error("Warning : Unknown direction is declerad")
 
   end
 
@@ -80,6 +77,9 @@ function UᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="x")
     end
 
     #Compute the diffusion term  - νk^2 u_i
+    uᵢ = direction == "x" ? vars.ux : direction == "y" ? vars.uy : vars.uz;
+    uᵢh = vars.nonlinh1;
+    mul!(uᵢh, grid.rfftplan, uᵢ); 
     @. ∂uᵢh∂t += -grid.Krsq*params.ν*uᵢh;
 
     # hyperdiffusion term
@@ -93,16 +93,10 @@ end
 
 function HDcalcN_advection!(N, sol, t, clock, vars, params, grid)
 
-  #Update V  Fourier Conponment
-  @. vars.uxh = sol[:, :, :, params.ux_ind];
-  @. vars.uyh = sol[:, :, :, params.uy_ind];
-  @. vars.uzh = sol[:, :, :, params.uz_ind];
-
-
   #Update V + B Real Conponment
-  ldiv!(vars.ux, grid.rfftplan, deepcopy(vars.uxh))
-  ldiv!(vars.uy, grid.rfftplan, deepcopy(vars.uyh))
-  ldiv!(vars.uz, grid.rfftplan, deepcopy(vars.uzh))
+  ldiv!(vars.ux, grid.rfftplan, deepcopy(@view sol[:, :, :, params.ux_ind]));
+  ldiv!(vars.uy, grid.rfftplan, deepcopy(@view sol[:, :, :, params.uy_ind]));
+  ldiv!(vars.uz, grid.rfftplan, deepcopy(@view sol[:, :, :, params.uz_ind]));
 
   #Update V Advection
   UᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="x")
@@ -112,23 +106,4 @@ function HDcalcN_advection!(N, sol, t, clock, vars, params, grid)
   return nothing
 end
 
-function HDupdatevars!(prob)
-  vars, grid, sol, params = prob.vars, prob.grid, prob.sol, prob.params
-
-  dealias!(sol, grid)
-
-  #Update V Fourier Conponment
-  @. vars.uxh = sol[:, :, :, params.ux_ind];
-  @. vars.uyh = sol[:, :, :, params.uy_ind];
-  @. vars.uzh = sol[:, :, :, params.uz_ind];
-
-
-  #Update V + B Real Conponment
-  ldiv!(vars.ux, grid.rfftplan, deepcopy(vars.uxh)) # deepcopy() since inverse real-fft destroys its input
-  ldiv!(vars.uy, grid.rfftplan, deepcopy(vars.uyh)) # deepcopy() since inverse real-fft destroys its input
-  ldiv!(vars.uz, grid.rfftplan, deepcopy(vars.uzh)) # deepcopy() since inverse real-fft destroys its input
-
-  return nothing
-end
-
 end

src/Solver/HDSolver_VP.jl

@@ -29,7 +29,6 @@ function UᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="x")
   	k⁻²  = grid.invKrsq;
     U₀   = params.U₀x;
     uᵢ   = vars.ux;
-  	uᵢh  = vars.uxh;
   	∂uᵢh∂t = @view N[:,:,:,params.ux_ind];
 
   elseif direction == "y"
@@ -39,7 +38,6 @@ function UᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="x")
   	k⁻²  = grid.invKrsq;
     U₀   = params.U₀y;
     uᵢ   = vars.uy;
-  	uᵢh  = vars.uyh;
   	∂uᵢh∂t = @view N[:,:,:,params.uy_ind];
 
   elseif direction == "z"
@@ -49,12 +47,11 @@ function UᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="x")
   	k⁻²  = grid.invKrsq;
     U₀   = params.U₀z;
     uᵢ   = vars.uz;
-  	uᵢh  = vars.uzh;
   	∂uᵢh∂t = @view N[:,:,:,params.uz_ind];
 
   else
 
-  	@warn "Warning : Unknown direction is declerad"
+  	error("Warning : Unknown direction is declerad")
 
   end
 
@@ -98,6 +95,9 @@ function UᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="x")
     end
 
     #Compute the diffusion term  - νk^2 u_i
+    uᵢ = direction == "x" ? vars.ux : direction == "y" ? vars.uy : vars.uz;
+    uᵢh = vars.nonlinh1;
+    mul!(uᵢh, grid.rfftplan, uᵢ); 
     @. ∂uᵢh∂t += -grid.Krsq*params.ν*uᵢh;
 
     # hyperdiffusion term
@@ -111,15 +111,10 @@ end
 
 function HDcalcN_advection!(N, sol, t, clock, vars, params, grid)
 
-  #Update V + B Fourier Conponment
-  @. vars.uxh = sol[:, :, :, params.ux_ind];
-  @. vars.uyh = sol[:, :, :, params.uy_ind];
-  @. vars.uzh = sol[:, :, :, params.uz_ind];
-
-  #Update V + B Real Conponment
-  ldiv!(vars.ux, grid.rfftplan, deepcopy(vars.uxh))
-  ldiv!(vars.uy, grid.rfftplan, deepcopy(vars.uyh))
-  ldiv!(vars.uz, grid.rfftplan, deepcopy(vars.uzh))
+  #Update V Real Conponment
+  ldiv!(vars.ux, grid.rfftplan, deepcopy(@view sol[:, :, :, params.ux_ind]))
+  ldiv!(vars.uy, grid.rfftplan, deepcopy(@view sol[:, :, :, params.uy_ind]))
+  ldiv!(vars.uz, grid.rfftplan, deepcopy(@view sol[:, :, :, params.uz_ind]))
 
   #Update V Advection
   UᵢUpdate!(N, sol, t, clock, vars, params, grid;direction="x")