From d246a0c5fdd7668c155b61f6b174e1b541a09b22 Mon Sep 17 00:00:00 2001 From: VanillaBrooks Date: Sun, 30 Apr 2023 11:28:02 -0700 Subject: [PATCH 01/14] Files generated by PkgTemplates PkgTemplates version: 0.7.34 --- .github/workflows/CI.yml | 36 ++++++++++++++++++++++++++++++ .github/workflows/CompatHelper.yml | 16 +++++++++++++ .github/workflows/TagBot.yml | 15 +++++++++++++ .gitignore | 1 + LICENSE | 21 +++++++++++++++++ Project.toml | 13 +++++++++++ README.md | 3 +++ src/spectralGPU.jl | 5 +++++ test/runtests.jl | 6 +++++ 9 files changed, 116 insertions(+) create mode 100644 .github/workflows/CI.yml create mode 100644 .github/workflows/CompatHelper.yml create mode 100644 .github/workflows/TagBot.yml create mode 100644 .gitignore create mode 100644 LICENSE create mode 100644 Project.toml create mode 100644 README.md create mode 100644 src/spectralGPU.jl create mode 100644 test/runtests.jl diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml new file mode 100644 index 0000000..68c126b --- /dev/null +++ b/.github/workflows/CI.yml @@ -0,0 +1,36 @@ +name: CI +on: + push: + branches: + - main + tags: ['*'] + pull_request: +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: only if it is a pull request build. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - '1.0' + - '1.8' + - 'nightly' + os: + - ubuntu-latest + arch: + - x64 + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v1 + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml new file mode 100644 index 0000000..cba9134 --- /dev/null +++ b/.github/workflows/CompatHelper.yml @@ -0,0 +1,16 @@ +name: CompatHelper +on: + schedule: + - cron: 0 0 * * * + workflow_dispatch: +jobs: + CompatHelper: + runs-on: ubuntu-latest + steps: + - name: Pkg.add("CompatHelper") + run: julia -e 'using Pkg; Pkg.add("CompatHelper")' + - name: CompatHelper.main() + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }} + run: julia -e 'using CompatHelper; CompatHelper.main()' diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml new file mode 100644 index 0000000..f49313b --- /dev/null +++ b/.github/workflows/TagBot.yml @@ -0,0 +1,15 @@ +name: TagBot +on: + issue_comment: + types: + - created + workflow_dispatch: +jobs: + TagBot: + if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' + runs-on: ubuntu-latest + steps: + - uses: JuliaRegistries/TagBot@v1 + with: + token: ${{ secrets.GITHUB_TOKEN }} + ssh: ${{ secrets.DOCUMENTER_KEY }} diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..b067edd --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +/Manifest.toml diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..5ecf4b1 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2023 brooks + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..ecc2b91 --- /dev/null +++ b/Project.toml @@ -0,0 +1,13 @@ +name = "spectralGPU" +uuid = "1f941c8d-9269-406b-be03-2a5a5e88bc1c" +authors = ["brooks"] +version = "1.0.0-DEV" + +[compat] +julia = "1" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] diff --git a/README.md b/README.md new file mode 100644 index 0000000..7714a85 --- /dev/null +++ b/README.md @@ -0,0 +1,3 @@ +# spectralGPU + +[![Build Status](https://github.com/vanillabrooks/spectralGPU.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/vanillabrooks/spectralGPU.jl/actions/workflows/CI.yml?query=branch%3Amain) diff --git a/src/spectralGPU.jl b/src/spectralGPU.jl new file mode 100644 index 0000000..0564f2b --- /dev/null +++ b/src/spectralGPU.jl @@ -0,0 +1,5 @@ +module spectralGPU + +# Write your package code here. + +end diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..0e7a4e6 --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,6 @@ +using spectralGPU +using Test + +@testset "spectralGPU.jl" begin + # Write your tests here. +end From 87c9563863581761bb9b1fac7b1d077c7cd12f73 Mon Sep 17 00:00:00 2001 From: VanillaBrooks Date: Sun, 30 Apr 2023 14:26:33 -0700 Subject: [PATCH 02/14] initial tests of fourier transforms etc --- .../fft_debug.py-checkpoint.ipynb | 107 ++ jupyter/fft_debug.py.ipynb | 225 +++ pluto/fft_debug | 1377 +++++++++++++++++ src/fft.jl | 24 + src/markers.jl | 11 + src/mesh.jl | 81 + src/spectralGPU.jl | 9 +- test/runtests.jl | 11 +- 8 files changed, 1842 insertions(+), 3 deletions(-) create mode 100644 jupyter/.ipynb_checkpoints/fft_debug.py-checkpoint.ipynb create mode 100644 jupyter/fft_debug.py.ipynb create mode 100644 pluto/fft_debug create mode 100644 src/fft.jl create mode 100644 src/markers.jl create mode 100644 src/mesh.jl diff --git a/jupyter/.ipynb_checkpoints/fft_debug.py-checkpoint.ipynb b/jupyter/.ipynb_checkpoints/fft_debug.py-checkpoint.ipynb new file mode 100644 index 0000000..2282b18 --- /dev/null +++ b/jupyter/.ipynb_checkpoints/fft_debug.py-checkpoint.ipynb @@ -0,0 +1,107 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "3d802773", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from numpy import fft" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "2a67eb8f", + "metadata": {}, + "outputs": [], + "source": [ + "N = 8\n", + "rank = 1\n", + "Np = 8\n", + "\n", + "kx = fft.fftfreq(N, 1./N)\n", + "ky = kx[rank*Np:(rank+1)*Np].copy()\n", + "kz = kx[:(N//2+1)].copy()\n", + "kz[-1] *= -1\n", + "\n", + "#(3, N, Np, Np)\n", + "K = np.array(np.meshgrid(kx, ky, kz, indexing='ij'), dtype=int)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "b329139c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([], shape=(3, 8, 0, 5), dtype=int64)" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "K" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "0e0fba90", + "metadata": {}, + "outputs": [], + "source": [ + "def fftn_single(u, fu):\n", + " fu [:] = rfftn (u , axes =( 0 ,1 , 2 ) )\n", + "\n", + "def ifftn_mpi(fu, u):\n", + " u[:] = irfftn ( fu , axes =( 0 ,1 , 2 ) )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "44cea490", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f26bb8c7", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/jupyter/fft_debug.py.ipynb b/jupyter/fft_debug.py.ipynb new file mode 100644 index 0000000..ca58391 --- /dev/null +++ b/jupyter/fft_debug.py.ipynb @@ -0,0 +1,225 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 54, + "id": "c5574ad8", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from numpy import fft\n", + "from math import pi\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "id": "b92b2e5d", + "metadata": {}, + "outputs": [], + "source": [ + "N = 8\n", + "rank = 0\n", + "Np = 8\n", + "\n", + "kn = N//2 + 1\n", + "\n", + "kx = fft.fftfreq(N, 1./N)\n", + "ky = kx[rank*Np:(rank+1)*Np].copy()\n", + "kz = kx[:(N//2+1)].copy()\n", + "kz[-1] *= -1\n", + "\n", + "#(3, N, Np, Np)\n", + "K = np.array(np.meshgrid(kx, ky, kz, indexing='ij'), dtype=int)\n", + "\n", + "X = np.mgrid[rank*Np:(rank+1)*Np, :N, :N].astype(float)*2*pi/N" + ] + }, + { + "cell_type": "code", + "execution_count": 56, + "id": "e23269fc", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(8,)\n", + "(8,)\n", + "(5,)\n" + ] + } + ], + "source": [ + "print(kx.shape)\n", + "print(ky.shape)\n", + "print(kz.shape)\n", + "K;" + ] + }, + { + "cell_type": "code", + "execution_count": 68, + "id": "f3e91211", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "5.497787143782138" + ] + }, + "execution_count": 68, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.max(X[0])" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "id": "f37638b8", + "metadata": {}, + "outputs": [], + "source": [ + "def fftn_single(u, fu):\n", + " fu[:] = fft.rfftn (u , axes =( 0 ,1 , 2 ) )\n", + "\n", + "def ifftn_single(fu, u):\n", + " u[:] = fft.irfftn ( fu , axes =( 0 ,1 , 2 ) )" + ] + }, + { + "cell_type": "code", + "execution_count": 58, + "id": "acc98409", + "metadata": {}, + "outputs": [], + "source": [ + "u = np.sin(X[0])\n", + "uback = np.zeros((N, N, N))\n", + "uhat = np.zeros((N, N, kn), dtype = complex)" + ] + }, + { + "cell_type": "code", + "execution_count": 59, + "id": "487cfa3b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0j\n", + "(8, 8, 5)\n" + ] + } + ], + "source": [ + "fftn_single(u, uhat);\n", + "ifftn_single(uhat, uback);\n", + "\n", + "print(uhat.sum())\n", + "print(uhat.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": 64, + "id": "61c1500f", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 64, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhsAAAGiCAYAAABOCgSdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAX5UlEQVR4nO3dfWyV9f3w8U+lclCEKswqvSlI1J9PPMioc4BzPrI0SDTL2Fx86Ob8g6Uq2Jg43B+4bKNuyTa3OJvBDJMYxSwTxSWAkAlsWdgAx0+GBlGMdP5kRCctkHuHGzj3P/d6y7DOU/hw2vp6Jd/E6+S6uD65bPWdc67DVVUqlUoBAJDkpEoPAAD0b2IDAEglNgCAVGIDAEglNgCAVGIDAEglNgCAVGIDAEglNgCAVGIDAEhVVmycc845UVVVddRqbm7Omg8A6OOqy9l5w4YNcejQoa7tv/71r3H99dfHzJkzj/tgAED/UHUsD2KbM2dO/Pa3v43t27dHVVXV8ZwLAOgnynpn44MOHDgQTzzxRLS0tHxkaBSLxSgWi13bhw8fjn/84x8xfPhwgQIAfUSpVIq9e/dGXV1dnHRSmbd8lnro6aefLg0YMKD09ttvf+R+8+bNK0WEZVmWZVn9YLW3t5fdDD3+GOULX/hCDBw4MJ5//vmP3O/f39no6OiIUaNGxf/68bfipFMG9eTUH8uounfT/ux/ufqs19LP0Z2bhvx3xc7NifHs3gkVO/eLf/+v9HPs/J9PpZ9jYPvA9HMMae/Rf0LLMvTN4n/eKcnA1/6nYufmxDjwX3Ufa7+DB4vxx/U/iD179kRNTU1Z5+jRxyhvvfVWrF69Op555pn/uG+hUIhCoXDU6yedMig1NqoHH33O423QaSenn6M7pw3xreX+blCpcj9f1fvyf38yf///ZcCg/NgYMDA/NqqrK/eRc/VJ+deQyjpcXd7vYk9ugejR/7EWLVoUtbW1MX369J4cDgB8gpQdG4cPH45FixZFU1NTVFf3+P5SAOATouzYWL16dezcuTPuuOOOjHkAgH6m7Lcmpk2bFj28pxQA+ARylyEAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpyo6Nt99+O2699dYYPnx4nHrqqXHppZfGpk2bMmYDAPqB6nJ2fv/992Pq1Klx9dVXx/Lly6O2tjbeeOONOP3005PGAwD6urJi4wc/+EHU19fHokWLul4755xzPvKYYrEYxWKxa7uzs7O8CQGAPq2sj1GWLVsWDQ0NMXPmzKitrY2JEyfGwoULP/KY1tbWqKmp6Vr19fXHNDAA0LeUFRs7duyItra2OP/882PlypUxa9asuOeee2Lx4sXdHjN37tzo6OjoWu3t7cc8NADQd5T1Mcrhw4ejoaEh5s+fHxEREydOjK1bt0ZbW1vcfvvtH3pMoVCIQqFw7JMCAH1SWe9sjBgxIi6++OIjXrvoooti586dx3UoAKD/KCs2pk6dGtu2bTvitddeey1Gjx59XIcCAPqPsmLj3nvvjfXr18f8+fPj9ddfjyeffDIWLFgQzc3NWfMBAH1cWbFx2WWXxdKlS+Opp56KsWPHxne/+914+OGH45ZbbsmaDwDo48q6QTQi4oYbbogbbrghYxYAoB/ybBQAIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSlRUbDz74YFRVVR2xzj777KzZAIB+oLrcAy655JJYvXp11/aAAQOO60AAQP9SdmxUV1eX9W5GsViMYrHYtd3Z2VnuKQGAPqzseza2b98edXV1MWbMmLj55ptjx44dH7l/a2tr1NTUdK36+voeDwsA9D1lxcbll18eixcvjpUrV8bChQtj165dMWXKlHjvvfe6PWbu3LnR0dHRtdrb2495aACg7yjrY5TGxsaufx43blxMnjw5zj333Hj88cejpaXlQ48pFApRKBSObUoAoM86pq++Dh48OMaNGxfbt28/XvMAAP3MMcVGsViMV199NUaMGHG85gEA+pmyYuO+++6LtWvXxptvvhl/+tOf4ktf+lJ0dnZGU1NT1nwAQB9X1j0bf/vb3+KrX/1qvPvuu3HmmWfGZz/72Vi/fn2MHj06az4AoI8rKzaWLFmSNQcA0E95NgoAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkOqYYqO1tTWqqqpizpw5x2kcAKC/6XFsbNiwIRYsWBDjx48/nvMAAP1Mj2Jj3759ccstt8TChQvjjDPOON4zAQD9SI9io7m5OaZPnx7XXXfdf9y3WCxGZ2fnEQsA+OSoLveAJUuWxEsvvRQbNmz4WPu3trbGd77znbIHAwD6h7Le2Whvb4/Zs2fHE088EYMGDfpYx8ydOzc6Ojq6Vnt7e48GBQD6prLe2di0aVPs3r07Jk2a1PXaoUOHYt26dfHII49EsViMAQMGHHFMoVCIQqFwfKYFAPqcsmLj2muvjS1bthzx2te//vW48MIL4/777z8qNAAAyoqNIUOGxNixY494bfDgwTF8+PCjXgcAiPA3iAIAycr+Nsq/W7NmzXEYAwDor7yzAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQKqyYqOtrS3Gjx8fQ4cOjaFDh8bkyZNj+fLlWbMBAP1AWbExcuTIeOihh2Ljxo2xcePGuOaaa+LGG2+MrVu3Zs0HAPRx1eXsPGPGjCO2v//970dbW1usX78+LrnkkuM6GADQP5QVGx906NCh+PWvfx379++PyZMnd7tfsViMYrHYtd3Z2dnTUwIAfVDZN4hu2bIlTjvttCgUCjFr1qxYunRpXHzxxd3u39raGjU1NV2rvr7+mAYGAPqWsmPjggsuiM2bN8f69evjm9/8ZjQ1NcUrr7zS7f5z586Njo6OrtXe3n5MAwMAfUvZH6MMHDgwzjvvvIiIaGhoiA0bNsRPf/rT+MUvfvGh+xcKhSgUCsc2JQDQZx3z37NRKpWOuCcDAOCDynpn44EHHojGxsaor6+PvXv3xpIlS2LNmjWxYsWKrPkAgD6urNj4+9//Hrfddlu88847UVNTE+PHj48VK1bE9ddfnzUfANDHlRUbjz32WNYcAEA/5dkoAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApCorNlpbW+Oyyy6LIUOGRG1tbdx0002xbdu2rNkAgH6grNhYu3ZtNDc3x/r162PVqlVx8ODBmDZtWuzfvz9rPgCgj6suZ+cVK1Ycsb1o0aKora2NTZs2xZVXXnlcBwMA+oeyYuPfdXR0RETEsGHDut2nWCxGsVjs2u7s7DyWUwIAfUyPbxAtlUrR0tISV1xxRYwdO7bb/VpbW6OmpqZr1dfX9/SUAEAf1OPYuOuuu+Lll1+Op5566iP3mzt3bnR0dHSt9vb2np4SAOiDevQxyt133x3Lli2LdevWxciRIz9y30KhEIVCoUfDAQB9X1mxUSqV4u67746lS5fGmjVrYsyYMVlzAQD9RFmx0dzcHE8++WQ899xzMWTIkNi1a1dERNTU1MQpp5ySMiAA0LeVdc9GW1tbdHR0xFVXXRUjRozoWk8//XTWfABAH1f2xygAAOXwbBQAIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSlR0b69atixkzZkRdXV1UVVXFs88+mzAWANBflB0b+/fvjwkTJsQjjzySMQ8A0M9Ul3tAY2NjNDY2fuz9i8ViFIvFru3Ozs5yTwkA9GHp92y0trZGTU1N16qvr88+JQDQi6THxty5c6Ojo6Nrtbe3Z58SAOhFyv4YpVyFQiEKhUL2aQCAXspXXwGAVGIDAEhV9sco+/bti9dff71r+80334zNmzfHsGHDYtSoUcd1OACg7ys7NjZu3BhXX31113ZLS0tERDQ1NcWvfvWr4zYYANA/lB0bV111VZRKpYxZAIB+yD0bAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApOpRbDz66KMxZsyYGDRoUEyaNCl+//vfH++5AIB+ouzYePrpp2POnDnx7W9/O/7yl7/E5z73uWhsbIydO3dmzAcA9HHV5R7w4x//OL7xjW/EnXfeGRERDz/8cKxcuTLa2tqitbX1qP2LxWIUi8Wu7Y6OjoiIOPy//9nTmT+Wg/uL/3mnY/TPff8n/Rzd2Vd1uGLn5sSo5M/Xifj9yf5vQETEoX/m/54cOlBKP8fBg/n/Prpz0uEDFTs3J8bBgx/vd/FfP4elUg9+5ktlKBaLpQEDBpSeeeaZI16/5557SldeeeWHHjNv3rxSRFiWZVmW1Q/WG2+8UU46lEqlUqmsdzbefffdOHToUJx11llHvH7WWWfFrl27PvSYuXPnRktLS9f2nj17YvTo0bFz586oqakp5/T9WmdnZ9TX10d7e3sMHTq00uP0Kq5N91yb7rk23XNtuufadK+joyNGjRoVw4YNK/vYsj9GiYioqqo6YrtUKh312r8UCoUoFApHvV5TU+Nf5IcYOnSo69IN16Z7rk33XJvuuTbdc226d9JJ5X+3pKwjPvWpT8WAAQOOehdj9+7dR73bAQAQUWZsDBw4MCZNmhSrVq064vVVq1bFlClTjutgAED/UPbHKC0tLXHbbbdFQ0NDTJ48ORYsWBA7d+6MWbNmfazjC4VCzJs370M/Wvkkc12659p0z7XpnmvTPdeme65N947l2lSVevAdlkcffTR++MMfxjvvvBNjx46Nn/zkJ3HllVeWfXIAoP/rUWwAAHxcno0CAKQSGwBAKrEBAKQSGwBAqhMaGx5N/+HWrVsXM2bMiLq6uqiqqopnn3220iP1Cq2trXHZZZfFkCFDora2Nm666abYtm1bpcfqFdra2mL8+PFdf8vh5MmTY/ny5ZUeq9dpbW2NqqqqmDNnTqVH6RUefPDBqKqqOmKdffbZlR6rV3j77bfj1ltvjeHDh8epp54al156aWzatKnSY1XcOeecc9TPTFVVVTQ3N5f155yw2PBo+u7t378/JkyYEI888kilR+lV1q5dG83NzbF+/fpYtWpVHDx4MKZNmxb79++v9GgVN3LkyHjooYdi48aNsXHjxrjmmmvixhtvjK1bt1Z6tF5jw4YNsWDBghg/fnylR+lVLrnkknjnnXe61pYtWyo9UsW9//77MXXq1Dj55JNj+fLl8corr8SPfvSjOP300ys9WsVt2LDhiJ+Xf/2lnjNnzizvDyr70W099JnPfKY0a9asI1678MILS9/61rdO1Ah9QkSUli5dWukxeqXdu3eXIqK0du3aSo/SK51xxhmlX/7yl5Ueo1fYu3dv6fzzzy+tWrWq9PnPf740e/bsSo/UK8ybN680YcKESo/R69x///2lK664otJj9AmzZ88unXvuuaXDhw+XddwJeWfjwIEDsWnTppg2bdoRr0+bNi3++Mc/nogR6Ac6OjoiInr0xMH+7NChQ7FkyZLYv39/TJ48udLj9ArNzc0xffr0uO666yo9Sq+zffv2qKurizFjxsTNN98cO3bsqPRIFbds2bJoaGiImTNnRm1tbUycODEWLlxY6bF6nQMHDsQTTzwRd9xxR7cPX+3OCYmNnjyaHj6oVCpFS0tLXHHFFTF27NhKj9MrbNmyJU477bQoFAoxa9asWLp0aVx88cWVHqvilixZEi+99FK0trZWepRe5/LLL4/FixfHypUrY+HChbFr166YMmVKvPfee5UeraJ27NgRbW1tcf7558fKlStj1qxZcc8998TixYsrPVqv8uyzz8aePXvia1/7WtnH9ugR8z1VzqPp4YPuuuuuePnll+MPf/hDpUfpNS644ILYvHlz7NmzJ37zm99EU1NTrF279hMdHO3t7TF79ux44YUXYtCgQZUep9dpbGzs+udx48bF5MmT49xzz43HH388WlpaKjhZZR0+fDgaGhpi/vz5ERExceLE2Lp1a7S1tcXtt99e4el6j8ceeywaGxujrq6u7GNPyDsbHk3Psbj77rtj2bJl8eKLL8bIkSMrPU6vMXDgwDjvvPOioaEhWltbY8KECfHTn/600mNV1KZNm2L37t0xadKkqK6ujurq6li7dm387Gc/i+rq6jh06FClR+xVBg8eHOPGjYvt27dXepSKGjFixFGRftFFF/kCwwe89dZbsXr16rjzzjt7dPwJiQ2PpqcnSqVS3HXXXfHMM8/E7373uxgzZkylR+rVSqVSFIvFSo9RUddee21s2bIlNm/e3LUaGhrilltuic2bN8eAAQMqPWKvUiwW49VXX40RI0ZUepSKmjp16lFfq3/ttddi9OjRFZqo91m0aFHU1tbG9OnTe3T8CfsY5VgfTd+f7du3L15//fWu7TfffDM2b94cw4YNi1GjRlVwsspqbm6OJ598Mp577rkYMmRI1ztjNTU1ccopp1R4usp64IEHorGxMerr62Pv3r2xZMmSWLNmTaxYsaLSo1XUkCFDjrqnZ/DgwTF8+HD3+kTEfffdFzNmzIhRo0bF7t2743vf+150dnZGU1NTpUerqHvvvTemTJkS8+fPjy9/+cvx5z//ORYsWBALFiyo9Gi9wuHDh2PRokXR1NQU1dU9zIaEb8Z06+c//3lp9OjRpYEDB5Y+/elP+wrj//Piiy+WIuKo1dTUVOnRKurDrklElBYtWlTp0Srujjvu6PpdOvPMM0vXXntt6YUXXqj0WL2Sr77+f1/5yldKI0aMKJ188smlurq60he/+MXS1q1bKz1Wr/D888+Xxo4dWyoUCqULL7ywtGDBgkqP1GusXLmyFBGlbdu29fjP8Ih5ACCVZ6MAAKnEBgCQSmwAAKnEBgCQSmwAAKnEBgCQSmwAAKnEBgCQSmwAAKnEBgCQSmwAAKn+L7k4WqKTQDMIAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.contourf(u[:, :, 0].transpose())" + ] + }, + { + "cell_type": "code", + "execution_count": 67, + "id": "84b291bd", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 67, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhsAAAGiCAYAAABOCgSdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAX5UlEQVR4nO3dfWyV9f3w8U+lclCEKswqvSlI1J9PPMioc4BzPrI0SDTL2Fx86Ob8g6Uq2Jg43B+4bKNuyTa3OJvBDJMYxSwTxSWAkAlsWdgAx0+GBlGMdP5kRCctkHuHGzj3P/d6y7DOU/hw2vp6Jd/E6+S6uD65bPWdc67DVVUqlUoBAJDkpEoPAAD0b2IDAEglNgCAVGIDAEglNgCAVGIDAEglNgCAVGIDAEglNgCAVGIDAEhVVmycc845UVVVddRqbm7Omg8A6OOqy9l5w4YNcejQoa7tv/71r3H99dfHzJkzj/tgAED/UHUsD2KbM2dO/Pa3v43t27dHVVXV8ZwLAOgnynpn44MOHDgQTzzxRLS0tHxkaBSLxSgWi13bhw8fjn/84x8xfPhwgQIAfUSpVIq9e/dGXV1dnHRSmbd8lnro6aefLg0YMKD09ttvf+R+8+bNK0WEZVmWZVn9YLW3t5fdDD3+GOULX/hCDBw4MJ5//vmP3O/f39no6OiIUaNGxf/68bfipFMG9eTUH8uounfT/ux/ufqs19LP0Z2bhvx3xc7NifHs3gkVO/eLf/+v9HPs/J9PpZ9jYPvA9HMMae/Rf0LLMvTN4n/eKcnA1/6nYufmxDjwX3Ufa7+DB4vxx/U/iD179kRNTU1Z5+jRxyhvvfVWrF69Op555pn/uG+hUIhCoXDU6yedMig1NqoHH33O423QaSenn6M7pw3xreX+blCpcj9f1fvyf38yf///ZcCg/NgYMDA/NqqrK/eRc/VJ+deQyjpcXd7vYk9ugejR/7EWLVoUtbW1MX369J4cDgB8gpQdG4cPH45FixZFU1NTVFf3+P5SAOATouzYWL16dezcuTPuuOOOjHkAgH6m7Lcmpk2bFj28pxQA+ARylyEAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpyo6Nt99+O2699dYYPnx4nHrqqXHppZfGpk2bMmYDAPqB6nJ2fv/992Pq1Klx9dVXx/Lly6O2tjbeeOONOP3005PGAwD6urJi4wc/+EHU19fHokWLul4755xzPvKYYrEYxWKxa7uzs7O8CQGAPq2sj1GWLVsWDQ0NMXPmzKitrY2JEyfGwoULP/KY1tbWqKmp6Vr19fXHNDAA0LeUFRs7duyItra2OP/882PlypUxa9asuOeee2Lx4sXdHjN37tzo6OjoWu3t7cc8NADQd5T1Mcrhw4ejoaEh5s+fHxEREydOjK1bt0ZbW1vcfvvtH3pMoVCIQqFw7JMCAH1SWe9sjBgxIi6++OIjXrvoooti586dx3UoAKD/KCs2pk6dGtu2bTvitddeey1Gjx59XIcCAPqPsmLj3nvvjfXr18f8+fPj9ddfjyeffDIWLFgQzc3NWfMBAH1cWbFx2WWXxdKlS+Opp56KsWPHxne/+914+OGH45ZbbsmaDwDo48q6QTQi4oYbbogbbrghYxYAoB/ybBQAIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSlRUbDz74YFRVVR2xzj777KzZAIB+oLrcAy655JJYvXp11/aAAQOO60AAQP9SdmxUV1eX9W5GsViMYrHYtd3Z2VnuKQGAPqzseza2b98edXV1MWbMmLj55ptjx44dH7l/a2tr1NTUdK36+voeDwsA9D1lxcbll18eixcvjpUrV8bChQtj165dMWXKlHjvvfe6PWbu3LnR0dHRtdrb2495aACg7yjrY5TGxsaufx43blxMnjw5zj333Hj88cejpaXlQ48pFApRKBSObUoAoM86pq++Dh48OMaNGxfbt28/XvMAAP3MMcVGsViMV199NUaMGHG85gEA+pmyYuO+++6LtWvXxptvvhl/+tOf4ktf+lJ0dnZGU1NT1nwAQB9X1j0bf/vb3+KrX/1qvPvuu3HmmWfGZz/72Vi/fn2MHj06az4AoI8rKzaWLFmSNQcA0E95NgoAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkEpsAACpxAYAkOqYYqO1tTWqqqpizpw5x2kcAKC/6XFsbNiwIRYsWBDjx48/nvMAAP1Mj2Jj3759ccstt8TChQvjjDPOON4zAQD9SI9io7m5OaZPnx7XXXfdf9y3WCxGZ2fnEQsA+OSoLveAJUuWxEsvvRQbNmz4WPu3trbGd77znbIHAwD6h7Le2Whvb4/Zs2fHE088EYMGDfpYx8ydOzc6Ojq6Vnt7e48GBQD6prLe2di0aVPs3r07Jk2a1PXaoUOHYt26dfHII49EsViMAQMGHHFMoVCIQqFwfKYFAPqcsmLj2muvjS1bthzx2te//vW48MIL4/777z8qNAAAyoqNIUOGxNixY494bfDgwTF8+PCjXgcAiPA3iAIAycr+Nsq/W7NmzXEYAwDor7yzAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQCqxAQCkEhsAQKqyYqOtrS3Gjx8fQ4cOjaFDh8bkyZNj+fLlWbMBAP1AWbExcuTIeOihh2Ljxo2xcePGuOaaa+LGG2+MrVu3Zs0HAPRx1eXsPGPGjCO2v//970dbW1usX78+LrnkkuM6GADQP5QVGx906NCh+PWvfx379++PyZMnd7tfsViMYrHYtd3Z2dnTUwIAfVDZN4hu2bIlTjvttCgUCjFr1qxYunRpXHzxxd3u39raGjU1NV2rvr7+mAYGAPqWsmPjggsuiM2bN8f69evjm9/8ZjQ1NcUrr7zS7f5z586Njo6OrtXe3n5MAwMAfUvZH6MMHDgwzjvvvIiIaGhoiA0bNsRPf/rT+MUvfvGh+xcKhSgUCsc2JQDQZx3z37NRKpWOuCcDAOCDynpn44EHHojGxsaor6+PvXv3xpIlS2LNmjWxYsWKrPkAgD6urNj4+9//Hrfddlu88847UVNTE+PHj48VK1bE9ddfnzUfANDHlRUbjz32WNYcAEA/5dkoAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApCorNlpbW+Oyyy6LIUOGRG1tbdx0002xbdu2rNkAgH6grNhYu3ZtNDc3x/r162PVqlVx8ODBmDZtWuzfvz9rPgCgj6suZ+cVK1Ycsb1o0aKora2NTZs2xZVXXnlcBwMA+oeyYuPfdXR0RETEsGHDut2nWCxGsVjs2u7s7DyWUwIAfUyPbxAtlUrR0tISV1xxRYwdO7bb/VpbW6OmpqZr1dfX9/SUAEAf1OPYuOuuu+Lll1+Op5566iP3mzt3bnR0dHSt9vb2np4SAOiDevQxyt133x3Lli2LdevWxciRIz9y30KhEIVCoUfDAQB9X1mxUSqV4u67746lS5fGmjVrYsyYMVlzAQD9RFmx0dzcHE8++WQ899xzMWTIkNi1a1dERNTU1MQpp5ySMiAA0LeVdc9GW1tbdHR0xFVXXRUjRozoWk8//XTWfABAH1f2xygAAOXwbBQAIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSiQ0AIJXYAABSlR0b69atixkzZkRdXV1UVVXFs88+mzAWANBflB0b+/fvjwkTJsQjjzySMQ8A0M9Ul3tAY2NjNDY2fuz9i8ViFIvFru3Ozs5yTwkA9GHp92y0trZGTU1N16qvr88+JQDQi6THxty5c6Ojo6Nrtbe3Z58SAOhFyv4YpVyFQiEKhUL2aQCAXspXXwGAVGIDAEhV9sco+/bti9dff71r+80334zNmzfHsGHDYtSoUcd1OACg7ys7NjZu3BhXX31113ZLS0tERDQ1NcWvfvWr4zYYANA/lB0bV111VZRKpYxZAIB+yD0bAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApBIbAEAqsQEApOpRbDz66KMxZsyYGDRoUEyaNCl+//vfH++5AIB+ouzYePrpp2POnDnx7W9/O/7yl7/E5z73uWhsbIydO3dmzAcA9HHV5R7w4x//OL7xjW/EnXfeGRERDz/8cKxcuTLa2tqitbX1qP2LxWIUi8Wu7Y6OjoiIOPy//9nTmT+Wg/uL/3mnY/TPff8n/Rzd2Vd1uGLn5sSo5M/Xifj9yf5vQETEoX/m/54cOlBKP8fBg/n/Prpz0uEDFTs3J8bBgx/vd/FfP4elUg9+5ktlKBaLpQEDBpSeeeaZI16/5557SldeeeWHHjNv3rxSRFiWZVmW1Q/WG2+8UU46lEqlUqmsdzbefffdOHToUJx11llHvH7WWWfFrl27PvSYuXPnRktLS9f2nj17YvTo0bFz586oqakp5/T9WmdnZ9TX10d7e3sMHTq00uP0Kq5N91yb7rk23XNtuufadK+joyNGjRoVw4YNK/vYsj9GiYioqqo6YrtUKh312r8UCoUoFApHvV5TU+Nf5IcYOnSo69IN16Z7rk33XJvuuTbdc226d9JJ5X+3pKwjPvWpT8WAAQOOehdj9+7dR73bAQAQUWZsDBw4MCZNmhSrVq064vVVq1bFlClTjutgAED/UPbHKC0tLXHbbbdFQ0NDTJ48ORYsWBA7d+6MWbNmfazjC4VCzJs370M/Wvkkc12659p0z7XpnmvTPdeme65N947l2lSVevAdlkcffTR++MMfxjvvvBNjx46Nn/zkJ3HllVeWfXIAoP/rUWwAAHxcno0CAKQSGwBAKrEBAKQSGwBAqhMaGx5N/+HWrVsXM2bMiLq6uqiqqopnn3220iP1Cq2trXHZZZfFkCFDora2Nm666abYtm1bpcfqFdra2mL8+PFdf8vh5MmTY/ny5ZUeq9dpbW2NqqqqmDNnTqVH6RUefPDBqKqqOmKdffbZlR6rV3j77bfj1ltvjeHDh8epp54al156aWzatKnSY1XcOeecc9TPTFVVVTQ3N5f155yw2PBo+u7t378/JkyYEI888kilR+lV1q5dG83NzbF+/fpYtWpVHDx4MKZNmxb79++v9GgVN3LkyHjooYdi48aNsXHjxrjmmmvixhtvjK1bt1Z6tF5jw4YNsWDBghg/fnylR+lVLrnkknjnnXe61pYtWyo9UsW9//77MXXq1Dj55JNj+fLl8corr8SPfvSjOP300ys9WsVt2LDhiJ+Xf/2lnjNnzizvDyr70W099JnPfKY0a9asI1678MILS9/61rdO1Ah9QkSUli5dWukxeqXdu3eXIqK0du3aSo/SK51xxhmlX/7yl5Ueo1fYu3dv6fzzzy+tWrWq9PnPf740e/bsSo/UK8ybN680YcKESo/R69x///2lK664otJj9AmzZ88unXvuuaXDhw+XddwJeWfjwIEDsWnTppg2bdoRr0+bNi3++Mc/nogR6Ac6OjoiInr0xMH+7NChQ7FkyZLYv39/TJ48udLj9ArNzc0xffr0uO666yo9Sq+zffv2qKurizFjxsTNN98cO3bsqPRIFbds2bJoaGiImTNnRm1tbUycODEWLlxY6bF6nQMHDsQTTzwRd9xxR7cPX+3OCYmNnjyaHj6oVCpFS0tLXHHFFTF27NhKj9MrbNmyJU477bQoFAoxa9asWLp0aVx88cWVHqvilixZEi+99FK0trZWepRe5/LLL4/FixfHypUrY+HChbFr166YMmVKvPfee5UeraJ27NgRbW1tcf7558fKlStj1qxZcc8998TixYsrPVqv8uyzz8aePXvia1/7WtnH9ugR8z1VzqPp4YPuuuuuePnll+MPf/hDpUfpNS644ILYvHlz7NmzJ37zm99EU1NTrF279hMdHO3t7TF79ux44YUXYtCgQZUep9dpbGzs+udx48bF5MmT49xzz43HH388WlpaKjhZZR0+fDgaGhpi/vz5ERExceLE2Lp1a7S1tcXtt99e4el6j8ceeywaGxujrq6u7GNPyDsbHk3Psbj77rtj2bJl8eKLL8bIkSMrPU6vMXDgwDjvvPOioaEhWltbY8KECfHTn/600mNV1KZNm2L37t0xadKkqK6ujurq6li7dm387Gc/i+rq6jh06FClR+xVBg8eHOPGjYvt27dXepSKGjFixFGRftFFF/kCwwe89dZbsXr16rjzzjt7dPwJiQ2PpqcnSqVS3HXXXfHMM8/E7373uxgzZkylR+rVSqVSFIvFSo9RUddee21s2bIlNm/e3LUaGhrilltuic2bN8eAAQMqPWKvUiwW49VXX40RI0ZUepSKmjp16lFfq3/ttddi9OjRFZqo91m0aFHU1tbG9OnTe3T8CfsY5VgfTd+f7du3L15//fWu7TfffDM2b94cw4YNi1GjRlVwsspqbm6OJ598Mp577rkYMmRI1ztjNTU1ccopp1R4usp64IEHorGxMerr62Pv3r2xZMmSWLNmTaxYsaLSo1XUkCFDjrqnZ/DgwTF8+HD3+kTEfffdFzNmzIhRo0bF7t2743vf+150dnZGU1NTpUerqHvvvTemTJkS8+fPjy9/+cvx5z//ORYsWBALFiyo9Gi9wuHDh2PRokXR1NQU1dU9zIaEb8Z06+c//3lp9OjRpYEDB5Y+/elP+wrj//Piiy+WIuKo1dTUVOnRKurDrklElBYtWlTp0Srujjvu6PpdOvPMM0vXXntt6YUXXqj0WL2Sr77+f1/5yldKI0aMKJ188smlurq60he/+MXS1q1bKz1Wr/D888+Xxo4dWyoUCqULL7ywtGDBgkqP1GusXLmyFBGlbdu29fjP8Ih5ACCVZ6MAAKnEBgCQSmwAAKnEBgCQSmwAAKnEBgCQSmwAAKnEBgCQSmwAAKnEBgCQSmwAAKn+L7k4WqKTQDMIAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.contourf(uback[:, :, 0].transpose())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eae3f5fc", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pluto/fft_debug b/pluto/fft_debug new file mode 100644 index 0000000..f5fe187 --- /dev/null +++ b/pluto/fft_debug @@ -0,0 +1,1377 @@ +### A Pluto.jl notebook ### +# v0.19.22 + +using Markdown +using InteractiveUtils + +# ╔═╡ 74dccc96-e792-11ed-2b5a-b305a053025e +begin + using FFTW + using LazyGrids +end + +# ╔═╡ 13168eeb-afa0-4225-8e44-53b4d1b09874 +begin + using CairoMakie +end + +# ╔═╡ 3101cf7e-77fa-404f-a871-8ed4423885ff +function ingredients(path::String) + # this is from the Julia source code (evalfile in base/loading.jl) + # but with the modification that it returns the module instead of the last object + name = Symbol(basename(path)) + m = Module(name) + Core.eval(m, + Expr(:toplevel, + :(eval(x) = $(Expr(:core, :eval))($name, x)), + :(include(x) = $(Expr(:top, :include))($name, x)), + :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)), + :(include($path)))) + m +end + +# ╔═╡ c68bb44c-f56d-486b-92b2-55adc7c0517f +sgpu = ingredients("../src/spectralGPU.jl").spectralGPU + +# ╔═╡ 9c077253-7ece-464c-9c85-89ca1c7cfb9b +names(sgpu) + +# ╔═╡ b7efeb89-8de4-4f61-a2d2-d46364864e9e +const par = sgpu.markers.SingleThread() + +# ╔═╡ 3b36b955-c489-4673-b8c3-56998354fae1 +const N = 128 + +# ╔═╡ ba5d6ac1-727f-4be7-acd5-8c098d93d7ba +const K = sgpu.mesh.wavenumbers(N) + +# ╔═╡ 046c6800-08f9-4e69-9570-0c802ce62423 +const mesh = sgpu.mesh.new_mesh(N) + +# ╔═╡ 5421714b-baf3-4136-9236-cb17ff7be80a +u = sin.(mesh.x); + +# ╔═╡ 85011187-8daa-4064-b93c-4624fa3dfcca +u |> sum + +# ╔═╡ 59144084-5b0e-4d98-b3ca-52dafc657f0a +typeof(u) + +# ╔═╡ ef502d5e-cfca-425b-96bc-385cc8dbd18b +contourf(u[:, :, 1]) + +# ╔═╡ 92782cc9-bc97-47df-a173-8be22bca7ddd +begin + uhat = ComplexF64.(zeros(K.kn, N, N)); + uback = zeros(N, N, N); + sgpu.fft.fftn_mpi!(par, u, uhat); + sgpu.fft.ifftn_mpi!(par, K, uhat, uback); + + uhat |> sum +end + +# ╔═╡ 3ac91add-acdd-4ec7-91f1-084a35367877 +contourf(uback[:, :, 1]) + +# ╔═╡ 00000000-0000-0000-0000-000000000001 +PLUTO_PROJECT_TOML_CONTENTS = """ +[deps] +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +LazyGrids = "7031d0ef-c40d-4431-b2f8-61a8d2f650db" + +[compat] +CairoMakie = "~0.10.4" +FFTW = "~1.6.0" +LazyGrids = "~0.5.0" +""" + +# ╔═╡ 00000000-0000-0000-0000-000000000002 +PLUTO_MANIFEST_TOML_CONTENTS = """ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.8.5" +manifest_format = "2.0" +project_hash = "71a25de324d66428b531c134cc3204997729a01c" + +[[deps.AbstractFFTs]] +deps = ["ChainRulesCore", "LinearAlgebra"] +git-tree-sha1 = "16b6dbc4cf7caee4e1e75c49485ec67b667098a0" +uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c" +version = "1.3.1" + +[[deps.AbstractTrees]] +git-tree-sha1 = "faa260e4cb5aba097a73fab382dd4b5819d8ec8c" +uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" +version = "0.4.4" + +[[deps.Adapt]] +deps = ["LinearAlgebra", "Requires"] +git-tree-sha1 = "cc37d689f599e8df4f464b2fa3870ff7db7492ef" +uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +version = "3.6.1" + +[[deps.Animations]] +deps = ["Colors"] +git-tree-sha1 = "e81c509d2c8e49592413bfb0bb3b08150056c79d" +uuid = "27a7e980-b3e6-11e9-2bcd-0b925532e340" +version = "0.4.1" + +[[deps.ArgTools]] +uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" +version = "1.1.1" + +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" + +[[deps.Automa]] +deps = ["Printf", "ScanByte", "TranscodingStreams"] +git-tree-sha1 = "d50976f217489ce799e366d9561d56a98a30d7fe" +uuid = "67c07d97-cdcb-5c2c-af73-a7f9c32a568b" +version = "0.8.2" + +[[deps.AxisAlgorithms]] +deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"] +git-tree-sha1 = "66771c8d21c8ff5e3a93379480a2307ac36863f7" +uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950" +version = "1.0.1" + +[[deps.AxisArrays]] +deps = ["Dates", "IntervalSets", "IterTools", "RangeArrays"] +git-tree-sha1 = "1dd4d9f5beebac0c03446918741b1a03dc5e5788" +uuid = "39de3d68-74b9-583c-8d2d-e117c070f3a9" +version = "0.4.6" + +[[deps.Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[deps.Bzip2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2" +uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0" +version = "1.0.8+0" + +[[deps.CEnum]] +git-tree-sha1 = "eb4cb44a499229b3b8426dcfb5dd85333951ff90" +uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82" +version = "0.4.2" + +[[deps.CRC32c]] +uuid = "8bf52ea8-c179-5cab-976a-9e18b702a9bc" + +[[deps.Cairo]] +deps = ["Cairo_jll", "Colors", "Glib_jll", "Graphics", "Libdl", "Pango_jll"] +git-tree-sha1 = "d0b3f8b4ad16cb0a2988c6788646a5e6a17b6b1b" +uuid = "159f3aea-2a34-519c-b102-8c37f9878175" +version = "1.0.5" + +[[deps.CairoMakie]] +deps = ["Base64", "Cairo", "Colors", "FFTW", "FileIO", "FreeType", "GeometryBasics", "LinearAlgebra", "Makie", "SHA", "SnoopPrecompile"] +git-tree-sha1 = "2aba202861fd2b7603beb80496b6566491229855" +uuid = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +version = "0.10.4" + +[[deps.Cairo_jll]] +deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"] +git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2" +uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a" +version = "1.16.1+1" + +[[deps.Calculus]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad" +uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" +version = "0.5.1" + +[[deps.ChainRulesCore]] +deps = ["Compat", "LinearAlgebra", "SparseArrays"] +git-tree-sha1 = "c6d890a52d2c4d55d326439580c3b8d0875a77d9" +uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +version = "1.15.7" + +[[deps.ChangesOfVariables]] +deps = ["LinearAlgebra", "Test"] +git-tree-sha1 = "f84967c4497e0e1955f9a582c232b02847c5f589" +uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" +version = "0.1.7" + +[[deps.ColorBrewer]] +deps = ["Colors", "JSON", "Test"] +git-tree-sha1 = "61c5334f33d91e570e1d0c3eb5465835242582c4" +uuid = "a2cac450-b92f-5266-8821-25eda20663c8" +version = "0.4.0" + +[[deps.ColorSchemes]] +deps = ["ColorTypes", "ColorVectorSpace", "Colors", "FixedPointNumbers", "PrecompileTools", "Random"] +git-tree-sha1 = "be6ab11021cd29f0344d5c4357b163af05a48cba" +uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4" +version = "3.21.0" + +[[deps.ColorTypes]] +deps = ["FixedPointNumbers", "Random"] +git-tree-sha1 = "eb7f0f8307f71fac7c606984ea5fb2817275d6e4" +uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" +version = "0.11.4" + +[[deps.ColorVectorSpace]] +deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"] +git-tree-sha1 = "600cc5508d66b78aae350f7accdb58763ac18589" +uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4" +version = "0.9.10" + +[[deps.Colors]] +deps = ["ColorTypes", "FixedPointNumbers", "Reexport"] +git-tree-sha1 = "fc08e5930ee9a4e03f84bfb5211cb54e7769758a" +uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" +version = "0.12.10" + +[[deps.Compat]] +deps = ["Dates", "LinearAlgebra", "UUIDs"] +git-tree-sha1 = "7a60c856b9fa189eb34f5f8a6f6b5529b7942957" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "4.6.1" + +[[deps.CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "1.0.1+0" + +[[deps.ConstructionBase]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "89a9db8d28102b094992472d333674bd1a83ce2a" +uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" +version = "1.5.1" + +[[deps.Contour]] +git-tree-sha1 = "d05d9e7b7aedff4e5b51a029dced05cfb6125781" +uuid = "d38c429a-6771-53c6-b99e-75d170b6e991" +version = "0.6.2" + +[[deps.DataAPI]] +git-tree-sha1 = "e8119c1a33d267e16108be441a287a6981ba1630" +uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" +version = "1.14.0" + +[[deps.DataStructures]] +deps = ["Compat", "InteractiveUtils", "OrderedCollections"] +git-tree-sha1 = "d1fff3a548102f48987a52a2e0d114fa97d730f0" +uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +version = "0.18.13" + +[[deps.DataValueInterfaces]] +git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6" +uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464" +version = "1.0.0" + +[[deps.Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" + +[[deps.DensityInterface]] +deps = ["InverseFunctions", "Test"] +git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b" +uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d" +version = "0.4.0" + +[[deps.Distributed]] +deps = ["Random", "Serialization", "Sockets"] +uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" + +[[deps.Distributions]] +deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"] +git-tree-sha1 = "13027f188d26206b9e7b863036f87d2f2e7d013a" +uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" +version = "0.25.87" + +[[deps.DocStringExtensions]] +deps = ["LibGit2"] +git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d" +uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +version = "0.9.3" + +[[deps.Downloads]] +deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] +uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +version = "1.6.0" + +[[deps.DualNumbers]] +deps = ["Calculus", "NaNMath", "SpecialFunctions"] +git-tree-sha1 = "5837a837389fccf076445fce071c8ddaea35a566" +uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74" +version = "0.6.8" + +[[deps.EarCut_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "e3290f2d49e661fbd94046d7e3726ffcb2d41053" +uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5" +version = "2.2.4+0" + +[[deps.Expat_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "bad72f730e9e91c08d9427d5e8db95478a3c323d" +uuid = "2e619515-83b5-522b-bb60-26c02a35a201" +version = "2.4.8+0" + +[[deps.Extents]] +git-tree-sha1 = "5e1e4c53fa39afe63a7d356e30452249365fba99" +uuid = "411431e0-e8b7-467b-b5e0-f676ba4f2910" +version = "0.1.1" + +[[deps.FFMPEG]] +deps = ["FFMPEG_jll"] +git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8" +uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" +version = "0.4.1" + +[[deps.FFMPEG_jll]] +deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Pkg", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"] +git-tree-sha1 = "74faea50c1d007c85837327f6775bea60b5492dd" +uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5" +version = "4.4.2+2" + +[[deps.FFTW]] +deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"] +git-tree-sha1 = "f9818144ce7c8c41edf5c4c179c684d92aa4d9fe" +uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +version = "1.6.0" + +[[deps.FFTW_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "c6033cc3892d0ef5bb9cd29b7f2f0331ea5184ea" +uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a" +version = "3.3.10+0" + +[[deps.FileIO]] +deps = ["Pkg", "Requires", "UUIDs"] +git-tree-sha1 = "7be5f99f7d15578798f338f5433b6c432ea8037b" +uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" +version = "1.16.0" + +[[deps.FileWatching]] +uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" + +[[deps.FillArrays]] +deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"] +git-tree-sha1 = "fc86b4fd3eff76c3ce4f5e96e2fdfa6282722885" +uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" +version = "1.0.0" + +[[deps.FixedPointNumbers]] +deps = ["Statistics"] +git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc" +uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" +version = "0.8.4" + +[[deps.Fontconfig_jll]] +deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"] +git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03" +uuid = "a3f928ae-7b40-5064-980b-68af3947d34b" +version = "2.13.93+0" + +[[deps.Formatting]] +deps = ["Printf"] +git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8" +uuid = "59287772-0a20-5a39-b81b-1366585eb4c0" +version = "0.4.2" + +[[deps.FreeType]] +deps = ["CEnum", "FreeType2_jll"] +git-tree-sha1 = "cabd77ab6a6fdff49bfd24af2ebe76e6e018a2b4" +uuid = "b38be410-82b0-50bf-ab77-7b57e271db43" +version = "4.0.0" + +[[deps.FreeType2_jll]] +deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"] +git-tree-sha1 = "87eb71354d8ec1a96d4a7636bd57a7347dde3ef9" +uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7" +version = "2.10.4+0" + +[[deps.FreeTypeAbstraction]] +deps = ["ColorVectorSpace", "Colors", "FreeType", "GeometryBasics"] +git-tree-sha1 = "38a92e40157100e796690421e34a11c107205c86" +uuid = "663a7486-cb36-511b-a19d-713bb74d65c9" +version = "0.10.0" + +[[deps.FriBidi_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91" +uuid = "559328eb-81f9-559d-9380-de523a88c83c" +version = "1.0.10+0" + +[[deps.Future]] +deps = ["Random"] +uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" + +[[deps.GPUArraysCore]] +deps = ["Adapt"] +git-tree-sha1 = "1cd7f0af1aa58abc02ea1d872953a97359cb87fa" +uuid = "46192b85-c4d5-4398-a991-12ede77f4527" +version = "0.1.4" + +[[deps.GeoInterface]] +deps = ["Extents"] +git-tree-sha1 = "0eb6de0b312688f852f347171aba888658e29f20" +uuid = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" +version = "1.3.0" + +[[deps.GeometryBasics]] +deps = ["EarCut_jll", "GeoInterface", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"] +git-tree-sha1 = "659140c9375afa2f685e37c1a0b9c9a60ef56b40" +uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326" +version = "0.4.7" + +[[deps.Gettext_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"] +git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046" +uuid = "78b55507-aeef-58d4-861c-77aaff3498b1" +version = "0.21.0+0" + +[[deps.Glib_jll]] +deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Pkg", "Zlib_jll"] +git-tree-sha1 = "d3b3624125c1474292d0d8ed0f65554ac37ddb23" +uuid = "7746bdde-850d-59dc-9ae8-88ece973131d" +version = "2.74.0+2" + +[[deps.Graphics]] +deps = ["Colors", "LinearAlgebra", "NaNMath"] +git-tree-sha1 = "d61890399bc535850c4bf08e4e0d3a7ad0f21cbd" +uuid = "a2bd30eb-e257-5431-a919-1863eab51364" +version = "1.1.2" + +[[deps.Graphite2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011" +uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472" +version = "1.3.14+0" + +[[deps.GridLayoutBase]] +deps = ["GeometryBasics", "InteractiveUtils", "Observables"] +git-tree-sha1 = "678d136003ed5bceaab05cf64519e3f956ffa4ba" +uuid = "3955a311-db13-416c-9275-1d80ed98e5e9" +version = "0.9.1" + +[[deps.Grisu]] +git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2" +uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe" +version = "1.0.2" + +[[deps.HarfBuzz_jll]] +deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"] +git-tree-sha1 = "129acf094d168394e80ee1dc4bc06ec835e510a3" +uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566" +version = "2.8.1+1" + +[[deps.HypergeometricFunctions]] +deps = ["DualNumbers", "LinearAlgebra", "OpenLibm_jll", "SpecialFunctions"] +git-tree-sha1 = "432b5b03176f8182bd6841fbfc42c718506a2d5f" +uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a" +version = "0.3.15" + +[[deps.ImageAxes]] +deps = ["AxisArrays", "ImageBase", "ImageCore", "Reexport", "SimpleTraits"] +git-tree-sha1 = "c54b581a83008dc7f292e205f4c409ab5caa0f04" +uuid = "2803e5a7-5153-5ecf-9a86-9b4c37f5f5ac" +version = "0.6.10" + +[[deps.ImageBase]] +deps = ["ImageCore", "Reexport"] +git-tree-sha1 = "b51bb8cae22c66d0f6357e3bcb6363145ef20835" +uuid = "c817782e-172a-44cc-b673-b171935fbb9e" +version = "0.1.5" + +[[deps.ImageCore]] +deps = ["AbstractFFTs", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Graphics", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "Reexport"] +git-tree-sha1 = "acf614720ef026d38400b3817614c45882d75500" +uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534" +version = "0.9.4" + +[[deps.ImageIO]] +deps = ["FileIO", "IndirectArrays", "JpegTurbo", "LazyModules", "Netpbm", "OpenEXR", "PNGFiles", "QOI", "Sixel", "TiffImages", "UUIDs"] +git-tree-sha1 = "342f789fd041a55166764c351da1710db97ce0e0" +uuid = "82e4d734-157c-48bb-816b-45c225c6df19" +version = "0.6.6" + +[[deps.ImageMetadata]] +deps = ["AxisArrays", "ImageAxes", "ImageBase", "ImageCore"] +git-tree-sha1 = "36cbaebed194b292590cba2593da27b34763804a" +uuid = "bc367c6b-8a6b-528e-b4bd-a4b897500b49" +version = "0.9.8" + +[[deps.Imath_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "3d09a9f60edf77f8a4d99f9e015e8fbf9989605d" +uuid = "905a6f67-0a94-5f89-b386-d35d92009cd1" +version = "3.1.7+0" + +[[deps.IndirectArrays]] +git-tree-sha1 = "012e604e1c7458645cb8b436f8fba789a51b257f" +uuid = "9b13fd28-a010-5f03-acff-a1bbcff69959" +version = "1.0.0" + +[[deps.Inflate]] +git-tree-sha1 = "5cd07aab533df5170988219191dfad0519391428" +uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9" +version = "0.1.3" + +[[deps.IntelOpenMP_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "0cb9352ef2e01574eeebdb102948a58740dcaf83" +uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0" +version = "2023.1.0+0" + +[[deps.InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" + +[[deps.Interpolations]] +deps = ["Adapt", "AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"] +git-tree-sha1 = "721ec2cf720536ad005cb38f50dbba7b02419a15" +uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +version = "0.14.7" + +[[deps.IntervalSets]] +deps = ["Dates", "Random", "Statistics"] +git-tree-sha1 = "16c0cc91853084cb5f58a78bd209513900206ce6" +uuid = "8197267c-284f-5f27-9208-e0e47529a953" +version = "0.7.4" + +[[deps.InverseFunctions]] +deps = ["Test"] +git-tree-sha1 = "6667aadd1cdee2c6cd068128b3d226ebc4fb0c67" +uuid = "3587e190-3f89-42d0-90ee-14403ec27112" +version = "0.1.9" + +[[deps.IrrationalConstants]] +git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2" +uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" +version = "0.2.2" + +[[deps.Isoband]] +deps = ["isoband_jll"] +git-tree-sha1 = "f9b6d97355599074dc867318950adaa6f9946137" +uuid = "f1662d9f-8043-43de-a69a-05efc1cc6ff4" +version = "0.1.1" + +[[deps.IterTools]] +git-tree-sha1 = "fa6287a4469f5e048d763df38279ee729fbd44e5" +uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e" +version = "1.4.0" + +[[deps.IteratorInterfaceExtensions]] +git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856" +uuid = "82899510-4779-5014-852e-03e436cf321d" +version = "1.0.0" + +[[deps.JLLWrappers]] +deps = ["Preferences"] +git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.4.1" + +[[deps.JSON]] +deps = ["Dates", "Mmap", "Parsers", "Unicode"] +git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a" +uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +version = "0.21.4" + +[[deps.JpegTurbo]] +deps = ["CEnum", "FileIO", "ImageCore", "JpegTurbo_jll", "TOML"] +git-tree-sha1 = "106b6aa272f294ba47e96bd3acbabdc0407b5c60" +uuid = "b835a17e-a41a-41e7-81f0-2f016b05efe0" +version = "0.1.2" + +[[deps.JpegTurbo_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "6f2675ef130a300a112286de91973805fcc5ffbc" +uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8" +version = "2.1.91+0" + +[[deps.KernelDensity]] +deps = ["Distributions", "DocStringExtensions", "FFTW", "Interpolations", "StatsBase"] +git-tree-sha1 = "4a9513ad756e712177bd342ba6c022b515ed8d76" +uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b" +version = "0.6.6" + +[[deps.LAME_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c" +uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d" +version = "3.100.1+0" + +[[deps.LZO_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6" +uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac" +version = "2.10.1+0" + +[[deps.LaTeXStrings]] +git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996" +uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" +version = "1.3.0" + +[[deps.LazyArtifacts]] +deps = ["Artifacts", "Pkg"] +uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3" + +[[deps.LazyGrids]] +deps = ["Statistics"] +git-tree-sha1 = "f43d10fea7e448a60e92976bbd8bfbca7a6e5d09" +uuid = "7031d0ef-c40d-4431-b2f8-61a8d2f650db" +version = "0.5.0" + +[[deps.LazyModules]] +git-tree-sha1 = "a560dd966b386ac9ae60bdd3a3d3a326062d3c3e" +uuid = "8cdb02fc-e678-4876-92c5-9defec4f444e" +version = "0.3.1" + +[[deps.LibCURL]] +deps = ["LibCURL_jll", "MozillaCACerts_jll"] +uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" +version = "0.6.3" + +[[deps.LibCURL_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] +uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" +version = "7.84.0+0" + +[[deps.LibGit2]] +deps = ["Base64", "NetworkOptions", "Printf", "SHA"] +uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" + +[[deps.LibSSH2_jll]] +deps = ["Artifacts", "Libdl", "MbedTLS_jll"] +uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" +version = "1.10.2+0" + +[[deps.Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" + +[[deps.Libffi_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "0b4a5d71f3e5200a7dff793393e09dfc2d874290" +uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490" +version = "3.2.2+1" + +[[deps.Libgcrypt_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"] +git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae" +uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4" +version = "1.8.7+0" + +[[deps.Libgpg_error_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9" +uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8" +version = "1.42.0+0" + +[[deps.Libiconv_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "c7cb1f5d892775ba13767a87c7ada0b980ea0a71" +uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" +version = "1.16.1+2" + +[[deps.Libmount_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73" +uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9" +version = "2.35.0+0" + +[[deps.Libuuid_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066" +uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700" +version = "2.36.0+0" + +[[deps.LinearAlgebra]] +deps = ["Libdl", "libblastrampoline_jll"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + +[[deps.LogExpFunctions]] +deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"] +git-tree-sha1 = "0a1b7c2863e44523180fdb3146534e265a91870b" +uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" +version = "0.3.23" + +[[deps.Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" + +[[deps.MKL_jll]] +deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"] +git-tree-sha1 = "2ce8695e1e699b68702c03402672a69f54b8aca9" +uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7" +version = "2022.2.0+0" + +[[deps.MacroTools]] +deps = ["Markdown", "Random"] +git-tree-sha1 = "42324d08725e200c23d4dfb549e0d5d89dede2d2" +uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +version = "0.5.10" + +[[deps.Makie]] +deps = ["Animations", "Base64", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "Distributions", "DocStringExtensions", "Downloads", "FFMPEG", "FileIO", "FixedPointNumbers", "Formatting", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageIO", "InteractiveUtils", "IntervalSets", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MakieCore", "Markdown", "Match", "MathTeXEngine", "MiniQhull", "Observables", "OffsetArrays", "Packing", "PlotUtils", "PolygonOps", "Printf", "Random", "RelocatableFolders", "Setfield", "Showoff", "SignedDistanceFields", "SnoopPrecompile", "SparseArrays", "StableHashTraits", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun"] +git-tree-sha1 = "74657542dc85c3b72b8a5a9392d57713d8b7a999" +uuid = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" +version = "0.19.4" + +[[deps.MakieCore]] +deps = ["Observables"] +git-tree-sha1 = "9926529455a331ed73c19ff06d16906737a876ed" +uuid = "20f20a25-4f0e-4fdf-b5d1-57303727442b" +version = "0.6.3" + +[[deps.MappedArrays]] +git-tree-sha1 = "e8b359ef06ec72e8c030463fe02efe5527ee5142" +uuid = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900" +version = "0.4.1" + +[[deps.Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" + +[[deps.Match]] +git-tree-sha1 = "1d9bc5c1a6e7ee24effb93f175c9342f9154d97f" +uuid = "7eb4fadd-790c-5f42-8a69-bfa0b872bfbf" +version = "1.2.0" + +[[deps.MathTeXEngine]] +deps = ["AbstractTrees", "Automa", "DataStructures", "FreeTypeAbstraction", "GeometryBasics", "LaTeXStrings", "REPL", "RelocatableFolders", "Test", "UnicodeFun"] +git-tree-sha1 = "8f52dbaa1351ce4cb847d95568cb29e62a307d93" +uuid = "0a4f8689-d25c-4efe-a92b-7142dfc1aa53" +version = "0.5.6" + +[[deps.MbedTLS_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" +version = "2.28.0+0" + +[[deps.MiniQhull]] +deps = ["QhullMiniWrapper_jll"] +git-tree-sha1 = "9dc837d180ee49eeb7c8b77bb1c860452634b0d1" +uuid = "978d7f02-9e05-4691-894f-ae31a51d76ca" +version = "0.4.0" + +[[deps.Missings]] +deps = ["DataAPI"] +git-tree-sha1 = "f66bdc5de519e8f8ae43bdc598782d35a25b1272" +uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" +version = "1.1.0" + +[[deps.Mmap]] +uuid = "a63ad114-7e13-5084-954f-fe012c677804" + +[[deps.MosaicViews]] +deps = ["MappedArrays", "OffsetArrays", "PaddedViews", "StackViews"] +git-tree-sha1 = "7b86a5d4d70a9f5cdf2dacb3cbe6d251d1a61dbe" +uuid = "e94cdb99-869f-56ef-bcf0-1ae2bcbe0389" +version = "0.3.4" + +[[deps.MozillaCACerts_jll]] +uuid = "14a3606d-f60d-562e-9121-12d972cd8159" +version = "2022.2.1" + +[[deps.NaNMath]] +deps = ["OpenLibm_jll"] +git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4" +uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +version = "1.0.2" + +[[deps.Netpbm]] +deps = ["FileIO", "ImageCore", "ImageMetadata"] +git-tree-sha1 = "5ae7ca23e13855b3aba94550f26146c01d259267" +uuid = "f09324ee-3d7c-5217-9330-fc30815ba969" +version = "1.1.0" + +[[deps.NetworkOptions]] +uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" +version = "1.2.0" + +[[deps.Observables]] +git-tree-sha1 = "6862738f9796b3edc1c09d0890afce4eca9e7e93" +uuid = "510215fc-4207-5dde-b226-833fc4488ee2" +version = "0.5.4" + +[[deps.OffsetArrays]] +deps = ["Adapt"] +git-tree-sha1 = "82d7c9e310fe55aa54996e6f7f94674e2a38fcb4" +uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +version = "1.12.9" + +[[deps.Ogg_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "887579a3eb005446d514ab7aeac5d1d027658b8f" +uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051" +version = "1.3.5+1" + +[[deps.OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.20+0" + +[[deps.OpenEXR]] +deps = ["Colors", "FileIO", "OpenEXR_jll"] +git-tree-sha1 = "327f53360fdb54df7ecd01e96ef1983536d1e633" +uuid = "52e1d378-f018-4a11-a4be-720524705ac7" +version = "0.3.2" + +[[deps.OpenEXR_jll]] +deps = ["Artifacts", "Imath_jll", "JLLWrappers", "Libdl", "Zlib_jll"] +git-tree-sha1 = "a4ca623df1ae99d09bc9868b008262d0c0ac1e4f" +uuid = "18a262bb-aa17-5467-a713-aee519bc75cb" +version = "3.1.4+0" + +[[deps.OpenLibm_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "05823500-19ac-5b8b-9628-191a04bc5112" +version = "0.8.1+0" + +[[deps.OpenSSL_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "9ff31d101d987eb9d66bd8b176ac7c277beccd09" +uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" +version = "1.1.20+0" + +[[deps.OpenSpecFun_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" +uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" +version = "0.5.5+0" + +[[deps.Opus_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720" +uuid = "91d4177d-7536-5919-b921-800302f37372" +version = "1.3.2+0" + +[[deps.OrderedCollections]] +git-tree-sha1 = "d321bf2de576bf25ec4d3e4360faca399afca282" +uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +version = "1.6.0" + +[[deps.PCRE2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15" +version = "10.40.0+0" + +[[deps.PDMats]] +deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"] +git-tree-sha1 = "67eae2738d63117a196f497d7db789821bce61d1" +uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" +version = "0.11.17" + +[[deps.PNGFiles]] +deps = ["Base64", "CEnum", "ImageCore", "IndirectArrays", "OffsetArrays", "libpng_jll"] +git-tree-sha1 = "f809158b27eba0c18c269cf2a2be6ed751d3e81d" +uuid = "f57f5aa1-a3ce-4bc8-8ab9-96f992907883" +version = "0.3.17" + +[[deps.Packing]] +deps = ["GeometryBasics"] +git-tree-sha1 = "ec3edfe723df33528e085e632414499f26650501" +uuid = "19eb6ba3-879d-56ad-ad62-d5c202156566" +version = "0.5.0" + +[[deps.PaddedViews]] +deps = ["OffsetArrays"] +git-tree-sha1 = "0fac6313486baae819364c52b4f483450a9d793f" +uuid = "5432bcbf-9aad-5242-b902-cca2824c8663" +version = "0.5.12" + +[[deps.Pango_jll]] +deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "FriBidi_jll", "Glib_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "84a314e3926ba9ec66ac097e3635e270986b0f10" +uuid = "36c8627f-9965-5494-a995-c6b170f724f3" +version = "1.50.9+0" + +[[deps.Parsers]] +deps = ["Dates", "SnoopPrecompile"] +git-tree-sha1 = "478ac6c952fddd4399e71d4779797c538d0ff2bf" +uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" +version = "2.5.8" + +[[deps.Pixman_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29" +uuid = "30392449-352a-5448-841d-b1acce4e97dc" +version = "0.40.1+0" + +[[deps.Pkg]] +deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] +uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +version = "1.8.0" + +[[deps.PkgVersion]] +deps = ["Pkg"] +git-tree-sha1 = "f6cf8e7944e50901594838951729a1861e668cb8" +uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688" +version = "0.3.2" + +[[deps.PlotUtils]] +deps = ["ColorSchemes", "Colors", "Dates", "PrecompileTools", "Printf", "Random", "Reexport", "Statistics"] +git-tree-sha1 = "f92e1315dadf8c46561fb9396e525f7200cdc227" +uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043" +version = "1.3.5" + +[[deps.PolygonOps]] +git-tree-sha1 = "77b3d3605fc1cd0b42d95eba87dfcd2bf67d5ff6" +uuid = "647866c9-e3ac-4575-94e7-e3d426903924" +version = "0.1.2" + +[[deps.PrecompileTools]] +deps = ["Preferences"] +git-tree-sha1 = "2e47054ffe7d0a8872e977c0d09eb4b3d162ebde" +uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +version = "1.0.2" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.3.0" + +[[deps.Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" + +[[deps.ProgressMeter]] +deps = ["Distributed", "Printf"] +git-tree-sha1 = "d7a7aef8f8f2d537104f170139553b14dfe39fe9" +uuid = "92933f4c-e287-5a05-a399-4b506db050ca" +version = "1.7.2" + +[[deps.QOI]] +deps = ["ColorTypes", "FileIO", "FixedPointNumbers"] +git-tree-sha1 = "18e8f4d1426e965c7b532ddd260599e1510d26ce" +uuid = "4b34888f-f399-49d4-9bb3-47ed5cae4e65" +version = "1.0.0" + +[[deps.QhullMiniWrapper_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Qhull_jll"] +git-tree-sha1 = "607cf73c03f8a9f83b36db0b86a3a9c14179621f" +uuid = "460c41e3-6112-5d7f-b78c-b6823adb3f2d" +version = "1.0.0+1" + +[[deps.Qhull_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "238dd7e2cc577281976b9681702174850f8d4cbc" +uuid = "784f63db-0788-585a-bace-daefebcd302b" +version = "8.0.1001+0" + +[[deps.QuadGK]] +deps = ["DataStructures", "LinearAlgebra"] +git-tree-sha1 = "6ec7ac8412e83d57e313393220879ede1740f9ee" +uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +version = "2.8.2" + +[[deps.REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] +uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" + +[[deps.Random]] +deps = ["SHA", "Serialization"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[deps.RangeArrays]] +git-tree-sha1 = "b9039e93773ddcfc828f12aadf7115b4b4d225f5" +uuid = "b3c3ace0-ae52-54e7-9d0b-2c1406fd6b9d" +version = "0.3.2" + +[[deps.Ratios]] +deps = ["Requires"] +git-tree-sha1 = "6d7bb727e76147ba18eed998700998e17b8e4911" +uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439" +version = "0.4.4" + +[[deps.Reexport]] +git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" +uuid = "189a3867-3050-52da-a836-e630ba90ab69" +version = "1.2.2" + +[[deps.RelocatableFolders]] +deps = ["SHA", "Scratch"] +git-tree-sha1 = "90bc7a7c96410424509e4263e277e43250c05691" +uuid = "05181044-ff0b-4ac5-8273-598c1e38db00" +version = "1.0.0" + +[[deps.Requires]] +deps = ["UUIDs"] +git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" +uuid = "ae029012-a4dd-5104-9daa-d747884805df" +version = "1.3.0" + +[[deps.Rmath]] +deps = ["Random", "Rmath_jll"] +git-tree-sha1 = "f65dcb5fa46aee0cf9ed6274ccbd597adc49aa7b" +uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa" +version = "0.7.1" + +[[deps.Rmath_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "6ed52fdd3382cf21947b15e8870ac0ddbff736da" +uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f" +version = "0.4.0+0" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" + +[[deps.SIMD]] +deps = ["SnoopPrecompile"] +git-tree-sha1 = "8b20084a97b004588125caebf418d8cab9e393d1" +uuid = "fdea26ae-647d-5447-a871-4b548cad5224" +version = "3.4.4" + +[[deps.ScanByte]] +deps = ["Libdl", "SIMD"] +git-tree-sha1 = "2436b15f376005e8790e318329560dcc67188e84" +uuid = "7b38b023-a4d7-4c5e-8d43-3f3097f304eb" +version = "0.3.3" + +[[deps.Scratch]] +deps = ["Dates"] +git-tree-sha1 = "30449ee12237627992a99d5e30ae63e4d78cd24a" +uuid = "6c6a2e73-6563-6170-7368-637461726353" +version = "1.2.0" + +[[deps.Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[deps.Setfield]] +deps = ["ConstructionBase", "Future", "MacroTools", "StaticArraysCore"] +git-tree-sha1 = "e2cc6d8c88613c05e1defb55170bf5ff211fbeac" +uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46" +version = "1.1.1" + +[[deps.SharedArrays]] +deps = ["Distributed", "Mmap", "Random", "Serialization"] +uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" + +[[deps.Showoff]] +deps = ["Dates", "Grisu"] +git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de" +uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f" +version = "1.0.3" + +[[deps.SignedDistanceFields]] +deps = ["Random", "Statistics", "Test"] +git-tree-sha1 = "d263a08ec505853a5ff1c1ebde2070419e3f28e9" +uuid = "73760f76-fbc4-59ce-8f25-708e95d2df96" +version = "0.4.0" + +[[deps.SimpleTraits]] +deps = ["InteractiveUtils", "MacroTools"] +git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231" +uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d" +version = "0.9.4" + +[[deps.Sixel]] +deps = ["Dates", "FileIO", "ImageCore", "IndirectArrays", "OffsetArrays", "REPL", "libsixel_jll"] +git-tree-sha1 = "8fb59825be681d451c246a795117f317ecbcaa28" +uuid = "45858cf5-a6b0-47a3-bbea-62219f50df47" +version = "0.1.2" + +[[deps.SnoopPrecompile]] +deps = ["Preferences"] +git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c" +uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c" +version = "1.0.3" + +[[deps.Sockets]] +uuid = "6462fe0b-24de-5631-8697-dd941f90decc" + +[[deps.SortingAlgorithms]] +deps = ["DataStructures"] +git-tree-sha1 = "a4ada03f999bd01b3a25dcaa30b2d929fe537e00" +uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" +version = "1.1.0" + +[[deps.SparseArrays]] +deps = ["LinearAlgebra", "Random"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[deps.SpecialFunctions]] +deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] +git-tree-sha1 = "ef28127915f4229c971eb43f3fc075dd3fe91880" +uuid = "276daf66-3868-5448-9aa4-cd146d93841b" +version = "2.2.0" + +[[deps.StableHashTraits]] +deps = ["CRC32c", "Compat", "Dates", "SHA", "Tables", "TupleTools", "UUIDs"] +git-tree-sha1 = "0b8b801b8f03a329a4e86b44c5e8a7d7f4fe10a3" +uuid = "c5dd0088-6c3f-4803-b00e-f31a60c170fa" +version = "0.3.1" + +[[deps.StackViews]] +deps = ["OffsetArrays"] +git-tree-sha1 = "46e589465204cd0c08b4bd97385e4fa79a0c770c" +uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15" +version = "0.1.1" + +[[deps.StaticArrays]] +deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"] +git-tree-sha1 = "c262c8e978048c2b095be1672c9bee55b4619521" +uuid = "90137ffa-7385-5640-81b9-e52037218182" +version = "1.5.24" + +[[deps.StaticArraysCore]] +git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a" +uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" +version = "1.4.0" + +[[deps.Statistics]] +deps = ["LinearAlgebra", "SparseArrays"] +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[[deps.StatsAPI]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "45a7769a04a3cf80da1c1c7c60caf932e6f4c9f7" +uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" +version = "1.6.0" + +[[deps.StatsBase]] +deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] +git-tree-sha1 = "d1bf48bfcc554a3761a133fe3a9bb01488e06916" +uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +version = "0.33.21" + +[[deps.StatsFuns]] +deps = ["ChainRulesCore", "HypergeometricFunctions", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"] +git-tree-sha1 = "f625d686d5a88bcd2b15cd81f18f98186fdc0c9a" +uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c" +version = "1.3.0" + +[[deps.StructArrays]] +deps = ["Adapt", "DataAPI", "GPUArraysCore", "StaticArraysCore", "Tables"] +git-tree-sha1 = "521a0e828e98bb69042fec1809c1b5a680eb7389" +uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" +version = "0.6.15" + +[[deps.SuiteSparse]] +deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] +uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" + +[[deps.TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.0" + +[[deps.TableTraits]] +deps = ["IteratorInterfaceExtensions"] +git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39" +uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c" +version = "1.0.1" + +[[deps.Tables]] +deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"] +git-tree-sha1 = "1544b926975372da01227b382066ab70e574a3ec" +uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" +version = "1.10.1" + +[[deps.Tar]] +deps = ["ArgTools", "SHA"] +uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" +version = "1.10.1" + +[[deps.TensorCore]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6" +uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50" +version = "0.1.1" + +[[deps.Test]] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[[deps.TiffImages]] +deps = ["ColorTypes", "DataStructures", "DocStringExtensions", "FileIO", "FixedPointNumbers", "IndirectArrays", "Inflate", "Mmap", "OffsetArrays", "PkgVersion", "ProgressMeter", "UUIDs"] +git-tree-sha1 = "8621f5c499a8aa4aa970b1ae381aae0ef1576966" +uuid = "731e570b-9d59-4bfa-96dc-6df516fadf69" +version = "0.6.4" + +[[deps.TranscodingStreams]] +deps = ["Random", "Test"] +git-tree-sha1 = "0b829474fed270a4b0ab07117dce9b9a2fa7581a" +uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" +version = "0.9.12" + +[[deps.TriplotBase]] +git-tree-sha1 = "4d4ed7f294cda19382ff7de4c137d24d16adc89b" +uuid = "981d1d27-644d-49a2-9326-4793e63143c3" +version = "0.1.0" + +[[deps.TupleTools]] +git-tree-sha1 = "3c712976c47707ff893cf6ba4354aa14db1d8938" +uuid = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" +version = "1.3.0" + +[[deps.UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" + +[[deps.Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" + +[[deps.UnicodeFun]] +deps = ["REPL"] +git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf" +uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1" +version = "0.4.1" + +[[deps.WoodburyMatrices]] +deps = ["LinearAlgebra", "SparseArrays"] +git-tree-sha1 = "de67fa59e33ad156a590055375a30b23c40299d3" +uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6" +version = "0.5.5" + +[[deps.XML2_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"] +git-tree-sha1 = "93c41695bc1c08c46c5899f4fe06d6ead504bb73" +uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" +version = "2.10.3+0" + +[[deps.XSLT_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"] +git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a" +uuid = "aed1982a-8fda-507f-9586-7b0439959a61" +version = "1.1.34+0" + +[[deps.Xorg_libX11_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"] +git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527" +uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc" +version = "1.6.9+4" + +[[deps.Xorg_libXau_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e" +uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec" +version = "1.0.9+4" + +[[deps.Xorg_libXdmcp_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4" +uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05" +version = "1.1.3+4" + +[[deps.Xorg_libXext_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"] +git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3" +uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3" +version = "1.3.4+4" + +[[deps.Xorg_libXrender_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"] +git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96" +uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa" +version = "0.9.10+4" + +[[deps.Xorg_libpthread_stubs_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb" +uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74" +version = "0.1.0+3" + +[[deps.Xorg_libxcb_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"] +git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6" +uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b" +version = "1.13.0+3" + +[[deps.Xorg_xtrans_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845" +uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10" +version = "1.4.0+3" + +[[deps.Zlib_jll]] +deps = ["Libdl"] +uuid = "83775a58-1f1d-513f-b197-d71354ab007a" +version = "1.2.12+3" + +[[deps.isoband_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "51b5eeb3f98367157a7a12a1fb0aa5328946c03c" +uuid = "9a68df92-36a6-505f-a73e-abb412b6bfb4" +version = "0.2.3+0" + +[[deps.libaom_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "3a2ea60308f0996d26f1e5354e10c24e9ef905d4" +uuid = "a4ae2306-e953-59d6-aa16-d00cac43593b" +version = "3.4.0+0" + +[[deps.libass_jll]] +deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"] +git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47" +uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0" +version = "0.15.1+0" + +[[deps.libblastrampoline_jll]] +deps = ["Artifacts", "Libdl", "OpenBLAS_jll"] +uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.1.1+0" + +[[deps.libfdk_aac_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55" +uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280" +version = "2.0.2+0" + +[[deps.libpng_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"] +git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c" +uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f" +version = "1.6.38+0" + +[[deps.libsixel_jll]] +deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "libpng_jll"] +git-tree-sha1 = "d4f63314c8aa1e48cd22aa0c17ed76cd1ae48c3c" +uuid = "075b6546-f08a-558a-be8f-8157d0f608a5" +version = "1.10.3+0" + +[[deps.libvorbis_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"] +git-tree-sha1 = "b910cb81ef3fe6e78bf6acee440bda86fd6ae00c" +uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a" +version = "1.3.7+1" + +[[deps.nghttp2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" +version = "1.48.0+0" + +[[deps.p7zip_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" +version = "17.4.0+0" + +[[deps.x264_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2" +uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a" +version = "2021.5.5+0" + +[[deps.x265_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9" +uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76" +version = "3.5.0+0" +""" + +# ╔═╡ Cell order: +# ╠═74dccc96-e792-11ed-2b5a-b305a053025e +# ╠═13168eeb-afa0-4225-8e44-53b4d1b09874 +# ╠═3101cf7e-77fa-404f-a871-8ed4423885ff +# ╠═c68bb44c-f56d-486b-92b2-55adc7c0517f +# ╠═9c077253-7ece-464c-9c85-89ca1c7cfb9b +# ╠═b7efeb89-8de4-4f61-a2d2-d46364864e9e +# ╠═3b36b955-c489-4673-b8c3-56998354fae1 +# ╠═ba5d6ac1-727f-4be7-acd5-8c098d93d7ba +# ╠═046c6800-08f9-4e69-9570-0c802ce62423 +# ╠═5421714b-baf3-4136-9236-cb17ff7be80a +# ╠═85011187-8daa-4064-b93c-4624fa3dfcca +# ╠═59144084-5b0e-4d98-b3ca-52dafc657f0a +# ╠═ef502d5e-cfca-425b-96bc-385cc8dbd18b +# ╠═92782cc9-bc97-47df-a173-8be22bca7ddd +# ╠═3ac91add-acdd-4ec7-91f1-084a35367877 +# ╟─00000000-0000-0000-0000-000000000001 +# ╟─00000000-0000-0000-0000-000000000002 diff --git a/src/fft.jl b/src/fft.jl new file mode 100644 index 0000000..42f6770 --- /dev/null +++ b/src/fft.jl @@ -0,0 +1,24 @@ +module fft + +export fftn_mpi!, ifftn_mpi! + +using FFTW +using ..markers: AbstractParallel, ParallelMpi, SingleThread +using ..mesh: Wavenumbers + +fftn_mpi!(parallel::P, u, uhat) where P <: AbstractParallel = error("function is not not yet implemented. This is a bug") + +# single threaded forward FFT implementation for Base arrays +function fftn_mpi!(parallel::SingleThread, u::Array{Float64, 3}, uhat::Array{ComplexF64, 3}) + uhat[:, :, :] .= rfft(u, 1:3) +end + +ifftn_mpi!(parallel::P, wavenumbers::Wavenumbers, uhat, u) where P <: AbstractParallel = error("function is not not yet implemented. This is a bug") + +# single threaded inverse FFT implementation for Base arrays +function ifftn_mpi!(parallel::SingleThread, wavenumbers::Wavenumbers, uhat::Array{ComplexF64, 3}, u::Array{Float64, 3}) + u[:, :, :] .= irfft(uhat, wavenumbers.n, 1:3) +end + +# +end diff --git a/src/markers.jl b/src/markers.jl new file mode 100644 index 0000000..c529f9a --- /dev/null +++ b/src/markers.jl @@ -0,0 +1,11 @@ +module markers + +export AbstractParallel, ParallelMpi, SingleThread + +abstract type AbstractParallel end +struct ParallelMpi <: AbstractParallel end +struct SingleThread <: AbstractParallel end + + +# +end diff --git a/src/mesh.jl b/src/mesh.jl new file mode 100644 index 0000000..4fafeee --- /dev/null +++ b/src/mesh.jl @@ -0,0 +1,81 @@ +module mesh + +export Wavenumbers, wavenumbers +export Mesh, new_mesh + +using LazyGrids + +struct Wavenumbers + # number of x space values in each direction + n::Int + # number of modes in the x direction (N //2) +1 + # so for N = 64, this is 33 + kn::Int + kx::LazyGrids.GridAV{Int64, 1, 3} + ky::LazyGrids.GridAV{Int64, 2, 3} + kz::LazyGrids.GridAV{Int64, 3, 3} +end + +function wavenumbers(n::Int)::Wavenumbers + increasing_wavenumbers_positive = 0:(div(n,2)-1) + increasing_wavenumbers_negative = -(div(n,2)):-1 + + # fft frequencies replicatese numpy fft.fftfreq(N, 1/N) + ky_linear = vcat([increasing_wavenumbers_positive, increasing_wavenumbers_negative]...) + + # wave numbers in the z direction is identical to y direction + kz_linear = copy(ky_linear) + + # FFTW abbreviates the data in the x direction + # div(3,2) = 1, it is integer division + kn = div(n,2) + 1 + kx_linear = ky_linear[1:kn] + + kx, ky, kz = ndgrid(kx_linear, ky_linear, kz_linear) + + Wavenumbers(n, kn, kx, ky, kz) +end + +# make the Wavenumber type indexable +function Base.getindex(k_all::Wavenumbers, dim::Int)::LazyGrids.GridAV + if dim == 1 + return k_all.kx + elseif dim == 2 + return k_all.ky + elseif dim == 3 + return k_all.kz + else + error("dimension was not 1,2,3") + end +end + +struct Mesh + n::Int + x::LazyGrids.GridSL{Float64, 1, 3, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}} + y::LazyGrids.GridSL{Float64, 2, 3, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}} + z::LazyGrids.GridSL{Float64, 3, 3, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}} +end + +# nprange is range but for spectral things. It operates with the same spacing +# parameters as numpy's mgrid[] +# +# nprange(0, 1, 10) = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 +# range(0, 1, 10) = 0.0 0.111111 0.222222 0.333333 0.444444 0.555556 0.666667 0.777778 0.888889 1.0 +function nprange(start, _end, n) + step = (_end - start) / n + start:step:(_end - step) +end + +# create a mesh that is identical to numpy's mgrid[] spacing +function new_mesh(n) + xvals = nprange(0, 2pi, n) + yvals = copy(xvals) + zvals = copy(xvals) + + mesh_x, mesh_y, mesh_z = ndgrid(xvals, yvals, zvals); size(mesh_x) + + return Mesh(n, mesh_x, mesh_y, mesh_z) +end + +# +end diff --git a/src/spectralGPU.jl b/src/spectralGPU.jl index 0564f2b..af216ab 100644 --- a/src/spectralGPU.jl +++ b/src/spectralGPU.jl @@ -1,5 +1,12 @@ module spectralGPU -# Write your package code here. +include("./markers.jl") +include("./mesh.jl") +include("./fft.jl") +export fft +export mesh +export markers + +# end diff --git a/test/runtests.jl b/test/runtests.jl index 0e7a4e6..c7bce37 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,13 @@ -using spectralGPU +include("../src/spectralGPU.jl"); +using .spectralGPU: mesh, fft using Test -@testset "spectralGPU.jl" begin +@testset "mesh.jl" begin + # ensure this compiles + k = mesh.wavenumbers(32) + + @test k.kx[:, 1, 1] == k.kx[:, 2, 2] + @test k.ky[1, :, 1] == k.ky[2, :, 2] + @test k.kz[1, 1, :] == k.kz[2, 2, :] # Write your tests here. end From 92ce4fa6ccb8daeca228f9b0dca2bebb004d549f Mon Sep 17 00:00:00 2001 From: VanillaBrooks Date: Sat, 6 May 2023 16:47:19 -0700 Subject: [PATCH 03/14] tests for initial condition setup, curl, cross --- jupyter/fft_debug.py.ipynb | 82 +++++++++++++++++++++++++------------- pluto/fft_debug | 48 +++++++++++++++++++++- src/config.jl | 17 ++++++++ src/initial_condition.jl | 35 ++++++++++++++++ src/markers.jl | 11 +++++ src/mesh.jl | 2 +- src/solver.jl | 78 ++++++++++++++++++++++++++++++++++++ src/spectralGPU.jl | 8 ++++ src/state.jl | 44 ++++++++++++++++++++ test/runtests.jl | 58 ++++++++++++++++++++++++++- 10 files changed, 351 insertions(+), 32 deletions(-) create mode 100644 src/config.jl create mode 100644 src/initial_condition.jl create mode 100644 src/solver.jl create mode 100644 src/state.jl diff --git a/jupyter/fft_debug.py.ipynb b/jupyter/fft_debug.py.ipynb index ca58391..4395c05 100644 --- a/jupyter/fft_debug.py.ipynb +++ b/jupyter/fft_debug.py.ipynb @@ -2,8 +2,8 @@ "cells": [ { "cell_type": "code", - "execution_count": 54, - "id": "c5574ad8", + "execution_count": 1, + "id": "8f243277", "metadata": {}, "outputs": [], "source": [ @@ -15,8 +15,8 @@ }, { "cell_type": "code", - "execution_count": 55, - "id": "b92b2e5d", + "execution_count": 2, + "id": "1e929060", "metadata": {}, "outputs": [], "source": [ @@ -34,13 +34,20 @@ "#(3, N, Np, Np)\n", "K = np.array(np.meshgrid(kx, ky, kz, indexing='ij'), dtype=int)\n", "\n", - "X = np.mgrid[rank*Np:(rank+1)*Np, :N, :N].astype(float)*2*pi/N" + "X = np.mgrid[rank*Np:(rank+1)*Np, :N, :N].astype(float)*2*pi/N\n", + "\n", + "kmax_dealias = 2./3.*kn\n", + "dealias = np.array((abs(K[0]) < kmax_dealias)*(abs(K[1]) < kmax_dealias)*\n", + " (abs(K[2]) < kmax_dealias), dtype=bool)\n", + "\n", + "K2 = np.sum(K*K, 0, dtype=int)\n", + "K_over_K2 = K.astype(float) / np.where(K2 == 0, 1, K2).astype(float)" ] }, { "cell_type": "code", - "execution_count": 56, - "id": "e23269fc", + "execution_count": 3, + "id": "0d973ef2", "metadata": {}, "outputs": [ { @@ -62,29 +69,50 @@ }, { "cell_type": "code", - "execution_count": 68, - "id": "f3e91211", + "execution_count": 10, + "id": "e2caffd1", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ True, True, True, True, False, True, True, True])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "dealias[:, 0, 0]" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "5fc57a8e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "5.497787143782138" + "(8, 8, 5)" ] }, - "execution_count": 68, + "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "np.max(X[0])" + "dealias.shape" ] }, { "cell_type": "code", - "execution_count": 57, - "id": "f37638b8", + "execution_count": 5, + "id": "a2133510", "metadata": {}, "outputs": [], "source": [ @@ -97,8 +125,8 @@ }, { "cell_type": "code", - "execution_count": 58, - "id": "acc98409", + "execution_count": 6, + "id": "e4757549", "metadata": {}, "outputs": [], "source": [ @@ -109,8 +137,8 @@ }, { "cell_type": "code", - "execution_count": 59, - "id": "487cfa3b", + "execution_count": 7, + "id": "67dd846e", "metadata": {}, "outputs": [ { @@ -132,17 +160,17 @@ }, { "cell_type": "code", - "execution_count": 64, - "id": "61c1500f", + "execution_count": 8, + "id": "6ab5d3cf", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 64, + "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, @@ -163,17 +191,17 @@ }, { "cell_type": "code", - "execution_count": 67, - "id": "84b291bd", + "execution_count": 9, + "id": "bf76c6f6", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 67, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, @@ -195,7 +223,7 @@ { "cell_type": "code", "execution_count": null, - "id": "eae3f5fc", + "id": "3f5623fd", "metadata": {}, "outputs": [], "source": [] diff --git a/pluto/fft_debug b/pluto/fft_debug index f5fe187..30b93ae 100644 --- a/pluto/fft_debug +++ b/pluto/fft_debug @@ -40,11 +40,46 @@ names(sgpu) const par = sgpu.markers.SingleThread() # ╔═╡ 3b36b955-c489-4673-b8c3-56998354fae1 -const N = 128 +const N = 8 # ╔═╡ ba5d6ac1-727f-4be7-acd5-8c098d93d7ba const K = sgpu.mesh.wavenumbers(N) +# ╔═╡ 177df9a1-654f-4fcc-b4e2-c23e0c25a071 +const K2 = K[1].^2 + K[2].^2 + K[3].^2; + +# ╔═╡ 73edb142-7502-4556-830c-e4d253514929 +K.ky + +# ╔═╡ 09e454e4-2bc0-4515-ae4c-0051142394a9 +begin + increasing_wavenumbers_positive = 0:(div(N,2)-1) + increasing_wavenumbers_negative = -(div(N,2)):-1 + + # fft frequencies replicatese numpy fft.fftfreq(N, 1/N) + ky_linear = vcat([increasing_wavenumbers_positive, increasing_wavenumbers_negative]...) + kz_linear = copy(ky_linear) + kx_linear = ky_linear[1:5] +end + +# ╔═╡ ed52b05d-c1ec-4e5c-aaaa-5931d52399a5 +kx, ky, kz = ndgrid(kx_linear, ky_linear, kz_linear) + +# ╔═╡ c42cb77a-d195-47fe-872b-9a2fba864800 +ky ./ Array(kz) + +# ╔═╡ d907dfdb-1d0a-4b47-8722-30fd51fb3fed +kmax_dealias = 2. / 3. *K.kn + +# ╔═╡ 8ff60bd9-ab26-45ba-a097-045652b64697 +# dealias = np.array((abs(K[0]) < kmax_dealias)*(abs(K[1]) < kmax_dealias)* (abs(K[2]) < kmax_dealias), dtype=bool) + +# ╔═╡ 13fd1af0-b585-4d08-b13e-9c952bf04432 +dealias = ((abs.(kx) .< kmax_dealias) .* (abs.(ky) .< kmax_dealias) .* (abs.(kz) .< kmax_dealias)) + +# ╔═╡ 77626510-bff9-406d-a71b-289fcf3edcab +BitArray + # ╔═╡ 046c6800-08f9-4e69-9570-0c802ce62423 const mesh = sgpu.mesh.new_mesh(N) @@ -92,7 +127,7 @@ PLUTO_MANIFEST_TOML_CONTENTS = """ julia_version = "1.8.5" manifest_format = "2.0" -project_hash = "71a25de324d66428b531c134cc3204997729a01c" +project_hash = "182cf8de945373c65ed154e4107d9883bb4932e7" [[deps.AbstractFFTs]] deps = ["ChainRulesCore", "LinearAlgebra"] @@ -1366,6 +1401,15 @@ version = "3.5.0+0" # ╠═b7efeb89-8de4-4f61-a2d2-d46364864e9e # ╠═3b36b955-c489-4673-b8c3-56998354fae1 # ╠═ba5d6ac1-727f-4be7-acd5-8c098d93d7ba +# ╠═177df9a1-654f-4fcc-b4e2-c23e0c25a071 +# ╠═73edb142-7502-4556-830c-e4d253514929 +# ╠═09e454e4-2bc0-4515-ae4c-0051142394a9 +# ╠═ed52b05d-c1ec-4e5c-aaaa-5931d52399a5 +# ╠═c42cb77a-d195-47fe-872b-9a2fba864800 +# ╠═d907dfdb-1d0a-4b47-8722-30fd51fb3fed +# ╠═8ff60bd9-ab26-45ba-a097-045652b64697 +# ╠═13fd1af0-b585-4d08-b13e-9c952bf04432 +# ╠═77626510-bff9-406d-a71b-289fcf3edcab # ╠═046c6800-08f9-4e69-9570-0c802ce62423 # ╠═5421714b-baf3-4136-9236-cb17ff7be80a # ╠═85011187-8daa-4064-b93c-4624fa3dfcca diff --git a/src/config.jl b/src/config.jl new file mode 100644 index 0000000..995effd --- /dev/null +++ b/src/config.jl @@ -0,0 +1,17 @@ +module config + +export Config, create_config + +struct Config + ν::Float64 + re::Float64 + N::Int +end + +function create_config(N::Int, re::Float64)::Config + ν = 1 / re + + return Config(ν, re, N) +end + +end diff --git a/src/initial_condition.jl b/src/initial_condition.jl new file mode 100644 index 0000000..fae9f0a --- /dev/null +++ b/src/initial_condition.jl @@ -0,0 +1,35 @@ +module initial_condition + +using ..mesh: Mesh +using ..fft: fftn_mpi! +using ..markers: AbstractInitialCondition, AbstractParallel, TaylorGreen + +# catch-all for initial condition setup +function setup_initial_condition( + parallel::P, + ic::I, + mesh::Mesh, + U::AbstractArray{Float64, 4}, + U_hat::AbstractArray{ComplexF64, 4} +) where P <: AbstractParallel where I <: AbstractInitialCondition + error("unimplemented initial condition for type " * typeof(ic)) +end + +# Taylor-Green +function setup_initial_condition( + parallel::P, + ic::TaylorGreen, + mesh::Mesh, + U::AbstractArray{Float64, 4}, + U_hat::AbstractArray{ComplexF64, 4} +) where P <: AbstractParallel + U[:, :, :, 1] .= sin.(mesh.x) .* cos.(mesh.y) .* cos.(mesh.z) + U[:, :, :, 2] .= - cos.(mesh.x) .* sin.(mesh.y) .* cos.(mesh.z) + U[:, :, :, 3] .= 0. + + for i in 1:3 + fftn_mpi!(parallel, U[:, :, :, i], U_hat[:, :, :, i]) + end +end + +end diff --git a/src/markers.jl b/src/markers.jl index c529f9a..2143722 100644 --- a/src/markers.jl +++ b/src/markers.jl @@ -1,11 +1,22 @@ module markers +# +# Parallization (MPI, Single thread) +# export AbstractParallel, ParallelMpi, SingleThread abstract type AbstractParallel end struct ParallelMpi <: AbstractParallel end struct SingleThread <: AbstractParallel end +# +# Initial conditions +# +export AbstractInitialCondition, TaylorGreen, LoadInitialCondition + +abstract type AbstractInitialCondition end +struct TaylorGreen <: AbstractInitialCondition end +struct LoadInitialCondition <: AbstractInitialCondition end # end diff --git a/src/mesh.jl b/src/mesh.jl index 4fafeee..aa79dc8 100644 --- a/src/mesh.jl +++ b/src/mesh.jl @@ -67,7 +67,7 @@ function nprange(start, _end, n) end # create a mesh that is identical to numpy's mgrid[] spacing -function new_mesh(n) +function new_mesh(n::Int) xvals = nprange(0, 2pi, n) yvals = copy(xvals) zvals = copy(xvals) diff --git a/src/solver.jl b/src/solver.jl new file mode 100644 index 0000000..6958b3a --- /dev/null +++ b/src/solver.jl @@ -0,0 +1,78 @@ +module solver + +export curl!, cross!, compute_rhs! + +using ..markers: AbstractParallel +using ..fft: ifftn_mpi!, fftn_mpi! +using ..mesh: Wavenumbers +using ..state: State +using ..config: Config + +# calculate the curl of a fourier space array `input` and store the result in the x-space array `out` +function curl!( + parallel::P, + K::Wavenumbers, + input::Array{ComplexF64, 4} + ; + out::Array{Float64, 4} +) where P <: AbstractParallel + j = complex(0, 1) + ifftn_mpi!(parallel, K, j*(K[1].*input[:, :, :, 2] .- K[2].*input[:, :, :, 1]), out[:, :, :, 3]) + ifftn_mpi!(parallel, K, j*(K[3].*input[:, :, :, 1] .- K[1].*input[:, :, :, 3]), out[:, :, :, 2]) + ifftn_mpi!(parallel, K, j*(K[2].*input[:, :, :, 3] .- K[3].*input[:, :, :, 2]), out[:, :, :, 1]) + + nothing +end + +# take the cross product of X-Space a,b vectors +# stores the result in (fourier space) `out` and then return the +function cross!( + parallel::P, + a::Array{Float64, 4}, + b::Array{Float64, 4} + ; + out::Array{ComplexF64, 4} +) where P <: AbstractParallel + fftn_mpi!(parallel, a[:, :, :, 2].*b[:, :, :, 3] .- a[:, :, :, 3].*b[:, :, :, 2], out[:, :, :, 1]) + fftn_mpi!(parallel, a[:, :, :, 3].*b[:, :, :, 1] .- a[:, :, :, 1].*b[:, :, :, 3], out[:, :, :, 2]) + fftn_mpi!(parallel, a[:, :, :, 1].*b[:, :, :, 2] .- a[:, :, :, 2].*b[:, :, :, 1], out[:, :, :, 3]) + + nothing +end + +function compute_rhs!( + rk_step::Int, + parallel::P, + K::Wavenumbers, + config::Config, + U::Array{Float64, 4}, + U_hat::Array{ComplexF64, 4}, + state::State +) where P <: AbstractParallel + + if rk_step > 0 + for i in 1:3 + ifftn_mpi!(U_hat[:, :, :, i], U[:, :, :, i]) + end + end + + curl!(parallel, K, U_hat, out = state.curl) + #cross!(U, state.curl; out = state.dU) + #state.dU .*= dealias + + #state.P_hat[:, :, :] .= 0. + 0j; + #sum(P_hat, state.dU .* state.K_over_K², dims=0) + + ## TODO: make any adjustments to the forcing vector + # + ## add Pressure term to dU/dt + #state.dU -= state.P_hat * K + + ## dU is now the Eulerian term / Nonlinear term in the NSE + ## now calculate the diffusion term + + #state.dU -= config.ν * state.K² * U_hat +end + +# +end diff --git a/src/spectralGPU.jl b/src/spectralGPU.jl index af216ab..b9c9e0d 100644 --- a/src/spectralGPU.jl +++ b/src/spectralGPU.jl @@ -3,10 +3,18 @@ module spectralGPU include("./markers.jl") include("./mesh.jl") include("./fft.jl") +include("./initial_condition.jl") +include("./state.jl") +include("./config.jl") + +include("./solver.jl") export fft export mesh export markers +export state +export config +export initial_condition # end diff --git a/src/state.jl b/src/state.jl new file mode 100644 index 0000000..799e405 --- /dev/null +++ b/src/state.jl @@ -0,0 +1,44 @@ +module state + +using ..mesh: Wavenumbers + +export State, create_state + +struct State + dU::Array{ComplexF64, 4} + curl::Array{Float64, 4} + dealias::BitArray{3} + P_hat::Array{ComplexF64, 3} + K²::Array{Int, 3} + K_over_K²::Array{Float64, 4} +end + +function create_state(N::Int, K::Wavenumbers)::State + dU = ComplexF64.(zeros(K.kn, N, N, 3)) + curl = zeros(N, N, N, 3) + P_hat = ComplexF64.(zeros(K.kn, N, N)) + + K² = K[1].^2 + K[2].^2 + K[3].^3 + + kmax_dealias = 2/3 * K.kn + dealias = (abs.(K.kx) .< kmax_dealias) .* (abs.(K.ky) .< kmax_dealias) .* (abs.(K.kz) .< kmax_dealias) + + K_over_K² = Float64.(zeros(K.kn, N, N, 3)) + for i in 1:3 + K_over_K²[:, :, :, i] .= Float64.(K[i]) + end + + K_over_K² ./= K² + + State( + dU, + curl, + dealias, + P_hat, + K², + K_over_K² + ) +end + +# +end diff --git a/test/runtests.jl b/test/runtests.jl index c7bce37..b18d32a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,5 @@ include("../src/spectralGPU.jl"); -using .spectralGPU: mesh, fft +using .spectralGPU: mesh, fft, markers, initial_condition, state, config, solver using Test @testset "mesh.jl" begin @@ -9,5 +9,59 @@ using Test @test k.kx[:, 1, 1] == k.kx[:, 2, 2] @test k.ky[1, :, 1] == k.ky[2, :, 2] @test k.kz[1, 1, :] == k.kz[2, 2, :] - # Write your tests here. +end + +@testset "initial_condition.jl" begin + parallel = markers.SingleThread() + N = 64 + + K = mesh.wavenumbers(N) + U = zeros(N, N, N, 3) + U_hat = ComplexF64.(zeros(K.kn, N, N, 3)) + + msh = mesh.new_mesh(N) + + # taylor green + @test begin + ic = markers.TaylorGreen() + initial_condition.setup_initial_condition(parallel, ic, msh, U, U_hat) + true + end +end + +@testset "solver.jl" begin + parallel = markers.SingleThread() + N = 64 + re = 40. + + K = mesh.wavenumbers(N) + st = state.create_state(N, K) + msh = mesh.new_mesh(N) + cfg = config.create_config(N, re) + + U = zeros(N, N, N, 3) + U_hat = ComplexF64.(zeros(K.kn, N, N, 3)) + + # curl + @test begin + #solver.curl!(K, U_hat; out = st.curl[:, :, :, :]) + solver.curl!(parallel, K, U_hat; out=st.curl) + true + end + + # cross + @test begin + solver.cross!(parallel, U, st.curl; out = st.dU) + true + end + + #solver.compute_rhs!( + # 0, + # parallel, + # K, + # cfg, + # U, + # U_hat, + # st + #) end From 86a9ef305dc7352600b5619264392746ca376b6a Mon Sep 17 00:00:00 2001 From: VanillaBrooks Date: Sat, 6 May 2023 17:01:16 -0700 Subject: [PATCH 04/14] compute_rhs! solver tests passing --- src/solver.jl | 32 ++++++++++++++++++++------------ src/state.jl | 5 ++++- test/runtests.jl | 23 ++++++++++++++--------- 3 files changed, 38 insertions(+), 22 deletions(-) diff --git a/src/solver.jl b/src/solver.jl index 6958b3a..121f4f0 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -40,6 +40,12 @@ function cross!( nothing end +function wavenumber_product!(arr::Array{ComplexF64, 3}, K::Wavenumbers; out::Array{ComplexF64, 4}) + out[:, :, :, 1] .= arr .* K.kx + out[:, :, :, 2] .= arr .* K.ky + out[:, :, :, 3] .= arr .* K.kz +end + function compute_rhs!( rk_step::Int, parallel::P, @@ -56,22 +62,24 @@ function compute_rhs!( end end - curl!(parallel, K, U_hat, out = state.curl) - #cross!(U, state.curl; out = state.dU) - #state.dU .*= dealias + curl!(parallel, K, U_hat; out = state.curl) + cross!(parallel, U, state.curl; out = state.dU) + state.dU .*= state.dealias - #state.P_hat[:, :, :] .= 0. + 0j; - #sum(P_hat, state.dU .* state.K_over_K², dims=0) + state.P_hat[:, :, :] .= complex(0., 0); + # compute P_hat = sum(dU * K_over_K²): + sum!(state.P_hat, state.dU .* state.K_over_K²) - ## TODO: make any adjustments to the forcing vector - # - ## add Pressure term to dU/dt - #state.dU -= state.P_hat * K + # TODO: make any adjustments to the forcing vector + + # add Pressure term to dU/dt + wavenumber_product!(state.P_hat, K; out = state.wavenumber_product_tmp) + state.dU .-= state.wavenumber_product_tmp - ## dU is now the Eulerian term / Nonlinear term in the NSE - ## now calculate the diffusion term + # dU is now the Eulerian term / Nonlinear term in the NSE + # now calculate the diffusion term - #state.dU -= config.ν * state.K² * U_hat + state.dU .-= config.ν .* state.K² .* U_hat end # diff --git a/src/state.jl b/src/state.jl index 799e405..b7cfe65 100644 --- a/src/state.jl +++ b/src/state.jl @@ -11,12 +11,14 @@ struct State P_hat::Array{ComplexF64, 3} K²::Array{Int, 3} K_over_K²::Array{Float64, 4} + wavenumber_product_tmp::Array{ComplexF64, 4} end function create_state(N::Int, K::Wavenumbers)::State dU = ComplexF64.(zeros(K.kn, N, N, 3)) curl = zeros(N, N, N, 3) P_hat = ComplexF64.(zeros(K.kn, N, N)) + wavenumber_product_tmp = ComplexF64.(zeros(K.kn, N, N, 3)) K² = K[1].^2 + K[2].^2 + K[3].^3 @@ -36,7 +38,8 @@ function create_state(N::Int, K::Wavenumbers)::State dealias, P_hat, K², - K_over_K² + K_over_K², + wavenumber_product_tmp ) end diff --git a/test/runtests.jl b/test/runtests.jl index b18d32a..73f5404 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -55,13 +55,18 @@ end true end - #solver.compute_rhs!( - # 0, - # parallel, - # K, - # cfg, - # U, - # U_hat, - # st - #) + # main solver call + @test begin + solver.compute_rhs!( + 0, + parallel, + K, + cfg, + U, + U_hat, + st + ) + true + end + end From f65c4d18dcc6b998be0b9386d495c5b3e6c70195 Mon Sep 17 00:00:00 2001 From: VanillaBrooks Date: Sat, 6 May 2023 23:43:56 -0700 Subject: [PATCH 05/14] fix fft, ifft copies causing issues --- src/config.jl | 21 ++++++++-- src/fft.jl | 11 ++++-- src/initial_condition.jl | 10 ++--- src/integrate.jl | 59 +++++++++++++++++++++++++++++ src/solver.jl | 11 +++--- src/spectralGPU.jl | 1 + src/state.jl | 17 ++++++++- test/runtests.jl | 82 ++++++++++++++++++++++++++++++++++++++-- 8 files changed, 189 insertions(+), 23 deletions(-) create mode 100644 src/integrate.jl diff --git a/src/config.jl b/src/config.jl index 995effd..7837af3 100644 --- a/src/config.jl +++ b/src/config.jl @@ -1,17 +1,32 @@ module config -export Config, create_config +export Config, create_config, taylor_green_validation struct Config ν::Float64 re::Float64 N::Int + time::Float64 + dt::Float64 end -function create_config(N::Int, re::Float64)::Config +function create_config(N::Int, re::Float64, time::Float64)::Config ν = 1 / re - return Config(ν, re, N) + dt = 0.01 + + return Config(ν, re, N, time, dt) +end + +function taylor_green_validation()::Config + N = 64 + re = 1600 + ν = 1 / re + + t = 0.1 + dt = 0.01 + + return Config(ν, re, N, t, dt) end end diff --git a/src/fft.jl b/src/fft.jl index 42f6770..81b5cf3 100644 --- a/src/fft.jl +++ b/src/fft.jl @@ -6,17 +6,20 @@ using FFTW using ..markers: AbstractParallel, ParallelMpi, SingleThread using ..mesh: Wavenumbers -fftn_mpi!(parallel::P, u, uhat) where P <: AbstractParallel = error("function is not not yet implemented. This is a bug") +fftn_mpi!(parallel::P, u, uhat) where P <: AbstractParallel = error("function is not not yet implemented for $(typeof(uhat)), $(typeof(u)). This is a bug") # single threaded forward FFT implementation for Base arrays -function fftn_mpi!(parallel::SingleThread, u::Array{Float64, 3}, uhat::Array{ComplexF64, 3}) +function fftn_mpi!(parallel::SingleThread, u::XARR, uhat::SubArray{ComplexF64, 3}) where XARR <: AbstractArray{Float64, 3} uhat[:, :, :] .= rfft(u, 1:3) + + nothing end -ifftn_mpi!(parallel::P, wavenumbers::Wavenumbers, uhat, u) where P <: AbstractParallel = error("function is not not yet implemented. This is a bug") +ifftn_mpi!(parallel::P, wavenumbers::Wavenumbers, uhat, u) where P <: AbstractParallel = error("function is not not yet implemented for $(typeof(uhat)), $(typeof(u)). This is a bug") # single threaded inverse FFT implementation for Base arrays -function ifftn_mpi!(parallel::SingleThread, wavenumbers::Wavenumbers, uhat::Array{ComplexF64, 3}, u::Array{Float64, 3}) +function ifftn_mpi!(parallel::SingleThread, wavenumbers::Wavenumbers, uhat::FARR, u::SubArray{Float64, 3} +) where FARR <: AbstractArray{ComplexF64, 3} u[:, :, :] .= irfft(uhat, wavenumbers.n, 1:3) end diff --git a/src/initial_condition.jl b/src/initial_condition.jl index fae9f0a..0391586 100644 --- a/src/initial_condition.jl +++ b/src/initial_condition.jl @@ -9,8 +9,8 @@ function setup_initial_condition( parallel::P, ic::I, mesh::Mesh, - U::AbstractArray{Float64, 4}, - U_hat::AbstractArray{ComplexF64, 4} + U::Array{Float64, 4}, + U_hat::Array{ComplexF64, 4} ) where P <: AbstractParallel where I <: AbstractInitialCondition error("unimplemented initial condition for type " * typeof(ic)) end @@ -20,14 +20,14 @@ function setup_initial_condition( parallel::P, ic::TaylorGreen, mesh::Mesh, - U::AbstractArray{Float64, 4}, - U_hat::AbstractArray{ComplexF64, 4} + U::Array{Float64, 4}, + U_hat::Array{ComplexF64, 4} ) where P <: AbstractParallel U[:, :, :, 1] .= sin.(mesh.x) .* cos.(mesh.y) .* cos.(mesh.z) U[:, :, :, 2] .= - cos.(mesh.x) .* sin.(mesh.y) .* cos.(mesh.z) U[:, :, :, 3] .= 0. - for i in 1:3 + @views for i in 1:3 fftn_mpi!(parallel, U[:, :, :, i], U_hat[:, :, :, i]) end end diff --git a/src/integrate.jl b/src/integrate.jl new file mode 100644 index 0000000..bdeb653 --- /dev/null +++ b/src/integrate.jl @@ -0,0 +1,59 @@ +module Integrate + +using ..markers: AbstractParallel +using ..fft: ifftn_mpi!, fftn_mpi! +using ..mesh: Wavenumbers +using ..state: State +using ..config: Config +using ..solver: compute_rhs! + +function integrate( + parallel::P, + K::Wavenumbers, + config::Config, + state::State, + U::Array{Float64, 4}, + U_hat::Array{ComplexF64, 4} +) where P<: AbstractParallel + t = 0 + tstep = 0 + + # https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#Examples + # a_j is equal to the canonical b_i + # b_i is equal to the canonical a_j + a = [1/6, 1/3, 1/3, 1/6] + b = [0.5, 0.5, 1.] + + while t < config.time - 1e-8 + println("stepping"); + + t += config.dt + tstep += 1 + + state.U_hat₁[:, :, :, :] .= U_hat[:, :, :, :] + state.U_hat₀[:, :, :, :] .= U_hat[:, :, :, :] + + for rk_step in 1:4 + compute_rhs!(rk_step, parallel, K, config, U, U_hat, state) + + println("du abs change: ", sum(abs.(state.dU))) + + if rk_step < 4 + U_hat[:, :, :, :] .= state.U_hat₀ .+ b[rk_step] * config.dt * state.dU + end + + state.U_hat₁ .+= a[rk_step] * config.dt .* state.dU + end + + U_hat[:, :, :, :] .= state.U_hat₁ + + # update U from the rk time integrated U_hat + @views for i in 1:3 + ifftn_mpi!(parallel, K, U_hat[:, :, :, i], U[:, :, :, i]) + end + end +end + + +# +end diff --git a/src/solver.jl b/src/solver.jl index 121f4f0..9b25f78 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -9,7 +9,7 @@ using ..state: State using ..config: Config # calculate the curl of a fourier space array `input` and store the result in the x-space array `out` -function curl!( +@views function curl!( parallel::P, K::Wavenumbers, input::Array{ComplexF64, 4} @@ -17,6 +17,7 @@ function curl!( out::Array{Float64, 4} ) where P <: AbstractParallel j = complex(0, 1) + ifftn_mpi!(parallel, K, j*(K[1].*input[:, :, :, 2] .- K[2].*input[:, :, :, 1]), out[:, :, :, 3]) ifftn_mpi!(parallel, K, j*(K[3].*input[:, :, :, 1] .- K[1].*input[:, :, :, 3]), out[:, :, :, 2]) ifftn_mpi!(parallel, K, j*(K[2].*input[:, :, :, 3] .- K[3].*input[:, :, :, 2]), out[:, :, :, 1]) @@ -26,7 +27,7 @@ end # take the cross product of X-Space a,b vectors # stores the result in (fourier space) `out` and then return the -function cross!( +@views function cross!( parallel::P, a::Array{Float64, 4}, b::Array{Float64, 4} @@ -56,9 +57,9 @@ function compute_rhs!( state::State ) where P <: AbstractParallel - if rk_step > 0 - for i in 1:3 - ifftn_mpi!(U_hat[:, :, :, i], U[:, :, :, i]) + if rk_step != 1 + @views for i in 1:3 + ifftn_mpi!(parallel, K, U_hat[:, :, :, i], U[:, :, :, i]) end end diff --git a/src/spectralGPU.jl b/src/spectralGPU.jl index b9c9e0d..c2cd6c2 100644 --- a/src/spectralGPU.jl +++ b/src/spectralGPU.jl @@ -8,6 +8,7 @@ include("./state.jl") include("./config.jl") include("./solver.jl") +include("./integrate.jl") export fft export mesh diff --git a/src/state.jl b/src/state.jl index b7cfe65..24c3182 100644 --- a/src/state.jl +++ b/src/state.jl @@ -6,6 +6,8 @@ export State, create_state struct State dU::Array{ComplexF64, 4} + U_hat₁::Array{ComplexF64, 4} + U_hat₀::Array{ComplexF64, 4} curl::Array{Float64, 4} dealias::BitArray{3} P_hat::Array{ComplexF64, 3} @@ -16,11 +18,20 @@ end function create_state(N::Int, K::Wavenumbers)::State dU = ComplexF64.(zeros(K.kn, N, N, 3)) + U_hat₁ = ComplexF64.(zeros(K.kn, N, N, 3)) + U_hat₀ = ComplexF64.(zeros(K.kn, N, N, 3)) + curl = zeros(N, N, N, 3) P_hat = ComplexF64.(zeros(K.kn, N, N)) wavenumber_product_tmp = ComplexF64.(zeros(K.kn, N, N, 3)) - K² = K[1].^2 + K[2].^2 + K[3].^3 + K² = K[1].^2 + K[2].^2 + K[3].^2 + + # create a K² array that does not have any zeros + # this is for the K/K² calculation that follows. If this + # term has any zeros then it will produce NaNs in the division + K²_nonzero = copy(K²) + K²_nonzero[K² .== 0] .= 1 kmax_dealias = 2/3 * K.kn dealias = (abs.(K.kx) .< kmax_dealias) .* (abs.(K.ky) .< kmax_dealias) .* (abs.(K.kz) .< kmax_dealias) @@ -30,10 +41,12 @@ function create_state(N::Int, K::Wavenumbers)::State K_over_K²[:, :, :, i] .= Float64.(K[i]) end - K_over_K² ./= K² + K_over_K² ./= K²_nonzero State( dU, + U_hat₁, + U_hat₀, curl, dealias, P_hat, diff --git a/test/runtests.jl b/test/runtests.jl index 73f5404..e4216e0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,5 @@ include("../src/spectralGPU.jl"); -using .spectralGPU: mesh, fft, markers, initial_condition, state, config, solver +using .spectralGPU: mesh, fft, markers, initial_condition, state, config, solver, Integrate using Test @testset "mesh.jl" begin @@ -25,7 +25,14 @@ end @test begin ic = markers.TaylorGreen() initial_condition.setup_initial_condition(parallel, ic, msh, U, U_hat) - true + u_hat_sum = sum(abs.(U_hat)) + + if u_hat_sum == 0 + println("uhat sum was zero: ", u_hat_sum) + false + else + true + end end end @@ -37,7 +44,7 @@ end K = mesh.wavenumbers(N) st = state.create_state(N, K) msh = mesh.new_mesh(N) - cfg = config.create_config(N, re) + cfg = config.create_config(N, re, 1.0) U = zeros(N, N, N, 3) U_hat = ComplexF64.(zeros(K.kn, N, N, 3)) @@ -58,7 +65,7 @@ end # main solver call @test begin solver.compute_rhs!( - 0, + 2, parallel, K, cfg, @@ -70,3 +77,70 @@ end end end + +@testset "integrate.jl" begin + parallel = markers.SingleThread() + N = 64 + re = 40. + time = 0.05 + + K = mesh.wavenumbers(N) + st = state.create_state(N, K) + msh = mesh.new_mesh(N) + cfg = config.create_config(N, re, time) + + U = zeros(N, N, N, 3) + U_hat = ComplexF64.(zeros(K.kn, N, N, 3)) + + # main solver call + @test begin + Integrate.integrate( + parallel, + K, + cfg, + st, + U, + U_hat, + ) + true + end +end + +@testset "integrate.jl - checked" begin + parallel = markers.SingleThread() + N = 64 + + K = mesh.wavenumbers(N) + st = state.create_state(N, K) + msh = mesh.new_mesh(N) + cfg = config.taylor_green_validation() + ic = markers.TaylorGreen() + + U = zeros(N, N, N, 3) + U_hat = ComplexF64.(zeros(K.kn, N, N, 3)) + initial_condition.setup_initial_condition(parallel, ic, msh, U, U_hat) + + u_sum = sum(abs.(U)) + println("sum of all values in U is ", u_sum); + u_hat_sum = sum(abs.(U_hat)) + println("sum of all values in U_hat is ", u_hat_sum); + + # main solver call + @test begin + Integrate.integrate( + parallel, + K, + cfg, + st, + U, + U_hat, + ) + + u_sum = sum(abs.(U)) + println("sum of all values in U is ", u_sum); + + k = (1/2) * sum(U .* U) * (1 / N)^3 + println("k is ", k) + round(k - 0.124953117517; digits=7) == 0 + end +end From 5246ada3f92efadb3468aeb4853469cd5614cb4a Mon Sep 17 00:00:00 2001 From: VanillaBrooks Date: Tue, 9 May 2023 18:36:54 -0700 Subject: [PATCH 06/14] add julia dependencies to project.toml --- .gitignore | 1 + Project.toml | 5 +++++ base.apptainer | 20 ++++++++++++++++++++ justfile | 5 +++++ spectralGPU.apptainer | 11 +++++++++++ 5 files changed, 42 insertions(+) create mode 100644 base.apptainer create mode 100644 justfile create mode 100644 spectralGPU.apptainer diff --git a/.gitignore b/.gitignore index b067edd..dcc0adb 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,2 @@ /Manifest.toml +*.sif diff --git a/Project.toml b/Project.toml index ecc2b91..7ad35e8 100644 --- a/Project.toml +++ b/Project.toml @@ -3,6 +3,11 @@ uuid = "1f941c8d-9269-406b-be03-2a5a5e88bc1c" authors = ["brooks"] version = "1.0.0-DEV" +[deps] +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" +LazyGrids = "7031d0ef-c40d-4431-b2f8-61a8d2f650db" + [compat] julia = "1" diff --git a/base.apptainer b/base.apptainer new file mode 100644 index 0000000..6c7f1d1 --- /dev/null +++ b/base.apptainer @@ -0,0 +1,20 @@ +Bootstrap: docker +From: ubuntu:20.04 + +%files + ./ /spectralGPU + +%post + rm /bin/sh && ln -s /bin/bash /bin/sh + + apt update -y && apt upgrade -y + apt install -y curl libssl-dev + + curl -fsSL https://install.julialang.org | sh -s -- --yes + + . /root/.profile + + julia --version + +%apprun distribute + echo "hello!" diff --git a/justfile b/justfile new file mode 100644 index 0000000..3f6d177 --- /dev/null +++ b/justfile @@ -0,0 +1,5 @@ +base: + apptainer build base.sif ./base.apptainer + +build: + apptainer build spectralGPU.sif ./spectralGPU.apptainer diff --git a/spectralGPU.apptainer b/spectralGPU.apptainer new file mode 100644 index 0000000..4e5d0f0 --- /dev/null +++ b/spectralGPU.apptainer @@ -0,0 +1,11 @@ +Bootstrap: localimage +From: base.sif + +%files + ./ /spectralGPU + +%post + julia -e 'import Pkg; Pkg.add("FFTW"); Pkg.add("LazyGrids")' + +%apprun distribute + echo "hello!" From e787fa08815f4357889952056712493704b43301 Mon Sep 17 00:00:00 2001 From: VanillaBrooks Date: Tue, 9 May 2023 18:52:14 -0700 Subject: [PATCH 07/14] move tests to single threaded tests --- justfile | 4 ++-- spectralGPU.apptainer | 4 ++++ test/{runtests.jl => cpu_singlethread.jl} | 0 3 files changed, 6 insertions(+), 2 deletions(-) rename test/{runtests.jl => cpu_singlethread.jl} (100%) diff --git a/justfile b/justfile index 3f6d177..9506690 100644 --- a/justfile +++ b/justfile @@ -1,5 +1,5 @@ base: - apptainer build base.sif ./base.apptainer + sudo -E apptainer build base.sif ./base.apptainer build: - apptainer build spectralGPU.sif ./spectralGPU.apptainer + sudo -E apptainer build spectralGPU.sif ./spectralGPU.apptainer diff --git a/spectralGPU.apptainer b/spectralGPU.apptainer index 4e5d0f0..9616b7b 100644 --- a/spectralGPU.apptainer +++ b/spectralGPU.apptainer @@ -5,6 +5,10 @@ From: base.sif ./ /spectralGPU %post + # ensure julia is in the $PATH variable + . /root/.profile + + julia --version julia -e 'import Pkg; Pkg.add("FFTW"); Pkg.add("LazyGrids")' %apprun distribute diff --git a/test/runtests.jl b/test/cpu_singlethread.jl similarity index 100% rename from test/runtests.jl rename to test/cpu_singlethread.jl From 85b11fda568b0f00687aac9c8aa1c2d4f07d3e45 Mon Sep 17 00:00:00 2001 From: VanillaBrooks Date: Tue, 9 May 2023 21:54:14 -0700 Subject: [PATCH 08/14] CUDA cross products, curl, fft, ifft --- justfile | 4 ++ src/fft.jl | 25 +++++++-- src/initial_condition.jl | 21 ++++---- src/markers.jl | 5 +- src/solver.jl | 28 +++++----- test/cpu_singlethread.jl | 8 +-- test/cuda.jl | 112 +++++++++++++++++++++++++++++++++++++++ 7 files changed, 170 insertions(+), 33 deletions(-) create mode 100644 test/cuda.jl diff --git a/justfile b/justfile index 9506690..824fb26 100644 --- a/justfile +++ b/justfile @@ -3,3 +3,7 @@ base: build: sudo -E apptainer build spectralGPU.sif ./spectralGPU.apptainer + +test: + julia ./test/cpu_singlethread.jl + julia ./test/cuda.jl diff --git a/src/fft.jl b/src/fft.jl index 81b5cf3..adfe68c 100644 --- a/src/fft.jl +++ b/src/fft.jl @@ -3,14 +3,23 @@ module fft export fftn_mpi!, ifftn_mpi! using FFTW -using ..markers: AbstractParallel, ParallelMpi, SingleThread +using CUDA.CUFFT +using CUDA +using ..markers: AbstractParallel, ParallelMpi, SingleThreadCPU, SingleThreadGPU using ..mesh: Wavenumbers fftn_mpi!(parallel::P, u, uhat) where P <: AbstractParallel = error("function is not not yet implemented for $(typeof(uhat)), $(typeof(u)). This is a bug") # single threaded forward FFT implementation for Base arrays -function fftn_mpi!(parallel::SingleThread, u::XARR, uhat::SubArray{ComplexF64, 3}) where XARR <: AbstractArray{Float64, 3} - uhat[:, :, :] .= rfft(u, 1:3) +function fftn_mpi!(parallel::SingleThreadCPU, u::XARR, uhat::SubArray{ComplexF64, 3}) where XARR <: AbstractArray{Float64, 3} + uhat[:, :, :] .= FFTW.rfft(u, 1:3) + + nothing +end + +# single threaded forward FFT implementation for CUDA arrays +function fftn_mpi!(parallel::SingleThreadGPU, u::CuArray{Float64, 3}, uhat::CuArray{ComplexF64, 3}) + uhat[:, :, :] .= CUFFT.rfft(u, 1:3) nothing end @@ -18,9 +27,15 @@ end ifftn_mpi!(parallel::P, wavenumbers::Wavenumbers, uhat, u) where P <: AbstractParallel = error("function is not not yet implemented for $(typeof(uhat)), $(typeof(u)). This is a bug") # single threaded inverse FFT implementation for Base arrays -function ifftn_mpi!(parallel::SingleThread, wavenumbers::Wavenumbers, uhat::FARR, u::SubArray{Float64, 3} +function ifftn_mpi!(parallel::SingleThreadCPU, wavenumbers::Wavenumbers, uhat::FARR, u::SubArray{Float64, 3} ) where FARR <: AbstractArray{ComplexF64, 3} - u[:, :, :] .= irfft(uhat, wavenumbers.n, 1:3) + u[:, :, :] .= FFTW.irfft(uhat, wavenumbers.n, 1:3) +end + +# single threaded inverse FFT implementation for CUDA arrays +function ifftn_mpi!(parallel::SingleThreadGPU, wavenumbers::Wavenumbers, uhat::CuArray{ComplexF64, 3}, u::CuArray{Float64, 3} +) + u[:, :, :] .= CUFFT.irfft(uhat, wavenumbers.n, 1:3) end # diff --git a/src/initial_condition.jl b/src/initial_condition.jl index 0391586..98574af 100644 --- a/src/initial_condition.jl +++ b/src/initial_condition.jl @@ -4,14 +4,16 @@ using ..mesh: Mesh using ..fft: fftn_mpi! using ..markers: AbstractInitialCondition, AbstractParallel, TaylorGreen +using CUDA + # catch-all for initial condition setup function setup_initial_condition( parallel::P, ic::I, mesh::Mesh, - U::Array{Float64, 4}, - U_hat::Array{ComplexF64, 4} -) where P <: AbstractParallel where I <: AbstractInitialCondition + U::XARR, + U_hat::FARR, +) where P <: AbstractParallel where I <: AbstractInitialCondition where XARR <: AbstractArray{Float64, 4} where FARR <: AbstractArray{ComplexF64, 4} error("unimplemented initial condition for type " * typeof(ic)) end @@ -20,12 +22,13 @@ function setup_initial_condition( parallel::P, ic::TaylorGreen, mesh::Mesh, - U::Array{Float64, 4}, - U_hat::Array{ComplexF64, 4} -) where P <: AbstractParallel - U[:, :, :, 1] .= sin.(mesh.x) .* cos.(mesh.y) .* cos.(mesh.z) - U[:, :, :, 2] .= - cos.(mesh.x) .* sin.(mesh.y) .* cos.(mesh.z) - U[:, :, :, 3] .= 0. + U::XARR, + U_hat::FARR, +) where P <: AbstractParallel where XARR <: AbstractArray{Float64, 4} where FARR <: AbstractArray{ComplexF64, 4} + @CUDA.sync U[:, :, :, 1] .= sin.(mesh.x) .* cos.(mesh.y) .* cos.(mesh.z) + @CUDA.sync U[:, :, :, 2] .= cos.(mesh.x) .* sin.(mesh.y) .* cos.(mesh.z) + @CUDA.sync U[:, :, :, 2] *= -1 + @CUDA.sync U[:, :, :, 3] .= 0. @views for i in 1:3 fftn_mpi!(parallel, U[:, :, :, i], U_hat[:, :, :, i]) diff --git a/src/markers.jl b/src/markers.jl index 2143722..58c43e3 100644 --- a/src/markers.jl +++ b/src/markers.jl @@ -3,11 +3,12 @@ module markers # # Parallization (MPI, Single thread) # -export AbstractParallel, ParallelMpi, SingleThread +export AbstractParallel, ParallelMpi, SingleThreadCPU abstract type AbstractParallel end struct ParallelMpi <: AbstractParallel end -struct SingleThread <: AbstractParallel end +struct SingleThreadCPU <: AbstractParallel end +struct SingleThreadGPU <: AbstractParallel end # # Initial conditions diff --git a/src/solver.jl b/src/solver.jl index 9b25f78..7b1d08a 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -8,19 +8,21 @@ using ..mesh: Wavenumbers using ..state: State using ..config: Config +using CUDA + # calculate the curl of a fourier space array `input` and store the result in the x-space array `out` @views function curl!( parallel::P, K::Wavenumbers, - input::Array{ComplexF64, 4} + input::FARRAY ; - out::Array{Float64, 4} -) where P <: AbstractParallel + out::XARRAY +) where P <: AbstractParallel where FARRAY <: AbstractArray{ComplexF64, 4} where XARRAY <: AbstractArray{Float64, 4} j = complex(0, 1) - ifftn_mpi!(parallel, K, j*(K[1].*input[:, :, :, 2] .- K[2].*input[:, :, :, 1]), out[:, :, :, 3]) - ifftn_mpi!(parallel, K, j*(K[3].*input[:, :, :, 1] .- K[1].*input[:, :, :, 3]), out[:, :, :, 2]) - ifftn_mpi!(parallel, K, j*(K[2].*input[:, :, :, 3] .- K[3].*input[:, :, :, 2]), out[:, :, :, 1]) + @CUDA.sync ifftn_mpi!(parallel, K, j*(K[1].*input[:, :, :, 2] .- K[2].*input[:, :, :, 1]), out[:, :, :, 3]) + @CUDA.sync ifftn_mpi!(parallel, K, j*(K[3].*input[:, :, :, 1] .- K[1].*input[:, :, :, 3]), out[:, :, :, 2]) + @CUDA.sync ifftn_mpi!(parallel, K, j*(K[2].*input[:, :, :, 3] .- K[3].*input[:, :, :, 2]), out[:, :, :, 1]) nothing end @@ -29,14 +31,14 @@ end # stores the result in (fourier space) `out` and then return the @views function cross!( parallel::P, - a::Array{Float64, 4}, - b::Array{Float64, 4} + a::XARRAY, + b::XARRAY, ; - out::Array{ComplexF64, 4} -) where P <: AbstractParallel - fftn_mpi!(parallel, a[:, :, :, 2].*b[:, :, :, 3] .- a[:, :, :, 3].*b[:, :, :, 2], out[:, :, :, 1]) - fftn_mpi!(parallel, a[:, :, :, 3].*b[:, :, :, 1] .- a[:, :, :, 1].*b[:, :, :, 3], out[:, :, :, 2]) - fftn_mpi!(parallel, a[:, :, :, 1].*b[:, :, :, 2] .- a[:, :, :, 2].*b[:, :, :, 1], out[:, :, :, 3]) + out::FARRAY +) where P <: AbstractParallel where FARRAY <: AbstractArray{ComplexF64, 4} where XARRAY <: AbstractArray{Float64, 4} + @CUDA.sync fftn_mpi!(parallel, a[:, :, :, 2].*b[:, :, :, 3] .- a[:, :, :, 3].*b[:, :, :, 2], out[:, :, :, 1]) + @CUDA.sync fftn_mpi!(parallel, a[:, :, :, 3].*b[:, :, :, 1] .- a[:, :, :, 1].*b[:, :, :, 3], out[:, :, :, 2]) + @CUDA.sync fftn_mpi!(parallel, a[:, :, :, 1].*b[:, :, :, 2] .- a[:, :, :, 2].*b[:, :, :, 1], out[:, :, :, 3]) nothing end diff --git a/test/cpu_singlethread.jl b/test/cpu_singlethread.jl index e4216e0..73d8503 100644 --- a/test/cpu_singlethread.jl +++ b/test/cpu_singlethread.jl @@ -12,7 +12,7 @@ using Test end @testset "initial_condition.jl" begin - parallel = markers.SingleThread() + parallel = markers.SingleThreadCPU() N = 64 K = mesh.wavenumbers(N) @@ -37,7 +37,7 @@ end end @testset "solver.jl" begin - parallel = markers.SingleThread() + parallel = markers.SingleThreadCPU() N = 64 re = 40. @@ -79,7 +79,7 @@ end end @testset "integrate.jl" begin - parallel = markers.SingleThread() + parallel = markers.SingleThreadCPU() N = 64 re = 40. time = 0.05 @@ -107,7 +107,7 @@ end end @testset "integrate.jl - checked" begin - parallel = markers.SingleThread() + parallel = markers.SingleThreadCPU() N = 64 K = mesh.wavenumbers(N) diff --git a/test/cuda.jl b/test/cuda.jl new file mode 100644 index 0000000..0c39464 --- /dev/null +++ b/test/cuda.jl @@ -0,0 +1,112 @@ +include("../src/spectralGPU.jl"); +using .spectralGPU: mesh, fft, markers, initial_condition, state, config, solver, Integrate +using Test +using CUDA + + +@testset "cuda fft" begin + parallel = markers.SingleThreadGPU() + N = 64 + + K = mesh.wavenumbers(N) + + U = CuArray(zeros(N, N, N, 3)) + U_hat = CuArray(ComplexF64.(zeros(K.kn, N, N, 3))) + + @test begin + @views fft.fftn_mpi!(parallel, U[:, :, :, 1], U_hat[:, :, :, 1]) + true + end + + @test begin + @views fft.ifftn_mpi!(parallel, K, U_hat[:, :, :, 1], U[:, :, :, 1]) + true + end +end + + +@testset "cuda fft initial condition" begin + parallel = markers.SingleThreadGPU() + N = 64 + + K = mesh.wavenumbers(N) + + U = CuArray(zeros(N, N, N, 3)) + U_inverse = CuArray(zeros(N, N, N, 3)) + U_hat = CuArray(ComplexF64.(zeros(K.kn, N, N, 3))) + + ic = markers.TaylorGreen() + msh = mesh.new_mesh(N) + + initial_condition.setup_initial_condition(parallel, ic, msh, U, U_hat) + + u_hat_sum = sum(abs.(U_hat)) + + # ensure the velocity component has been populated + @test begin + U_norm = sum(abs.(U)) + U_norm > 1000 + end + + # ensure forward FFT works well for CUDA arrays + @test begin + if u_hat_sum == 0 + println("uhat sum was zero: ", u_hat_sum) + false + else + true + end + end + + # ensure the inverse FFT works equally well for CUDA arrays + @test begin + @views for i in 1:3 + fft.ifftn_mpi!(parallel, K, U_hat[:, :, :, i], U_inverse[:, :, :, i]) + end + + diff = U - U_inverse + l1_error = sum(abs.(diff)) + l1_error < 1e-5 + end +end + +#@testset "solver.jl" begin +# parallel = markers.SingleThreadGPU() +# N = 64 +# re = 40. +# +# K = mesh.wavenumbers(N) +# st = state.create_state(N, K) +# msh = mesh.new_mesh(N) +# cfg = config.create_config(N, re, 1.0) +# +# U = CuArray(zeros(N, N, N, 3)) +# U_hat = CuArray(ComplexF64.(zeros(K.kn, N, N, 3))) +# +# # curl +# @test begin +# solver.curl!(parallel, K, U_hat; out=st.curl) +# true +# end +# +# # cross +# @test begin +# solver.cross!(parallel, U, st.curl; out = st.dU) +# true +# end +# +# # main solver call +# @test begin +# solver.compute_rhs!( +# 2, +# parallel, +# K, +# cfg, +# U, +# U_hat, +# st +# ) +# true +# end +# +#end From bc5b9924f658b1293ffaf52dadfe3425c2a2c46b Mon Sep 17 00:00:00 2001 From: VanillaBrooks Date: Tue, 9 May 2023 23:51:35 -0700 Subject: [PATCH 09/14] all solver.jl tests passing w/ cuda, new WavenumbersGPU and StateGPU --- src/fft.jl | 4 +-- src/markers.jl | 10 ++++++ src/mesh.jl | 43 +++++++++++++++++++++-- src/solver.jl | 83 ++++++++++++++++++++++++++++++++++++++------ src/state.jl | 47 +++++++++++++++++++++---- test/cuda.jl | 93 +++++++++++++++++++++++++++----------------------- 6 files changed, 216 insertions(+), 64 deletions(-) diff --git a/src/fft.jl b/src/fft.jl index adfe68c..2a2fbdf 100644 --- a/src/fft.jl +++ b/src/fft.jl @@ -6,7 +6,7 @@ using FFTW using CUDA.CUFFT using CUDA using ..markers: AbstractParallel, ParallelMpi, SingleThreadCPU, SingleThreadGPU -using ..mesh: Wavenumbers +using ..mesh: Wavenumbers, WavenumbersGPU fftn_mpi!(parallel::P, u, uhat) where P <: AbstractParallel = error("function is not not yet implemented for $(typeof(uhat)), $(typeof(u)). This is a bug") @@ -33,7 +33,7 @@ function ifftn_mpi!(parallel::SingleThreadCPU, wavenumbers::Wavenumbers, uhat::F end # single threaded inverse FFT implementation for CUDA arrays -function ifftn_mpi!(parallel::SingleThreadGPU, wavenumbers::Wavenumbers, uhat::CuArray{ComplexF64, 3}, u::CuArray{Float64, 3} +function ifftn_mpi!(parallel::SingleThreadGPU, wavenumbers::WavenumbersGPU, uhat::CuArray{ComplexF64, 3}, u::CuArray{Float64, 3} ) u[:, :, :] .= CUFFT.irfft(uhat, wavenumbers.n, 1:3) end diff --git a/src/markers.jl b/src/markers.jl index 58c43e3..3c3098d 100644 --- a/src/markers.jl +++ b/src/markers.jl @@ -19,5 +19,15 @@ abstract type AbstractInitialCondition end struct TaylorGreen <: AbstractInitialCondition end struct LoadInitialCondition <: AbstractInitialCondition end +# State +export AbstractState + +abstract type AbstractState end + +# Wavenumbers +export AbstractWavenumbers + +abstract type AbstractWavenumbers end + # end diff --git a/src/mesh.jl b/src/mesh.jl index aa79dc8..ef03318 100644 --- a/src/mesh.jl +++ b/src/mesh.jl @@ -1,11 +1,14 @@ module mesh -export Wavenumbers, wavenumbers +using ..markers: AbstractWavenumbers + +export Wavenumbers, wavenumbers, WavenumbersGPU, wavenumbers_gpu export Mesh, new_mesh +using CUDA using LazyGrids -struct Wavenumbers +struct Wavenumbers <: AbstractWavenumbers # number of x space values in each direction n::Int # number of modes in the x direction (N //2) +1 @@ -16,6 +19,17 @@ struct Wavenumbers kz::LazyGrids.GridAV{Int64, 3, 3} end +struct WavenumbersGPU <: AbstractWavenumbers + # number of x space values in each direction + n::Int + # number of modes in the x direction (N //2) +1 + # so for N = 64, this is 33 + kn::Int + kx::CuArray{Int64, 3} + ky::CuArray{Int64, 3} + kz::CuArray{Int64, 3} +end + function wavenumbers(n::Int)::Wavenumbers increasing_wavenumbers_positive = 0:(div(n,2)-1) increasing_wavenumbers_negative = -(div(n,2)):-1 @@ -36,6 +50,18 @@ function wavenumbers(n::Int)::Wavenumbers Wavenumbers(n, kn, kx, ky, kz) end +function wavenumbers_gpu(n::Int)::WavenumbersGPU + K = wavenumbers(n) + + return WavenumbersGPU( + K.n, + K.kn, + CuArray(K.kx), + CuArray(K.ky), + CuArray(K.kz), + ) +end + # make the Wavenumber type indexable function Base.getindex(k_all::Wavenumbers, dim::Int)::LazyGrids.GridAV if dim == 1 @@ -49,6 +75,19 @@ function Base.getindex(k_all::Wavenumbers, dim::Int)::LazyGrids.GridAV end end +# make the Wavenumber type indexable +function Base.getindex(k_all::WavenumbersGPU, dim::Int)::CuArray{Int64, 3} + if dim == 1 + return k_all.kx + elseif dim == 2 + return k_all.ky + elseif dim == 3 + return k_all.kz + else + error("dimension was not 1,2,3") + end +end + struct Mesh n::Int x::LazyGrids.GridSL{Float64, 1, 3, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}} diff --git a/src/solver.jl b/src/solver.jl index 7b1d08a..a503c92 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -2,10 +2,10 @@ module solver export curl!, cross!, compute_rhs! -using ..markers: AbstractParallel +using ..markers: AbstractParallel, AbstractState, AbstractWavenumbers using ..fft: ifftn_mpi!, fftn_mpi! -using ..mesh: Wavenumbers -using ..state: State +using ..mesh: Wavenumbers, WavenumbersGPU +using ..state: State, StateGPU using ..config: Config using CUDA @@ -14,28 +14,59 @@ using CUDA @views function curl!( parallel::P, K::Wavenumbers, - input::FARRAY + input::Array{ComplexF64, 4} ; - out::XARRAY -) where P <: AbstractParallel where FARRAY <: AbstractArray{ComplexF64, 4} where XARRAY <: AbstractArray{Float64, 4} + out::Array{Float64, 4} +) where P <: AbstractParallel j = complex(0, 1) - @CUDA.sync ifftn_mpi!(parallel, K, j*(K[1].*input[:, :, :, 2] .- K[2].*input[:, :, :, 1]), out[:, :, :, 3]) + ifftn_mpi!(parallel, K, j*(K[1].*input[:, :, :, 2] .- K[2].*input[:, :, :, 1]), out[:, :, :, 3]) + ifftn_mpi!(parallel, K, j*(K[3].*input[:, :, :, 1] .- K[1].*input[:, :, :, 3]), out[:, :, :, 2]) + ifftn_mpi!(parallel, K, j*(K[2].*input[:, :, :, 3] .- K[3].*input[:, :, :, 2]), out[:, :, :, 1]) + + nothing +end + +@views function curl!( + parallel::P, + K::WavenumbersGPU, + input::CuArray{ComplexF64, 4} + ; + out::CuArray{Float64, 4} +) where P <: AbstractParallel + j = complex(0, 1) + + @CUDA.sync ifftn_mpi!(parallel, K,j * (CuArray(K[1]).*input[:, :, :, 2] .- CuArray(K[2]).*input[:, :, :, 1]), out[:, :, :, 3]) @CUDA.sync ifftn_mpi!(parallel, K, j*(K[3].*input[:, :, :, 1] .- K[1].*input[:, :, :, 3]), out[:, :, :, 2]) @CUDA.sync ifftn_mpi!(parallel, K, j*(K[2].*input[:, :, :, 3] .- K[3].*input[:, :, :, 2]), out[:, :, :, 1]) nothing end + # take the cross product of X-Space a,b vectors # stores the result in (fourier space) `out` and then return the @views function cross!( parallel::P, - a::XARRAY, - b::XARRAY, + a::Array{Float64, 4}, + b::Array{Float64, 4}, + ; + out::Array{ComplexF64, 4} +) where P <: AbstractParallel + fftn_mpi!(parallel, a[:, :, :, 2].*b[:, :, :, 3] .- a[:, :, :, 3].*b[:, :, :, 2], out[:, :, :, 1]) + fftn_mpi!(parallel, a[:, :, :, 3].*b[:, :, :, 1] .- a[:, :, :, 1].*b[:, :, :, 3], out[:, :, :, 2]) + fftn_mpi!(parallel, a[:, :, :, 1].*b[:, :, :, 2] .- a[:, :, :, 2].*b[:, :, :, 1], out[:, :, :, 3]) + + nothing +end + +@views function cross!( + parallel::P, + a::CuArray{Float64, 4}, + b::CuArray{Float64, 4}, ; - out::FARRAY -) where P <: AbstractParallel where FARRAY <: AbstractArray{ComplexF64, 4} where XARRAY <: AbstractArray{Float64, 4} + out::CuArray{ComplexF64, 4} +) where P <: AbstractParallel @CUDA.sync fftn_mpi!(parallel, a[:, :, :, 2].*b[:, :, :, 3] .- a[:, :, :, 3].*b[:, :, :, 2], out[:, :, :, 1]) @CUDA.sync fftn_mpi!(parallel, a[:, :, :, 3].*b[:, :, :, 1] .- a[:, :, :, 1].*b[:, :, :, 3], out[:, :, :, 2]) @CUDA.sync fftn_mpi!(parallel, a[:, :, :, 1].*b[:, :, :, 2] .- a[:, :, :, 2].*b[:, :, :, 1], out[:, :, :, 3]) @@ -49,6 +80,24 @@ function wavenumber_product!(arr::Array{ComplexF64, 3}, K::Wavenumbers; out::Arr out[:, :, :, 3] .= arr .* K.kz end +function wavenumber_product!(arr::CuArray{ComplexF64, 3}, K::WavenumbersGPU; out::CuArray{ComplexF64, 4}) + out[:, :, :, 1] .= arr .* K.kx + out[:, :, :, 2] .= arr .* K.ky + out[:, :, :, 3] .= arr .* K.kz +end + +function compute_rhs!( + rk_step::Int, + parallel::P, + K::WavenumbersGPU, + config::Config, + U::CuArray{Float64, 4}, + U_hat::CuArray{ComplexF64, 4}, + state::StateGPU +) where P <: AbstractParallel + __compute_rhs!(rk_step, parallel, K, config, U, U_hat, state) +end + function compute_rhs!( rk_step::Int, parallel::P, @@ -58,6 +107,18 @@ function compute_rhs!( U_hat::Array{ComplexF64, 4}, state::State ) where P <: AbstractParallel + __compute_rhs!(rk_step, parallel, K, config, U, U_hat, state) +end + +function __compute_rhs!( + rk_step::Int, + parallel::P, + K::WAVE, + config::Config, + U::ARRAY, + U_hat::FARRAY, + state::STATE +) where P <: AbstractParallel where FARRAY <: AbstractArray{ComplexF64, 4} where ARRAY <: AbstractArray{Float64, 4} where STATE <: AbstractState where WAVE <: AbstractWavenumbers if rk_step != 1 @views for i in 1:3 diff --git a/src/state.jl b/src/state.jl index 24c3182..a70eb36 100644 --- a/src/state.jl +++ b/src/state.jl @@ -1,10 +1,13 @@ module state using ..mesh: Wavenumbers +using ..markers: AbstractState, AbstractWavenumbers -export State, create_state +using CUDA -struct State +export State, create_state, create_state_gpu + +struct State <: AbstractState dU::Array{ComplexF64, 4} U_hat₁::Array{ComplexF64, 4} U_hat₀::Array{ComplexF64, 4} @@ -16,7 +19,19 @@ struct State wavenumber_product_tmp::Array{ComplexF64, 4} end -function create_state(N::Int, K::Wavenumbers)::State +struct StateGPU <: AbstractState + dU::CuArray{ComplexF64, 4} + U_hat₁::CuArray{ComplexF64, 4} + U_hat₀::CuArray{ComplexF64, 4} + curl::CuArray{Float64, 4} + dealias::CuArray{Int8, 3} + P_hat::CuArray{ComplexF64, 3} + K²::CuArray{Int, 3} + K_over_K²::CuArray{Float64, 4} + wavenumber_product_tmp::CuArray{ComplexF64, 4} +end + +function create_state(N::Int, K::WAVE)::State where WAVE <: AbstractWavenumbers dU = ComplexF64.(zeros(K.kn, N, N, 3)) U_hat₁ = ComplexF64.(zeros(K.kn, N, N, 3)) U_hat₀ = ComplexF64.(zeros(K.kn, N, N, 3)) @@ -25,20 +40,20 @@ function create_state(N::Int, K::Wavenumbers)::State P_hat = ComplexF64.(zeros(K.kn, N, N)) wavenumber_product_tmp = ComplexF64.(zeros(K.kn, N, N, 3)) - K² = K[1].^2 + K[2].^2 + K[3].^2 + K² = Array(K[1].^2 + K[2].^2 + K[3].^2) # create a K² array that does not have any zeros # this is for the K/K² calculation that follows. If this # term has any zeros then it will produce NaNs in the division - K²_nonzero = copy(K²) + K²_nonzero = Array(copy(K²)) K²_nonzero[K² .== 0] .= 1 kmax_dealias = 2/3 * K.kn - dealias = (abs.(K.kx) .< kmax_dealias) .* (abs.(K.ky) .< kmax_dealias) .* (abs.(K.kz) .< kmax_dealias) + dealias = (abs.(Array(K.kx)) .< kmax_dealias) .* (abs.(Array(K.ky)) .< kmax_dealias) .* (abs.(Array(K.kz)) .< kmax_dealias) K_over_K² = Float64.(zeros(K.kn, N, N, 3)) for i in 1:3 - K_over_K²[:, :, :, i] .= Float64.(K[i]) + K_over_K²[:, :, :, i] .= Float64.(Array(K[i])) end K_over_K² ./= K²_nonzero @@ -56,5 +71,23 @@ function create_state(N::Int, K::Wavenumbers)::State ) end +function create_state_gpu(N::Int, K::WAVE)::StateGPU where WAVE <: AbstractWavenumbers + cpu_state = create_state(N, K) + dealias::Array{Int8, 3} = Array(cpu_state.dealias) + + return StateGPU( + # + CuArray(cpu_state.dU), + CuArray(cpu_state.U_hat₁), + CuArray(cpu_state.U_hat₀), + CuArray(cpu_state.curl), + CuArray(dealias), + CuArray(cpu_state.P_hat), + CuArray(cpu_state.K²), + CuArray(cpu_state.K_over_K²), + CuArray(cpu_state.wavenumber_product_tmp), + ) +end + # end diff --git a/test/cuda.jl b/test/cuda.jl index 0c39464..a3405a6 100644 --- a/test/cuda.jl +++ b/test/cuda.jl @@ -3,12 +3,21 @@ using .spectralGPU: mesh, fft, markers, initial_condition, state, config, solver using Test using CUDA +@testset "state.jl" begin + N = 64 + K = mesh.wavenumbers_gpu(N) + @test begin + st::state.StateGPU = state.create_state_gpu(N, K) + true + end +end + @testset "cuda fft" begin parallel = markers.SingleThreadGPU() N = 64 - K = mesh.wavenumbers(N) + K = mesh.wavenumbers_gpu(N) U = CuArray(zeros(N, N, N, 3)) U_hat = CuArray(ComplexF64.(zeros(K.kn, N, N, 3))) @@ -29,7 +38,7 @@ end parallel = markers.SingleThreadGPU() N = 64 - K = mesh.wavenumbers(N) + K = mesh.wavenumbers_gpu(N) U = CuArray(zeros(N, N, N, 3)) U_inverse = CuArray(zeros(N, N, N, 3)) @@ -70,43 +79,43 @@ end end end -#@testset "solver.jl" begin -# parallel = markers.SingleThreadGPU() -# N = 64 -# re = 40. -# -# K = mesh.wavenumbers(N) -# st = state.create_state(N, K) -# msh = mesh.new_mesh(N) -# cfg = config.create_config(N, re, 1.0) -# -# U = CuArray(zeros(N, N, N, 3)) -# U_hat = CuArray(ComplexF64.(zeros(K.kn, N, N, 3))) -# -# # curl -# @test begin -# solver.curl!(parallel, K, U_hat; out=st.curl) -# true -# end -# -# # cross -# @test begin -# solver.cross!(parallel, U, st.curl; out = st.dU) -# true -# end -# -# # main solver call -# @test begin -# solver.compute_rhs!( -# 2, -# parallel, -# K, -# cfg, -# U, -# U_hat, -# st -# ) -# true -# end -# -#end +@testset "solver.jl" begin + parallel = markers.SingleThreadGPU() + N = 64 + re = 40. + + K = mesh.wavenumbers_gpu(N) + st = state.create_state_gpu(N, K) + msh = mesh.new_mesh(N) + cfg = config.create_config(N, re, 1.0) + + U = CuArray(zeros(N, N, N, 3)) + U_hat = CuArray(ComplexF64.(zeros(K.kn, N, N, 3))) + + # curl + @test begin + solver.curl!(parallel, K, U_hat; out=st.curl) + true + end + + # cross + @test begin + solver.cross!(parallel, U, st.curl; out = st.dU) + true + end + + # main solver call + @test begin + solver.compute_rhs!( + 2, + parallel, + K, + cfg, + U, + U_hat, + st + ) + true + end + +end From ca8ee5cbdfdd7817f49f7ec809078baff4316fb7 Mon Sep 17 00:00:00 2001 From: VanillaBrooks Date: Tue, 9 May 2023 23:59:02 -0700 Subject: [PATCH 10/14] RK integration with cuda --- justfile | 2 +- src/integrate.jl | 30 +++++++++++++++++++++++++++--- test/cuda.jl | 39 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 67 insertions(+), 4 deletions(-) diff --git a/justfile b/justfile index 824fb26..ecf682f 100644 --- a/justfile +++ b/justfile @@ -5,5 +5,5 @@ build: sudo -E apptainer build spectralGPU.sif ./spectralGPU.apptainer test: - julia ./test/cpu_singlethread.jl + #julia ./test/cpu_singlethread.jl julia ./test/cuda.jl diff --git a/src/integrate.jl b/src/integrate.jl index bdeb653..f2f9dfa 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -1,12 +1,14 @@ module Integrate -using ..markers: AbstractParallel +using ..markers: AbstractParallel, AbstractState, AbstractWavenumbers using ..fft: ifftn_mpi!, fftn_mpi! -using ..mesh: Wavenumbers -using ..state: State +using ..mesh: Wavenumbers, WavenumbersGPU +using ..state: State, StateGPU using ..config: Config using ..solver: compute_rhs! +using CUDA + function integrate( parallel::P, K::Wavenumbers, @@ -15,6 +17,28 @@ function integrate( U::Array{Float64, 4}, U_hat::Array{ComplexF64, 4} ) where P<: AbstractParallel + __integrate(parallel, K, config, state, U, U_hat) +end + +function integrate( + parallel::P, + K::WavenumbersGPU, + config::Config, + state::StateGPU, + U::CuArray{Float64, 4}, + U_hat::CuArray{ComplexF64, 4} +) where P<: AbstractParallel + __integrate(parallel, K, config, state, U, U_hat) +end + +function __integrate( + parallel::P, + K::WAVE, + config::Config, + state::STATE, + U::XARRAY, + U_hat::FARRAY, +) where P<: AbstractParallel where FARRAY <: AbstractArray{ComplexF64, 4} where XARRAY <: AbstractArray{Float64, 4} where STATE <: AbstractState where WAVE <: AbstractWavenumbers t = 0 tstep = 0 diff --git a/test/cuda.jl b/test/cuda.jl index a3405a6..d9d61fa 100644 --- a/test/cuda.jl +++ b/test/cuda.jl @@ -119,3 +119,42 @@ end end end + +@testset "integrate.jl - checked" begin + parallel = markers.SingleThreadGPU() + N = 64 + + K = mesh.wavenumbers_gpu(N) + st = state.create_state_gpu(N, K) + msh = mesh.new_mesh(N) + cfg = config.taylor_green_validation() + ic = markers.TaylorGreen() + + U = CuArray(zeros(N, N, N, 3)) + U_hat = CuArray(ComplexF64.(zeros(K.kn, N, N, 3))) + initial_condition.setup_initial_condition(parallel, ic, msh, U, U_hat) + + u_sum = sum(abs.(U)) + println("sum of all values in U is ", u_sum); + u_hat_sum = sum(abs.(U_hat)) + println("sum of all values in U_hat is ", u_hat_sum); + + # main solver call + @test begin + Integrate.integrate( + parallel, + K, + cfg, + st, + U, + U_hat, + ) + + u_sum = sum(abs.(U)) + println("sum of all values in U is ", u_sum); + + k = (1/2) * sum(U .* U) * (1 / N)^3 + println("k is ", k) + round(k - 0.124953117517; digits=7) == 0 + end +end From 91eb4a32d5aa81ed4719544ce031858fe4dba9f2 Mon Sep 17 00:00:00 2001 From: VanillaBrooks Date: Wed, 10 May 2023 01:02:03 -0700 Subject: [PATCH 11/14] benchmarks --- benches/cpu_gpu.jl | 65 ++++++++++++++++++++++++++++++++++++++++++++++ justfile | 3 +++ src/config.jl | 35 +++++++++++++++++++------ src/integrate.jl | 28 ++++++++++---------- src/markers.jl | 5 ++++ src/solver.jl | 14 +++++----- 6 files changed, 122 insertions(+), 28 deletions(-) create mode 100644 benches/cpu_gpu.jl diff --git a/benches/cpu_gpu.jl b/benches/cpu_gpu.jl new file mode 100644 index 0000000..a537d0d --- /dev/null +++ b/benches/cpu_gpu.jl @@ -0,0 +1,65 @@ +include("../src/spectralGPU.jl"); +using .spectralGPU: mesh, fft, markers, initial_condition, state, config, solver, Integrate + +using CUDA +using BenchmarkTools + +######### +######### GPU benchmark +######### + +N = 256 +begin + local parallel = markers.SingleThreadGPU() + + local K = mesh.wavenumbers_gpu(N) + local st = state.create_state_gpu(N, K) + local msh = mesh.new_mesh(N) + local cfg = config.taylor_green_validation() + local ic = markers.TaylorGreen() + + local U = CuArray(zeros(N, N, N, 3)) + local U_hat = CuArray(ComplexF64.(zeros(K.kn, N, N, 3))) + initial_condition.setup_initial_condition(parallel, ic, msh, U, U_hat) + + r = @benchmark Integrate.integrate( + $parallel, + $K, + $cfg, + $st, + $U, + $U_hat, + ) + + show(r) + println("\n"); +end + +######### +######### CPU benchmark +######### + +begin + local parallel = markers.SingleThreadCPU() + + local K = mesh.wavenumbers(N) + local st = state.create_state(N, K) + local msh = mesh.new_mesh(N) + local cfg = config.taylor_green_validation() + local ic = markers.TaylorGreen() + + local U = zeros(N, N, N, 3) + local U_hat = ComplexF64.(zeros(K.kn, N, N, 3)) + initial_condition.setup_initial_condition(parallel, ic, msh, U, U_hat) + + r = @benchmark Integrate.integrate( + $parallel, + $K, + $cfg, + $st, + $U, + $U_hat, + ) + + show(r) +end diff --git a/justfile b/justfile index ecf682f..3dc0758 100644 --- a/justfile +++ b/justfile @@ -7,3 +7,6 @@ build: test: #julia ./test/cpu_singlethread.jl julia ./test/cuda.jl + +bench: + julia ./benches/cpu_gpu.jl diff --git a/src/config.jl b/src/config.jl index 7837af3..29df6e7 100644 --- a/src/config.jl +++ b/src/config.jl @@ -1,32 +1,51 @@ module config -export Config, create_config, taylor_green_validation +using ..markers: AbstractConfig -struct Config +export Config, create_config, taylor_green_validation, calculate_dt + +struct Config <: AbstractConfig + ν::Float64 + re::Float64 + N::Int + time::Float64 +end + +struct ConfigValidation <: AbstractConfig ν::Float64 re::Float64 N::Int time::Float64 - dt::Float64 end function create_config(N::Int, re::Float64, time::Float64)::Config ν = 1 / re - dt = 0.01 + CFL_NUM = 0.75 - return Config(ν, re, N, time, dt) + return Config(ν, re, N, time) end -function taylor_green_validation()::Config +function taylor_green_validation()::ConfigValidation N = 64 re = 1600 ν = 1 / re t = 0.1 - dt = 0.01 - return Config(ν, re, N, t, dt) + return ConfigValidation(ν, re, N, t) +end + +function calculate_dt(velocity::XARRAY, cfg::Config)::Float64 where XARRAY <: AbstractArray{Float64, 4} + max_velo = maximum(velocity) + CFL = 0.75 + + return CFL / (cfg.N * max_velo) +end + +function calculate_dt(velocity::XARRAY, cfg::ConfigValidation)::Float64 where XARRAY <: AbstractArray{Float64, 4} + return 0.01 end +# end diff --git a/src/integrate.jl b/src/integrate.jl index f2f9dfa..3dc174f 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -1,10 +1,10 @@ module Integrate -using ..markers: AbstractParallel, AbstractState, AbstractWavenumbers +using ..markers: AbstractParallel, AbstractState, AbstractWavenumbers, AbstractConfig using ..fft: ifftn_mpi!, fftn_mpi! using ..mesh: Wavenumbers, WavenumbersGPU using ..state: State, StateGPU -using ..config: Config +using ..config: Config, calculate_dt using ..solver: compute_rhs! using CUDA @@ -12,33 +12,33 @@ using CUDA function integrate( parallel::P, K::Wavenumbers, - config::Config, + config::CFG, state::State, U::Array{Float64, 4}, U_hat::Array{ComplexF64, 4} -) where P<: AbstractParallel +) where P<: AbstractParallel where CFG <: AbstractConfig __integrate(parallel, K, config, state, U, U_hat) end function integrate( parallel::P, K::WavenumbersGPU, - config::Config, + config::CFG, state::StateGPU, U::CuArray{Float64, 4}, U_hat::CuArray{ComplexF64, 4} -) where P<: AbstractParallel +) where P<: AbstractParallel where CFG <: AbstractConfig __integrate(parallel, K, config, state, U, U_hat) end function __integrate( parallel::P, K::WAVE, - config::Config, + config::CFG, state::STATE, U::XARRAY, U_hat::FARRAY, -) where P<: AbstractParallel where FARRAY <: AbstractArray{ComplexF64, 4} where XARRAY <: AbstractArray{Float64, 4} where STATE <: AbstractState where WAVE <: AbstractWavenumbers +) where P<: AbstractParallel where FARRAY <: AbstractArray{ComplexF64, 4} where XARRAY <: AbstractArray{Float64, 4} where STATE <: AbstractState where WAVE <: AbstractWavenumbers where CFG <: AbstractConfig t = 0 tstep = 0 @@ -48,10 +48,12 @@ function __integrate( a = [1/6, 1/3, 1/3, 1/6] b = [0.5, 0.5, 1.] + dt = calculate_dt(U, config) + while t < config.time - 1e-8 - println("stepping"); + #println("stepping"); - t += config.dt + t += dt tstep += 1 state.U_hat₁[:, :, :, :] .= U_hat[:, :, :, :] @@ -60,13 +62,13 @@ function __integrate( for rk_step in 1:4 compute_rhs!(rk_step, parallel, K, config, U, U_hat, state) - println("du abs change: ", sum(abs.(state.dU))) + #println("du abs change: ", sum(abs.(state.dU))) if rk_step < 4 - U_hat[:, :, :, :] .= state.U_hat₀ .+ b[rk_step] * config.dt * state.dU + U_hat[:, :, :, :] .= state.U_hat₀ .+ b[rk_step] * dt * state.dU end - state.U_hat₁ .+= a[rk_step] * config.dt .* state.dU + state.U_hat₁ .+= a[rk_step] * dt .* state.dU end U_hat[:, :, :, :] .= state.U_hat₁ diff --git a/src/markers.jl b/src/markers.jl index 3c3098d..8955224 100644 --- a/src/markers.jl +++ b/src/markers.jl @@ -29,5 +29,10 @@ export AbstractWavenumbers abstract type AbstractWavenumbers end +# Config +export AbstractConfig + +abstract type AbstractConfig end + # end diff --git a/src/solver.jl b/src/solver.jl index a503c92..edc0822 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -2,7 +2,7 @@ module solver export curl!, cross!, compute_rhs! -using ..markers: AbstractParallel, AbstractState, AbstractWavenumbers +using ..markers: AbstractParallel, AbstractState, AbstractWavenumbers, AbstractConfig using ..fft: ifftn_mpi!, fftn_mpi! using ..mesh: Wavenumbers, WavenumbersGPU using ..state: State, StateGPU @@ -90,11 +90,11 @@ function compute_rhs!( rk_step::Int, parallel::P, K::WavenumbersGPU, - config::Config, + config::CFG, U::CuArray{Float64, 4}, U_hat::CuArray{ComplexF64, 4}, state::StateGPU -) where P <: AbstractParallel +) where P <: AbstractParallel where CFG <: AbstractConfig __compute_rhs!(rk_step, parallel, K, config, U, U_hat, state) end @@ -102,11 +102,11 @@ function compute_rhs!( rk_step::Int, parallel::P, K::Wavenumbers, - config::Config, + config::CFG, U::Array{Float64, 4}, U_hat::Array{ComplexF64, 4}, state::State -) where P <: AbstractParallel +) where P <: AbstractParallel where CFG <: AbstractConfig __compute_rhs!(rk_step, parallel, K, config, U, U_hat, state) end @@ -114,11 +114,11 @@ function __compute_rhs!( rk_step::Int, parallel::P, K::WAVE, - config::Config, + config::CFG, U::ARRAY, U_hat::FARRAY, state::STATE -) where P <: AbstractParallel where FARRAY <: AbstractArray{ComplexF64, 4} where ARRAY <: AbstractArray{Float64, 4} where STATE <: AbstractState where WAVE <: AbstractWavenumbers +) where P <: AbstractParallel where FARRAY <: AbstractArray{ComplexF64, 4} where ARRAY <: AbstractArray{Float64, 4} where STATE <: AbstractState where WAVE <: AbstractWavenumbers where CFG <: AbstractConfig if rk_step != 1 @views for i in 1:3 From da588f43368f7e79d4a4c32f3b279fabe959e510 Mon Sep 17 00:00:00 2001 From: VanillaBrooks Date: Wed, 10 May 2023 13:42:52 -0700 Subject: [PATCH 12/14] configuration type level validation parameters, move more things to gpu --- benches/cpu_gpu.jl | 14 ++++++------ justfile | 2 +- src/config.jl | 38 ++++++++++++++----------------- src/integrate.jl | 22 +++++++----------- src/markers.jl | 2 ++ src/solver.jl | 49 ++++++++++++++++++++-------------------- src/spectralGPU.jl | 2 +- src/state.jl | 33 +++++++++++++++++++++++---- test/cpu_singlethread.jl | 13 +++++------ test/cuda.jl | 9 ++++---- 10 files changed, 100 insertions(+), 84 deletions(-) diff --git a/benches/cpu_gpu.jl b/benches/cpu_gpu.jl index a537d0d..ece90e1 100644 --- a/benches/cpu_gpu.jl +++ b/benches/cpu_gpu.jl @@ -8,14 +8,14 @@ using BenchmarkTools ######### GPU benchmark ######### -N = 256 +N = 128 begin local parallel = markers.SingleThreadGPU() local K = mesh.wavenumbers_gpu(N) - local st = state.create_state_gpu(N, K) - local msh = mesh.new_mesh(N) local cfg = config.taylor_green_validation() + local st = state.create_state_gpu(N, K, cfg) + local msh = mesh.new_mesh(N) local ic = markers.TaylorGreen() local U = CuArray(zeros(N, N, N, 3)) @@ -31,7 +31,7 @@ begin $U_hat, ) - show(r) + display(r) println("\n"); end @@ -43,9 +43,9 @@ begin local parallel = markers.SingleThreadCPU() local K = mesh.wavenumbers(N) - local st = state.create_state(N, K) - local msh = mesh.new_mesh(N) local cfg = config.taylor_green_validation() + local st = state.create_state(N, K, cfg) + local msh = mesh.new_mesh(N) local ic = markers.TaylorGreen() local U = zeros(N, N, N, 3) @@ -61,5 +61,5 @@ begin $U_hat, ) - show(r) + display(r) end diff --git a/justfile b/justfile index 3dc0758..e7e1c85 100644 --- a/justfile +++ b/justfile @@ -5,7 +5,7 @@ build: sudo -E apptainer build spectralGPU.sif ./spectralGPU.apptainer test: - #julia ./test/cpu_singlethread.jl + julia ./test/cpu_singlethread.jl julia ./test/cuda.jl bench: diff --git a/src/config.jl b/src/config.jl index 29df6e7..5cf3104 100644 --- a/src/config.jl +++ b/src/config.jl @@ -1,49 +1,45 @@ module config -using ..markers: AbstractConfig +using ..markers: AbstractConfig, ProductionConfig, ValidationConfig export Config, create_config, taylor_green_validation, calculate_dt -struct Config <: AbstractConfig - ν::Float64 - re::Float64 - N::Int - time::Float64 -end +using CUDA -struct ConfigValidation <: AbstractConfig +struct Config{C <: AbstractConfig} ν::Float64 re::Float64 N::Int time::Float64 + mode::C end -function create_config(N::Int, re::Float64, time::Float64)::Config +function create_config(N::Int, re::Float64, time::Float64)::Config{ProductionConfig} ν = 1 / re - - CFL_NUM = 0.75 - - return Config(ν, re, N, time) + return Config(ν, re, N, time, ProductionConfig()) end -function taylor_green_validation()::ConfigValidation +function taylor_green_validation()::Config{ValidationConfig} N = 64 - re = 1600 + re = 1600. ν = 1 / re t = 0.1 - return ConfigValidation(ν, re, N, t) + return Config(ν, re, N, t, ValidationConfig()) end -function calculate_dt(velocity::XARRAY, cfg::Config)::Float64 where XARRAY <: AbstractArray{Float64, 4} - max_velo = maximum(velocity) - CFL = 0.75 +cfl_calc(max_velo::Float64, N::Int) = 0.75 / (N * max_velo) - return CFL / (cfg.N * max_velo) +function calculate_dt(velocity::XARRAY, cfg::Config{ProductionConfig})::Float64 where XARRAY <: AbstractArray{Float64, 4} + max_velo = maximum(velocity) + return cfl_calc(max_velo, cfg.N) end -function calculate_dt(velocity::XARRAY, cfg::ConfigValidation)::Float64 where XARRAY <: AbstractArray{Float64, 4} +# +# Validation dt cases +# +function calculate_dt(velocity::XARRAY, cfg::Config{ValidationConfig})::Float64 where XARRAY <: AbstractArray{Float64, 4} return 0.01 end diff --git a/src/integrate.jl b/src/integrate.jl index 3dc174f..ddb59d7 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -12,7 +12,7 @@ using CUDA function integrate( parallel::P, K::Wavenumbers, - config::CFG, + config::Config{CFG}, state::State, U::Array{Float64, 4}, U_hat::Array{ComplexF64, 4} @@ -23,7 +23,7 @@ end function integrate( parallel::P, K::WavenumbersGPU, - config::CFG, + config::Config{CFG}, state::StateGPU, U::CuArray{Float64, 4}, U_hat::CuArray{ComplexF64, 4} @@ -32,9 +32,9 @@ function integrate( end function __integrate( - parallel::P, - K::WAVE, - config::CFG, + parallel::P, + K::WAVE, + config::Config{CFG}, state::STATE, U::XARRAY, U_hat::FARRAY, @@ -42,12 +42,6 @@ function __integrate( t = 0 tstep = 0 - # https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#Examples - # a_j is equal to the canonical b_i - # b_i is equal to the canonical a_j - a = [1/6, 1/3, 1/3, 1/6] - b = [0.5, 0.5, 1.] - dt = calculate_dt(U, config) while t < config.time - 1e-8 @@ -60,15 +54,15 @@ function __integrate( state.U_hat₀[:, :, :, :] .= U_hat[:, :, :, :] for rk_step in 1:4 - compute_rhs!(rk_step, parallel, K, config, U, U_hat, state) + compute_rhs!(rk_step, parallel, K, U, U_hat, state) #println("du abs change: ", sum(abs.(state.dU))) if rk_step < 4 - U_hat[:, :, :, :] .= state.U_hat₀ .+ b[rk_step] * dt * state.dU + U_hat[:, :, :, :] .= state.U_hat₀ .+ state.b[rk_step] * dt * state.dU end - state.U_hat₁ .+= a[rk_step] * dt .* state.dU + state.U_hat₁ .+= state.a[rk_step] * dt .* state.dU end U_hat[:, :, :, :] .= state.U_hat₁ diff --git a/src/markers.jl b/src/markers.jl index 8955224..46cd687 100644 --- a/src/markers.jl +++ b/src/markers.jl @@ -33,6 +33,8 @@ abstract type AbstractWavenumbers end export AbstractConfig abstract type AbstractConfig end +struct ProductionConfig <: AbstractConfig end +struct ValidationConfig <: AbstractConfig end # end diff --git a/src/solver.jl b/src/solver.jl index edc0822..1d8d749 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -11,13 +11,13 @@ using ..config: Config using CUDA # calculate the curl of a fourier space array `input` and store the result in the x-space array `out` -@views function curl!( +@views function __curl!( parallel::P, - K::Wavenumbers, - input::Array{ComplexF64, 4} + K::WAVE, + input::FARRAY ; - out::Array{Float64, 4} -) where P <: AbstractParallel + out::XARRAY +) where P <: AbstractParallel where FARRAY <: AbstractArray{ComplexF64, 4} where XARRAY <: AbstractArray{Float64, 4} where WAVE <: AbstractWavenumbers j = complex(0, 1) ifftn_mpi!(parallel, K, j*(K[1].*input[:, :, :, 2] .- K[2].*input[:, :, :, 1]), out[:, :, :, 3]) @@ -27,6 +27,16 @@ using CUDA nothing end +@views function curl!( + parallel::P, + K::Wavenumbers, + input::Array{ComplexF64, 4} + ; + out::Array{Float64, 4} +) where P <: AbstractParallel + __curl!(parallel, K, input, out = out); +end + @views function curl!( parallel::P, K::WavenumbersGPU, @@ -34,13 +44,7 @@ end ; out::CuArray{Float64, 4} ) where P <: AbstractParallel - j = complex(0, 1) - - @CUDA.sync ifftn_mpi!(parallel, K,j * (CuArray(K[1]).*input[:, :, :, 2] .- CuArray(K[2]).*input[:, :, :, 1]), out[:, :, :, 3]) - @CUDA.sync ifftn_mpi!(parallel, K, j*(K[3].*input[:, :, :, 1] .- K[1].*input[:, :, :, 3]), out[:, :, :, 2]) - @CUDA.sync ifftn_mpi!(parallel, K, j*(K[2].*input[:, :, :, 3] .- K[3].*input[:, :, :, 2]), out[:, :, :, 1]) - - nothing + __curl!(parallel, K, input, out = out); end @@ -67,9 +71,9 @@ end ; out::CuArray{ComplexF64, 4} ) where P <: AbstractParallel - @CUDA.sync fftn_mpi!(parallel, a[:, :, :, 2].*b[:, :, :, 3] .- a[:, :, :, 3].*b[:, :, :, 2], out[:, :, :, 1]) - @CUDA.sync fftn_mpi!(parallel, a[:, :, :, 3].*b[:, :, :, 1] .- a[:, :, :, 1].*b[:, :, :, 3], out[:, :, :, 2]) - @CUDA.sync fftn_mpi!(parallel, a[:, :, :, 1].*b[:, :, :, 2] .- a[:, :, :, 2].*b[:, :, :, 1], out[:, :, :, 3]) + fftn_mpi!(parallel, a[:, :, :, 2].*b[:, :, :, 3] .- a[:, :, :, 3].*b[:, :, :, 2], out[:, :, :, 1]) + fftn_mpi!(parallel, a[:, :, :, 3].*b[:, :, :, 1] .- a[:, :, :, 1].*b[:, :, :, 3], out[:, :, :, 2]) + fftn_mpi!(parallel, a[:, :, :, 1].*b[:, :, :, 2] .- a[:, :, :, 2].*b[:, :, :, 1], out[:, :, :, 3]) nothing end @@ -90,35 +94,32 @@ function compute_rhs!( rk_step::Int, parallel::P, K::WavenumbersGPU, - config::CFG, U::CuArray{Float64, 4}, U_hat::CuArray{ComplexF64, 4}, state::StateGPU -) where P <: AbstractParallel where CFG <: AbstractConfig - __compute_rhs!(rk_step, parallel, K, config, U, U_hat, state) +) where P <: AbstractParallel + __compute_rhs!(rk_step, parallel, K, U, U_hat, state) end function compute_rhs!( rk_step::Int, parallel::P, K::Wavenumbers, - config::CFG, U::Array{Float64, 4}, U_hat::Array{ComplexF64, 4}, state::State -) where P <: AbstractParallel where CFG <: AbstractConfig - __compute_rhs!(rk_step, parallel, K, config, U, U_hat, state) +) where P <: AbstractParallel + __compute_rhs!(rk_step, parallel, K, U, U_hat, state) end function __compute_rhs!( rk_step::Int, parallel::P, K::WAVE, - config::CFG, U::ARRAY, U_hat::FARRAY, state::STATE -) where P <: AbstractParallel where FARRAY <: AbstractArray{ComplexF64, 4} where ARRAY <: AbstractArray{Float64, 4} where STATE <: AbstractState where WAVE <: AbstractWavenumbers where CFG <: AbstractConfig +) where P <: AbstractParallel where FARRAY <: AbstractArray{ComplexF64, 4} where ARRAY <: AbstractArray{Float64, 4} where STATE <: AbstractState where WAVE <: AbstractWavenumbers if rk_step != 1 @views for i in 1:3 @@ -143,7 +144,7 @@ function __compute_rhs!( # dU is now the Eulerian term / Nonlinear term in the NSE # now calculate the diffusion term - state.dU .-= config.ν .* state.K² .* U_hat + state.dU .-= state.ν .* state.K² .* U_hat end # diff --git a/src/spectralGPU.jl b/src/spectralGPU.jl index c2cd6c2..5133ac5 100644 --- a/src/spectralGPU.jl +++ b/src/spectralGPU.jl @@ -4,8 +4,8 @@ include("./markers.jl") include("./mesh.jl") include("./fft.jl") include("./initial_condition.jl") -include("./state.jl") include("./config.jl") +include("./state.jl") include("./solver.jl") include("./integrate.jl") diff --git a/src/state.jl b/src/state.jl index a70eb36..0f36c28 100644 --- a/src/state.jl +++ b/src/state.jl @@ -1,7 +1,8 @@ module state using ..mesh: Wavenumbers -using ..markers: AbstractState, AbstractWavenumbers +using ..markers: AbstractState, AbstractWavenumbers, AbstractConfig +using ..config: Config using CUDA @@ -17,6 +18,9 @@ struct State <: AbstractState K²::Array{Int, 3} K_over_K²::Array{Float64, 4} wavenumber_product_tmp::Array{ComplexF64, 4} + ν::Float64 + a::Vector{Float64} + b::Vector{Float64} end struct StateGPU <: AbstractState @@ -29,9 +33,14 @@ struct StateGPU <: AbstractState K²::CuArray{Int, 3} K_over_K²::CuArray{Float64, 4} wavenumber_product_tmp::CuArray{ComplexF64, 4} + ν::CuVector{Float64} + # RK integration constants + a::Vector{Float64} + # RK integration constants + b::Vector{Float64} end -function create_state(N::Int, K::WAVE)::State where WAVE <: AbstractWavenumbers +function create_state(N::Int, K::WAVE, config::Config{CFG})::State where WAVE <: AbstractWavenumbers where CFG <: AbstractConfig dU = ComplexF64.(zeros(K.kn, N, N, 3)) U_hat₁ = ComplexF64.(zeros(K.kn, N, N, 3)) U_hat₀ = ComplexF64.(zeros(K.kn, N, N, 3)) @@ -58,6 +67,12 @@ function create_state(N::Int, K::WAVE)::State where WAVE <: AbstractWavenumbers K_over_K² ./= K²_nonzero + # https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#Examples + # a_j is equal to the canonical b_i + # b_i is equal to the canonical a_j + a = [1/6, 1/3, 1/3, 1/6] + b = [0.5, 0.5, 1.] + State( dU, U_hat₁, @@ -67,12 +82,15 @@ function create_state(N::Int, K::WAVE)::State where WAVE <: AbstractWavenumbers P_hat, K², K_over_K², - wavenumber_product_tmp + wavenumber_product_tmp, + config.ν, + a, + b ) end -function create_state_gpu(N::Int, K::WAVE)::StateGPU where WAVE <: AbstractWavenumbers - cpu_state = create_state(N, K) +function create_state_gpu(N::Int, K::WAVE, config::Config{CFG})::StateGPU where WAVE <: AbstractWavenumbers where CFG <: AbstractConfig + cpu_state = create_state(N, K, config) dealias::Array{Int8, 3} = Array(cpu_state.dealias) return StateGPU( @@ -86,6 +104,11 @@ function create_state_gpu(N::Int, K::WAVE)::StateGPU where WAVE <: AbstractWaven CuArray(cpu_state.K²), CuArray(cpu_state.K_over_K²), CuArray(cpu_state.wavenumber_product_tmp), + CuVector([cpu_state.ν]), + #CuVector(cpu_state.a), + #CuVector(cpu_state.b), + cpu_state.a, + cpu_state.b, ) end diff --git a/test/cpu_singlethread.jl b/test/cpu_singlethread.jl index 73d8503..7a4619d 100644 --- a/test/cpu_singlethread.jl +++ b/test/cpu_singlethread.jl @@ -42,9 +42,9 @@ end re = 40. K = mesh.wavenumbers(N) - st = state.create_state(N, K) - msh = mesh.new_mesh(N) cfg = config.create_config(N, re, 1.0) + st = state.create_state(N, K, cfg) + msh = mesh.new_mesh(N) U = zeros(N, N, N, 3) U_hat = ComplexF64.(zeros(K.kn, N, N, 3)) @@ -68,7 +68,6 @@ end 2, parallel, K, - cfg, U, U_hat, st @@ -85,9 +84,9 @@ end time = 0.05 K = mesh.wavenumbers(N) - st = state.create_state(N, K) - msh = mesh.new_mesh(N) cfg = config.create_config(N, re, time) + st = state.create_state(N, K, cfg) + msh = mesh.new_mesh(N) U = zeros(N, N, N, 3) U_hat = ComplexF64.(zeros(K.kn, N, N, 3)) @@ -111,9 +110,9 @@ end N = 64 K = mesh.wavenumbers(N) - st = state.create_state(N, K) - msh = mesh.new_mesh(N) cfg = config.taylor_green_validation() + st = state.create_state(N, K, cfg) + msh = mesh.new_mesh(N) ic = markers.TaylorGreen() U = zeros(N, N, N, 3) diff --git a/test/cuda.jl b/test/cuda.jl index d9d61fa..15b7f85 100644 --- a/test/cuda.jl +++ b/test/cuda.jl @@ -6,8 +6,9 @@ using CUDA @testset "state.jl" begin N = 64 K = mesh.wavenumbers_gpu(N) + cfg = config.create_config(N, 40., 0.00001) @test begin - st::state.StateGPU = state.create_state_gpu(N, K) + st::state.StateGPU = state.create_state_gpu(N, K, cfg) true end end @@ -85,7 +86,8 @@ end re = 40. K = mesh.wavenumbers_gpu(N) - st = state.create_state_gpu(N, K) + cfg = config.create_config(N, 40., 0.00001) + st = state.create_state_gpu(N, K, cfg) msh = mesh.new_mesh(N) cfg = config.create_config(N, re, 1.0) @@ -110,7 +112,6 @@ end 2, parallel, K, - cfg, U, U_hat, st @@ -125,9 +126,9 @@ end N = 64 K = mesh.wavenumbers_gpu(N) - st = state.create_state_gpu(N, K) msh = mesh.new_mesh(N) cfg = config.taylor_green_validation() + st = state.create_state_gpu(N, K, cfg) ic = markers.TaylorGreen() U = CuArray(zeros(N, N, N, 3)) From 605420f8bed94771a0cd1e5a21d0320370713bc6 Mon Sep 17 00:00:00 2001 From: VanillaBrooks Date: Wed, 10 May 2023 14:38:04 -0700 Subject: [PATCH 13/14] cross, curl, fft benchmarks --- benches/cpu_gpu.jl | 2 + benches/cross.jl | 39 ++++++++++++++++++ benches/curl.jl | 37 +++++++++++++++++ benches/fft.jl | 99 ++++++++++++++++++++++++++++++++++++++++++++++ justfile | 3 ++ src/solver.jl | 32 +++++++++------ 6 files changed, 200 insertions(+), 12 deletions(-) create mode 100644 benches/cross.jl create mode 100644 benches/curl.jl create mode 100644 benches/fft.jl diff --git a/benches/cpu_gpu.jl b/benches/cpu_gpu.jl index ece90e1..754beba 100644 --- a/benches/cpu_gpu.jl +++ b/benches/cpu_gpu.jl @@ -31,6 +31,7 @@ begin $U_hat, ) + println("GPU full solver integration") display(r) println("\n"); end @@ -61,5 +62,6 @@ begin $U_hat, ) + println("CPU full solver integration") display(r) end diff --git a/benches/cross.jl b/benches/cross.jl new file mode 100644 index 0000000..c45df19 --- /dev/null +++ b/benches/cross.jl @@ -0,0 +1,39 @@ +include("../src/spectralGPU.jl"); +using .spectralGPU: mesh, fft, markers, initial_condition, state, config, solver, Integrate + +using CUDA +using BenchmarkTools + +N = 64 + +begin + local parallel = markers.SingleThreadGPU() + local re = 40. + local K = mesh.wavenumbers_gpu(N) + local cfg = config.create_config(N, re, 0.00001) + local st = state.create_state_gpu(N, K, cfg) + local U = CuArray(zeros(N, N, N, 3)) + local U_hat = CuArray(ComplexF64.(zeros(K.kn, N, N, 3))) + + r = @benchmark solver.cross!($parallel, $U, $st.curl; out = $st.dU) + + println("GPU cross") + display(r) + println("") +end + +begin + local parallel = markers.SingleThreadCPU() + local re = 40. + local K = mesh.wavenumbers(N) + local cfg = config.create_config(N, re, 0.00001) + local st = state.create_state(N, K, cfg) + local U = zeros(N, N, N, 3) + local U_hat = ComplexF64.(zeros(K.kn, N, N, 3)) + + r = @benchmark solver.cross!($parallel, $U, $st.curl; out = $st.dU) + + println("CPU cross") + display(r) + println("") +end diff --git a/benches/curl.jl b/benches/curl.jl new file mode 100644 index 0000000..9f526db --- /dev/null +++ b/benches/curl.jl @@ -0,0 +1,37 @@ +include("../src/spectralGPU.jl"); +using .spectralGPU: mesh, fft, markers, initial_condition, state, config, solver, Integrate + +using CUDA +using BenchmarkTools + +N = 64 + +begin + local parallel = markers.SingleThreadGPU() + local re = 40. + local K = mesh.wavenumbers_gpu(N) + local cfg = config.create_config(N, re, 0.00001) + local st = state.create_state_gpu(N, K, cfg) + local U_hat = CuArray(ComplexF64.(zeros(K.kn, N, N, 3))) + + r = @benchmark solver.curl!($parallel, $K, $U_hat; out = $st.curl) + + println("GPU curl") + display(r) + println("") +end + +begin + local parallel = markers.SingleThreadCPU() + local re = 40. + local K = mesh.wavenumbers(N) + local cfg = config.create_config(N, re, 0.00001) + local st = state.create_state(N, K, cfg) + local U_hat = ComplexF64.(zeros(K.kn, N, N, 3)) + + r = @benchmark solver.curl!($parallel, $K, $U_hat; out = $st.curl) + + println("CPU curl") + display(r) + println("") +end diff --git a/benches/fft.jl b/benches/fft.jl new file mode 100644 index 0000000..e612b1f --- /dev/null +++ b/benches/fft.jl @@ -0,0 +1,99 @@ +include("../src/spectralGPU.jl"); +using .spectralGPU: mesh, fft, markers, initial_condition, state, config, solver, Integrate + +using CUDA +using BenchmarkTools + +######### +######### GPU (forward) benchmark +######### + +N = 64 + +begin + parallel = markers.SingleThreadGPU() + + K = mesh.wavenumbers_gpu(N) + + U = CuArray(zeros(N, N, N, 3)) + U_hat = CuArray(ComplexF64.(zeros(K.kn, N, N, 3))) + + msh = mesh.new_mesh(N) + ic = markers.TaylorGreen() + initial_condition.setup_initial_condition(parallel, ic, msh, U, U_hat) + + r = @benchmark @views fft.fftn_mpi!($parallel, $U[:, :, :, 1], $U_hat[:, :, :, 1]) + + println("gpu forward FFT") + display(r) + println("") +end + +######### +######### CPU (forward) benchmark +######### + +begin + parallel = markers.SingleThreadCPU() + + K = mesh.wavenumbers(N) + + U = zeros(N, N, N, 3) + U_hat = ComplexF64.(zeros(K.kn, N, N, 3)) + + msh = mesh.new_mesh(N) + ic = markers.TaylorGreen() + initial_condition.setup_initial_condition(parallel, ic, msh, U, U_hat) + + r = @benchmark @views fft.fftn_mpi!($parallel, $U[:, :, :, 1], $U_hat[:, :, :, 1]) + + println("cpu forward FFT") + display(r) + println("") +end + +######### +######### GPU (inverse) benchmark +######### + +begin + parallel = markers.SingleThreadGPU() + + K = mesh.wavenumbers_gpu(N) + + U = CuArray(zeros(N, N, N, 3)) + U_hat = CuArray(ComplexF64.(zeros(K.kn, N, N, 3))) + + msh = mesh.new_mesh(N) + ic = markers.TaylorGreen() + initial_condition.setup_initial_condition(parallel, ic, msh, U, U_hat) + + r = @benchmark @views fft.ifftn_mpi!($parallel, $K, $U_hat[:, :, :, 1], $U[:, :, :, 1]) + + println("gpu inverse FFT") + display(r) + println("") +end + +######### +######### CPU (inverse) benchmark +######### + +begin + parallel = markers.SingleThreadCPU() + + K = mesh.wavenumbers(N) + + U = zeros(N, N, N, 3) + U_hat = ComplexF64.(zeros(K.kn, N, N, 3)) + + msh = mesh.new_mesh(N) + ic = markers.TaylorGreen() + initial_condition.setup_initial_condition(parallel, ic, msh, U, U_hat) + + r = @benchmark @views fft.ifftn_mpi!($parallel, $K, $U_hat[:, :, :, 1], $U[:, :, :, 1]) + + println("cpu inverse FFT") + display(r) + println("") +end diff --git a/justfile b/justfile index e7e1c85..6db8691 100644 --- a/justfile +++ b/justfile @@ -9,4 +9,7 @@ test: julia ./test/cuda.jl bench: + #julia ./benches/fft.jl + #julia ./benches/curl.jl + #julia ./benches/cross.jl julia ./benches/cpu_gpu.jl diff --git a/src/solver.jl b/src/solver.jl index 1d8d749..e6d6381 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -27,7 +27,7 @@ using CUDA nothing end -@views function curl!( +function curl!( parallel::P, K::Wavenumbers, input::Array{ComplexF64, 4} @@ -37,7 +37,7 @@ end __curl!(parallel, K, input, out = out); end -@views function curl!( +function curl!( parallel::P, K::WavenumbersGPU, input::CuArray{ComplexF64, 4} @@ -50,32 +50,40 @@ end # take the cross product of X-Space a,b vectors # stores the result in (fourier space) `out` and then return the -@views function cross!( +@views function __cross!( + parallel::P, + a::XARRAY, + b::XARRAY + ; + out::FARRAY +) where P <: AbstractParallel where FARRAY <: AbstractArray{ComplexF64, 4} where XARRAY <: AbstractArray{Float64, 4} + fftn_mpi!(parallel, a[:, :, :, 2].*b[:, :, :, 3] .- a[:, :, :, 3].*b[:, :, :, 2], out[:, :, :, 1]) + fftn_mpi!(parallel, a[:, :, :, 3].*b[:, :, :, 1] .- a[:, :, :, 1].*b[:, :, :, 3], out[:, :, :, 2]) + fftn_mpi!(parallel, a[:, :, :, 1].*b[:, :, :, 2] .- a[:, :, :, 2].*b[:, :, :, 1], out[:, :, :, 3]) + + nothing +end + +function cross!( parallel::P, a::Array{Float64, 4}, b::Array{Float64, 4}, ; out::Array{ComplexF64, 4} ) where P <: AbstractParallel - fftn_mpi!(parallel, a[:, :, :, 2].*b[:, :, :, 3] .- a[:, :, :, 3].*b[:, :, :, 2], out[:, :, :, 1]) - fftn_mpi!(parallel, a[:, :, :, 3].*b[:, :, :, 1] .- a[:, :, :, 1].*b[:, :, :, 3], out[:, :, :, 2]) - fftn_mpi!(parallel, a[:, :, :, 1].*b[:, :, :, 2] .- a[:, :, :, 2].*b[:, :, :, 1], out[:, :, :, 3]) + __cross!(parallel, a, b, out = out); nothing end -@views function cross!( +function cross!( parallel::P, a::CuArray{Float64, 4}, b::CuArray{Float64, 4}, ; out::CuArray{ComplexF64, 4} ) where P <: AbstractParallel - fftn_mpi!(parallel, a[:, :, :, 2].*b[:, :, :, 3] .- a[:, :, :, 3].*b[:, :, :, 2], out[:, :, :, 1]) - fftn_mpi!(parallel, a[:, :, :, 3].*b[:, :, :, 1] .- a[:, :, :, 1].*b[:, :, :, 3], out[:, :, :, 2]) - fftn_mpi!(parallel, a[:, :, :, 1].*b[:, :, :, 2] .- a[:, :, :, 2].*b[:, :, :, 1], out[:, :, :, 3]) - - nothing + __cross!(parallel, a, b, out = out); end function wavenumber_product!(arr::Array{ComplexF64, 3}, K::Wavenumbers; out::Array{ComplexF64, 4}) From 31cd0805689fcd63cecf4f61802bb76c8eb2548f Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Thu, 11 May 2023 01:21:07 +0000 Subject: [PATCH 14/14] CompatHelper: add new compat entry for CUDA at version 4, (keep existing compat) --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 7ad35e8..c5f6515 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" LazyGrids = "7031d0ef-c40d-4431-b2f8-61a8d2f650db" [compat] +CUDA = "4" julia = "1" [extras]