diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json index 9bffab0..1eaffcb 100644 --- a/.devcontainer/devcontainer.json +++ b/.devcontainer/devcontainer.json @@ -1,42 +1,25 @@ +// For format details, see https://aka.ms/devcontainer.json. For config options, see the +// README at: https://github.com/devcontainers/templates/tree/main/src/python { - "name": "python dev", - "image": "mcr.microsoft.com/devcontainers/python:3.12", - "mounts": [ - "source=${localEnv:HOME}/.ssh,target=/home/vscode/.ssh,type=bind,consistency=cached" - ], - "customizations":{ + "name": "HydraxMPM", + // Or use a Dockerfile or Docker Compose file. More info: https://containers.dev/guide/dockerfile + "image": "mcr.microsoft.com/devcontainers/python:1-3.11-bullseye", + "features": {}, + + // Use 'postCreateCommand' to run commands after the container is created. + "postCreateCommand": "./.devcontainer/postCreateCommand.sh", + + // Configure tool-specific properties. + "customizations": { "vscode": { - "extensions":[ - "ms-python.python", - "ms-python.vscode-pylance", - "GitHub.copilot", - "oderwat.indent-rainbow", - "aaron-bond.better-comments", - "njpwerner.autodocstring", - "spywhere.guides", - "charliermarsh.ruff", - "streetsidesoftware.code-spell-checker", - "ms-toolsai.jupyter", - "littlefoxteam.vscode-python-test-adapter", - "VisualStudioExptTeam.vscodeintellicode" - ], + "extensions": ["ms-python.python", "editorconfig.editorconfig"], "settings": { - "python.pythonPath": ".venv/bin/python", - "editor.rulers": [120], - "python.testing.unittestArgs": [ - "-v", - "-s", - "./tests", - "-p", - "test_*.py" - ], - "python.testing.pytestEnabled": false, - "python.testing.unittestEnabled": true, - "editor.defaultFormatter": "charliermarsh.ruff", - "editor.formatOnSave": true + "python.testing.pytestArgs": ["tests"], + "python.testing.unittestEnabled": false, + "python.testing.pytestEnabled": true, + "python.defaultInterpreterPath": "/workspaces/Shongololo/.venv/bin/python", + "python.testing.pytestPath": "/workspaces/Shongololo/.venv/bin/pytest" } } - }, - "postCreateCommand": "curl -sSL https://install.python-poetry.org | python3 - && poetry install", - -} \ No newline at end of file + } +} diff --git a/.devcontainer/postCreateCommand.sh b/.devcontainer/postCreateCommand.sh new file mode 100755 index 0000000..0f40854 --- /dev/null +++ b/.devcontainer/postCreateCommand.sh @@ -0,0 +1,10 @@ +#! /usr/bin/env bash + +# Install uv +curl -LsSf https://astral.sh/uv/install.sh | sh + +# Install Dependencies +uv sync + +# Install pre-commit hooks +uv run pre-commit install --install-hooks diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9a1a546..5a50301 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -1,8 +1,8 @@ -name: ci +name: ci on: push: branches: - - master + - master - main permissions: contents: write @@ -26,11 +26,11 @@ jobs: run: uv python install 3.12 - name: Sync UV run: uv sync --all-extras - - run: echo "cache_id=$(date --utc '+%V')" >> $GITHUB_ENV + - run: echo "cache_id=$(date --utc '+%V')" >> $GITHUB_ENV - uses: actions/cache@v4 with: key: mkdocs-material-${{ env.cache_id }} path: .cache restore-keys: | mkdocs-material- - - run: uv run mkdocs gh-deploy --force \ No newline at end of file + - run: uv run mkdocs gh-deploy --force diff --git a/.github/workflows/run_tests.yml b/.github/workflows/run_tests.yml index 1f80bb1..a0a0c37 100644 --- a/.github/workflows/run_tests.yml +++ b/.github/workflows/run_tests.yml @@ -15,7 +15,7 @@ jobs: strategy: matrix: python-version: ["3.11", "3.12"] - + steps: - uses: actions/checkout@v4 diff --git a/.gitignore b/.gitignore index dec60ef..3b75fd0 100644 --- a/.gitignore +++ b/.gitignore @@ -11,6 +11,12 @@ __pycache__/ *.py[cod] *$py.class +core.* +uv.lock +imgui.ini +polyscope.ini +settings.json + flowdex/ # C extensions diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..827bb81 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,20 @@ +repos: + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: "v4.4.0" + hooks: + - id: check-case-conflict + - id: check-merge-conflict + - id: check-toml + - id: end-of-file-fixer + - id: check-added-large-files + - id: check-yaml + - id: debug-statements + - id: trailing-whitespace + # - id: check-docstring-first + + # - repo: https://github.com/astral-sh/ruff-pre-commit + # rev: "v0.6.3" + # hooks: + # - id: ruff + # args: [--exit-non-zero-on-fix] + # - id: ruff-format diff --git a/.vscode/settings.json b/.vscode/settings.json index 9b38853..a3a1838 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -4,4 +4,4 @@ ], "python.testing.unittestEnabled": false, "python.testing.pytestEnabled": true -} \ No newline at end of file +} diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..d9dcf5b --- /dev/null +++ b/Makefile @@ -0,0 +1,74 @@ +.PHONY: install +install: ## Install the virtual environment and install the pre-commit hooks + @echo "๐Ÿš€ Creating virtual environment using uv" + @uv sync + @uv run pre-commit install + +.PHONY: check +check: ## Run code quality tools. + @echo "๐Ÿš€ Checking lock file consistency with 'pyproject.toml'" + @uv lock --locked + @echo "๐Ÿš€ Linting code: Running pre-commit" + @uv run pre-commit run -a + # @echo "๐Ÿš€ Static type checking: Running mypy" + # @uv run mypy + @echo "๐Ÿš€ Checking for obsolete dependencies: Running deptry" + @uv run deptry . + +.PHONY: test +test: ## Test the code with pytest + @echo "๐Ÿš€ Testing code: Running pytest" + @uv run python -m pytest --cov --cov-config=pyproject.toml --cov-report=xml + + +# EXAMPLES +.PHONY: examples +examples: ## run all examples + @uv run examples/bvp/two_sphere_impact/two_sphere_impact.py + + +.PHONY: verification +verification: ## run all examples + @echo "๐Ÿ“ˆ Plotting verification tests" + @echo "MCC TRX CD CU '" + @uv run verification/mcc/trx_cu_cd.py + + + + +.PHONY: help +help: + @uv run python -c "import re; \ + [[print(f'\033[36m{m[0]:<20}\033[0m {m[1]}') for m in re.findall(r'^([a-zA-Z_-]+):.*?## (.*)$$', open(makefile).read(), re.M)] for makefile in ('$(MAKEFILE_LIST)').strip().split()]" + + +.DEFAULT_GOAL := help + + + + +# .PHONY: build +# build: clean-build ## Build wheel file +# @echo "๐Ÿš€ Creating wheel file" +# @uvx --from build pyproject-build --installer uv + +# .PHONY: clean-build +# clean-build: ## Clean build artifacts +# @echo "๐Ÿš€ Removing build artifacts" +# @uv run python -c "import shutil; import os; shutil.rmtree('dist') if os.path.exists('dist') else None" + +# .PHONY: publish +# publish: ## Publish a release to PyPI. +# @echo "๐Ÿš€ Publishing." +# @uvx twine upload --repository-url https://upload.pypi.org/legacy/ dist/* + +# .PHONY: build-and-publish +# build-and-publish: build publish ## Build and publish. + +# .PHONY: docs-test +# docs-test: ## Test if documentation can be built without warnings or errors +# @uv run mkdocs build -s + +# .PHONY: docs +# docs: ## Build and serve the documentation +# @uv run mkdocs serve \ No newline at end of file diff --git a/README.md b/README.md index 1c970db..a6cf7a1 100644 --- a/README.md +++ b/README.md @@ -1,27 +1,109 @@ +

-A Material Point Method solver developed with [JAX](https://jax.readthedocs.io/en/latest/quickstart.html) +
+

+ HydraxMPM Logo +

+
+

HydraxMPM

+

+ A JAX-powered Material Point Method & Single Integration Point simulation environment for granular materials +

+

+ Develop ยท Diagnose ยท Simulate +

+

+ ๐Ÿ“– Documentation (coming Soon) | โš™๏ธ Installation | ๐Ÿค Contributing (coming Soon) + +

-This code is still in pre-alpha, i.e., all models may not be feature-complete. +HydraxMPM integrates the **Material Point Method (MPM)** solver for large-scale granular dynamics simulations and **Single Integration Point (SIP)** testing , within one environment. Built on JAX, it leverages automatic differentiation and hardware acceleration (CPU/GPU/TPU) for research and development of numerical models capturing solid-like and fluid-like behavior of granular materials. -## Current features -Basic shape functions (linear, cubic) -Solvers (USL, APIC) -Materials ( Drucker Prager, Linear isotropic Elastic, Modified Cam Clay, Newtonian Fluid, $\mu (I)$ rheology) -Forces (Rigid body contact, gravity, slip and no slip boundaries) -## Installation instructions -- Install uv [here](https://docs.astral.sh/uv/getting-started/installation/) -- Clone repository `git clone git@github.com:GrainLearning/HydraxMPM.git && cd HydraxMPM` -- Install dependencies `uv sync` -- Run an example, e.g., `uv run examples/dambreak/dambreak.py`. Output is found in the `./examples/dambreak/` directory. +

+ + + + Animation demonstrating HydraxMPM simulation (light/dark mode versions) + +

-## Contributors: -Retief Lubbe, Hongyang Cheng -University of Twente, SMM -## Acknowledgements + + + + + +
+ + + MPM simulation SS (light mode) + + + + MPM simulation KE (light mode) +
+ + + + + +## Capabilities + +* ๐Ÿ”ฌ **Diagnose Constitutive Models:** Perform controlled SIP tests (triaxial, shear, and compression) for advanced constitutive model analysis. +* โ›ฐ๏ธ **Simulate Large-Scale Processes:** Model complex, large-deformation granular processes (e.g., landslides) using MPM. +* โš–๏ธ **Validate Across Scales:** Compare solid-l ike and fluid-like model behavior at both element (SIP) and system (MPM) levels. +* โˆ‡ **Gradient-Aware:** Utilize automatic differentiation for sensitivity analysis, inverse problems, and optimization. + +## Key Features + +* **Unified MPM & SIP:** Shared API facilitates rapid prototyping and validation. +* **High Performance:** JAX backend with JIT compilation. +* **Differentiable:** Enables advanced gradient-based studies. +* **Modular:** Designed for extensibility in research settings. +* **Solvers & Schemes:** Explicit MPM (USL) with FLIP/PIC, APIC, AFLIP transfer; Linear, Quadratic, Cubic B-spline basis functions. +* **Available Models:** Drucker-Prager, Modified Cam-Clay, Newtonian Fluid, Incompressible $\mu (I)$ rheology. +* **SIP Tests:** Triaxial (Drained/Undrained), Constant Pressure/Volume Shear, Isotropic Compression. +* **Contact & Boundaries:** Rigid body contact (penalty-based), slip/no-slip conditions. +* **Time Stepping & Stability:** Fixed and adaptive time stepping with Courantโ€“Friedrichsโ€“Lewy (CFL) condition. + + +## Installation + +1. **Install uv:** Follow instructions [here](https://docs.astral.sh/uv/getting-started/installation/). +2. **Clone & Install Dependencies:** + ```bash + git clone https://github.com/GrainLearning/HydraxMPM.git && cd HydraxMPM + uv sync + ``` +3. **Run Example:** + ```bash + uv run ./examples/dambreak/dambreak.py + ``` + *(Output in `./examples/dambreak/output`)* + + + +## ๐Ÿ‘ฅ Contributors: + +* Retief Lubbe (Soil Micro Mechanics group / University of Twente) +* Hongyang Cheng (Soil Micro Mechanics group / University of Twente) + +## ๐Ÿ™ Acknowledgements This research is part of the project TUSAIL [Training in Upscaling Particle Systems: Advancing Industry across Length-scales](https://tusail.eu)โ€ฏโ€ฏandโ€ฏhas received funding from the European Horizon2020 Framework Programme for research, technological development and demonstration under grant agreement ID 955661. +> [!WARNING] +> This is a research software under active development (pre-alpha). APIs and functionality are subject to change without notice. + + + + diff --git a/TODO.md b/TODO.md new file mode 100644 index 0000000..20db199 --- /dev/null +++ b/TODO.md @@ -0,0 +1,101 @@ + +## Core + + + +## Documentation + + +### How to guides + +- [ ] How to install + +- [ ] How to create a custom constitutive model + + -[ ] Example MU I / Phi I + +- [ ] How to create a custom element benchmark + +- [ ] How to visualize results + +- [ ] How to calibrate model... + +### Tutorials +- [ ] Granular column collapse + +### References + +- [ ] Examples +- [ ] Unit tests +- [ ] 3D 2D + + +### Explaination +- [ ] Motivation +- [ ] JAX +- [ ] Philsophy +- [ ] Model backgroud + +- [ ] Remove support for init by density.. + - [ ] Add documentation note, it can be done to overwrite init_state method... + +### Supported laws +- [ ] Drucker-Prager +- [ ] Modified Cam Clay +- [ ] Mu I rheology + + +## Unoffical +- [ ] UH model** + +- Fix post update for energies +- ADD FAQ on how to use CPU GPU os.env flat + + + +- Add error handler in mpm solver to output data if throws error + - After cleanup +- write some standard checks, unit tests, to see if models are functioning properly + - TRX compression + - ISO extension +- Make a list of experimental validation simulations +- Make file to run all checks + +- Open questions + - How to organize plots? worth automated way + + +- Add note that large velocity gradients may be due to, a large spatial variation in the materila points velocity.. + - particles moving opposite direction, or shearing.. + Excessively large gradients can destabilize simulations (e.g., causing unphysical stresses or grid artifacts) + + + +- [ ] Add gifs + - [ ] SIP test moving 4 models + - [ ] MPM column collapse (4 models) + +- [ ] Badges + - [ ] license + - [ ] build + - [ ] docs + + + +- ADD License in files +- Later +-[ ] Add detailed table of models in documentation[] +- Mention the code realizies on Equionox +- Prelimnary is knowledge of jax + + + + +- MCC change R to OCR + \ No newline at end of file diff --git a/benchmarks/benchmark_usl.py b/benchmarks/benchmark_usl.py deleted file mode 100644 index 2b5acc1..0000000 --- a/benchmarks/benchmark_usl.py +++ /dev/null @@ -1,84 +0,0 @@ -import jax -import jax.numpy as jnp -import pymudokon as pm -import pyperf - - -def run_p2g_batch(particles, nodes, shapefunctions): - nodes = pm.solvers.usl.p2g_batch( - nodes=nodes, - particles=particles, - shapefunctions=shapefunctions, - dt=0.1, - ) - return nodes - - -def run_p2g(particles, nodes, shapefunctions): - nodes = pm.solvers.usl.p2g( - nodes=nodes, - particles=particles, - shapefunctions=shapefunctions, - dt=0.1, - ) - return nodes - - -def create_system(num_particles, cell_size): - key = jax.random.key(0) - positions = jax.random.uniform(key, (num_particles, 2)) - - particles = pm.Particles.create(positions) - nodes = pm.Nodes.create( - origin=jnp.array([0.0, 0.0]), - end=jnp.array([1.0, 1.0]), - node_spacing=cell_size, - ) - shapefunctions = pm.LinearShapeFunction.create(num_particles, 2) - - shapefunctions = shapefunctions.get_interactions(particles, nodes) - - shapefunctions = shapefunctions.calculate_shapefunction(nodes) - - return particles, nodes, shapefunctions - - -runner = pyperf.Runner(loops=2) - -runner.warmups = 2 # Number of warm-up runs -runner.samples = 2 # Number of benchmark runs - - -run_p2g_batch_jitted = jax.jit(run_p2g_batch) - -run_p2g_jitted = jax.jit(run_p2g) - - -for num_particles in [10000]: - for cell_size in [0.001]: - particles, nodes, shapefunctions = create_system(num_particles, cell_size) - runner.bench_func( - f"p2g_batch/{num_particles}/{cell_size}", - lambda: jax.block_until_ready( - run_p2g_batch_jitted(particles, nodes, shapefunctions) - ), - ) - - -for num_particles in [10000]: - for cell_size in [0.001]: - particles, nodes, shapefunctions = create_system(num_particles, cell_size) - runner.bench_func( - f"p2g/{num_particles}/{cell_size}", - lambda: jax.block_until_ready( - run_p2g_jitted(particles, nodes, shapefunctions) - ), - ) - -# runner = pyperf.Runner() -# for size_pow10 in range(5): -# size = pow(10, size_pow10) - -# in1 = [4.2 for i in range(size)] -# in2 = [2.4 for i in range(size)] -# runner.bench_func(f'sum_append/{size}', lambda: vector_sum_append(in1, in2)) diff --git a/docs/_static/extra.css b/docs/_static/extra.css index a9847c7..0ee8983 100644 --- a/docs/_static/extra.css +++ b/docs/_static/extra.css @@ -3,7 +3,7 @@ .doc-symbol-function::after { content: "function"; } - + .doc-symbol-method::after { content: "method"; } @@ -22,14 +22,14 @@ --md-accent-fg-color: #A5F3FC; --md-accent-bg-color: hsla(0, 0%, 100%, 1); --md-accent-bg-color--light: hsla(0, 0%, 100%, 0.7); - + } /* Indentation. */ div.doc-contents:not(.first) { padding-left: 30px; border-left: .05rem solid var(--md-typeset-table-color); - + } /* Mark external links as such. */ @@ -221,9 +221,9 @@ h2 { /* background-color: #A21CAF; */ /* border: 1px solid #A21CAF; */ /* background-color: #A21CAF; */ - + /* padding: 4px 20px; */ - + color: #A5F3FC; font-weight: bold; line-height: 1.2; diff --git a/docs/_static/hydraxmpm.png b/docs/_static/hydraxmpm.png new file mode 100644 index 0000000..173c7f3 Binary files /dev/null and b/docs/_static/hydraxmpm.png differ diff --git a/docs/_static/mpm_models_ke_dark.png b/docs/_static/mpm_models_ke_dark.png new file mode 100644 index 0000000..fa5e33a Binary files /dev/null and b/docs/_static/mpm_models_ke_dark.png differ diff --git a/docs/_static/mpm_models_ke_light.png b/docs/_static/mpm_models_ke_light.png new file mode 100644 index 0000000..8488eb5 Binary files /dev/null and b/docs/_static/mpm_models_ke_light.png differ diff --git a/docs/_static/mpm_models_ss_dark.png b/docs/_static/mpm_models_ss_dark.png new file mode 100644 index 0000000..7992b06 Binary files /dev/null and b/docs/_static/mpm_models_ss_dark.png differ diff --git a/docs/_static/mpm_models_ss_light.png b/docs/_static/mpm_models_ss_light.png new file mode 100644 index 0000000..84e2e48 Binary files /dev/null and b/docs/_static/mpm_models_ss_light.png differ diff --git a/docs/_static/sip_animation_dark.gif b/docs/_static/sip_animation_dark.gif new file mode 100644 index 0000000..8043c2b Binary files /dev/null and b/docs/_static/sip_animation_dark.gif differ diff --git a/docs/_static/sip_animation_light.gif b/docs/_static/sip_animation_light.gif new file mode 100644 index 0000000..58baeff Binary files /dev/null and b/docs/_static/sip_animation_light.gif differ diff --git a/docs/advance/custom_models.md b/docs/advance/custom_models.md index 878b126..a59c06a 100644 --- a/docs/advance/custom_models.md +++ b/docs/advance/custom_models.md @@ -1 +1 @@ -# Custom models \ No newline at end of file +# Custom models diff --git a/docs/api/config.md b/docs/api/config.md index 79b536a..e5e1b8a 100644 --- a/docs/api/config.md +++ b/docs/api/config.md @@ -1,30 +1,10 @@ -# Configuration +# Configuration -A configuration data object is passed down to all `HydraxMPM` objects. This pre-allocates necessary memory for simulations. +A configuration data object is passed down to `HydraxMPM` solvers. The same configuration is used for both single integration point solver and material point method solver. -Two config types are available: +## API +::: solvers.config.Config -* **`MPMConfig`**: For Material Point Method (MPM) simulations. -* **`IPConfig`**: For single integration point element tests. - - - -!!! tip "Modifying config" - - Lorem ipsum dolor sit amet, consectetur adipiscing elit. Nulla et euismod - nulla. Curabitur feugiat, tortor non consequat finibus, justo purus auctor - massa, nec semper lorem quam in massa. - ---- - -## MPMConfig -::: config.mpm_config.MPMConfig - - -## IPConfig -::: config.mpm_config.MPMConfig - - -**Integration Point Config** +::: solvers.mpm_solver.MPMSolver diff --git a/docs/api/materials.md b/docs/api/materials.md index 0d22206..8895655 100644 --- a/docs/api/materials.md +++ b/docs/api/materials.md @@ -50,12 +50,12 @@ $$ $$ -State boundary layer +State boundary layer Distance between NCL and current state during yielding $$ -\ln v_m = \ln Z +\ln v_m = \ln Z - \lambda \ln \left( \frac{p + p_s}{1 + p_s} \right) - (\lambda - \kappa) \ln \left( 1+ \frac{m^2}{M^2} \right) $$ --> diff --git a/docs/api/math.md b/docs/api/math.md index e8d8752..75eabe8 100644 --- a/docs/api/math.md +++ b/docs/api/math.md @@ -23,4 +23,4 @@ pressure_stack = pm.get_pressure_stack(stress_stack) # [1000,1000] [Pa] and shap ``` -::: utils.math_helpers --> \ No newline at end of file +::: utils.math_helpers --> diff --git a/docs/api/nodes.md b/docs/api/nodes.md index e69de29..5bf6847 100644 --- a/docs/api/nodes.md +++ b/docs/api/nodes.md @@ -0,0 +1,5 @@ +## Grid +::: grid.grid.Grid + + diff --git a/docs/api/particles.md b/docs/api/particles.md index 1a9fa36..1bc27ec 100644 --- a/docs/api/particles.md +++ b/docs/api/particles.md @@ -1,29 +1,2 @@ - +## API +::: material_points.material_points.MaterialPoints diff --git a/docs/api/shapefunctions.md b/docs/api/shapefunctions.md index 11ea8e5..8662159 100644 --- a/docs/api/shapefunctions.md +++ b/docs/api/shapefunctions.md @@ -6,7 +6,7 @@ The shapefunction quantites are stored in this class. Two shapefunctions are currently implemented, namely, cubic splines and linear shape functions [^1]. !!! note - It is important to specify the `shapefuction_type` in `MPMConfig` before creating the shapefunction class. + It is important to specify the `shapefuction_type` in `MPMConfig` before creating the shapefunction class. ::: shapefunctions.shapefunctions @@ -20,7 +20,7 @@ ```py import hydraxmpm as hdx - + config = hdx.MPMConfig( shapefunction_type="cubic" ) @@ -29,7 +29,7 @@ ``` - + @@ -43,7 +43,7 @@ ```py import hydraxmpm as hdx - + config = hdx.MPMConfig( shapefunction="linear" ) @@ -57,4 +57,4 @@ ------------ **References:** -[^1]: De Vaucorbeil, Alban, et al. 'Material point method after 25 years: theory, implementation, and applications.' --> \ No newline at end of file +[^1]: De Vaucorbeil, Alban, et al. 'Material point method after 25 years: theory, implementation, and applications.' --> diff --git a/docs/basics.md b/docs/basics.md index e512586..b875158 100644 --- a/docs/basics.md +++ b/docs/basics.md @@ -7,4 +7,4 @@ HydraxMPM follows an array of structures design. That is, a set of scalars has a shape `(number of material points,)`, vectors `(number of material points, dimension)`, tensors `(number of material points, dimension)`. These arrays are typically denoted as stacks (e.g., `mass_stack`, `velocity_stack`). Some arrays are always defined in 3D, e.g., the deformation gradient `F_stack`, which has a shape `(number of material points,3)`. -- Objects are defined as [dataclasses](https://docs.python.org/3/library/dataclasses.html), common practice in JAX \ No newline at end of file +- Objects are defined as [dataclasses](https://docs.python.org/3/library/dataclasses.html), common practice in JAX diff --git a/docs/extra.css b/docs/extra.css index 2985cd1..83df50a 100644 --- a/docs/extra.css +++ b/docs/extra.css @@ -2,7 +2,7 @@ div.doc-contents:not(.first) { padding-left: 25px; border-left: .05rem solid var(--md-typeset-table-color); - + } /* Mark external links as such. */ @@ -32,9 +32,9 @@ a.autorefs-external:hover::after { .doc-symbol-function::after { content: "function"; } - + .doc-symbol-method::after { content: "method"; } - -} \ No newline at end of file + +} diff --git a/docs/how-tos/initialize_material_points.md b/docs/how-tos/initialize_material_points.md index 50257a6..fd28e48 100644 --- a/docs/how-tos/initialize_material_points.md +++ b/docs/how-tos/initialize_material_points.md @@ -1,3 +1,3 @@ -This is empty and sad, wait until its happier :( \ No newline at end of file +This is empty and sad, wait until its happier :( diff --git a/docs/index.md b/docs/index.md index ea5a51e..ee1e196 100644 --- a/docs/index.md +++ b/docs/index.md @@ -21,4 +21,4 @@ ## Whats in the box? -- Plane strain \ No newline at end of file +- Plane strain diff --git a/docs/other/animated_sip.py b/docs/other/animated_sip.py new file mode 100644 index 0000000..32d52a2 --- /dev/null +++ b/docs/other/animated_sip.py @@ -0,0 +1,233 @@ +import os + + +os.environ["CUDA_VISIBLE_DEVICES"] = "" # Hide GPUs from all libraries (including JAX) +os.environ["JAX_PLATFORMS"] = "cpu" + +import jax +import jax.numpy as jnp + +import hydraxmpm as hdx + +import matplotlib.pyplot as plt + + +import matplotlib.pyplot as plt +import matplotlib as mpl + +import scienceplots + +plt.style.use(["science", "no-latex"]) + + +import equinox as eqx +import os + +jax.config.update("jax_enable_x64", True) + + +dir_path = os.path.dirname(os.path.realpath(__file__)) + + +mpl.rcParams["lines.linewidth"] = 2 + +cycle = plt.rcParams["axes.prop_cycle"].by_key()["color"] + + +p_0 = 1000.0 + +fric_angle = 19.8 +fric_angle_rad = jnp.deg2rad(fric_angle) +rho_0 = 2650.0 +rho_p = 2650.0 + +K = 7e5 # [Pa] + +# matching Mohr-Coulomb criterion under triaxial extension conditions +mu = ( + 6 * jnp.sin(fric_angle_rad) / (jnp.sqrt(3) * (3 + jnp.sin(jnp.deg2rad(fric_angle)))) +) + + +models = ( + hdx.ModifiedCamClay( + nu=0.3, + M=mu * jnp.sqrt(3), + lam=0.0058, + kap=0.0012, + R=1.0, + rho_0=rho_0, + rho_p=rho_p, + settings=dict(atol=1e-3), + other=dict(label="Modified Cam-Clay", zorder=-1, color=cycle[0]), + ), + hdx.DruckerPrager( + nu=0.3, + K=K, + mu_1=mu, + rho_0=rho_0, + rho_p=rho_p, + other=dict(label="Drucker-Prager", zorder=-1, color=cycle[2]), + ), + hdx.MuI_incompressible( + other=dict(label="$ \\mu (I)$-rheology", ls="-", zorder=-1, color=cycle[1]), + mu_s=mu, + mu_d=1.9, + I_0=0.279, + K=K, + d=0.00125, + rho_p=rho_p, + rho_0=rho_0, + ), + hdx.NewtonFluid( + other=dict(label="Newton Fluid", ls="-", zorder=0, color=cycle[4]), + K=K, + viscosity=0.002, + alpha=7.0, + rho_0=rho_0, + ), +) + + +dgamma_dt_start = 0.5 +dgamma_dt_end = 50.0 + +print(dgamma_dt_start, dgamma_dt_end) + +x_slow = dgamma_dt_start * 2 +x_fast = dgamma_dt_end * 2 + +num_steps_s1 = 25000 +x1_stack = jnp.ones(num_steps_s1) * x_slow + + +num_steps_s2 = 5000 +x2_stack = jnp.linspace(x_slow, x_fast, num_steps_s2) + + +num_steps_s3 = 15000 +x3_stack = jnp.ones(num_steps_s3) * x_fast + + +x_total_stack = jnp.concat((x1_stack, x2_stack, x3_stack)) + +sip_benchmarks = ( + hdx.S_CD( + deps_xy_dt=jnp.concat((x1_stack, x2_stack, x3_stack)), + # deps_xy_dt=4.0, + num_steps=x_total_stack.shape[0], + p0=p_0, + init_material_points=True, + other=dict(type="S_CD"), + ), +) + + +fig, axes = plt.subplots( + ncols=3, + figsize=(8, 3), + dpi=200, + layout="constrained", +) +dt = 0.00001 + +for mi, model in enumerate(models): + for ie, sip_benchmark in enumerate(sip_benchmarks): + solver = hdx.SIPSolver( + material_points=hdx.MaterialPoints( + p_stack=jnp.array([p_0]), + ), + output_dict=("p_stack", "q_stack", "eps_v_stack", "specific_volume_stack"), + constitutive_law=model, + sip_benchmarks=sip_benchmark, + ) + t_stack = jnp.arange(0, sip_benchmark.num_steps) * dt + + solver = solver.setup() + + (p_stack, q_stack, eps_v_stack, specific_volume_stack) = solver.run(dt=dt) + + if ie == 0: + label = model.other.get("label", "") + else: + label = None + + print(eps_v_stack) + p_stack = p_stack + q_stack = q_stack + + hdx.make_plot( + axes.flat[0], + p_stack, + q_stack, + xlabel="$p$ [Pa]", + ylabel="$q$ [Pa]", + color=model.other.get("color", "black"), + label=label, + ls=sip_benchmark.other.get("ls", "-"), + zorder=model.other.get("zorder", 1), + start_end_markersize=7, + ) + + hdx.make_plot( + axes.flat[1], + t_stack, + q_stack / p_stack, + xlabel="$t$ [s]", + ylabel="$q/p$ [-]", + color=model.other["color"], + ls="-", + zorder=model.other["zorder"], + start_end_markersize=7, + ) + + hdx.make_plot( + axes.flat[2], + t_stack, + eps_v_stack, + xlabel="$t$ [s]", + ylabel="$\\varepsilon_v$ [-]", + color=model.other["color"], + ls="-", + zorder=model.other["zorder"], + start_end_markersize=7, + ) + print(f"{model.other['label']} {sip_benchmark.other['type']} done..") + + +# plot CSL +p_aux = jnp.arange(0, 1500, 10) +q_csl = models[0].CSL_q_p(p_aux) + +axes.flat[0].plot(p_aux, q_csl, "r-", lw=1.0, zorder=-1, label="CSL") + +t_aux = jnp.arange(0, 0.5, 0.01) +axes.flat[1].plot(t_aux, jnp.ones_like(t_aux) * models[0].M, "r-", lw=1.0, zorder=-3) + + +# set limits +axes.flat[0].margins(0.05) +axes.flat[1].margins(0.05) +axes.flat[2].margins(0.05) + + +axes.flat[0].grid(True) +axes.flat[1].grid(True) +axes.flat[2].grid(True) + + +def create_legend(fig, ncols=5): + lines_labels = [ax.get_legend_handles_labels() for ax in fig.axes] + lines, labels = [sum(lol, []) for lol in zip(*lines_labels)] + fig.legend(lines, labels, ncols=ncols, loc="outside lower center") + return fig + + +create_legend(fig, 5) + +for i, label in enumerate(["(a)", "(b)", "(c)"]): + axes.flat[i].set_title(label, y=0, pad=-35, verticalalignment="top") +# fig.savefig(dir_path + "/plots/sip_pressure_control.png") + + +# ax = plt.gca() # get axis handle \ No newline at end of file diff --git a/docs/tutorials/1_granular_column.md b/docs/tutorials/1_granular_column.md index c43e978..4a0a6d5 100644 --- a/docs/tutorials/1_granular_column.md +++ b/docs/tutorials/1_granular_column.md @@ -16,7 +16,7 @@ It is recommended that you follow along with the tutorial. Stuck or in a rush? t ## Learning objectives -By the end of this tutorial you will be able to +By the end of this tutorial you will be able to - Understand the overall code style of a HydraxMPM script @@ -26,23 +26,23 @@ Create a project directory containing the driver script and two sub-folders with ``` # driver script -/tutorial/t1_granular_column.py +/tutorial/t1_granular_column.py # output of gravity pack -/tutorial/output/t1_pack/ +/tutorial/output/t1_pack/ # output of collapse -/tutorial/output/t2_collapse +/tutorial/output/t2_collapse ``` The output is stored as [vtk](https://docs.vtk.org/en/latest/design_documents/VTKFileFormats.html) files in the `/tutorial/output/t1_pack/` and `/tutorial/output/t2_collapse/` folders, respectively. ## Step 2: import modules -Import HydraxMPM and supporting the JAX dependencies. +Import HydraxMPM and supporting the JAX dependencies. ```python {hl_lines="3"} ---8<-- "t1_granular_column.py:2:5" +# --8<-- "t1_granular_column.py:2:5" ``` @@ -50,11 +50,11 @@ HydraxMPM is top-level package. It contains all the bits and pieces under the sa ## Step 3: create points representing initial granular column -Create rectangular column a rectangular column of material points. Particles are spaced evenly given a cell size and the number of particles per cell (in one direction). We pad so that material points do not touch the boundary. +Create rectangular column a rectangular column of material points. Particles are spaced evenly given a cell size and the number of particles per cell (in one direction). We pad so that material points do not touch the boundary. -```python +```python ---8<-- "t1_granular_column.py:8:23" +# --8<-- "t1_granular_column.py:8:23" ``` @@ -62,9 +62,9 @@ Create rectangular column a rectangular column of material points. Particles are There are several ways of initializing material points. See the [how-to initialize material points](/how-tos/initialize_material_points) . ??? Tip - [Matplotlib](https://matplotlib.org/) or [Pyvista](https://docs.pyvista.org/) may be used visualize initial positions. To visualize the particles copy and paste the following (note plt.show works only when GUI is active). + [Matplotlib](https://matplotlib.org/) or [Pyvista](https://docs.pyvista.org/) may be used visualize initial positions. To visualize the particles copy and paste the following (note plt.show works only when GUI is active). ```python - + import matplotlib.pyplot as plt plt.scatter(*position_stack.T ) # unpack x and y and plot @@ -75,9 +75,9 @@ Create rectangular column a rectangular column of material points. Particles are This is a juicy sandwich of all common general simulation parameters. -```python +```python ---8<-- "t1_granular_column.py:24:37" +# --8<-- "t1_granular_column.py:24:37" ``` Ok... it seems like a lot, but actually pretty straight forward: @@ -92,7 +92,7 @@ This is a juicy sandwich of all common general simulation parameters. - `shapefunction` the interpolation type from particle to grid. Two shapefunctions are supported, either `linear` or `cubic`. Cubic is better in most cases. - `num_steps` the total iteration count - `store_every` output every nth step - - `default_gpu_id` if you are working on a shared GPU workstation, this is the parameter you change to avoid making the other person(s) angry! Run `nvidia-smi`, in your terminal to find the id of an empty GPU. + - `default_gpu_id` if you are working on a shared GPU workstation, this is the parameter you change to avoid making the other person(s) angry! Run `nvidia-smi`, in your terminal to find the id of an empty GPU. - `dt` is a constant time step - `file=__file__` this records the path of your driver script, which is important to save relative output in the correct folder. @@ -107,9 +107,9 @@ This is a juicy sandwich of all common general simulation parameters. Lets see the summary of the config -```python +```python ---8<-- "t1_granular_column.py:38:38" +# --8<-- "t1_granular_column.py:38:38" ``` If all went well you should get the following output: @@ -134,28 +134,28 @@ total time: 12.000000000000002 Lets create a material with some initial bulk density and a particle density -```python + ---8<-- "t1_granular_column.py:40:44" +# --8<-- "t1_granular_column.py:40:44" -``` + Now the material dataclass -```python + ---8<-- "t1_granular_column.py:45:56" +# --8<-- "t1_granular_column.py:45:56" -``` + Lets break this down: - `hdx.ModifiedCamClay` is the Modified Cam Clay material class - We always pass `config` along when creating HydraxMPM dataclasses. -- The only non-parameter input is the reference solid volume fraction. +- The only non-parameter input is the reference solid volume fraction. -!!! Tip +!!! Tip Materials in HydraxMPM are normally initialize either with a reference pressure or reference solid volume fraction. Given one, we can find the other. ??? Tip @@ -165,66 +165,66 @@ Lets break this down: ## Step 5: create the `Nodes` and `Particles` data classes -We initialize particle stresses to match the known pressure state, given a predefined solid-volume fraction above. +We initialize particle stresses to match the known pressure state, given a predefined solid-volume fraction above. -Here we finally make use of JAX [vmap](https://jax.readthedocs.io/en/latest/_autosummary/jax.vmap.html), to get the stress tensor. +Here we finally make use of JAX [vmap](https://jax.readthedocs.io/en/latest/_autosummary/jax.vmap.html), to get the stress tensor. -```python + ---8<-- "t1_granular_column.py:58:61" +# --8<-- "t1_granular_column.py:58:61" -``` + Pass all positions and stresses to the `Particles` dataclass. -```python + ---8<-- "t1_granular_column.py:63:65" +# --8<-- "t1_granular_column.py:63:65" -``` + We can create background grid nodes via the config. -```python + ---8<-- "t1_granular_column.py:67:67" +# --8<-- "t1_granular_column.py:67:67" -``` + The `discretize` function determines initial particle volume by dividing the number of particles in a cell by the cell size. -```python +```python ---8<-- "t1_granular_column.py:69:72" +# --8<-- "t1_granular_column.py:69:72" ``` ## Step 6: create the `Gravity` and `Domain` -Gravity is slowly ramped up +Gravity is slowly ramped up -```python +```python ---8<-- "t1_granular_column.py:73:77" +# --8<-- "t1_granular_column.py:73:77" ``` Creating the outside domain box. -```python + ## Step 7: create the `Solver` The solver determines how the background grid and material points interact. -```python + ??? Tip We recommend using the `hdx.ASFLIP` for granular materials @@ -235,11 +235,11 @@ The solver determines how the background grid and material points interact. - See [available callback functions]() -```python + @@ -247,19 +247,19 @@ The solver determines how the background grid and material points interact. -```python + ## Step 10: Granular column collapse -```python + diff --git a/docs/usage/index.md b/docs/usage/index.md index 56b1d8b..1c54e5d 100644 --- a/docs/usage/index.md +++ b/docs/usage/index.md @@ -4,7 +4,7 @@ ## Sampling particles -