v 0.2.0 Update

LICENSE

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2022 : Ka Wai HO
+Copyright (c) 2023 : Ka Wai HO
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

Project.toml

@@ -1,7 +1,7 @@
 name = "MHDFlows"
 uuid = "1d939cba-ab73-4bc0-975c-87d4c856e1f9"
 authors = ["Ka Wai HO <woodyho13@gmail.com> "]
-version = "0.1.4"
+version = "0.2.0"
 
 [deps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
@@ -10,18 +10,15 @@ FourierFlows = "2aec4490-903f-5c70-9b11-9bed06a700e1"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
+MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+PyPlot =  "d330b81b-6aea-500a-939a-2ce795aea3ee"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 
 [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"
-TimeOutputs = "^0.5.21"
+julia = "^1.7.3"

README.md

@@ -3,18 +3,20 @@
 
 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:
 
-1. 2D incompressible HD/MHD simulation (periodic boundary)
-2. 3D incompressible HD/MHD simulation (periodic boundary)
+1. 2/3D incompressible HD/MHD simulation (periodic boundary)
 3. Incompressible  HD/MHD simulation with volume penalization method
-4. Passive Dye Tracer (Experimental Feature)
+4. Isothermal compressible  HD/MHD simulation (periodic boundary)
+5. 2/3D Electron magnetohydrodynamic simulation (periodic boundary)
+6. Passive Dye Tracer (Experimental Feature)
 
 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.4
+v 0.2.0
+note : v 0.2.0 will be the final major update before the multi-gpu version release 
 
 ## Installation Guide & compatibility 
-The current version is tested on v1.5.3/1.7.3/1.8.2 version.
+The current version is tested on v1.7.3/1.8.2/1.9.0 version.
 
 Currently, you have two way of installing MHDFlows.jl
 
@@ -24,7 +26,7 @@ Currently, you have two way of installing MHDFlows.jl
 
    ```julia
    julia>
-   (v1.7) pkg> add MHDFlows
+   (v1.8) pkg> add MHDFlows
    ```
 
 
@@ -37,13 +39,14 @@ The MHD Solver could either run on CPU or GPU. The scalability is same as Fourie
 
 **Memory usage**
 
-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.
+For GPU users, here are some useful numbers of memory requirement for choosing the resolution of the simulation with RK4/ LSRK4 method. You may end up getting higher resolution for the same memory.
 
-| Memory Size | Maximum Resolution ( $N^3$ )    |
-| ----------- | ------------------------------ |
+| Memory Size | Maximum Resolution ( $N^3$ )  |
+| ----------- | ------------------------------|
 | 6 GB        | $256^3$ (pure MHD simulation) |
-| 10 GB       | $300^3$ (pure MHD simulation) |
-| 40 GB       | $512^3$ (pure MHD simulation) |
+| 10 GB       | $320^3$ (pure MHD simulation) |
+| 24 GB       | $512^3$ (pure MHD simulation) |
+| 80 GB       | $700^3$ (pure MHD simulation) |
 
 **Speed**
 
@@ -57,18 +60,18 @@ Environment: WSL2 in Win11 (Ubuntu 18.04/20.04 LTS through jupyter-lab)
 
 | Spec CPU/GPU                | $32^3$ | $64^3$ | $128^3$ | $256^3$ |
 | --------------------------- | ------ | ------ | ------- | ------- |
-| AMD Ryzen 7 5800x 8 threads | 0.139s | 0.178s | 0.764s  | 7.025s  |
-| NVIDIA RTX 3080 10GB        | 0.016s | 0.018s | 0.038s  | 0.211s  |
+| AMD Ryzen 7 5800x 8 threads | 0.040s | 0.074s | 0.490S  | 7.025s  |
+| NVIDIA RTX 3080 10GB        | 0.016s | 0.018s | 0.023s  | 0.152s  |
 
 **MHD** (Taylor Green Vortex, T = Float32)
 
 | Spec CPU/GPU                | $32^3$ | $64^3$ | $128^3$ | $256^3$ |
 | --------------------------- | ------ | ------ | ------- | ------- |
-| AMD Ryzen 7 5800x 8 threads | 0.12s  | 0.220s | 1.49s   | 14.50s  |
-| NVIDIA RTX 3080 10GB        | 0.019s | 0.019s | 0.050s  | 0.39 s  |
+| AMD Ryzen 7 5800x 8 threads | 0.049s | 0.180s | 1.490s  | 14.50s  |
+| NVIDIA RTX 3080 10GB        | 0.013s | 0.012s | 0.037s  | 0.271s  |
 
 ## 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. 
+Few examples were set up to illustrate the workflow of using this package. [Check out](https://github.com/MHDFlows/MHDFlows-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).

src/DyeModule.jl

@@ -106,25 +106,25 @@ function DyeEqn!(N, sol, t, clock, vars, params, grid)
     # ∂ρₖ/∂t = ∑ᵢ -im*kᵢ*(ρₖ vₖᵢ)
 
     # Initialization of rho
-    @. N*=0;
+    @. N*=0
 
     for (uᵢ,kᵢ) ∈ zip([vars.ux,vars.uy,vars.uz],[grid.kr,grid.l,grid.m])
          # Initialization
         @. vars.nonlin1 *= 0;
-        ρuᵢ  = vars.nonlin1;
-        ρuᵢh = vars.nonlinh1;
+        ρuᵢ  = vars.nonlin1
+        ρuᵢh = vars.nonlinh1
 
         # get back the updated sol in real space using fft
         ldiv!(ρuᵢ, grid.rfftplan, deepcopy(sol))  
 
         # Pre-Calculation in Real Space
-        @. ρuᵢ *= uᵢ;
+        @. ρuᵢ *= uᵢ
 
         # Fourier transform 
-        mul!(ρuᵢh, grid.rfftplan, ρuᵢ);
+        mul!(ρuᵢh, grid.rfftplan, ρuᵢ)
 
         # Perform the actual calculation
-        @. N += -im*kᵢ*ρuᵢh;
+        @. N += -im*kᵢ*ρuᵢh
     end
     return nothing
 end

src/MHDFlows.jl

@@ -2,15 +2,18 @@ module MHDFlows
 
 using 
   CUDA,
+  FastBroadcast,
   Statistics,
   Reexport,
   DocStringExtensions,
   HDF5,
   FFTW,
-  ProgressMeter
+  ProgressMeter,
+  TimerOutputs
 
 @reexport using FourierFlows
 
+using Random: rand!
 using LinearAlgebra: mul!, ldiv!
 import Base: show, summary
 
@@ -21,29 +24,51 @@ abstract type MHDVars <: AbstractVars end
 include("DyeModule.jl")
 include("Problems.jl")
 include("pgen.jl")
-include("Solver/VPSolver.jl")
-include("Solver/HDSolver.jl")
-include("Solver/MHDSolver.jl")
+# Data Structure
+include("Structure/datastructure.jl")
+include("Structure/HDParams.jl")
+include("Structure/HDVars.jl")
+include("Structure/MHDParams.jl")
+include("Structure/MHDVars.jl")
+# Solver
+include("Solver/VPSolver.jl");
+include("Solver/HDSolver.jl");
+include("Solver/MHDSolver.jl");
+include("Solver/ShearingBox.jl");
+include("Solver/HDSolver_Compessible.jl");
+include("Solver/MHDSolver_Compessible.jl");
+# integrator related
 include("DiagnosticWrapper.jl")
 include("integrator.jl")
-include("datastructure.jl")
+
+# timestepper
+include("timestepper/timestepper.jl")
+
+#utils
 include("utils/utils.jl");
 include("utils/VectorCalculus.jl")
 include("utils/MHDAnalysis.jl")
 include("utils/GeometryFunction.jl")
 include("utils/IC.jl")
 include("utils/UserInterface.jl")
-
+include("utils/TurbStatTool.jl")
 
 #pgen module
+include("pgen/A99ForceDriving_GPU.jl")
 include("pgen/A99ForceDriving.jl")
 include("pgen/TaylorGreenDynamo.jl")
 include("pgen/NegativeDamping.jl")
 
+DivVCorrection! = VPSolver.DivVCorrection!;
+DivBCorrection! = VPSolver.DivBCorrection!;
+
 export Problem,           
        TimeIntegrator!,
        Restart!,
+       DivVCorrection!,
+       DivBCorrection!,
        Cylindrical_Mask_Function,
+       DivFreeSpectraMap,
        SetUpProblemIC!,
        readMHDFlows,
        Curl,            

src/Problems.jl

@@ -3,42 +3,31 @@
 # ----------
 
 
-"""
-    struct Equation{T, TL, G<:AbstractFloat}
-    
-The equation to be solved `∂u/∂t = L*u + N(u)`. Array `L` includes the coefficients
-of the linear term `L*u` and `calcN!` is a function which computes the nonlinear
-term `N(u)`. The struct also includes the problem's `grid` and the float type of the
-state vector (and consequently of `N(u)`).
-$(TYPEDFIELDS)
-"""
-struct Equation{T, TL, G<:AbstractFloat}
-    "array with coefficient for the linear part of the equation"
-       L :: TL
-    "function that computes the nonlinear part of the equation"
-  calcN! :: Function
-    "the grid"
-    grid :: AbstractGrid{G}
-    "the dimensions of `L`"
-    dims :: Tuple
-    "the float type for the state vector"
-       T :: T # eltype or tuple of eltypes of sol and N
-end
-
 """
     Equation(L, calcN!, grid; dims=supersize(L), T=nothing)
     
 The equation constructor from the array `L` of the coefficients of the linear term, the function 
 `calcN!` that computes the nonlinear term and the `grid` for the problem.
 """
-function Equation(L, calcN!, grid::AbstractGrid{G}; dims=supersize(L), T=nothing) where G
+function Setup_Equation(calcN!, grid::AbstractGrid{G}; T=nothing, Nl = 3) where G
+  dims = tuple(size(grid.Krsq)...,Nl)
   T = T == nothing ? T = cxtype(G) : T
-
-  return Equation(L, calcN!, grid, dims, T)
+  #Compatibility to FourierFlows.Equation 
+  L = 0
+  return FourierFlows.Equation(L, calcN!, grid; dims=dims)
 end
 
 CheckON(Flag_equal_to_True::Bool) = Flag_equal_to_True ? string("ON") : string("OFF")
 
+function CheckON(FlagB::Bool, FlagE::Bool) 
+  if FlagB 
+    FlagE ? string("ON (EMHD)") : string("ON (Ideal MHD)")
+  else
+    string("OFF")
+  end
+
+end
+
 function CheckDye(dye::Dye)
     if dye.dyeflag == true
         return string("ON, at prob.dye")
@@ -79,8 +68,14 @@ $(TYPEDFIELDS)
 struct Flag
     "Magnetic Field"
      b :: Bool
+     "EMHD Field"
+     e :: Bool
     "Volume Penalization"
     vp :: Bool
+    "Compressibility"
+     c :: Bool
+     "Shear"
+     s :: Bool
 end
 
 """
@@ -89,15 +84,15 @@ end
 A problem that represents a partial differential equation.
 $(TYPEDFIELDS)
 """
-struct MHDFlowsProblem{T, A<:AbstractArray, Tg<:AbstractFloat, TL, Dye, usr_foo} <: AbstractProblem
+struct MHDFlowsProblem{T, A<:AbstractArray, Tg<:AbstractFloat, TL, Dye, usr_foo, AbstractGrid} <: AbstractProblem
     "the state vector"
           sol :: A
     "the problem's slock"
         clock :: FourierFlows.Clock{Tg}
     "the equation"
           eqn :: FourierFlows.Equation{T, TL, Tg}
     "the grid"
-         grid :: AbstractGrid{Tg}
+         grid :: AbstractGrid
     "the variables"
          vars :: AbstractVars
     "the parameters"
@@ -122,36 +117,43 @@ to the time-stepper constructor.
 """
 function MHDFLowsProblem(eqn::FourierFlows.Equation, stepper, dt, grid::AbstractGrid{T}, 
                  vars=EmptyVars, params=EmptyParams, dev::Device=CPU(); 
-                 BFlag = false, VPFlag = false, DyeFlag = false, usr_func = [],
+                 BFlag = false, EFlag = false, VPFlag = false, CFlag = false, SFlag = false, DyeFlag = false, usr_func = [],
                  stepperkwargs...) where T
 
   clock = FourierFlows.Clock{T}(dt, 0, 0)
+  if EFlag && stepper == "HM89"
+  #  timestepper = eSSPIFRK3TimeStepper(eqn, dev) #For SFlag
+    timestepper = HM89TimeStepper(eqn, dev) #For EFlag
+  else
+    timestepper = FourierFlows.TimeStepper(stepper, eqn, dt, dev)
+  end
+  sol = zeros(dev, eqn.T, eqn.dims)
 
-  timestepper = FourierFlows.TimeStepper(stepper, eqn, dt, dev);
-
-  sol = zeros(dev, eqn.T, eqn.dims);
+  flag = Flag(BFlag, EFlag, VPFlag, CFlag, SFlag)
 
-  flag = Flag(BFlag, VPFlag);
+  dye = DyeContructer(dev, DyeFlag, grid)
 
-  dye = DyeContructer(dev, DyeFlag, grid);
+  usr_func = length(usr_func) == 0 ? [nothingfunction] : usr_func
 
-  usr_func = length(usr_func) == 0 ? [nothingfunction] : usr_func;
   return MHDFlowsProblem(sol, clock, eqn, grid, vars, params, timestepper, flag, usr_func, dye)
+
 end
 
 show(io::IO, problem::MHDFlowsProblem) =
     print(io, "MHDFlows Problem\n",
-    	        "  │    Funtions\n",
-    		      "  │     ├──────── B-field: "*CheckON(problem.flag.b),'\n',
-    		      "  ├─────├────── VP Method: "*CheckON(problem.flag.vp),'\n',
-    		      "  │     ├──────────── Dye: "*CheckDye(problem.dye),'\n',
-    		      "  │     └── user function: "*CheckFunction(problem.usr_func),'\n',
-    		      "  │                        ",'\n',
-              "  │     Features           ",'\n',  
-              "  │     ├─────────── grid: grid (on " * string(typeof(problem.grid.device)) * ")", '\n',
-              "  │     ├───── parameters: params", '\n',
-              "  │     ├────── variables: vars", '\n',
-              "  └─────├─── state vector: sol", '\n',
-              "        ├─────── equation: eqn", '\n',
-              "        ├────────── clock: clock", '\n',
-              "        └──── timestepper: ", string(nameof(typeof(problem.timestepper))))
+    	    "  │    Funtions\n",
+          "  │     ├ Compressibility: "*CheckON(problem.flag.c),'\n',
+          "  │     ├──────── B-field: "*CheckON(problem.flag.b,problem.flag.e),'\n',
+          "  │     ├────────── Shear: "*CheckON(problem.flag.s),'\n',
+    		  "  ├─────├────── VP Method: "*CheckON(problem.flag.vp),'\n',
+    		  "  │     ├──────────── Dye: "*CheckDye(problem.dye),'\n',
+    		  "  │     └── user function: "*CheckFunction(problem.usr_func),'\n',
+    		  "  │                        ",'\n',
+          "  │     Features           ",'\n',  
+          "  │     ├─────────── grid: grid (on " * string(typeof(problem.grid.device)) * ")", '\n',
+          "  │     ├───── parameters: params", '\n',
+          "  │     ├────── variables: vars", '\n',
+          "  └─────├─── state vector: sol", '\n',
+          "        ├─────── equation: eqn", '\n',
+          "        ├────────── clock: clock", '\n',
+          "        └──── timestepper: ", string(nameof(typeof(problem.timestepper))))