diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md new file mode 100644 index 0000000..ba1dbb5 --- /dev/null +++ b/.github/pull_request_template.md @@ -0,0 +1,16 @@ +*We set up this PR template to confirm that you have considered all the necessary aspects of your contribution so that we are on the right track. Please fill out the template below to the best of your ability. If you have any questions, please feel free to get in touch.* + +## Proposed Changes +*Give a brief overview of your contribution here in a few sentences. For bug fixes, please detail the bug and include a test case. For a new feature, please clearly describe the design, its rationale, the possible alternatives considered. Note that each PR should address only one issue* + +## Related Work +*Resolve any issues (bug fix or feature request), note any related PRs, or mention interactions with the work of others, if any.* + +## Checklist +*Put an x by all that apply. You can fill this out after submitting the PR.* + +- [ ] I am submitting my contribution to the development branch. +- [ ] My contribution is commented and uses the code style of the project (see .JuliaFormatter.toml). +- [ ] I have added a unit test to verify my contribution, if necessary. +- [ ] I have updated appropriate documentation (README, documentation, tutorials, etc.), if applicable. +- [ ] I have added an entry in NEWS.md, if necessary. \ No newline at end of file diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 8b57e61..28c458b 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -32,8 +32,7 @@ jobs: arch: - x64 steps: - - run: | - pip3 install meshio scipy --user + - run: pip3 install scipy --user - uses: actions/checkout@v2 - uses: julia-actions/setup-julia@v1 with: @@ -56,6 +55,8 @@ jobs: - uses: julia-actions/julia-runtest@latest continue-on-error: ${{ matrix.version == 'nightly' }} - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v1 + - uses: codecov/codecov-action@v4 with: - file: lcov.info \ No newline at end of file + file: lcov.info + token: ${{ secrets.CODECOV_TOKEN }} + fail_ci_if_error: true \ No newline at end of file diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml index 07a31f6..45b30ed 100644 --- a/.github/workflows/format.yml +++ b/.github/workflows/format.yml @@ -2,6 +2,7 @@ name: FormatPR on: schedule: - cron: '0 0 * * *' +permissions: write-all jobs: build: runs-on: ubuntu-latest @@ -28,4 +29,4 @@ jobs: - name: Check outputs run: | echo "Pull Request Number - ${{ steps.cpr.outputs.pull-request-number }}" - echo "Pull Request URL - ${{ steps.cpr.outputs.pull-request-url }}" \ No newline at end of file + echo "Pull Request URL - ${{ steps.cpr.outputs.pull-request-url }}" diff --git a/Project.toml b/Project.toml index 6064281..8fc424c 100644 --- a/Project.toml +++ b/Project.toml @@ -7,24 +7,39 @@ version = "0.1.0" AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" -JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" -Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -MAT = "23992714-dd62-5051-b70f-ba57cb901cac" -MKL_jll = "856f044c-d86e-5d09-b602-aeab76dc8ba7" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" -MPICH_jll = "7cb0a576-ebde-5e09-9194-50597f1243b4" P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" -ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" -SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -Tar = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" + +[compat] +AbstractTrees = "0.4" +CairoMakie = "0.12, 0.13, 0.14" +DocStringExtensions = "0.9" +JLD2 = "0.4, 0.5" +LinearAlgebra = "<0.0.1, 1" +MPI = "0.20" +P4est = "0.4" +Parameters = "0.12" +PrecompileTools = "^1" +PythonCall = "0.9" +RecipesBase = "^1" +Reexport = "^1" +SparseArrays = "1.10" +SpecialFunctions = "^2" +StaticArrays = "1.9" +julia = "^1" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] diff --git a/README.md b/README.md index 9a9b2cb..aab3e56 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,111 @@ # KitAMR.jl -Large-scale distributed computational fluid dynamics computing with adaptive mesh refinement + +![CI](https://img.shields.io/github/actions/workflow/status/vavrines/KitBase.jl/ci.yml?branch=main) +[![codecov](https://img.shields.io/codecov/c/github/CFDML/KitAMR.jl)](https://codecov.io/gh/CFDML/KitAMR.jl) + +## Motivation +**KitAMR.jl** is a distributed adaptive Cartesian grid solver for kinetic equations, developed based on [P4est.jl](https://github.com/trixi-framework/P4est.jl). Its goal is to achieve large-scale parallel solving of 2D and 3D flows across all regimes, leveraging GPUs, differentiable programming, and machine learning to enhance solving efficiency. + +Currently, the development is in its initial stages, with the following functionalities implemented: +- 2D flow +- Square domain grid generation +- Geometric adaptation of physical space grids +- Dynamic adaptation of physical space grids +- Dynamic adaptation of velocity space grids +- Load balancing based on physical space grid density +- VanLeer reconstruction +- UGKS fluxes +- Post-processing + +## Theory +### Kinetic Methods +In kinetic theory, methods simulate the macroscopic motion of fluids by describing the evolution of distribution functions of fluid particles in velocity space over time. Such methods include direct simulation Monte Carlo (DSMC), discrete velocity method (DVM), lattice Boltzmann equation (LBE), gas-kinetic scheme, semi-Lagrangian method, implicit-explicit (IMEX) method, and others. +### Distribution Function +In kinetic methods, the state of the fluid is described by a distribution function $f(\mathbf{x},\mathbf{u},t)$. The distribution function represents the number of particles at a given moment in time at a specific physical space point in a particular velocity space element. Its normalization property is + +$$ +\rho(\mathbf{x},t)=\iiint_{-\infty}^{\infty}f(\mathbf{x},\mathbf{u},t)d\mathbf{u} +$$ + +, where $\rho(\mathbf{x},t)$ is the density of flows at $\mathbf{x}$ in time $t$. At higher temperatures, molecular internal degrees of freedom are excited, and the distribution function will depend on these degrees of freedom as well. For simplicity, we do not currently consider this situation. + +
+ Image 1 +
+With the distribution function, macroscopic conservative variables can be caculated as follows: + +$$ +\begin{split} +&\rho = \iiint_{-\infty}^{\infty}f(\mathbf{x},\mathbf{u},t)d\mathbf{u}\\ +&\rho\mathbf U = \iiint_{-\infty}^{\infty}\mathbf uf(\mathbf{x},\mathbf{u},t)d\mathbf{u}\\ +&\rho E = \iiint_{-\infty}^{\infty} \frac 12 \mathbf u^2f(\mathbf{x},\mathbf{u},t)d\mathbf{u}\\ +&\mathbf q=\iiint_{-\infty}^{\infty} \mathbf u\frac 12 \mathbf u^2f(\mathbf{x},\mathbf{u},t)d\mathbf{u} +\end{split} +$$ + +### Boltzmann Equation +The evolution of the distribution function follows the Boltzmann equation. + +$$ +\frac{\partial f}{\partial t}+\mathbf{u}\cdot\frac{\partial f}{\partial \mathbf{x}}+\mathbf{F}\cdot\frac{\partial f}{\partial \mathbf{u}}=\iiint_{-\infty}^{\infty}\int_0^{4\pi}(f^{\*}f_1^{\*}-ff_1)c_r\sigma d\Omega dc_1 +$$ + + +The Boltzmann equation holds a central position in kinetic theory, which is an integro-differential equation, with a four-fold integral term on the right-hand side. The complexity introduced by this term makes solving the Boltzmann equation challenging. Therefore, as a simplification, kinetic methods often solve its model equations instead. +### Shakhov Model Equation +The simplified equation obtained after simplifying the collision term is known as the Shakhov model equation: + +$$ +\begin{cases} +\frac{\partial f}{\partial t}+\mathbf{u}\cdot\frac{\partial f}{\partial \mathbf{x}}+\mathbf{F}\cdot\frac{\partial f}{\partial \mathbf{u}}=\frac{g^+-f}{\tau}\\ +g^+=g[1+\frac45(1-\mathrm{Pr})\lambda^2\frac{\mathbf{u}\cdot\mathbf{q}}{\rho}(2\lambda \mathbf{u}^2-5)] +\end{cases} +$$ + +, Where $\mathrm {Pr}$ is the Prandtl number, $\mathbf{q}$ is the heat flux, and $\tau$ is the relaxation time typically given through phenomenological models to match experimental data. The Shakhov model can provide results close to the Boltzmann equation for mass, momentum, energy, and heat fluxes in situations where flow velocities are not very high and deviations from equilibrium are not very large. And the current solver is based on Shakhov model. + +## Computing Methods +### Finite Volume Method (FVM) +The underlying discretization of the Boltzmann equation is based on finite volume method, with its discrete form as follows: + +$$ +\{\bar U_j\Omega_j\}^{n+1} = +\{\bar U_j \Omega_j\}^n-\Delta t\sum_{\partial\Omega}\mathbf{F}^{\*}\cdot\Delta\mathbf{S}+\Delta t \bar{Q}_j\Omega_j +$$ + +, where $\bar U_j^n\coloneqq\frac{1}{\Omega_j}\int_{\Omega_j}Ud\Omega_j|^n$, $\bar Q_j^n\coloneqq\frac{1}{\Omega_j}\int_{\Omega_j}Qd\Omega_j$ are the cell-averaged conservative variable and source. And $\bar Q_j^{\*}$ and $\mathbf{F}^{\*}$ are respectively cell- and time-averaged sources and numerical flux, which are defined as $\bar Q_j^{\*}\coloneqq\frac{1}{\Delta t}\int_n^{n+1}\bar Q_jdt$ and $\mathbf F^{\*}\cdot \Delta \mathbf{S}\coloneqq \frac{1}{\Delta t}\int_n^{n+1}\mathbf F\cdot\Delta \mathbf{S}dt$. +It is an exact relation for the time evolution of the space averaged conservative avriables over cell $j$ from time step $n$ to $n+1$. And the numerical flux $\mathbf F^{\*}$ identifies completely a scheme by the way it approximates the time-averaged physical flux along each cell face. +### Numerical Flux +The numerical scheme that the solver currently supports is UGKS, with interface variables obtained from VanLeer reconstruction. For specific implementation details, please refer to the relevant literature. + +### Discrete Velocity Method +DVM is one of the popular approaches for solving rarefied flow problems. In this method, the particle velocity space is discretized into a finite set of points and the numerical quadrature rule is utilized to approximate the integration of moments. Since the particle velocity space is discretized into a finite set of points, the continuum Boltzmann equation is reduced to the corresponding discrete velocity Boltzmann equation (DVBE). A great variety of algorithms have been developed, including unified gas-kinetic scheme (UGKS), discrete unified gas-kinetic scheme (DUGKS), semi-Lagrangian method, etc. Overall, most of the above methods can be applied from free molecular regime to continuum regime. But in order to make the quadrature error to be small enough, a large number of discrete velocity points are usually required. In particular, for fluid flows near continuum regime, the computational cost of DVM is much larger than those traditional CFD methods based on the Navier-Stokes equation. + +### Adaptive Mesh Refinement (AMR) +Considering the above limitations of DVM, we adopt Adaptive Mesh Refinement (AMR) to improve solving efficiency. AMR reallocates computational resources based on flow features, balancing efficiency and accuracy. + +
+ Image 1 +
+ +
+ Image 1 + Image 2 +
+ + +The current solver discretizes with tree-based Cartesian grids, adapting the mesh in both physical and velocity space simultaneously, according to characteristics including velocity field gradients, boundary complexity, and energy proportion at phase points. Users can also customize the criteria for grid refinement or coarsening according to their own needs. + +### Cartesian Grids +Cartesian grids offer advantages such as high grid quality and automation but struggle with complex boundary shapes. In the subsequent development of the solver, we will introduce Immersed Boundary Method (IBM) to address this challenge. + +## Installation +The current solver does not yet provide a user interface. If you'd like to try it out, we recommend cloning the code repository from Git to your local machine and then using the Julia package manager to install its dependencies in the directory where the code resides: +```julia +julia> ] +(v1.10) pkg> activate . +(//your_path//) pkg> instantiate +``` +Soon, we will provide a complete user interface with customizable boundaries, and we will register it in the [Julia package registry](https://github.com/JuliaRegistries/General). Stay tuned! +## Current Usage +Currently, flow properties, initial contitions and boundary conditions are defined in `Global_Data` in `/src/types.jl`. It should be noted specifically that the boundary follows the sequence of -x, +x, -y, +y. The maximum level of the refinement is defined in `/src/abstract.jl` by `DVM_PS_MAXLEVEL` and `DVM_VS_MAXLEVEL` respectively. Grid refinement and coarsening behaviors are defined in `adaptive.jl` and `vs_adaptive.jl`. Feel free to modify them according to your needs. \ No newline at end of file diff --git a/src/DVM.jl b/deprecated/DVM.jl similarity index 83% rename from src/DVM.jl rename to deprecated/DVM.jl index 19b7e48..36fd9ea 100644 --- a/src/DVM.jl +++ b/deprecated/DVM.jl @@ -1,12 +1,15 @@ -using Pkg -Pkg.activate(@__DIR__) +#using Pkg +#Pkg.activate(@__DIR__) using MPI using P4est -using P4estTypes +#using P4estTypes +include("../lib/P4estTypes/src/P4estTypes.jl") +using .P4estTypes + using StaticArrays using Parameters using SpecialFunctions -using Makie +#using Makie using CairoMakie using PythonCall scipy = pyimport("scipy") diff --git a/src/flux.jl b/deprecated/flux.jl similarity index 100% rename from src/flux.jl rename to deprecated/flux.jl diff --git a/example/cavity_sim.jl b/example/cavity_sim.jl new file mode 100644 index 0000000..7d48de0 --- /dev/null +++ b/example/cavity_sim.jl @@ -0,0 +1,90 @@ +""" +To execute the code (assuming 16 cores are going to be used): +- from the main directory: mpiexecjl -n 16 --project=. julia example/cavity_sim.jl +- from the file directory: mpiexecjl -n 16 --project=.. julia cavity_sim.jl +""" + +using KitAMR, MPI, CairoMakie + +MPI.Init() +ps4est, DVM_data = KitAMR.init() +if MPI.Comm_rank(MPI.COMM_WORLD) == 0 + KitAMR.save_set(DVM_data) +end +KitAMR.save_result(DVM_data) +while true + if DVM_data.global_data.adapt_step == 10 + KitAMR.update_slope!(DVM_data) + KitAMR.update_gradmax!(DVM_data) + KitAMR.PS_refine!(ps4est) + KitAMR.PS_coarsen!(ps4est) + KitAMR.PS_balance!(ps4est) + if DVM_data.global_data.partition_step == 30 + KitAMR.PS_partition!(ps4est, DVM_data) + DVM_data.global_data.partition_step = 0 + end + if DVM_data.global_data.vs_adapt_step == 50 + KitAMR.vs_refine!(DVM_data) + KitAMR.vs_coarsen!(DVM_data) + DVM_data.global_data.vs_adapt_step = 0 + end + KitAMR.update_ghost!(ps4est, DVM_data) + KitAMR.update_neighbor!(ps4est, DVM_data) + KitAMR.update_faces!(ps4est, DVM_data) + DVM_data.global_data.adapt_step = 0 + end + KitAMR.update_Δt!(DVM_data) # 798.125 μs (41515 allocations: 806.91 KiB) + KitAMR.update_slope!(DVM_data) # 1.304 s (22645678 allocations: 1.43 GiB) -> final: 336.005 ms (32818 allocations: 1.35 MiB) + KitAMR.DVM_data_exchange!(ps4est, DVM_data, Val(2)) # 0.121942 + KitAMR.update_flux!(DVM_data) # 4.437 s (174354747 allocations: 5.30 GiB)->3.268 s (75782427 allocations: 4.84 GiB) -> final: 540.574 ms (1301847 allocations: 1.42 GiB)->411.685 ms (1262009 allocations: 1.51 GiB) + KitAMR.update_volume!(DVM_data) # 1.262 s (59874178 allocations: 1.55 GiB) ->469.264 ms (9070978 allocations: 623.28 MiB)-> final 50.184 ms (248053 allocations: 158.25 MiB) + KitAMR.DVM_data_exchange!(ps4est, DVM_data, Val(1)) + DVM_data.global_data.adapt_step += 1 + DVM_data.global_data.vs_adapt_step += 1 + DVM_data.global_data.partition_step += 1 + KitAMR.save_result(DVM_data) + if MPI.Comm_rank(MPI.COMM_WORLD) == 0 + @show DVM_data.global_data.gas.sim_time + pp = PointerWrapper(ps4est) + @show pp.local_num_quadrants[] + end + DVM_data.global_data.gas.sim_time > DVM_data.data_set.end_time && break + MPI.Barrier(MPI.COMM_WORLD) + if MPI.Comm_rank(MPI.COMM_WORLD) == MPI.Comm_size(MPI.COMM_WORLD) - 1 + ps_data = DVM_data.trees.data[end][end] + @show ps_data.vs_data.vs_num + end +end +KitAMR.save_VS_final(DVM_data) +meshes = KitAMR.collect_mesh(DVM_data) +solutions = KitAMR.collect_solution(ps4est, DVM_data) + +if MPI.Comm_rank(MPI.COMM_WORLD) == 0 + x, y = KitAMR.reshape_mesh(meshes) + f = Figure() + ax = Axis(f[1, 1]) + KitAMR.mesh_plot(x, y, ax) + save("mesh.png", f) + + x, y, variable = KitAMR.reshape_solutions(solutions, DVM_data.global_data, :prim, 4) + variable = 1 ./ variable + f = Figure() + ax = Axis(f[1, 1]) + co = contourf!(x, y, variable) + Colorbar(f[1, 2], co) + save("test_vs_adaptive.png", f) +end + +KitAMR.finalize!(ps4est, DVM_data) + +phase_cells = 0 +for i in eachindex(DVM_data.trees.data) + for j in eachindex(DVM_data.trees.data[i]) + global phase_cells + phase_cells += DVM_data.trees.data[i][j].vs_data.vs_num + end +end +phase_cells = MPI.Allreduce(phase_cells, +, MPI.COMM_WORLD) +if MPI.Comm_rank(MPI.COMM_WORLD) == 0 + @show phase_cells +end diff --git a/example/main_simulation.jl b/example/main_simulation.jl deleted file mode 100644 index 2d009bb..0000000 --- a/example/main_simulation.jl +++ /dev/null @@ -1,83 +0,0 @@ -include("DVM.jl") -MPI.Init() -ps4est, DVM_data = init() -if MPI.Comm_rank(MPI.COMM_WORLD) == 0 - save_set(DVM_data) -end -save_result(DVM_data) -while true - if DVM_data.global_data.adapt_step == 10 - update_slope!(DVM_data) - update_gradmax!(DVM_data) - PS_refine!(ps4est) - PS_coarsen!(ps4est) - PS_balance!(ps4est) - if DVM_data.global_data.partition_step == 30 - PS_partition!(ps4est, DVM_data) - DVM_data.global_data.partition_step = 0 - end - if DVM_data.global_data.vs_adapt_step == 50 - vs_refine!(DVM_data) - vs_coarsen!(DVM_data) - DVM_data.global_data.vs_adapt_step = 0 - end - update_ghost!(ps4est, DVM_data) - update_neighbor!(ps4est, DVM_data) - update_faces!(ps4est, DVM_data) - DVM_data.global_data.adapt_step = 0 - end - update_Δt!(DVM_data) # 798.125 μs (41515 allocations: 806.91 KiB) - update_slope!(DVM_data) # 1.304 s (22645678 allocations: 1.43 GiB) -> final: 336.005 ms (32818 allocations: 1.35 MiB) - DVM_data_exchange!(ps4est, DVM_data, Val(2)) # 0.121942 - update_flux!(DVM_data) # 4.437 s (174354747 allocations: 5.30 GiB)->3.268 s (75782427 allocations: 4.84 GiB) -> final: 540.574 ms (1301847 allocations: 1.42 GiB)->411.685 ms (1262009 allocations: 1.51 GiB) - update_volume!(DVM_data) # 1.262 s (59874178 allocations: 1.55 GiB) ->469.264 ms (9070978 allocations: 623.28 MiB)-> final 50.184 ms (248053 allocations: 158.25 MiB) - DVM_data_exchange!(ps4est, DVM_data, Val(1)) - DVM_data.global_data.adapt_step += 1 - DVM_data.global_data.vs_adapt_step += 1 - DVM_data.global_data.partition_step += 1 - save_result(DVM_data) - if MPI.Comm_rank(MPI.COMM_WORLD) == 0 - @show DVM_data.global_data.gas.sim_time - pp = PointerWrapper(ps4est) - @show pp.local_num_quadrants[] - end - DVM_data.global_data.gas.sim_time > DVM_data.data_set.end_time && break - MPI.Barrier(MPI.COMM_WORLD) - if MPI.Comm_rank(MPI.COMM_WORLD) == MPI.Comm_size(MPI.COMM_WORLD) - 1 - ps_data = DVM_data.trees.data[end][end] - @show ps_data.vs_data.vs_num - end -end -save_VS_final(DVM_data) -meshes = collect_mesh(DVM_data) -solutions = collect_solution(ps4est, DVM_data) -# write_result!(ps4est,DVM_data,"write_test") -if MPI.Comm_rank(MPI.COMM_WORLD) == 0 - x, y = reshape_mesh(meshes) - f = Figure() - ax = Axis(f[1, 1]) - mesh_plot(x, y, ax) - save("mesh.png", f) - - x, y, variable = reshape_solutions(solutions, DVM_data.global_data, :prim, 4) - variable = 1 ./ variable - f = Figure() - ax = Axis(f[1, 1]) - co = contourf!(x, y, variable) - Colorbar(f[1, 2], co) - save("test_vs_adaptive.png", f) -end - -finalize!(ps4est, DVM_data) -# # mpiexecjl -n 6 --project=. julia ./main.jl 2592.10s user 308.70s system 552% cpu 8:45.11 total -phase_cells = 0 -for i in eachindex(DVM_data.trees.data) - for j in eachindex(DVM_data.trees.data[i]) - global phase_cells - phase_cells += DVM_data.trees.data[i][j].vs_data.vs_num - end -end -phase_cells = MPI.Allreduce(phase_cells, +, MPI.COMM_WORLD) -if MPI.Comm_rank(MPI.COMM_WORLD) == 0 - @show phase_cells -end diff --git a/lib/P4estTypes/src/pxest.jl b/lib/P4estTypes/src/pxest.jl index f482902..d37960e 100644 --- a/lib/P4estTypes/src/pxest.jl +++ b/lib/P4estTypes/src/pxest.jl @@ -1051,6 +1051,7 @@ The keyword arguments (`kw...`) for the balancing are: + `P4estTypes.CONNECT_EDGE(Val(8))`: enforce across face and edge. + `P4estTypes.CONNECT_CORNER(Val(4))`: enforce across face and corner. + `P4estTypes.CONNECT_CORNER(Val(8))`: enforce across face, edge, and corner. + - `init = nothing`: callback function with prototype `init(forest, treeid, quadrant)` called for each quadrant created to initialized the user data. diff --git a/precompile.sh b/precompile.sh new file mode 100644 index 0000000..5fc0265 --- /dev/null +++ b/precompile.sh @@ -0,0 +1,2 @@ +echo Triggering precompilation +julia --project=. -e 'using Pkg; Pkg.precompile(); println("Precompilation complete");' \ No newline at end of file diff --git a/src/KitAMR.jl b/src/KitAMR.jl index 1af49dc..5b81b2a 100644 --- a/src/KitAMR.jl +++ b/src/KitAMR.jl @@ -1,22 +1,20 @@ module KitAMR +using CairoMakie using MPI -using P4est -using StaticArrays +using JLD2 +using LinearAlgebra using Parameters using SpecialFunctions -using CairoMakie -using LinearAlgebra -using HDF5 -using JLD2 +using StaticArrays +using PythonCall + +using Reexport +@reexport using P4est include("../lib/P4estTypes/src/P4estTypes.jl") using .P4estTypes -using PythonCall -scipy = pyimport("scipy") -np = pyimport("numpy") - include("abstract.jl") include("types.jl") include("model.jl") @@ -37,4 +35,11 @@ include("partition.jl") include("postprocess.jl") include("finalize.jl") +const np = Ref{Py}() +const scipy = Ref{Py}() +function __init__() + np[] = pyimport("numpy") + scipy[] = pyimport("scipy") +end + end # module diff --git a/src/postprocess.jl b/src/postprocess.jl index 7078d17..23a8a72 100644 --- a/src/postprocess.jl +++ b/src/postprocess.jl @@ -74,12 +74,12 @@ function reshape_solutions( midpoints[i, :] .= solutions[i].midpoint z[i] = getfield(solutions[i], field)[index] end - itp = scipy.interpolate.CloughTocher2DInterpolator(midpoints, z) + itp = scipy[].interpolate.CloughTocher2DInterpolator(midpoints, z) dx = (xmax - xmin) / Nx / 2^DVM_PS_MAXLEVEL dy = (ymax - ymin) / Ny / 2^DVM_PS_MAXLEVEL X = collect(xmin+EPS+dx/2:dx:xmax-EPS-dx/2) Y = collect(ymin+EPS+dy/2:dy:ymax-EPS-dy/2) - pyZ = np.zeros((length(X), length(Y))) + pyZ = np[].zeros((length(X), length(Y))) for i in eachindex(X) for j in eachindex(Y) pyZ[i-1, j-1] = itp(X[i], Y[j])[] @@ -495,12 +495,12 @@ function reshape_solutions_vs(midpoint::AM, df::AM, global_data::Global_Data) Nx, Ny = global_data.vs_trees_num xmin, xmax, ymin, ymax = global_data.quadrature z = @view(df[:, 1]) - itp = scipy.interpolate.CloughTocher2DInterpolator(midpoint, PyArray(z)) + itp = scipy[].interpolate.CloughTocher2DInterpolator(midpoint, PyArray(z)) dx = (xmax - xmin) / Nx / 2^DVM_VS_MAXLEVEL dy = (ymax - ymin) / Ny / 2^DVM_VS_MAXLEVEL X = collect(xmin+EPS+dx/2:dx:xmax-EPS-dx/2) Y = collect(ymin+EPS+dy/2:dy:ymax-EPS-dy/2) - pyZ = np.zeros((length(X), length(Y))) + pyZ = np[].zeros((length(X), length(Y))) for i in eachindex(X) for j in eachindex(Y) pyZ[i-1, j-1] = itp(X[i], Y[j])[] diff --git a/src/types.jl b/src/types.jl index f575ddf..e913e35 100644 --- a/src/types.jl +++ b/src/types.jl @@ -228,7 +228,7 @@ struct Data_Set VS_end_time::Float64 end function Data_Set() - return Data_Set(0.1, 0.05, 3.0, 0.1) + return Data_Set(0.1, 0.05, 0.02, 0.1) end mutable struct DVM_Data global_data::Global_Data diff --git a/test/runtests.jl b/test/runtests.jl index 083c51a..d942feb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1 +1,40 @@ -using KitAMR +using KitAMR, MPI + +MPI.Init() +ps4est, DVM_data = KitAMR.init(); +if MPI.Comm_rank(MPI.COMM_WORLD) == 0 + KitAMR.save_set(DVM_data) +end + +# minimal solution algorithm +begin + KitAMR.update_slope!(DVM_data) + KitAMR.update_gradmax!(DVM_data) + KitAMR.PS_refine!(ps4est) + KitAMR.PS_coarsen!(ps4est) + KitAMR.PS_balance!(ps4est) + KitAMR.PS_partition!(ps4est, DVM_data) + KitAMR.update_ghost!(ps4est, DVM_data) + KitAMR.update_neighbor!(ps4est, DVM_data) + KitAMR.update_faces!(ps4est, DVM_data) + KitAMR.update_Δt!(DVM_data) + KitAMR.update_slope!(DVM_data) + KitAMR.DVM_data_exchange!(ps4est, DVM_data, Val(2)) + KitAMR.update_flux!(DVM_data) + KitAMR.update_volume!(DVM_data) + KitAMR.DVM_data_exchange!(ps4est, DVM_data, Val(1)) + KitAMR.save_result(DVM_data) + MPI.Barrier(MPI.COMM_WORLD) +end + +KitAMR.save_VS_final(DVM_data) +meshes = KitAMR.collect_mesh(DVM_data) +solutions = KitAMR.collect_solution(ps4est, DVM_data) + +# skip plotting +if MPI.Comm_rank(MPI.COMM_WORLD) == 0 + x, y = KitAMR.reshape_mesh(meshes) + x, y, variable = KitAMR.reshape_solutions(solutions, DVM_data.global_data, :prim, 4) +end + +KitAMR.finalize!(ps4est, DVM_data)