diff --git a/CITATION.cff b/CITATION.cff index 23cbac3..20b0311 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -1,32 +1,29 @@ cff-version: 1.2.0 -message: "If you use this software, please cite it as below." +message: "If you find MaterialPointSolver.jl useful or have used it in your research, please cite it as follows." authors: -- family-names: "Huo" - given-names: "Zenan" - orcid: "https://orcid.org/0000-0000-0000-0000" -- family-names: "Bot" - given-names: "Hew" - orcid: "https://orcid.org/0000-0000-0000-0000" -title: "My Research Software" -version: 2.0.4 -doi: 10.5281/zenodo.1234 -date-released: 2017-12-18 -url: "https://github.com/github-linguist/linguist" -preferred-citation: - type: article - authors: - - family-names: "Lisa" - given-names: "Mona" - orcid: "https://orcid.org/0000-0000-0000-0000" - - family-names: "Bot" - given-names: "Hew" - orcid: "https://orcid.org/0000-0000-0000-0000" - doi: "10.0000/00000" - journal: "Journal Title" - month: 9 - start: 1 # First page number - end: 10 # Last page number - title: "My awesome research software" - issue: 1 - volume: 1 - year: 2021 \ No newline at end of file + - family-names: Huo + given-names: Zenan + - family-names: Alkhimenkov + given-names: Yury + - family-names: Jaboyedoff + given-names: Michel + - family-names: Podladchikov + given-names: Yury + - family-names: Räss + given-names: Ludovic + - family-names: Wyser + given-names: Emmanuel + - family-names: Mei + given-names: Gang +title: "A high-performance backend-agnostic Material Point Method solver in Julia" +journal: "Computers and Geotechnics" +volume: 183 +start: 107189 # 对应 pages 的起始页码 +year: 2025 +doi: "10.1016/j.compgeo.2025.107189" +url: "https://www.sciencedirect.com/science/article/pii/S0266352X25001387" +keywords: + - "MPM" + - "Julia language" + - "Heterogeneous computing" + - "Effective memory throughput" \ No newline at end of file diff --git a/Project.toml b/Project.toml index 35f8ca0..9cbf725 100644 --- a/Project.toml +++ b/Project.toml @@ -1,27 +1,32 @@ name = "MaterialPointSolver" uuid = "dc4aa359-7d89-437d-91bb-f4b330c522c1" authors = ["Zenan Huo "] -version = "0.4.5" +version = "0.4.6" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" +Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" SysInfo = "90a7ee08-a23f-48b9-9006-0e0e2a9e4608" [compat] Adapt = "4.0" BenchmarkTools = "1.5" +CodecZlib = "0.7" DelimitedFiles = "1.9" HDF5 = "0.17" KernelAbstractions = "0.9" ProgressMeter = "1.10" +StatsBase = "0.34" Suppressor = "0.2" SysInfo = "0.3" julia = "1.11" @@ -31,12 +36,14 @@ AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" Metal = "dde4c033-4e86-420c-a63e-0dd931031962" oneAPI = "8f75cd03-7ff8-4ecb-9b8f-daf728133b1b" +WGLMakie = "276b4fcb-3e11-5398-bf8b-a0c2d153d008" [extensions] MaterialPointSolverAMDGPUExt = "AMDGPU" MaterialPointSolverCUDAExt = "CUDA" MaterialPointSolverMetalExt = "Metal" MaterialPointSolveroneAPIExt = "oneAPI" +MaterialPointSolverDebugExt = "WGLMakie" [extras] AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" @@ -46,6 +53,8 @@ oneAPI = "8f75cd03-7ff8-4ecb-9b8f-daf728133b1b" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +WGLMakie = "276b4fcb-3e11-5398-bf8b-a0c2d153d008" [targets] test = ["Test", "BenchmarkTools", "DelimitedFiles"] \ No newline at end of file diff --git a/README.md b/README.md index 94fac8b..2657286 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@ [![CI](https://github.com/LandslideSIM/MaterialPointSolver.jl/actions/workflows/ci.yml/badge.svg)](https://github.com/LandslideSIM/MaterialPointSolver.jl/actions/workflows/ci.yml) [![Stable](https://img.shields.io/badge/docs-stable-blue.svg?logo=quicklook)](https://landslidesim.github.io/MaterialPointSolver.jl/stable/) [![Dev](https://img.shields.io/badge/docs-dev-red.svg?logo=quicklook)](https://landslidesim.github.io/MaterialPointSolver.jl/dev/) -[![Version](https://img.shields.io/badge/version-v0.4.5-pink)]() +[![Version](https://img.shields.io/badge/version-v0.4.6-pink)]() [![](https://img.shields.io/badge/NVIDIA-CUDA-green.svg?logo=nvidia)](https://developer.nvidia.com/cuda-toolkit) [![](https://img.shields.io/badge/AMD-ROCm-red.svg?logo=amd)](https://www.amd.com/en/products/software/rocm.html) @@ -14,7 +14,9 @@ This package provides a high-performance, backend-agnostic implementation of the Material Point Method (MPM) using the Julia Language. It is lightweight and user-friendly, allowing efficient execution on various hardware accelerators with a single codebase. Please check here for the documentation.

- +> If you have a GPU from Intel® and want to try on it, please [contact me 📧](mailto:zenan.huo@outlook.com). + + ## Installation ⚙️ @@ -65,31 +67,44 @@ julia> ] - ✅ one-click switch between `FP64` and `FP32` - ✅ user-defined algorithms/extensions at any level +There is a `debug` model can be used to make sure the simulation is working as expect. It's also can be used for in-situ visualization in VSCode. + + + ## Citation 🔥 If you find `MaterialPointSolver.jl` useful or have used it in your research, please cite it as follows: ```bib -@article{index, - title={Here is the title}, - author={authors}, - journal={journal}, - year={year} +@article{HUO2025107189, +title = {A high-performance backend-agnostic Material Point Method solver in Julia}, +journal = {Computers and Geotechnics}, +volume = {183}, +pages = {107189}, +year = {2025}, +issn = {0266-352X}, +doi = {https://doi.org/10.1016/j.compgeo.2025.107189}, +url = {https://www.sciencedirect.com/science/article/pii/S0266352X25001387}, +author = {Zenan Huo and Yury Alkhimenkov and Michel Jaboyedoff and Yury Podladchikov and Ludovic Räss and Emmanuel Wyser and Gang Mei}, +keywords = {MPM, Julia language, Heterogeneous computing, Effective memory throughput} } ``` > [!CAUTION] > This is the latest version of `MaterialPointSover.jl`, if you want to see the examples in the paper, please move to [https://github.com/LandslideSIM/Archive_MaterialPointSolver.jl_paper](https://github.com/LandslideSIM/Archive_MaterialPointSolver.jl_paper). +> [!TIP] +> After the article was published, we released many new features and achieved significant performance improvements. We are currently actively working on completing the documentation. If possible, you can directly review the source code. + ## Acknowledgement 👍 This project is sponsored by [Risk Group | Université de Lausanne](https://wp.unil.ch/risk/) and [China Scholarship Council [中国国家留学基金管理委员会]](https://www.csc.edu.cn/). ## MPM ➕ Julia -* [[package]: elastoPlasm.jl](https://github.com/ewyser/elastoPlasm.jl) is fully witten in Julia, it solves explicit elasto-plastic problems within a finite deformation framework. +* [📦 [package]: elastoPlasm.jl](https://github.com/ewyser/elastoPlasm.jl) is fully witten in Julia, it solves explicit elasto-plastic problems within a finite deformation framework. -* [[package]: Tesserae.jl](https://github.com/KeitaNakamura/Tesserae.jl) is an MPM-related Julia package, it provides some useful functions that can be used for MPM, such as convenient macros for transferring data between grids and particles. +* [📦 [package]: Tesserae.jl](https://github.com/KeitaNakamura/Tesserae.jl) is an MPM-related Julia package, it provides some useful functions that can be used for MPM, such as convenient macros for transferring data between grids and particles. -* [[code]: MPM-Julia](https://github.com/vinhphunguyen/MPM-Julia) is the code for the paper: Sinai, V.P. Nguyen, C.T. Nguyen and S. Bordas. Programming the Material Point Method in Julia. Advances in Engineering Software,105: 17--29, 2017. +* [✍️ [code]: MPM-Julia](https://github.com/vinhphunguyen/MPM-Julia) is the code for the paper: Sinai, V.P. Nguyen, C.T. Nguyen and S. Bordas. Programming the Material Point Method in Julia. Advances in Engineering Software,105: 17--29, 2017. -* [[code]: jump](https://github.com/vinhphunguyen/jump) is for the theory of the MPM described in the book 'The Material Point Method: Theory, Implementations and Applications (Scientific Computation) 1st ed. 2023 Edition'. [https://link.springer.com/book/10.1007/978-3-031-24070-6](https://link.springer.com/book/10.1007/978-3-031-24070-6) +* [[✍️ code]: jump](https://github.com/vinhphunguyen/jump) is for the theory of the MPM described in the book 'The Material Point Method: Theory, Implementations and Applications (Scientific Computation) 1st ed. 2023 Edition'. [https://link.springer.com/book/10.1007/978-3-031-24070-6](https://link.springer.com/book/10.1007/978-3-031-24070-6) diff --git a/docs/Project.toml b/docs/Project.toml index 3ae66de..ab3476d 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -2,10 +2,12 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365" MaterialPointSolver = "dc4aa359-7d89-437d-91bb-f4b330c522c1" +WGLMakie = "276b4fcb-3e11-5398-bf8b-a0c2d153d008" [sources] MaterialPointSolver = {path = ".."} [compat] Documenter = "1.8" -DocumenterVitepress = "0.1" \ No newline at end of file +DocumenterVitepress = "0.1" +WGLMakie = "0.11" \ No newline at end of file diff --git a/docs/assets/debug.gif b/docs/assets/debug.gif new file mode 100644 index 0000000..bf15a03 Binary files /dev/null and b/docs/assets/debug.gif differ diff --git a/docs/make.jl b/docs/make.jl index 2121f7a..a8e3e1f 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,4 +1,4 @@ -using Documenter, DocumenterVitepress, MaterialPointSolver +using Documenter, DocumenterVitepress, MaterialPointSolver, WGLMakie makedocs( modules = [MaterialPointSolver], @@ -23,8 +23,7 @@ makedocs( "interface/kernelfunctions.md", "interface/mpmprocess.md", "interface/workflow.md", - "interface/import.md", - "interface/export.md" + "interface/imexport.md", ], "Tutorials" => Any[ "tutorials/example1.md", @@ -38,8 +37,9 @@ makedocs( "advancedtopics/customextension.md", "advancedtopics/customfriction.md" ], - "utils.md", - "showcases.md" + "Useful Tools" => Any[ + "utils/debug.md", + ] ], warnonly = [:missing_docs, :cross_references], ) diff --git a/docs/src/interface/export.md b/docs/src/interface/export.md deleted file mode 100644 index 9c00155..0000000 --- a/docs/src/interface/export.md +++ /dev/null @@ -1 +0,0 @@ -# Export simulation results \ No newline at end of file diff --git a/docs/src/interface/imexport.md b/docs/src/interface/imexport.md new file mode 100644 index 0000000..9e6244a --- /dev/null +++ b/docs/src/interface/imexport.md @@ -0,0 +1,21 @@ +# Import/Export Simulation Results + +We use Julia's built-in [`Serialization`](https://docs.julialang.org/en/v1/stdlib/Serialization/#Serialization) functionality to implement the import and export of models. Sometimes the model file size is too large, so we use CodecZlib for compression (optional), and automatically determine whether it is a compressed file when reading. + + +!!! note + + We use the file suffix `.mpm` to denote model files generated by ***MaterialPointSolver.jl***. + +```@docs +exportmodel( + args :: DeviceArgs{T1, T2}, + grid :: DeviceGrid{T1, T2}, + mp :: DeviceParticle{T1, T2}, + attr :: DeviceProperty{T1, T2}, + bc ::DeviceVBoundary{T1, T2}; + model_file::String = "temp.mpm", + compress ::Bool = true +) where {T1, T2} +importmodel(mpm_file::String) +``` \ No newline at end of file diff --git a/docs/src/interface/import.md b/docs/src/interface/import.md deleted file mode 100644 index 61d9424..0000000 --- a/docs/src/interface/import.md +++ /dev/null @@ -1 +0,0 @@ -# Import computing models \ No newline at end of file diff --git a/docs/src/utils.md b/docs/src/utils.md deleted file mode 100644 index a77562a..0000000 --- a/docs/src/utils.md +++ /dev/null @@ -1 +0,0 @@ -# Useful Tools \ No newline at end of file diff --git a/docs/src/utils/debug.md b/docs/src/utils/debug.md new file mode 100644 index 0000000..1c96129 --- /dev/null +++ b/docs/src/utils/debug.md @@ -0,0 +1,113 @@ +# Debug Mode (in-situ visualization) + +We have provided a mode that can display simulation results in real-time, which we call "debug mode." It is also backend-agnostic and happens to enable in-situ visualization. The original intent of this design was to allow users to immediately view the computation results while running a simulation. If the results do not develop as expected, the computation can be terminated at any time for a quick code review. + +Note that this feature will only be enabled after the user explicitly calls WGLMakie: + +```@example debug +using MaterialPointSolver +using WGLMakie +``` + +!!! note + + - We implement this functionality through WGLMakie.jl, which theoretically supports computations on a remote headless server. However, as WGLMakie.jl is still in an experimental stage, we recommend launching Debug mode only on local machine for now. + - We prefer to always use the plot panel provided by the Julia Language plugin in VSCode to display results; please enable the plot panel in the julia-vscode plugin. + - The current maximum number of particles displayed is 500,000, but this can be customized by the user. + +## API + +In normal mode, we assume that you have already instantiated and obtained `args`, `grid`, `mp`, `attr`, and `bc`. + +:::tabs + +== normal mode + +```julia +materialpointsolver!(args, grid, mp, attr, bc) +``` + +== debug mode + +```julia +materialpointsolver!(args, grid, mp, attr, bc, debug) +``` + +::: + +where the type of `debug` is `DebugConfig`: + +```@docs +DebugConfig +``` + +The current `DebugPlot` type has only one field, `plot::DebugPlot`, but more features may be added in the future. + +The type of `plot` is `DebugPlot`: + +```@docs +DebugPlot +``` + +Here, we may need to explain the meaning of the calculate field in `DebugPlot`. Sometimes, the values we want to display might be derived from multiple attributes, such as displacement, so we need to define the computation rules. + +!!! note + + - The return type of the user-defined `calculate` function must currently be a vector with the same length as the number of material points, and the input parameters of this custom function should always be: (i) `grid`, (ii) `mp`, (iii) `attr`, (iv) `bc`, and (v) vid, with the order fixed and unchangeable. + - Users can directly use `debug=DebugConfig()`, which will apply the default debug configuration, visualizing the color based on the density of the material points. + + +## Examples + +If you want to visualize the equivalent plastic strain, then： + +```@example debug +plotfunc(grid, mp, attr, bc, vid) = Array{Float32, 1}(mp.ϵq) +plotconfig = DebugPlot(cbname="equivalent plastic strain", calculate=plotfunc) +debug = DebugConfig(plot=plotconfig) +``` + +```@raw html + +``` + +If you want to visualize the displacement： + +```@example debug +function tmpf(grid, mp, attr, bc, vid) + ξ0 = Array{Float32, 2}(mp.ξ0[vid, :]) + ξ1 = Array{Float32, 2}(mp.ξ[vid, :] ) + + val = @. sqrt((ξ1[:, 1] - ξ0[:, 1])^2 .+ + (ξ1[:, 2] - ξ0[:, 2])^2 .+ + (ξ1[:, 3] - ξ0[:, 3])^2) + return val +end +plotconfig = DebugPlot(pointsize=1, axis=true, axfontsize=10, calculate=tmpf, cbname="disp") +debug = DebugConfig(plotconfig) +``` + +```@raw html + +``` + +If you need a smooth visualization, try to decrease the `pointsize` or `pointnum`: + + +```@example debug +function tmpf(grid, mp, attr, bc, vid) + ξ0 = Array{Float32, 2}(mp.ξ0[vid, :]) + ξ1 = Array{Float32, 2}(mp.ξ[vid, :] ) + + val = @. sqrt((ξ1[:, 1] - ξ0[:, 1])^2 .+ + (ξ1[:, 2] - ξ0[:, 2])^2 .+ + (ξ1[:, 3] - ξ0[:, 3])^2) + return val +end +plotconfig = DebugPlot(pointsize=1f-2, axis=true, axfontsize=10, calculate=tmpf, cbname="disp", pointnum=1000) # [!code focus] +debug = DebugConfig(plotconfig) +``` + +```@raw html + +``` \ No newline at end of file diff --git a/docs/src/utils/figures/debug1.png b/docs/src/utils/figures/debug1.png new file mode 100644 index 0000000..e3782eb Binary files /dev/null and b/docs/src/utils/figures/debug1.png differ diff --git a/docs/src/utils/figures/debug2.png b/docs/src/utils/figures/debug2.png new file mode 100644 index 0000000..ce5d0ec Binary files /dev/null and b/docs/src/utils/figures/debug2.png differ diff --git a/docs/src/utils/figures/debug3.png b/docs/src/utils/figures/debug3.png new file mode 100644 index 0000000..7f0a016 Binary files /dev/null and b/docs/src/utils/figures/debug3.png differ diff --git a/ext/CUDAExt/devicehelpfunc_cuda.jl b/ext/CUDAExt/devicehelpfunc_cuda.jl index 8894185..09ac4cc 100644 --- a/ext/CUDAExt/devicehelpfunc_cuda.jl +++ b/ext/CUDAExt/devicehelpfunc_cuda.jl @@ -40,8 +40,7 @@ function device2host!( args :: DeviceArgs{T1, T2}, mp ::DeviceParticle{T1, T2}, dev_mp ::DeviceParticle{T1, T2}, - ::Val{:CUDA}; - verbose::Bool=false + ::Val{:CUDA} ) where {T1, T2} copyto!(mp.σm , dev_mp.σm ) copyto!(mp.vs , dev_mp.vs ) @@ -60,11 +59,53 @@ function device2host!( copyto!(mp.ϵijw, dev_mp.ϵijw); copyto!(mp.n , dev_mp.n ); ) : nothing - if verbose == true + + return nothing +end + +function device2host!( + grid :: DeviceGrid{T1, T2}, + mp :: DeviceParticle{T1, T2}, + attr :: DeviceProperty{T1, T2}, + bc ::DeviceVBoundary{T1, T2}, + dev_grid:: DeviceGrid{T1, T2}, + dev_mp :: DeviceParticle{T1, T2}, + dev_attr:: DeviceProperty{T1, T2}, + dev_bc ::DeviceVBoundary{T1, T2}, + ::Val{:CUDA}; + verbose ::Bool=true +) where {T1, T2} + # download grid data from device + @inbounds for field in fieldnames(typeof(grid)) + if getfield(grid, field) isa AbstractArray + copyto!(getfield(grid, field), getfield(dev_grid, field)) + end + end + # download mp data from device + @inbounds for field in fieldnames(typeof(mp)) + if getfield(mp, field) isa AbstractArray + copyto!(getfield(mp, field), getfield(dev_mp, field)) + end + end + # download attr data from device + @inbounds for field in fieldnames(typeof(attr)) + if getfield(attr, field) isa AbstractArray + copyto!(getfield(attr, field), getfield(dev_attr, field)) + end + end + # download bc data from device + @inbounds for field in fieldnames(typeof(bc)) + if getfield(bc, field) isa AbstractArray + copyto!(getfield(bc, field), getfield(dev_bc, field)) + end + end + + if verbose dev_id = CUDA.device().handle content = "downloading from :CUDA [$(dev_id)] → host" println("\e[1;31m[▼ I/O:\e[0m \e[0;31m$(content)\e[0m") end + return nothing end diff --git a/ext/DebugExt/plotutils.jl b/ext/DebugExt/plotutils.jl new file mode 100644 index 0000000..01828b0 --- /dev/null +++ b/ext/DebugExt/plotutils.jl @@ -0,0 +1,47 @@ +function gettheme() + custom_dark = theme_dark() + custom_dark.textcolor = :white + custom_dark.linecolor = :white + custom_dark.backgroundcolor = "#181818" + return custom_dark +end + +function getbox(grid::DeviceGrid2D{T1, T2}) where {T1, T2} + vmin = minimum(grid.ξ, dims=1) + vmax = maximum(grid.ξ, dims=1) + + # 直接创建最终的8个Point2f元组数组（4条边 × 2个点） + return Point2f[ + (vmin[1], vmin[2]), (vmax[1], vmin[2]), # 底边 + (vmax[1], vmin[2]), (vmax[1], vmax[2]), # 右边 + (vmax[1], vmax[2]), (vmin[1], vmax[2]), # 顶边 + (vmin[1], vmax[2]), (vmin[1], vmin[2]) # 左边 + ] +end + +function getbox(grid::DeviceGrid3D{T1, T2}) where {T1, T2} + vmin = minimum(grid.ξ, dims=1) + vmax = maximum(grid.ξ, dims=1) + + # 直接创建24个Point3f元组数组（12条边 × 2个点） + return Point3f[ + # 底面 (4条边) + (vmin[1], vmin[2], vmin[3]), (vmax[1], vmin[2], vmin[3]), + (vmax[1], vmin[2], vmin[3]), (vmax[1], vmax[2], vmin[3]), + (vmax[1], vmax[2], vmin[3]), (vmin[1], vmax[2], vmin[3]), + (vmin[1], vmax[2], vmin[3]), (vmin[1], vmin[2], vmin[3]), + + # 顶面 (4条边) + (vmin[1], vmin[2], vmax[3]), (vmax[1], vmin[2], vmax[3]), + (vmax[1], vmin[2], vmax[3]), (vmax[1], vmax[2], vmax[3]), + (vmax[1], vmax[2], vmax[3]), (vmin[1], vmax[2], vmax[3]), + (vmin[1], vmax[2], vmax[3]), (vmin[1], vmin[2], vmax[3]), + + # 侧边 (4条边) + (vmin[1], vmin[2], vmin[3]), (vmin[1], vmin[2], vmax[3]), + (vmax[1], vmin[2], vmin[3]), (vmax[1], vmin[2], vmax[3]), + (vmax[1], vmax[2], vmin[3]), (vmax[1], vmax[2], vmax[3]), + (vmin[1], vmax[2], vmin[3]), (vmin[1], vmax[2], vmax[3]) + ] +end + diff --git a/ext/MaterialPointSolverDebugExt.jl b/ext/MaterialPointSolverDebugExt.jl new file mode 100644 index 0000000..d4fa61d --- /dev/null +++ b/ext/MaterialPointSolverDebugExt.jl @@ -0,0 +1,123 @@ +module MaterialPointSolverDebugExt + +using MaterialPointSolver +using KernelAbstractions +using WGLMakie + +import MaterialPointSolver: materialpointsolver!, info_print, submit_work!, perf, procedure! +import MaterialPointSolver: DeviceArgs, DeviceGrid, DeviceParticle, DeviceProperty, + DeviceVBoundary, DebugConfig +import MaterialPointSolver: initmpstatus!, progressinfo, host2device, updatepb!, KAsync, + getBackend, clean_device! +import StatsBase: mean, sample + +include(joinpath(@__DIR__, "DebugExt/plotutils.jl")) + +""" + materialpointsolver!(args::DeviceArgs, grid::DeviceGrid, mp::DeviceParticle, + attr::DeviceProperty, bc::DeviceVBoundary, debug::DebugConfig; + workflow::F=procedure!) + +Description: +--- +This function is the main entry point for the material point method (MPM) solver. Here we +added `debug` to enable an in-situ visulaization of the simulation. This can help to check +the problem of the code. Note that it is just used for debug not for production mode. +""" +function materialpointsolver!( + args :: DeviceArgs{T1, T2}, + grid :: DeviceGrid{T1, T2}, + mp :: DeviceParticle{T1, T2}, + attr :: DeviceProperty{T1, T2}, + bc ::DeviceVBoundary{T1, T2}, + debug ::DebugConfig; + workflow::F=procedure! +) where {T1, T2, F<:Function} + info_print(args, grid, mp) # terminal info + submit_work!(args, grid, mp, attr, bc, workflow, debug) # MPM solver + perf(args) # performance summary + return nothing +end + +@views function submit_work!( + args :: DeviceArgs{T1, T2}, + grid :: DeviceGrid{T1, T2}, + mp :: DeviceParticle{T1, T2}, + attr :: DeviceProperty{T1, T2}, + bc ::DeviceVBoundary{T1, T2}, + workflow::F, + debug ::DebugConfig +) where {T1, T2, F<:Function} + initmpstatus!(CPU())(ndrange=mp.np, grid, mp, Val(args.basis)) + # variables setup for the simulation + Ti = T2(0.0) + pc = Ref{T1}(0) + pb = progressinfo(args, "solving") + ΔT = args.ΔT + #ΔT = args.time_step==:auto ? cfl(args, grid, mp, attr, Val(args.coupling)) : args.ΔT + dev_grid, dev_mp, dev_attr, dev_bc = host2device(grid, mp, attr, bc, Val(args.device)) + + # plot setup + vid = mp.np > debug.plot.pointnum ? sample(1:mp.np, debug.plot.pointnum, replace=false) : (1:mp.np) + getvattr = debug.plot.calculate(grid, mp, attr, bc, vid) + pdims = size(mp.ξ, 2) + vtime = Observable(Ti) + vpts = Observable(Array{Float32, 2}(mp.ξ[vid, :])) + vattr = Observable(Array{Float32, 1}(getvattr)) + + set_theme!(gettheme()) + fig = Figure() + ax = LScene(fig[1, 1], show_axis=debug.plot.axis) + if debug.plot.axis + canvas = ax.scene.plots[1] + canvas.ticks[:textcolor] = :white + canvas.ticks[:fontsize ] = debug.plot.axfontsize + canvas.frame[:axiscolor] = "#818181" + canvas.names[:textcolor] = :white + canvas.names[:fontsize ] = debug.plot.axfontsize + end + Label(fig[1, 1, Top()], @lift("Simulation Time: $(lpad(round($vtime, digits=2), 8)) s"), + padding = (0, 0, 5, 0)) + Label(fig[1, 1, Bottom()], "Debug mode: ON") + linesegments!(ax, getbox(grid), linewidth=0) + plt = scatter!(ax, vpts, color=vattr, markersize=debug.plot.pointsize, + colormap=debug.plot.colormap, transparency=false, depthsorting=false, marker=:rect) + if debug.plot.colorbar + Colorbar(fig[1, 2], plt, label=debug.plot.cbname, spinewidth=0, + tickformat=debug.plot.tickformat) + end + pdims==2 && cam2d!(ax.scene) + display(fig) + + # main part: debug ON + debug_id = T1(1) # debug index + debug_switch = T1(0) # debug step + args.start_time = time() + while Ti < args.Ttol + workflow(args, dev_grid, dev_mp, dev_attr, dev_bc, ΔT, Ti, + Val(args.coupling), Val(args.scheme)) + ΔT = args.time_step == :auto ? args.αT * reduce(min, dev_mp.cfl) : args.ΔT + Ti += ΔT + debug_switch += 1 + args.iter_num += 1 + updatepb!(pc, Ti, args.Ttol, pb) + if (debug_switch == args.hdf5_step) || (debug_switch == T1(0)) + copyto!(vpts[], dev_mp.ξ[vid, :]) + copyto!(vattr[], debug.plot.calculate(dev_grid, dev_mp, dev_attr, dev_bc, vid)) + notify(vpts) + notify(vattr) + vtime[] = Ti + yield() + debug_switch = T1(0) + debug_id += T1(1) + end + end + args.end_time = time() + + KAsync(getBackend(Val(args.device))) + device2host!(grid, mp, attr, bc, dev_grid, dev_mp, dev_attr, dev_bc, Val(args.device)) + clean_device!(dev_grid, dev_mp, dev_attr, dev_bc, Val(args.device)) + return nothing +end + +end \ No newline at end of file diff --git a/ext/MetalExt/devicehelpfunc_metal.jl b/ext/MetalExt/devicehelpfunc_metal.jl index 5d92b96..c553f6b 100644 --- a/ext/MetalExt/devicehelpfunc_metal.jl +++ b/ext/MetalExt/devicehelpfunc_metal.jl @@ -11,6 +11,7 @@ | 3. clean_gpu! [2D & 3D] | | 4. Tpeak | | 5. getBackend | +| 6. getArray | +==========================================================================================# function host2device( @@ -38,8 +39,7 @@ function device2host!( args :: DeviceArgs{T1, T2}, mp ::DeviceParticle{T1, T2}, dev_mp ::DeviceParticle{T1, T2}, - ::Val{:Metal}; - verbose::Bool=false + ::Val{:Metal} ) where {T1, T2} copyto!(mp.σm , dev_mp.σm ) copyto!(mp.vs , dev_mp.vs ) @@ -58,10 +58,52 @@ function device2host!( copyto!(mp.ϵijw, dev_mp.ϵijw); copyto!(mp.n , dev_mp.n ); ) : nothing - if verbose == true + + return nothing +end + +function device2host!( + grid :: DeviceGrid{T1, T2}, + mp :: DeviceParticle{T1, T2}, + attr :: DeviceProperty{T1, T2}, + bc ::DeviceVBoundary{T1, T2}, + dev_grid:: DeviceGrid{T1, T2}, + dev_mp :: DeviceParticle{T1, T2}, + dev_attr:: DeviceProperty{T1, T2}, + dev_bc ::DeviceVBoundary{T1, T2}, + ::Val{:Metal}; + verbose ::Bool=true +) where {T1, T2} + # download grid data from device + @inbounds for field in fieldnames(typeof(grid)) + if getfield(grid, field) isa AbstractArray + copyto!(getfield(grid, field), getfield(dev_grid, field)) + end + end + # download mp data from device + @inbounds for field in fieldnames(typeof(mp)) + if getfield(mp, field) isa AbstractArray + copyto!(getfield(mp, field), getfield(dev_mp, field)) + end + end + # download attr data from device + @inbounds for field in fieldnames(typeof(attr)) + if getfield(attr, field) isa AbstractArray + copyto!(getfield(attr, field), getfield(dev_attr, field)) + end + end + # download bc data from device + @inbounds for field in fieldnames(typeof(bc)) + if getfield(bc, field) isa AbstractArray + copyto!(getfield(bc, field), getfield(dev_bc, field)) + end + end + + if verbose content = "downloading from :Metal [1] → host" println("\e[1;31m[▼ I/O:\e[0m \e[0;31m$(content)\e[0m") end + return nothing end diff --git a/src/MaterialPointSolver.jl b/src/MaterialPointSolver.jl index 7f74b47..3fc20e2 100644 --- a/src/MaterialPointSolver.jl +++ b/src/MaterialPointSolver.jl @@ -13,6 +13,8 @@ module MaterialPointSolver using Adapt, BenchmarkTools, Dates, DelimitedFiles, HDF5, KernelAbstractions, Printf, ProgressMeter, SysInfo# PrecompileTools +import Serialization: serialize, deserialize +import CodecZlib: ZlibCompressorStream, ZlibDecompressorStream import KernelAbstractions.synchronize as KAsync import KernelAbstractions.Extras: @unroll as @KAunroll import Suppressor.@suppress as @MPSsuppress diff --git a/src/solver.jl b/src/solver.jl index d5fec02..5861e54 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -134,7 +134,7 @@ This function will start to run the MPM solver. args.end_time = time() end KAsync(getBackend(Val(args.device))) - device2host!(args, mp, dev_mp, Val(args.device), verbose=true) + device2host!(grid, mp, attr, bc, dev_grid, dev_mp, dev_attr, dev_bc, Val(args.device)) clean_device!(dev_grid, dev_mp, dev_attr, dev_bc, Val(args.device)) return nothing end \ No newline at end of file diff --git a/src/toolkits/devicehelpfunc.jl b/src/toolkits/devicehelpfunc.jl index 54a8b46..24b76a1 100644 --- a/src/toolkits/devicehelpfunc.jl +++ b/src/toolkits/devicehelpfunc.jl @@ -6,16 +6,18 @@ | Programmer : Zenan Huo | | Start Date : 01/01/2022 | | Affiliation: Risk Group, UNIL-ISTE | -| Functions : 1. host2device [2D & 3D] | -| 2. device2host! [2D & 3D] | -| 3. clean_gpu! [2D & 3D] | -| 4. memcpy! | -| 5. Tpeak | -| 6. getBackend | -| 7. getArray | +| Functions : 1. host2device [2D & 3D] | +| 2. device2host_p! [2D & 3D] | +| 3. device2host! [2D & 3D] | +| 4. clean_gpu! [2D & 3D] | +| 5. memcpy! | +| 6. Tpeak | +| 7. getBackend | +| 8. getArray | +==========================================================================================# -export host2device, device2host!, clean_device!, memcpy!, Tpeak, getBackend, getArray +export host2device, device2host_p!, device2host!, clean_device!, memcpy!, Tpeak, getBackend, + getArray """ host2device(grid::DeviceGrid{T1, T2}, mp::DeviceParticle{T1, T2}, @@ -37,18 +39,43 @@ end """ device2host!(args::DeviceArgs{T1, T2}, mp::DeviceParticle{T1, T2}, - dev_mp::GPUDeviceParticle{T1, T2}, ::Val{:CPU/:CUDA/...}; verbose::Bool=false) + dev_mp::GPUDeviceParticle{T1, T2}, ::Val{:CPU}) Description: --- -Transfer data from device to host. (2D & 3D) +Transfer data partialy from device to host. (2D & 3D) """ function device2host!( args :: DeviceArgs{T1, T2}, mp ::DeviceParticle{T1, T2}, dev_mp ::DeviceParticle{T1, T2}, - ::Val{:CPU}; - verbose::Bool=false + ::Val{:CPU} +) where {T1, T2} + return nothing +end + +""" + device2host!(grid::DeviceArgs{T1, T2}, mp::DeviceParticle{T1, T2}, + attr::DeviceProperty{T1, T2}, bc::DeviceVBoundary{T1, T2}, + dev_grid::GPUDeviceGrid{T1, T2}, dev_mp::GPUDeviceParticle{T1, T2}, + dev_attr::GPUDeviceParticle{T1, T2}, dev_bc::GPUDeviceVBoundary{T1, T2}, + ::Val{:CPU}; verbose::Bool=false) + +Description: +--- +Transfer full data from device to host. (2D & 3D) +""" +function device2host!( + grid :: DeviceGrid{T1, T2}, + mp :: DeviceParticle{T1, T2}, + attr :: DeviceProperty{T1, T2}, + bc ::DeviceVBoundary{T1, T2}, + dev_grid:: DeviceGrid{T1, T2}, + dev_mp :: DeviceParticle{T1, T2}, + dev_attr:: DeviceProperty{T1, T2}, + dev_bc ::DeviceVBoundary{T1, T2}, + ::Val{:CPU}; + verbose ::Bool=true ) where {T1, T2} return nothing end diff --git a/src/toolkits/modelinfo.jl b/src/toolkits/modelinfo.jl index 07b5dc4..c753051 100644 --- a/src/toolkits/modelinfo.jl +++ b/src/toolkits/modelinfo.jl @@ -8,9 +8,14 @@ | Affiliation: Risk Group, UNIL-ISTE | | Functions : 01. check_datasize | | 02. @memcheck | +| 03. exportmodel | +| 04. importmodel | +==========================================================================================# -export check_datasize, @memcheck +export check_datasize +export @memcheck +export exportmodel +export importmodel function check_datasize(args:: DeviceArgs{T1, T2}, grid:: DeviceGrid{T1, T2}, @@ -57,3 +62,105 @@ macro memcheck(expr) @info "💾 $data_str" # print info end end + +""" + exportmodel(args, grid, mp, attr, bc; model_file="temp.mpm", compress=true) + +Description: +--- +Export model data to a compressed or uncompressed file (.mpm). If user does not identify the +file path, the `.mpm` file will save in the project directory with the project name. +""" +function exportmodel( + args :: DeviceArgs{T1, T2}, + grid :: DeviceGrid{T1, T2}, + mp :: DeviceParticle{T1, T2}, + attr :: DeviceProperty{T1, T2}, + bc ::DeviceVBoundary{T1, T2}; + model_file::String = "temp.mpm", + compress ::Bool = true +) where {T1, T2} + # get model file path and timestamp + model_file == "temp.mpm" && (model_file = joinpath( + args.project_path, args.project_name, args.project_name * ".mpm")) + ctime = now() + timestamp = Int64(floor(datetime2unix(ctime))) + + # print info for the user + model_status = compress ? "compressed model" : "uncompressed model" + println("\e[1;33m[ Export:\e[0m", " $(model_status) data (.mpm)") + println(" ", "\e[1;33m⦿\e[0m", " time stamp: $(ctime)") + + # save model data + # 1) compress header, 2) timestamp, 3) args, 4) grid, 5) mp, 6) attr, 7) bc + open(model_file, "w") do io + # compress header（0x01 means compressed, 0x00 means uncompressed） + compress ? write(io, UInt8(0x01)) : write(io, UInt8(0x00)) + # write timestamp + write(io, timestamp) + cstream = compress ? ZlibCompressorStream(io) : io + serialize(cstream, args) + println(" ", "\e[1;33m⦿\e[0m", " args: ", "\e[0;33mdone\e[0m") + serialize(cstream, grid) + println(" ", "\e[1;33m⦿\e[0m", " grid: ", "\e[0;33mdone\e[0m") + serialize(cstream, mp ) + println(" ", "\e[1;33m⦿\e[0m", " mp: ", "\e[0;33mdone\e[0m") + serialize(cstream, attr) + println(" ", "\e[1;33m⦿\e[0m", " attr: ", "\e[0;33mdone\e[0m") + serialize(cstream, bc ) + println(" ", "\e[1;33m⦿\e[0m", " bc: ", "\e[0;33mdone\e[0m") + close(cstream) + end + + # print info + modelsize = round(filesize(model_file) / 1024^3, digits=2) + println(" ", "\e[1;33m⦿\e[0m", " file size: ", "$(modelsize) GiB") + println(" ", "\e[1;33m⦿\e[0m", " file path: ", model_file) + return nothing +end + +""" + importmodel(mpm_file::String) + +Description: +--- +Import model data from a compressed or uncompressed file (.mpm). This function will detect +wheather the model file is compressed or not, and then read the model data. +""" +function importmodel(mpm_file::String) + # input file check + isfile(mpm_file) || throw(ArgumentError("file not found: $mpm_file")) + splitext(mpm_file)[end] == ".mpm" || throw(ArgumentError( + "invalid file format: $(splitext(mpm_file)[end])")) + + # read model data + args, grid, mp, attr, bc = open(mpm_file, "r") do io + # get comrression header + compression_flag = read(io, UInt8) + compress = compression_flag == 0x01 ? true : false + # print info for the user + model_status = compress ? "compressed model" : "uncompressed model" + println("\e[1;33m[ Import:\e[0m", " $(model_status) data (.mpm)") + # get the timestamp + timestamp = read(io, Int64) + creation_time = unix2datetime(Float64(timestamp)) + println(" ", "\e[1;33m⦿\e[0m", " time stamp: $(creation_time)") + # read model data + cstream = compress ? ZlibDecompressorStream(io) : io + args = deserialize(cstream) + println(" ", "\e[1;33m⦿\e[0m", " args: ", "\e[0;33mdone\e[0m") + grid = deserialize(cstream) + println(" ", "\e[1;33m⦿\e[0m", " grid: ", "\e[0;33mdone\e[0m") + mp = deserialize(cstream) + println(" ", "\e[1;33m⦿\e[0m", " mp: ", "\e[0;33mdone\e[0m") + attr = deserialize(cstream) + println(" ", "\e[1;33m⦿\e[0m", " attr: ", "\e[0;33mdone\e[0m") + bc = deserialize(cstream) + println(" ", "\e[1;33m⦿\e[0m", " bc: ", "\e[0;33mdone\e[0m") + close(cstream) + println() + args, grid, mp, attr, bc + end + + return args, grid, mp, attr, bc +end \ No newline at end of file diff --git a/src/type.jl b/src/type.jl index aa9bc39..1054080 100644 --- a/src/type.jl +++ b/src/type.jl @@ -12,4 +12,5 @@ include(joinpath(@__DIR__, "types/modelargs.jl")) include(joinpath(@__DIR__, "types/grid.jl" )) include(joinpath(@__DIR__, "types/particle.jl" )) include(joinpath(@__DIR__, "types/property.jl" )) -include(joinpath(@__DIR__, "types/boundary.jl" )) \ No newline at end of file +include(joinpath(@__DIR__, "types/boundary.jl" )) +include(joinpath(@__DIR__, "types/debug.jl" )) \ No newline at end of file diff --git a/src/types/debug.jl b/src/types/debug.jl new file mode 100644 index 0000000..1b2f9ce --- /dev/null +++ b/src/types/debug.jl @@ -0,0 +1,69 @@ +#==========================================================================================+ +| MaterialPointSolver.jl: High-performance MPM Solver for Geomechanics | ++------------------------------------------------------------------------------------------+ +| File Name : debug.jl | +| Description: Type system for debug configurations in MaterialPointSolver.jl | +| Programmer : Zenan Huo | +| Start Date : 01/01/2022 | +| Affiliation: Risk Group, UNIL-ISTE | +| License : MIT License | ++==========================================================================================# + +export DebugPlot +export DebugConfig + +def_func(grid, mp, attr, bc, vid) = Array{Float32}(mp.ρs) + +""" + DebugPlot + +Description: +--- +This struct is used to store the WGLMakie configurations for the in-situ visualization of +the simulation. + +Users need to provice `attr` and `cbname` to specify the attribute to plot and the name of +the colorbar, even you choose `axis = false`. + +- `axis::Bool = false` : Whether to show the axis +- `colorbar::Bool = true` : Whether to show the colorbar +- `colormap::Symbol = :readblue` : The colormap for the plot +- `pointsize::Real = 1f-2` : The size of the points +- `pointnum::Int = 500000` : The number of points to plot +- `tickformat::String = "{:9.2f}"` : The format of the ticks +- `calculate::Function = def_func` : The function to calculate the value to plot +- `axfontsize::Real = 24` : The font size of the axis +- `cbfontsize::Real = 24` : The font size of the colorbar +- `cbname::String` : The name of the colorbar +""" +@kwdef struct DebugPlot + axis ::Bool = false + colorbar ::Bool = true + colormap ::Symbol = :redblue + pointsize ::Real = 1f-2 + pointnum ::Int = 500000 + tickformat::String = "{:9.2f}" + calculate ::Function = def_func + axfontsize::Real = 24 + cbfontsize::Real = 24 + cbname ::String +end + + +def_debug = DebugPlot(cbname="density") + +""" + DebugConfig + +Description: +--- +This struct is used to store the debug configurations for the in-situ visualization of the +simulation. + +- `plot`::DebugPlot : The WGLMakie configurations for the debug mode +- ... : More debug configurations can be added in the future +""" +@kwdef struct DebugConfig + plot::DebugPlot=def_debug + # here can put more debug configurations +end \ No newline at end of file diff --git a/src/types/modelargs.jl b/src/types/modelargs.jl index 6848c05..8d2e09c 100644 --- a/src/types/modelargs.jl +++ b/src/types/modelargs.jl @@ -64,7 +64,7 @@ function UserArgs2D(; Ttol, Te=0, ΔT, time_step=:fixed, FLIP=1, PIC=0, constitu T2 = ϵ=="FP32" ? Float32 : Float64 # project default value folderdir = joinpath(abspath(project_path), project_name) - mkpath(folderdir); rm(folderdir, recursive=true, force=true); mkpath(folderdir) + isdir(folderdir) || mkpath(folderdir) cop_set = [:OS, :TS] bas_set = [:uGIMP, :linear, :bspline2, :bspline3] dev_set = [:CPU, :CUDA, :ROCm, :oneAPI, :Metal] @@ -150,7 +150,7 @@ function UserArgs3D(; Ttol, Te=0, ΔT, time_step=:fixed, FLIP=1, PIC=0, constitu T2 = ϵ=="FP32" ? Float32 : Float64 # project default value folderdir = joinpath(abspath(project_path), project_name) - mkpath(folderdir); rm(folderdir, recursive=true, force=true); mkpath(folderdir) + isdir(folderdir) || mkpath(folderdir) cop_set = [:OS, :TS] bas_set = [:uGIMP, :linear, :bspline2, :bspline3] dev_set = [:CPU, :CUDA, :ROCm, :oneAPI, :Metal] diff --git a/test/validation/druckerprager/3d_druckerprager.jl b/test/validation/druckerprager/3d_druckerprager.jl index 49c9d0c..434d9d9 100644 --- a/test/validation/druckerprager/3d_druckerprager.jl +++ b/test/validation/druckerprager/3d_druckerprager.jl @@ -8,12 +8,14 @@ | Affiliation: Risk Group, UNIL-ISTE | +==========================================================================================# +using WGLMakie using MaterialPointSolver using MaterialPointGenerator -using CairoMakie -using CUDA +using Metal +WGLMakie.activate!(; resize_to=:parent) +backend_name = :Metal -MaterialPointSolver.warmup(Val(:CUDA)) +MaterialPointSolver.warmup(Val(backend_name)) # 0.01000: 8000 pts # 0.00660: 27000 pts @@ -21,24 +23,24 @@ MaterialPointSolver.warmup(Val(:CUDA)) # 0.00400: 125000 pts | grid: x[-0.008, 0.808] y[-0.008, 0.056] z[-0.08, 0.108] # 0.00333: 216000 pts # 0.00250: 512000 pts -init_grid_space_x = 0.0025 -init_grid_space_y = 0.0025 -init_grid_space_z = 0.0025 +init_grid_space_x = 0.005 +init_grid_space_y = 0.005 +init_grid_space_z = 0.005 init_grid_range_x = [-0.02, 0.07] init_grid_range_y = [-0.02, 0.75] init_grid_range_z = [-0.02, 0.12] init_mp_in_space = 2 -init_T = 1 +init_T = 0.6 init_ρs = 2650 init_ν = 0.3 init_Ks = 7e5 init_Es = init_Ks * (3 * (1 - 2 * init_ν)) init_Gs = init_Es / (2 * (1 + init_ν)) init_ΔT = 0.5 * init_grid_space_x / sqrt(init_Es / init_ρs) -init_step = floor(init_T / init_ΔT / 50) +init_step = floor(init_T / init_ΔT / 100) init_ϕ = deg2rad(19.8) -init_FP = "FP64" -init_basis = :uGIMP +init_FP = "FP32" +init_basis = :bspline2 init_NIC = 27 # args setup @@ -54,9 +56,9 @@ args = UserArgs3D( hdf5 = false, hdf5_step = init_step, MVL = false, - device = :CUDA, + device = backend_name, coupling = :OS, - scheme = :MUSL, + scheme = :MLS, gravity = -9.8, ζs = 0, project_name = "3d_druckerprager", @@ -89,14 +91,16 @@ pts = meshbuilder(0 + dx / 2 : dx : 0.05 - dx / 2, 0 + dz / 2 : dz : 0.10 - dz / 2) mpρs = ones(size(pts, 1)) * init_ρs mp = UserParticle3D( - ϵ = init_FP, - phase = 1, - NIC = init_NIC, - dx = dx, - dy = dy, - dz = dz, - ξ = pts, - ρs = mpρs + ϵ = init_FP, + phase = 1, + NIC = init_NIC, + dx = dx, + dy = dy, + dz = dz, + ξ = pts, + ρs = mpρs, + affine = true, + mls = true ) # property setup @@ -139,19 +143,30 @@ bc = UserVBoundary3D( ) # solver setup -materialpointsolver!(args, grid, mp, attr, bc) +function tmpf(grid, mp, attr, bc, vid) + ξ0 = Array{Float32, 2}(mp.ξ0[vid, :]) + ξ1 = Array{Float32, 2}(mp.ξ[vid, :] ) -let - figfont = MaterialPointSolver.tnr - fig = Figure(size=(1200, 700), fonts=(; regular=figfont, bold=figfont), fontsize=30) - ax = Axis3(fig[1, 1], xlabel=L"x\ (m)", ylabel=L"y\ (m)", zlabel=L"z\ (m)", - aspect=:data, azimuth=0.2*π, elevation=0.1*π, xlabeloffset=60, zlabeloffset=80, - protrusions=100, xticks=(0:0.04:0.04), height=450, width=950) - pl1 = scatter!(ax, mp.ξ, color=log10.(mp.ϵq.+1), colormap=:jet, markersize=3, - colorrange=(0, 1)) - Colorbar(fig[1, 1], limits=(0, 1), colormap=:jet, size=16, ticks=0:0.5:1, spinewidth=0, - label=L"log_{10}(\epsilon_{II}+1)", vertical=false, tellwidth=false, width=200, - halign=:right, valign=:top, flipaxis=false) - display(fig) + val = @. sqrt((ξ1[:, 1] - ξ0[:, 1])^2 .+ + (ξ1[:, 2] - ξ0[:, 2])^2 .+ + (ξ1[:, 3] - ξ0[:, 3])^2) + return val end -rm(joinpath(abspath(args.project_path), args.project_name), recursive=true, force=true) \ No newline at end of file +plotconfig = DebugPlot(pointsize=1, axis=true, axfontsize=10, calculate=tmpf, cbname="disp", pointnum=1000) +debug = DebugConfig(plotconfig) +materialpointsolver!(args, grid, mp, attr, bc, debug) + +# let +# figfont = MaterialPointSolver.tnr +# fig = Figure(size=(1200, 700), fonts=(; regular=figfont, bold=figfont), fontsize=30) +# ax = Axis3(fig[1, 1], xlabel=L"x\ (m)", ylabel=L"y\ (m)", zlabel=L"z\ (m)", +# aspect=:data, azimuth=0.2*π, elevation=0.1*π, xlabeloffset=60, zlabeloffset=80, +# protrusions=100, xticks=(0:0.04:0.04), height=450, width=950) +# pl1 = scatter!(ax, mp.ξ, color=log10.(mp.ϵq.+1), colormap=:jet, markersize=3, +# colorrange=(0, 1)) +# Colorbar(fig[1, 1], limits=(0, 1), colormap=:jet, size=16, ticks=0:0.5:1, spinewidth=0, +# label=L"log_{10}(\epsilon_{II}+1)", vertical=false, tellwidth=false, width=200, +# halign=:right, valign=:top, flipaxis=false) +# display(fig) +# end +# rm(joinpath(abspath(args.project_path), args.project_name), recursive=true, force=true) \ No newline at end of file