From 7beded2ad02cd22d2b3ae5ac8c73b9d6da18ca48 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Simonis?= Date: Fri, 30 Sep 2022 18:31:19 +0200 Subject: [PATCH 01/11] Add elastic-tube-1d rust solvers --- elastic-tube-1d/fluid-rust/.gitignore | 1 + elastic-tube-1d/fluid-rust/Cargo.lock | 421 +++++++++++++++++++++++++ elastic-tube-1d/fluid-rust/Cargo.toml | 12 + elastic-tube-1d/fluid-rust/clean.sh | 7 + elastic-tube-1d/fluid-rust/run.sh | 4 + elastic-tube-1d/fluid-rust/src/main.rs | 403 +++++++++++++++++++++++ elastic-tube-1d/solid-rust/.gitignore | 1 + elastic-tube-1d/solid-rust/Cargo.lock | 274 ++++++++++++++++ elastic-tube-1d/solid-rust/Cargo.toml | 11 + elastic-tube-1d/solid-rust/clean.sh | 6 + elastic-tube-1d/solid-rust/run.sh | 4 + elastic-tube-1d/solid-rust/src/main.rs | 114 +++++++ 12 files changed, 1258 insertions(+) create mode 100644 elastic-tube-1d/fluid-rust/.gitignore create mode 100644 elastic-tube-1d/fluid-rust/Cargo.lock create mode 100644 elastic-tube-1d/fluid-rust/Cargo.toml create mode 100755 elastic-tube-1d/fluid-rust/clean.sh create mode 100755 elastic-tube-1d/fluid-rust/run.sh create mode 100644 elastic-tube-1d/fluid-rust/src/main.rs create mode 100644 elastic-tube-1d/solid-rust/.gitignore create mode 100644 elastic-tube-1d/solid-rust/Cargo.lock create mode 100644 elastic-tube-1d/solid-rust/Cargo.toml create mode 100755 elastic-tube-1d/solid-rust/clean.sh create mode 100755 elastic-tube-1d/solid-rust/run.sh create mode 100644 elastic-tube-1d/solid-rust/src/main.rs diff --git a/elastic-tube-1d/fluid-rust/.gitignore b/elastic-tube-1d/fluid-rust/.gitignore new file mode 100644 index 000000000..2f7896d1d --- /dev/null +++ b/elastic-tube-1d/fluid-rust/.gitignore @@ -0,0 +1 @@ +target/ diff --git a/elastic-tube-1d/fluid-rust/Cargo.lock b/elastic-tube-1d/fluid-rust/Cargo.lock new file mode 100644 index 000000000..a62a2e950 --- /dev/null +++ b/elastic-tube-1d/fluid-rust/Cargo.lock @@ -0,0 +1,421 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 3 + +[[package]] +name = "ansi_term" +version = "0.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d52a9bb7ec0cf484c551830a7ce27bd20d67eac647e1befb56b0be4ee39a55d2" +dependencies = [ + "winapi", +] + +[[package]] +name = "approx" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" +dependencies = [ + "num-traits", +] + +[[package]] +name = "atty" +version = "0.2.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d9b39be18770d11421cdb1b9947a45dd3f37e93092cbf377614828a319d5fee8" +dependencies = [ + "hermit-abi", + "libc", + "winapi", +] + +[[package]] +name = "autocfg" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa" + +[[package]] +name = "bitflags" +version = "1.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" + +[[package]] +name = "bytemuck" +version = "1.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2f5715e491b5a1598fc2bef5a606847b5dc1d48ea625bd3c02c00de8285591da" + +[[package]] +name = "cc" +version = "1.0.73" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2fff2a6927b3bb87f9595d67196a70493f627687a71d87a0d692242c33f58c11" + +[[package]] +name = "clap" +version = "2.34.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a0610544180c38b88101fecf2dd634b174a62eef6946f84dfc6a7127512b381c" +dependencies = [ + "ansi_term", + "atty", + "bitflags", + "strsim", + "textwrap", + "unicode-width", + "vec_map", +] + +[[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 = "cxx" +version = "1.0.78" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "19f39818dcfc97d45b03953c1292efc4e80954e1583c4aa770bac1383e2310a4" +dependencies = [ + "cc", + "cxxbridge-flags", + "cxxbridge-macro", + "link-cplusplus", +] + +[[package]] +name = "cxx-build" +version = "1.0.78" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3e580d70777c116df50c390d1211993f62d40302881e54d4b79727acb83d0199" +dependencies = [ + "cc", + "codespan-reporting", + "once_cell", + "proc-macro2", + "quote", + "scratch", + "syn", +] + +[[package]] +name = "cxxbridge-flags" +version = "1.0.78" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "56a46460b88d1cec95112c8c363f0e2c39afdb237f60583b0b36343bf627ea9c" + +[[package]] +name = "cxxbridge-macro" +version = "1.0.78" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "747b608fecf06b0d72d440f27acc99288207324b793be2c17991839f3d4995ea" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "fluid-rust" +version = "0.1.0" +dependencies = [ + "nalgebra", + "precice", +] + +[[package]] +name = "hermit-abi" +version = "0.1.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "62b467343b94ba476dcb2500d242dadbb39557df889310ac77c5d99100aaac33" +dependencies = [ + "libc", +] + +[[package]] +name = "libc" +version = "0.2.139" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "201de327520df007757c1f0adce6e827fe8562fbc28bfd9c15571c66ca1f5f79" + +[[package]] +name = "link-cplusplus" +version = "1.0.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9272ab7b96c9046fbc5bc56c06c117cb639fe2d509df0c421cad82d2915cf369" +dependencies = [ + "cc", +] + +[[package]] +name = "matrixmultiply" +version = "0.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "add85d4dd35074e6fedc608f8c8f513a3548619a9024b751949ef0e8e45a4d84" +dependencies = [ + "rawpointer", +] + +[[package]] +name = "nalgebra" +version = "0.31.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e9e0a04ce089f9401aac565c740ed30c46291260f27d4911fdbaa6ca65fa3044" +dependencies = [ + "approx", + "matrixmultiply", + "nalgebra-macros", + "num-complex", + "num-rational", + "num-traits", + "simba", + "typenum", +] + +[[package]] +name = "nalgebra-macros" +version = "0.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "01fcc0b8149b4632adc89ac3b7b31a12fb6099a0317a4eb2ebff574ef7de7218" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "num-complex" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7ae39348c8bc5fbd7f40c727a9925f03517afd2ab27d46702108b6a7e5414c19" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-integer" +version = "0.1.45" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "225d3389fb3509a24c93f5c29eb6bde2586b98d9f016636dff58d7c6f7569cd9" +dependencies = [ + "autocfg", + "num-traits", +] + +[[package]] +name = "num-rational" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0638a1c9d0a3c0914158145bc76cff373a75a627e6ecbfb71cbe6f453a5a19b0" +dependencies = [ + "autocfg", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "578ede34cf02f8924ab9447f50c28075b4d3e5b269972345e7e0372b38c6cdcd" +dependencies = [ + "autocfg", +] + +[[package]] +name = "once_cell" +version = "1.15.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e82dad04139b71a90c080c8463fe0dc7902db5192d939bd0950f074d014339e1" + +[[package]] +name = "paste" +version = "1.0.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b1de2e551fb905ac83f73f7aedf2f0cb4a0da7e35efa24a202a936269f1f18e1" + +[[package]] +name = "pkg-config" +version = "0.3.25" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1df8c4ec4b0627e53bdf214615ad287367e482558cf84b109250b37464dc03ae" + +[[package]] +name = "precice" +version = "2.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dec31fcf7edff2ec658584c7010a28da5f1e3d87b42b4f08a17df79d63305d91" +dependencies = [ + "clap", + "cxx", + "cxx-build", + "pkg-config", + "semver", +] + +[[package]] +name = "proc-macro2" +version = "1.0.46" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "94e2ef8dbfc347b10c094890f778ee2e36ca9bb4262e86dc99cd217e35f3470b" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "quote" +version = "1.0.21" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bbe448f377a7d6961e30f5955f9b8d106c3f5e449d493ee1b125c1d43c2b5179" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "rawpointer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" + +[[package]] +name = "safe_arch" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "794821e4ccb0d9f979512f9c1973480123f9bd62a90d74ab0f9426fcf8f4a529" +dependencies = [ + "bytemuck", +] + +[[package]] +name = "scratch" +version = "1.0.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9c8132065adcfd6e02db789d9285a0deb2f3fcb04002865ab67d5fb103533898" + +[[package]] +name = "semver" +version = "1.0.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "58bc9567378fc7690d6b2addae4e60ac2eeea07becb2c64b9f218b53865cba2a" + +[[package]] +name = "simba" +version = "0.7.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c48e45e5961033db030b56ad67aef22e9c908c493a6e8348c0a0f6b93433cd77" +dependencies = [ + "approx", + "num-complex", + "num-traits", + "paste", + "wide", +] + +[[package]] +name = "strsim" +version = "0.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8ea5119cdb4c55b55d432abb513a0429384878c15dde60cc77b1c99de1a95a6a" + +[[package]] +name = "syn" +version = "1.0.101" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e90cde112c4b9690b8cbe810cba9ddd8bc1d7472e2cae317b69e9438c1cba7d2" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "termcolor" +version = "1.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bab24d30b911b2376f3a13cc2cd443142f0c81dda04c118693e35b3835757755" +dependencies = [ + "winapi-util", +] + +[[package]] +name = "textwrap" +version = "0.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d326610f408c7a4eb6f51c37c330e496b08506c9457c9d34287ecc38809fb060" +dependencies = [ + "unicode-width", +] + +[[package]] +name = "typenum" +version = "1.15.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dcf81ac59edc17cc8697ff311e8f5ef2d99fcbd9817b34cec66f90b6c3dfd987" + +[[package]] +name = "unicode-ident" +version = "1.0.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dcc811dc4066ac62f84f11307873c4850cb653bfa9b1719cee2bd2204a4bc5dd" + +[[package]] +name = "unicode-width" +version = "0.1.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c0edd1e5b14653f783770bce4a4dabb4a5108a5370a5f5d8cfe8710c361f6c8b" + +[[package]] +name = "vec_map" +version = "0.8.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f1bddf1187be692e79c5ffeab891132dfb0f236ed36a43c7ed39f1165ee20191" + +[[package]] +name = "wide" +version = "0.7.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b3aba2d1dac31ac7cae82847ac5b8be822aee8f99a4e100f279605016b185c5f" +dependencies = [ + "bytemuck", + "safe_arch", +] + +[[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.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "70ec6ce85bb158151cae5e5c87f95a8e97d2c0c4b001223f33a334e3ce5de178" +dependencies = [ + "winapi", +] + +[[package]] +name = "winapi-x86_64-pc-windows-gnu" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" diff --git a/elastic-tube-1d/fluid-rust/Cargo.toml b/elastic-tube-1d/fluid-rust/Cargo.toml new file mode 100644 index 000000000..c490d7a09 --- /dev/null +++ b/elastic-tube-1d/fluid-rust/Cargo.toml @@ -0,0 +1,12 @@ +[package] +name = "fluid-rust" +version = "0.1.0" +edition = "2021" +publish = false +repository = "https://github.com/precice/tutorials" +description = "The fluid participant of the elastic-tube 1D tutorial" +license = "LGPL-3.0-or-later" + +[dependencies] +"precice" = "2.5" +"nalgebra" = "0.31" diff --git a/elastic-tube-1d/fluid-rust/clean.sh b/elastic-tube-1d/fluid-rust/clean.sh new file mode 100755 index 000000000..ff2895ffa --- /dev/null +++ b/elastic-tube-1d/fluid-rust/clean.sh @@ -0,0 +1,7 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +rm -rvf ./output/*.vtk +clean_precice_logs . diff --git a/elastic-tube-1d/fluid-rust/run.sh b/elastic-tube-1d/fluid-rust/run.sh new file mode 100755 index 000000000..d9103ab27 --- /dev/null +++ b/elastic-tube-1d/fluid-rust/run.sh @@ -0,0 +1,4 @@ +#!/bin/sh +set -e -u + +cargo run --release ../precice-config.xml diff --git a/elastic-tube-1d/fluid-rust/src/main.rs b/elastic-tube-1d/fluid-rust/src/main.rs new file mode 100644 index 000000000..254d55f6d --- /dev/null +++ b/elastic-tube-1d/fluid-rust/src/main.rs @@ -0,0 +1,403 @@ +use nalgebra as na; +use std::env; +use std::process::ExitCode; + +use std::fs; +use std::fs::File; +use std::io::{BufWriter, Write}; +use std::path; + +fn write_vtk( + file: &str, + dx: f64, + velocity: &[f64], + pressure: &[f64], + cross_section_length: &[f64], +) -> std::io::Result<()> { + let filepath = path::Path::new(file); + if !filepath.parent().unwrap().exists() { + fs::create_dir(filepath.parent().unwrap())?; + } + + // let mut f = File::create(file).expect("Unable to open vtk file"); + let file = File::create(file).expect("Unable to open vtk file"); + let mut f = BufWriter::new(file); + + let n_points = velocity.len(); + + write!(f, "# vtk DataFile Version 2.0\n\n")?; + write!(f, "ASCII\n\n")?; + write!(f, "DATASET UNSTRUCTURED_GRID\n\n")?; + write!(f, "POINTS {n_points} float\n\n")?; + for i in 0..n_points { + write!( + f, + "{:.16e} 0.0000000000000000e+00 0.0000000000000000e+00\n", + i as f64 * dx + )?; + } + writeln!(f)?; + + write!(f, "POINT_DATA {n_points}\n\n")?; + + write!(f, "VECTORS velocity float\n")?; + for v in velocity { + write!( + f, + "{:.16e} 0.0000000000000000e+00 0.0000000000000000e+00\n", + v + )?; + } + writeln!(f)?; + + write!(f, "SCALARS pressure float\n")?; + write!(f, "LOOKUP_TABLE default\n")?; + for v in pressure { + write!(f, "{:.16e}\n", v)?; + } + writeln!(f)?; + + write!(f, "SCALARS diameter float\n")?; + write!(f, "LOOKUP_TABLE default\n")?; + for v in cross_section_length { + write!(f, "{:.16e}\n", v)?; + } + writeln!(f)?; + Ok(()) +} + +fn fluid_compute_solution( + velocity_old: &[f64], + pressure_old: &[f64], + cross_section_length_old: &[f64], + cross_section_length: &[f64], + t: f64, + n: usize, + kappa: f64, + tau: f64, + velocity: &mut [f64], + pressure: &mut [f64], +) { + // Initial guess + pressure.copy_from_slice(&pressure_old[..]); + + const E: f64 = 10000_f64; + //const C_MK2 : f64 = E / std::f64::consts::FRAC_2_SQRT_PI; + let c_mk2: f64 = E / 2.0 * std::f64::consts::PI.sqrt(); + + // lhs = 2*N+2 + + const ALPHA: f64 = 0.0; + const L: f64 = 10.0; + let dx = L / kappa; + + let equations = 2 * n + 2; + + let mut k = 0; + loop { + k += 1; + + let res = { + let mut res = na::DVector::from_element(equations, 0.0); + for i in 1..n { + let [p_l, p_c, p_r]: [f64; 3] = pressure[i - 1..i + 2].try_into().unwrap(); + let [v_l, v_c, v_r]: [f64; 3] = velocity[i - 1..i + 2].try_into().unwrap(); + let [c_l, c_c, c_r]: [f64; 3] = + cross_section_length[i - 1..i + 2].try_into().unwrap(); + + /* Momentum */ + res[i] = (velocity_old[i] * cross_section_length_old[i] - v_c * c_c) * dx / tau + + 0.25 * (-c_r * v_c * v_r - c_c * v_c * v_r) + + 0.25 + * (-c_r * v_c * v_c - c_c * v_c * v_c + c_c * v_l * v_c + c_l * v_l * v_c) + + 0.25 * (c_l * v_l * v_l + c_c * v_l * v_l) + + 0.25 + * (c_l * p_l + c_c * p_l - c_l * p_c + c_r * p_c - c_c * p_r - c_r * p_r); + + /* Continuity */ + res[i + n + 1] = (cross_section_length_old[i] - c_c) * dx / tau + + 0.25 + * (c_l * v_l + c_c * v_l + c_l * v_c - c_r * v_c - c_c * v_r - c_r * v_r) + + ALPHA * (p_l - 2.0 * p_c + p_r); + } + /* Boundary */ + + /* Velocity Inlet is prescribed */ + let u0 = 10.0; + let ampl = 3.0; + let frequency = 10.0; + let t_shift = 0.0; + let velocity_in = + u0 + ampl * (frequency * (t + tau + t_shift) * std::f64::consts::PI).sin(); + res[0] = velocity_in - velocity[0]; + + /* Pressure Inlet is linearly interpolated */ + res[n + 1] = -pressure[0] + 2.0 * pressure[1] - pressure[2]; + + /* Velocity Outlet is linearly interpolated */ + res[n] = -velocity[n] + 2.0 * velocity[n - 1] - velocity[n - 2]; + + /* Pressure Outlet is "non-reflecting" */ + let tmp2 = + (c_mk2 - pressure_old[n] / 2.0).sqrt() - (velocity[n] - velocity_old[n]) / 4.0; + res[2 * n + 1] = -pressure[n] + 2.0 * (c_mk2 - tmp2.powi(2)); + res + }; + + // compute norm of residual + let norm: f64 = { + let norm1 = res.map(|x: f64| x * x).sum().sqrt(); + let norm2 = pressure + .iter() + .chain(velocity.iter()) + .map(|x| x * x) + .sum::() + .sqrt(); + norm1 / norm2 + }; + + const TOL: f64 = 1e-10; + const MAX_ITER: usize = 1000; + if norm < TOL && k > 1 { + //println!("Success at k={} res={} tol={}", k, norm, TOL); + break; + } else if k > MAX_ITER { + panic!("Nonlinear Solver break, iterations: {k}, residual norm: {norm}"); + } else { + //println!("Continuing k={} res={:.10} tol={:.10}", k, norm, TOL); + } + + /* Initializing the the LHS i.e. Left Hand Side */ + let lhs = { + let mut lhs = na::DMatrix::from_element(equations, equations, 0.0); + for i in 1..n { + let [v_l, v_c, v_r]: [f64; 3] = velocity[i - 1..i + 2].try_into().unwrap(); + let [c_l, c_c, c_r]: [f64; 3] = + cross_section_length[i - 1..i + 2].try_into().unwrap(); + + // Momentum, Velocity + lhs[(i, i - 1)] = + 0.25 * (-2.0 * c_l * v_l - 2.0 * c_c * v_l - c_c * v_c - c_l * v_c); + + lhs[(i, i)] = c_c * dx / tau + + 0.25 + * (c_r * v_r + c_c * v_r + 2.0 * c_r * v_c + 2.0 * c_c * v_c + - c_c * v_l + - c_l * v_l); + lhs[(i, i + 1)] = 0.25 * (c_r * v_c + c_c * v_c); + + // Momentum, Pressure + lhs[(i, n + 1 + i - 1)] = -0.25 * (c_l + c_c); + lhs[(i, n + 1 + i)] = 0.25 * (c_l - c_r); + lhs[(i, n + 1 + i + 1)] = 0.25 * (c_c + c_r); + + // Continuity, Velocity + lhs[(i + n + 1, i - 1)] = -0.25 * (c_l + c_c); + lhs[(i + n + 1, i)] = -0.25 * (c_l - c_r); + lhs[(i + n + 1, i + 1)] = 0.25 * (c_c + c_r); + + // Continuity, Pressure + lhs[(i + n + 1, n + 1 + i - 1)] = -ALPHA; + lhs[(i + n + 1, n + 1 + i)] = 2.0 * ALPHA; + lhs[(i + n + 1, n + 1 + i + 1)] = -ALPHA; + } + + /* Boundary */ + + // Velocity Inlet is prescribed + lhs[(0, 0)] = 1.0; + // Pressure Inlet is linearly interpolated + lhs[(n + 1, n + 1)] = 1.0; + lhs[(n + 1, n + 2)] = -2.0; + lhs[(n + 1, n + 3)] = 1.0; + // Velocity Outlet is linearly interpolated + lhs[(n, n)] = 1.0; + lhs[(n, n - 1)] = -2.0; + lhs[(n, n - 2)] = 1.0; + // Pressure Outlet is non-Reflecting + lhs[(2 * n + 1, 2 * n + 1)] = 1.0; + lhs[(2 * n + 1, n)] = + -((c_mk2 - pressure_old[n] / 2.0).sqrt() - (velocity[n] - velocity_old[n]) / 4.0); + lhs + }; + + let sol = lhs + .lu() + .solve(&res) + .expect("Linear Solver did not converge!"); + + velocity + .iter_mut() + .zip(sol.as_slice()[..n + 1].iter()) + .for_each(|p| *p.0 += p.1); + pressure + .iter_mut() + .zip(sol.as_slice()[n + 1..].iter()) + .for_each(|p| *p.0 += p.1); + } +} + +fn main() -> ExitCode { + println!("Starting Fluid Solver..."); + + let args: Vec<_> = env::args().collect(); + + if args.len() != 2 { + println!("Fluid: Usage: {} ", args[0]); + return ExitCode::from(1); + } + + let config = &args[1]; + + const DOMAIN_SIZE: usize = 100; + const CHUNK_SIZE: usize = DOMAIN_SIZE + 1; + + let mut interface = precice::new("Fluid", &config, 0, 1); + + println!("preCICE configured..."); + + let dimensions = interface.get_dimensions(); + let mesh_id = interface.get_mesh_id("Fluid-Nodes-Mesh"); + let cross_section_length_id = interface.get_data_id("CrossSectionLength", mesh_id); + let pressure_id = interface.get_data_id("Pressure", mesh_id); + + const KAPPA: f64 = 100_f64; + const L: f64 = 10_f64; + const DX: f64 = L / KAPPA; + + const R0: f64 = 2.0 * std::f64::consts::FRAC_2_SQRT_PI; + const A0: f64 = R0 * R0 * std::f64::consts::PI; + const U0: f64 = 10.0; + const AMPL: f64 = 3.0; + const FREQUENCY: f64 = 10.0; + const T_SHIFT: f64 = 0.0; + const P0: f64 = 0.0; + let vel_in0: f64 = U0 + AMPL * (FREQUENCY * T_SHIFT * std::f64::consts::PI).sin(); + + let mut pressure: Vec = vec![P0; CHUNK_SIZE as usize]; + let mut cross_section_length: Vec = vec![A0; CHUNK_SIZE as usize]; + let mut velocity: Vec = vec![vel_in0; CHUNK_SIZE as usize]; + + let mut pressure_old = pressure.clone(); + + const CELLWIDTH: f64 = L / DOMAIN_SIZE as f64; + let grid_size = CHUNK_SIZE * dimensions as usize; + let grid: Vec = { + let mut v: Vec = vec![0_f64; grid_size]; + for i in 0..CHUNK_SIZE { + let idx = i * dimensions as usize; + v[idx] = i as f64 * CELLWIDTH as f64; + } + v + }; + + let vertex_ids = { + let mut ids = vec![-1; CHUNK_SIZE]; + interface + .pin_mut() + .set_mesh_vertices(mesh_id, &grid[..], &mut ids[..]); + ids + }; + + println!("Initializing preCICE..."); + + let mut t = 0.0; + let dt = interface.pin_mut().initialize(); + + if interface.is_action_required(&precice::action_write_initial_data()) { + interface.pin_mut().write_block_scalar_data( + cross_section_length_id, + &vertex_ids[..], + &pressure[..], + ); + interface + .pin_mut() + .mark_action_fulfilled(&precice::action_write_initial_data()); + } + + interface.pin_mut().initialize_data(); + + interface.read_block_scalar_data( + cross_section_length_id, + &vertex_ids[..], + &mut cross_section_length[..], + ); + + let mut cross_section_length_old = cross_section_length.clone(); + + // initialize such that mass conservation is fulfilled + let mut velocity_old = { + let csl0 = cross_section_length[0]; + cross_section_length + .iter() + .map(|csl| vel_in0 * csl0 / csl) + .collect::>() + }; + + let mut out_counter = 0; + + while interface.is_coupling_ongoing() { + if interface.is_action_required(&precice::action_write_iteration_checkpoint()) { + interface + .pin_mut() + .mark_action_fulfilled(&precice::action_write_iteration_checkpoint()); + } + + fluid_compute_solution( + &velocity_old[..], + &pressure_old[..], + &cross_section_length_old[..], + &cross_section_length[..], + t, + DOMAIN_SIZE, + KAPPA, + dt, + &mut velocity[..], + &mut pressure[..], + ); + + interface + .pin_mut() + .write_block_scalar_data(pressure_id, &vertex_ids[..], &pressure[..]); + + interface.pin_mut().advance(dt); + + interface.read_block_scalar_data( + cross_section_length_id, + &vertex_ids[..], + &mut cross_section_length[..], + ); + + if interface.is_action_required(&precice::action_read_iteration_checkpoint()) { + // i.e. fluid not yet converged + interface + .pin_mut() + .mark_action_fulfilled(&precice::action_read_iteration_checkpoint()); + + //pressure.copy_from_slice(&pressure_old[..]); + //velocity.copy_from_slice(&velocity_old[..]); + } else { + t += dt; + pressure_old.copy_from_slice(&pressure[..]); + velocity_old.copy_from_slice(&velocity[..]); + cross_section_length_old.copy_from_slice(&cross_section_length[..]); + + let filename = format!("./output/out_fluid{out_counter}.vtk"); + println!("writing timestep at t={} to {}", t, filename); + write_vtk( + &filename, + DX, + &velocity_old[..], + &pressure_old[..], + &cross_section_length_old[..], + ) + .expect("Unable to write the vtk file"); + out_counter += 1; + } + } + + println!("Exiting FluidSolver at t={}", t); + interface.pin_mut().finalize(); + + return ExitCode::SUCCESS; +} diff --git a/elastic-tube-1d/solid-rust/.gitignore b/elastic-tube-1d/solid-rust/.gitignore new file mode 100644 index 000000000..2f7896d1d --- /dev/null +++ b/elastic-tube-1d/solid-rust/.gitignore @@ -0,0 +1 @@ +target/ diff --git a/elastic-tube-1d/solid-rust/Cargo.lock b/elastic-tube-1d/solid-rust/Cargo.lock new file mode 100644 index 000000000..064988cbe --- /dev/null +++ b/elastic-tube-1d/solid-rust/Cargo.lock @@ -0,0 +1,274 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 3 + +[[package]] +name = "ansi_term" +version = "0.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d52a9bb7ec0cf484c551830a7ce27bd20d67eac647e1befb56b0be4ee39a55d2" +dependencies = [ + "winapi", +] + +[[package]] +name = "atty" +version = "0.2.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d9b39be18770d11421cdb1b9947a45dd3f37e93092cbf377614828a319d5fee8" +dependencies = [ + "hermit-abi", + "libc", + "winapi", +] + +[[package]] +name = "bitflags" +version = "1.3.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" + +[[package]] +name = "cc" +version = "1.0.79" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "50d30906286121d95be3d479533b458f87493b30a4b5f79a607db8f5d11aa91f" + +[[package]] +name = "clap" +version = "2.34.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a0610544180c38b88101fecf2dd634b174a62eef6946f84dfc6a7127512b381c" +dependencies = [ + "ansi_term", + "atty", + "bitflags", + "strsim", + "textwrap", + "unicode-width", + "vec_map", +] + +[[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 = "cxx" +version = "1.0.91" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "86d3488e7665a7a483b57e25bdd90d0aeb2bc7608c8d0346acf2ad3f1caf1d62" +dependencies = [ + "cc", + "cxxbridge-flags", + "cxxbridge-macro", + "link-cplusplus", +] + +[[package]] +name = "cxx-build" +version = "1.0.91" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "48fcaf066a053a41a81dfb14d57d99738b767febb8b735c3016e469fac5da690" +dependencies = [ + "cc", + "codespan-reporting", + "once_cell", + "proc-macro2", + "quote", + "scratch", + "syn", +] + +[[package]] +name = "cxxbridge-flags" +version = "1.0.91" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a2ef98b8b717a829ca5603af80e1f9e2e48013ab227b68ef37872ef84ee479bf" + +[[package]] +name = "cxxbridge-macro" +version = "1.0.91" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "086c685979a698443656e5cf7856c95c642295a38599f12fb1ff76fb28d19892" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "hermit-abi" +version = "0.1.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "62b467343b94ba476dcb2500d242dadbb39557df889310ac77c5d99100aaac33" +dependencies = [ + "libc", +] + +[[package]] +name = "libc" +version = "0.2.139" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "201de327520df007757c1f0adce6e827fe8562fbc28bfd9c15571c66ca1f5f79" + +[[package]] +name = "link-cplusplus" +version = "1.0.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ecd207c9c713c34f95a097a5b029ac2ce6010530c7b49d7fea24d977dede04f5" +dependencies = [ + "cc", +] + +[[package]] +name = "once_cell" +version = "1.17.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b7e5500299e16ebb147ae15a00a942af264cf3688f47923b8fc2cd5858f23ad3" + +[[package]] +name = "pkg-config" +version = "0.3.26" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6ac9a59f73473f1b8d852421e59e64809f025994837ef743615c6d0c5b305160" + +[[package]] +name = "precice" +version = "2.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "dec31fcf7edff2ec658584c7010a28da5f1e3d87b42b4f08a17df79d63305d91" +dependencies = [ + "clap", + "cxx", + "cxx-build", + "pkg-config", + "semver", +] + +[[package]] +name = "proc-macro2" +version = "1.0.51" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5d727cae5b39d21da60fa540906919ad737832fe0b1c165da3a34d6548c849d6" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "quote" +version = "1.0.23" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8856d8364d252a14d474036ea1358d63c9e6965c8e5c1885c18f73d70bff9c7b" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "scratch" +version = "1.0.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ddccb15bcce173023b3fedd9436f882a0739b8dfb45e4f6b6002bee5929f61b2" + +[[package]] +name = "semver" +version = "1.0.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "58bc9567378fc7690d6b2addae4e60ac2eeea07becb2c64b9f218b53865cba2a" + +[[package]] +name = "solid-rust" +version = "0.1.0" +dependencies = [ + "precice", +] + +[[package]] +name = "strsim" +version = "0.8.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8ea5119cdb4c55b55d432abb513a0429384878c15dde60cc77b1c99de1a95a6a" + +[[package]] +name = "syn" +version = "1.0.107" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1f4064b5b16e03ae50984a5a8ed5d4f8803e6bc1fd170a3cda91a1be4b18e3f5" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "termcolor" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "be55cf8942feac5c765c2c993422806843c9a9a45d4d5c407ad6dd2ea95eb9b6" +dependencies = [ + "winapi-util", +] + +[[package]] +name = "textwrap" +version = "0.11.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d326610f408c7a4eb6f51c37c330e496b08506c9457c9d34287ecc38809fb060" +dependencies = [ + "unicode-width", +] + +[[package]] +name = "unicode-ident" +version = "1.0.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "84a22b9f218b40614adcb3f4ff08b703773ad44fa9423e4e0d346d5db86e4ebc" + +[[package]] +name = "unicode-width" +version = "0.1.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c0edd1e5b14653f783770bce4a4dabb4a5108a5370a5f5d8cfe8710c361f6c8b" + +[[package]] +name = "vec_map" +version = "0.8.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f1bddf1187be692e79c5ffeab891132dfb0f236ed36a43c7ed39f1165ee20191" + +[[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.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "70ec6ce85bb158151cae5e5c87f95a8e97d2c0c4b001223f33a334e3ce5de178" +dependencies = [ + "winapi", +] + +[[package]] +name = "winapi-x86_64-pc-windows-gnu" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" diff --git a/elastic-tube-1d/solid-rust/Cargo.toml b/elastic-tube-1d/solid-rust/Cargo.toml new file mode 100644 index 000000000..e74a65c3c --- /dev/null +++ b/elastic-tube-1d/solid-rust/Cargo.toml @@ -0,0 +1,11 @@ +[package] +name = "solid-rust" +version = "0.1.0" +edition = "2021" +publish = false +repository = "https://github.com/precice/tutorials" +description = "The solid participant of the elastic-tube 1D tutorial" +license = "LGPL-3.0-or-later" + +[dependencies] +"precice" = "2.5" diff --git a/elastic-tube-1d/solid-rust/clean.sh b/elastic-tube-1d/solid-rust/clean.sh new file mode 100755 index 000000000..7330c8aaa --- /dev/null +++ b/elastic-tube-1d/solid-rust/clean.sh @@ -0,0 +1,6 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_precice_logs . diff --git a/elastic-tube-1d/solid-rust/run.sh b/elastic-tube-1d/solid-rust/run.sh new file mode 100755 index 000000000..d9103ab27 --- /dev/null +++ b/elastic-tube-1d/solid-rust/run.sh @@ -0,0 +1,4 @@ +#!/bin/sh +set -e -u + +cargo run --release ../precice-config.xml diff --git a/elastic-tube-1d/solid-rust/src/main.rs b/elastic-tube-1d/solid-rust/src/main.rs new file mode 100644 index 000000000..b3d858e37 --- /dev/null +++ b/elastic-tube-1d/solid-rust/src/main.rs @@ -0,0 +1,114 @@ +use std::env; +use std::process::ExitCode; + +fn solid_compute_solution(pressure: &[f64], cross_section_length: &mut [f64]) { + assert!(pressure.len() == cross_section_length.len()); + const E: f64 = 10000.0; + const C_MK2: f64 = E / std::f64::consts::FRAC_2_SQRT_PI; + const PRESSURE0: f64 = 0.0; + let new: Vec<_> = pressure + .iter() + .map(|p| ((PRESSURE0 - 2.0 * C_MK2) / (p - 2.0 * C_MK2)).powi(2)) + .collect(); + cross_section_length.copy_from_slice(&new[..]); +} + +fn main() -> ExitCode { + println!("Starting Solid Solver..."); + + let args: Vec<_> = env::args().collect(); + + if args.len() != 2 { + println!("Fluid: Usage: {} ", args[0]); + return ExitCode::from(1); + } + + let config = &args[1]; + + const DOMAIN_SIZE: usize = 100; + const CHUNK_SIZE: usize = DOMAIN_SIZE + 1; + + let mut interface = precice::new("Solid", &config, 0, 1); + + println!("preCICE configured..."); + + let dimensions = interface.get_dimensions(); + let mesh_id = interface.get_mesh_id("Solid-Nodes-Mesh"); + let cross_section_length_id = interface.get_data_id("CrossSectionLength", mesh_id); + let pressure_id = interface.get_data_id("Pressure", mesh_id); + + let mut pressure: Vec = vec![0.0; CHUNK_SIZE as usize]; + let mut cross_section_length: Vec = vec![1.0; CHUNK_SIZE as usize]; + + let grid_size = CHUNK_SIZE * dimensions as usize; + let grid: Vec = { + let mut v: Vec = vec![0_f64; grid_size]; + for i in 0..CHUNK_SIZE - 1 { + for j in 0..(dimensions as usize) - 1 { + let idx = i * dimensions as usize + j; + v[idx] = (i * (1 - j)) as f64; + } + } + v + }; + + let vertex_ids = { + let mut ids = vec![-1; CHUNK_SIZE]; + interface + .pin_mut() + .set_mesh_vertices(mesh_id, &grid[..], &mut ids[..]); + ids + }; + + println!("Initializing preCICE..."); + + let mut t = 0.0; + let dt = interface.pin_mut().initialize(); + + if interface.is_action_required(&precice::action_write_initial_data()) { + interface.pin_mut().write_block_scalar_data( + cross_section_length_id, + &vertex_ids[..], + &cross_section_length[..], + ); + interface + .pin_mut() + .mark_action_fulfilled(&precice::action_write_initial_data()); + } + + interface.pin_mut().initialize_data(); + + while interface.is_coupling_ongoing() { + if interface.is_action_required(&precice::action_write_iteration_checkpoint()) { + interface + .pin_mut() + .mark_action_fulfilled(&precice::action_write_iteration_checkpoint()); + } + + interface.read_block_scalar_data(pressure_id, &vertex_ids[..], &mut pressure[..]); + + solid_compute_solution(&pressure, &mut cross_section_length); + + interface.pin_mut().write_block_scalar_data( + cross_section_length_id, + &vertex_ids[..], + &cross_section_length[..], + ); + + interface.pin_mut().advance(dt); + + if interface.is_action_required(&precice::action_read_iteration_checkpoint()) { + // i.e. fluid not yet converged + interface + .pin_mut() + .mark_action_fulfilled(&precice::action_read_iteration_checkpoint()); + } else { + t += dt; + } + } + + println!("Exiting SolidSolver at t={}", t); + interface.pin_mut().finalize(); + + return ExitCode::SUCCESS; +} From 7e7fa26facf7258900a0bfe758b9fd45faee749f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Simonis?= Date: Tue, 12 Dec 2023 17:36:10 +0100 Subject: [PATCH 02/11] Port eleastic-tube-1d --- elastic-tube-1d/fluid-rust/Cargo.toml | 2 +- elastic-tube-1d/fluid-rust/run.sh | 2 +- elastic-tube-1d/fluid-rust/src/main.rs | 67 +++--- elastic-tube-1d/solid-rust/Cargo.lock | 274 ------------------------- elastic-tube-1d/solid-rust/Cargo.toml | 2 +- elastic-tube-1d/solid-rust/src/main.rs | 62 +++--- 6 files changed, 65 insertions(+), 344 deletions(-) delete mode 100644 elastic-tube-1d/solid-rust/Cargo.lock diff --git a/elastic-tube-1d/fluid-rust/Cargo.toml b/elastic-tube-1d/fluid-rust/Cargo.toml index c490d7a09..b7e82a6e5 100644 --- a/elastic-tube-1d/fluid-rust/Cargo.toml +++ b/elastic-tube-1d/fluid-rust/Cargo.toml @@ -8,5 +8,5 @@ description = "The fluid participant of the elastic-tube 1D tutorial" license = "LGPL-3.0-or-later" [dependencies] -"precice" = "2.5" +"precice" = "3" "nalgebra" = "0.31" diff --git a/elastic-tube-1d/fluid-rust/run.sh b/elastic-tube-1d/fluid-rust/run.sh index d9103ab27..79a022fa3 100755 --- a/elastic-tube-1d/fluid-rust/run.sh +++ b/elastic-tube-1d/fluid-rust/run.sh @@ -1,4 +1,4 @@ #!/bin/sh set -e -u -cargo run --release ../precice-config.xml +cargo run --release ../precice-config.xml diff --git a/elastic-tube-1d/fluid-rust/src/main.rs b/elastic-tube-1d/fluid-rust/src/main.rs index 254d55f6d..1fcbb81de 100644 --- a/elastic-tube-1d/fluid-rust/src/main.rs +++ b/elastic-tube-1d/fluid-rust/src/main.rs @@ -252,14 +252,14 @@ fn main() -> ExitCode { const DOMAIN_SIZE: usize = 100; const CHUNK_SIZE: usize = DOMAIN_SIZE + 1; - let mut interface = precice::new("Fluid", &config, 0, 1); + let mut participant = precice::new("Fluid", &config, 0, 1); println!("preCICE configured..."); - let dimensions = interface.get_dimensions(); - let mesh_id = interface.get_mesh_id("Fluid-Nodes-Mesh"); - let cross_section_length_id = interface.get_data_id("CrossSectionLength", mesh_id); - let pressure_id = interface.get_data_id("Pressure", mesh_id); + let mesh_name = "Fluid-Nodes-Mesh"; + let dimensions = participant.get_mesh_dimensions(mesh_name); + assert!(participant.get_data_dimensions(mesh_name, "CrossSectionLength") == 1); + assert!(participant.get_data_dimensions(mesh_name, "Pressure") == 1); const KAPPA: f64 = 100_f64; const L: f64 = 10_f64; @@ -293,33 +293,30 @@ fn main() -> ExitCode { let vertex_ids = { let mut ids = vec![-1; CHUNK_SIZE]; - interface + participant .pin_mut() - .set_mesh_vertices(mesh_id, &grid[..], &mut ids[..]); + .set_mesh_vertices(mesh_name, &grid[..], &mut ids[..]); ids }; println!("Initializing preCICE..."); - let mut t = 0.0; - let dt = interface.pin_mut().initialize(); - - if interface.is_action_required(&precice::action_write_initial_data()) { - interface.pin_mut().write_block_scalar_data( - cross_section_length_id, + if participant.pin_mut().requires_initial_data() { + participant.pin_mut().write_data( + mesh_name, + "Pressure", &vertex_ids[..], &pressure[..], ); - interface - .pin_mut() - .mark_action_fulfilled(&precice::action_write_initial_data()); } - interface.pin_mut().initialize_data(); + participant.pin_mut().initialize(); - interface.read_block_scalar_data( - cross_section_length_id, + participant.read_data( + mesh_name, + "CrossSectionLength", &vertex_ids[..], + 0.0, &mut cross_section_length[..], ); @@ -336,13 +333,15 @@ fn main() -> ExitCode { let mut out_counter = 0; - while interface.is_coupling_ongoing() { - if interface.is_action_required(&precice::action_write_iteration_checkpoint()) { - interface - .pin_mut() - .mark_action_fulfilled(&precice::action_write_iteration_checkpoint()); + let mut t = 0.0; + + while participant.is_coupling_ongoing() { + if participant.pin_mut().requires_writing_checkpoint() { + // do nothing } + let dt = participant.get_max_time_step_size(); + fluid_compute_solution( &velocity_old[..], &pressure_old[..], @@ -356,23 +355,23 @@ fn main() -> ExitCode { &mut pressure[..], ); - interface + participant .pin_mut() - .write_block_scalar_data(pressure_id, &vertex_ids[..], &pressure[..]); + .write_data(mesh_name, "Pressure", &vertex_ids[..], &pressure[..]); + + participant.pin_mut().advance(dt); - interface.pin_mut().advance(dt); - interface.read_block_scalar_data( - cross_section_length_id, + participant.read_data( + mesh_name, + "CrossSectionLength", &vertex_ids[..], + participant.get_max_time_step_size(), &mut cross_section_length[..], ); - if interface.is_action_required(&precice::action_read_iteration_checkpoint()) { + if participant.pin_mut().requires_reading_checkpoint() { // i.e. fluid not yet converged - interface - .pin_mut() - .mark_action_fulfilled(&precice::action_read_iteration_checkpoint()); //pressure.copy_from_slice(&pressure_old[..]); //velocity.copy_from_slice(&velocity_old[..]); @@ -397,7 +396,7 @@ fn main() -> ExitCode { } println!("Exiting FluidSolver at t={}", t); - interface.pin_mut().finalize(); + participant.pin_mut().finalize(); return ExitCode::SUCCESS; } diff --git a/elastic-tube-1d/solid-rust/Cargo.lock b/elastic-tube-1d/solid-rust/Cargo.lock deleted file mode 100644 index 064988cbe..000000000 --- a/elastic-tube-1d/solid-rust/Cargo.lock +++ /dev/null @@ -1,274 +0,0 @@ -# This file is automatically @generated by Cargo. -# It is not intended for manual editing. -version = 3 - -[[package]] -name = "ansi_term" -version = "0.12.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d52a9bb7ec0cf484c551830a7ce27bd20d67eac647e1befb56b0be4ee39a55d2" -dependencies = [ - "winapi", -] - -[[package]] -name = "atty" -version = "0.2.14" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d9b39be18770d11421cdb1b9947a45dd3f37e93092cbf377614828a319d5fee8" -dependencies = [ - "hermit-abi", - "libc", - "winapi", -] - -[[package]] -name = "bitflags" -version = "1.3.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" - -[[package]] -name = "cc" -version = "1.0.79" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "50d30906286121d95be3d479533b458f87493b30a4b5f79a607db8f5d11aa91f" - -[[package]] -name = "clap" -version = "2.34.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a0610544180c38b88101fecf2dd634b174a62eef6946f84dfc6a7127512b381c" -dependencies = [ - "ansi_term", - "atty", - "bitflags", - "strsim", - "textwrap", - "unicode-width", - "vec_map", -] - -[[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 = "cxx" -version = "1.0.91" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "86d3488e7665a7a483b57e25bdd90d0aeb2bc7608c8d0346acf2ad3f1caf1d62" -dependencies = [ - "cc", - "cxxbridge-flags", - "cxxbridge-macro", - "link-cplusplus", -] - -[[package]] -name = "cxx-build" -version = "1.0.91" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "48fcaf066a053a41a81dfb14d57d99738b767febb8b735c3016e469fac5da690" -dependencies = [ - "cc", - "codespan-reporting", - "once_cell", - "proc-macro2", - "quote", - "scratch", - "syn", -] - -[[package]] -name = "cxxbridge-flags" -version = "1.0.91" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a2ef98b8b717a829ca5603af80e1f9e2e48013ab227b68ef37872ef84ee479bf" - -[[package]] -name = "cxxbridge-macro" -version = "1.0.91" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "086c685979a698443656e5cf7856c95c642295a38599f12fb1ff76fb28d19892" -dependencies = [ - "proc-macro2", - "quote", - "syn", -] - -[[package]] -name = "hermit-abi" -version = "0.1.19" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "62b467343b94ba476dcb2500d242dadbb39557df889310ac77c5d99100aaac33" -dependencies = [ - "libc", -] - -[[package]] -name = "libc" -version = "0.2.139" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "201de327520df007757c1f0adce6e827fe8562fbc28bfd9c15571c66ca1f5f79" - -[[package]] -name = "link-cplusplus" -version = "1.0.8" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ecd207c9c713c34f95a097a5b029ac2ce6010530c7b49d7fea24d977dede04f5" -dependencies = [ - "cc", -] - -[[package]] -name = "once_cell" -version = "1.17.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b7e5500299e16ebb147ae15a00a942af264cf3688f47923b8fc2cd5858f23ad3" - -[[package]] -name = "pkg-config" -version = "0.3.26" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6ac9a59f73473f1b8d852421e59e64809f025994837ef743615c6d0c5b305160" - -[[package]] -name = "precice" -version = "2.5.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "dec31fcf7edff2ec658584c7010a28da5f1e3d87b42b4f08a17df79d63305d91" -dependencies = [ - "clap", - "cxx", - "cxx-build", - "pkg-config", - "semver", -] - -[[package]] -name = "proc-macro2" -version = "1.0.51" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5d727cae5b39d21da60fa540906919ad737832fe0b1c165da3a34d6548c849d6" -dependencies = [ - "unicode-ident", -] - -[[package]] -name = "quote" -version = "1.0.23" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8856d8364d252a14d474036ea1358d63c9e6965c8e5c1885c18f73d70bff9c7b" -dependencies = [ - "proc-macro2", -] - -[[package]] -name = "scratch" -version = "1.0.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ddccb15bcce173023b3fedd9436f882a0739b8dfb45e4f6b6002bee5929f61b2" - -[[package]] -name = "semver" -version = "1.0.16" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "58bc9567378fc7690d6b2addae4e60ac2eeea07becb2c64b9f218b53865cba2a" - -[[package]] -name = "solid-rust" -version = "0.1.0" -dependencies = [ - "precice", -] - -[[package]] -name = "strsim" -version = "0.8.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8ea5119cdb4c55b55d432abb513a0429384878c15dde60cc77b1c99de1a95a6a" - -[[package]] -name = "syn" -version = "1.0.107" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1f4064b5b16e03ae50984a5a8ed5d4f8803e6bc1fd170a3cda91a1be4b18e3f5" -dependencies = [ - "proc-macro2", - "quote", - "unicode-ident", -] - -[[package]] -name = "termcolor" -version = "1.2.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "be55cf8942feac5c765c2c993422806843c9a9a45d4d5c407ad6dd2ea95eb9b6" -dependencies = [ - "winapi-util", -] - -[[package]] -name = "textwrap" -version = "0.11.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d326610f408c7a4eb6f51c37c330e496b08506c9457c9d34287ecc38809fb060" -dependencies = [ - "unicode-width", -] - -[[package]] -name = "unicode-ident" -version = "1.0.6" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "84a22b9f218b40614adcb3f4ff08b703773ad44fa9423e4e0d346d5db86e4ebc" - -[[package]] -name = "unicode-width" -version = "0.1.10" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c0edd1e5b14653f783770bce4a4dabb4a5108a5370a5f5d8cfe8710c361f6c8b" - -[[package]] -name = "vec_map" -version = "0.8.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f1bddf1187be692e79c5ffeab891132dfb0f236ed36a43c7ed39f1165ee20191" - -[[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.5" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "70ec6ce85bb158151cae5e5c87f95a8e97d2c0c4b001223f33a334e3ce5de178" -dependencies = [ - "winapi", -] - -[[package]] -name = "winapi-x86_64-pc-windows-gnu" -version = "0.4.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" diff --git a/elastic-tube-1d/solid-rust/Cargo.toml b/elastic-tube-1d/solid-rust/Cargo.toml index e74a65c3c..c9de9e33d 100644 --- a/elastic-tube-1d/solid-rust/Cargo.toml +++ b/elastic-tube-1d/solid-rust/Cargo.toml @@ -8,4 +8,4 @@ description = "The solid participant of the elastic-tube 1D tutorial" license = "LGPL-3.0-or-later" [dependencies] -"precice" = "2.5" +"precice" = "3" diff --git a/elastic-tube-1d/solid-rust/src/main.rs b/elastic-tube-1d/solid-rust/src/main.rs index b3d858e37..298b74c6d 100644 --- a/elastic-tube-1d/solid-rust/src/main.rs +++ b/elastic-tube-1d/solid-rust/src/main.rs @@ -28,14 +28,14 @@ fn main() -> ExitCode { const DOMAIN_SIZE: usize = 100; const CHUNK_SIZE: usize = DOMAIN_SIZE + 1; - let mut interface = precice::new("Solid", &config, 0, 1); + let mut participant = precice::new("Solid", &config, 0, 1); println!("preCICE configured..."); - let dimensions = interface.get_dimensions(); - let mesh_id = interface.get_mesh_id("Solid-Nodes-Mesh"); - let cross_section_length_id = interface.get_data_id("CrossSectionLength", mesh_id); - let pressure_id = interface.get_data_id("Pressure", mesh_id); + let mesh_name = "Solid-Nodes-Mesh"; + let dimensions = participant.get_mesh_dimensions(mesh_name); + assert!(participant.get_data_dimensions(mesh_name, "CrossSectionLength") == 1); + assert!(participant.get_data_dimensions(mesh_name, "Pressure") == 1); let mut pressure: Vec = vec![0.0; CHUNK_SIZE as usize]; let mut cross_section_length: Vec = vec![1.0; CHUNK_SIZE as usize]; @@ -54,61 +54,57 @@ fn main() -> ExitCode { let vertex_ids = { let mut ids = vec![-1; CHUNK_SIZE]; - interface + participant .pin_mut() - .set_mesh_vertices(mesh_id, &grid[..], &mut ids[..]); + .set_mesh_vertices(mesh_name, &grid[..], &mut ids[..]); ids }; - println!("Initializing preCICE..."); - - let mut t = 0.0; - let dt = interface.pin_mut().initialize(); - - if interface.is_action_required(&precice::action_write_initial_data()) { - interface.pin_mut().write_block_scalar_data( - cross_section_length_id, + if participant.pin_mut().requires_initial_data() { + participant.pin_mut().write_data( + mesh_name, + "CrossSectionLength", &vertex_ids[..], &cross_section_length[..], ); - interface - .pin_mut() - .mark_action_fulfilled(&precice::action_write_initial_data()); } - interface.pin_mut().initialize_data(); - while interface.is_coupling_ongoing() { - if interface.is_action_required(&precice::action_write_iteration_checkpoint()) { - interface - .pin_mut() - .mark_action_fulfilled(&precice::action_write_iteration_checkpoint()); + println!("Initializing preCICE..."); + + participant.pin_mut().initialize(); + + let mut t = 0.0; + while participant.is_coupling_ongoing() { + if participant.pin_mut().requires_writing_checkpoint() { + // no nothing } - interface.read_block_scalar_data(pressure_id, &vertex_ids[..], &mut pressure[..]); + let dt = participant.get_max_time_step_size(); + + participant.read_data(mesh_name, "Pressure", &vertex_ids[..], dt, &mut pressure[..]); solid_compute_solution(&pressure, &mut cross_section_length); - interface.pin_mut().write_block_scalar_data( - cross_section_length_id, + participant.pin_mut().write_data( + mesh_name, + "CrossSectionLength", &vertex_ids[..], &cross_section_length[..], ); - interface.pin_mut().advance(dt); + participant.pin_mut().advance(dt); - if interface.is_action_required(&precice::action_read_iteration_checkpoint()) { + if participant.pin_mut().requires_reading_checkpoint() { // i.e. fluid not yet converged - interface - .pin_mut() - .mark_action_fulfilled(&precice::action_read_iteration_checkpoint()); + // do nothing } else { t += dt; } } println!("Exiting SolidSolver at t={}", t); - interface.pin_mut().finalize(); + participant.pin_mut().finalize(); return ExitCode::SUCCESS; } From bc499ef013cf60e08d52b516e54067805820dbe2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Simonis?= Date: Wed, 13 Dec 2023 13:14:40 +0100 Subject: [PATCH 03/11] Migrate to wrapper --- elastic-tube-1d/fluid-rust/src/main.rs | 25 +++++++++++-------------- elastic-tube-1d/solid-rust/src/main.rs | 23 +++++++++++------------ 2 files changed, 22 insertions(+), 26 deletions(-) diff --git a/elastic-tube-1d/fluid-rust/src/main.rs b/elastic-tube-1d/fluid-rust/src/main.rs index 1fcbb81de..76fb7bbd1 100644 --- a/elastic-tube-1d/fluid-rust/src/main.rs +++ b/elastic-tube-1d/fluid-rust/src/main.rs @@ -6,6 +6,7 @@ use std::fs; use std::fs::File; use std::io::{BufWriter, Write}; use std::path; +use precice; fn write_vtk( file: &str, @@ -252,7 +253,7 @@ fn main() -> ExitCode { const DOMAIN_SIZE: usize = 100; const CHUNK_SIZE: usize = DOMAIN_SIZE + 1; - let mut participant = precice::new("Fluid", &config, 0, 1); + let mut participant = precice::Participant::new("Fluid", &config, 0, 1); println!("preCICE configured..."); @@ -293,16 +294,14 @@ fn main() -> ExitCode { let vertex_ids = { let mut ids = vec![-1; CHUNK_SIZE]; - participant - .pin_mut() - .set_mesh_vertices(mesh_name, &grid[..], &mut ids[..]); + participant.set_mesh_vertices(mesh_name, &grid[..], &mut ids[..]); ids }; println!("Initializing preCICE..."); - if participant.pin_mut().requires_initial_data() { - participant.pin_mut().write_data( + if participant.requires_initial_data() { + participant.write_data( mesh_name, "Pressure", &vertex_ids[..], @@ -310,7 +309,7 @@ fn main() -> ExitCode { ); } - participant.pin_mut().initialize(); + participant.initialize(); participant.read_data( mesh_name, @@ -336,7 +335,7 @@ fn main() -> ExitCode { let mut t = 0.0; while participant.is_coupling_ongoing() { - if participant.pin_mut().requires_writing_checkpoint() { + if participant.requires_writing_checkpoint() { // do nothing } @@ -355,11 +354,9 @@ fn main() -> ExitCode { &mut pressure[..], ); - participant - .pin_mut() - .write_data(mesh_name, "Pressure", &vertex_ids[..], &pressure[..]); + participant.write_data(mesh_name, "Pressure", &vertex_ids[..], &pressure[..]); - participant.pin_mut().advance(dt); + participant.advance(dt); participant.read_data( @@ -370,7 +367,7 @@ fn main() -> ExitCode { &mut cross_section_length[..], ); - if participant.pin_mut().requires_reading_checkpoint() { + if participant.requires_reading_checkpoint() { // i.e. fluid not yet converged //pressure.copy_from_slice(&pressure_old[..]); @@ -396,7 +393,7 @@ fn main() -> ExitCode { } println!("Exiting FluidSolver at t={}", t); - participant.pin_mut().finalize(); + participant.finalize(); return ExitCode::SUCCESS; } diff --git a/elastic-tube-1d/solid-rust/src/main.rs b/elastic-tube-1d/solid-rust/src/main.rs index 298b74c6d..eb4b752cf 100644 --- a/elastic-tube-1d/solid-rust/src/main.rs +++ b/elastic-tube-1d/solid-rust/src/main.rs @@ -1,5 +1,6 @@ use std::env; use std::process::ExitCode; +use precice; fn solid_compute_solution(pressure: &[f64], cross_section_length: &mut [f64]) { assert!(pressure.len() == cross_section_length.len()); @@ -28,7 +29,7 @@ fn main() -> ExitCode { const DOMAIN_SIZE: usize = 100; const CHUNK_SIZE: usize = DOMAIN_SIZE + 1; - let mut participant = precice::new("Solid", &config, 0, 1); + let mut participant = precice::Participant::new("Solid", &config, 0, 1); println!("preCICE configured..."); @@ -54,14 +55,12 @@ fn main() -> ExitCode { let vertex_ids = { let mut ids = vec![-1; CHUNK_SIZE]; - participant - .pin_mut() - .set_mesh_vertices(mesh_name, &grid[..], &mut ids[..]); + participant.set_mesh_vertices(mesh_name, &grid[..], &mut ids[..]); ids }; - if participant.pin_mut().requires_initial_data() { - participant.pin_mut().write_data( + if participant.requires_initial_data() { + participant.write_data( mesh_name, "CrossSectionLength", &vertex_ids[..], @@ -72,11 +71,11 @@ fn main() -> ExitCode { println!("Initializing preCICE..."); - participant.pin_mut().initialize(); + participant.initialize(); let mut t = 0.0; while participant.is_coupling_ongoing() { - if participant.pin_mut().requires_writing_checkpoint() { + if participant.requires_writing_checkpoint() { // no nothing } @@ -86,16 +85,16 @@ fn main() -> ExitCode { solid_compute_solution(&pressure, &mut cross_section_length); - participant.pin_mut().write_data( + participant.write_data( mesh_name, "CrossSectionLength", &vertex_ids[..], &cross_section_length[..], ); - participant.pin_mut().advance(dt); + participant.advance(dt); - if participant.pin_mut().requires_reading_checkpoint() { + if participant.requires_reading_checkpoint() { // i.e. fluid not yet converged // do nothing } else { @@ -104,7 +103,7 @@ fn main() -> ExitCode { } println!("Exiting SolidSolver at t={}", t); - participant.pin_mut().finalize(); + participant.finalize(); return ExitCode::SUCCESS; } From 041c192ab087110d435f57360744ad751a44308a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Simonis?= Date: Wed, 13 Dec 2023 13:16:15 +0100 Subject: [PATCH 04/11] Format rust solvers --- elastic-tube-1d/fluid-rust/src/main.rs | 10 ++-------- elastic-tube-1d/solid-rust/src/main.rs | 11 ++++++++--- 2 files changed, 10 insertions(+), 11 deletions(-) diff --git a/elastic-tube-1d/fluid-rust/src/main.rs b/elastic-tube-1d/fluid-rust/src/main.rs index 76fb7bbd1..84c096f04 100644 --- a/elastic-tube-1d/fluid-rust/src/main.rs +++ b/elastic-tube-1d/fluid-rust/src/main.rs @@ -2,11 +2,11 @@ use nalgebra as na; use std::env; use std::process::ExitCode; +use precice; use std::fs; use std::fs::File; use std::io::{BufWriter, Write}; use std::path; -use precice; fn write_vtk( file: &str, @@ -301,12 +301,7 @@ fn main() -> ExitCode { println!("Initializing preCICE..."); if participant.requires_initial_data() { - participant.write_data( - mesh_name, - "Pressure", - &vertex_ids[..], - &pressure[..], - ); + participant.write_data(mesh_name, "Pressure", &vertex_ids[..], &pressure[..]); } participant.initialize(); @@ -358,7 +353,6 @@ fn main() -> ExitCode { participant.advance(dt); - participant.read_data( mesh_name, "CrossSectionLength", diff --git a/elastic-tube-1d/solid-rust/src/main.rs b/elastic-tube-1d/solid-rust/src/main.rs index eb4b752cf..831ca4071 100644 --- a/elastic-tube-1d/solid-rust/src/main.rs +++ b/elastic-tube-1d/solid-rust/src/main.rs @@ -1,6 +1,6 @@ +use precice; use std::env; use std::process::ExitCode; -use precice; fn solid_compute_solution(pressure: &[f64], cross_section_length: &mut [f64]) { assert!(pressure.len() == cross_section_length.len()); @@ -68,7 +68,6 @@ fn main() -> ExitCode { ); } - println!("Initializing preCICE..."); participant.initialize(); @@ -81,7 +80,13 @@ fn main() -> ExitCode { let dt = participant.get_max_time_step_size(); - participant.read_data(mesh_name, "Pressure", &vertex_ids[..], dt, &mut pressure[..]); + participant.read_data( + mesh_name, + "Pressure", + &vertex_ids[..], + dt, + &mut pressure[..], + ); solid_compute_solution(&pressure, &mut cross_section_length); From 6c29feaf723bbcef4e36434791c14e8fd18ed974 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Simonis?= Date: Thu, 11 Jan 2024 10:36:39 +0100 Subject: [PATCH 05/11] Delete the lock --- elastic-tube-1d/fluid-rust/Cargo.lock | 421 -------------------------- 1 file changed, 421 deletions(-) delete mode 100644 elastic-tube-1d/fluid-rust/Cargo.lock diff --git a/elastic-tube-1d/fluid-rust/Cargo.lock b/elastic-tube-1d/fluid-rust/Cargo.lock deleted file mode 100644 index a62a2e950..000000000 --- a/elastic-tube-1d/fluid-rust/Cargo.lock +++ /dev/null @@ -1,421 +0,0 @@ -# This file is automatically @generated by Cargo. -# It is not intended for manual editing. -version = 3 - -[[package]] -name = "ansi_term" -version = "0.12.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d52a9bb7ec0cf484c551830a7ce27bd20d67eac647e1befb56b0be4ee39a55d2" -dependencies = [ - "winapi", -] - -[[package]] -name = "approx" -version = "0.5.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" -dependencies = [ - "num-traits", -] - -[[package]] -name = "atty" -version = "0.2.14" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d9b39be18770d11421cdb1b9947a45dd3f37e93092cbf377614828a319d5fee8" -dependencies = [ - "hermit-abi", - "libc", - "winapi", -] - -[[package]] -name = "autocfg" -version = "1.1.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa" - -[[package]] -name = "bitflags" -version = "1.3.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" - -[[package]] -name = "bytemuck" -version = "1.12.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2f5715e491b5a1598fc2bef5a606847b5dc1d48ea625bd3c02c00de8285591da" - -[[package]] -name = "cc" -version = "1.0.73" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2fff2a6927b3bb87f9595d67196a70493f627687a71d87a0d692242c33f58c11" - -[[package]] -name = "clap" -version = "2.34.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a0610544180c38b88101fecf2dd634b174a62eef6946f84dfc6a7127512b381c" -dependencies = [ - "ansi_term", - "atty", - "bitflags", - "strsim", - "textwrap", - "unicode-width", - "vec_map", -] - -[[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 = "cxx" -version = "1.0.78" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "19f39818dcfc97d45b03953c1292efc4e80954e1583c4aa770bac1383e2310a4" -dependencies = [ - "cc", - "cxxbridge-flags", - "cxxbridge-macro", - "link-cplusplus", -] - -[[package]] -name = "cxx-build" -version = "1.0.78" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3e580d70777c116df50c390d1211993f62d40302881e54d4b79727acb83d0199" -dependencies = [ - "cc", - "codespan-reporting", - "once_cell", - "proc-macro2", - "quote", - "scratch", - "syn", -] - -[[package]] -name = "cxxbridge-flags" -version = "1.0.78" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "56a46460b88d1cec95112c8c363f0e2c39afdb237f60583b0b36343bf627ea9c" - -[[package]] -name = "cxxbridge-macro" -version = "1.0.78" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "747b608fecf06b0d72d440f27acc99288207324b793be2c17991839f3d4995ea" -dependencies = [ - "proc-macro2", - "quote", - "syn", -] - -[[package]] -name = "fluid-rust" -version = "0.1.0" -dependencies = [ - "nalgebra", - "precice", -] - -[[package]] -name = "hermit-abi" -version = "0.1.19" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "62b467343b94ba476dcb2500d242dadbb39557df889310ac77c5d99100aaac33" -dependencies = [ - "libc", -] - -[[package]] -name = "libc" -version = "0.2.139" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "201de327520df007757c1f0adce6e827fe8562fbc28bfd9c15571c66ca1f5f79" - -[[package]] -name = "link-cplusplus" -version = "1.0.7" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9272ab7b96c9046fbc5bc56c06c117cb639fe2d509df0c421cad82d2915cf369" -dependencies = [ - "cc", -] - -[[package]] -name = "matrixmultiply" -version = "0.3.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "add85d4dd35074e6fedc608f8c8f513a3548619a9024b751949ef0e8e45a4d84" -dependencies = [ - "rawpointer", -] - -[[package]] -name = "nalgebra" -version = "0.31.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e9e0a04ce089f9401aac565c740ed30c46291260f27d4911fdbaa6ca65fa3044" -dependencies = [ - "approx", - "matrixmultiply", - "nalgebra-macros", - "num-complex", - "num-rational", - "num-traits", - "simba", - "typenum", -] - -[[package]] -name = "nalgebra-macros" -version = "0.1.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "01fcc0b8149b4632adc89ac3b7b31a12fb6099a0317a4eb2ebff574ef7de7218" -dependencies = [ - "proc-macro2", - "quote", - "syn", -] - -[[package]] -name = "num-complex" -version = "0.4.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7ae39348c8bc5fbd7f40c727a9925f03517afd2ab27d46702108b6a7e5414c19" -dependencies = [ - "num-traits", -] - -[[package]] -name = "num-integer" -version = "0.1.45" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "225d3389fb3509a24c93f5c29eb6bde2586b98d9f016636dff58d7c6f7569cd9" -dependencies = [ - "autocfg", - "num-traits", -] - -[[package]] -name = "num-rational" -version = "0.4.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0638a1c9d0a3c0914158145bc76cff373a75a627e6ecbfb71cbe6f453a5a19b0" -dependencies = [ - "autocfg", - "num-integer", - "num-traits", -] - -[[package]] -name = "num-traits" -version = "0.2.15" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "578ede34cf02f8924ab9447f50c28075b4d3e5b269972345e7e0372b38c6cdcd" -dependencies = [ - "autocfg", -] - -[[package]] -name = "once_cell" -version = "1.15.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e82dad04139b71a90c080c8463fe0dc7902db5192d939bd0950f074d014339e1" - -[[package]] -name = "paste" -version = "1.0.9" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b1de2e551fb905ac83f73f7aedf2f0cb4a0da7e35efa24a202a936269f1f18e1" - -[[package]] -name = "pkg-config" -version = "0.3.25" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1df8c4ec4b0627e53bdf214615ad287367e482558cf84b109250b37464dc03ae" - -[[package]] -name = "precice" -version = "2.5.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "dec31fcf7edff2ec658584c7010a28da5f1e3d87b42b4f08a17df79d63305d91" -dependencies = [ - "clap", - "cxx", - "cxx-build", - "pkg-config", - "semver", -] - -[[package]] -name = "proc-macro2" -version = "1.0.46" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "94e2ef8dbfc347b10c094890f778ee2e36ca9bb4262e86dc99cd217e35f3470b" -dependencies = [ - "unicode-ident", -] - -[[package]] -name = "quote" -version = "1.0.21" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "bbe448f377a7d6961e30f5955f9b8d106c3f5e449d493ee1b125c1d43c2b5179" -dependencies = [ - "proc-macro2", -] - -[[package]] -name = "rawpointer" -version = "0.2.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" - -[[package]] -name = "safe_arch" -version = "0.6.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "794821e4ccb0d9f979512f9c1973480123f9bd62a90d74ab0f9426fcf8f4a529" -dependencies = [ - "bytemuck", -] - -[[package]] -name = "scratch" -version = "1.0.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9c8132065adcfd6e02db789d9285a0deb2f3fcb04002865ab67d5fb103533898" - -[[package]] -name = "semver" -version = "1.0.16" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "58bc9567378fc7690d6b2addae4e60ac2eeea07becb2c64b9f218b53865cba2a" - -[[package]] -name = "simba" -version = "0.7.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c48e45e5961033db030b56ad67aef22e9c908c493a6e8348c0a0f6b93433cd77" -dependencies = [ - "approx", - "num-complex", - "num-traits", - "paste", - "wide", -] - -[[package]] -name = "strsim" -version = "0.8.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8ea5119cdb4c55b55d432abb513a0429384878c15dde60cc77b1c99de1a95a6a" - -[[package]] -name = "syn" -version = "1.0.101" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e90cde112c4b9690b8cbe810cba9ddd8bc1d7472e2cae317b69e9438c1cba7d2" -dependencies = [ - "proc-macro2", - "quote", - "unicode-ident", -] - -[[package]] -name = "termcolor" -version = "1.1.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "bab24d30b911b2376f3a13cc2cd443142f0c81dda04c118693e35b3835757755" -dependencies = [ - "winapi-util", -] - -[[package]] -name = "textwrap" -version = "0.11.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d326610f408c7a4eb6f51c37c330e496b08506c9457c9d34287ecc38809fb060" -dependencies = [ - "unicode-width", -] - -[[package]] -name = "typenum" -version = "1.15.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "dcf81ac59edc17cc8697ff311e8f5ef2d99fcbd9817b34cec66f90b6c3dfd987" - -[[package]] -name = "unicode-ident" -version = "1.0.4" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "dcc811dc4066ac62f84f11307873c4850cb653bfa9b1719cee2bd2204a4bc5dd" - -[[package]] -name = "unicode-width" -version = "0.1.10" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c0edd1e5b14653f783770bce4a4dabb4a5108a5370a5f5d8cfe8710c361f6c8b" - -[[package]] -name = "vec_map" -version = "0.8.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f1bddf1187be692e79c5ffeab891132dfb0f236ed36a43c7ed39f1165ee20191" - -[[package]] -name = "wide" -version = "0.7.4" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b3aba2d1dac31ac7cae82847ac5b8be822aee8f99a4e100f279605016b185c5f" -dependencies = [ - "bytemuck", - "safe_arch", -] - -[[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.5" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "70ec6ce85bb158151cae5e5c87f95a8e97d2c0c4b001223f33a334e3ce5de178" -dependencies = [ - "winapi", -] - -[[package]] -name = "winapi-x86_64-pc-windows-gnu" -version = "0.4.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" From 01370aa5912e918c6a07d0d6754b7406eefe0bad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Simonis?= Date: Thu, 11 Jan 2024 14:06:36 +0100 Subject: [PATCH 06/11] Fix et1d fluid rust output naming --- elastic-tube-1d/fluid-rust/src/main.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/elastic-tube-1d/fluid-rust/src/main.rs b/elastic-tube-1d/fluid-rust/src/main.rs index 84c096f04..bd49f463f 100644 --- a/elastic-tube-1d/fluid-rust/src/main.rs +++ b/elastic-tube-1d/fluid-rust/src/main.rs @@ -372,7 +372,7 @@ fn main() -> ExitCode { velocity_old.copy_from_slice(&velocity[..]); cross_section_length_old.copy_from_slice(&cross_section_length[..]); - let filename = format!("./output/out_fluid{out_counter}.vtk"); + let filename = format!("./output/out_fluid_{out_counter}.vtk"); println!("writing timestep at t={} to {}", t, filename); write_vtk( &filename, From fa9e15b2ae82ea629fa607332c8236314f685c34 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Simonis?= Date: Thu, 11 Jan 2024 14:07:04 +0100 Subject: [PATCH 07/11] Add et1d rust fluid to plot --- elastic-tube-1d/plot-diameter.sh | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/elastic-tube-1d/plot-diameter.sh b/elastic-tube-1d/plot-diameter.sh index f4d31bb82..8d1b17c4e 100755 --- a/elastic-tube-1d/plot-diameter.sh +++ b/elastic-tube-1d/plot-diameter.sh @@ -16,3 +16,10 @@ if [ -n "$(ls -A ./fluid-python/output/*.vtk 2>/dev/null)" ]; then else echo "No results to plot from fluid-python." fi + +# Plot diameter from fluid-rust +if [ -n "$(ls -A ./fluid-rust/output/*.vtk 2>/dev/null)" ]; then + python3 plot-vtk.py diameter fluid-rust/output/out_fluid_ & +else + echo "No results to plot from fluid-python." +fi From c6790f3992349c71bb279c00b6c637bfd01ec6dd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Simonis?= Date: Thu, 11 Jan 2024 14:07:18 +0100 Subject: [PATCH 08/11] Fix coordinates in et1d solid rust --- elastic-tube-1d/solid-rust/src/main.rs | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/elastic-tube-1d/solid-rust/src/main.rs b/elastic-tube-1d/solid-rust/src/main.rs index 831ca4071..f85f22f1d 100644 --- a/elastic-tube-1d/solid-rust/src/main.rs +++ b/elastic-tube-1d/solid-rust/src/main.rs @@ -28,6 +28,7 @@ fn main() -> ExitCode { const DOMAIN_SIZE: usize = 100; const CHUNK_SIZE: usize = DOMAIN_SIZE + 1; + const TUBE_LENGTH: f64 = 10.0; let mut participant = precice::Participant::new("Solid", &config, 0, 1); @@ -44,11 +45,9 @@ fn main() -> ExitCode { let grid_size = CHUNK_SIZE * dimensions as usize; let grid: Vec = { let mut v: Vec = vec![0_f64; grid_size]; + const DX : f64 = TUBE_LENGTH / DOMAIN_SIZE as f64; for i in 0..CHUNK_SIZE - 1 { - for j in 0..(dimensions as usize) - 1 { - let idx = i * dimensions as usize + j; - v[idx] = (i * (1 - j)) as f64; - } + v[i * dimensions as usize] = DX * i as f64; } v }; From 565b485170e859a661d9bdc3bfda1fa778967428 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Simonis?= Date: Fri, 2 Feb 2024 13:40:09 +0100 Subject: [PATCH 09/11] Ignore Cargo.locks --- elastic-tube-1d/.gitignore | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/elastic-tube-1d/.gitignore b/elastic-tube-1d/.gitignore index 584556938..71c464a3c 100644 --- a/elastic-tube-1d/.gitignore +++ b/elastic-tube-1d/.gitignore @@ -1,9 +1,7 @@ -fluid-cpp/build/ -fluid-cpp/output/*.vtk -solid-cpp/build/ -fluid-python/__pycache__/ -fluid-python/output/*.vtk -solid-python/__pycache__/ +*-cpp/build/ +*-python/__pycache__/ +*-rust/Cargo.lock +fluid-*/output/*.vtk *.o *.log From e8a82c80aea9ff8c47f8834e9f686f6f5731c2e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Simonis?= Date: Thu, 15 Feb 2024 10:54:08 +0100 Subject: [PATCH 10/11] Cleanup et1d rust solvers --- elastic-tube-1d/fluid-rust/src/main.rs | 243 +---------------------- elastic-tube-1d/fluid-rust/src/solver.rs | 172 ++++++++++++++++ elastic-tube-1d/fluid-rust/src/utils.rs | 62 ++++++ elastic-tube-1d/solid-rust/src/main.rs | 16 +- elastic-tube-1d/solid-rust/src/solver.rs | 11 + 5 files changed, 253 insertions(+), 251 deletions(-) create mode 100644 elastic-tube-1d/fluid-rust/src/solver.rs create mode 100644 elastic-tube-1d/fluid-rust/src/utils.rs create mode 100644 elastic-tube-1d/solid-rust/src/solver.rs diff --git a/elastic-tube-1d/fluid-rust/src/main.rs b/elastic-tube-1d/fluid-rust/src/main.rs index bd49f463f..c4d83b3fa 100644 --- a/elastic-tube-1d/fluid-rust/src/main.rs +++ b/elastic-tube-1d/fluid-rust/src/main.rs @@ -1,242 +1,9 @@ -use nalgebra as na; +use precice; use std::env; use std::process::ExitCode; -use precice; -use std::fs; -use std::fs::File; -use std::io::{BufWriter, Write}; -use std::path; - -fn write_vtk( - file: &str, - dx: f64, - velocity: &[f64], - pressure: &[f64], - cross_section_length: &[f64], -) -> std::io::Result<()> { - let filepath = path::Path::new(file); - if !filepath.parent().unwrap().exists() { - fs::create_dir(filepath.parent().unwrap())?; - } - - // let mut f = File::create(file).expect("Unable to open vtk file"); - let file = File::create(file).expect("Unable to open vtk file"); - let mut f = BufWriter::new(file); - - let n_points = velocity.len(); - - write!(f, "# vtk DataFile Version 2.0\n\n")?; - write!(f, "ASCII\n\n")?; - write!(f, "DATASET UNSTRUCTURED_GRID\n\n")?; - write!(f, "POINTS {n_points} float\n\n")?; - for i in 0..n_points { - write!( - f, - "{:.16e} 0.0000000000000000e+00 0.0000000000000000e+00\n", - i as f64 * dx - )?; - } - writeln!(f)?; - - write!(f, "POINT_DATA {n_points}\n\n")?; - - write!(f, "VECTORS velocity float\n")?; - for v in velocity { - write!( - f, - "{:.16e} 0.0000000000000000e+00 0.0000000000000000e+00\n", - v - )?; - } - writeln!(f)?; - - write!(f, "SCALARS pressure float\n")?; - write!(f, "LOOKUP_TABLE default\n")?; - for v in pressure { - write!(f, "{:.16e}\n", v)?; - } - writeln!(f)?; - - write!(f, "SCALARS diameter float\n")?; - write!(f, "LOOKUP_TABLE default\n")?; - for v in cross_section_length { - write!(f, "{:.16e}\n", v)?; - } - writeln!(f)?; - Ok(()) -} - -fn fluid_compute_solution( - velocity_old: &[f64], - pressure_old: &[f64], - cross_section_length_old: &[f64], - cross_section_length: &[f64], - t: f64, - n: usize, - kappa: f64, - tau: f64, - velocity: &mut [f64], - pressure: &mut [f64], -) { - // Initial guess - pressure.copy_from_slice(&pressure_old[..]); - - const E: f64 = 10000_f64; - //const C_MK2 : f64 = E / std::f64::consts::FRAC_2_SQRT_PI; - let c_mk2: f64 = E / 2.0 * std::f64::consts::PI.sqrt(); - - // lhs = 2*N+2 - - const ALPHA: f64 = 0.0; - const L: f64 = 10.0; - let dx = L / kappa; - - let equations = 2 * n + 2; - - let mut k = 0; - loop { - k += 1; - - let res = { - let mut res = na::DVector::from_element(equations, 0.0); - for i in 1..n { - let [p_l, p_c, p_r]: [f64; 3] = pressure[i - 1..i + 2].try_into().unwrap(); - let [v_l, v_c, v_r]: [f64; 3] = velocity[i - 1..i + 2].try_into().unwrap(); - let [c_l, c_c, c_r]: [f64; 3] = - cross_section_length[i - 1..i + 2].try_into().unwrap(); - - /* Momentum */ - res[i] = (velocity_old[i] * cross_section_length_old[i] - v_c * c_c) * dx / tau - + 0.25 * (-c_r * v_c * v_r - c_c * v_c * v_r) - + 0.25 - * (-c_r * v_c * v_c - c_c * v_c * v_c + c_c * v_l * v_c + c_l * v_l * v_c) - + 0.25 * (c_l * v_l * v_l + c_c * v_l * v_l) - + 0.25 - * (c_l * p_l + c_c * p_l - c_l * p_c + c_r * p_c - c_c * p_r - c_r * p_r); - - /* Continuity */ - res[i + n + 1] = (cross_section_length_old[i] - c_c) * dx / tau - + 0.25 - * (c_l * v_l + c_c * v_l + c_l * v_c - c_r * v_c - c_c * v_r - c_r * v_r) - + ALPHA * (p_l - 2.0 * p_c + p_r); - } - /* Boundary */ - - /* Velocity Inlet is prescribed */ - let u0 = 10.0; - let ampl = 3.0; - let frequency = 10.0; - let t_shift = 0.0; - let velocity_in = - u0 + ampl * (frequency * (t + tau + t_shift) * std::f64::consts::PI).sin(); - res[0] = velocity_in - velocity[0]; - - /* Pressure Inlet is linearly interpolated */ - res[n + 1] = -pressure[0] + 2.0 * pressure[1] - pressure[2]; - - /* Velocity Outlet is linearly interpolated */ - res[n] = -velocity[n] + 2.0 * velocity[n - 1] - velocity[n - 2]; - - /* Pressure Outlet is "non-reflecting" */ - let tmp2 = - (c_mk2 - pressure_old[n] / 2.0).sqrt() - (velocity[n] - velocity_old[n]) / 4.0; - res[2 * n + 1] = -pressure[n] + 2.0 * (c_mk2 - tmp2.powi(2)); - res - }; - - // compute norm of residual - let norm: f64 = { - let norm1 = res.map(|x: f64| x * x).sum().sqrt(); - let norm2 = pressure - .iter() - .chain(velocity.iter()) - .map(|x| x * x) - .sum::() - .sqrt(); - norm1 / norm2 - }; - - const TOL: f64 = 1e-10; - const MAX_ITER: usize = 1000; - if norm < TOL && k > 1 { - //println!("Success at k={} res={} tol={}", k, norm, TOL); - break; - } else if k > MAX_ITER { - panic!("Nonlinear Solver break, iterations: {k}, residual norm: {norm}"); - } else { - //println!("Continuing k={} res={:.10} tol={:.10}", k, norm, TOL); - } - - /* Initializing the the LHS i.e. Left Hand Side */ - let lhs = { - let mut lhs = na::DMatrix::from_element(equations, equations, 0.0); - for i in 1..n { - let [v_l, v_c, v_r]: [f64; 3] = velocity[i - 1..i + 2].try_into().unwrap(); - let [c_l, c_c, c_r]: [f64; 3] = - cross_section_length[i - 1..i + 2].try_into().unwrap(); - - // Momentum, Velocity - lhs[(i, i - 1)] = - 0.25 * (-2.0 * c_l * v_l - 2.0 * c_c * v_l - c_c * v_c - c_l * v_c); - - lhs[(i, i)] = c_c * dx / tau - + 0.25 - * (c_r * v_r + c_c * v_r + 2.0 * c_r * v_c + 2.0 * c_c * v_c - - c_c * v_l - - c_l * v_l); - lhs[(i, i + 1)] = 0.25 * (c_r * v_c + c_c * v_c); - - // Momentum, Pressure - lhs[(i, n + 1 + i - 1)] = -0.25 * (c_l + c_c); - lhs[(i, n + 1 + i)] = 0.25 * (c_l - c_r); - lhs[(i, n + 1 + i + 1)] = 0.25 * (c_c + c_r); - - // Continuity, Velocity - lhs[(i + n + 1, i - 1)] = -0.25 * (c_l + c_c); - lhs[(i + n + 1, i)] = -0.25 * (c_l - c_r); - lhs[(i + n + 1, i + 1)] = 0.25 * (c_c + c_r); - - // Continuity, Pressure - lhs[(i + n + 1, n + 1 + i - 1)] = -ALPHA; - lhs[(i + n + 1, n + 1 + i)] = 2.0 * ALPHA; - lhs[(i + n + 1, n + 1 + i + 1)] = -ALPHA; - } - - /* Boundary */ - - // Velocity Inlet is prescribed - lhs[(0, 0)] = 1.0; - // Pressure Inlet is linearly interpolated - lhs[(n + 1, n + 1)] = 1.0; - lhs[(n + 1, n + 2)] = -2.0; - lhs[(n + 1, n + 3)] = 1.0; - // Velocity Outlet is linearly interpolated - lhs[(n, n)] = 1.0; - lhs[(n, n - 1)] = -2.0; - lhs[(n, n - 2)] = 1.0; - // Pressure Outlet is non-Reflecting - lhs[(2 * n + 1, 2 * n + 1)] = 1.0; - lhs[(2 * n + 1, n)] = - -((c_mk2 - pressure_old[n] / 2.0).sqrt() - (velocity[n] - velocity_old[n]) / 4.0); - lhs - }; - - let sol = lhs - .lu() - .solve(&res) - .expect("Linear Solver did not converge!"); - - velocity - .iter_mut() - .zip(sol.as_slice()[..n + 1].iter()) - .for_each(|p| *p.0 += p.1); - pressure - .iter_mut() - .zip(sol.as_slice()[n + 1..].iter()) - .for_each(|p| *p.0 += p.1); - } -} +mod solver; +mod utils; fn main() -> ExitCode { println!("Starting Fluid Solver..."); @@ -336,7 +103,7 @@ fn main() -> ExitCode { let dt = participant.get_max_time_step_size(); - fluid_compute_solution( + solver::fluid_compute_solution( &velocity_old[..], &pressure_old[..], &cross_section_length_old[..], @@ -374,7 +141,7 @@ fn main() -> ExitCode { let filename = format!("./output/out_fluid_{out_counter}.vtk"); println!("writing timestep at t={} to {}", t, filename); - write_vtk( + utils::write_vtk( &filename, DX, &velocity_old[..], diff --git a/elastic-tube-1d/fluid-rust/src/solver.rs b/elastic-tube-1d/fluid-rust/src/solver.rs new file mode 100644 index 000000000..8f57cc3f2 --- /dev/null +++ b/elastic-tube-1d/fluid-rust/src/solver.rs @@ -0,0 +1,172 @@ +use nalgebra::{DMatrix, DVector}; + +pub fn fluid_compute_solution( + velocity_old: &[f64], + pressure_old: &[f64], + cross_section_length_old: &[f64], + cross_section_length: &[f64], + t: f64, + n: usize, + kappa: f64, + tau: f64, + velocity: &mut [f64], + pressure: &mut [f64], +) { + // Initial guess + pressure.copy_from_slice(&pressure_old[..]); + + const E: f64 = 10000_f64; + //const C_MK2 : f64 = E / std::f64::consts::FRAC_2_SQRT_PI; + let c_mk2: f64 = E / 2.0 * std::f64::consts::PI.sqrt(); + + // lhs = 2*N+2 + + const ALPHA: f64 = 0.0; + const L: f64 = 10.0; + let dx = L / kappa; + + let equations = 2 * n + 2; + + let mut k = 0; + loop { + k += 1; + + let res = { + let mut res = DVector::from_element(equations, 0.0); + for i in 1..n { + let [p_l, p_c, p_r]: [f64; 3] = pressure[i - 1..i + 2].try_into().unwrap(); + let [v_l, v_c, v_r]: [f64; 3] = velocity[i - 1..i + 2].try_into().unwrap(); + let [c_l, c_c, c_r]: [f64; 3] = + cross_section_length[i - 1..i + 2].try_into().unwrap(); + + /* Momentum */ + res[i] = (velocity_old[i] * cross_section_length_old[i] - v_c * c_c) * dx / tau + + 0.25 * (-c_r * v_c * v_r - c_c * v_c * v_r) + + 0.25 + * (-c_r * v_c * v_c - c_c * v_c * v_c + c_c * v_l * v_c + c_l * v_l * v_c) + + 0.25 * (c_l * v_l * v_l + c_c * v_l * v_l) + + 0.25 + * (c_l * p_l + c_c * p_l - c_l * p_c + c_r * p_c - c_c * p_r - c_r * p_r); + + /* Continuity */ + res[i + n + 1] = (cross_section_length_old[i] - c_c) * dx / tau + + 0.25 + * (c_l * v_l + c_c * v_l + c_l * v_c - c_r * v_c - c_c * v_r - c_r * v_r) + + ALPHA * (p_l - 2.0 * p_c + p_r); + } + /* Boundary */ + + /* Velocity Inlet is prescribed */ + let u0 = 10.0; + let ampl = 3.0; + let frequency = 10.0; + let t_shift = 0.0; + let velocity_in = + u0 + ampl * (frequency * (t + tau + t_shift) * std::f64::consts::PI).sin(); + res[0] = velocity_in - velocity[0]; + + /* Pressure Inlet is linearly interpolated */ + res[n + 1] = -pressure[0] + 2.0 * pressure[1] - pressure[2]; + + /* Velocity Outlet is linearly interpolated */ + res[n] = -velocity[n] + 2.0 * velocity[n - 1] - velocity[n - 2]; + + /* Pressure Outlet is "non-reflecting" */ + let tmp2 = + (c_mk2 - pressure_old[n] / 2.0).sqrt() - (velocity[n] - velocity_old[n]) / 4.0; + res[2 * n + 1] = -pressure[n] + 2.0 * (c_mk2 - tmp2.powi(2)); + res + }; + + // compute norm of residual + let norm: f64 = { + let norm1 = res.map(|x: f64| x * x).sum().sqrt(); + let norm2 = pressure + .iter() + .chain(velocity.iter()) + .map(|x| x * x) + .sum::() + .sqrt(); + norm1 / norm2 + }; + + const TOL: f64 = 1e-10; + const MAX_ITER: usize = 1000; + if norm < TOL && k > 1 { + //println!("Success at k={} res={} tol={}", k, norm, TOL); + break; + } else if k > MAX_ITER { + panic!("Nonlinear Solver break, iterations: {k}, residual norm: {norm}"); + } else { + //println!("Continuing k={} res={:.10} tol={:.10}", k, norm, TOL); + } + + /* Initializing the the LHS i.e. Left Hand Side */ + let lhs = { + let mut lhs = DMatrix::from_element(equations, equations, 0.0); + for i in 1..n { + let [v_l, v_c, v_r]: [f64; 3] = velocity[i - 1..i + 2].try_into().unwrap(); + let [c_l, c_c, c_r]: [f64; 3] = + cross_section_length[i - 1..i + 2].try_into().unwrap(); + + // Momentum, Velocity + lhs[(i, i - 1)] = + 0.25 * (-2.0 * c_l * v_l - 2.0 * c_c * v_l - c_c * v_c - c_l * v_c); + + lhs[(i, i)] = c_c * dx / tau + + 0.25 + * (c_r * v_r + c_c * v_r + 2.0 * c_r * v_c + 2.0 * c_c * v_c + - c_c * v_l + - c_l * v_l); + lhs[(i, i + 1)] = 0.25 * (c_r * v_c + c_c * v_c); + + // Momentum, Pressure + lhs[(i, n + 1 + i - 1)] = -0.25 * (c_l + c_c); + lhs[(i, n + 1 + i)] = 0.25 * (c_l - c_r); + lhs[(i, n + 1 + i + 1)] = 0.25 * (c_c + c_r); + + // Continuity, Velocity + lhs[(i + n + 1, i - 1)] = -0.25 * (c_l + c_c); + lhs[(i + n + 1, i)] = -0.25 * (c_l - c_r); + lhs[(i + n + 1, i + 1)] = 0.25 * (c_c + c_r); + + // Continuity, Pressure + lhs[(i + n + 1, n + 1 + i - 1)] = -ALPHA; + lhs[(i + n + 1, n + 1 + i)] = 2.0 * ALPHA; + lhs[(i + n + 1, n + 1 + i + 1)] = -ALPHA; + } + + /* Boundary */ + + // Velocity Inlet is prescribed + lhs[(0, 0)] = 1.0; + // Pressure Inlet is linearly interpolated + lhs[(n + 1, n + 1)] = 1.0; + lhs[(n + 1, n + 2)] = -2.0; + lhs[(n + 1, n + 3)] = 1.0; + // Velocity Outlet is linearly interpolated + lhs[(n, n)] = 1.0; + lhs[(n, n - 1)] = -2.0; + lhs[(n, n - 2)] = 1.0; + // Pressure Outlet is non-Reflecting + lhs[(2 * n + 1, 2 * n + 1)] = 1.0; + lhs[(2 * n + 1, n)] = + -((c_mk2 - pressure_old[n] / 2.0).sqrt() - (velocity[n] - velocity_old[n]) / 4.0); + lhs + }; + + let sol = lhs + .lu() + .solve(&res) + .expect("Linear Solver did not converge!"); + + velocity + .iter_mut() + .zip(sol.as_slice()[..n + 1].iter()) + .for_each(|p| *p.0 += p.1); + pressure + .iter_mut() + .zip(sol.as_slice()[n + 1..].iter()) + .for_each(|p| *p.0 += p.1); + } +} diff --git a/elastic-tube-1d/fluid-rust/src/utils.rs b/elastic-tube-1d/fluid-rust/src/utils.rs new file mode 100644 index 000000000..f593a6ee4 --- /dev/null +++ b/elastic-tube-1d/fluid-rust/src/utils.rs @@ -0,0 +1,62 @@ +use std::fs::{create_dir, File}; +use std::io::{BufWriter, Write}; +use std::path; + +pub fn write_vtk( + file: &str, + dx: f64, + velocity: &[f64], + pressure: &[f64], + cross_section_length: &[f64], +) -> std::io::Result<()> { + let filepath = path::Path::new(file); + if !filepath.parent().unwrap().exists() { + create_dir(filepath.parent().unwrap())?; + } + + // let mut f = File::create(file).expect("Unable to open vtk file"); + let file = File::create(file).expect("Unable to open vtk file"); + let mut f = BufWriter::new(file); + + let n_points = velocity.len(); + + write!(f, "# vtk DataFile Version 2.0\n\n")?; + write!(f, "ASCII\n\n")?; + write!(f, "DATASET UNSTRUCTURED_GRID\n\n")?; + write!(f, "POINTS {n_points} float\n\n")?; + for i in 0..n_points { + write!( + f, + "{:.16e} 0.0000000000000000e+00 0.0000000000000000e+00\n", + i as f64 * dx + )?; + } + writeln!(f)?; + + write!(f, "POINT_DATA {n_points}\n\n")?; + + write!(f, "VECTORS velocity float\n")?; + for v in velocity { + write!( + f, + "{:.16e} 0.0000000000000000e+00 0.0000000000000000e+00\n", + v + )?; + } + writeln!(f)?; + + write!(f, "SCALARS pressure float\n")?; + write!(f, "LOOKUP_TABLE default\n")?; + for v in pressure { + write!(f, "{:.16e}\n", v)?; + } + writeln!(f)?; + + write!(f, "SCALARS diameter float\n")?; + write!(f, "LOOKUP_TABLE default\n")?; + for v in cross_section_length { + write!(f, "{:.16e}\n", v)?; + } + writeln!(f)?; + Ok(()) +} diff --git a/elastic-tube-1d/solid-rust/src/main.rs b/elastic-tube-1d/solid-rust/src/main.rs index f85f22f1d..c7a2417fd 100644 --- a/elastic-tube-1d/solid-rust/src/main.rs +++ b/elastic-tube-1d/solid-rust/src/main.rs @@ -2,17 +2,7 @@ use precice; use std::env; use std::process::ExitCode; -fn solid_compute_solution(pressure: &[f64], cross_section_length: &mut [f64]) { - assert!(pressure.len() == cross_section_length.len()); - const E: f64 = 10000.0; - const C_MK2: f64 = E / std::f64::consts::FRAC_2_SQRT_PI; - const PRESSURE0: f64 = 0.0; - let new: Vec<_> = pressure - .iter() - .map(|p| ((PRESSURE0 - 2.0 * C_MK2) / (p - 2.0 * C_MK2)).powi(2)) - .collect(); - cross_section_length.copy_from_slice(&new[..]); -} +mod solver; fn main() -> ExitCode { println!("Starting Solid Solver..."); @@ -45,7 +35,7 @@ fn main() -> ExitCode { let grid_size = CHUNK_SIZE * dimensions as usize; let grid: Vec = { let mut v: Vec = vec![0_f64; grid_size]; - const DX : f64 = TUBE_LENGTH / DOMAIN_SIZE as f64; + const DX: f64 = TUBE_LENGTH / DOMAIN_SIZE as f64; for i in 0..CHUNK_SIZE - 1 { v[i * dimensions as usize] = DX * i as f64; } @@ -87,7 +77,7 @@ fn main() -> ExitCode { &mut pressure[..], ); - solid_compute_solution(&pressure, &mut cross_section_length); + solver::solid_compute_solution(&pressure, &mut cross_section_length); participant.write_data( mesh_name, diff --git a/elastic-tube-1d/solid-rust/src/solver.rs b/elastic-tube-1d/solid-rust/src/solver.rs new file mode 100644 index 000000000..2a37a227c --- /dev/null +++ b/elastic-tube-1d/solid-rust/src/solver.rs @@ -0,0 +1,11 @@ +pub fn solid_compute_solution(pressure: &[f64], cross_section_length: &mut [f64]) { + assert!(pressure.len() == cross_section_length.len()); + const E: f64 = 10000.0; + const C_MK2: f64 = E / std::f64::consts::FRAC_2_SQRT_PI; + const PRESSURE0: f64 = 0.0; + let new: Vec<_> = pressure + .iter() + .map(|p| ((PRESSURE0 - 2.0 * C_MK2) / (p - 2.0 * C_MK2)).powi(2)) + .collect(); + cross_section_length.copy_from_slice(&new[..]); +} From f81ee9ac58a22606bb542dcc3c03d3f046879336 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20Simonis?= Date: Thu, 15 Feb 2024 11:08:06 +0100 Subject: [PATCH 11/11] Apply cargo linting hints --- elastic-tube-1d/fluid-rust/src/main.rs | 13 ++++++------- elastic-tube-1d/fluid-rust/src/solver.rs | 2 +- elastic-tube-1d/fluid-rust/src/utils.rs | 22 +++++++++++----------- elastic-tube-1d/solid-rust/src/main.rs | 9 ++++----- 4 files changed, 22 insertions(+), 24 deletions(-) diff --git a/elastic-tube-1d/fluid-rust/src/main.rs b/elastic-tube-1d/fluid-rust/src/main.rs index c4d83b3fa..f9d472f4a 100644 --- a/elastic-tube-1d/fluid-rust/src/main.rs +++ b/elastic-tube-1d/fluid-rust/src/main.rs @@ -1,4 +1,3 @@ -use precice; use std::env; use std::process::ExitCode; @@ -20,7 +19,7 @@ fn main() -> ExitCode { const DOMAIN_SIZE: usize = 100; const CHUNK_SIZE: usize = DOMAIN_SIZE + 1; - let mut participant = precice::Participant::new("Fluid", &config, 0, 1); + let mut participant = precice::Participant::new("Fluid", config, 0, 1); println!("preCICE configured..."); @@ -42,9 +41,9 @@ fn main() -> ExitCode { const P0: f64 = 0.0; let vel_in0: f64 = U0 + AMPL * (FREQUENCY * T_SHIFT * std::f64::consts::PI).sin(); - let mut pressure: Vec = vec![P0; CHUNK_SIZE as usize]; - let mut cross_section_length: Vec = vec![A0; CHUNK_SIZE as usize]; - let mut velocity: Vec = vec![vel_in0; CHUNK_SIZE as usize]; + let mut pressure: Vec = vec![P0; CHUNK_SIZE]; + let mut cross_section_length: Vec = vec![A0; CHUNK_SIZE]; + let mut velocity: Vec = vec![vel_in0; CHUNK_SIZE]; let mut pressure_old = pressure.clone(); @@ -54,7 +53,7 @@ fn main() -> ExitCode { let mut v: Vec = vec![0_f64; grid_size]; for i in 0..CHUNK_SIZE { let idx = i * dimensions as usize; - v[idx] = i as f64 * CELLWIDTH as f64; + v[idx] = i as f64 * CELLWIDTH; } v }; @@ -156,5 +155,5 @@ fn main() -> ExitCode { println!("Exiting FluidSolver at t={}", t); participant.finalize(); - return ExitCode::SUCCESS; + ExitCode::SUCCESS } diff --git a/elastic-tube-1d/fluid-rust/src/solver.rs b/elastic-tube-1d/fluid-rust/src/solver.rs index 8f57cc3f2..ed5f4534b 100644 --- a/elastic-tube-1d/fluid-rust/src/solver.rs +++ b/elastic-tube-1d/fluid-rust/src/solver.rs @@ -13,7 +13,7 @@ pub fn fluid_compute_solution( pressure: &mut [f64], ) { // Initial guess - pressure.copy_from_slice(&pressure_old[..]); + pressure.copy_from_slice(pressure_old); const E: f64 = 10000_f64; //const C_MK2 : f64 = E / std::f64::consts::FRAC_2_SQRT_PI; diff --git a/elastic-tube-1d/fluid-rust/src/utils.rs b/elastic-tube-1d/fluid-rust/src/utils.rs index f593a6ee4..3ad0c8070 100644 --- a/elastic-tube-1d/fluid-rust/src/utils.rs +++ b/elastic-tube-1d/fluid-rust/src/utils.rs @@ -25,9 +25,9 @@ pub fn write_vtk( write!(f, "DATASET UNSTRUCTURED_GRID\n\n")?; write!(f, "POINTS {n_points} float\n\n")?; for i in 0..n_points { - write!( + writeln!( f, - "{:.16e} 0.0000000000000000e+00 0.0000000000000000e+00\n", + "{:.16e} 0.0000000000000000e+00 0.0000000000000000e+00", i as f64 * dx )?; } @@ -35,27 +35,27 @@ pub fn write_vtk( write!(f, "POINT_DATA {n_points}\n\n")?; - write!(f, "VECTORS velocity float\n")?; + writeln!(f, "VECTORS velocity float")?; for v in velocity { - write!( + writeln!( f, - "{:.16e} 0.0000000000000000e+00 0.0000000000000000e+00\n", + "{:.16e} 0.0000000000000000e+00 0.0000000000000000e+00", v )?; } writeln!(f)?; - write!(f, "SCALARS pressure float\n")?; - write!(f, "LOOKUP_TABLE default\n")?; + writeln!(f, "SCALARS pressure float")?; + writeln!(f, "LOOKUP_TABLE default")?; for v in pressure { - write!(f, "{:.16e}\n", v)?; + writeln!(f, "{:.16e}", v)?; } writeln!(f)?; - write!(f, "SCALARS diameter float\n")?; - write!(f, "LOOKUP_TABLE default\n")?; + writeln!(f, "SCALARS diameter float")?; + writeln!(f, "LOOKUP_TABLE default")?; for v in cross_section_length { - write!(f, "{:.16e}\n", v)?; + writeln!(f, "{:.16e}", v)?; } writeln!(f)?; Ok(()) diff --git a/elastic-tube-1d/solid-rust/src/main.rs b/elastic-tube-1d/solid-rust/src/main.rs index c7a2417fd..cf091a832 100644 --- a/elastic-tube-1d/solid-rust/src/main.rs +++ b/elastic-tube-1d/solid-rust/src/main.rs @@ -1,4 +1,3 @@ -use precice; use std::env; use std::process::ExitCode; @@ -20,7 +19,7 @@ fn main() -> ExitCode { const CHUNK_SIZE: usize = DOMAIN_SIZE + 1; const TUBE_LENGTH: f64 = 10.0; - let mut participant = precice::Participant::new("Solid", &config, 0, 1); + let mut participant = precice::Participant::new("Solid", config, 0, 1); println!("preCICE configured..."); @@ -29,8 +28,8 @@ fn main() -> ExitCode { assert!(participant.get_data_dimensions(mesh_name, "CrossSectionLength") == 1); assert!(participant.get_data_dimensions(mesh_name, "Pressure") == 1); - let mut pressure: Vec = vec![0.0; CHUNK_SIZE as usize]; - let mut cross_section_length: Vec = vec![1.0; CHUNK_SIZE as usize]; + let mut pressure: Vec = vec![0.0; CHUNK_SIZE]; + let mut cross_section_length: Vec = vec![1.0; CHUNK_SIZE]; let grid_size = CHUNK_SIZE * dimensions as usize; let grid: Vec = { @@ -99,5 +98,5 @@ fn main() -> ExitCode { println!("Exiting SolidSolver at t={}", t); participant.finalize(); - return ExitCode::SUCCESS; + ExitCode::SUCCESS }