diff --git a/.gitignore b/.gitignore index 95d5975..aeecb12 100644 --- a/.gitignore +++ b/.gitignore @@ -40,6 +40,19 @@ docs/theory/solver_2x_sdf/out/ # ноутбук kbc_lbm.ipynb их встраивает, и без них он нечитаем на свежем клоне # (перегенерация требует GPU-сервера). См. заметку о них в README ветки research. +# --------------------------------------------------------------------------- +# Rust (docs/theory/2d_solver — решатель KBC на Rust) +# --------------------------------------------------------------------------- +# Каталог сборки cargo: объектники, зависимости, инкрементальный кэш. +target/ +# Cargo.lock для бинарного крейта обычно коммитят (воспроизводимость прогонов) — +# он НЕ игнорируется намеренно. + +# Артефакты прогонов решателя: гифки, ряды, сводки. Воспроизводятся запуском. +docs/theory/2d_solver/*.gif +docs/theory/2d_solver/*.csv +docs/theory/2d_solver/out/ + # --------------------------------------------------------------------------- # Редакторы и ОС # --------------------------------------------------------------------------- diff --git a/docs/theory/2d_solver/Cargo.lock b/docs/theory/2d_solver/Cargo.lock new file mode 100644 index 0000000..40f9cd9 --- /dev/null +++ b/docs/theory/2d_solver/Cargo.lock @@ -0,0 +1,1308 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 3 + +[[package]] +name = "android_system_properties" +version = "0.1.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ae221649c9976a6f6c56ae1facf410f3ddb33cc661c4b7b61020a912d4237fbc" +dependencies = [ + "libc", +] + +[[package]] +name = "anstream" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "824a212faf96e9acacdbd09febd34438f8f711fb84e09a8916013cd7815ca28d" +dependencies = [ + "anstyle", + "anstyle-parse", + "anstyle-query", + "anstyle-wincon", + "colorchoice", + "is_terminal_polyfill", + "utf8parse", +] + +[[package]] +name = "anstyle" +version = "1.0.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "940b3a0ca603d1eade50a4846a2afffd5ef57a9feac2c0e2ec2e14f9ead76000" + +[[package]] +name = "anstyle-parse" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "52ce7f38b242319f7cabaa6813055467063ecdc9d355bbb4ce0c68908cd8130e" +dependencies = [ + "utf8parse", +] + +[[package]] +name = "anstyle-query" +version = "1.1.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "40c48f72fd53cd289104fc64099abca73db4166ad86ea0b4341abe65af83dadc" +dependencies = [ + "windows-sys", +] + +[[package]] +name = "anstyle-wincon" +version = "3.0.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "291e6a250ff86cd4a820112fb8898808a366d8f9f58ce16d1f538353ad55747d" +dependencies = [ + "anstyle", + "once_cell_polyfill", + "windows-sys", +] + +[[package]] +name = "arrayvec" +version = "0.7.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d3fb67a6e08acf24fdeccbac2cb6ac4305825bd1f117462e0e6f2f193345ad56" + +[[package]] +name = "ash" +version = "0.38.0+1.3.281" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0bb44936d800fea8f016d7f2311c6a4f97aebd5dc86f09906139ec848cf3a46f" +dependencies = [ + "libloading", +] + +[[package]] +name = "bit-set" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f0481a0e032742109b1133a095184ee93d88f3dc9e0d28a5d033dc77a073f44f" +dependencies = [ + "bit-vec", +] + +[[package]] +name = "bit-vec" +version = "0.7.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d2c54ff287cfc0a34f38a6b832ea1bd8e448a330b3e40a50859e6488bee07f22" + +[[package]] +name = "bitflags" +version = "1.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" + +[[package]] +name = "bitflags" +version = "2.13.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b588b76d00fde79687d7646a9b5bdf3cc0f655e0bbd080335a95d7e96f3587da" + +[[package]] +name = "block" +version = "0.1.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0d8c1fef690941d3e7788d328517591fecc684c084084702d6ff1641e993699a" + +[[package]] +name = "bumpalo" +version = "3.20.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72f5acc6cb2ba439de613abc23857ec3d78374d8ed5ac84e9d11336e87da8649" + +[[package]] +name = "bytemuck" +version = "1.25.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95832e849adfb21180ccb6826a99da14e5d266ae5c2e668e1602cf234f153797" +dependencies = [ + "bytemuck_derive", +] + +[[package]] +name = "bytemuck_derive" +version = "1.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fc0e56a716f1e132ff6bf4bdac1c944a3fcdc1cae65f70a4a2a1ac3b401d2d1f" +dependencies = [ + "proc-macro2", + "quote", + "syn 3.0.3", +] + +[[package]] +name = "cfg-if" +version = "1.0.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801" + +[[package]] +name = "cfg_aliases" +version = "0.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fd16c4719339c4530435d38e511904438d07cce7950afa3718a84ac36c10e89e" + +[[package]] +name = "clap" +version = "4.6.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "473c7e07f409a8d772161724aa8db6a765a2532a70f9667eeb7b49d3d02fbdca" +dependencies = [ + "clap_builder", + "clap_derive", +] + +[[package]] +name = "clap_builder" +version = "4.6.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7b48fea5a88e9ae728a2dcbedbfc0e730f7d60da42e1cb049a83c9fb8b789889" +dependencies = [ + "anstream", + "anstyle", + "clap_lex", + "strsim", +] + +[[package]] +name = "clap_derive" +version = "4.6.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d012d2b9d65aca7f18f4d9878a045bc17899bba951561ba5ec3c2ba1eed9a061" +dependencies = [ + "heck", + "proc-macro2", + "quote", + "syn 3.0.3", +] + +[[package]] +name = "clap_lex" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c8d4a3bb8b1e0c1050499d1815f5ab16d04f0959b233085fb31653fbfc9d98f9" + +[[package]] +name = "codespan-reporting" +version = "0.11.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3538270d33cc669650c4b093848450d380def10c331d38c768e34cac80576e6e" +dependencies = [ + "termcolor", + "unicode-width", +] + +[[package]] +name = "color_quant" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3d7b894f5411737b7867f4827955924d7c254fc9f4d91a6aad6b097804b1018b" + +[[package]] +name = "colorchoice" +version = "1.0.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1d07550c9036bf2ae0c684c4297d503f838287c83c53686d05370d0e139ae570" + +[[package]] +name = "com" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7e17887fd17353b65b1b2ef1c526c83e26cd72e74f598a8dc1bee13a48f3d9f6" +dependencies = [ + "com_macros", +] + +[[package]] +name = "com_macros" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d375883580a668c7481ea6631fc1a8863e33cc335bf56bfad8d7e6d4b04b13a5" +dependencies = [ + "com_macros_support", + "proc-macro2", + "syn 1.0.109", +] + +[[package]] +name = "com_macros_support" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ad899a1087a9296d5644792d7cb72b8e34c1bec8e7d4fbc002230169a6e8710c" +dependencies = [ + "proc-macro2", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "core-foundation" +version = "0.9.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91e195e091a93c46f7102ec7818a2aa394e1e1771c3ab4825963fa03e45afb8f" +dependencies = [ + "core-foundation-sys", + "libc", +] + +[[package]] +name = "core-foundation-sys" +version = "0.8.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "773648b94d0e5d620f64f280777445740e61fe701025087ec8b57f45c791888b" + +[[package]] +name = "core-graphics-types" +version = "0.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "45390e6114f68f718cc7a830514a96f903cccd70d02a8f6d9f643ac4ba45afaf" +dependencies = [ + "bitflags 1.3.2", + "core-foundation", + "libc", +] + +[[package]] +name = "crossbeam-deque" +version = "0.8.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5181e0de7b61eb03a81e347d6dd8797bae9da5146707b51077e2d71a54ec0ceb" +dependencies = [ + "crossbeam-epoch", + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-epoch" +version = "0.9.20" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2d6914041f254d6e9176c01941b21115dcfb7089e55135a35411081bd106ef3f" +dependencies = [ + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-utils" +version = "0.8.22" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "61803da095bee82a81bb1a452ecc25d3b2f1416d1897eb86430c6159ef717c17" + +[[package]] +name = "d3d12" +version = "22.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bdbd1f579714e3c809ebd822c81ef148b1ceaeb3d535352afc73fd0c4c6a0017" +dependencies = [ + "bitflags 2.13.1", + "libloading", + "winapi", +] + +[[package]] +name = "document-features" +version = "0.2.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d4b8a88685455ed29a21542a33abd9cb6510b6b129abadabdcef0f4c55bc8f61" +dependencies = [ + "litrs", +] + +[[package]] +name = "either" +version = "1.17.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9e5e8f6c15a24b9a3ee5efec809ccd006d3b30e8b3bb63c39af737c7f87daa1d" + +[[package]] +name = "equivalent" +version = "1.0.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "877a4ace8713b0bcf2a4e7eec82529c029f1d0619886d18145fea96c3ffe5c0f" + +[[package]] +name = "foldhash" +version = "0.1.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d9c4f5dac5e15c24eb999c26181a6ca40b39fe946cbe4c263c7209467bc83af2" + +[[package]] +name = "foreign-types" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d737d9aa519fb7b749cbc3b962edcf310a8dd1f4b67c91c4f83975dbdd17d965" +dependencies = [ + "foreign-types-macros", + "foreign-types-shared", +] + +[[package]] +name = "foreign-types-macros" +version = "0.2.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ea5190182e6915eb873ddbc16e23b711b6eb1f9c00a0d0a3a91b5f6228475225" +dependencies = [ + "proc-macro2", + "quote", + "syn 3.0.3", +] + +[[package]] +name = "foreign-types-shared" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "aa9a19cbb55df58761df49b23516a86d432839add4af60fc256da840f66ed35b" + +[[package]] +name = "futures-core" +version = "0.3.34" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "92d699e522242e69e3003b94ecc1f960f3a5e015aa7c5d7486e65ad01dd94f5e" + +[[package]] +name = "futures-task" +version = "0.3.34" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cd417de3d1d015fc3bfd2b1ea46dfc7bab72ef86f1cc7cc9c78e728b34a6d1fd" + +[[package]] +name = "futures-util" +version = "0.3.34" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0d50a92467f8ba5dd6e3ee5d4bd04d73ab2e4e1c44474a0674821dfce14b79bc" +dependencies = [ + "futures-core", + "futures-task", + "pin-project-lite", + "slab", +] + +[[package]] +name = "gif" +version = "0.13.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4ae047235e33e2829703574b54fdec96bfbad892062d97fed2f76022287de61b" +dependencies = [ + "color_quant", + "weezl", +] + +[[package]] +name = "gl_generator" +version = "0.14.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1a95dfc23a2b4a9a2f5ab41d194f8bfda3cabec42af4e39f08c339eb2a0c124d" +dependencies = [ + "khronos_api", + "log", + "xml-rs", +] + +[[package]] +name = "glow" +version = "0.13.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bd348e04c43b32574f2de31c8bb397d96c9fcfa1371bd4ca6d8bdc464ab121b1" +dependencies = [ + "js-sys", + "slotmap", + "wasm-bindgen", + "web-sys", +] + +[[package]] +name = "glutin_wgl_sys" +version = "0.6.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2c4ee00b289aba7a9e5306d57c2d05499b2e5dc427f84ac708bd2c090212cf3e" +dependencies = [ + "gl_generator", +] + +[[package]] +name = "gpu-alloc" +version = "0.6.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "45cf04b2726f02df5508c6de726acdc90cdf97ac771a9a0ffd8ba10a6e696bf9" +dependencies = [ + "bitflags 2.13.1", + "gpu-alloc-types", +] + +[[package]] +name = "gpu-alloc-types" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b2bbed164dd10ed526c2e4fe3e721ca4a71c61730e5aafac6844b417b3227058" +dependencies = [ + "bitflags 2.13.1", +] + +[[package]] +name = "gpu-allocator" +version = "0.26.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fdd4240fc91d3433d5e5b0fc5b67672d771850dc19bbee03c1381e19322803d7" +dependencies = [ + "log", + "presser", + "thiserror", + "winapi", + "windows", +] + +[[package]] +name = "gpu-descriptor" +version = "0.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b89c83349105e3732062a895becfc71a8f921bb71ecbbdd8ff99263e3b53a0ca" +dependencies = [ + "bitflags 2.13.1", + "gpu-descriptor-types", + "hashbrown 0.15.5", +] + +[[package]] +name = "gpu-descriptor-types" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fdf242682df893b86f33a73828fb09ca4b2d3bb6cc95249707fc684d27484b91" +dependencies = [ + "bitflags 2.13.1", +] + +[[package]] +name = "hashbrown" +version = "0.15.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9229cfe53dfd69f0609a49f65461bd93001ea1ef889cd5529dd176593f5338a1" +dependencies = [ + "foldhash", +] + +[[package]] +name = "hashbrown" +version = "0.17.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ed5909b6e89a2db4456e54cd5f673791d7eca6732202bbf2a9cc504fe2f9b84a" + +[[package]] +name = "hassle-rs" +version = "0.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "af2a7e73e1f34c48da31fb668a907f250794837e08faa144fd24f0b8b741e890" +dependencies = [ + "bitflags 2.13.1", + "com", + "libc", + "libloading", + "thiserror", + "widestring", + "winapi", +] + +[[package]] +name = "heck" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" + +[[package]] +name = "hexf-parse" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dfa686283ad6dd069f105e5ab091b04c62850d3e4cf5d67debad1933f55023df" + +[[package]] +name = "indexmap" +version = "2.14.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d466e9454f08e4a911e14806c24e16fba1b4c121d1ea474396f396069cf949d9" +dependencies = [ + "equivalent", + "hashbrown 0.17.1", +] + +[[package]] +name = "is_terminal_polyfill" +version = "1.70.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a6cb138bb79a146c1bd460005623e142ef0181e3d0219cb493e02f7d08a35695" + +[[package]] +name = "jni-sys" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "41a652e1f9b6e0275df1f15b32661cf0d4b78d4d87ddec5e0c3c20f097433258" +dependencies = [ + "jni-sys 0.4.1", +] + +[[package]] +name = "jni-sys" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c6377a88cb3910bee9b0fa88d4f42e1d2da8e79915598f65fb0c7ee14c878af2" +dependencies = [ + "jni-sys-macros", +] + +[[package]] +name = "jni-sys-macros" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "38c0b942f458fe50cdac086d2f946512305e5631e720728f2a61aabcd47a6264" +dependencies = [ + "quote", + "syn 2.0.119", +] + +[[package]] +name = "js-sys" +version = "0.3.104" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0e0c1080212aad755ea003d18543e8768dd432c48819efd73a7bf1e39b7a5a3a" +dependencies = [ + "cfg-if", + "futures-util", + "wasm-bindgen", +] + +[[package]] +name = "kbc2d" +version = "0.1.0" +dependencies = [ + "bytemuck", + "clap", + "gif", + "pollster", + "rayon", + "wgpu", +] + +[[package]] +name = "khronos-egl" +version = "6.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6aae1df220ece3c0ada96b8153459b67eebe9ae9212258bb0134ae60416fdf76" +dependencies = [ + "libc", + "libloading", + "pkg-config", +] + +[[package]] +name = "khronos_api" +version = "3.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e2db585e1d738fc771bf08a151420d3ed193d9d895a36df7f6f8a9456b911ddc" + +[[package]] +name = "libc" +version = "0.2.189" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3eaf3ede3fee6db1a4c2ee091bf8a8b4dccdc6d17f656fb07896ee72867612f2" + +[[package]] +name = "libloading" +version = "0.8.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d7c4b02199fee7c5d21a5ae7d8cfa79a6ef5bb2fc834d6e9058e89c825efdc55" +dependencies = [ + "cfg-if", + "windows-link", +] + +[[package]] +name = "litrs" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "11d3d7f243d5c5a8b9bb5d6dd2b1602c0cb0b9db1621bafc7ed66e35ff9fe092" + +[[package]] +name = "lock_api" +version = "0.4.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "224399e74b87b5f3557511d98dff8b14089b3dadafcab6bb93eab67d3aace965" +dependencies = [ + "scopeguard", +] + +[[package]] +name = "log" +version = "0.4.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0ceec5bc11778974d1bcb055b18002eba7f4b3518b6a0081b3af5f21666da9ad" + +[[package]] +name = "malloc_buf" +version = "0.0.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "62bb907fe88d54d8d9ce32a3cceab4218ed2f6b7d35617cafe9adf84e43919cb" +dependencies = [ + "libc", +] + +[[package]] +name = "metal" +version = "0.29.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7ecfd3296f8c56b7c1f6fbac3c71cefa9d78ce009850c45000015f206dc7fa21" +dependencies = [ + "bitflags 2.13.1", + "block", + "core-graphics-types", + "foreign-types", + "log", + "objc", + "paste", +] + +[[package]] +name = "naga" +version = "22.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8bd5a652b6faf21496f2cfd88fc49989c8db0825d1f6746b1a71a6ede24a63ad" +dependencies = [ + "arrayvec", + "bit-set", + "bitflags 2.13.1", + "cfg_aliases", + "codespan-reporting", + "hexf-parse", + "indexmap", + "log", + "rustc-hash", + "spirv", + "termcolor", + "thiserror", + "unicode-xid", +] + +[[package]] +name = "ndk-sys" +version = "0.5.0+25.2.9519653" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8c196769dd60fd4f363e11d948139556a344e79d451aeb2fa2fd040738ef7691" +dependencies = [ + "jni-sys 0.3.1", +] + +[[package]] +name = "objc" +version = "0.2.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "915b1b472bc21c53464d6c8461c9d3af805ba1ef837e1cac254428f4a77177b1" +dependencies = [ + "malloc_buf", +] + +[[package]] +name = "once_cell" +version = "1.21.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9f7c3e4beb33f85d45ae3e3a1792185706c8e16d043238c593331cc7cd313b50" + +[[package]] +name = "once_cell_polyfill" +version = "1.70.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "384b8ab6d37215f3c5301a95a4accb5d64aa607f1fcb26a11b5303878451b4fe" + +[[package]] +name = "parking_lot" +version = "0.12.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "93857453250e3077bd71ff98b6a65ea6621a19bb0f559a85248955ac12c45a1a" +dependencies = [ + "lock_api", + "parking_lot_core", +] + +[[package]] +name = "parking_lot_core" +version = "0.9.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2621685985a2ebf1c516881c026032ac7deafcda1a2c9b7850dc81e3dfcb64c1" +dependencies = [ + "cfg-if", + "libc", + "redox_syscall", + "smallvec", + "windows-link", +] + +[[package]] +name = "paste" +version = "1.0.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" + +[[package]] +name = "pin-project-lite" +version = "0.2.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a89322df9ebe1c1578d689c92318e070967d1042b512afbe49518723f4e6d5cd" + +[[package]] +name = "pkg-config" +version = "0.3.34" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f6b464fbc74e149a392436b17d523f769e057cb6877f6a5c4618bc6f11800548" + +[[package]] +name = "pollster" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "22686f4785f02a4fcc856d3b3bb19bf6c8160d103f7a99cc258bddd0251dc7f2" + +[[package]] +name = "presser" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e8cf8e6a8aa66ce33f63993ffc4ea4271eb5b0530a9002db8455ea6050c77bfa" + +[[package]] +name = "proc-macro2" +version = "1.0.107" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "985e7ec9bb745e6ce6535b544d84d6cd6f7ad8bd711c398938ae983b91a766d9" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "profiling" +version = "1.0.18" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3d595e54a326bc53c1c197b32d295e14b169e3cfeaa8dc82b529f947fba6bcf5" + +[[package]] +name = "quote" +version = "1.0.47" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fbf4db142a473a8d80c26bbf18454ed458bf8d26c8219c331daecfdbd079001" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "range-alloc" +version = "0.1.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ca45419789ae5a7899559e9512e58ca889e41f04f1f2445e9f4b290ceccd1d08" + +[[package]] +name = "raw-window-handle" +version = "0.6.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "20675572f6f24e9e76ef639bc5552774ed45f1c30e2951e1e99c59888861c539" + +[[package]] +name = "rayon" +version = "1.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fb39b166781f92d482534ef4b4b1b2568f42613b53e5b6c160e24cfbfa30926d" +dependencies = [ + "either", + "rayon-core", +] + +[[package]] +name = "rayon-core" +version = "1.13.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "22e18b0f0062d30d4230b2e85ff77fdfe4326feb054b9783a3460d8435c8ab91" +dependencies = [ + "crossbeam-deque", + "crossbeam-utils", +] + +[[package]] +name = "redox_syscall" +version = "0.5.18" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ed2bf2547551a7053d6fdfafda3f938979645c44812fbfcda098faae3f1a362d" +dependencies = [ + "bitflags 2.13.1", +] + +[[package]] +name = "renderdoc-sys" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "19b30a45b0cd0bcca8037f3d0dc3421eaf95327a17cad11964fb8179b4fc4832" + +[[package]] +name = "rustc-hash" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08d43f7aa6b08d49f382cde6a7982047c3426db949b1424bc4b7ec9ae12c6ce2" + +[[package]] +name = "rustversion" +version = "1.0.23" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cf54715a573b99ac80df0bc206da022bcd442c974952c7b9720069370852e21f" + +[[package]] +name = "scopeguard" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49" + +[[package]] +name = "slab" +version = "0.4.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0c790de23124f9ab44544d7ac05d60440adc586479ce501c1d6d7da3cd8c9cf5" + +[[package]] +name = "slotmap" +version = "1.1.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bdd58c3c93c3d278ca835519292445cb4b0d4dc59ccfdf7ceadaab3f8aeb4038" +dependencies = [ + "version_check", +] + +[[package]] +name = "smallvec" +version = "1.15.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8ed6a63f02c8539c91a8685a86f4099661ba3da017932f6ebbea6de3f0fa7c90" + +[[package]] +name = "spirv" +version = "0.3.0+sdk-1.3.268.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "eda41003dc44290527a59b13432d4a0379379fa074b70174882adfbdfd917844" +dependencies = [ + "bitflags 2.13.1", +] + +[[package]] +name = "static_assertions" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a2eb9349b6444b326872e140eb1cf5e7c522154d69e7a0ffb0fb81c06b37543f" + +[[package]] +name = "strsim" +version = "0.11.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7da8b5736845d9f2fcb837ea5d9e2628564b3b043a70948a3f0b778838c5fb4f" + +[[package]] +name = "syn" +version = "1.0.109" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72b64191b275b66ffe2469e8af2c1cfe3bafa67b529ead792a6d0160888b4237" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "syn" +version = "2.0.119" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "872831b642d1a07999a962a351ed35b955ea2cfc8f3862091e2a240a84f17297" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "syn" +version = "3.0.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "53e9bae58849f64dfa4f5d5ae372c8341f7305f82a3868709269343628b659a3" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "termcolor" +version = "1.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "06794f8f6c5c898b3275aebefa6b8a1cb24cd2c6c79397ab15774837a0bc5755" +dependencies = [ + "winapi-util", +] + +[[package]] +name = "thiserror" +version = "1.0.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b6aaf5339b578ea85b50e080feb250a3e8ae8cfcdff9a461c9ec2904bc923f52" +dependencies = [ + "thiserror-impl", +] + +[[package]] +name = "thiserror-impl" +version = "1.0.69" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4fee6c4efc90059e10f81e6d42c60a18f76588c3d74cb83a0b242a2b6c7504c1" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.119", +] + +[[package]] +name = "unicode-ident" +version = "1.0.24" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6e4313cd5fcd3dad5cafa179702e2b244f760991f45397d14d4ebf38247da75" + +[[package]] +name = "unicode-width" +version = "0.1.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7dd6e30e90baa6f72411720665d41d89b9a3d039dc45b8faea1ddd07f617f6af" + +[[package]] +name = "unicode-xid" +version = "0.2.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ebc1c04c71510c7f702b52b7c350734c9ff1295c464a03335b00bb84fc54f853" + +[[package]] +name = "utf8parse" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "06abde3611657adf66d383f00b093d7faecc7fa57071cce2578660c9f1010821" + +[[package]] +name = "version_check" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" + +[[package]] +name = "wasm-bindgen" +version = "0.2.127" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1b70935747edd64d89de3efa29d73789b806c15798f8e7dca4d8ac356b50ce70" +dependencies = [ + "cfg-if", + "once_cell", + "rustversion", + "wasm-bindgen-macro", + "wasm-bindgen-shared", +] + +[[package]] +name = "wasm-bindgen-futures" +version = "0.4.77" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6b7777d5cc23d0e91404e53ce2d5e8ec7acae3026b16233dba62cd3246457950" +dependencies = [ + "js-sys", + "wasm-bindgen", +] + +[[package]] +name = "wasm-bindgen-macro" +version = "0.2.127" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "77775f8f3f7217702089053b94958f8f54061a3f663417df76e19cbdcca29bc1" +dependencies = [ + "quote", + "wasm-bindgen-macro-support", +] + +[[package]] +name = "wasm-bindgen-macro-support" +version = "0.2.127" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e11d33f857dc2fb11b8bc75aee111aa9cbeb12cd9f25efd3d4c2a3dd4e235284" +dependencies = [ + "bumpalo", + "proc-macro2", + "quote", + "syn 2.0.119", + "wasm-bindgen-shared", +] + +[[package]] +name = "wasm-bindgen-shared" +version = "0.2.127" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7ef64dbcc55df09c7e5a46182d181c2cfa3e925f3da937ea764728b4bbb9dcbf" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "web-sys" +version = "0.3.104" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c435338968042f4f59a557f690a253676d47ce13ceb55d70100e7facf6620a30" +dependencies = [ + "js-sys", + "wasm-bindgen", +] + +[[package]] +name = "weezl" +version = "0.1.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a28ac98ddc8b9274cb41bb4d9d4d5c425b6020c50c46f25559911905610b4a88" + +[[package]] +name = "wgpu" +version = "22.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e1d1c4ba43f80542cf63a0a6ed3134629ae73e8ab51e4b765a67f3aa062eb433" +dependencies = [ + "arrayvec", + "cfg_aliases", + "document-features", + "js-sys", + "log", + "naga", + "parking_lot", + "profiling", + "raw-window-handle", + "smallvec", + "static_assertions", + "wasm-bindgen", + "wasm-bindgen-futures", + "web-sys", + "wgpu-core", + "wgpu-hal", + "wgpu-types", +] + +[[package]] +name = "wgpu-core" +version = "22.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0348c840d1051b8e86c3bcd31206080c5e71e5933dabd79be1ce732b0b2f089a" +dependencies = [ + "arrayvec", + "bit-vec", + "bitflags 2.13.1", + "cfg_aliases", + "document-features", + "indexmap", + "log", + "naga", + "once_cell", + "parking_lot", + "profiling", + "raw-window-handle", + "rustc-hash", + "smallvec", + "thiserror", + "wgpu-hal", + "wgpu-types", +] + +[[package]] +name = "wgpu-hal" +version = "22.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f6bbf4b4de8b2a83c0401d9e5ae0080a2792055f25859a02bf9be97952bbed4f" +dependencies = [ + "android_system_properties", + "arrayvec", + "ash", + "bit-set", + "bitflags 2.13.1", + "block", + "cfg_aliases", + "core-graphics-types", + "d3d12", + "glow", + "glutin_wgl_sys", + "gpu-alloc", + "gpu-allocator", + "gpu-descriptor", + "hassle-rs", + "js-sys", + "khronos-egl", + "libc", + "libloading", + "log", + "metal", + "naga", + "ndk-sys", + "objc", + "once_cell", + "parking_lot", + "profiling", + "range-alloc", + "raw-window-handle", + "renderdoc-sys", + "rustc-hash", + "smallvec", + "thiserror", + "wasm-bindgen", + "web-sys", + "wgpu-types", + "winapi", +] + +[[package]] +name = "wgpu-types" +version = "22.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bc9d91f0e2c4b51434dfa6db77846f2793149d8e73f800fa2e41f52b8eac3c5d" +dependencies = [ + "bitflags 2.13.1", + "js-sys", + "web-sys", +] + +[[package]] +name = "widestring" +version = "1.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72069c3113ab32ab29e5584db3c6ec55d416895e60715417b5b883a357c3e471" + +[[package]] +name = "winapi" +version = "0.3.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5c839a674fcd7a98952e593242ea400abe93992746761e38641405d28b00f419" +dependencies = [ + "winapi-i686-pc-windows-gnu", + "winapi-x86_64-pc-windows-gnu", +] + +[[package]] +name = "winapi-i686-pc-windows-gnu" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ac3b87c63620426dd9b991e5ce0329eff545bccbbb34f3be09ff6fb6ab51b7b6" + +[[package]] +name = "winapi-util" +version = "0.1.11" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c2a7b1c03c876122aa43f3020e6c3c3ee5c05081c9a00739faf7503aeba10d22" +dependencies = [ + "windows-sys", +] + +[[package]] +name = "winapi-x86_64-pc-windows-gnu" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" + +[[package]] +name = "windows" +version = "0.52.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e48a53791691ab099e5e2ad123536d0fff50652600abaf43bbf952894110d0be" +dependencies = [ + "windows-core", + "windows-targets", +] + +[[package]] +name = "windows-core" +version = "0.52.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "33ab640c8d7e35bf8ba19b884ba838ceb4fba93a4e8c65a9059d08afcfc683d9" +dependencies = [ + "windows-targets", +] + +[[package]] +name = "windows-link" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f0805222e57f7521d6a62e36fa9163bc891acd422f971defe97d64e70d0a4fe5" + +[[package]] +name = "windows-sys" +version = "0.61.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ae137229bcbd6cdf0f7b80a31df61766145077ddf49416a728b02cb3921ff3fc" +dependencies = [ + "windows-link", +] + +[[package]] +name = "windows-targets" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9b724f72796e036ab90c1021d4780d4d3d648aca59e491e6b98e725b84e99973" +dependencies = [ + "windows_aarch64_gnullvm", + "windows_aarch64_msvc", + "windows_i686_gnu", + "windows_i686_gnullvm", + "windows_i686_msvc", + "windows_x86_64_gnu", + "windows_x86_64_gnullvm", + "windows_x86_64_msvc", +] + +[[package]] +name = "windows_aarch64_gnullvm" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "32a4622180e7a0ec044bb555404c800bc9fd9ec262ec147edd5989ccd0c02cd3" + +[[package]] +name = "windows_aarch64_msvc" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "09ec2a7bb152e2252b53fa7803150007879548bc709c039df7627cabbd05d469" + +[[package]] +name = "windows_i686_gnu" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8e9b5ad5ab802e97eb8e295ac6720e509ee4c243f69d781394014ebfe8bbfa0b" + +[[package]] +name = "windows_i686_gnullvm" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0eee52d38c090b3caa76c563b86c3a4bd71ef1a819287c19d586d7334ae8ed66" + +[[package]] +name = "windows_i686_msvc" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "240948bc05c5e7c6dabba28bf89d89ffce3e303022809e73deaefe4f6ec56c66" + +[[package]] +name = "windows_x86_64_gnu" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "147a5c80aabfbf0c7d901cb5895d1de30ef2907eb21fbbab29ca94c5b08b1a78" + +[[package]] +name = "windows_x86_64_gnullvm" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "24d5b23dc417412679681396f2b49f3de8c1473deb516bd34410872eff51ed0d" + +[[package]] +name = "windows_x86_64_msvc" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "589f6da84c646204747d1270a2a5661ea66ed1cced2631d546fdfb155959f9ec" + +[[package]] +name = "xml-rs" +version = "0.8.29" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e450f9b2ed1dff33c94c12589a87338689467b9c4f5d8a5710bd09a847d2c8a7" diff --git a/docs/theory/2d_solver/Cargo.toml b/docs/theory/2d_solver/Cargo.toml new file mode 100644 index 0000000..84b3f40 --- /dev/null +++ b/docs/theory/2d_solver/Cargo.toml @@ -0,0 +1,35 @@ +[package] +name = "kbc2d" +version = "0.1.0" +edition = "2021" +rust-version = "1.75" +description = "Двумерный решатель LBM D2Q9 с энтропийным оператором столкновения KBC (модель D / N1)" + +[[bin]] +name = "kbc2d" +path = "src/main.rs" + +[dependencies] +clap = { version = "4", features = ["derive"] } +rayon = "1" +gif = "0.13" +bytemuck = { version = "1", features = ["derive"] } +wgpu = { version = "22", optional = true } +pollster = { version = "0.3", optional = true } + +[features] +default = ["gpu"] +# GPU-бэкенд на wgpu (Vulkan/DX12/Metal). Отключается: --no-default-features +gpu = ["dep:wgpu", "dep:pollster"] + +[profile.release] +opt-level = 3 +lto = "thin" +codegen-units = 1 +panic = "abort" + +# Отдельный профиль для прогонов с проверками переполнения индексов +[profile.release-dbg] +inherits = "release" +debug = true +panic = "unwind" diff --git a/docs/theory/2d_solver/README.md b/docs/theory/2d_solver/README.md new file mode 100644 index 0000000..f02454c --- /dev/null +++ b/docs/theory/2d_solver/README.md @@ -0,0 +1,226 @@ +# kbc2d — двумерный решатель LBM D2Q9 с энтропийным столкновением KBC + +Переписанный на Rust решатель обтекания тела в канале. Физика — та же, что в +`docs/theory/solver_2x_sdf` (python/CuPy), но собранная заново: с тестами против формул +первоисточников, двумя взаимозаменяемыми бэкендами и анимацией, привязанной к физическому +времени потока. + +Модель столкновения — **KBC D** по таблице I работы Bösch, Chikatamarla, Karlin, *Entropic +Multi-Relaxation Models for Simulation of Fluid Turbulence* (arXiv:1507.02509); в трёхмерных +работах тех же авторов она называется **KBC-N1**. Сдвиговая часть `s` несёт натуральные +моменты {N, Π_xy}, всё остальное (T, Q_xyy, Q_yxx, A) уходит в `h`. Статьи лежат в +`docs/origins/`; ссылки на формулы в коде даны по нумерации 2D-статьи. + +## Устройство + +Пять файлов по ролям — каждый отвечает ровно за одно: + +| файл | роль | +|---|---| +| `src/math.rs` | **математика решателя.** Решётка D2Q9, энтропийное равновесие в product-form, проектор на сдвиг, стабилизатор γ, столкновение, Zou–He, геометрия тел через SDF, перевод единиц, спектральная диагностика. Всё узловое и чистое; здесь же тесты против формул статьи. | +| `src/cpu.rs` | **бэкенд под процессор.** Раскладка AoS, обход сетки, rayon, сборка Bouzidi-линков, связка уровней AMR. Физику берёт из `math`. | +| `src/gpu.rs` | **бэкенд под видеокарту.** wgpu + WGSL (Vulkan/DX12/Metal), раскладка SoA. Физика построчно повторяет `math.rs` на f32; топологию (маски, линки, рамка патча) не дублирует, а берёт из `cpu`. | +| `src/main.rs` | **запуск и оркестрирование.** Разбор параметров, сборка постановки, цикл по шагам, живой вывод и итоговый отчёт, выгрузка рядов в CSV. Здесь же контракт `Spec` / `StepRec` / `FieldKind`, общий для обоих бэкендов. | +| `src/gif.rs` | **создание гифок.** Тайминг относительно физического времени, палитры, нормировка, служебная надпись, кодирование. | + +## Сборка и запуск + +Нужен Rust 1.75+. + +```sh +cargo build --release # с GPU-бэкендом +cargo build --release --no-default-features # только CPU (без wgpu) +cargo test --release # 21 быстрый тест +cargo test --release -- --ignored # плюс эталонный сдвиговый слой (~3 с) +``` + +Пример: цилиндр Re=150, гифка завихренности в реальном времени. + +```sh +./target/release/kbc2d --shape cylinder --size 24 --re 150 \ + --nx 480 --ny 240 --steps 40000 --sponge-len 32 \ + --gif wake.gif --gif-field vorticity --verbose full +``` + +`--help` показывает все ключи, разбитые по группам: Физика, Сетка, Тело, Время, Схема, +Анимация, Вывод. + +## Синхронизация анимации с физическим временем + +Требование: гифка идёт с той же скоростью, что и настоящий поток, независимо от того, с какой +скоростью считает машина. Реальная производительность в тайминг не входит вообще. + +В LBM скорость самой решётки жёстко равна c = δx/δt = 1, поэтому шаг по времени однозначно +определяется тем, какую **решёточную** скорость `u_lat` мы назначаем физическому потоку: + +``` +δt = u_lat · δx / u_phys [с/шаг] +шагов в секунду = u_phys / (u_lat · δx) +задержка кадра = (шагов на кадр) · δt / playback +``` + +**Про пример из постановки.** «30 м/с, ячейка 0.1 м ⇒ 300 шагов/с» — это арифметика +δt = δx/u_phys, то есть `u_lat = 1`: поток проходит ровно ячейку за шаг. Формула +воспроизводится буквально ключом `--u-lat 1.0` (тест `units_time_scaling` это проверяет), но +физически такой режим негоден: Ma = u_lat/c_s = √3 ≈ 1.73, сверхзвук, разложение +Чепмена–Энскога не работает. Поэтому по умолчанию `u_lat = 0.05` (Ma ≈ 0.087), и те же 30 м/с +при ячейке 0.1 м дают 6000 шагов/с. Синхронность гифки выдерживается в обоих случаях — меняется +только, сколько шагов приходится на кадр. + +**Две тонкости формата GIF**, обе разобраны: + +1. *Задержка хранится в сотых долях секунды.* Точная физическая задержка почти никогда не целая: + 30 кадр/с — это 3⅓ сотых. Покадровое округление до 3 дало бы анимацию на 11% быстрее + реальности, и уход копился бы линейно (к тысячному кадру — 3.3 секунды). Поэтому задержки + выдаёт накопитель `DelayDither`: суммарное время кадров отслеживает точное физическое с + точностью до одной сотой (3, 3, 4, 3, 3, 4, …), ошибка ограничена ±5 мс и не растёт. +2. *Меньше одной сотой не бывает.* При физически корректном `u_lat` кадр каждые 10 шагов — это + 600 кадр/с, чего формат не умеет. Поэтому режим по умолчанию `--gif-every auto` подбирает + шаг сам под `--gif-fps` (30) так, чтобы получилось ровно реальное время. Если шаг задан + жёстко и задержка не представима, программа печатает фактический коэффициент расхождения и + конкретный совет, как починить. + +`--gif-speed 0.1` даёт замедление в 10 раз (тоже точно, через тот же накопитель). + +## Что параметризовано + +- **поток**: скорость (м/с), направление, число Рейнольдса, решёточная скорость (число Маха); +- **сетка**: размеры домена и размер ячейки в метрах (плотность сетки), коэффициент измельчения + вложенного патча и его границы; +- **тело**: семь форм на выбор — `cylinder`, `square`, `diamond`, `ellipse`, `naca`, `triangle`, + `plate` — плюс характерный размер, относительная толщина, угол атаки и положение; +- **время**: число шагов, длина разгона, амплитуда и длительность стартового возмущения; +- **схема**: оператор столкновения (`kbc`/`bgk`), режим выхода, поглощающая губка, бэкенд, + число потоков; +- **анимация**: файл, поле (`speed`/`vorticity`/`density`/`gamma`), палитра, масштаб, шаг кадра, + частота, скорость воспроизведения, диапазон нормировки; +- **вывод**: период живых строк, три уровня подробности, число окон в отчёте о сходимости, CSV. + +## Что печатает отчёт + +*Шапка* — вся постановка с производными величинами: δt, шагов на секунду, Ma, τ, физическая +вязкость, блокировка канала, геометрия патча, полный план тайминга анимации. + +*Живой вывод* — шаг, физическое время, ⟨ρ⟩, max|u|, Cd, Cl, скорость счёта и ETA; на уровне +`full` дополнительно ⟨γ⟩ с размахом, доля вырожденных узлов, доля узлов с ξ < 0, MLUPS и +отношение скорости счёта к реальному времени. + +*Итог* — установившийся режим (St, ⟨Cd⟩, rms Cl, ⟨Cm⟩) сырой и с поправкой на блокировку, рядом +литературные значения для цилиндра; таблица сходимости по окнам с вердиктом о дрейфе массы и +насыщении; разбор стабилизатора γ; производительность. + +## Состояние проверки + +**Ядро схемы проверено против формул статей** (`cargo test`, 21 тест + 1 длинный): + +- проектор Δs, выписанный аналитически из представления популяций через натуральные моменты + (ур. 10), совпадает с матричным `M⁻¹·diag(…,1,1,…)·M` до 1e-13; он идемпотентен и не несёт + ни массы, ни импульса; +- равновесие в product-form сохраняет ρ и ρu до 1e-13; +- γ из замкнутой оценки (ур. 17) — корень условия критической точки энтропии (ур. 15): невязка + при γ\* более чем в 20 раз меньше, чем при γ\*±1; +- при γ = 2 схема совпадает с LBGK поточечно; +- **сдвиговые моменты релаксируют ровно с 2β при любой γ** (проверено при β = 0.3, 0.6, 0.95) — + это и есть гарантия того, что стабилизатор не трогает вязкость; +- затухание сдвиговой волны даёт ν из ур. (5) с погрешностью < 1% при τ = 0.6 и τ = 1.0; +- Zou–He ставит ровно заданные скорость на входе и плотность на выходе; +- SDF всех семи форм: знак верен внутри и снаружи, |∇φ| = 1 ± 0.05 на контрольном кольце. + +**Сквозная сверка с эталоном.** Дважды периодический сдвиговый слой — один из трёх бенчмарков +2D-статьи (N=128, Re=30000, u₀=0.04, κ=80, δ=0.05, одно конвективное время). Отношение +энстрофии к начальной сходится с fp64-эталоном питоновского решателя **0.6035**. Тест +чувствителен именно к тому, что важно: на испорченном (абсолютном) пороге вырожденности γ тот +же прогон давал 0.6599, то есть +9.3%, а чистый LBGK при этих параметрах разваливается. + +**Паритет бэкендов.** CPU (f64) и GPU (f32) на одной постановке совпадают до 4–5 значащих +цифр шаг в шаг: ⟨ρ⟩ 1.04933 против 1.04934, Cd 2.339 против 2.338, ⟨γ⟩ 1.2645 против 1.2646. +На Intel Iris Xe GPU даёт ≈82 MLUPS против ≈18 MLUPS у процессора. + +**Согласованность уровней AMR.** Один и тот же случай, посчитанный с патчем ×2 и вовсе без +измельчения (`--refine 1`), даёт St 0.1951 против 0.1970 и ⟨Cd⟩ 1.956 против 1.943 — расхождение +в пределах 1–3%. То есть связка уровней (рамка, подшаги, рестрикция) не вносит систематики. + +**Формы тел ведут себя физично.** Один короткий прогон на каждую форму (260×130, размер 20, +угол атаки 12°) даёт ожидаемый порядок сопротивления: обтекаемые — профиль 0.44, эллипс 0.45, +пластина 0.60; промежуточные — цилиндр 1.33, ромб 1.65; тупые — квадрат 2.38, треугольник 2.72. +Момент ⟨Cm⟩ при этом ≈ 0 **только** у круга (−0.0002), которому угол атаки безразличен, а у +несимметричных под углом тел он ненулевой (профиль +0.31, эллипс +0.15, пластина +0.13) — +то есть и плечо, и знак момента считаются осмысленно. + +### Сверка канального случая с питоновским решателем + +Постановка ровно та же, что в `solver_2x_sdf` по умолчанию: 174×90, D=16, cx=40, Re=150, +U=0.07, патч ×2 на [16,140]×[13,77], жёсткий ноль поперечной скорости на выходе, 100 000 шагов. +Слева — что даёт этот решатель, справа — что записано в README питоновского. + +| величина | kbc2d | solver_2x_sdf | | +|---|---|---|---| +| ⟨ρ⟩ | 1.001 | 1.001 | совпало | +| ⟨Cm⟩ | +0.00001 | −0.00003 | оба ≈ 0 — симметрия считывания в порядке | +| rms Cl (с попр.) | 0.571 | 0.557 | расхождение 2.5% | +| ⟨Cd⟩ (сырой) | 1.956 | 1.799 | **на 9% выше** | +| St·(1−β) | 0.161 | 0.183 | **на 12% ниже** | + +Расхождение по Cd и St не объяснено. Что про него известно: + +- оно **не** от связки уровней: без измельчения картина та же (см. выше); +- оно **не** от ядра схемы: сдвиговый слой из статьи воспроизводится с эталоном, вязкость точна + до 1%, сдвиговые моменты релаксируют ровно с 2β; +- после поправки на блокировку **этот** Cd = 1.32 ближе к литературным 1.33 (0.6%), чем + питоновский 1.22 (8.6% ниже), а вот питоновский St ближе к литературным 0.183; +- rms Cl у обоих решателей вдвое выше литературной ~0.3 — это известная особенность постановки, + разобранная в `solver_2x_sdf/README.md` (продольные границы завышают давленческие амплитуды). + +Развести это можно только прогоном обоих кодов бок о бок на одной машине; питоновская версия +требует GPU и CuPy, которых на машине сборки нет, поэтому сравнение сделано с числами, +записанными в её README. + +### Известные отличия от питоновского решателя + +- **Внутри тела столкновение не считается.** Эти популяции фиктивны: Bouzidi перекрывает всё, + что могло бы прийти из тела в жидкость, так что на физику они не влияют. Питоновская версия + считает их наравне со всеми. Побочный эффект — статистика γ здесь собирается строго по + жидкости. +- **Рестрикция дополнительно пропускает узлы, у которых тонкий узел-источник лежит внутри + тела.** В питоновской версии маска строится только по грубому уровню. Случай краевой (тело на + обоих уровнях одно и то же), но здесь он закрыт явно: из фиктивного узла в жидкий переносить + нечего. +- **Точность.** Процессорный бэкенд работает в f64 (питоновский по умолчанию в f32, + переключается переменной `AMR_FP64`); GPU-бэкенд — в f32, потому что в WGSL нет двойной + точности. + +### Порог вырожденности γ + +`GREL = 1e-8` — **относительный** порог, доля от ⟨Δ|Δ⟩, а не абсолютный. Это принципиально: +знаменатель ⟨Δh|Δh⟩ квадратичен по неравновесию и физически мал (~1e-7…1e-9 в развитом следе), +поэтому абсолютный порог срабатывает на подавляющем большинстве узлов и молча подменяет γ на 2, +то есть гонит чистый LBGK вместо KBC. Доля вырожденных узлов печатается в отчёте — на исправном +пороге она обязана быть ~0. + +Единственное исключение — самый первый шаг: поле в точности равно равновесию, Δ ≡ 0, и порог +честно срабатывает на всех узлах. На γ это не влияет, потому что она умножается на Δh = 0. + +## Акустика канала + +Пара «вход по скорости / выход по давлению» — недодемпфированный акустический резонатор: +затухание продольной моды идёт как ν(π/Nx)², то есть на длинном домене её почти ничто не гасит. +Это разобрано в факторном исследовании питоновского решателя (`solver_2x_sdf/README.md`, +эксперимент №1): длинные домены без демпфера дают смещённый режим ⟨ρ⟩ ≈ 1.66 или развал счёта. + +Отсюда два решения в умолчаниях: + +- `--outlet extrapolate` стоит по умолчанию (в том исследовании — лучший вариант по всем + метрикам: жёсткий ноль поперечной скорости отражает вихри дорожки обратно к телу); +- при `nx ≥ 250` и выключенной губке программа печатает предупреждение с готовым рецептом + (`--sponge-len 32`). + +Мгновенное ⟨ρ⟩ при этом всё равно осциллирует вокруг единицы — это сама акустическая мода. +Смотреть надо на **оконные средние** в таблице сходимости: именно они должны стоять на 1.000. + +## Дальше + +- σ·n-кросс-чек силы (интеграл тензора напряжений по контуру) как независимая проверка GMEM — + в питоновской версии есть, здесь пока нет; +- подвижные и вращающиеся тела: GMEM уже записан в галилей-инвариантной форме и принимает + скорость стенки на линке, но подача этой скорости не подключена; +- несколько тел одновременно (сейчас — одно тело на выбор из семи форм). diff --git a/docs/theory/2d_solver/src/cpu.rs b/docs/theory/2d_solver/src/cpu.rs new file mode 100644 index 0000000..c213d8a --- /dev/null +++ b/docs/theory/2d_solver/src/cpu.rs @@ -0,0 +1,935 @@ +//! БЭКЕНД ПОД ПРОЦЕССОР. Раскладка памяти, обход сетки, AMR-связка уровней и параллелизм +//! через rayon. Вся физика берётся из `math` — здесь её нет, есть только организация счёта. +//! +//! Раскладка: AoS, `Vec<[R; Q]>`, узел = y·nx + x. Столкновение — самая дорогая часть шага +//! (sqrt на узел) и в AoS оно тривиально распараллеливается и идеально ложится в кэш; +//! перенос при этом собирает 9 значений с 9 разных узлов, что дешевле, чем кажется, т.к. +//! соседи по x лежат рядом. GPU-бэкенд использует SoA — там важнее коалесцированный доступ. + +use rayon::prelude::*; + +use crate::math::{self, Body, Kbc, Link, LinkKind, R, CX, CY, OPP, Q}; +use crate::{Collision, FieldKind, Spec, StepRec}; + +// ───────────────────────────────────────────────────────────────────────────── +// Геометрия уровня +// ───────────────────────────────────────────────────────────────────────────── + +/// Маски, SDF и предсобранные Bouzidi-линки одного уровня сетки. +pub struct Geom { + pub nx: usize, + pub ny: usize, + pub solid: Vec, + /// Линки «жидкий → твёрдый» с готовой геометрией пересечения. + pub links: Vec, + /// Плечи для момента: координаты узлов относительно центра тела. + pub body_cx: R, + pub body_cy: R, +} + +impl Geom { + /// Собрать геометрию уровня по телу. `solid` определяется знаком SDF: φ ≤ 0 — тело. + pub fn build(nx: usize, ny: usize, body: &Body) -> Self { + let n = nx * ny; + let mut phi = vec![0.0; n]; + let mut solid = vec![false; n]; + for y in 0..ny { + for x in 0..nx { + let p = body.sdf(x as R, y as R); + phi[y * nx + x] = p; + solid[y * nx + x] = p <= 0.0; + } + } + let links = build_links(nx, ny, &solid, &phi); + Geom { nx, ny, solid, links, body_cx: body.cx, body_cy: body.cy } + } +} + +/// Индекс соседа по направлению i с периодическим заворотом (перенос сделан через roll, +/// поэтому решётка периодична по обеим осям; ГУ снимают заворот там, где он не физичен). +#[inline] +fn nb(nx: usize, ny: usize, x: usize, y: usize, i: usize, sign: i32) -> usize { + let xs = (x as i32 + sign * CX[i]).rem_euclid(nx as i32) as usize; + let ys = (y as i32 + sign * CY[i]).rem_euclid(ny as i32) as usize; + ys * nx + xs +} + +/// Для каждого жидкого узла и каждого направления, упирающегося в тело, решаем, какая из трёх +/// формул Bouzidi применима, и считаем долю пересечения q из SDF. Делается один раз. +fn build_links(nx: usize, ny: usize, solid: &[bool], phi: &[R]) -> Vec { + let mut links = Vec::new(); + for y in 0..ny { + for x in 0..nx { + let node = y * nx + x; + if solid[node] { + continue; + } + for i in 1..Q { + let s = nb(nx, ny, x, y, i, 1); // сосед по +c_i + if !solid[s] { + continue; + } + let far = nb(nx, ny, x, y, i, -1); // «дальний» сосед x_f − c_i + let q = math::bouzidi_q(phi[node], phi[s]); + let kind = if q >= 0.5 { + LinkKind::Far + } else if !solid[far] { + LinkKind::Near + } else { + LinkKind::Simple + }; + links.push(Link { + node: node as u32, + far: far as u32, + i: i as u8, + ib: OPP[i] as u8, + kind, + q, + body: true, + }); + } + } + } + links +} + +// ───────────────────────────────────────────────────────────────────────────── +// Уровень сетки +// ───────────────────────────────────────────────────────────────────────────── + +/// Один уровень: поля популяций до/после столкновения, β и геометрия. +pub struct Level { + pub nx: usize, + pub ny: usize, + /// Текущее состояние (после переноса и ГУ). + pub f: Vec<[R; Q]>, + /// Состояние после столкновения — источник переноса и вход Bouzidi/GMEM. + pub post: Vec<[R; Q]>, + /// Энтропийный стабилизатор поузлово (для диагностики и картинки). + pub gamma: Vec, + /// β по столбцам x (губка делает его полем); длина nx. + pub beta: Vec, + pub geom: Geom, +} + +impl Level { + fn new(nx: usize, ny: usize, beta: Vec, geom: Geom) -> Self { + let f0 = math::feq(1.0, 0.0, 0.0); + let n = nx * ny; + Level { nx, ny, f: vec![f0; n], post: vec![f0; n], gamma: vec![2.0; n], beta, geom } + } + + /// Столкновение по всем жидким узлам. Внутри тела не считаем: эти популяции фиктивны + /// (ГУ Bouzidi перезаписывает всё, что могло бы прийти из тела в жидкость), а счёт там + /// только жжёт такты и способен родить NaN при экстремальных режимах. + fn collide(&mut self, op: Collision) -> Stats { + let nx = self.nx; + let beta = &self.beta; + let solid = &self.geom.solid; + let post = &mut self.post; + let gamma = &mut self.gamma; + let f = &self.f; + + post.par_iter_mut() + .zip(gamma.par_iter_mut()) + .enumerate() + .map(|(node, (p, g))| { + *p = f[node]; + if solid[node] { + *g = 2.0; + return Stats::EMPTY; + } + let b = beta[node % nx]; + let k: Kbc = match op { + Collision::Kbc => math::collide_node(p, b), + Collision::Bgk => math::collide_node_bgk(p, b), + }; + *g = k.gamma; + Stats::of(k, b) + }) + .reduce(|| Stats::EMPTY, Stats::merge) + } + + /// Перенос f_i(x, t+1) = f_i(x − c_i, t), периодический по обеим осям. + fn stream(&mut self) { + let (nx, ny) = (self.nx, self.ny); + let post = &self.post; + self.f + .par_chunks_mut(nx) + .enumerate() + .for_each(|(y, row)| { + for (x, cell) in row.iter_mut().enumerate() { + for i in 0..Q { + let xs = (x as i32 - CX[i]).rem_euclid(nx as i32) as usize; + let ys = (y as i32 - CY[i]).rem_euclid(ny as i32) as usize; + cell[i] = post[ys * nx + xs][i]; + } + } + }); + let _ = ny; + } + + /// Интерполированный отскок Bouzidi по всем линкам тела. + /// `f` — поле ПОСЛЕ переноса (его правим), `post` — ПОСЛЕ столкновения (до переноса). + fn apply_bouzidi(&mut self) { + for l in &self.geom.links { + let node = l.node as usize; + let i = l.i as usize; + let ib = l.ib as usize; + let fi = self.post[node][i]; + let v = match l.kind { + LinkKind::Near => 2.0 * l.q * fi + (1.0 - 2.0 * l.q) * self.post[l.far as usize][i], + LinkKind::Far => { + let h = 1.0 / (2.0 * l.q); + h * fi + (1.0 - h) * self.post[node][ib] + } + LinkKind::Simple => fi, + }; + self.f[node][ib] = v; + } + } + + /// FREE-SLIP (зеркальные) стенки канала сверху и снизу: касательный импульс сохраняется, + /// нормальный заворачивается, масса сохраняется. Применять ПОСЛЕ переноса и ПЕРЕД Zou–He. + fn free_slip_walls(&mut self) { + let (nx, ny) = (self.nx, self.ny); + for x in 0..nx { + let b = &mut self.f[x]; + b[2] = b[4]; + b[5] = b[8]; + b[6] = b[7]; + } + let top = (ny - 1) * nx; + for x in 0..nx { + let t = &mut self.f[top + x]; + t[4] = t[2]; + t[7] = t[6]; + t[8] = t[5]; + } + } + + /// Zou–He вход/выход по всему столбцу, включая угловые узлы. + /// + /// ПОРЯДОК КРИТИЧЕН: только ПОСЛЕ free_slip_walls и обязательно на ВЕСЬ столбец. Перенос + /// сделан через roll и потому периодичен по обеим осям: заворот по y снимают стенки, + /// заворот по x — этот ГУ. Если оставить ряды 0 и ny−1 без Zou–He, в них популяции с + /// c_x = +1 на входе приходят прямо из столбца выхода, и вход с выходом оказываются + /// физически связаны в четырёх узлах. + fn channel_bc(&mut self, ux_in: R, uy_in: R, rho_out: R, outlet_extrapolate: bool) { + let (nx, ny) = (self.nx, self.ny); + for y in 0..ny { + math::zou_he_inlet(&mut self.f[y * nx], ux_in, uy_in); + } + for y in 0..ny { + let uy_out = if outlet_extrapolate { + // нуль-градиент поперечной скорости: берём u_y с предвыходного столбца + let c = &self.f[y * nx + nx - 2]; + let s: R = c.iter().sum(); + ((c[2] + c[5] + c[6]) - (c[4] + c[7] + c[8])) / s + } else { + 0.0 + }; + math::zou_he_outlet(&mut self.f[y * nx + nx - 1], rho_out, uy_out); + } + } + + /// Сила и момент на теле по GMEM, суммой по линкам тела. + fn force(&self) -> (R, R, R) { + let (mut fx, mut fy, mut tz) = (0.0, 0.0, 0.0); + let nx = self.nx; + for l in &self.geom.links { + if !l.body { + continue; + } + let node = l.node as usize; + let i = l.i as usize; + let (dfx, dfy) = math::gmem_link(i, self.post[node][i], self.f[node][l.ib as usize], 0.0, 0.0); + fx += dfx; + fy += dfy; + // плечо — до ТОЧКИ ПЕРЕСЕЧЕНИЯ линка со стенкой r_w = r_f + q·c_i, а не до узла: + // узел дал бы ошибку плеча до целой ячейки, а q уже посчитан + let rx = (node % nx) as R - self.geom.body_cx + l.q * CX[i] as R; + let ry = (node / nx) as R - self.geom.body_cy + l.q * CY[i] as R; + tz += rx * dfy - ry * dfx; + } + (fx, fy, tz) + } + + #[inline] + pub fn macros_at(&self, node: usize) -> (R, R, R) { + math::macros(&self.f[node]) + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Статистика KBC за шаг +// ───────────────────────────────────────────────────────────────────────────── + +/// Сводка по стабилизатору γ за один проход столкновения. Собирается редукцией rayon. +#[derive(Clone, Copy, Debug)] +pub struct Stats { + pub n: u64, + pub gsum: R, + pub gmin: R, + pub gmax: R, + /// Узлы, где ⟨Δh|Δh⟩ выродилось и γ подменён на 2 (локально — чистый LBGK). + pub degenerate: u64, + /// Узлы с отрицательной объёмной вязкостью ξ = c_s²(1/(γβ) − ½) — локальное антизатухание. + pub xi_negative: u64, +} + +impl Stats { + pub const EMPTY: Stats = Stats { + n: 0, + gsum: 0.0, + gmin: R::INFINITY, + gmax: R::NEG_INFINITY, + degenerate: 0, + xi_negative: 0, + }; + #[inline] + fn of(k: Kbc, beta: R) -> Stats { + Stats { + n: 1, + gsum: k.gamma, + gmin: k.gamma, + gmax: k.gamma, + degenerate: k.degenerate as u64, + xi_negative: (math::xi_of_gamma(k.gamma, beta) < 0.0) as u64, + } + } + fn merge(a: Stats, b: Stats) -> Stats { + Stats { + n: a.n + b.n, + gsum: a.gsum + b.gsum, + gmin: a.gmin.min(b.gmin), + gmax: a.gmax.max(b.gmax), + degenerate: a.degenerate + b.degenerate, + xi_negative: a.xi_negative + b.xi_negative, + } + } + pub fn gamma_mean(&self) -> R { + if self.n == 0 { + 0.0 + } else { + self.gsum / self.n as R + } + } + pub fn degenerate_frac(&self) -> R { + if self.n == 0 { + 0.0 + } else { + self.degenerate as R / self.n as R + } + } + pub fn xi_negative_frac(&self) -> R { + if self.n == 0 { + 0.0 + } else { + self.xi_negative as R / self.n as R + } + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// AMR-патч +// ───────────────────────────────────────────────────────────────────────────── + +/// Узел рамки тонкого уровня с готовым билинейным стенсилем с грубого уровня. +/// Публичен: GPU-бэкенд переиспользует ту же топологию патча, а не строит её заново. +#[derive(Clone, Copy)] +pub struct Ghost { + pub fine: u32, + pub c00: u32, + pub c10: u32, + pub c01: u32, + pub c11: u32, + pub tx: R, + pub ty: R, +} + +/// Связка «грубый L0 ↔ тонкий L1»: на каждый шаг L0 тонкий уровень делает r подшагов. +/// Рамка патча заполняется интерполяцией с L0 (билинейно по пространству и линейно по +/// времени между состояниями «до» и «после» шага L0), внутренность после подшагов +/// проецируется обратно на L0. +/// +/// Неравновесная часть при смене уровня масштабируется: f^neq ∝ τ·δt, поэтому +/// коэффициент грубый→тонкий равен R01 = τ_f/(r·τ_c), обратно — 1/R01. +pub struct Patch { + pub ax: usize, + pub bx: usize, + pub ay: usize, + pub by: usize, + pub r: usize, + pub nfx: usize, + pub r01: R, + ghosts: Vec, + /// Пары (узел L0, узел L1) для рестрикции — только внутренние жидкие узлы перекрытия. + restrict: Vec<(u32, u32)>, +} + +impl Patch { + /// Топология рамки и рестрикции. Публична: GPU-бэкенд загружает эти же списки в буферы, + /// вместо того чтобы строить их заново. В сборке без GPU читателей у них нет. + #[cfg_attr(not(feature = "gpu"), allow(dead_code))] + pub fn ghosts(&self) -> &[Ghost] { + &self.ghosts + } + #[cfg_attr(not(feature = "gpu"), allow(dead_code))] + pub fn restrict_pairs(&self) -> &[(u32, u32)] { + &self.restrict + } + + pub fn new(spec: &Spec, coarse: &Geom, fine_solid: &[bool], r01: R) -> Patch { + let (ax, bx, ay, by) = spec.patch.expect("патч запрошен без границ"); + let r = spec.refine; + let nfx = r * (bx - ax) + 1; + let nfy = r * (by - ay) + 1; + let cnx = coarse.nx; + + // рамка: один ряд по периметру тонкого поля + let mut ghosts = Vec::with_capacity(2 * (nfx + nfy)); + let push = |gx: usize, gy: usize, out: &mut Vec| { + // координата тонкого узла в системе L0 + let fx = ax as R + gx as R / r as R; + let fy = ay as R + gy as R / r as R; + let x0 = fx.floor() as usize; + let y0 = fy.floor() as usize; + let x1 = (x0 + 1).min(cnx - 1); + let y1 = (y0 + 1).min(coarse.ny - 1); + out.push(Ghost { + fine: (gy * nfx + gx) as u32, + c00: (y0 * cnx + x0) as u32, + c10: (y0 * cnx + x1) as u32, + c01: (y1 * cnx + x0) as u32, + c11: (y1 * cnx + x1) as u32, + tx: fx - x0 as R, + ty: fy - y0 as R, + }); + }; + for gx in 0..nfx { + push(gx, 0, &mut ghosts); + push(gx, nfy - 1, &mut ghosts); + } + for gy in 1..nfy - 1 { + push(0, gy, &mut ghosts); + push(nfx - 1, gy, &mut ghosts); + } + + // рестрикция: каждый r-й тонкий узел внутренности, только там, где на L0 жидкость + let mut restrict = Vec::new(); + for cy in ay + 1..by { + for cx in ax + 1..bx { + let cnode = cy * cnx + cx; + if coarse.solid[cnode] { + continue; + } + let fnode = ((cy - ay) * r) * nfx + (cx - ax) * r; + if fine_solid[fnode] { + continue; + } + restrict.push((cnode as u32, fnode as u32)); + } + } + // рамка — ровно периметр тонкого поля: два ряда по nfx плюс два столбца без углов + assert_eq!(ghosts.len(), 2 * nfx + 2 * (nfy - 2)); + Patch { ax, bx, ay, by, r, nfx, r01, ghosts, restrict } + } + + /// Ghost-значения из грубого поля: равновесие по интерполированным ρ, u плюс + /// масштабированная неравновесная часть. + fn ghost_values(&self, coarse: &[[R; Q]], out: &mut Vec<[R; Q]>) { + out.clear(); + out.reserve(self.ghosts.len()); + for g in &self.ghosts { + let w00 = (1.0 - g.tx) * (1.0 - g.ty); + let w10 = g.tx * (1.0 - g.ty); + let w01 = (1.0 - g.tx) * g.ty; + let w11 = g.tx * g.ty; + // Интерполируются ОТДЕЛЬНО макропеременные и отдельно неравновесная часть: + // рамка = feq(интерполированные ρ, u) + R01·интерполированная neq. Так сделано в + // эталонном python-решателе, и это не то же самое, что интерполировать сами + // популяции: u = interp(ρu)/interp(ρ) отличается от interp(u) во втором порядке, + // и расщепление на eq/neq тогда тоже смещается. Разница мала поузлово, но рамка + // задаёт весь обмен между уровнями каждый шаг, так что копится. + let mut rho = 0.0; + let mut ux = 0.0; + let mut uy = 0.0; + let mut neq = [0.0; Q]; + for (w, c) in [ + (w00, g.c00 as usize), + (w10, g.c10 as usize), + (w01, g.c01 as usize), + (w11, g.c11 as usize), + ] { + let cf = &coarse[c]; + let (r, x, y) = math::macros(cf); + rho += w * r; + ux += w * x; + uy += w * y; + let fe = math::feq(r, x, y); + for i in 0..Q { + neq[i] += w * (cf[i] - fe[i]); + } + } + let fe = math::feq(rho, ux, uy); + let mut v = [0.0; Q]; + for i in 0..Q { + v[i] = fe[i] + self.r01 * neq[i]; + } + out.push(v); + } + } + + /// Записать рамку тонкого поля линейной комбинацией ghost-значений «до» и «после». + fn fill(&self, fine: &mut [[R; Q]], old: &[[R; Q]], new: &[[R; Q]], w: R) { + for (k, g) in self.ghosts.iter().enumerate() { + let dst = &mut fine[g.fine as usize]; + for i in 0..Q { + dst[i] = (1.0 - w) * old[k][i] + w * new[k][i]; + } + } + } + + /// Спроецировать внутренность тонкого уровня обратно на грубый, масштабируя neq на 1/R01. + fn restrict_to(&self, fine: &[[R; Q]], coarse: &mut [[R; Q]]) { + let rfc = 1.0 / self.r01; + for &(cnode, fnode) in &self.restrict { + let ff = &fine[fnode as usize]; + let (rho, ux, uy) = math::macros(ff); + let fe = math::feq(rho, ux, uy); + let dst = &mut coarse[cnode as usize]; + for i in 0..Q { + dst[i] = fe[i] + rfc * (ff[i] - fe[i]); + } + } + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Симуляция +// ───────────────────────────────────────────────────────────────────────────── + +pub struct Sim { + pub spec: Spec, + pub l0: Level, + pub l1: Option, + pub patch: Option, + /// Состояние L0 до столкновения — «старый» край для временной интерполяции рамки. + pre: Vec<[R; Q]>, + gh_old: Vec<[R; Q]>, + gh_new: Vec<[R; Q]>, + fluid_count: R, + probe_node: usize, + probe_on_fine: bool, + pub step_index: u64, +} + +impl Sim { + pub fn new(spec: Spec) -> Sim { + let (nx, ny) = (spec.nx, spec.ny); + let geom0 = Geom::build(nx, ny, &spec.body); + + // губка: плавный рост вязкости в последних sponge_len столбцах. Канал + // «скорость-вход + давление-выход» — недодемпфированный акустический резонатор, + // губка гасит и вихри, и акустику до прихода на выход. + let beta0: Vec = (0..nx) + .map(|x| { + if spec.sponge_len == 0 { + return spec.beta0; + } + let start = nx - 1 - spec.sponge_len; + let s = math::smoothstep((x as R - start as R) / spec.sponge_len as R); + let nu = math::nu_of_beta(spec.beta0); + math::beta_of_nu(nu * (1.0 + (spec.sponge_mult - 1.0) * s)) + }) + .collect(); + + let l0 = Level::new(nx, ny, beta0, geom0); + let fluid_count = l0.geom.solid.iter().filter(|s| !**s).count() as R; + + let (l1, patch) = if spec.refine > 1 { + let (ax, bx, ay, by) = spec.patch.expect("refine > 1 требует патч"); + let r = spec.refine; + let body1 = spec.body.refined(r as R, ax as R, ay as R); + let nfx = r * (bx - ax) + 1; + let nfy = r * (by - ay) + 1; + let geom1 = Geom::build(nfx, nfy, &body1); + // τ_f = r(τ_c − ½) + ½ ⇒ одинаковая ν на обоих уровнях + let tau0 = 1.0 / (2.0 * spec.beta0); + let tau1 = r as R * (tau0 - 0.5) + 0.5; + let beta1 = 1.0 / (2.0 * tau1); + let r01 = tau1 / (r as R * tau0); + let patch = Patch::new(&spec, &l0.geom, &geom1.solid, r01); + (Some(Level::new(nfx, nfy, vec![beta1; nfx], geom1)), Some(patch)) + } else { + (None, None) + }; + + // зонд следа: берём с тонкой сетки, если точка внутри патча (меньше численного + // размытия вихрей → чище спектр и St), иначе с грубой + let (px, py) = spec.probe; + let (probe_node, probe_on_fine) = match &patch { + Some(p) if px >= p.ax && px <= p.bx && py >= p.ay && py <= p.by => { + (((py - p.ay) * p.r) * p.nfx + (px - p.ax) * p.r, true) + } + _ => (py * nx + px, false), + }; + + let n = nx * ny; + Sim { + spec, + l0, + l1, + patch, + pre: vec![[0.0; Q]; n], + gh_old: Vec::new(), + gh_new: Vec::new(), + fluid_count, + probe_node, + probe_on_fine, + step_index: 0, + } + } + + /// Скорость на входе в момент t: разгон smoothstep плюс окно поперечного возмущения. + /// + /// Разгон гасит импульсный старт (мгновенное включение входа шлёт по домену ударную волну). + /// Возмущение — короткий поперечный импульс, сбивающий симметрию: без него дорожка Кармана + /// заводится только на численном шуме и стартует на порядок позже. + fn inlet(&self, t: u64) -> (R, R) { + let sp = &self.spec; + let u = sp.units.u_lat * math::smoothstep(t as R / sp.ramp.max(1) as R); + let (s, c) = sp.flow_angle.sin_cos(); + let mut uy = u * s; + if sp.pert_dur > 0 && t >= sp.ramp && t < sp.ramp + sp.pert_dur { + let ph = (t - sp.ramp) as R / sp.pert_dur as R; + uy += sp.pert_amp * sp.units.u_lat * (std::f64::consts::PI * ph).sin(); + } + (u * c, uy) + } + + pub fn step(&mut self) -> StepRec { + let t = self.step_index; + let (ux_in, uy_in) = self.inlet(t); + let sp_collision = self.spec.collision; + let outlet_extrap = self.spec.outlet_extrapolate; + + // ── уровень 0 ── + self.pre.copy_from_slice(&self.l0.f); + let stats0 = self.l0.collide(sp_collision); + self.l0.stream(); + self.l0.apply_bouzidi(); + self.l0.free_slip_walls(); + self.l0.channel_bc(ux_in, uy_in, 1.0, outlet_extrap); + + // ── уровень 1: r подшагов с временной интерполяцией рамки ── + let (mut fx, mut fy, mut tz) = (0.0, 0.0, 0.0); + if let (Some(l1), Some(p)) = (self.l1.as_mut(), self.patch.as_ref()) { + p.ghost_values(&self.pre, &mut self.gh_old); + p.ghost_values(&self.l0.f, &mut self.gh_new); + for s in 0..p.r { + l1.collide(sp_collision); + l1.stream(); + l1.apply_bouzidi(); + // силу снимаем на КАЖДОМ подшаге и усредняем — мгновенное значение на + // последнем подшаге даёт лишний шум в рядах при том же среднем + let (a, b, c) = l1.force(); + fx += a; + fy += b; + tz += c; + let w = (s + 1) as R / p.r as R; + p.fill(&mut l1.f, &self.gh_old, &self.gh_new, w); + } + let inv = 1.0 / p.r as R; + fx *= inv; + fy *= inv; + tz *= inv; + p.restrict_to(&l1.f, &mut self.l0.f); + } else { + let (a, b, c) = self.l0.force(); + fx = a; + fy = b; + tz = c; + } + + // ── диагностика ── + let probe = if self.probe_on_fine { + self.l1.as_ref().unwrap().macros_at(self.probe_node) + } else { + self.l0.macros_at(self.probe_node) + }; + + let solid = &self.l0.geom.solid; + let (rho_sum, max_u) = self + .l0 + .f + .par_iter() + .enumerate() + .filter(|(n, _)| !solid[*n]) + .map(|(_, c)| { + let (r, ux, uy) = math::macros(c); + (r, (ux * ux + uy * uy).sqrt()) + }) + .reduce(|| (0.0, 0.0), |a, b| (a.0 + b.0, a.1.max(b.1))); + + self.step_index += 1; + StepRec { + step: t, + fx, + fy, + tz, + uy_probe: probe.2, + rho_mean: rho_sum / self.fluid_count, + max_u, + gamma_mean: stats0.gamma_mean(), + gamma_min: stats0.gmin, + gamma_max: stats0.gmax, + degenerate_frac: stats0.degenerate_frac(), + xi_negative_frac: stats0.xi_negative_frac(), + } + } + + /// Характерный размер тела в единицах того уровня, где снимается сила. + pub fn force_ref_size(&self) -> R { + match &self.patch { + Some(p) => self.spec.body.d * p.r as R, + None => self.spec.body.d, + } + } + + /// Поле для картинки: (значение, маска тела) на сетке L0. + pub fn sample_field(&self, kind: FieldKind) -> (Vec, &[bool]) { + let (nx, ny) = (self.l0.nx, self.l0.ny); + let mut out = vec![0.0; nx * ny]; + match kind { + FieldKind::Speed => { + for n in 0..nx * ny { + let (_, ux, uy) = self.l0.macros_at(n); + out[n] = (ux * ux + uy * uy).sqrt(); + } + } + FieldKind::Vorticity => { + // ω = ∂u_y/∂x − ∂u_x/∂y, центральные разности с заворотом по краям + let mut ux = vec![0.0; nx * ny]; + let mut uy = vec![0.0; nx * ny]; + for n in 0..nx * ny { + let (_, a, b) = self.l0.macros_at(n); + ux[n] = a; + uy[n] = b; + } + for y in 0..ny { + for x in 0..nx { + let xp = (x + 1).min(nx - 1); + let xm = x.saturating_sub(1); + let yp = (y + 1).min(ny - 1); + let ym = y.saturating_sub(1); + let dvdx = (uy[y * nx + xp] - uy[y * nx + xm]) / (xp - xm).max(1) as R; + let dudy = (ux[yp * nx + x] - ux[ym * nx + x]) / (yp - ym).max(1) as R; + out[y * nx + x] = dvdx - dudy; + } + } + } + FieldKind::Density => { + for n in 0..nx * ny { + out[n] = self.l0.macros_at(n).0; + } + } + FieldKind::Gamma => out.copy_from_slice(&self.l0.gamma), + } + (out, &self.l0.geom.solid) + } + + /// Есть ли в поле NaN/inf — признак развала счёта. + pub fn is_finite(&self) -> bool { + self.l0.f.par_iter().all(|c| c.iter().all(|v| v.is_finite())) + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Тесты бэкенда +// ───────────────────────────────────────────────────────────────────────────── + +#[cfg(test)] +mod tests { + use super::*; + + fn bare_level(nx: usize, ny: usize) -> Level { + Level::new( + nx, + ny, + vec![0.5; nx], + Geom { + nx, + ny, + solid: vec![false; nx * ny], + links: Vec::new(), + body_cx: 0.0, + body_cy: 0.0, + }, + ) + } + + /// Голый периодический шаг: столкновение + перенос, без единого ГУ. + /// Нужен, чтобы отделить ядро схемы от граничных условий. + fn periodic_step(f: &mut [[R; Q]], tmp: &mut [[R; Q]], nx: usize, ny: usize, beta: R) { + for (n, c) in f.iter().enumerate() { + tmp[n] = *c; + math::collide_node(&mut tmp[n], beta); + } + for y in 0..ny { + for x in 0..nx { + for i in 0..Q { + let xs = (x as i32 - CX[i]).rem_euclid(nx as i32) as usize; + let ys = (y as i32 - CY[i]).rem_euclid(ny as i32) as usize; + f[y * nx + x][i] = tmp[ys * nx + xs][i]; + } + } + } + } + + /// Перенос обязан сдвигать каждую популяцию ровно на её c_i, с заворотом. + #[test] + fn streaming_shifts_by_lattice_velocity() { + let (nx, ny) = (7usize, 5usize); + for i in 0..Q { + let mut lvl = bare_level(nx, ny); + for c in lvl.post.iter_mut() { + *c = [0.0; Q]; + } + lvl.post[2 * nx + 3][i] = 1.0; // дельта в (x=3, y=2) + lvl.stream(); + let xd = (3 + CX[i]).rem_euclid(nx as i32) as usize; + let yd = (2 + CY[i]).rem_euclid(ny as i32) as usize; + assert!( + (lvl.f[yd * nx + xd][i] - 1.0).abs() < 1e-15, + "направление {i}: масса не пришла в ({xd},{yd})" + ); + let total: R = lvl.f.iter().map(|c| c[i]).sum(); + assert!((total - 1.0).abs() < 1e-15, "направление {i}: масса не сохранилась"); + } + } + + /// Однородное равновесие — неподвижная точка схемы: не должно никуда уехать. + #[test] + fn uniform_equilibrium_is_fixed_point() { + let (nx, ny) = (16usize, 16usize); + let mut f = vec![math::feq(1.0, 0.03, -0.01); nx * ny]; + let mut tmp = f.clone(); + let f0 = f[0]; + for _ in 0..50 { + periodic_step(&mut f, &mut tmp, nx, ny, 0.9); + } + for c in &f { + for i in 0..Q { + assert!((c[i] - f0[i]).abs() < 1e-12, "однородное равновесие поехало"); + } + } + } + + /// Затухание сдвиговой волны обязано идти с ν = c_s²(1/(2β) − ½), формула (5). + /// Это прямая проверка того, что стабилизатор γ НЕ трогает вязкость. + #[test] + fn shear_wave_decays_at_prescribed_viscosity() { + let n = 48usize; + for &tau in &[0.6, 1.0] { + let beta = 1.0 / (2.0 * tau); + let nu = math::nu_of_beta(beta); + let k = 2.0 * std::f64::consts::PI / n as R; + let amp = 0.01; + let mut f = vec![[0.0; Q]; n * n]; + for y in 0..n { + for x in 0..n { + f[y * n + x] = math::feq(1.0, amp * (k * y as R).sin(), 0.0); + } + } + let mut tmp = f.clone(); + let mode = |f: &[[R; Q]]| -> R { + let mut s = 0.0; + for y in 0..n { + for x in 0..n { + s += math::macros(&f[y * n + x]).1 * (k * y as R).sin(); + } + } + (2.0 * s / (n * n) as R).abs() + }; + let a0 = mode(&f); + let steps = 1500; + for _ in 0..steps { + periodic_step(&mut f, &mut tmp, n, n, beta); + } + let a1 = mode(&f); + let nu_measured = -(a1 / a0).ln() / (k * k * steps as R); + let err = (nu_measured - nu).abs() / nu; + assert!(err < 1e-2, "τ={tau}: ν измеренная {nu_measured:.6e} vs заданная {nu:.6e}"); + } + } + + /// СКВОЗНАЯ СВЕРКА С ЭТАЛОНОМ. Дважды периодический сдвиговый слой — один из трёх + /// бенчмарков 2D-статьи (разд. VIII). Постановка: N=128, Re=30000, u0=0.04, κ=80, δ=0.05, + /// одно конвективное время t_c = N/u0 шагов. + /// + /// Эталон 0.6035 — отношение энстрофии к начальной на fp64-версии python-решателя + /// (docs/theory/solver_2x_sdf) после исправления порога вырожденности γ. Тот же прогон на + /// испорченном (абсолютном) пороге давал 0.6599, то есть +9.3%, а чистый LBGK при этих + /// параметрах разваливается. Тест ловит и «схема молча выродилась в LBGK», и «схема + /// считает не то». + #[test] + #[ignore = "долгий (3200 шагов на 128²); запуск: cargo test --release -- --ignored"] + fn doubly_periodic_shear_layer_matches_reference() { + let n = 128usize; + let (u0, kappa, delta, re) = (0.04, 80.0, 0.05, 30000.0); + let nu = u0 * n as R / re; + let beta = math::beta_of_nu(nu); + let tc = (n as R / u0) as usize; + + let mut f = vec![[0.0; Q]; n * n]; + for y in 0..n { + let yy = y as R / n as R; + let ux = if y as R <= n as R / 2.0 { + u0 * (kappa * (yy - 0.25)).tanh() + } else { + u0 * (kappa * (0.75 - yy)).tanh() + }; + for x in 0..n { + let uy = delta * u0 * (2.0 * std::f64::consts::PI * (x as R / n as R + 0.25)).sin(); + f[y * n + x] = math::feq(1.0, ux, uy); + } + } + let mut tmp = f.clone(); + + let enstrophy = |f: &[[R; Q]]| -> R { + let mut ux = vec![0.0; n * n]; + let mut uy = vec![0.0; n * n]; + for k in 0..n * n { + let (_, a, b) = math::macros(&f[k]); + ux[k] = a; + uy[k] = b; + } + let mut s = 0.0; + for y in 0..n { + for x in 0..n { + let (xp, xm) = ((x + 1) % n, (x + n - 1) % n); + let (yp, ym) = ((y + 1) % n, (y + n - 1) % n); + let w = 0.5 * (uy[y * n + xp] - uy[y * n + xm]) + - 0.5 * (ux[yp * n + x] - ux[ym * n + x]); + s += w * w; + } + } + s / (n * n) as R + }; + + let e0 = enstrophy(&f); + for _ in 0..=tc { + periodic_step(&mut f, &mut tmp, n, n, beta); + } + assert!(f.iter().all(|c| c.iter().all(|v| v.is_finite())), "счёт развалился"); + let ratio = enstrophy(&f) / e0; + assert!( + (ratio - 0.6035).abs() < 0.012, + "энстрофия/начальная = {ratio:.4}, эталон fp64 python-решателя 0.6035 \ + (испорченный порог γ давал 0.6599)" + ); + } +} + diff --git a/docs/theory/2d_solver/src/gif.rs b/docs/theory/2d_solver/src/gif.rs new file mode 100644 index 0000000..72a373b --- /dev/null +++ b/docs/theory/2d_solver/src/gif.rs @@ -0,0 +1,588 @@ +//! БЛОК СОЗДАНИЯ ГИФОК: перевод поля в индексированный кадр, палитры, служебная надпись и — +//! главное — СИНХРОНИЗАЦИЯ АНИМАЦИИ С ФИЗИЧЕСКИМ ВРЕМЕНЕМ ПОТОКА. +//! +//! Требование: гифка идёт с той же скоростью, что и настоящий поток, независимо от того, с +//! какой скоростью считает машина. Реальная производительность (шагов в секунду wall-clock) +//! в тайминг не входит вообще — берётся только физическая длительность шага δt из `Units`: +//! +//! δt = u_lat·δx/u_phys с/шаг +//! кадр каждые S шагов ⇒ между кадрами S·δt секунд физического времени +//! задержка кадра = S·δt / playback с (playback = 1 ⇒ реальное время) +//! +//! Ограничение формата: задержка в GIF хранится в СОТЫХ долях секунды (u16), т.е. сетка 10 мс, +//! и почти все декодеры считают задержку 0–1 сс как «как можно быстрее». Поэтому реальное время +//! достижимо не при любом S: при δt = 1.7e-4 с (30 м/с, ячейка 0.1 м, u_lat = 0.05) кадр каждые +//! 10 шагов — это 600 кадр/с, чего GIF не умеет. Отсюда два режима: +//! +//! * `--gif-every N` — шаг задан жёстко; если реальное время не представимо, задержка +//! зажимается и в отчёт печатается фактический коэффициент замедления; +//! * `--gif-every auto` — шаг подбирается так, чтобы реальное время получилось точно при +//! целевой частоте `--gif-fps` (по умолчанию 30). Это режим по +//! умолчанию: синхронность важнее, чем конкретное число шагов на кадр. + +use std::fs::File; +use std::io::BufWriter; + +use crate::math::{Units, R}; + +// ───────────────────────────────────────────────────────────────────────────── +// Тайминг +// ───────────────────────────────────────────────────────────────────────────── + +/// Как кадры анимации ложатся на физическое время. +#[derive(Clone, Copy, Debug)] +pub struct GifPlan { + /// Кадр каждые столько шагов симуляции. + pub stride: u64, + /// ТОЧНАЯ задержка кадра в сотых долях секунды — вообще говоря дробная. + pub exact_delay_cs: R, + /// Физическое время между соседними кадрами, с. + pub frame_dt: R, + /// Запрошенная скорость воспроизведения (1 = реальное время, 0.1 = 10× замедление). + pub requested_playback: R, + /// Фактическая, в среднем по анимации. + pub actual_playback: R, + /// Средняя частота кадров получившейся гифки. + pub fps: R, + /// Шаг подбирался автоматически. + pub auto: bool, + /// Точная задержка нецелая ⇒ задержки кадров чередуются (см. `DelayDither`). + pub dithered: bool, + /// Точная задержка меньше одной сотой ⇒ формат физически не тянет, время сжато. + pub clamped: bool, +} + +impl GifPlan { + /// `stride`: Some(N) — жёстко N шагов на кадр; None — подобрать под `target_fps`. + pub fn new(units: &Units, stride: Option, target_fps: R, playback: R) -> GifPlan { + let dt = units.dt; + let playback = if playback > 0.0 { playback } else { 1.0 }; + let auto = stride.is_none(); + + let stride = match stride { + Some(s) => s.max(1), + None => { + // хотим кадр раз в 1/target_fps секунды экранного времени, что при скорости + // playback отвечает playback/target_fps секундам физического времени + ((playback / target_fps / dt).round() as i64).max(1) as u64 + } + }; + + let frame_dt = stride as R * dt; + let exact_delay_cs = 100.0 * frame_dt / playback; + // 0 сотых большинство декодеров трактует как «как можно быстрее», поэтому пол — единица + let clamped = exact_delay_cs < 1.0; + let mean_delay_cs = if clamped { 1.0 } else { exact_delay_cs }; + + GifPlan { + stride, + exact_delay_cs, + frame_dt, + requested_playback: playback, + actual_playback: frame_dt / (mean_delay_cs / 100.0), + fps: 100.0 / mean_delay_cs, + auto, + dithered: !clamped && (exact_delay_cs - exact_delay_cs.round()).abs() > 1e-9, + clamped, + } + } + + /// Отношение фактической скорости воспроизведения к запрошенной (1.0 = точно). + pub fn sync_error(&self) -> R { + self.actual_playback / self.requested_playback + } + + /// Синхронность выдержана: средняя задержка совпадает с точной. + pub fn is_synced(&self) -> bool { + (self.sync_error() - 1.0).abs() < 0.02 + } + + /// Совет, как починить рассинхрон, если формат его не тянет. + pub fn advice(&self, units: &Units) -> Option { + // разойтись с запрошенной скоростью можно только одним способом: точная задержка + // не влезла в минимальную единицу формата и была зажата + if !self.clamped { + return None; + } + let good = ((0.01 * self.requested_playback / units.dt).ceil() as i64).max(1); + Some(format!( + "кадр требует задержки {:.4} сотых — меньше минимальной единицы формата, поэтому \ + анимация идёт в {:.2}× быстрее запрошенного. Починить: --gif-every auto (подберёт \ + шаг сам), либо --gif-every {} (минимальный шаг с представимой задержкой), либо \ + --gif-speed {:.4} (замедление вместо реального времени)", + self.exact_delay_cs, + self.sync_error(), + good, + self.requested_playback / self.sync_error() + )) + } +} + +/// Раскладка дробной задержки по целым сотым без накопления ошибки. +/// +/// Формат GIF хранит задержку кадра целым числом сотых долей секунды. Точная физическая +/// задержка почти никогда не целая: 30 кадр/с — это 3⅓ сотых. Округление каждого кадра +/// по отдельности даёт систематический уход (3 вместо 3.333 — анимация на 11% быстрее +/// реальности, и расхождение копится линейно). Поэтому задержки выдаются так, чтобы +/// НАКОПЛЕННОЕ время кадров всегда совпадало с накопленным физическим с точностью до одной +/// сотой: для 3.333 получается 3, 3, 4, 3, 3, 4, … Ошибка ограничена ±5 мс и не растёт. +#[derive(Clone, Copy, Debug, Default)] +pub struct DelayDither { + exact_cum: R, + emitted_cum: i64, +} + +impl DelayDither { + /// Задержка очередного кадра в сотых долях секунды. + pub fn next(&mut self, exact_delay_cs: R) -> u16 { + self.exact_cum += exact_delay_cs; + let want = self.exact_cum.round() as i64; + let d = (want - self.emitted_cum).max(1); + self.emitted_cum += d; + d.clamp(1, u16::MAX as i64) as u16 + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Палитры +// ───────────────────────────────────────────────────────────────────────────── + +/// Индексы служебных цветов в палитре (после градиента). +const NCOLORS: usize = 250; +const IDX_BODY: u8 = 250; +const IDX_PATCH: u8 = 251; +const IDX_TEXT: u8 = 252; +const IDX_TEXT_BG: u8 = 253; + +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum ColorMap { + /// Последовательная, воспринимаемо-равномерная — для |u| и ρ. + Viridis, + /// Яркая последовательная с широким динамическим диапазоном. + Turbo, + /// Расходящаяся сине-бело-красная — для знакопеременных полей (завихренность, γ−2). + CoolWarm, +} + +impl ColorMap { + pub fn from_str(s: &str) -> Option { + Some(match s.to_ascii_lowercase().as_str() { + "viridis" => ColorMap::Viridis, + "turbo" => ColorMap::Turbo, + "coolwarm" | "rdbu" => ColorMap::CoolWarm, + _ => return None, + }) + } + pub const ALL: [&'static str; 3] = ["viridis", "turbo", "coolwarm"]; + + /// Опорные точки градиента; между ними линейная интерполяция в sRGB. + fn anchors(&self) -> &'static [[u8; 3]] { + match self { + ColorMap::Viridis => &[ + [68, 1, 84], + [72, 40, 120], + [62, 74, 137], + [49, 104, 142], + [38, 130, 142], + [31, 158, 137], + [53, 183, 121], + [109, 205, 89], + [180, 222, 44], + [253, 231, 37], + ], + ColorMap::Turbo => &[ + [48, 18, 59], + [70, 107, 227], + [40, 187, 236], + [49, 242, 153], + [138, 252, 61], + [210, 226, 27], + [254, 165, 45], + [239, 89, 17], + [180, 27, 1], + [122, 4, 3], + ], + // Расходящаяся палитра Морланда. Опорных точек НЕЧЁТНОЕ число: только тогда + // середина диапазона (для знакопеременного поля — ноль) попадает ровно на + // нейтральный серый, а не между двумя соседними точками. + ColorMap::CoolWarm => &[ + [59, 76, 192], + [98, 130, 234], + [141, 176, 254], + [184, 208, 249], + [221, 221, 221], + [245, 196, 173], + [244, 154, 123], + [222, 96, 77], + [180, 4, 38], + ], + } + } + + /// Палитра GIF: NCOLORS ступеней градиента плюс служебные цвета. + fn palette(&self) -> Vec { + let a = self.anchors(); + let mut p = vec![0u8; 256 * 3]; + for k in 0..NCOLORS { + let t = k as R / (NCOLORS - 1) as R * (a.len() - 1) as R; + let i = (t.floor() as usize).min(a.len() - 2); + let fr = t - i as R; + for c in 0..3 { + p[k * 3 + c] = (a[i][c] as R * (1.0 - fr) + a[i + 1][c] as R * fr).round() as u8; + } + } + let set = |p: &mut Vec, idx: u8, rgb: [u8; 3]| { + let o = idx as usize * 3; + p[o] = rgb[0]; + p[o + 1] = rgb[1]; + p[o + 2] = rgb[2]; + }; + set(&mut p, IDX_BODY, [25, 25, 28]); + set(&mut p, IDX_PATCH, [255, 255, 255]); + set(&mut p, IDX_TEXT, [255, 255, 255]); + set(&mut p, IDX_TEXT_BG, [0, 0, 0]); + p + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Микрошрифт 5×7 для служебной надписи +// ───────────────────────────────────────────────────────────────────────────── + +/// Битовая маска глифа: 7 строк по 5 бит (старший бит — левый пиксель). +fn glyph(c: char) -> [u8; 7] { + match c.to_ascii_uppercase() { + '0' => [0x0E, 0x11, 0x13, 0x15, 0x19, 0x11, 0x0E], + '1' => [0x04, 0x0C, 0x04, 0x04, 0x04, 0x04, 0x0E], + '2' => [0x0E, 0x11, 0x01, 0x02, 0x04, 0x08, 0x1F], + '3' => [0x1F, 0x02, 0x04, 0x02, 0x01, 0x11, 0x0E], + '4' => [0x02, 0x06, 0x0A, 0x12, 0x1F, 0x02, 0x02], + '5' => [0x1F, 0x10, 0x1E, 0x01, 0x01, 0x11, 0x0E], + '6' => [0x06, 0x08, 0x10, 0x1E, 0x11, 0x11, 0x0E], + '7' => [0x1F, 0x01, 0x02, 0x04, 0x08, 0x08, 0x08], + '8' => [0x0E, 0x11, 0x11, 0x0E, 0x11, 0x11, 0x0E], + '9' => [0x0E, 0x11, 0x11, 0x0F, 0x01, 0x02, 0x0C], + 'A' => [0x0E, 0x11, 0x11, 0x1F, 0x11, 0x11, 0x11], + 'B' => [0x1E, 0x11, 0x11, 0x1E, 0x11, 0x11, 0x1E], + 'C' => [0x0E, 0x11, 0x10, 0x10, 0x10, 0x11, 0x0E], + 'D' => [0x1C, 0x12, 0x11, 0x11, 0x11, 0x12, 0x1C], + 'E' => [0x1F, 0x10, 0x10, 0x1E, 0x10, 0x10, 0x1F], + 'F' => [0x1F, 0x10, 0x10, 0x1E, 0x10, 0x10, 0x10], + 'G' => [0x0E, 0x11, 0x10, 0x17, 0x11, 0x11, 0x0F], + 'H' => [0x11, 0x11, 0x11, 0x1F, 0x11, 0x11, 0x11], + 'I' => [0x0E, 0x04, 0x04, 0x04, 0x04, 0x04, 0x0E], + 'K' => [0x11, 0x12, 0x14, 0x18, 0x14, 0x12, 0x11], + 'L' => [0x10, 0x10, 0x10, 0x10, 0x10, 0x10, 0x1F], + 'M' => [0x11, 0x1B, 0x15, 0x15, 0x11, 0x11, 0x11], + 'N' => [0x11, 0x19, 0x15, 0x13, 0x11, 0x11, 0x11], + 'O' => [0x0E, 0x11, 0x11, 0x11, 0x11, 0x11, 0x0E], + 'P' => [0x1E, 0x11, 0x11, 0x1E, 0x10, 0x10, 0x10], + 'R' => [0x1E, 0x11, 0x11, 0x1E, 0x14, 0x12, 0x11], + 'S' => [0x0F, 0x10, 0x10, 0x0E, 0x01, 0x01, 0x1E], + 'T' => [0x1F, 0x04, 0x04, 0x04, 0x04, 0x04, 0x04], + 'U' => [0x11, 0x11, 0x11, 0x11, 0x11, 0x11, 0x0E], + 'V' => [0x11, 0x11, 0x11, 0x11, 0x11, 0x0A, 0x04], + 'W' => [0x11, 0x11, 0x11, 0x15, 0x15, 0x1B, 0x11], + 'X' => [0x11, 0x11, 0x0A, 0x04, 0x0A, 0x11, 0x11], + 'Y' => [0x11, 0x11, 0x0A, 0x04, 0x04, 0x04, 0x04], + 'Z' => [0x1F, 0x01, 0x02, 0x04, 0x08, 0x10, 0x1F], + '.' => [0x00, 0x00, 0x00, 0x00, 0x00, 0x0C, 0x0C], + ',' => [0x00, 0x00, 0x00, 0x00, 0x0C, 0x04, 0x08], + ':' => [0x00, 0x0C, 0x0C, 0x00, 0x0C, 0x0C, 0x00], + '-' => [0x00, 0x00, 0x00, 0x1F, 0x00, 0x00, 0x00], + '+' => [0x00, 0x04, 0x04, 0x1F, 0x04, 0x04, 0x00], + '=' => [0x00, 0x00, 0x1F, 0x00, 0x1F, 0x00, 0x00], + '/' => [0x01, 0x02, 0x02, 0x04, 0x08, 0x08, 0x10], + '%' => [0x19, 0x1A, 0x02, 0x04, 0x08, 0x0B, 0x13], + '(' => [0x02, 0x04, 0x08, 0x08, 0x08, 0x04, 0x02], + ')' => [0x08, 0x04, 0x02, 0x02, 0x02, 0x04, 0x08], + 'Ч' => [0x11, 0x11, 0x11, 0x0F, 0x01, 0x01, 0x01], + _ => [0; 7], // пробел и всё незнакомое + } +} + +/// Нарисовать строку в индексированный буфер (координаты — левый верхний угол). +fn draw_text(buf: &mut [u8], w: usize, h: usize, x0: usize, y0: usize, s: &str, scale: usize) { + let gw = 6 * scale; // 5 пикселей глифа + 1 пробел + // подложка, чтобы текст читался на любом фоне + let tw = s.chars().count() * gw; + for yy in y0.saturating_sub(scale)..(y0 + 7 * scale + scale).min(h) { + for xx in x0.saturating_sub(scale)..(x0 + tw + scale).min(w) { + buf[yy * w + xx] = IDX_TEXT_BG; + } + } + for (k, ch) in s.chars().enumerate() { + let g = glyph(ch); + for (row, bits) in g.iter().enumerate() { + for col in 0..5 { + if bits & (1 << (4 - col)) == 0 { + continue; + } + for sy in 0..scale { + for sx in 0..scale { + let px = x0 + k * gw + col * scale + sx; + let py = y0 + row * scale + sy; + if px < w && py < h { + buf[py * w + px] = IDX_TEXT; + } + } + } + } + } + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Диапазон нормировки +// ───────────────────────────────────────────────────────────────────────────── + +/// Диапазон значений, отображаемый на палитру. Фиксируется на весь прогон: плавающая +/// автонормировка делает анимацию нечитаемой — цвет перестаёт что-либо значить. +#[derive(Clone, Copy, Debug)] +pub struct Range { + pub lo: R, + pub hi: R, +} + +impl Range { + #[inline] + fn index(&self, v: R) -> u8 { + if !v.is_finite() { + return 0; + } + let t = ((v - self.lo) / (self.hi - self.lo)).clamp(0.0, 1.0); + (t * (NCOLORS - 1) as R).round() as u8 + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Запись гифки +// ───────────────────────────────────────────────────────────────────────────── + +/// Что подписывать на кадре. +pub struct Hud { + pub u_phys: R, + pub field: &'static str, + pub show: bool, +} + +pub struct GifWriter { + enc: gif::Encoder>, + pub plan: GifPlan, + nx: usize, + ny: usize, + scale: usize, + w: usize, + h: usize, + range: Range, + hud: Hud, + patch: Option<(usize, usize, usize, usize)>, + dither: DelayDither, + pub frames: u32, +} + +impl GifWriter { + #[allow(clippy::too_many_arguments)] + pub fn create( + path: &str, + nx: usize, + ny: usize, + scale: usize, + cmap: ColorMap, + range: Range, + plan: GifPlan, + hud: Hud, + patch: Option<(usize, usize, usize, usize)>, + ) -> std::io::Result { + let scale = scale.max(1); + let w = nx * scale; + let h = ny * scale; + let file = BufWriter::new(File::create(path)?); + let mut enc = gif::Encoder::new(file, w as u16, h as u16, &cmap.palette()) + .map_err(std::io::Error::other)?; + enc.set_repeat(gif::Repeat::Infinite).map_err(std::io::Error::other)?; + Ok(GifWriter { + enc, + plan, + nx, + ny, + scale, + w, + h, + range, + hud, + patch, + dither: DelayDither::default(), + frames: 0, + }) + } + + /// Добавить кадр. `field` — значения по узлам L0, `solid` — маска тела, + /// `t_phys` — физическое время этого кадра (с), `cd` — текущий коэффициент сопротивления. + pub fn push( + &mut self, + field: &[R], + solid: &[bool], + t_phys: R, + cd: Option, + ) -> std::io::Result<()> { + let mut buf = vec![0u8; self.w * self.h]; + // строка 0 изображения — это ВЕРХ, а y = 0 решётки — низ канала: переворачиваем + for y in 0..self.ny { + let src = (self.ny - 1 - y) * self.nx; + for x in 0..self.nx { + let idx = if solid[src + x] { IDX_BODY } else { self.range.index(field[src + x]) }; + for sy in 0..self.scale { + let row = (y * self.scale + sy) * self.w + x * self.scale; + for sx in 0..self.scale { + buf[row + sx] = idx; + } + } + } + } + if let Some((ax, bx, ay, by)) = self.patch { + self.draw_patch_outline(&mut buf, ax, bx, ay, by); + } + if self.hud.show { + let line = match cd { + Some(c) => format!( + "T={:.4}S U={:.1}M/S CD={:.2}", + t_phys, self.hud.u_phys, c + ), + None => format!("T={:.4}S U={:.1}M/S", t_phys, self.hud.u_phys), + }; + let sc = (self.scale.min(3)).max(1); + draw_text(&mut buf, self.w, self.h, 2 * sc, 2 * sc, &line, sc); + draw_text(&mut buf, self.w, self.h, 2 * sc, self.h - 9 * sc, self.hud.field, sc); + } + + let mut frame = gif::Frame::from_indexed_pixels(self.w as u16, self.h as u16, buf, None); + // задержка берётся из накопителя: сумма задержек кадров отслеживает точное + // физическое время, а не округляется покадрово (см. DelayDither) + frame.delay = self.dither.next(self.plan.exact_delay_cs.max(1.0)); + self.enc.write_frame(&frame).map_err(std::io::Error::other)?; + self.frames += 1; + Ok(()) + } + + /// Контур области измельчения — видно, где считает тонкая сетка. + fn draw_patch_outline(&self, buf: &mut [u8], ax: usize, bx: usize, ay: usize, by: usize) { + let flip = |y: usize| self.ny - 1 - y.min(self.ny - 1); + let mut put = |x: usize, y: usize| { + if x < self.nx && y < self.ny { + let py = flip(y) * self.scale; + let px = x * self.scale; + for sy in 0..self.scale { + for sx in 0..self.scale { + let o = (py + sy) * self.w + px + sx; + if o < buf.len() { + buf[o] = IDX_PATCH; + } + } + } + } + }; + for x in ax..=bx.min(self.nx - 1) { + put(x, ay); + put(x, by); + } + for y in ay..=by.min(self.ny - 1) { + put(ax, y); + put(bx, y); + } + } + + pub fn finish(self) -> std::io::Result { + Ok(self.frames) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// Пример из постановки: 30 м/с, ячейка 0.1 м, u_lat = 1 ⇒ 300 шагов/с; кадр каждые + /// 10 шагов ⇒ 30 кадр/с. Точная задержка 3⅓ сотых — нецелая, поэтому включается дизеринг, + /// и в среднем гифка идёт ровно в реальном времени. + #[test] + fn spec_example_is_realtime() { + let u = Units::new(0.1, 30.0, 1.0, 150.0, 16.0); + assert!((u.steps_per_second() - 300.0).abs() < 1e-9); + let p = GifPlan::new(&u, Some(10), 30.0, 1.0); + assert!((p.exact_delay_cs - 10.0 / 3.0).abs() < 1e-9, "{}", p.exact_delay_cs); + assert!(p.dithered && !p.clamped); + assert!((p.fps - 30.0).abs() < 1e-9); + assert!(p.is_synced(), "коэффициент {}", p.sync_error()); + } + + /// Дизеринг обязан вести НАКОПЛЕННОЕ время кадров вплотную к физическому, без ухода. + /// Покадровое округление 3.333 → 3 дало бы к тысячному кадру уход в 3.3 секунды. + #[test] + fn dither_tracks_exact_time_without_drift() { + let exact = 10.0 / 3.0; + let mut d = DelayDither::default(); + let mut sum = 0i64; + for k in 1..=1000 { + sum += d.next(exact) as i64; + let want = exact * k as R; + assert!( + (sum as R - want).abs() <= 0.5 + 1e-9, + "кадр {k}: накоплено {sum} сотых против точных {want:.3}" + ); + } + // а наивное округление к этому моменту отстало бы на треть сотой на кадр + assert!((sum as R - exact * 1000.0).abs() < 1.0); + assert!((3000 - sum).abs() > 300, "дизеринг обязан отличаться от покадрового округления"); + } + + /// Тайминг обязан зависеть ТОЛЬКО от физики, не от того, как быстро считает машина: + /// один и тот же δt даёт одну и ту же задержку, вдвое больший шаг — вдвое большую. + #[test] + fn timing_ignores_wall_clock() { + let u = Units::new(0.1, 30.0, 1.0, 150.0, 16.0); + let a = GifPlan::new(&u, Some(10), 30.0, 1.0); + let b = GifPlan::new(&u, Some(10), 30.0, 1.0); + assert!((a.exact_delay_cs - b.exact_delay_cs).abs() < 1e-12); + let c = GifPlan::new(&u, Some(20), 30.0, 1.0); + assert!((c.exact_delay_cs - 2.0 * a.exact_delay_cs).abs() < 1e-12); + assert!(c.is_synced() && a.is_synced()); + } + + /// Авто-режим обязан вытянуть реальное время при физически корректном u_lat, где жёсткий + /// шаг «каждые 10» потребовал бы 600 кадр/с и формат бы этого не вытянул. + #[test] + fn auto_stride_restores_realtime() { + let u = Units::new(0.1, 30.0, 0.05, 150.0, 16.0); + assert!((u.steps_per_second() - 6000.0).abs() < 1e-6); + let fixed = GifPlan::new(&u, Some(10), 30.0, 1.0); + assert!(fixed.clamped, "0.167 сотых обязаны упереться в пол формата"); + assert!(!fixed.is_synced(), "жёсткий шаг не должен был попасть в реальное время"); + assert!(fixed.advice(&u).is_some(), "в таком режиме обязан быть совет"); + let auto = GifPlan::new(&u, None, 30.0, 1.0); + assert_eq!(auto.stride, 200); // 1/30 с при 6000 шаг/с + assert!(auto.is_synced(), "коэффициент {}", auto.sync_error()); + assert!((auto.fps - 30.0).abs() < 1e-9); + } + + /// Замедление: playback = 0.1 растягивает то же физическое время ровно в 10 раз. + #[test] + fn slow_motion_scales_delay() { + let u = Units::new(0.1, 30.0, 1.0, 150.0, 16.0); + let fast = GifPlan::new(&u, Some(10), 30.0, 1.0); + let slow = GifPlan::new(&u, Some(10), 30.0, 0.1); + assert!((slow.exact_delay_cs - 10.0 * fast.exact_delay_cs).abs() < 1e-9); + assert!((slow.actual_playback - 0.1).abs() < 1e-9); + assert!(slow.is_synced()); + } + + /// Палитра обязана быть ровно 256 цветов по 3 байта и не затирать служебные индексы. + #[test] + fn palette_layout() { + for cm in [ColorMap::Viridis, ColorMap::Turbo, ColorMap::CoolWarm] { + let p = cm.palette(); + assert_eq!(p.len(), 768); + let o = IDX_BODY as usize * 3; + assert_eq!(&p[o..o + 3], &[25, 25, 28]); + } + } +} diff --git a/docs/theory/2d_solver/src/gpu.rs b/docs/theory/2d_solver/src/gpu.rs new file mode 100644 index 0000000..e0d037f --- /dev/null +++ b/docs/theory/2d_solver/src/gpu.rs @@ -0,0 +1,1409 @@ +//! БЭКЕНД ПОД ВИДЕОКАРТУ на wgpu (Vulkan / DX12 / Metal). +//! +//! Физика построчно повторяет `math.rs`, но на f32: в WGSL нет двойной точности. Это законно +//! именно потому, что порог вырожденности γ в схеме ОТНОСИТЕЛЬНЫЙ (`GREL`, доля от ⟨Δ|Δ⟩), а не +//! абсолютный — иначе в f32 он срабатывал бы на подавляющем большинстве узлов и схема молча +//! выродилась бы в LBGK. Топология задачи (маски, Bouzidi-линки, рамка и рестрикция патча) не +//! дублируется: она строится теми же процедурами из `cpu`, а сюда загружается готовыми списками. +//! +//! Раскладка полей — SoA: `f[i*n + node]`. В отличие от AoS на процессоре, здесь важен +//! коалесцированный доступ: соседние потоки читают соседние узлы одной и той же популяции. +//! +//! Ограничение реализации: за каждый шаг делается одно чтение 48-байтового буфера итогов +//! (сила, зонд, статистика). Это синхронизация с GPU на каждом шаге; она заметно ограничивает +//! частоту шагов на мелких сетках, зато `step()` возвращает те же величины, что и процессорный +//! бэкенд, шаг в шаг. На крупных сетках стоимость самой синхронизации теряется на фоне счёта. + +use std::borrow::Cow; + +use bytemuck::{Pod, Zeroable}; +use wgpu::util::DeviceExt; + +use crate::cpu; +use crate::math::{self, R}; +use crate::{Collision, FieldKind, Spec, StepRec}; + +const WG: u32 = 64; + +// ───────────────────────────────────────────────────────────────────────────── +// Структуры, разделяемые с шейдером +// ───────────────────────────────────────────────────────────────────────────── + +#[repr(C)] +#[derive(Clone, Copy, Pod, Zeroable)] +struct LevelParams { + n: u32, + nx: u32, + ny: u32, + nlinks: u32, + bcx: f32, + bcy: f32, + probe_node: u32, + /// бит 0 — писать частичные суммы статистики; бит 1 — на этом уровне лежит зонд + flags: u32, +} + +#[repr(C)] +#[derive(Clone, Copy, Pod, Zeroable)] +struct Dyn { + ux_in: f32, + uy_in: f32, + rho_out: f32, + _pad0: f32, + outlet_extrap: u32, + nparts: u32, + refine: u32, + collision: u32, +} + +#[repr(C)] +#[derive(Clone, Copy, Pod, Zeroable)] +struct Substep { + idx: u32, + _p: [u32; 3], +} + +#[repr(C)] +#[derive(Clone, Copy, Pod, Zeroable)] +struct AmrParams { + nghost: u32, + nrestrict: u32, + ccount: u32, + fcount: u32, + r01: f32, + rfc: f32, + w: f32, + _pad: f32, +} + +#[repr(C)] +#[derive(Clone, Copy, Pod, Zeroable)] +struct GLink { + node: u32, + far: u32, + i: u32, + ib: u32, + kind: u32, + q: f32, + _p: [u32; 2], +} + +#[repr(C)] +#[derive(Clone, Copy, Pod, Zeroable)] +struct GGhost { + fine: u32, + c00: u32, + c10: u32, + c01: u32, + c11: u32, + tx: f32, + ty: f32, + _p: u32, +} + +/// Итоги шага, читаемые на хост. Порядок обязан совпадать с `results` в шейдере. +#[repr(C)] +#[derive(Clone, Copy, Pod, Zeroable, Default, Debug)] +struct Results { + fx: f32, + fy: f32, + tz: f32, + uy_probe: f32, + rho_sum: f32, + max_u: f32, + g_sum: f32, + g_min: f32, + g_max: f32, + cnt: f32, + degen: f32, + xineg: f32, +} + +// ───────────────────────────────────────────────────────────────────────────── +// Шейдер +// ───────────────────────────────────────────────────────────────────────────── + +const SHADER: &str = r#" +// ───── структуры (обязаны совпадать с gpu.rs) ───── +struct LevelParams { n:u32, nx:u32, ny:u32, nlinks:u32, bcx:f32, bcy:f32, probe_node:u32, flags:u32 }; +struct Dyn { ux_in:f32, uy_in:f32, rho_out:f32, p0:f32, outlet_extrap:u32, nparts:u32, refine:u32, collision:u32 }; +struct Substep { idx:u32, p0:u32, p1:u32, p2:u32 }; +struct AmrParams { nghost:u32, nrestrict:u32, ccount:u32, fcount:u32, r01:f32, rfc:f32, w:f32, pad:f32 }; +struct GLink { node:u32, far:u32, i:u32, ib:u32, kind:u32, q:f32, p0:u32, p1:u32 }; +struct GGhost { fine:u32, c00:u32, c10:u32, c01:u32, c11:u32, tx:f32, ty:f32, p0:u32 }; +struct Partial { rho:f32, maxu:f32, gsum:f32, gmin:f32, gmax:f32, cnt:f32, degen:f32, xineg:f32 }; + +@group(0) @binding(0) var f : array; +@group(0) @binding(1) var post : array; +@group(0) @binding(2) var gam : array; +@group(0) @binding(3) var solid: array; +@group(0) @binding(4) var links: array; +@group(0) @binding(5) var beta : array; +@group(0) @binding(6) var parts: array; +@group(0) @binding(7) var P : LevelParams; + +@group(1) @binding(0) var D : Dyn; +@group(1) @binding(1) var S : Substep; +// results[0..11] — итоги шага, results[16 + s*4 + c] — сила подшага s. +// Сила и итоги держатся в ОДНОМ буфере намеренно: связка уровня использует 7 storage-биндингов, +// и отдельный буфер сил вывел бы конвейер за предел 8, гарантируемый лимитами по умолчанию. +@group(1) @binding(2) var results : array; + +@group(2) @binding(0) var amr_src: array; +@group(2) @binding(1) var amr_dst: array; +@group(2) @binding(2) var ghosts : array; +@group(2) @binding(3) var gh_a : array; +@group(2) @binding(4) var gh_b : array; +@group(2) @binding(5) var rest : array>; +@group(2) @binding(6) var A : AmrParams; + +const GREL: f32 = 1e-8; +const CS2 : f32 = 0.3333333333; + +// ───── физика: построчный перенос math.rs ───── + +// Энтропийное равновесие в product-form (точный максимизатор энтропии при заданных rho, rho*u) +fn feq9(rho: f32, ux0: f32, uy0: f32) -> array { + let ux = clamp(ux0, -0.95, 0.95); + let uy = clamp(uy0, -0.95, 0.95); + let sx = sqrt(1.0 + 3.0*ux*ux); + let sy = sqrt(1.0 + 3.0*uy*uy); + let base = rho * (2.0 - sx) * (2.0 - sy); + let qx = (2.0*ux + sx) / (1.0 - ux); + let qy = (2.0*uy + sy) / (1.0 - uy); + let ix = 1.0 / qx; + let iy = 1.0 / qy; + return array( + base*(4.0/9.0), + base*(1.0/9.0)*qx, base*(1.0/9.0)*qy, base*(1.0/9.0)*ix, base*(1.0/9.0)*iy, + base*(1.0/36.0)*qx*qy, base*(1.0/36.0)*ix*qy, base*(1.0/36.0)*ix*iy, base*(1.0/36.0)*qx*iy); +} + +fn macros9(fv: array) -> vec3 { + let rho = fv[0]+fv[1]+fv[2]+fv[3]+fv[4]+fv[5]+fv[6]+fv[7]+fv[8]; + let mx = fv[1]+fv[5]+fv[8]-fv[3]-fv[6]-fv[7]; + let my = fv[2]+fv[5]+fv[6]-fv[4]-fv[7]-fv[8]; + return vec3(rho, mx/rho, my/rho); +} + +// Проекция неравновесия на сдвиговую часть модели KBC D: моменты N = M20-M02 и Pi_xy = M11 +fn shift9(d: array) -> array { + let a = 0.25*(d[1] + d[3] - d[2] - d[4]); + let b = 0.25*(d[5] + d[7] - d[6] - d[8]); + return array(0.0, a, -a, a, -a, b, -b, b, -b); +} + +fn load9(base: ptr>, off: u32, n: u32) { + for (var i = 0u; i < 9u; i = i + 1u) { (*base)[i] = f[i*n + off]; } +} + +// Возвращает (пост-столкновительные популяции, gamma, признак вырождения) +struct CollOut { fv: array, gamma: f32, degen: f32 }; + +fn collide9(fin: array, b: f32, op: u32) -> CollOut { + var out: CollOut; + // WGSL разрешает переменный индекс только по памяти (var), а не по значению (let/параметр), + // поэтому всё, что индексируется в цикле, кладётся в var + var fv = fin; + let m = macros9(fin); + var fe = feq9(m.x, m.y, m.z); + if (op == 1u) { // LBGK: то же самое при gamma = 2 + var g: array; + for (var i = 0u; i < 9u; i = i + 1u) { g[i] = fv[i] + 2.0*b*(fe[i] - fv[i]); } + out.fv = g; out.gamma = 2.0; out.degen = 0.0; + return out; + } + var d: array; + for (var i = 0u; i < 9u; i = i + 1u) { d[i] = fv[i] - fe[i]; } + var ds = shift9(d); + var num = 0.0; var den = 0.0; var nrm = 0.0; + for (var i = 0u; i < 9u; i = i + 1u) { + let inv = 1.0 / fe[i]; + let dh = d[i] - ds[i]; + num = num + ds[i]*dh*inv; + den = den + dh*dh*inv; + nrm = nrm + d[i]*d[i]*inv; + } + // относительный порог: den квадратична по неравновесию и физически мала + let ok = den > GREL*nrm; + let binv = 1.0 / b; + var gm = 2.0; + if (ok) { gm = binv - (2.0 - binv)*num/den; } + var g: array; + for (var i = 0u; i < 9u; i = i + 1u) { + let dh = d[i] - ds[i]; + g[i] = fv[i] - b*(2.0*ds[i] + gm*dh); + } + out.fv = g; out.gamma = gm; + out.degen = select(1.0, 0.0, ok); + return out; +} + +fn zou_he_inlet(fin: array, ux: f32, uy: f32) -> array { + var g = fin; + let rho = (g[0] + g[2] + g[4] + 2.0*(g[3] + g[6] + g[7])) / (1.0 - ux); + let dd = 0.5*(g[2] - g[4]); + g[1] = g[3] + (2.0/3.0)*rho*ux; + g[5] = g[7] - dd + (1.0/6.0)*rho*ux + 0.5*rho*uy; + g[8] = g[6] + dd + (1.0/6.0)*rho*ux - 0.5*rho*uy; + return g; +} + +fn zou_he_outlet(fin: array, rho_out: f32, uy_out: f32) -> array { + var g = fin; + let ux = (g[0] + g[2] + g[4] + 2.0*(g[1] + g[5] + g[8])) / rho_out - 1.0; + let dd = 0.5*(g[2] - g[4]); + g[3] = g[1] - (2.0/3.0)*rho_out*ux; + g[6] = g[8] - dd - (1.0/6.0)*rho_out*ux + 0.5*rho_out*uy_out; + g[7] = g[5] + dd - (1.0/6.0)*rho_out*ux - 0.5*rho_out*uy_out; + return g; +} + +// ───── ядра ───── + +@compute @workgroup_size(64) +fn k_collide(@builtin(global_invocation_id) gid: vec3) { + let nd = gid.x; + if (nd >= P.n) { return; } + var fv: array; + load9(&fv, nd, P.n); + if (solid[nd] != 0u) { + // внутри тела не считаем: популяции там фиктивны, Bouzidi всё равно их перекрывает + for (var i = 0u; i < 9u; i = i + 1u) { post[i*P.n + nd] = fv[i]; } + gam[nd] = 2.0; + return; + } + var r = collide9(fv, beta[nd % P.nx], D.collision); + for (var i = 0u; i < 9u; i = i + 1u) { post[i*P.n + nd] = r.fv[i]; } + gam[nd] = r.gamma; +} + +@compute @workgroup_size(64) +fn k_stream(@builtin(global_invocation_id) gid: vec3) { + let nd = gid.x; + if (nd >= P.n) { return; } + var cx = array(0, 1, 0, -1, 0, 1, -1, -1, 1); + var cy = array(0, 0, 1, 0, -1, 1, 1, -1, -1); + let x = i32(nd % P.nx); + let y = i32(nd / P.nx); + let nxi = i32(P.nx); + let nyi = i32(P.ny); + for (var i = 0u; i < 9u; i = i + 1u) { + var xs = (x - cx[i]) % nxi; if (xs < 0) { xs = xs + nxi; } + var ys = (y - cy[i]) % nyi; if (ys < 0) { ys = ys + nyi; } + f[i*P.n + nd] = post[i*P.n + u32(ys*nxi + xs)]; + } +} + +@compute @workgroup_size(64) +fn k_bouzidi(@builtin(global_invocation_id) gid: vec3) { + let k = gid.x; + if (k >= P.nlinks) { return; } + let L = links[k]; + let fi = post[L.i*P.n + L.node]; + var v = fi; + if (L.kind == 0u) { // q < 1/2, есть дальний жидкий сосед + v = 2.0*L.q*fi + (1.0 - 2.0*L.q)*post[L.i*P.n + L.far]; + } else if (L.kind == 1u) { // q >= 1/2 + let h = 1.0/(2.0*L.q); + v = h*fi + (1.0 - h)*post[L.ib*P.n + L.node]; + } + f[L.ib*P.n + L.node] = v; +} + +@compute @workgroup_size(64) +fn k_walls(@builtin(global_invocation_id) gid: vec3) { + let x = gid.x; + if (x >= P.nx) { return; } + // зеркальное отражение: касательный импульс сохраняется, нормальный заворачивается + f[2u*P.n + x] = f[4u*P.n + x]; + f[5u*P.n + x] = f[8u*P.n + x]; + f[6u*P.n + x] = f[7u*P.n + x]; + let t = (P.ny - 1u)*P.nx + x; + f[4u*P.n + t] = f[2u*P.n + t]; + f[7u*P.n + t] = f[6u*P.n + t]; + f[8u*P.n + t] = f[5u*P.n + t]; +} + +@compute @workgroup_size(64) +fn k_channel(@builtin(global_invocation_id) gid: vec3) { + let y = gid.x; + if (y >= P.ny) { return; } + // вход: скоростной Zou-He по ВСЕМУ столбцу, включая угловые узлы + let a = y*P.nx; + var fv: array; + load9(&fv, a, P.n); + var gi = zou_he_inlet(fv, D.ux_in, D.uy_in); + for (var i = 0u; i < 9u; i = i + 1u) { f[i*P.n + a] = gi[i]; } + // выход: давление-Zou-He + let b = y*P.nx + P.nx - 1u; + var gv: array; + load9(&gv, b, P.n); + var uyo = 0.0; + if (D.outlet_extrap != 0u) { + let c = y*P.nx + P.nx - 2u; + var cv: array; + load9(&cv, c, P.n); + let s = cv[0]+cv[1]+cv[2]+cv[3]+cv[4]+cv[5]+cv[6]+cv[7]+cv[8]; + uyo = ((cv[2]+cv[5]+cv[6]) - (cv[4]+cv[7]+cv[8])) / s; + } + var go = zou_he_outlet(gv, D.rho_out, uyo); + for (var i = 0u; i < 9u; i = i + 1u) { f[i*P.n + b] = go[i]; } +} + +// сила и момент по GMEM, редукция одной рабочей группой +var wfx: array; +var wfy: array; +var wtz: array; + +@compute @workgroup_size(256) +fn k_force(@builtin(local_invocation_id) lid: vec3) { + let t = lid.x; + var cx = array(0.0, 1.0, 0.0, -1.0, 0.0, 1.0, -1.0, -1.0, 1.0); + var cy = array(0.0, 0.0, 1.0, 0.0, -1.0, 1.0, 1.0, -1.0, -1.0); + var sx = 0.0; var sy = 0.0; var sz = 0.0; + var k = t; + loop { + if (k >= P.nlinks) { break; } + let L = links[k]; + let fp = post[L.i*P.n + L.node]; + let fb = f[L.ib*P.n + L.node]; + let dfx = cx[L.i]*fp + cx[L.i]*fb; + let dfy = cy[L.i]*fp + cy[L.i]*fb; + sx = sx + dfx; + sy = sy + dfy; + // плечо до точки пересечения линка со стенкой, а не до узла + let rx = f32(L.node % P.nx) - P.bcx + L.q*cx[L.i]; + let ry = f32(L.node / P.nx) - P.bcy + L.q*cy[L.i]; + sz = sz + rx*dfy - ry*dfx; + k = k + 256u; + } + wfx[t] = sx; wfy[t] = sy; wtz[t] = sz; + workgroupBarrier(); + var s = 128u; + loop { + if (s == 0u) { break; } + if (t < s) { + wfx[t] = wfx[t] + wfx[t + s]; + wfy[t] = wfy[t] + wfy[t + s]; + wtz[t] = wtz[t] + wtz[t + s]; + } + workgroupBarrier(); + s = s >> 1u; + } + if (t == 0u) { + results[16u + S.idx*4u + 0u] = wfx[0]; + results[16u + S.idx*4u + 1u] = wfy[0]; + results[16u + S.idx*4u + 2u] = wtz[0]; + } +} + +// частичные суммы диагностики: одна Partial на рабочую группу +var wp: array; + +@compute @workgroup_size(64) +fn k_stats1(@builtin(global_invocation_id) gid: vec3, + @builtin(local_invocation_id) lid: vec3, + @builtin(workgroup_id) wid: vec3) { + let nd = gid.x; + let t = lid.x; + var p: Partial; + p.rho = 0.0; p.maxu = 0.0; p.gsum = 0.0; p.gmin = 1e30; p.gmax = -1e30; + p.cnt = 0.0; p.degen = 0.0; p.xineg = 0.0; + if (nd < P.n && solid[nd] == 0u) { + var fv: array; + load9(&fv, nd, P.n); + let m = macros9(fv); + p.rho = m.x; + p.maxu = sqrt(m.y*m.y + m.z*m.z); + let g = gam[nd]; + p.gsum = g; p.gmin = g; p.gmax = g; p.cnt = 1.0; + // вырожденный узел помечен ровно gamma = 2 (см. k_collide) + p.degen = select(0.0, 1.0, g == 2.0); + // объёмная вязкость модели D: xi = cs^2 (1/(gamma*beta) - 1/2) + let b = beta[nd % P.nx]; + let xi = CS2*(1.0/(g*b) - 0.5); + p.xineg = select(0.0, 1.0, xi < 0.0); + } + // зонд следа снимается с того уровня, на котором он лежит + if ((P.flags & 2u) != 0u && nd == P.probe_node) { + var fv: array; + load9(&fv, nd, P.n); + let m = macros9(fv); + results[3] = m.z; + } + wp[t] = p; + workgroupBarrier(); + var s = 32u; + loop { + if (s == 0u) { break; } + if (t < s) { + let o = wp[t + s]; + wp[t].rho = wp[t].rho + o.rho; + wp[t].maxu = max(wp[t].maxu, o.maxu); + wp[t].gsum = wp[t].gsum + o.gsum; + wp[t].gmin = min(wp[t].gmin, o.gmin); + wp[t].gmax = max(wp[t].gmax, o.gmax); + wp[t].cnt = wp[t].cnt + o.cnt; + wp[t].degen = wp[t].degen + o.degen; + wp[t].xineg = wp[t].xineg + o.xineg; + } + workgroupBarrier(); + s = s >> 1u; + } + if (t == 0u && (P.flags & 1u) != 0u) { parts[wid.x] = wp[0]; } +} + +// Размер группы обязан совпадать с длиной wp: при 256 потоках на 64 ячейки четыре потока +// писали бы в одну и ту же ячейку и три четверти частичных сумм терялись бы молча. +@compute @workgroup_size(64) +fn k_stats2(@builtin(local_invocation_id) lid: vec3) { + let t = lid.x; + var p: Partial; + p.rho = 0.0; p.maxu = 0.0; p.gsum = 0.0; p.gmin = 1e30; p.gmax = -1e30; + p.cnt = 0.0; p.degen = 0.0; p.xineg = 0.0; + var k = t; + loop { + if (k >= D.nparts) { break; } + let o = parts[k]; + p.rho = p.rho + o.rho; + p.maxu = max(p.maxu, o.maxu); + p.gsum = p.gsum + o.gsum; + p.gmin = min(p.gmin, o.gmin); + p.gmax = max(p.gmax, o.gmax); + p.cnt = p.cnt + o.cnt; + p.degen = p.degen + o.degen; + p.xineg = p.xineg + o.xineg; + k = k + 64u; + } + wp[t] = p; + workgroupBarrier(); + // сведение 256 потоков через 64 ячейки делаем последовательно нулевым потоком: + // объём данных крошечный, а корректность важнее пары микросекунд + if (t == 0u) { + var q: Partial; + q.rho = 0.0; q.maxu = 0.0; q.gsum = 0.0; q.gmin = 1e30; q.gmax = -1e30; + q.cnt = 0.0; q.degen = 0.0; q.xineg = 0.0; + for (var j = 0u; j < 64u; j = j + 1u) { + let o = wp[j]; + q.rho = q.rho + o.rho; + q.maxu = max(q.maxu, o.maxu); + q.gsum = q.gsum + o.gsum; + q.gmin = min(q.gmin, o.gmin); + q.gmax = max(q.gmax, o.gmax); + q.cnt = q.cnt + o.cnt; + q.degen = q.degen + o.degen; + q.xineg = q.xineg + o.xineg; + } + var fx = 0.0; var fy = 0.0; var tz = 0.0; + for (var s = 0u; s < D.refine; s = s + 1u) { + fx = fx + results[16u + s*4u + 0u]; + fy = fy + results[16u + s*4u + 1u]; + tz = tz + results[16u + s*4u + 2u]; + } + let inv = 1.0 / f32(D.refine); + results[0] = fx*inv; + results[1] = fy*inv; + results[2] = tz*inv; + results[4] = q.rho; + results[5] = q.maxu; + results[6] = q.gsum; + results[7] = q.gmin; + results[8] = q.gmax; + results[9] = q.cnt; + results[10] = q.degen; + results[11] = q.xineg; + } +} + +// ───── AMR ───── + +@compute @workgroup_size(64) +fn k_ghost(@builtin(global_invocation_id) gid: vec3) { + let k = gid.x; + if (k >= A.nghost) { return; } + let g = ghosts[k]; + let w00 = (1.0 - g.tx)*(1.0 - g.ty); + let w10 = g.tx*(1.0 - g.ty); + let w01 = (1.0 - g.tx)*g.ty; + let w11 = g.tx*g.ty; + // отдельно интерполируются макропеременные, отдельно неравновесная часть (как в cpu.rs) + var ws = array(w00, w10, w01, w11); + var cs = array(g.c00, g.c10, g.c01, g.c11); + var rho = 0.0; var ux = 0.0; var uy = 0.0; + var neq = array(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0); + for (var j = 0u; j < 4u; j = j + 1u) { + var cf: array; + for (var i = 0u; i < 9u; i = i + 1u) { cf[i] = amr_src[i*A.ccount + cs[j]]; } + let m = macros9(cf); + rho = rho + ws[j]*m.x; + ux = ux + ws[j]*m.y; + uy = uy + ws[j]*m.z; + var fe0 = feq9(m.x, m.y, m.z); + for (var i = 0u; i < 9u; i = i + 1u) { neq[i] = neq[i] + ws[j]*(cf[i] - fe0[i]); } + } + var fe = feq9(rho, ux, uy); + // неравновесная часть масштабируется: f^neq ~ tau*dt + for (var i = 0u; i < 9u; i = i + 1u) { amr_dst[k*9u + i] = fe[i] + A.r01*neq[i]; } +} + +@compute @workgroup_size(64) +fn k_fill(@builtin(global_invocation_id) gid: vec3) { + let k = gid.x; + if (k >= A.nghost) { return; } + let g = ghosts[k]; + // временная интерполяция рамки между состояниями L0 «до» и «после» шага + for (var i = 0u; i < 9u; i = i + 1u) { + amr_dst[i*A.fcount + g.fine] = (1.0 - A.w)*gh_a[k*9u + i] + A.w*gh_b[k*9u + i]; + } +} + +@compute @workgroup_size(64) +fn k_restrict(@builtin(global_invocation_id) gid: vec3) { + let k = gid.x; + if (k >= A.nrestrict) { return; } + let pr = rest[k]; + var ff: array; + for (var i = 0u; i < 9u; i = i + 1u) { ff[i] = amr_src[i*A.fcount + pr.y]; } + let m = macros9(ff); + var fe = feq9(m.x, m.y, m.z); + for (var i = 0u; i < 9u; i = i + 1u) { + amr_dst[i*A.ccount + pr.x] = fe[i] + A.rfc*(ff[i] - fe[i]); + } +} +"#; + +// ───────────────────────────────────────────────────────────────────────────── +// Уровень на устройстве +// ───────────────────────────────────────────────────────────────────────────── + +struct GpuLevel { + nx: usize, + ny: usize, + n: usize, + f: wgpu::Buffer, + /// Стабилизатор γ поузлово — нужен для картинки; забирается с устройства как есть. + gam: wgpu::Buffer, + bind: wgpu::BindGroup, + /// число рабочих групп редукции статистики (оно же длина буфера частичных сумм) + parts_count: u32, + nlinks: u32, +} + +impl GpuLevel { + fn link_groups(&self) -> u32 { + ceil_div(self.nlinks.max(1), WG) + } +} + +fn soa_equilibrium(n: usize) -> Vec { + let fe = math::feq(1.0, 0.0, 0.0); + let mut v = vec![0.0f32; 9 * n]; + for i in 0..9 { + for k in 0..n { + v[i * n + k] = fe[i] as f32; + } + } + v +} + +/// Залить уровень на устройство: поля, маски, линки, β и групповая привязка. +fn make_level( + device: &wgpu::Device, + layout: &wgpu::BindGroupLayout, + nx: usize, + ny: usize, + geom: &cpu::Geom, + beta: &[f32], + probe_node: u32, + flags: u32, +) -> GpuLevel { + let n = nx * ny; + let init = soa_equilibrium(n); + let mkf = |label: &str| { + device.create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: Some(label), + contents: bytemuck::cast_slice(&init), + usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC + | wgpu::BufferUsages::COPY_DST, + }) + }; + let f = mkf("f"); + let post = mkf("post"); + let gam = device.create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: Some("gamma"), + contents: bytemuck::cast_slice(&vec![2.0f32; n]), + usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC, + }); + let solid: Vec = geom.solid.iter().map(|s| *s as u32).collect(); + let solid_buf = storage_init(device, "solid", bytemuck::cast_slice(&solid)); + let links: Vec = geom + .links + .iter() + .map(|l| GLink { + node: l.node, + far: l.far, + i: l.i as u32, + ib: l.ib as u32, + kind: match l.kind { + math::LinkKind::Near => 0, + math::LinkKind::Far => 1, + math::LinkKind::Simple => 2, + }, + q: l.q as f32, + _p: [0; 2], + }) + .collect(); + let links_buf = storage_init(device, "links", bytemuck::cast_slice(&links)); + let beta_buf = storage_init(device, "beta", bytemuck::cast_slice(beta)); + + let parts_count = ceil_div(n as u32, WG); + let parts = storage(device, "parts", (parts_count as u64) * 32); + let params = device.create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: Some("level params"), + contents: bytemuck::bytes_of(&LevelParams { + n: n as u32, + nx: nx as u32, + ny: ny as u32, + nlinks: links.len() as u32, + bcx: geom.body_cx as f32, + bcy: geom.body_cy as f32, + probe_node, + flags, + }), + usage: wgpu::BufferUsages::UNIFORM, + }); + + let bind = device.create_bind_group(&wgpu::BindGroupDescriptor { + label: Some("level"), + layout, + entries: &[ + bind(0, &f), + bind(1, &post), + bind(2, &gam), + bind(3, &solid_buf), + bind(4, &links_buf), + bind(5, &beta_buf), + bind(6, &parts), + bind(7, ¶ms), + ], + }); + GpuLevel { nx, ny, n, f, gam, bind, parts_count, nlinks: links.len() as u32 } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Симуляция +// ───────────────────────────────────────────────────────────────────────────── + +pub struct Sim { + spec: Spec, + device: wgpu::Device, + queue: wgpu::Queue, + adapter_name: String, + + l0: GpuLevel, + l1: Option, + pre: wgpu::Buffer, + + dyn_buf: wgpu::Buffer, + results_buf: wgpu::Buffer, + staging: wgpu::Buffer, + /// по одной группе на подшаг — различаются только индексом подшага в `Substep` + dyn_binds: Vec, + + amr: Option, + + pipes: Pipes, + readback: wgpu::Buffer, + + fluid_count: R, + step_index: u64, + d_ref: R, +} + +struct AmrRes { + nghost: u32, + nrestrict: u32, + bg_ghost_old: wgpu::BindGroup, + bg_ghost_new: wgpu::BindGroup, + bg_fill: Vec, + bg_restrict: wgpu::BindGroup, +} + +struct Pipes { + collide: wgpu::ComputePipeline, + stream: wgpu::ComputePipeline, + bouzidi: wgpu::ComputePipeline, + walls: wgpu::ComputePipeline, + channel: wgpu::ComputePipeline, + force: wgpu::ComputePipeline, + stats1: wgpu::ComputePipeline, + stats2: wgpu::ComputePipeline, + ghost: wgpu::ComputePipeline, + fill: wgpu::ComputePipeline, + restrict: wgpu::ComputePipeline, + empty_bg: wgpu::BindGroup, +} + +impl Sim { + pub fn new(spec: Spec) -> Result { + let instance = wgpu::Instance::new(wgpu::InstanceDescriptor { + backends: wgpu::Backends::PRIMARY, + ..Default::default() + }); + let adapter = pollster::block_on(instance.request_adapter(&wgpu::RequestAdapterOptions { + power_preference: wgpu::PowerPreference::HighPerformance, + compatible_surface: None, + force_fallback_adapter: false, + })) + .ok_or("подходящий GPU-адаптер не найден (нужен Vulkan, DX12 или Metal)")?; + let info = adapter.get_info(); + let adapter_name = format!("{} ({:?})", info.name, info.backend); + + let (device, queue) = pollster::block_on(adapter.request_device( + &wgpu::DeviceDescriptor { + label: Some("kbc2d"), + required_features: wgpu::Features::empty(), + // берём то, что реально умеет адаптер: связке уровня нужно 8 storage-биндингов, + // downlevel-профиль даёт всего 4 + required_limits: adapter.limits(), + memory_hints: wgpu::MemoryHints::Performance, + }, + None, + )) + .map_err(|e| format!("не удалось получить устройство: {e}"))?; + + // ── топология берётся из процессорного бэкенда, а не строится заново ── + let geom0 = cpu::Geom::build(spec.nx, spec.ny, &spec.body); + let fluid_count = geom0.solid.iter().filter(|s| !**s).count() as R; + + let beta0: Vec = (0..spec.nx) + .map(|x| { + if spec.sponge_len == 0 { + return spec.beta0 as f32; + } + let start = spec.nx - 1 - spec.sponge_len; + let s = math::smoothstep((x as R - start as R) / spec.sponge_len as R); + let nu = math::nu_of_beta(spec.beta0); + math::beta_of_nu(nu * (1.0 + (spec.sponge_mult - 1.0) * s)) as f32 + }) + .collect(); + + let module = device.create_shader_module(wgpu::ShaderModuleDescriptor { + label: Some("kbc2d.wgsl"), + source: wgpu::ShaderSource::Wgsl(Cow::Borrowed(SHADER)), + }); + + let bgl_level = level_layout(&device); + let bgl_dyn = dyn_layout(&device); + let bgl_amr = amr_layout(&device); + let bgl_empty = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor { + label: Some("empty"), + entries: &[], + }); + + let pl_level = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor { + label: Some("level"), + bind_group_layouts: &[&bgl_level, &bgl_dyn], + push_constant_ranges: &[], + }); + let pl_amr = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor { + label: Some("amr"), + bind_group_layouts: &[&bgl_empty, &bgl_empty, &bgl_amr], + push_constant_ranges: &[], + }); + + let mk = |layout: &wgpu::PipelineLayout, entry: &str| { + device.create_compute_pipeline(&wgpu::ComputePipelineDescriptor { + label: Some(entry), + layout: Some(layout), + module: &module, + entry_point: entry, + compilation_options: Default::default(), + cache: None, + }) + }; + + // ── общие буферы ── + let dyn_buf = device.create_buffer(&wgpu::BufferDescriptor { + label: Some("dyn"), + size: std::mem::size_of::() as u64, + usage: wgpu::BufferUsages::UNIFORM | wgpu::BufferUsages::COPY_DST, + mapped_at_creation: false, + }); + // 16 слотов итогов + 8 подшагов по 4 числа на силу + let results_buf = device.create_buffer(&wgpu::BufferDescriptor { + label: Some("results"), + size: (16 + 8 * 4) * 4, + usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC, + mapped_at_creation: false, + }); + let staging = device.create_buffer(&wgpu::BufferDescriptor { + label: Some("staging"), + size: std::mem::size_of::() as u64, + usage: wgpu::BufferUsages::MAP_READ | wgpu::BufferUsages::COPY_DST, + mapped_at_creation: false, + }); + + let nsub = spec.refine.max(1); + let dyn_binds: Vec = (0..nsub) + .map(|s| { + let sb = device.create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: Some("substep"), + contents: bytemuck::bytes_of(&Substep { idx: s as u32, _p: [0; 3] }), + usage: wgpu::BufferUsages::UNIFORM, + }); + device.create_bind_group(&wgpu::BindGroupDescriptor { + label: Some("dyn"), + layout: &bgl_dyn, + entries: &[bind(0, &dyn_buf), bind(1, &sb), bind(2, &results_buf)], + }) + }) + .collect(); + + // ── зонд ── + let (px, py) = spec.probe; + let probe_on_fine = matches!(spec.patch, Some((ax, bx, ay, by)) + if px >= ax && px <= bx && py >= ay && py <= by) + && spec.refine > 1; + + // ── уровень 0 ── + let probe0 = if probe_on_fine { 0 } else { py * spec.nx + px }; + let flags0 = 1 | if probe_on_fine { 0 } else { 2 }; + let l0 = make_level( + &device, + &bgl_level, + spec.nx, + spec.ny, + &geom0, + &beta0, + probe0 as u32, + flags0, + ); + + let pre = device.create_buffer(&wgpu::BufferDescriptor { + label: Some("pre"), + size: (9 * l0.n * 4) as u64, + usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_DST, + mapped_at_creation: false, + }); + + // ── уровень 1 и AMR ── + let mut l1 = None; + let mut amr = None; + if spec.refine > 1 { + let (ax, bx, ay, by) = spec.patch.ok_or("refine > 1 требует патч")?; + let r = spec.refine; + let body1 = spec.body.refined(r as R, ax as R, ay as R); + let nfx = r * (bx - ax) + 1; + let nfy = r * (by - ay) + 1; + let geom1 = cpu::Geom::build(nfx, nfy, &body1); + let tau0 = 1.0 / (2.0 * spec.beta0); + let tau1 = r as R * (tau0 - 0.5) + 0.5; + let beta1 = vec![(1.0 / (2.0 * tau1)) as f32; nfx]; + let r01 = tau1 / (r as R * tau0); + let patch = cpu::Patch::new(&spec, &geom0, &geom1.solid, r01); + + let probe1 = if probe_on_fine { ((py - ay) * r) * nfx + (px - ax) * r } else { 0 }; + let flags1 = if probe_on_fine { 2 } else { 0 }; + let lvl1 = + make_level(&device, &bgl_level, nfx, nfy, &geom1, &beta1, probe1 as u32, flags1); + + let ghosts: Vec = patch + .ghosts() + .iter() + .map(|g| GGhost { + fine: g.fine, + c00: g.c00, + c10: g.c10, + c01: g.c01, + c11: g.c11, + tx: g.tx as f32, + ty: g.ty as f32, + _p: 0, + }) + .collect(); + let rest: Vec<[u32; 2]> = + patch.restrict_pairs().iter().map(|&(c, fi)| [c, fi]).collect(); + + let ghost_bytes = (ghosts.len() * 9 * 4) as u64; + let gh_old = storage(&device, "gh_old", ghost_bytes); + let gh_new = storage(&device, "gh_new", ghost_bytes); + let ghosts_buf = storage_init(&device, "ghosts", bytemuck::cast_slice(&ghosts)); + let rest_buf = storage_init(&device, "rest", bytemuck::cast_slice(&rest)); + let dummy = storage(&device, "dummy", 16); + + let base = AmrParams { + nghost: ghosts.len() as u32, + nrestrict: rest.len() as u32, + ccount: l0.n as u32, + fcount: lvl1.n as u32, + r01: r01 as f32, + rfc: (1.0 / r01) as f32, + w: 0.0, + _pad: 0.0, + }; + let mk_amr_bg = |src: &wgpu::Buffer, + dst: &wgpu::Buffer, + a: &wgpu::Buffer, + b: &wgpu::Buffer, + params: AmrParams| { + let ub = device.create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: Some("amr params"), + contents: bytemuck::bytes_of(¶ms), + usage: wgpu::BufferUsages::UNIFORM, + }); + device.create_bind_group(&wgpu::BindGroupDescriptor { + label: Some("amr"), + layout: &bgl_amr, + entries: &[ + bind(0, src), + bind(1, dst), + bind(2, &ghosts_buf), + bind(3, a), + bind(4, b), + bind(5, &rest_buf), + bind(6, &ub), + ], + }) + }; + + let bg_ghost_old = mk_amr_bg(&pre, &gh_old, &dummy, &dummy, base); + let bg_ghost_new = mk_amr_bg(&l0.f, &gh_new, &dummy, &dummy, base); + let bg_fill: Vec = (0..r) + .map(|s| { + let mut p = base; + p.w = ((s + 1) as R / r as R) as f32; + mk_amr_bg(&dummy, &lvl1.f, &gh_old, &gh_new, p) + }) + .collect(); + let bg_restrict = mk_amr_bg(&lvl1.f, &l0.f, &dummy, &dummy, base); + + amr = Some(AmrRes { + nghost: ghosts.len() as u32, + nrestrict: rest.len() as u32, + bg_ghost_old, + bg_ghost_new, + bg_fill, + bg_restrict, + }); + l1 = Some(lvl1); + } + + let pipes = Pipes { + collide: mk(&pl_level, "k_collide"), + stream: mk(&pl_level, "k_stream"), + bouzidi: mk(&pl_level, "k_bouzidi"), + walls: mk(&pl_level, "k_walls"), + channel: mk(&pl_level, "k_channel"), + force: mk(&pl_level, "k_force"), + stats1: mk(&pl_level, "k_stats1"), + stats2: mk(&pl_level, "k_stats2"), + ghost: mk(&pl_amr, "k_ghost"), + fill: mk(&pl_amr, "k_fill"), + restrict: mk(&pl_amr, "k_restrict"), + empty_bg: device.create_bind_group(&wgpu::BindGroupDescriptor { + label: Some("empty"), + layout: &bgl_empty, + entries: &[], + }), + }; + + let readback = device.create_buffer(&wgpu::BufferDescriptor { + label: Some("readback"), + size: (9 * l0.n * 4) as u64, + usage: wgpu::BufferUsages::MAP_READ | wgpu::BufferUsages::COPY_DST, + mapped_at_creation: false, + }); + + let d_ref = match spec.patch { + Some(_) if spec.refine > 1 => spec.body.d * spec.refine as R, + _ => spec.body.d, + }; + + Ok(Sim { + spec, + device, + queue, + adapter_name, + l0, + l1, + pre, + dyn_buf, + results_buf, + staging, + dyn_binds, + amr, + pipes, + readback, + fluid_count, + step_index: 0, + d_ref, + }) + } + + pub fn name(&self) -> &'static str { + // сигнатура бэкенда фиксирована; конкретный адаптер печатается отдельно + Box::leak(format!("GPU: {} , f32", self.adapter_name).into_boxed_str()) + } + + pub fn force_ref_size(&self) -> R { + self.d_ref + } + + fn inlet(&self, t: u64) -> (R, R) { + let sp = &self.spec; + let u = sp.units.u_lat * math::smoothstep(t as R / sp.ramp.max(1) as R); + let (s, c) = sp.flow_angle.sin_cos(); + let mut uy = u * s; + if sp.pert_dur > 0 && t >= sp.ramp && t < sp.ramp + sp.pert_dur { + let ph = (t - sp.ramp) as R / sp.pert_dur as R; + uy += sp.pert_amp * sp.units.u_lat * (std::f64::consts::PI * ph).sin(); + } + (u * c, uy) + } + + pub fn step(&mut self) -> StepRec { + let t = self.step_index; + let (ux_in, uy_in) = self.inlet(t); + let refine = self.spec.refine.max(1); + self.queue.write_buffer( + &self.dyn_buf, + 0, + bytemuck::bytes_of(&Dyn { + ux_in: ux_in as f32, + uy_in: uy_in as f32, + rho_out: 1.0, + _pad0: 0.0, + outlet_extrap: self.spec.outlet_extrapolate as u32, + nparts: self.l0.parts_count, + refine: refine as u32, + collision: (self.spec.collision == Collision::Bgk) as u32, + }), + ); + + let mut enc = self + .device + .create_command_encoder(&wgpu::CommandEncoderDescriptor { label: Some("step") }); + + // «до» — состояние L0 перед столкновением, нужно как старый край рамки патча + if self.amr.is_some() { + enc.copy_buffer_to_buffer(&self.l0.f, 0, &self.pre, 0, (9 * self.l0.n * 4) as u64); + } + + { + let mut p = enc.begin_compute_pass(&wgpu::ComputePassDescriptor { + label: Some("L0"), + timestamp_writes: None, + }); + p.set_bind_group(0, &self.l0.bind, &[]); + p.set_bind_group(1, &self.dyn_binds[0], &[]); + let ncell = ceil_div(self.l0.n as u32, WG); + p.set_pipeline(&self.pipes.collide); + p.dispatch_workgroups(ncell, 1, 1); + p.set_pipeline(&self.pipes.stream); + p.dispatch_workgroups(ncell, 1, 1); + p.set_pipeline(&self.pipes.bouzidi); + p.dispatch_workgroups(self.l0.link_groups(), 1, 1); + // порядок обязателен: стенки снимают заворот по y, затем Zou-He — по x, на весь столбец + p.set_pipeline(&self.pipes.walls); + p.dispatch_workgroups(ceil_div(self.l0.nx as u32, WG), 1, 1); + p.set_pipeline(&self.pipes.channel); + p.dispatch_workgroups(ceil_div(self.l0.ny as u32, WG), 1, 1); + } + + if let (Some(l1), Some(a)) = (self.l1.as_ref(), self.amr.as_ref()) { + { + let mut p = enc.begin_compute_pass(&wgpu::ComputePassDescriptor { + label: Some("ghost"), + timestamp_writes: None, + }); + p.set_bind_group(0, &self.pipes.empty_bg, &[]); + p.set_bind_group(1, &self.pipes.empty_bg, &[]); + p.set_pipeline(&self.pipes.ghost); + p.set_bind_group(2, &a.bg_ghost_old, &[]); + p.dispatch_workgroups(ceil_div(a.nghost, WG), 1, 1); + p.set_bind_group(2, &a.bg_ghost_new, &[]); + p.dispatch_workgroups(ceil_div(a.nghost, WG), 1, 1); + } + let ncell1 = ceil_div(l1.n as u32, WG); + let nlink1 = l1.link_groups(); + for s in 0..refine { + { + let mut p = enc.begin_compute_pass(&wgpu::ComputePassDescriptor { + label: Some("L1"), + timestamp_writes: None, + }); + p.set_bind_group(0, &l1.bind, &[]); + p.set_bind_group(1, &self.dyn_binds[s], &[]); + p.set_pipeline(&self.pipes.collide); + p.dispatch_workgroups(ncell1, 1, 1); + p.set_pipeline(&self.pipes.stream); + p.dispatch_workgroups(ncell1, 1, 1); + p.set_pipeline(&self.pipes.bouzidi); + p.dispatch_workgroups(nlink1, 1, 1); + // силу снимаем на каждом подшаге, усредняется она в k_stats2 + p.set_pipeline(&self.pipes.force); + p.dispatch_workgroups(1, 1, 1); + } + { + let mut p = enc.begin_compute_pass(&wgpu::ComputePassDescriptor { + label: Some("fill"), + timestamp_writes: None, + }); + p.set_bind_group(0, &self.pipes.empty_bg, &[]); + p.set_bind_group(1, &self.pipes.empty_bg, &[]); + p.set_bind_group(2, &a.bg_fill[s], &[]); + p.set_pipeline(&self.pipes.fill); + p.dispatch_workgroups(ceil_div(a.nghost, WG), 1, 1); + } + } + { + let mut p = enc.begin_compute_pass(&wgpu::ComputePassDescriptor { + label: Some("restrict"), + timestamp_writes: None, + }); + p.set_bind_group(0, &self.pipes.empty_bg, &[]); + p.set_bind_group(1, &self.pipes.empty_bg, &[]); + p.set_bind_group(2, &a.bg_restrict, &[]); + p.set_pipeline(&self.pipes.restrict); + p.dispatch_workgroups(ceil_div(a.nrestrict, WG), 1, 1); + } + } else { + let mut p = enc.begin_compute_pass(&wgpu::ComputePassDescriptor { + label: Some("force L0"), + timestamp_writes: None, + }); + p.set_bind_group(0, &self.l0.bind, &[]); + p.set_bind_group(1, &self.dyn_binds[0], &[]); + p.set_pipeline(&self.pipes.force); + p.dispatch_workgroups(1, 1, 1); + } + + { + let mut p = enc.begin_compute_pass(&wgpu::ComputePassDescriptor { + label: Some("stats"), + timestamp_writes: None, + }); + p.set_bind_group(1, &self.dyn_binds[0], &[]); + p.set_bind_group(0, &self.l0.bind, &[]); + p.set_pipeline(&self.pipes.stats1); + p.dispatch_workgroups(self.l0.parts_count, 1, 1); + if let Some(l1) = self.l1.as_ref() { + // на тонком уровне stats1 нужен только чтобы снять зонд (флаг записи сумм снят) + p.set_bind_group(0, &l1.bind, &[]); + p.dispatch_workgroups(ceil_div(l1.n as u32, WG), 1, 1); + p.set_bind_group(0, &self.l0.bind, &[]); + } + p.set_pipeline(&self.pipes.stats2); + p.dispatch_workgroups(1, 1, 1); + } + + enc.copy_buffer_to_buffer( + &self.results_buf, + 0, + &self.staging, + 0, + std::mem::size_of::() as u64, + ); + self.queue.submit(Some(enc.finish())); + + let res: Results = self.read_staging(); + self.step_index += 1; + + let cnt = if res.cnt > 0.0 { res.cnt as R } else { self.fluid_count }; + StepRec { + step: t, + fx: res.fx as R, + fy: res.fy as R, + tz: res.tz as R, + uy_probe: res.uy_probe as R, + rho_mean: res.rho_sum as R / self.fluid_count, + max_u: res.max_u as R, + gamma_mean: res.g_sum as R / cnt, + gamma_min: res.g_min as R, + gamma_max: res.g_max as R, + degenerate_frac: res.degen as R / cnt, + xi_negative_frac: res.xineg as R / cnt, + } + } + + fn read_staging(&self) -> Results { + let slice = self.staging.slice(..); + let (tx, rx) = std::sync::mpsc::channel(); + slice.map_async(wgpu::MapMode::Read, move |r| { + let _ = tx.send(r); + }); + self.device.poll(wgpu::Maintain::Wait); + let out = match rx.recv() { + Ok(Ok(())) => { + let data = slice.get_mapped_range(); + *bytemuck::from_bytes::(&data[..std::mem::size_of::()]) + } + _ => Results::default(), + }; + self.staging.unmap(); + out + } + + /// Скачать популяции L0 на хост (нужно для кадров и проверки на NaN). + fn download_l0(&self) -> Vec { + self.download(&self.l0.f, (9 * self.l0.n * 4) as u64) + } + + /// Скачать произвольный буфер уровня. Буфер приёма выделен один раз при сборке под самый + /// большой запрос (популяции L0); для меньших читается только его начало. + fn download(&self, src: &wgpu::Buffer, bytes: u64) -> Vec { + let rb = &self.readback; + let mut enc = self + .device + .create_command_encoder(&wgpu::CommandEncoderDescriptor { label: Some("dl") }); + enc.copy_buffer_to_buffer(src, 0, rb, 0, bytes); + self.queue.submit(Some(enc.finish())); + + let slice = rb.slice(..); + let (tx, rx) = std::sync::mpsc::channel(); + slice.map_async(wgpu::MapMode::Read, move |r| { + let _ = tx.send(r); + }); + self.device.poll(wgpu::Maintain::Wait); + let out = match rx.recv() { + Ok(Ok(())) => { + let m = slice.get_mapped_range(); + bytemuck::cast_slice::(&m[..bytes as usize]).to_vec() + } + _ => vec![f32::NAN; (bytes / 4) as usize], + }; + rb.unmap(); + out + } + + pub fn is_finite(&self) -> bool { + self.download_l0().iter().all(|v| v.is_finite()) + } + + pub fn sample_field(&self, kind: FieldKind) -> (Vec, Vec) { + let (nx, ny, n) = (self.l0.nx, self.l0.ny, self.l0.n); + let raw = self.download_l0(); + let node = |k: usize| -> [R; 9] { + std::array::from_fn(|i| raw[i * n + k] as R) + }; + let mut out = vec![0.0; n]; + match kind { + FieldKind::Speed => { + for k in 0..n { + let (_, ux, uy) = math::macros(&node(k)); + out[k] = (ux * ux + uy * uy).sqrt(); + } + } + FieldKind::Density => { + for k in 0..n { + out[k] = math::macros(&node(k)).0; + } + } + FieldKind::Vorticity => { + let mut ux = vec![0.0; n]; + let mut uy = vec![0.0; n]; + for k in 0..n { + let (_, a, b) = math::macros(&node(k)); + ux[k] = a; + uy[k] = b; + } + for y in 0..ny { + for x in 0..nx { + let xp = (x + 1).min(nx - 1); + let xm = x.saturating_sub(1); + let yp = (y + 1).min(ny - 1); + let ym = y.saturating_sub(1); + let dvdx = (uy[y * nx + xp] - uy[y * nx + xm]) / (xp - xm).max(1) as R; + let dudy = (ux[yp * nx + x] - ux[ym * nx + x]) / (yp - ym).max(1) as R; + out[y * nx + x] = dvdx - dudy; + } + } + } + FieldKind::Gamma => { + // берём ровно то, что посчитало устройство: пересчёт на хосте с базовым β + // соврал бы всюду, где β — поле (включённая губка) + let g = self.download(&self.l0.gam, (n * 4) as u64); + for k in 0..n { + out[k] = g[k] as R; + } + } + } + let solid = cpu::Geom::build(nx, ny, &self.spec.body).solid; + (out, solid) + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Мелкие помощники +// ───────────────────────────────────────────────────────────────────────────── + +fn ceil_div(a: u32, b: u32) -> u32 { + (a + b - 1) / b +} + +fn bind(binding: u32, buf: &wgpu::Buffer) -> wgpu::BindGroupEntry<'_> { + wgpu::BindGroupEntry { binding, resource: buf.as_entire_binding() } +} + +fn storage(device: &wgpu::Device, label: &str, size: u64) -> wgpu::Buffer { + device.create_buffer(&wgpu::BufferDescriptor { + label: Some(label), + size: size.max(16), + usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_DST, + mapped_at_creation: false, + }) +} + +fn storage_init(device: &wgpu::Device, label: &str, data: &[u8]) -> wgpu::Buffer { + let pad; + let data = if data.is_empty() { + pad = [0u8; 16]; + &pad[..] + } else { + data + }; + device.create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: Some(label), + contents: data, + usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_DST, + }) +} + +fn ro(binding: u32) -> wgpu::BindGroupLayoutEntry { + entry(binding, wgpu::BufferBindingType::Storage { read_only: true }) +} +fn rw(binding: u32) -> wgpu::BindGroupLayoutEntry { + entry(binding, wgpu::BufferBindingType::Storage { read_only: false }) +} +fn un(binding: u32) -> wgpu::BindGroupLayoutEntry { + entry(binding, wgpu::BufferBindingType::Uniform) +} +fn entry(binding: u32, ty: wgpu::BufferBindingType) -> wgpu::BindGroupLayoutEntry { + wgpu::BindGroupLayoutEntry { + binding, + visibility: wgpu::ShaderStages::COMPUTE, + ty: wgpu::BindingType::Buffer { ty, has_dynamic_offset: false, min_binding_size: None }, + count: None, + } +} + +fn level_layout(device: &wgpu::Device) -> wgpu::BindGroupLayout { + device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor { + label: Some("level"), + entries: &[rw(0), rw(1), rw(2), ro(3), ro(4), ro(5), rw(6), un(7)], + }) +} +fn dyn_layout(device: &wgpu::Device) -> wgpu::BindGroupLayout { + device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor { + label: Some("dyn"), + entries: &[un(0), un(1), rw(2)], + }) +} +fn amr_layout(device: &wgpu::Device) -> wgpu::BindGroupLayout { + device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor { + label: Some("amr"), + entries: &[ro(0), rw(1), ro(2), ro(3), ro(4), ro(5), un(6)], + }) +} diff --git a/docs/theory/2d_solver/src/main.rs b/docs/theory/2d_solver/src/main.rs new file mode 100644 index 0000000..964d883 --- /dev/null +++ b/docs/theory/2d_solver/src/main.rs @@ -0,0 +1,852 @@ +//! ОСНОВНОЙ ФАЙЛ: разбор параметров, сборка задачи, цикл по шагам, вывод. +//! +//! Здесь же живёт контракт между оркестратором и бэкендами: `Spec` (что считать), +//! `StepRec` (что бэкенд отдаёт за шаг) и `FieldKind` (что рисовать). Бэкенды (cpu.rs, +//! gpu.rs) реализуют один и тот же интерфейс и взаимозаменяемы флагом `--backend`. +//! +//! Раскладка задачи — канал с обтекаемым телом: вход слева (Zou–He по скорости), выход +//! справа (Zou–He по давлению), сверху и снизу зеркальные (free-slip) стенки, тело — +//! no-slip через SDF + интерполированный отскок Bouzidi. + +mod cpu; +mod gif; +mod math; +#[cfg(feature = "gpu")] +mod gpu; + +use std::io::Write as _; +use std::time::Instant; + +use clap::Parser; + +use math::{Body, ShapeKind, Units, R}; + +// ───────────────────────────────────────────────────────────────────────────── +// Контракт оркестратор ↔ бэкенд +// ───────────────────────────────────────────────────────────────────────────── + +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum Collision { + /// Энтропийный KBC-D (он же N1) — основной оператор модели. + Kbc, + /// Классический LBGK — для сравнительных прогонов; это ровно KBC при γ ≡ 2. + Bgk, +} + +/// Что рисовать в анимации. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum FieldKind { + Speed, + Vorticity, + Density, + Gamma, +} + +impl FieldKind { + pub fn label(&self) -> &'static str { + match self { + FieldKind::Speed => "SPEED", + FieldKind::Vorticity => "VORTICITY", + FieldKind::Density => "DENSITY", + FieldKind::Gamma => "GAMMA", + } + } +} + +/// Запись за один шаг — сырьё для всего отчёта. +#[derive(Clone, Copy, Debug, Default)] +pub struct StepRec { + pub step: u64, + pub fx: R, + pub fy: R, + pub tz: R, + pub uy_probe: R, + pub rho_mean: R, + pub max_u: R, + pub gamma_mean: R, + pub gamma_min: R, + pub gamma_max: R, + pub degenerate_frac: R, + pub xi_negative_frac: R, +} + +/// Полностью разрешённая постановка задачи в решёточных единицах. +#[derive(Clone)] +pub struct Spec { + pub nx: usize, + pub ny: usize, + pub body: Body, + pub refine: usize, + /// Границы патча измельчения в координатах L0: (ax, bx, ay, by). + pub patch: Option<(usize, usize, usize, usize)>, + pub units: Units, + pub beta0: R, + pub steps: u64, + pub ramp: u64, + /// Направление набегающего потока, рад (0 — вдоль +x). + pub flow_angle: R, + pub pert_amp: R, + pub pert_dur: u64, + pub outlet_extrapolate: bool, + pub sponge_len: usize, + pub sponge_mult: R, + pub collision: Collision, + /// Узел зонда следа в координатах L0. + pub probe: (usize, usize), +} + +// ───────────────────────────────────────────────────────────────────────────── +// Параметры запуска +// ───────────────────────────────────────────────────────────────────────────── + +#[derive(Parser, Debug)] +#[command( + name = "kbc2d", + about = "Двумерный решатель LBM D2Q9 с энтропийным столкновением KBC (модель D/N1)", + long_about = "Обтекание тела в канале: KBC-D по Bösch/Chikatamarla/Karlin, SDF+Bouzidi на \ +теле, Zou–He вход/выход, вложенный патч измельчения. Анимация синхронизируется с ФИЗИЧЕСКИМ \ +временем потока, а не со скоростью счёта." +)] +struct Cli { + // ── физика потока ── + /// Скорость набегающего потока, м/с + #[arg(long, default_value_t = 30.0, help_heading = "Физика")] + u_phys: R, + /// Размер ячейки (плотность сетки), м + #[arg(long, default_value_t = 0.1, help_heading = "Физика")] + dx: R, + /// Решёточная скорость потока (ячеек за шаг). Задаёт δt = u_lat·dx/u_phys. + /// 0.05 ⇒ Ma≈0.09. Значение 1.0 воспроизводит арифметику «30 м/с ÷ 0.1 м = 300 шагов/с», + /// но это Ma≈1.73 — за пределами применимости LBM. + #[arg(long, default_value_t = 0.05, help_heading = "Физика")] + u_lat: R, + /// Число Рейнольдса по характерному размеру тела + #[arg(long, default_value_t = 150.0, help_heading = "Физика")] + re: R, + /// Направление потока, градусы (0 — вдоль канала) + #[arg(long, default_value_t = 0.0, help_heading = "Физика")] + flow_angle: R, + + // ── сетка ── + /// Размер домена по x, ячеек + #[arg(long, default_value_t = 480, help_heading = "Сетка")] + nx: usize, + /// Размер домена по y, ячеек + #[arg(long, default_value_t = 240, help_heading = "Сетка")] + ny: usize, + /// Коэффициент измельчения вложенного патча (1 — патча нет) + #[arg(long, default_value_t = 2, help_heading = "Сетка")] + refine: usize, + /// Границы патча в координатах L0: ax,bx,ay,by (по умолчанию строятся вокруг тела) + #[arg(long, help_heading = "Сетка")] + patch: Option, + + // ── тело ── + /// Форма обтекаемого тела + #[arg(long, default_value = "cylinder", + value_parser = ShapeKind::ALL, help_heading = "Тело")] + shape: String, + /// Характерный размер тела (диаметр/сторона/хорда), ячеек + #[arg(long, default_value_t = 24.0, help_heading = "Тело")] + size: R, + /// Относительная толщина для ellipse / naca / plate + #[arg(long, default_value_t = 0.3, help_heading = "Тело")] + thickness: R, + /// Угол атаки тела, градусы + #[arg(long, default_value_t = 0.0, help_heading = "Тело")] + body_angle: R, + /// Положение центра тела по x, ячеек (по умолчанию nx/4) + #[arg(long, help_heading = "Тело")] + body_x: Option, + /// Положение центра тела по y, ячеек (по умолчанию ny/2) + #[arg(long, help_heading = "Тело")] + body_y: Option, + + // ── время ── + /// Число шагов симуляции + #[arg(long, default_value_t = 20000, help_heading = "Время")] + steps: u64, + /// Длина smoothstep-разгона входа, шагов + #[arg(long, default_value_t = 1000, help_heading = "Время")] + ramp: u64, + /// Амплитуда стартового поперечного импульса (доля от U); 0 — выключить + #[arg(long, default_value_t = 0.1, help_heading = "Время")] + pert_amp: R, + /// Длительность стартового импульса, шагов + #[arg(long, default_value_t = 300, help_heading = "Время")] + pert_dur: u64, + + // ── численная схема ── + /// Оператор столкновения + #[arg(long, default_value = "kbc", value_parser = ["kbc", "bgk"], help_heading = "Схема")] + collision: String, + /// Поперечная скорость на выходе: extrapolate — нуль-градиент (умолчание: в факторном + /// исследовании эталонного решателя это лучший вариант по всем метрикам), zero — жёсткий + /// ноль (классический Zou–He; отражает вихри дорожки назад к телу и завышает rms Cl) + #[arg(long, default_value = "extrapolate", value_parser = ["zero", "extrapolate"], + help_heading = "Схема")] + outlet: String, + /// Длина поглощающей губки перед выходом, столбцов (0 — выключена) + #[arg(long, default_value_t = 0, help_heading = "Схема")] + sponge_len: usize, + /// Во сколько раз губка поднимает вязкость + #[arg(long, default_value_t = 30.0, help_heading = "Схема")] + sponge_mult: R, + /// Бэкенд вычислений + #[arg(long, default_value = "cpu", value_parser = ["cpu", "gpu"], help_heading = "Схема")] + backend: String, + /// Число потоков CPU (0 — по числу ядер) + #[arg(long, default_value_t = 0, help_heading = "Схема")] + threads: usize, + + // ── анимация ── + /// Файл гифки (не указан — анимация не пишется) + #[arg(long, help_heading = "Анимация")] + gif: Option, + /// Шагов на кадр: число либо auto (подобрать под --gif-fps в реальном времени) + #[arg(long, default_value = "auto", help_heading = "Анимация")] + gif_every: String, + /// Целевая частота кадров для режима auto + #[arg(long, default_value_t = 30.0, help_heading = "Анимация")] + gif_fps: R, + /// Скорость воспроизведения: 1 — реальное время, 0.1 — замедление в 10 раз + #[arg(long, default_value_t = 1.0, help_heading = "Анимация")] + gif_speed: R, + /// Что рисовать + #[arg(long, default_value = "vorticity", + value_parser = ["speed", "vorticity", "density", "gamma"], help_heading = "Анимация")] + gif_field: String, + /// Палитра + #[arg(long, default_value = "coolwarm", + value_parser = gif::ColorMap::ALL, help_heading = "Анимация")] + gif_cmap: String, + /// Целочисленное увеличение картинки + #[arg(long, default_value_t = 2, help_heading = "Анимация")] + gif_scale: usize, + /// Диапазон нормировки цвета: lo,hi (по умолчанию — по скорости потока) + #[arg(long, help_heading = "Анимация")] + gif_range: Option, + /// Не рисовать служебную надпись на кадре + #[arg(long, help_heading = "Анимация")] + no_hud: bool, + + // ── вывод ── + /// Период строк живого отчёта, шагов (0 — выключить) + #[arg(long, default_value_t = 1000, help_heading = "Вывод")] + report_every: u64, + /// Подробность: quiet — только итог, normal — живой отчёт, full — плюс разбор KBC + #[arg(long, default_value = "normal", value_parser = ["quiet", "normal", "full"], + help_heading = "Вывод")] + verbose: String, + /// Число окон в отчёте о сходимости + #[arg(long, default_value_t = 10, help_heading = "Вывод")] + windows: usize, + /// Выгрузить временные ряды в CSV + #[arg(long, help_heading = "Вывод")] + csv: Option, +} + +// ───────────────────────────────────────────────────────────────────────────── +// Сборка постановки +// ───────────────────────────────────────────────────────────────────────────── + +fn parse_pair(s: &str, what: &str) -> Result<(R, R), String> { + let p: Vec<&str> = s.split(',').collect(); + if p.len() != 2 { + return Err(format!("{what}: ожидалось два числа через запятую, получено «{s}»")); + } + let a = p[0].trim().parse::().map_err(|e| format!("{what}: {e}"))?; + let b = p[1].trim().parse::().map_err(|e| format!("{what}: {e}"))?; + Ok((a, b)) +} + +fn build_spec(cli: &Cli) -> Result { + if cli.nx < 16 || cli.ny < 16 { + return Err("сетка меньше 16×16 не имеет смысла".into()); + } + if cli.u_lat <= 0.0 || cli.u_lat >= 1.0 { + return Err(format!("--u-lat обязан лежать в (0,1), получено {}", cli.u_lat)); + } + if cli.u_phys <= 0.0 || cli.dx <= 0.0 || cli.re <= 0.0 { + return Err("--u-phys, --dx и --re обязаны быть положительными".into()); + } + if cli.refine == 0 { + return Err("--refine обязан быть ≥ 1".into()); + } + + let kind = ShapeKind::from_str(&cli.shape).ok_or("неизвестная форма")?; + let cx = cli.body_x.unwrap_or(cli.nx as R / 4.0); + let cy = cli.body_y.unwrap_or(cli.ny as R / 2.0); + let body = Body::new(kind, cx, cy, cli.size, cli.body_angle, cli.thickness); + + let units = Units::new(cli.dx, cli.u_phys, cli.u_lat, cli.re, cli.size); + let beta0 = math::beta_of_nu(units.nu_lat); + let tau0 = 1.0 / (2.0 * beta0); + if tau0 <= 0.5 { + return Err(format!( + "τ = {tau0:.4} ≤ ½: вязкость неположительна. Подними --u-lat или опусти --re" + )); + } + + // патч измельчения: по умолчанию охватывает тело и ближний след + let patch = if cli.refine > 1 { + let d = cli.size; + let (ax, bx, ay, by) = match &cli.patch { + Some(s) => { + let p: Vec<&str> = s.split(',').collect(); + if p.len() != 4 { + return Err("--patch: ожидалось ax,bx,ay,by".into()); + } + let v: Result, _> = + p.iter().map(|t| t.trim().parse::()).collect(); + let v = v.map_err(|e| format!("--patch: {e}"))?; + (v[0], v[1], v[2], v[3]) + } + None => ( + (cx - 1.5 * d).round().max(2.0) as usize, + (cx + 6.25 * d).round().min(cli.nx as R - 3.0) as usize, + (cy - 2.0 * d).round().max(1.0) as usize, + (cy + 2.0 * d).round().min(cli.ny as R - 2.0) as usize, + ), + }; + // патч обязан лежать строго внутри домена и НЕ накрывать зеркальные ряды стенок: + // на них ГУ работает по всему ряду, а рамка патча их бы перезаписала + if !(1 <= ay && ay < by && by <= cli.ny - 2) { + return Err(format!( + "патч по y [{ay},{by}] обязан лежать внутри (0,{}) и не трогать ряды стенок; \ + при теле {d} ячеек минимальное ny ≈ {}", + cli.ny - 1, + (4.0 * d) as usize + 6 + )); + } + if !(2 <= ax && ax < bx && bx <= cli.nx - 2) { + return Err(format!("патч по x [{ax},{bx}] обязан лежать внутри (1,{})", cli.nx - 1)); + } + if cli.sponge_len > 0 && cli.nx - 1 - cli.sponge_len <= bx { + return Err(format!( + "губка (последние {} столбцов) накрывает патч (bx={bx})", + cli.sponge_len + )); + } + Some((ax, bx, ay, by)) + } else { + None + }; + + let probe = ( + ((cx + 3.0 * cli.size).round() as usize).min(cli.nx - 2), + (cy.round() as usize).min(cli.ny - 2), + ); + + Ok(Spec { + nx: cli.nx, + ny: cli.ny, + body, + refine: cli.refine, + patch, + units, + beta0, + steps: cli.steps, + ramp: cli.ramp.max(1), + flow_angle: cli.flow_angle * std::f64::consts::PI / 180.0, + pert_amp: cli.pert_amp, + pert_dur: cli.pert_dur, + outlet_extrapolate: cli.outlet == "extrapolate", + sponge_len: cli.sponge_len, + sponge_mult: cli.sponge_mult, + collision: if cli.collision == "bgk" { Collision::Bgk } else { Collision::Kbc }, + probe, + }) +} + +// ───────────────────────────────────────────────────────────────────────────── +// Бэкенд +// ───────────────────────────────────────────────────────────────────────────── + +enum Backend { + Cpu(Box), + #[cfg(feature = "gpu")] + Gpu(Box), +} + +impl Backend { + fn step(&mut self) -> StepRec { + match self { + Backend::Cpu(s) => s.step(), + #[cfg(feature = "gpu")] + Backend::Gpu(s) => s.step(), + } + } + fn sample_field(&mut self, k: FieldKind) -> (Vec, Vec) { + match self { + Backend::Cpu(s) => { + let (v, m) = s.sample_field(k); + (v, m.to_vec()) + } + #[cfg(feature = "gpu")] + Backend::Gpu(s) => s.sample_field(k), + } + } + fn is_finite(&mut self) -> bool { + match self { + Backend::Cpu(s) => s.is_finite(), + #[cfg(feature = "gpu")] + Backend::Gpu(s) => s.is_finite(), + } + } + fn force_ref_size(&self) -> R { + match self { + Backend::Cpu(s) => s.force_ref_size(), + #[cfg(feature = "gpu")] + Backend::Gpu(s) => s.force_ref_size(), + } + } + fn name(&self) -> &'static str { + match self { + Backend::Cpu(_) => "CPU (rayon, f64)", + #[cfg(feature = "gpu")] + Backend::Gpu(s) => s.name(), + } + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// main +// ───────────────────────────────────────────────────────────────────────────── + +fn main() { + let cli = Cli::parse(); + if let Err(e) = run(cli) { + eprintln!("\nОШИБКА: {e}"); + std::process::exit(1); + } +} + +fn run(cli: Cli) -> Result<(), String> { + if cli.threads > 0 { + rayon::ThreadPoolBuilder::new() + .num_threads(cli.threads) + .build_global() + .map_err(|e| e.to_string())?; + } + let spec = build_spec(&cli)?; + let quiet = cli.verbose == "quiet"; + let full = cli.verbose == "full"; + + // ── анимация: план тайминга считается ДО прогона и печатается в шапке ── + let stride = if cli.gif_every == "auto" { + None + } else { + Some(cli.gif_every.parse::().map_err(|e| format!("--gif-every: {e}"))?) + }; + let plan = gif::GifPlan::new(&spec.units, stride, cli.gif_fps, cli.gif_speed); + let field_kind = match cli.gif_field.as_str() { + "speed" => FieldKind::Speed, + "density" => FieldKind::Density, + "gamma" => FieldKind::Gamma, + _ => FieldKind::Vorticity, + }; + + if !quiet { + print_header(&spec, &cli, &plan, field_kind); + } + + // ── бэкенд ── + let mut back = match cli.backend.as_str() { + "gpu" => { + #[cfg(feature = "gpu")] + { + Backend::Gpu(Box::new(gpu::Sim::new(spec.clone())?)) + } + #[cfg(not(feature = "gpu"))] + { + return Err("сборка без поддержки GPU (--no-default-features)".into()); + } + } + _ => Backend::Cpu(Box::new(cpu::Sim::new(spec.clone()))), + }; + if !quiet { + println!("бэкенд: {}", back.name()); + } + + // ── писатель гифки ── + let mut writer = match &cli.gif { + Some(path) => { + let cmap = gif::ColorMap::from_str(&cli.gif_cmap).ok_or("неизвестная палитра")?; + let range = match &cli.gif_range { + Some(s) => { + let (lo, hi) = parse_pair(s, "--gif-range")?; + gif::Range { lo, hi } + } + None => default_range(field_kind, spec.units.u_lat, spec.units.d_lat), + }; + let hud = gif::Hud { + u_phys: spec.units.u_phys, + field: field_kind.label(), + show: !cli.no_hud, + }; + Some( + gif::GifWriter::create( + path, + spec.nx, + spec.ny, + cli.gif_scale, + cmap, + range, + plan, + hud, + spec.patch, + ) + .map_err(|e| format!("не удалось создать {path}: {e}"))?, + ) + } + None => None, + }; + + // ── цикл ── + let d_ref = back.force_ref_size(); + let mut recs: Vec = Vec::with_capacity(spec.steps as usize + 1); + let t_start = Instant::now(); + let mut blew_up = false; + let nodes_per_step = (spec.nx * spec.ny + + spec + .patch + .map(|(ax, bx, ay, by)| { + spec.refine * (spec.refine * (bx - ax) + 1) * (spec.refine * (by - ay) + 1) + }) + .unwrap_or(0)) as f64; + + for t in 0..=spec.steps { + let rec = back.step(); + recs.push(rec); + + if let Some(w) = writer.as_mut() { + if t % plan.stride == 0 { + let (field, solid) = back.sample_field(field_kind); + let (cd, _, _) = + math::coefficients(rec.fx, rec.fy, rec.tz, spec.units.u_lat, d_ref); + w.push(&field, &solid, spec.units.time_of_step(t), Some(cd)) + .map_err(|e| format!("запись кадра: {e}"))?; + } + } + + if cli.report_every > 0 && !quiet && t % cli.report_every == 0 { + live_line(&spec, &rec, d_ref, t_start.elapsed().as_secs_f64(), nodes_per_step, full); + } + + // развал счёта: проверяем редко, проверка стоит полного прохода по полю + if t % 500 == 0 && !back.is_finite() { + eprintln!("\nСЧЁТ РАЗВАЛИЛСЯ на шаге {t}: в поле появились NaN/inf."); + blew_up = true; + break; + } + } + let wall = t_start.elapsed().as_secs_f64(); + + if let Some(w) = writer { + let frames = w.frames; + w.finish().map_err(|e| format!("закрытие гифки: {e}"))?; + if !quiet { + println!("\nгифка: {} — кадров {frames}", cli.gif.as_ref().unwrap()); + } + } + + if let Some(path) = &cli.csv { + write_csv(path, &recs, &spec, d_ref).map_err(|e| format!("CSV: {e}"))?; + if !quiet { + println!("ряды: {path}"); + } + } + + final_report(&spec, &cli, &recs, d_ref, wall, nodes_per_step, blew_up, &plan, full); + Ok(()) +} + +fn default_range(k: FieldKind, u_lat: R, d_lat: R) -> gif::Range { + match k { + FieldKind::Speed => gif::Range { lo: 0.0, hi: 1.7 * u_lat }, + FieldKind::Vorticity => { + let w = 5.0 * u_lat / d_lat; + gif::Range { lo: -w, hi: w } + } + FieldKind::Density => gif::Range { lo: 0.985, hi: 1.015 }, + FieldKind::Gamma => gif::Range { lo: 0.0, hi: 4.0 }, + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Вывод +// ───────────────────────────────────────────────────────────────────────────── + +fn print_header(spec: &Spec, cli: &Cli, plan: &gif::GifPlan, field: FieldKind) { + let u = &spec.units; + let tau = 1.0 / (2.0 * spec.beta0); + let blockage = spec.body.frontal_extent() / spec.ny as R; + + println!("╔══════════════════════════════════════════════════════════════════════════╗"); + println!("║ KBC-2D — обтекание тела в канале, D2Q9 + энтропийное столкновение KBC-D ║"); + println!("╚══════════════════════════════════════════════════════════════════════════╝"); + + println!("\n── Физика ──────────────────────────────────────────────────────────────────"); + println!(" скорость потока {:>12.4} м/с направление {:>6.1}°", u.u_phys, cli.flow_angle); + println!(" размер ячейки {:>12.5} м домен {:.2} × {:.2} м", + u.dx, spec.nx as R * u.dx, spec.ny as R * u.dx); + println!(" тело {:>12} {:.3} м ({:.0} ячеек), угол атаки {:.1}°", + format!("{:?}", spec.body.kind).to_lowercase(), u.d_phys(), u.d_lat, cli.body_angle); + println!(" Рейнольдс {:>12.1} вязкость {:.4e} м²/с", u.re, u.nu_phys()); + println!(" Мах (реш.) {:>12.4} u_lat = {:.4}, τ = {:.4}", u.mach(), u.u_lat, tau); + println!(" блокировка канала {:>12.3} (габарит тела / ширина канала)", blockage); + + println!("\n── Время ───────────────────────────────────────────────────────────────────"); + println!(" шаг по времени {:>12.6e} с {:.1} шагов на секунду физ. времени", + u.dt, u.steps_per_second()); + println!(" всего шагов {:>12} = {:.4} с физ. времени = {:.1} конв. времён D/U", + spec.steps, u.time_of_step(spec.steps), u.convective_times(spec.steps)); + println!(" разгон / возмущение {:>12} {} шагов импульса амплитудой {:.2}·U", + spec.ramp, spec.pert_dur, spec.pert_amp); + + println!("\n── Сетка ───────────────────────────────────────────────────────────────────"); + println!(" L0 {:>12} {} узлов", format!("{}×{}", spec.nx, spec.ny), + spec.nx * spec.ny); + match spec.patch { + Some((ax, bx, ay, by)) => { + let (nfx, nfy) = (spec.refine * (bx - ax) + 1, spec.refine * (by - ay) + 1); + let tau1 = spec.refine as R * (tau - 0.5) + 0.5; + println!(" L1 (×{}) {:>12} область x∈[{ax},{bx}] y∈[{ay},{by}], τ₁={:.4}", + spec.refine, format!("{nfx}×{nfy}"), tau1); + println!(" {:>12} {} подшагов на шаг L0", "", spec.refine); + } + None => println!(" L1 {:>12} (измельчение выключено)", "нет"), + } + println!(" оператор {:>12}", if spec.collision == Collision::Kbc { "KBC-D" } else { "LBGK" }); + println!(" выход по u_y {:>12} губка {} столбцов ×{:.0}", + if spec.outlet_extrapolate { "extrapolate" } else { "zero" }, + spec.sponge_len, spec.sponge_mult); + + // Пара «вход по скорости / выход по давлению» — недодемпфированный акустический резонатор: + // затухание продольной моды ~ ν(π/Nx)², то есть на длинном домене её почти ничто не гасит, + // и стартовый транзиент способен раскачать её до смещённого режима или до развала счёта. + if spec.sponge_len == 0 && spec.nx >= 250 { + println!(" + ⚠ домен длинный ({} столбцов), губка выключена. Продольная акустическая", spec.nx); + println!(" мода затухает как ν(π/Nx)² и на такой длине почти не гасится: возможен уход"); + println!(" ⟨ρ⟩ от единицы (режим смещён, числа несопоставимы) вплоть до развала счёта."); + println!(" Смотрите ⟨ρ⟩ в отчёте о сходимости; лечится ключом --sponge-len 32."); + } + + if cli.gif.is_some() { + println!("\n── Анимация (синхронизация с физическим временем) ───────────────────────────"); + println!(" поле {:>12} палитра {}", field.label().to_lowercase(), cli.gif_cmap); + println!(" шаг кадра {:>12} {}", plan.stride, + if plan.auto { "подобран автоматически" } else { "задан явно" }); + println!(" между кадрами {:>12.6e} с физ. времени", plan.frame_dt); + println!(" задержка кадра {:>12.4} сотых {:.2} кадр/с{}", plan.exact_delay_cs, + plan.fps, + if plan.dithered { " (задержки чередуются — точное время без ухода)" } else { "" }); + println!(" воспроизведение {:>12.4}× от реального времени (запрошено {:.4}×)", + plan.actual_playback, plan.requested_playback); + match plan.advice(&spec.units) { + Some(a) => println!(" ⚠ {a}"), + None if plan.is_synced() => { + println!(" ✓ синхронно с физическим временем (расхождение {:.4}%)", + (plan.sync_error() - 1.0).abs() * 100.0) + } + None => println!(" ⚠ расхождение с реальным временем {:.3}×", plan.sync_error()), + } + } + println!(); +} + +fn live_line(spec: &Spec, r: &StepRec, d_ref: R, wall: f64, nodes: f64, full: bool) { + let (cd, cl, _) = math::coefficients(r.fx, r.fy, r.tz, spec.units.u_lat, d_ref); + let done = (r.step + 1) as f64; + let sps = done / wall.max(1e-9); + let eta = (spec.steps + 1 - r.step.min(spec.steps)) as f64 / sps.max(1e-9); + let mlups = sps * nodes / 1e6; + // во сколько раз счёт быстрее реального времени + let rt = sps / spec.units.steps_per_second(); + let head = format!( + "шаг {:>8} t={:>9.4}с ⟨ρ⟩={:.5} max|u|={:.4} Cd={:>7.3} Cl={:>7.3}", + r.step, + spec.units.time_of_step(r.step), + r.rho_mean, + r.max_u, + cd, + cl + ); + if full { + println!( + "{head}\n ⟨γ⟩={:.3} [{:.2}…{:.2}] вырожд={:.2}% ξ<0={:.2}% \ + {:.0} шаг/с ({:.1} MLUPS, {:.2}× реального) ETA {:.0}с", + r.gamma_mean, + r.gamma_min, + r.gamma_max, + r.degenerate_frac * 100.0, + r.xi_negative_frac * 100.0, + sps, + mlups, + rt, + eta + ); + } else { + println!("{head} {:.0} шаг/с ETA {:.0}с", sps, eta); + } + std::io::stdout().flush().ok(); +} + +#[allow(clippy::too_many_arguments)] +fn final_report( + spec: &Spec, + cli: &Cli, + recs: &[StepRec], + d_ref: R, + wall: f64, + nodes: f64, + blew_up: bool, + plan: &gif::GifPlan, + full: bool, +) { + let u = &spec.units; + let n = recs.len(); + if n < 8 { + println!("слишком мало шагов для отчёта"); + return; + } + let h = n / 2; // установившийся режим — вторая половина ряда + + let cd: Vec = recs.iter().map(|r| 2.0 * r.fx / (u.u_lat * u.u_lat * d_ref)).collect(); + let cl: Vec = recs.iter().map(|r| 2.0 * r.fy / (u.u_lat * u.u_lat * d_ref)).collect(); + let cm: Vec = + recs.iter().map(|r| 2.0 * r.tz / (u.u_lat * u.u_lat * d_ref * d_ref)).collect(); + let uy: Vec = recs.iter().map(|r| r.uy_probe).collect(); + let rho: Vec = recs.iter().map(|r| r.rho_mean).collect(); + + let (cd_m, _) = math::mean_std(&cd[h..]); + let (_, cl_rms) = math::mean_std(&cl[h..]); + let (cm_m, _) = math::mean_std(&cm[h..]); + let (_, uy_rms) = math::mean_std(&uy[h..]); + let (st, _, _) = math::strouhal(&uy, u.d_lat, u.u_lat, 8192); + + println!("\n╔══════════════════════════════════════════════════════════════════════════╗"); + println!("║ ИТОГОВЫЙ ОТЧЁТ ║"); + println!("╚══════════════════════════════════════════════════════════════════════════╝"); + if blew_up { + println!("⚠ ПРОГОН ОБОРВАН: счёт развалился. Числа ниже относятся к тому, что успело сойтись."); + } + + // ── установившийся режим ── + let beta = spec.body.frontal_extent() / spec.ny as R; + println!("\n── Установившийся режим (вторая половина ряда, шаги {}–{}) ──────────────────", h, n - 1); + println!(" {:<22}{:>10}{:>10}{:>10}{:>10}", "величина", "St", "⟨Cd⟩", "rms Cl", "⟨Cm⟩"); + println!(" {:<22}{:>10.4}{:>10.4}{:>10.4}{:>10.5}", "как посчитано", st, cd_m, cl_rms, cm_m); + // поправка на стеснение канала: St кинематическое ⇒ ×(1−β); Cd и Cl ~ скорость² ⇒ ×(1−β)² + println!( + " {:<22}{:>10.4}{:>10.4}{:>10.4}{:>10}", + "с поправ. на блокир.", + st * (1.0 - beta), + cd_m * (1.0 - beta).powi(2), + cl_rms * (1.0 - beta).powi(2), + "—" + ); + if spec.body.kind == ShapeKind::Cylinder && u.re > 100.0 && u.re < 200.0 { + println!(" {:<22}{:>10.4}{:>10.4}{:>10}{:>10.1}", "литература Re≈150", 0.183, 1.330, "~0.30", 0.0); + } + println!(" rms u_y в зонде = {uy_rms:.5} (узел {:?}); блокировка β = {beta:.3}", spec.probe); + println!(" поправка: St×(1−β), Cd и rms Cl ×(1−β)². Литература — безграничный цилиндр;"); + println!(" ⟨Cm⟩≈0 для симметричного тела под нулевым углом — это контроль симметрии считывания."); + + // ── сходимость ── + let nw = cli.windows.max(2).min(n / 2); + let w = n / nw; + println!("\n── Сходимость по времени (окна по {w} шагов) ────────────────────────────────"); + println!(" {:>4}{:>16}{:>10}{:>10}{:>10}{:>10}{:>10}", "окно", "шаги", "⟨ρ⟩", "⟨Cd⟩", "⟨Cd⟩/ρ", "rms Cl", "rms u_y"); + let mut rows = Vec::new(); + for k in 0..nw { + let s = k * w; + let e = if k + 1 == nw { n } else { (k + 1) * w }; + let (rm, _) = math::mean_std(&rho[s..e]); + let (cdm, _) = math::mean_std(&cd[s..e]); + let (_, clr) = math::mean_std(&cl[s..e]); + let (_, uyr) = math::mean_std(&uy[s..e]); + let cdn = if rm != 0.0 { cdm / rm } else { R::NAN }; + rows.push((rm, cdm, cdn, clr, uyr)); + println!(" {:>4}{:>16}{:>10.5}{:>10.4}{:>10.4}{:>10.4}{:>10.5}", + k + 1, format!("{s}–{e}"), rm, cdm, cdn, clr, uyr); + } + let (rm1, _, cdn1, clr1, _) = rows[nw - 1]; + let (rm0, _, _, clr0, _) = rows[nw - 2]; + print!(" вердикт: "); + if (rm1 - 1.0).abs() > 0.02 { + // «стабильно, но не там»: система села на смещённую ветвь (акустическая накачка массы), + // на ней Cd/Cl несопоставимы с литературой + print!("⟨ρ⟩={rm1:.4} — РЕЖИМ СМЕЩЁН, прогон невалиден для сравнения; "); + } else if (rm1 - rm0).abs() > 1e-3 { + print!("⟨ρ⟩={rm1:.4}, Δ за окно {:+.5} — ДРЕЙФ МАССЫ (правь ГУ выхода); ", rm1 - rm0); + } else { + print!("масса стабильна (⟨ρ⟩={rm1:.4}, Δ {:+.5}); ", rm1 - rm0); + } + if clr1 - clr0 > 1e-3 { + println!("rms Cl ещё растёт ({:+.4}) — не насыщено; дрейф-устойчивый Cd = {cdn1:.3}", clr1 - clr0); + } else { + println!("rms Cl насыщен ({:+.4}); дрейф-устойчивый Cd = {cdn1:.3}", clr1 - clr0); + } + + // ── диагностика KBC ── + if spec.collision == Collision::Kbc { + let (g_m, g_s) = math::mean_std(&recs[h..].iter().map(|r| r.gamma_mean).collect::>()); + let gmin = recs[h..].iter().map(|r| r.gamma_min).fold(R::INFINITY, R::min); + let gmax = recs[h..].iter().map(|r| r.gamma_max).fold(R::NEG_INFINITY, R::max); + let (deg, _) = math::mean_std(&recs[h..].iter().map(|r| r.degenerate_frac).collect::>()); + let (xin, _) = math::mean_std(&recs[h..].iter().map(|r| r.xi_negative_frac).collect::>()); + println!("\n── Энтропийный стабилизатор γ (формула (17) статьи) ─────────────────────────"); + println!(" ⟨γ⟩ = {g_m:.4} ± {g_s:.4} размах по узлам [{gmin:.3} … {gmax:.3}]"); + println!(" вырожденных узлов (⟨Δh|Δh⟩→0, γ подменён на 2 = локальный LBGK): {:.3}%", deg * 100.0); + println!(" узлов с отрицательной объёмной вязкостью ξ = c_s²(1/(γβ) − ½): {:.2}%", xin * 100.0); + if full { + println!(" Пояснение. γ ≡ 2 — это в точности LBGK, поэтому доля вырожденных узлов есть"); + println!(" доля решётки, где KBC ничего не даёт; при относительном пороге она обязана быть"); + println!(" ~0 (абсолютный порог загонял её в 77–99% и молча превращал схему в LBGK)."); + println!(" ξ по формуле (57) у моделей B и D зависит от γ, то есть от точки и времени;"); + println!(" ξ < 0 означает локальное антизатухание акустики. У моделей A и C ξ = ν всегда."); + } + } + + // ── производительность и синхронность ── + let sps = n as f64 / wall.max(1e-9); + println!("\n── Производительность ──────────────────────────────────────────────────────"); + println!(" время счёта {wall:.1} с {sps:.0} шаг/с {:.1} MLUPS", sps * nodes / 1e6); + println!(" физического времени просчитано {:.4} с ⇒ счёт идёт в {:.3}× от реального времени", + u.time_of_step(n as u64 - 1), sps / u.steps_per_second()); + if cli.gif.is_some() { + println!(" гифка воспроизводится в {:.4}× от реального времени — это НЕ зависит от скорости", + plan.actual_playback); + println!(" счёта: задержка кадра берётся из δt, а не из wall-clock."); + } + println!(); +} + +fn write_csv(path: &str, recs: &[StepRec], spec: &Spec, d_ref: R) -> std::io::Result<()> { + let mut f = std::io::BufWriter::new(std::fs::File::create(path)?); + writeln!(f, "step,t_phys_s,cd,cl,cm,uy_probe,rho_mean,max_u,gamma_mean,gamma_min,gamma_max,degenerate_frac,xi_negative_frac")?; + let s = spec.units.u_lat * spec.units.u_lat * d_ref; + for r in recs { + writeln!( + f, + "{},{:.9e},{:.6e},{:.6e},{:.6e},{:.6e},{:.9},{:.6e},{:.6},{:.6},{:.6},{:.6},{:.6}", + r.step, + spec.units.time_of_step(r.step), + 2.0 * r.fx / s, + 2.0 * r.fy / s, + 2.0 * r.tz / (s * d_ref), + r.uy_probe, + r.rho_mean, + r.max_u, + r.gamma_mean, + r.gamma_min, + r.gamma_max, + r.degenerate_frac, + r.xi_negative_frac + )?; + } + Ok(()) +} diff --git a/docs/theory/2d_solver/src/math.rs b/docs/theory/2d_solver/src/math.rs new file mode 100644 index 0000000..c6ee540 --- /dev/null +++ b/docs/theory/2d_solver/src/math.rs @@ -0,0 +1,1013 @@ +//! МАТЕМАТИКА РЕШАТЕЛЯ — всё, что не зависит от бэкенда и от способа хранения полей. +//! +//! Модель: D2Q9, энтропийный оператор столкновения KBC в варианте, который в 2D-статье +//! (Bösch, Chikatamarla, Karlin, «Entropic Multi-Relaxation Models for Simulation of Fluid +//! Turbulence», arXiv:1507.02509, табл. I) называется **KBC D**, а в 3D-работах — **KBC-N1**: +//! сдвиговая часть s несёт натуральные моменты {N, Π_xy}, всё остальное (T, Q_xyy, Q_yxx, A) +//! уходит в h. Ссылки на формулы ниже — по нумерации 2D-статьи. +//! +//! Все функции здесь узловые и чистые: принимают/возвращают [R; Q]. Бэкенды (cpu.rs, gpu.rs) +//! отвечают только за раскладку памяти, обход сетки и параллелизм. WGSL-ядра в gpu.rs — это +//! построчный перенос этого файла на f32; при правке физики надо править оба места. + +#![allow(clippy::needless_range_loop)] + +use std::f64::consts::PI; + +/// Точность CPU-пути. GPU-путь работает в f32 (в WGSL нет f64) — см. заметку о GREL ниже. +pub type R = f64; + +// ───────────────────────────────────────────────────────────────────────────── +// Решётка D2Q9 +// ───────────────────────────────────────────────────────────────────────────── + +pub const Q: usize = 9; + +/// Нумерация направлений: 0:(0,0) 1:E 2:N 3:W 4:S 5:NE 6:NW 7:SW 8:SE. +pub const CX: [i32; Q] = [0, 1, 0, -1, 0, 1, -1, -1, 1]; +pub const CY: [i32; Q] = [0, 0, 1, 0, -1, 1, 1, -1, -1]; +/// Индекс противоположного направления (c_ī = −c_i). +pub const OPP: [usize; Q] = [0, 3, 4, 1, 2, 7, 8, 5, 6]; +pub const W: [R; Q] = [ + 4.0 / 9.0, + 1.0 / 9.0, + 1.0 / 9.0, + 1.0 / 9.0, + 1.0 / 9.0, + 1.0 / 36.0, + 1.0 / 36.0, + 1.0 / 36.0, + 1.0 / 36.0, +]; +/// Квадрат скорости звука решётки. +pub const CS2: R = 1.0 / 3.0; + +/// Порог вырожденности знаменателя γ — ОТНОСИТЕЛЬНЫЙ (доля от ‖Δ‖² = ⟨Δ|Δ⟩), не абсолютный. +/// +/// den = ⟨Δh|Δh⟩ квадратична по неравновесию и физически мала (~1e-7…1e-9 в развитом следе). +/// Абсолютный порог 1e-6 в f32 срабатывал бы на 77–99% узлов и молча подменял γ на 2, т.е. гнал +/// чистый LBGK вместо KBC (в python-версии это измерено: энстрофия сдвигового слоя Re=3e4 +/// уходила на +9%). den — сумма неотрицательных слагаемых, её ОТНОСИТЕЛЬНАЯ точность ~eps типа, +/// поэтому одна и та же дробь годится и для f32, и для f64. +pub const GREL: R = 1e-8; + +// ───────────────────────────────────────────────────────────────────────────── +// Равновесие и макропеременные +// ───────────────────────────────────────────────────────────────────────────── + +/// Энтропийное равновесие в product-form — точный максимизатор энтропии (3) при заданных ρ, ρu, +/// а не усечённый по скорости полином (4): +/// +/// feq_i = W_i·ρ·Π_d (2 − √(1+3u_d²)) · ((2u_d + √(1+3u_d²))/(1 − u_d))^{c_id} +/// +/// Показатели c_id ∈ {−1,0,+1}, поэтому возведение в степень не нужно — берём q или 1/q. +#[inline] +pub fn feq(rho: R, ux: R, uy: R) -> [R; Q] { + // клип скорости: product-form расходится при |u|→1 (делитель 1−u) + let ux = ux.clamp(-0.95, 0.95); + let uy = uy.clamp(-0.95, 0.95); + let sx = (1.0 + 3.0 * ux * ux).sqrt(); + let sy = (1.0 + 3.0 * uy * uy).sqrt(); + let base = rho * (2.0 - sx) * (2.0 - sy); + let qx = (2.0 * ux + sx) / (1.0 - ux); + let qy = (2.0 * uy + sy) / (1.0 - uy); + let ix = 1.0 / qx; + let iy = 1.0 / qy; + [ + base * W[0], + base * W[1] * qx, + base * W[2] * qy, + base * W[3] * ix, + base * W[4] * iy, + base * W[5] * qx * qy, + base * W[6] * ix * qy, + base * W[7] * ix * iy, + base * W[8] * qx * iy, + ] +} + +/// Макропеременные: ρ = Σ f_i, u = (Σ c_i f_i)/ρ. +#[inline] +pub fn macros(f: &[R; Q]) -> (R, R, R) { + let rho = f.iter().sum::(); + let mx = f[1] + f[5] + f[8] - f[3] - f[6] - f[7]; + let my = f[2] + f[5] + f[6] - f[4] - f[7] - f[8]; + (rho, mx / rho, my / rho) +} + +// ───────────────────────────────────────────────────────────────────────────── +// KBC: проектор на сдвиг, энтропийный стабилизатор, столкновение +// ───────────────────────────────────────────────────────────────────────────── + +/// Δs = Ps·Δ — проекция неравновесия на сдвиговую часть модели KBC D: натуральные моменты +/// N = M20 − M02 и Π_xy = M11. +/// +/// Аналитическая форма, прямо из представления популяций через натуральные моменты (10): +/// вклад N сидит только в f(σ,0) (+ρN/4) и f(0,λ) (−ρN/4), вклад Π_xy — только в f(σ,λ) +/// (+σλ·ρΠ_xy/4). Поскольку n_i и p_i линейны по своим моментам, а ρ у f и feq одинакова, +/// разности берутся покомпонентно. Эквивалентно матричному Ps = M⁻¹·diag(0,0,0,0,1,1,0,0,0)·M +/// (проверяется тестом `projector_matches_moment_matrix`), но без матрицы 9×9. +#[inline] +pub fn project_shift(d: &[R; Q]) -> [R; Q] { + // ρΔN = Σ (c_x² − c_y²)Δ_i и ρΔΠ_xy = Σ c_x c_y Δ_i + let dn = d[1] + d[3] - d[2] - d[4]; + let dp = d[5] + d[7] - d[6] - d[8]; + let a = 0.25 * dn; + let b = 0.25 * dp; + [0.0, a, -a, a, -a, b, -b, b, -b] +} + +/// Что вернуло столкновение помимо новых популяций — нужно для диагностики KBC. +#[derive(Clone, Copy, Debug, Default)] +pub struct Kbc { + /// Энтропийный стабилизатор γ данного узла. + pub gamma: R, + /// true, если ⟨Δh|Δh⟩ выродилось и γ подменён на 2 (локально это в точности LBGK). + pub degenerate: bool, +} + +/// Столкновение KBC-D за один узел. f перезаписывается пост-столкновительным состоянием. +/// +/// Шаги дословно по статье: ρ,u → feq (3) → Δ = f − feq → Δs = Ps·Δ (13) → Δh = Δ − Δs → +/// γ* по замкнутой оценке (17) → релаксация к зеркальному состоянию (14): +/// +/// γ* = 1/β − (2 − 1/β)·⟨Δs|Δh⟩ / ⟨Δh|Δh⟩, ⟨X|Y⟩ = Σ X_iY_i/feq_i (16),(17) +/// f' = f − β(2Δs + γΔh) +/// +/// При γ = 2 совпадает с LBGK. Сдвиговые моменты релаксируют с точным 2β = 1/τ при любом γ, +/// поэтому кинематическая вязкость (5) от стабилизатора не зависит. +#[inline] +pub fn collide_node(f: &mut [R; Q], beta: R) -> Kbc { + let (rho, ux, uy) = macros(f); + let fe = feq(rho, ux, uy); + + let mut d = [0.0; Q]; + for i in 0..Q { + d[i] = f[i] - fe[i]; + } + let ds = project_shift(&d); + + let (mut num, mut den, mut nrm) = (0.0, 0.0, 0.0); + for i in 0..Q { + let inv = 1.0 / fe[i]; + let dh = d[i] - ds[i]; + num += ds[i] * dh * inv; + den += dh * dh * inv; + nrm += d[i] * d[i] * inv; + } + + let ok = den > GREL * nrm; + let binv = 1.0 / beta; + let gamma = if ok { binv - (2.0 - binv) * num / den } else { 2.0 }; + + for i in 0..Q { + let dh = d[i] - ds[i]; + f[i] -= beta * (2.0 * ds[i] + gamma * dh); + } + Kbc { gamma, degenerate: !ok } +} + +/// Столкновение LBGK — для сравнительных прогонов (`--collision bgk`). +/// Это ровно KBC при γ ≡ 2, т.е. зеркальное состояние (2). +#[inline] +pub fn collide_node_bgk(f: &mut [R; Q], beta: R) -> Kbc { + let (rho, ux, uy) = macros(f); + let fe = feq(rho, ux, uy); + for i in 0..Q { + f[i] += 2.0 * beta * (fe[i] - f[i]); + } + Kbc { gamma: 2.0, degenerate: false } +} + +/// Кинематическая вязкость по β, формула (5): ν = c_s²(1/(2β) − 1/2). +#[inline] +pub fn nu_of_beta(beta: R) -> R { + CS2 * (1.0 / (2.0 * beta) - 0.5) +} + +/// β по вязкости (обратно к (5)). +#[inline] +pub fn beta_of_nu(nu: R) -> R { + 1.0 / (2.0 * (nu / CS2 + 0.5)) +} + +/// ОБЪЁМНАЯ вязкость модели KBC D, формула (57): ξ = c_s²(1/(γβ) − 1/2). +/// +/// В отличие от моделей A и C (где ξ = ν), у D она зависит от узла и времени вместе с γ. +/// Отрицательная ξ означает локальное антизатухание акустической моды — поэтому в отчёте +/// печатается доля узлов с ξ < 0 (см. report в main.rs). +#[inline] +pub fn xi_of_gamma(gamma: R, beta: R) -> R { + CS2 * (1.0 / (gamma * beta) - 0.5) +} + +// ───────────────────────────────────────────────────────────────────────────── +// Единицы: физика ↔ решётка +// ───────────────────────────────────────────────────────────────────────────── + +/// Перевод между физическими и решёточными единицами. +/// +/// В LBM скорость самой решётки жёстко равна c = δx/δt = 1, поэтому шаг по времени +/// однозначно определяется тем, какую решёточную скорость `u_lat` мы назначаем физическому +/// потоку `u_phys`: +/// +/// δt = u_lat · δx / u_phys [с/шаг] +/// шагов в секунду = 1/δt = u_phys / (u_lat · δx) +/// +/// ВНИМАНИЕ к постановке «30 м/с, ячейка 0.1 м ⇒ 300 шагов/с»: эта арифметика (δt = δx/u_phys) +/// отвечает u_lat = 1, т.е. поток движется ровно на ячейку за шаг. Для LBM это далеко за +/// границей применимости: Ma = u_lat/c_s = √3 ≈ 1.73, сверхзвук, разложение Чепмена–Энскога +/// не работает. Значение воспроизводится буквально флагом `--u-lat 1.0`, но по умолчанию +/// стоит u_lat = 0.05 (Ma ≈ 0.087) — тогда 30 м/с при 0.1 м дают 6000 шагов/с. Гифка +/// синхронизируется с ФИЗИЧЕСКИМ временем в обоих случаях (см. gif.rs). +#[derive(Clone, Copy, Debug)] +pub struct Units { + /// Размер ячейки, м. + pub dx: R, + /// Скорость набегающего потока, м/с. + pub u_phys: R, + /// Она же в решёточных единицах (ячеек за шаг). + pub u_lat: R, + /// Физическая длительность одного шага, с. + pub dt: R, + /// Число Рейнольдса по характерному размеру тела. + pub re: R, + /// Характерный размер тела в ячейках. + pub d_lat: R, + /// Кинематическая вязкость в решёточных единицах. + pub nu_lat: R, +} + +impl Units { + pub fn new(dx: R, u_phys: R, u_lat: R, re: R, d_lat: R) -> Self { + let dt = u_lat * dx / u_phys; + let nu_lat = u_lat * d_lat / re; + Units { dx, u_phys, u_lat, dt, re, d_lat, nu_lat } + } + /// Шагов симуляции на секунду физического времени. + pub fn steps_per_second(&self) -> R { + 1.0 / self.dt + } + /// Физический размер тела, м. + pub fn d_phys(&self) -> R { + self.d_lat * self.dx + } + /// Физическая вязкость, м²/с (= ν_lat·δx²/δt). + pub fn nu_phys(&self) -> R { + self.nu_lat * self.dx * self.dx / self.dt + } + /// Число Маха относительно скорости звука решётки. + pub fn mach(&self) -> R { + self.u_lat / CS2.sqrt() + } + /// Физическое время, прошедшее за `steps` шагов, с. + pub fn time_of_step(&self, steps: u64) -> R { + steps as R * self.dt + } + /// Время в единицах «конвективных времён» D/U — привычная шкала для дорожки Кармана. + pub fn convective_times(&self, steps: u64) -> R { + steps as R * self.u_lat / self.d_lat + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Обтекаемые тела: SDF +// ───────────────────────────────────────────────────────────────────────────── + +/// Соглашение по SDF: φ > 0 в жидкости, φ < 0 внутри тела, |φ| ≈ расстояние до поверхности. +/// Это соглашение использует Bouzidi (доля пересечения q считается как φ/(φ − φ_соседа)), +/// поэтому важна не только правильность знака, но и то, что |∇φ| ≈ 1 у поверхности. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum ShapeKind { + /// Круговой цилиндр диаметра D. Эталон: есть литература по Cd/St. + Cylinder, + /// Квадрат со стороной D, грань перпендикулярна потоку. + Square, + /// Ромб: тот же квадрат, повёрнутый на 45° (диагональ D). + Diamond, + /// Эллипс: большая ось D (вдоль x тела), малая — ratio·D. + Ellipse, + /// Симметричный профиль NACA00xx, хорда D, относительная толщина ratio. + Naca, + /// Равносторонний треугольник, вершина против потока, сторона D. + Triangle, + /// Тонкая пластина: длина D, толщина ratio·D. + Plate, +} + +impl ShapeKind { + pub fn from_str(s: &str) -> Option { + Some(match s.to_ascii_lowercase().as_str() { + "cylinder" | "circle" | "круг" | "цилиндр" => ShapeKind::Cylinder, + "square" | "квадрат" => ShapeKind::Square, + "diamond" | "rhombus" | "ромб" => ShapeKind::Diamond, + "ellipse" | "эллипс" => ShapeKind::Ellipse, + "naca" | "airfoil" | "профиль" => ShapeKind::Naca, + "triangle" | "треугольник" => ShapeKind::Triangle, + "plate" | "пластина" => ShapeKind::Plate, + _ => return None, + }) + } + pub const ALL: [&'static str; 7] = + ["cylinder", "square", "diamond", "ellipse", "naca", "triangle", "plate"]; +} + +/// Тело: форма + положение + масштаб + поворот. Для многоугольных форм контур считается один раз. +#[derive(Clone, Debug)] +pub struct Body { + pub kind: ShapeKind, + /// Центр в решёточных координатах уровня, на котором тело строится. + pub cx: R, + pub cy: R, + /// Характерный размер (диаметр / сторона / хорда) в ячейках того же уровня. + pub d: R, + /// Поворот тела против часовой стрелки, рад (угол атаки). + pub angle: R, + /// Относительная толщина для Ellipse/Naca/Plate. + pub ratio: R, + /// Контур в системе тела (без поворота), для многоугольных форм. + poly: Vec<[R; 2]>, +} + +impl Body { + pub fn new(kind: ShapeKind, cx: R, cy: R, d: R, angle_deg: R, ratio: R) -> Self { + let angle = angle_deg * PI / 180.0; + let poly = build_polygon(kind, d, ratio); + Body { kind, cx, cy, d, angle, ratio, poly } + } + + /// Тот же контур, пересчитанный на уровень с измельчением `r` (координаты и размер ×r). + pub fn refined(&self, r: R, ox: R, oy: R) -> Body { + Body::new( + self.kind, + (self.cx - ox) * r, + (self.cy - oy) * r, + self.d * r, + self.angle * 180.0 / PI, + self.ratio, + ) + } + + /// Знаковое расстояние до поверхности тела в точке (x, y) решёточных координат. + pub fn sdf(&self, x: R, y: R) -> R { + // в систему тела: сдвиг к центру + обратный поворот + let (s, c) = self.angle.sin_cos(); + let px = x - self.cx; + let py = y - self.cy; + let bx = c * px + s * py; + let by = -s * px + c * py; + + match self.kind { + ShapeKind::Cylinder => (bx * bx + by * by).sqrt() - 0.5 * self.d, + ShapeKind::Ellipse => sd_ellipse(bx, by, 0.5 * self.d, 0.5 * self.d * self.ratio), + _ => sd_polygon(bx, by, &self.poly), + } + } + + /// Площадь сечения в ячейках² — нужна для отчёта (блокировка канала). + pub fn frontal_extent(&self) -> R { + match self.kind { + ShapeKind::Cylinder => self.d, + ShapeKind::Ellipse | ShapeKind::Naca | ShapeKind::Plate => { + // проекция повёрнутого габарита на ось y + let a = 0.5 * self.d; + let b = 0.5 * self.d * self.ratio; + 2.0 * ((a * self.angle.sin()).powi(2) + (b * self.angle.cos()).powi(2)).sqrt() + } + _ => { + let mut lo = R::INFINITY; + let mut hi = R::NEG_INFINITY; + let (s, c) = self.angle.sin_cos(); + for p in &self.poly { + let wy = s * p[0] + c * p[1]; + lo = lo.min(wy); + hi = hi.max(wy); + } + hi - lo + } + } + } +} + +/// Контур тела в его собственной системе координат (центр в нуле, без поворота). +fn build_polygon(kind: ShapeKind, d: R, ratio: R) -> Vec<[R; 2]> { + let h = 0.5 * d; + match kind { + ShapeKind::Cylinder | ShapeKind::Ellipse => Vec::new(), // аналитические + ShapeKind::Square => vec![[-h, -h], [h, -h], [h, h], [-h, h]], + ShapeKind::Diamond => vec![[h, 0.0], [0.0, h], [-h, 0.0], [0.0, -h]], + ShapeKind::Plate => { + let t = 0.5 * d * ratio; + vec![[-h, -t], [h, -t], [h, t], [-h, t]] + } + ShapeKind::Triangle => { + // равносторонний, вершина в +x (против потока), центр — в центроиде + let s = d; // сторона + let r_circ = s / 3.0_f64.sqrt(); // радиус описанной окружности + (0..3) + .map(|k| { + let a = 2.0 * PI * k as R / 3.0; + [r_circ * a.cos(), r_circ * a.sin()] + }) + .collect() + } + ShapeKind::Naca => naca_symmetric(d, ratio, 64), + } +} + +/// Симметричный профиль NACA00xx: y_t = 5t·c·(0.2969√ξ − 0.1260ξ − 0.3516ξ² + 0.2843ξ³ − 0.1036ξ⁴), +/// ξ = x/c. Хорда направлена по +x тела, начало отсчёта смещено так, чтобы центр вращения +/// (точка приложения угла атаки) был в c/4 — стандартная аэродинамическая четверть хорды. +fn naca_symmetric(chord: R, t: R, n: usize) -> Vec<[R; 2]> { + let yt = |xi: R| { + 5.0 * t + * chord + * (0.2969 * xi.sqrt() - 0.1260 * xi - 0.3516 * xi * xi + 0.2843 * xi.powi(3) + - 0.1036 * xi.powi(4)) + }; + // косинусное сгущение к носку и хвосту + let xs: Vec = (0..=n) + .map(|k| 0.5 * (1.0 - (PI * k as R / n as R).cos())) + .collect(); + let x0 = 0.25 * chord; // четверть хорды в нуле системы тела + let mut poly = Vec::with_capacity(2 * n); + for &xi in xs.iter() { + // нижняя поверхность, от носка к хвосту (обход по часовой → замкнётся против) + poly.push([xi * chord - x0, -yt(xi)]); + } + for &xi in xs.iter().rev().skip(1) { + poly.push([xi * chord - x0, yt(xi)]); + } + // профиль строился носком в −x; развернём, чтобы носок смотрел против потока (+x → навстречу) + for p in poly.iter_mut() { + p[0] = -p[0]; + } + poly +} + +/// Точное знаковое расстояние до простого многоугольника (Íñigo Quílez). +/// Знак берётся по правилу пересечений (winding), расстояние — минимум по рёбрам. +fn sd_polygon(px: R, py: R, v: &[[R; 2]]) -> R { + let n = v.len(); + debug_assert!(n >= 3); + let mut d = (px - v[0][0]).powi(2) + (py - v[0][1]).powi(2); + let mut s = 1.0; + let mut j = n - 1; + for i in 0..n { + let ex = v[j][0] - v[i][0]; + let ey = v[j][1] - v[i][1]; + let wx = px - v[i][0]; + let wy = py - v[i][1]; + let t = ((wx * ex + wy * ey) / (ex * ex + ey * ey)).clamp(0.0, 1.0); + let bx = wx - ex * t; + let by = wy - ey * t; + d = d.min(bx * bx + by * by); + let c1 = py >= v[i][1]; + let c2 = py < v[j][1]; + let c3 = ex * wy > ey * wx; + if (c1 && c2 && c3) || (!c1 && !c2 && !c3) { + s = -s; + } + j = i; + } + s * d.sqrt() +} + +/// ТОЧНОЕ знаковое расстояние до эллипса: ищем ближайшую точку поверхности итерациями по +/// параметру (метод неподвижной точки через эволюту), затем берём длину невязки со знаком. +/// +/// Дешёвое приближение k(k−1)/|∇k| здесь не годится: |∇φ| у него уходит до ~0.6 уже на +/// расстоянии порядка полуоси, а φ используется не только для знака, но и для доли +/// пересечения Bouzidi q = φ_f/(φ_f − φ_s) — там ошибка масштаба напрямую сдвигает стенку. +fn sd_ellipse(px: R, py: R, a: R, b: R) -> R { + let (sx, sy) = (px.signum(), py.signum()); + let (x, y) = (px.abs(), py.abs()); + // вырожденный центр: ближайшая точка — конец малой полуоси + if x < 1e-12 && y < 1e-12 { + return -a.min(b); + } + let inside = (x / a).powi(2) + (y / b).powi(2) <= 1.0; + + // старт в середине первого квадранта, итерации по (tx, ty) = (cos t, sin t) + let mut tx = std::f64::consts::FRAC_1_SQRT_2; + let mut ty = std::f64::consts::FRAC_1_SQRT_2; + for _ in 0..12 { + // центр кривизны (эволюта) в текущей точке + let ex = (a * a - b * b) * tx.powi(3) / a; + let ey = (b * b - a * a) * ty.powi(3) / b; + let rx = a * tx - ex; + let ry = b * ty - ey; + let qx = x - ex; + let qy = y - ey; + let r = (rx * rx + ry * ry).sqrt(); + let q = (qx * qx + qy * qy).sqrt(); + if q < 1e-300 { + break; + } + tx = ((qx * r / q + ex) / a).clamp(0.0, 1.0); + ty = ((qy * r / q + ey) / b).clamp(0.0, 1.0); + let n = (tx * tx + ty * ty).sqrt(); + if n < 1e-300 { + break; + } + tx /= n; + ty /= n; + } + let cx = a * tx; + let cy = b * ty; + let dist = ((x - cx).powi(2) + (y - cy).powi(2)).sqrt(); + let _ = (sx, sy); // знак по положению точки, а не по квадранту + if inside { + -dist + } else { + dist + } +} + +// ───────────────────────────────────────────────────────────────────────────── +// Граничные условия +// ───────────────────────────────────────────────────────────────────────────── + +/// Тип формулы Bouzidi для конкретного линка (решается один раз при сборке). +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum LinkKind { + /// q < ½ и «дальний» сосед x_f − c_i тоже жидкий: f_ī = 2q·f_i + (1−2q)·f_i(дальний). + Near, + /// q ≥ ½: f_ī = (1/2q)·f_i + (1 − 1/2q)·f_ī. + Far, + /// q < ½, но дальнего жидкого соседа нет — откат на простой отскок f_ī = f_i. + Simple, +} + +/// Линк «жидкий узел → твёрдый сосед» с предвычисленной геометрией пересечения. +#[derive(Clone, Copy, Debug)] +pub struct Link { + /// Плоский индекс жидкого узла (y·nx + x). + pub node: u32, + /// Плоский индекс «дальнего» соседа x_f − c_i (валиден только при LinkKind::Near). + pub far: u32, + /// Направление, уходящее в тело. + pub i: u8, + /// Противоположное направление (то, что вернётся в узел). + pub ib: u8, + pub kind: LinkKind, + /// Доля пересечения q = |x_f→стенка| / |x_f→x_solid| ∈ (0,1), из SDF, клип [0.02, 0.98]. + pub q: R, + /// Линк принадлежит обтекаемому телу (а не иной твёрдой поверхности) — по нему считается сила. + pub body: bool, +} + +/// Доля пересечения по SDF: φ_f / (φ_f − φ_solid). Клип отсекает вырожденные линки, +/// где узел лежит практически на поверхности (иначе 1/(2q) взрывается). +#[inline] +pub fn bouzidi_q(phi_fluid: R, phi_solid: R) -> R { + let den = phi_fluid - phi_solid; + if den.abs() < 1e-12 { + return 0.5; + } + (phi_fluid / den).clamp(0.02, 0.98) +} + +/// Zou–He, скоростной вход (западная грань): задаём u = (ux, uy), ρ плавает. +/// Достраиваем приходящие извне популяции 1, 5, 8 из известных 0, 2, 3, 4, 6, 7. +#[inline] +pub fn zou_he_inlet(f: &mut [R; Q], ux: R, uy: R) { + let rho = (f[0] + f[2] + f[4] + 2.0 * (f[3] + f[6] + f[7])) / (1.0 - ux); + let d = 0.5 * (f[2] - f[4]); + f[1] = f[3] + (2.0 / 3.0) * rho * ux; + f[5] = f[7] - d + (1.0 / 6.0) * rho * ux + 0.5 * rho * uy; + f[8] = f[6] + d + (1.0 / 6.0) * rho * ux - 0.5 * rho * uy; +} + +/// Zou–He, давление-выход (восточная грань): задаём ρ = rho_out, u_x плавает. +/// `uy_out` — поперечная скорость: 0 (классика) либо экстраполяция с предвыходного столбца. +/// +/// Зачем именно давление-выход: грубый zero-gradient по f не выпускает поток, масса и +/// противодавление копятся, поток глохнет. Давление-выход якорит массу, скоростной вход её +/// не пересоздаёт — дрейфа нет. +#[inline] +pub fn zou_he_outlet(f: &mut [R; Q], rho_out: R, uy_out: R) { + let ux = (f[0] + f[2] + f[4] + 2.0 * (f[1] + f[5] + f[8])) / rho_out - 1.0; + let d = 0.5 * (f[2] - f[4]); + f[3] = f[1] - (2.0 / 3.0) * rho_out * ux; + f[6] = f[8] - d - (1.0 / 6.0) * rho_out * ux + 0.5 * rho_out * uy_out; + f[7] = f[5] + d - (1.0 / 6.0) * rho_out * ux - 0.5 * rho_out * uy_out; +} + +// ───────────────────────────────────────────────────────────────────────────── +// Сила на теле +// ───────────────────────────────────────────────────────────────────────────── + +/// Вклад одного линка в силу по GMEM (галилей-инвариантный обмен импульсом, Wen et al. 2014): +/// +/// ΔF = (c_i − u_w)·f_i^{после столкновения} − (c_ī − u_w)·f_ī^{после стриминга}, c_ī = −c_i +/// +/// f_ī берётся из поля ПОСЛЕ стриминга — иначе ведущий симметричный член ~2ρW_i сокращается +/// и Cd систематически врёт. При u_w = 0 сводится к классическому Σ c_i(f_i + f_ī). +#[inline] +pub fn gmem_link(i: usize, f_post: R, f_back: R, uwx: R, uwy: R) -> (R, R) { + let cx = CX[i] as R; + let cy = CY[i] as R; + ((cx - uwx) * f_post + (cx + uwx) * f_back, (cy - uwy) * f_post + (cy + uwy) * f_back) +} + +/// Безразмерные коэффициенты. `d_ref` — характерный размер В ТЕХ ЖЕ единицах, что сила +/// (т.е. на том уровне сетки, где сила снималась). +#[inline] +pub fn coefficients(fx: R, fy: R, tz: R, u: R, d_ref: R) -> (R, R, R) { + let s = u * u * d_ref; + (2.0 * fx / s, 2.0 * fy / s, 2.0 * tz / (s * d_ref)) +} + +// ───────────────────────────────────────────────────────────────────────────── +// Спектральная диагностика (число Струхаля) +// ───────────────────────────────────────────────────────────────────────────── + +/// Комплексное БПФ по основанию 2, на месте. n обязан быть степенью двойки. +fn fft_pow2(re: &mut [R], im: &mut [R]) { + let n = re.len(); + debug_assert!(n.is_power_of_two() && im.len() == n); + // бит-реверс + let mut j = 0usize; + for i in 1..n { + let mut bit = n >> 1; + while j & bit != 0 { + j ^= bit; + bit >>= 1; + } + j |= bit; + if i < j { + re.swap(i, j); + im.swap(i, j); + } + } + let mut len = 2usize; + while len <= n { + let ang = -2.0 * PI / len as R; + let (ws, wc) = ang.sin_cos(); + let mut i = 0usize; + while i < n { + let (mut cr, mut ci) = (1.0, 0.0); + for k in 0..len / 2 { + let ur = re[i + k]; + let ui = im[i + k]; + let vr = re[i + k + len / 2] * cr - im[i + k + len / 2] * ci; + let vi = re[i + k + len / 2] * ci + im[i + k + len / 2] * cr; + re[i + k] = ur + vr; + im[i + k] = ui + vi; + re[i + k + len / 2] = ur - vr; + im[i + k + len / 2] = ui - vi; + let nr = cr * wc - ci * ws; + ci = cr * ws + ci * wc; + cr = nr; + let _ = k; + } + i += len; + } + len <<= 1; + } +} + +/// Число Струхаля St = f·D/U по ряду поперечной скорости в следе. +/// +/// Берётся вторая половина ряда (установившийся режим), убирается среднее, спектр считается +/// с нулевым дополнением до `pad`, вершина пика уточняется параболической интерполяцией по +/// трём точкам — иначе St «залипает» на сетке БПФ-бинов (шаг ≈ D/(U·pad)). +/// Возвращает (St, амплитудный спектр, шаг частоты). +pub fn strouhal(signal: &[R], d_lat: R, u_lat: R, pad: usize) -> (R, Vec, R) { + let pad = pad.next_power_of_two(); + if signal.len() < 8 { + return (0.0, Vec::new(), 0.0); + } + let tail = &signal[signal.len() / 2..]; + let mean = tail.iter().sum::() / tail.len() as R; + let mut re = vec![0.0; pad]; + let mut im = vec![0.0; pad]; + for (k, &v) in tail.iter().enumerate().take(pad) { + re[k] = v - mean; + } + fft_pow2(&mut re, &mut im); + let half = pad / 2; + let amp: Vec = (0..half).map(|k| (re[k] * re[k] + im[k] * im[k]).sqrt()).collect(); + if amp.len() <= 3 { + return (0.0, amp, 1.0 / pad as R); + } + let mut k = 1; + for j in 2..amp.len() - 1 { + if amp[j] > amp[k] { + k = j; + } + } + let delta = { + let (a0, a1, a2) = (amp[k - 1], amp[k], amp[k + 1]); + let den = a0 - 2.0 * a1 + a2; + if den != 0.0 { + (0.5 * (a0 - a2) / den).clamp(-0.5, 0.5) + } else { + 0.0 + } + }; + let f_peak = (k as R + delta) / pad as R; // циклов за шаг + (f_peak * d_lat / u_lat, amp, 1.0 / pad as R) +} + +/// Среднее и стандартное отклонение среза. +pub fn mean_std(v: &[R]) -> (R, R) { + if v.is_empty() { + return (0.0, 0.0); + } + let m = v.iter().sum::() / v.len() as R; + let var = v.iter().map(|x| (x - m) * (x - m)).sum::() / v.len() as R; + (m, var.sqrt()) +} + +/// smoothstep-разгон: гасит импульсный старт (иначе ударная волна от мгновенного включения входа). +#[inline] +pub fn smoothstep(t: R) -> R { + let t = t.clamp(0.0, 1.0); + t * t * (3.0 - 2.0 * t) +} + +// ───────────────────────────────────────────────────────────────────────────── +// Тесты: сверка с формулами статей +// ───────────────────────────────────────────────────────────────────────────── + +#[cfg(test)] +mod tests { + use super::*; + + /// Проектор Δs, построенный аналитически, обязан совпадать с матричным + /// Ps = M⁻¹·diag(…,1,1,…)·M, где M — базис натуральных моментов из (6)–(7). + #[test] + fn projector_matches_moment_matrix() { + // строки M: 1, cx, cy, 3(cx²+cy²)−2, cx²−cy², cx·cy, cx²cy, cx cy², cx²cy² + let mut m = [[0.0f64; Q]; Q]; + for i in 0..Q { + let (x, y) = (CX[i] as f64, CY[i] as f64); + m[0][i] = 1.0; + m[1][i] = x; + m[2][i] = y; + m[3][i] = 3.0 * (x * x + y * y) - 2.0; + m[4][i] = x * x - y * y; + m[5][i] = x * y; + m[6][i] = x * x * y; + m[7][i] = x * y * y; + m[8][i] = x * x * y * y; + } + let minv = invert9(&m); + // Ps = M⁻¹ · D · M, D оставляет строки 4 (N) и 5 (Π_xy) + let mut ps = [[0.0f64; Q]; Q]; + for a in 0..Q { + for b in 0..Q { + ps[a][b] = minv[a][4] * m[4][b] + minv[a][5] * m[5][b]; + } + } + // сверяем на случайных (детерминированных) векторах + let mut seed = 12345u64; + let mut rnd = || { + seed = seed.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + ((seed >> 33) as f64 / (1u64 << 31) as f64) - 1.0 + }; + for _ in 0..200 { + let d: [f64; Q] = std::array::from_fn(|_| rnd()); + let got = project_shift(&d); + for a in 0..Q { + let want: f64 = (0..Q).map(|b| ps[a][b] * d[b]).sum(); + assert!((got[a] - want).abs() < 1e-13, "i={a}: {} vs {}", got[a], want); + } + } + } + + /// Ps обязан быть идемпотентным (это проектор) и не трогать сохраняющиеся моменты. + #[test] + fn projector_is_idempotent_and_conservative() { + let mut seed = 999u64; + let mut rnd = || { + seed = seed.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + ((seed >> 33) as f64 / (1u64 << 31) as f64) - 1.0 + }; + for _ in 0..100 { + let d: [f64; Q] = std::array::from_fn(|_| rnd()); + let s1 = project_shift(&d); + let s2 = project_shift(&s1); + for i in 0..Q { + assert!((s1[i] - s2[i]).abs() < 1e-14); + } + // Δs не несёт массы и импульса + assert!(s1.iter().sum::().abs() < 1e-13); + let mx: f64 = (0..Q).map(|i| CX[i] as f64 * s1[i]).sum(); + let my: f64 = (0..Q).map(|i| CY[i] as f64 * s1[i]).sum(); + assert!(mx.abs() < 1e-13 && my.abs() < 1e-13); + } + } + + /// Равновесие обязано точно сохранять ρ и ρu (оно — максимизатор энтропии при этих связях). + #[test] + fn equilibrium_conserves_moments() { + for &(rho, ux, uy) in &[(1.0, 0.0, 0.0), (1.2, 0.08, -0.03), (0.9, -0.11, 0.07)] { + let fe = feq(rho, ux, uy); + let (r, x, y) = macros(&fe); + assert!((r - rho).abs() < 1e-13, "ρ: {r} vs {rho}"); + assert!((x - ux).abs() < 1e-13, "ux: {x} vs {ux}"); + assert!((y - uy).abs() < 1e-13, "uy: {y} vs {uy}"); + } + } + + /// γ из замкнутой формулы (17) обязана быть корнем условия критической точки (15) + /// с точностью до порядка разложения. Проверяем, что невязка (15) мала по сравнению + /// с её же масштабом при γ, отличающейся на единицу. + #[test] + fn gamma_is_root_of_entropy_condition() { + let beta = 1.0 / (2.0 * 0.6); + let f = perturbed_state(0.9); + let (rho, ux, uy) = macros(&f); + let fe = feq(rho, ux, uy); + let d: [f64; Q] = std::array::from_fn(|i| f[i] - fe[i]); + let ds = project_shift(&d); + let dh: [f64; Q] = std::array::from_fn(|i| d[i] - ds[i]); + + // невязка условия (15) + let resid = |g: f64| -> f64 { + (0..Q) + .map(|i| { + dh[i] * (1.0 + ((1.0 - beta * g) * dh[i] - (2.0 * beta - 1.0) * ds[i]) / fe[i]).ln() + }) + .sum() + }; + let mut fc = f; + let info = collide_node(&mut fc, beta); + let r0 = resid(info.gamma).abs(); + let r1 = resid(info.gamma + 1.0).abs().min(resid(info.gamma - 1.0).abs()); + assert!(r0 < 0.05 * r1, "невязка (15) при γ*: {r0:.3e}, на ±1 от него: {r1:.3e}"); + } + + /// При γ = 2 KBC обязан совпасть с LBGK поточечно. + #[test] + fn kbc_reduces_to_bgk_at_gamma_two() { + let beta = 0.7; + let f = perturbed_state(0.5); + // конструируем состояние с Δh ≡ 0: возмущаем только сдвиговые моменты + let (rho, ux, uy) = macros(&f); + let fe = feq(rho, ux, uy); + let d: [f64; Q] = std::array::from_fn(|i| f[i] - fe[i]); + let ds = project_shift(&d); + let mut pure_shear: [f64; Q] = std::array::from_fn(|i| fe[i] + ds[i]); + let mut bgk = pure_shear; + collide_node(&mut pure_shear, beta); + collide_node_bgk(&mut bgk, beta); + for i in 0..Q { + assert!((pure_shear[i] - bgk[i]).abs() < 1e-12, "i={i}"); + } + } + + /// Сдвиговые моменты обязаны релаксировать ровно с 2β при ЛЮБОЙ γ — именно это гарантирует, + /// что вязкость (5) задаётся только β, а стабилизатор её не трогает. + #[test] + fn shear_relaxes_at_two_beta_regardless_of_gamma() { + for &beta in &[0.3, 0.6, 0.95] { + let mut f = perturbed_state(0.8); + let (rho, ux, uy) = macros(&f); + let fe = feq(rho, ux, uy); + let d0: [f64; Q] = std::array::from_fn(|i| f[i] - fe[i]); + let s0 = project_shift(&d0); + collide_node(&mut f, beta); + let (r1, x1, y1) = macros(&f); + let fe1 = feq(r1, x1, y1); + let d1: [f64; Q] = std::array::from_fn(|i| f[i] - fe1[i]); + let s1 = project_shift(&d1); + for i in 0..Q { + let want = (1.0 - 2.0 * beta) * s0[i]; + assert!((s1[i] - want).abs() < 1e-10, "β={beta} i={i}: {} vs {want}", s1[i]); + } + } + } + + /// Столкновение обязано сохранять ρ и ρu (иначе оператор не консервативен). + #[test] + fn collision_conserves_mass_and_momentum() { + for &beta in &[0.25, 0.5, 0.99] { + let mut f = perturbed_state(0.7); + let before = macros(&f); + collide_node(&mut f, beta); + let after = macros(&f); + assert!((before.0 - after.0).abs() < 1e-13); + assert!((before.1 - after.1).abs() < 1e-12); + assert!((before.2 - after.2).abs() < 1e-12); + } + } + + /// Zou–He вход обязан ставить ровно заданную скорость. + #[test] + fn zou_he_inlet_sets_velocity() { + let mut f = perturbed_state(0.3); + zou_he_inlet(&mut f, 0.05, -0.01); + let (_, ux, uy) = macros(&f); + assert!((ux - 0.05).abs() < 1e-12, "ux={ux}"); + assert!((uy + 0.01).abs() < 1e-12, "uy={uy}"); + } + + /// Zou–He выход обязан ставить ровно заданную плотность. + #[test] + fn zou_he_outlet_sets_density() { + let mut f = perturbed_state(0.3); + zou_he_outlet(&mut f, 1.0, 0.0); + let (rho, _, _) = macros(&f); + assert!((rho - 1.0).abs() < 1e-12, "ρ={rho}"); + } + + /// SDF: знак внутри/снаружи и |∇φ| ≈ 1 у поверхности для каждой формы. + #[test] + fn sdf_sign_and_gradient() { + for kind in [ + ShapeKind::Cylinder, + ShapeKind::Square, + ShapeKind::Diamond, + ShapeKind::Ellipse, + ShapeKind::Naca, + ShapeKind::Triangle, + ShapeKind::Plate, + ] { + let b = Body::new(kind, 50.0, 50.0, 20.0, 15.0, 0.3); + assert!(b.sdf(50.0, 50.0) < 0.0, "{kind:?}: центр обязан быть внутри"); + assert!(b.sdf(200.0, 200.0) > 0.0, "{kind:?}: далёкая точка обязана быть снаружи"); + // градиент на кольце вокруг тела + let h = 1e-4; + for k in 0..24 { + let a = 2.0 * PI * k as f64 / 24.0; + let (x, y) = (50.0 + 25.0 * a.cos(), 50.0 + 25.0 * a.sin()); + let gx = (b.sdf(x + h, y) - b.sdf(x - h, y)) / (2.0 * h); + let gy = (b.sdf(x, y + h) - b.sdf(x, y - h)) / (2.0 * h); + let g = (gx * gx + gy * gy).sqrt(); + assert!((g - 1.0).abs() < 0.05, "{kind:?}: |∇φ|={g} при угле {a}"); + } + } + } + + /// Единицы: δt обязан удовлетворять δt = u_lat·δx/u_phys, а формула из ТЗ + /// (300 шагов/с при 30 м/с и 0.1 м) обязана воспроизводиться при u_lat = 1. + #[test] + fn units_time_scaling() { + let u = Units::new(0.1, 30.0, 1.0, 150.0, 16.0); + assert!((u.steps_per_second() - 300.0).abs() < 1e-9, "{}", u.steps_per_second()); + let u2 = Units::new(0.1, 30.0, 0.05, 150.0, 16.0); + assert!((u2.steps_per_second() - 6000.0).abs() < 1e-9); + // ν_phys = U·D/Re независимо от выбора u_lat + let want = 30.0 * (16.0 * 0.1) / 150.0; + assert!((u.nu_phys() - want).abs() < 1e-12); + assert!((u2.nu_phys() - want).abs() < 1e-12); + } + + /// БПФ: чистая синусоида обязана давать пик ровно на своей частоте. + #[test] + fn strouhal_recovers_known_frequency() { + let f0 = 0.017; // циклов за шаг + let sig: Vec = (0..4096).map(|k| (2.0 * PI * f0 * k as f64).sin()).collect(); + let (st, _, _) = strouhal(&sig, 16.0, 0.07, 8192); + let want = f0 * 16.0 / 0.07; + assert!((st - want).abs() / want < 2e-3, "St={st} vs {want}"); + } + + // ── вспомогательное ── + + /// Состояние с заметным неравновесием (масштаб `amp` от равновесия) — детерминированное. + fn perturbed_state(amp: f64) -> [f64; Q] { + let fe = feq(1.05, 0.06, -0.02); + let pat = [0.3, -0.7, 0.5, 0.2, -0.4, 0.9, -0.6, 0.15, -0.35]; + std::array::from_fn(|i| fe[i] * (1.0 + amp * 0.1 * pat[i])) + } + + /// Обращение матрицы 9×9 методом Гаусса–Жордана (только для теста). + fn invert9(a: &[[f64; Q]; Q]) -> [[f64; Q]; Q] { + let mut m = *a; + let mut inv = [[0.0f64; Q]; Q]; + for i in 0..Q { + inv[i][i] = 1.0; + } + for col in 0..Q { + let mut piv = col; + for r in col + 1..Q { + if m[r][col].abs() > m[piv][col].abs() { + piv = r; + } + } + m.swap(col, piv); + inv.swap(col, piv); + let d = m[col][col]; + assert!(d.abs() > 1e-12, "вырожденная матрица моментов"); + for c in 0..Q { + m[col][c] /= d; + inv[col][c] /= d; + } + for r in 0..Q { + if r != col { + let k = m[r][col]; + for c in 0..Q { + m[r][c] -= k * m[col][c]; + inv[r][c] -= k * inv[col][c]; + } + } + } + } + inv + } +}