diff --git a/.github/workflows/rustbca_compile_check.yml b/.github/workflows/rustbca_compile_check.yml index 584476c..382c9b8 100644 --- a/.github/workflows/rustbca_compile_check.yml +++ b/.github/workflows/rustbca_compile_check.yml @@ -59,12 +59,12 @@ jobs: cargo run --release 0D examples/boron_nitride_0D.toml ./target/release/RustBCA 0D examples/titanium_dioxide_0D.toml ./target/release/RustBCA 1D examples/layered_geometry_1D.toml - cat 2000.0eV_0.0001deg_He_TiO2_Al_Sisummary.output + echo "layered geometry 1D:"; cat 2000.0eV_0.0001deg_He_TiO2_Al_Sisummary.output ./target/release/RustBCA examples/layered_geometry.toml - cat 2000.0eV_0.0001deg_He_TiO2_Al_Sisummary.output + echo "layered geometry 2D:"; cat 2000.0eV_0.0001deg_He_TiO2_Al_Sisummary.output ./target/release/RustBCA SPHERE examples/boron_nitride_sphere.toml cargo run --release --features parry3d TRIMESH examples/tungsten_twist_trimesh.toml ./target/release/RustBCA examples/boron_nitride_wire.toml - cat boron_nitride_summary.output + echo "boron nitride 2D:"; cat boron_nitride_summary.output ./target/release/RustBCA HOMOGENEOUS2D examples/boron_nitride_wire_homogeneous.toml - cat boron_nitride_h_summary.output + echo "boron nitride 2D Homogeneous:"; cat boron_nitride_h_summary.output diff --git a/Cargo.toml b/Cargo.toml index a335ee2..dd45c67 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -5,6 +5,8 @@ default-run = "RustBCA" authors = ["Jon Drobny ", "Jon Drobny "] edition = "2021" + + [[bin]] name = "RustBCA" path = "src/main.rs" @@ -27,7 +29,8 @@ serde = { version = "1.0.163", features = ["derive"] } hdf5 = {version = "0.7.1", optional = true} rcpr = { git = "https://github.com/drobnyjt/rcpr", optional = true} ndarray = {version = "0.14.0", features = ["serde"], optional = true} -parry3d-f64 = {optional = true, version="0.2.0"} +parry3d-f64 = {optional = true, version="0.23.0"} +parry2d-f64 = {optional = false, version="0.23.0"} [dependencies.pyo3] version = "0.19.0" @@ -41,7 +44,7 @@ float-cmp = "0.8.0" lto = "fat" codegen-units = 1 opt-level = 3 -debug = false +debug = 1 [features] hdf5_input = ["hdf5"] diff --git a/README.md b/README.md index 8006e2c..7315a93 100644 --- a/README.md +++ b/README.md @@ -1,289 +1,289 @@ -# `RustBCA` - -`RustBCA` is a general-purpose, high-performance code for simulating -ion-material interactions including sputtering, reflection, and implantation -using the binary collision approximation ([BCA]), written in [Rust]! -RustBCA includes a standalone version and libraries for including -ion-material interactions in simulations written in C/C++, Python, -and Fortran. - -By discretizing the collision cascade into a sequence of binary collisions, -[BCA] codes can accurately and efficiently model the prompt interaction -between an energetic ion and a target material. This includes reflection, -implantation, and transmission of the incident ion, as well as sputtering -and displacement damage of the target. Generally, [BCA] codes can be -valid for incident ion energies between several eV/nucleon to -<1 GeV/nucleon. Improvements to RustBCA have expanded the regime -of validity for some quantities, such as reflection coefficients, below -1 eV/nucleon for some ion/target pairs. - -Check out the `RustBCA` [Wiki] for detailed information, installation -instructions, use cases, examples, and more. See the RustBCA paper at the -Journal of Open Source Software by clicking the badge below: - -[![DOI](https://joss.theoj.org/papers/10.21105/joss.03298/status.svg)](https://doi.org/10.21105/joss.03298) - -Selected citations of RustBCA as of 5/24/23: -* [Simulation of liquid lithium divertor geometry using SOLPS-ITER](https://doi.org/10.1109/TPS.2022.3166402), JD Lore et al. (2022) -* [Characterizing W sources in the all-W wall, all-RF WEST tokamak environment](https://doi.org/10.1088/1361-6587/ac8acc), CC Klepper et al. (2022) -* [hPIC2: A hardware-accelerated, hybrid particle-in-cell code for dynamic plasma-material interactions](https://doi.org/10.1016/j.cpc.2022.108569), LT Meredith et al. (2023) -* [Global sensitivity analysis of a coupled multiphysics model to predict surface evolution in fusion plasma–surface interactions](https://doi.org/10.1016/j.commatsci.2023.112229), P. Robbe et al. (2023) -* [Modeling the effect of nitrogen recycling on the erosion and leakage of tungsten impurities from the SAS-VW divertor in DIII-D during nitrogen gas injection](https://doi.org/10.1016/j.nme.2022.101254), MS Parsons et al. (2023) -* [Enabling attractive-repulsive potentials in binary-collision-approximation monte-carlo codes for ion-surface interactions](https://doi.org/10.1088/2053-1591/ad1262), J Drobny and D Curreli (2023) -* [Multi-physics modeling of tungsten collector probe samples during the WEST C4 He campaign](https://doi.org.10.1088/1741-4326/ad6c5b), A. Lasa et al. (2024) - -## Getting started - -The easiest way to get started is with the ergonomic Python functions. -These functions use the default RustBCA options detailed on the -[Input Files](https://github.com/lcpp-org/RustBCA/wiki/Standalone-Code:-Input-File) page, which are not universally applicable. -These examples use example material parameters located in -`scripts/materials.py` that should be verified before use. -Follow these steps to install, build, and run simple RustBCA simulations -for sputtering yields and reflection coefficients: -``` -git clone https://github.com/lcpp-org/rustbca -cd rustbca -python -m pip install . -python -Python 3.9.6 (tags/v3.9.6:db3ff76, Jun 28 2021, 15:26:21) [MSC v.1929 64 bit (AMD64)] on win32 -Type "help", "copyright", "credits" or "license" for more information. ->>> from libRustBCA import *; from scripts.materials import *; import numpy as np ->>> angle = 0.0 # deg ->>> energy = 1000.0 # eV ->>> num_samples = 10000 ->>> 1 < sputtering_yield(argon, tungsten, energy, angle, num_samples) < 1.1 # Y approx. 1.04 -True ->>> R_N, R_E = reflection_coefficient(argon, tungsten, energy, angle, num_samples) ->>> 0.3 < R_N < 0.4 # R_N approx. 0.35 -True ->>> 0.0 < R_E < 0.2 # R_E approx 0.1 -True -``` - -For those eager to get started with the standalone code, try running one of the examples in the -`RustBCA/examples` directory. Note that to automatically manipulate input files and reproduce -the plots located on the [Wiki], these may require some optional -[Python] packages (`matplotlib`, `numpy`, `scipy`, `shapely`, and `toml`). - -### H trajectories and collision cascades in boron nitride -First, run the example using: - -```shell -$ cargo run --release 0D examples/boron_nitride_0D.toml 2>/dev/null #suppress progress bar for automatic testing -Processing 10 ions... -Initializing with 4 threads... -Finished! -``` - -Afterwords, fire up your favourite [Python] interpreter -(e.g., [IPython]) and execute: - -```python -from scripts.rustbca import * -do_trajectory_plot("boron_nitride_") -``` - -### He implantation into a layered TiO2/Al/Si target - -First, run the example using: - -```shell -$ cargo run --release examples/layered_geometry.toml 2>/dev/null #suppress progress bar for automatic testing -Processing 10000 ions... -Initializing with 4 threads... -Finished! -``` - -Afterwords, fire up your favourite [Python] interpreter -(e.g., [IPython]) and execute: - -```python -import numpy as np -import matplotlib.pyplot as plt - -deposited_ions = np.genfromtxt( - "2000.0eV_0.0001deg_He_TiO2_Al_Sideposited.output", - delimiter=",", - names=["M", "Z", "x", "y", "z", "collisions"], -) - -plt.hist(deposited_ions["x"], bins=100) - -plt.show() -``` - -## Features - -The following features are implemented in `RustBCA`: - -* Ion-material interactions for all combinations of incident ion and target species. -* Infinite, homogeneous targets (Mesh0D), Layered, finite-depth inhomogeneous targets (Mesh1D), arbitrary 2D composition through a triangular mesh (Mesh2D), fast homogeneous 2D geometry (Homogeneous2D), homogeneous spherical geometry (Sphere), and homogeneous 3D triangular mesh geometry (TriMesh). -* Amorphous Solid/Liquid targets, Gaseous targets, and targets with both solid/liquid and gaseous elements -* Low energy (< 25 keV/nucleon) electronic stopping modes including: - * local (Oen-Robinson), - * nonlocal (Lindhard-Scharff), - * and equipartition -* Biersack-Varelas interpolation is also included for electronic stopping up to ~1 GeV/nucleon. Note that high energy physics beyond electronic stopping are not included, and that Biersack-Varelas may not be as accurate as other methods. -* Biersack-Haggmark treatment of high-energy free-flight paths between collisions can be included to greatly speed up high-energy simulations (i.e., by neglecting very small angle scattering). -* A wide range of interaction potentials are provided, including: - * the Kr-C, ZBL, Lenz-Jensen, and Moliere universal, screened-Coulomb potentials. - * the Lennard-Jones 12-6 and Morse attractive-repulsive potentials. -* Solving the distance-of-closest-approach problem is achieved using: - * the Newton-Raphson method for purely repulsive potentials, - * or, for attractive-repulsive potentials, an Adaptive Chebyshev Proxy Rootfinder with Automatic Subdivision algorithm and a polynomial root-finding algorithm are provided through [rcpr]. -* Multiple interaction potentials can be used in a single simulation for any number of potentials/species. - * For example, the He-W interaction can be specified using a Lennard-Jones 12-6 potential, while the W-W interaction can be defined using a Kr-C potential. -* The scattering integral can be calculated using: - * Gauss-Mehler quadrature, - * Gauss-Legendre quadrature, - * Mendenall-Weller quadrature, - * or the MAGIC algorithm (for certain screened Coulomb potentials only). -* Input files use the [TOML] format, making them both human-readable and easily parsable. -* RustBCA generates user-friendly, context-providing error messages, which help pinpoint the cause of errors and provide suggested fixes to the user. -* The simulation results are comma-delimited (`csv` format) and include: - * the energies and directions of emitted particles (reflected ions and sputtered atoms), - * the final positions of implanted ions, - * displacements, - * full trajectory tracking for both the incident ions and target atoms, - * and many other parameters such as position of origin of sputtered particles and energy loss along trajectories. -* Optionally, the code can produce energy-angle and implantation distributions when built with the `--features distributions` flag and disable space-intensive particle list output with `--features no_list_output`. -* Library functions for modeling ion reflection, implantation, and sputtering in C++/C, Python, and Fortran codes. - -## Installation - -Without optional features, `RustBCA` should compile with `cargo` alone on -Windows, MacOS, and Linux systems: - -``` -cargo build --release -``` - -will add an executable at `target/release/`. - -[HDF5] for particle list input has been tested on Windows, but version 1.10.6 must be used. - -#### Manual Dependences - -* [rustup], the [Rust] toolchain (includes `cargo`, the [Rust] package manager, `rustc`, the [Rust] compiler, and more). - -#### Automatic Dependencies - -* see [Cargo.toml](https://github.com/lcpp-org/RustBCA/blob/master/Cargo.toml) for a complete list of required and optional dependencies managed by `cargo`. - -#### Optional Dependencies - -* [HDF5] libraries -* For manipulating input files and running associated scripts, the following are suggested: - * [Python] 3.6+ - * [Python] libraries: `numpy`, `matplotlib`, `toml`, `shapely`, and `scipy`. - -### Detailed instructions for Ubuntu 18.04 LTS - -1. (Optional) Install Python 3.6+ (this comes natively in Ubuntu 18.04) -2. Install `curl`: -```bash -sudo apt-get install curl -``` -3. Install [rustup], the Rust toolchain (includes rustc, the compiler, and cargo, the package manager) from https://rustup.rs/ by running the following command and following on-screen instructions: -```bash -curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y -``` -4. (Optional) Install `pip` for [Python]: -```bash -sudo apt-get install python3-pip -``` -5. (Optional) Install [Python] libraries for making input files: -```bash -python3 -m pip install numpy matplotlib shapely scipy toml -``` -6. (Optional - should come with rustup) Install `cargo`: -```bash -sudo apt-get install cargo -``` -7. Build `RustBCA`: -```bash -git clone https://github.com/lcpp-org/RustBCA -cd RustBCA -cargo build --release -``` -8. (Optional) Build `RustBCA` with optional dependencies, `hdf5` and/or `rcpr`: -```bash -cargo build --release --features cpr_rootfinder,hdf5 -``` -9. `input.toml` is the input file - see the [Input File](https://github.com/lcpp-org/RustBCA/wiki/Standalone-Code:-Input-File) page for more information -10. Run the required tests using: -```bash -cargo test -``` -11. (Optional) Run the tests for the advanced rootfinder: -```bash -cargo test --features cpr_rootfinder -``` - -### Detailed instructions for Fedora 33 - -Most of the ingredients for building `RustBCA` and running the [Python] helper -scripts are available natively in the Fedora software repository, so the setup -is relatively painless. - -The [Rust] toolchain can be aquired using: - -```bash -sudo dnf install rust rust-src rust-std-static rust-analysis rust-gdb rust-lldb rustfmt -``` - -The (optional) [Python] packages can be obtained using: - -```bash -sudo dnf install python3-numpy python3-scipy python3-matplotlib python3-toml python3-shapely -``` - -or, alternatively, using `pip3`. - -Building `RustBCA` is straightforward, and can be done using: - -```bash -git clone https://github.com/lcpp-org/rustBCA -cd RustBCA -cargo build --release -``` - -with all of the explicit dependencies listed in `Cargo.toml` handled -automatically during the build. - -## Usage - -To use `RustBCA`, modify an `input.toml` file, which is used to configure each -simulation. -To run a simulation, execute: - -```bash -./RustBCA -``` - -with `input.toml` in the same directory as `RustBCA`. - -Alternatively, `RustBCA` accepts the name of a`.toml` input file as a single -command line argument. - -```bash -./RustBCA /path/to/input.toml -``` - -Additionally, `RustBCA` accepts an input file type (one of: `0D`, `1D`, `2D`, `TRIMESH`, `SPHERE` - see the wiki for more details): -```bash -./RustBCA 0D /path/to/input.toml -``` -**Warning: RustBCA defaults to the 2D triangular mesh input mode.** For more details, see [Input Files](https://github.com/lcpp-org/RustBCA/wiki/Standalone-Code:-Input-File). -Also have a look at the examples on the [Wiki] to see some examples of RustBCA input files. - -[BCA]: https://en.wikipedia.org/wiki/Binary_collision_approximation -[HDF5]: https://en.wikipedia.org/wiki/Hierarchical_Data_Format -[IPython]: https://en.wikipedia.org/wiki/IPython -[Python]: https://en.wikipedia.org/wiki/Python_(programming_language) -[rcpr]: https://github.com/drobnyjt/rcpr -[rustup]: https://rustup.rs -[Rust]: https://en.wikipedia.org/wiki/Rust_(programming_language) -[TOML]: https://en.wikipedia.org/wiki/TOML -[Wiki]: https://github.com/lcpp-org/RustBCA/wiki +# `RustBCA` + +`RustBCA` is a general-purpose, high-performance code for simulating +ion-material interactions including sputtering, reflection, and implantation +using the binary collision approximation ([BCA]), written in [Rust]! +RustBCA includes a standalone version and libraries for including +ion-material interactions in simulations written in C/C++, Python, +and Fortran. + +By discretizing the collision cascade into a sequence of binary collisions, +[BCA] codes can accurately and efficiently model the prompt interaction +between an energetic ion and a target material. This includes reflection, +implantation, and transmission of the incident ion, as well as sputtering +and displacement damage of the target. Generally, [BCA] codes can be +valid for incident ion energies between several eV/nucleon to +<1 GeV/nucleon. Improvements to RustBCA have expanded the regime +of validity for some quantities, such as reflection coefficients, below +1 eV/nucleon for some ion/target pairs. + +Check out the `RustBCA` [Wiki] for detailed information, installation +instructions, use cases, examples, and more. See the RustBCA paper at the +Journal of Open Source Software by clicking the badge below: + +[![DOI](https://joss.theoj.org/papers/10.21105/joss.03298/status.svg)](https://doi.org/10.21105/joss.03298) + +Selected citations of RustBCA as of 5/24/23: +* [Simulation of liquid lithium divertor geometry using SOLPS-ITER](https://doi.org/10.1109/TPS.2022.3166402), JD Lore et al. (2022) +* [Characterizing W sources in the all-W wall, all-RF WEST tokamak environment](https://doi.org/10.1088/1361-6587/ac8acc), CC Klepper et al. (2022) +* [hPIC2: A hardware-accelerated, hybrid particle-in-cell code for dynamic plasma-material interactions](https://doi.org/10.1016/j.cpc.2022.108569), LT Meredith et al. (2023) +* [Global sensitivity analysis of a coupled multiphysics model to predict surface evolution in fusion plasma–surface interactions](https://doi.org/10.1016/j.commatsci.2023.112229), P. Robbe et al. (2023) +* [Modeling the effect of nitrogen recycling on the erosion and leakage of tungsten impurities from the SAS-VW divertor in DIII-D during nitrogen gas injection](https://doi.org/10.1016/j.nme.2022.101254), MS Parsons et al. (2023) +* [Enabling attractive-repulsive potentials in binary-collision-approximation monte-carlo codes for ion-surface interactions](https://doi.org/10.1088/2053-1591/ad1262), J Drobny and D Curreli (2023) +* [Multi-physics modeling of tungsten collector probe samples during the WEST C4 He campaign](https://doi.org.10.1088/1741-4326/ad6c5b), A. Lasa et al. (2024) + +## Getting started + +The easiest way to get started is with the ergonomic Python functions. +These functions use the default RustBCA options detailed on the +[Input Files](https://github.com/lcpp-org/RustBCA/wiki/Standalone-Code:-Input-File) page, which are not universally applicable. +These examples use example material parameters located in +`scripts/materials.py` that should be verified before use. +Follow these steps to install, build, and run simple RustBCA simulations +for sputtering yields and reflection coefficients: +``` +git clone https://github.com/lcpp-org/rustbca +cd rustbca +python -m pip install . +python +Python 3.9.6 (tags/v3.9.6:db3ff76, Jun 28 2021, 15:26:21) [MSC v.1929 64 bit (AMD64)] on win32 +Type "help", "copyright", "credits" or "license" for more information. +>>> from libRustBCA import *; from scripts.materials import *; import numpy as np +>>> angle = 0.0 # deg +>>> energy = 1000.0 # eV +>>> num_samples = 10000 +>>> 1 < sputtering_yield(argon, tungsten, energy, angle, num_samples) < 1.1 # Y approx. 1.04 +True +>>> R_N, R_E = reflection_coefficient(argon, tungsten, energy, angle, num_samples) +>>> 0.3 < R_N < 0.4 # R_N approx. 0.35 +True +>>> 0.0 < R_E < 0.2 # R_E approx 0.1 +True +``` + +For those eager to get started with the standalone code, try running one of the examples in the +`RustBCA/examples` directory. Note that to automatically manipulate input files and reproduce +the plots located on the [Wiki], these may require some optional +[Python] packages (`matplotlib`, `numpy`, `scipy`, `shapely`, and `toml`). + +### H trajectories and collision cascades in boron nitride +First, run the example using: + +```shell +$ cargo run --release 0D examples/boron_nitride_0D.toml 2>/dev/null #suppress progress bar for automatic testing +Processing 10 ions... +Initializing with 4 threads... +Finished! +``` + +Afterwords, fire up your favourite [Python] interpreter +(e.g., [IPython]) and execute: + +```python +from scripts.rustbca import * +do_trajectory_plot("boron_nitride_") +``` + +### He implantation into a layered TiO2/Al/Si target + +First, run the example using: + +```shell +$ cargo run --release examples/layered_geometry.toml 2>/dev/null #suppress progress bar for automatic testing +Processing 10000 ions... +Initializing with 4 threads... +Finished! +``` + +Afterwords, fire up your favourite [Python] interpreter +(e.g., [IPython]) and execute: + +```python +import numpy as np +import matplotlib.pyplot as plt + +deposited_ions = np.genfromtxt( + "2000.0eV_0.0001deg_He_TiO2_Al_Sideposited.output", + delimiter=",", + names=["M", "Z", "x", "y", "z", "collisions"], +) + +plt.hist(deposited_ions["x"], bins=100) + +plt.show() +``` + +## Features + +The following features are implemented in `RustBCA`: + +* Ion-material interactions for all combinations of incident ion and target species. +* Infinite, homogeneous targets (Mesh0D), Layered, finite-depth inhomogeneous targets (Mesh1D), arbitrary 2D composition through a triangular mesh (Mesh2D), fast homogeneous 2D geometry (Homogeneous2D), homogeneous spherical geometry (Sphere), and homogeneous 3D triangular mesh geometry (TriMesh). +* Amorphous Solid/Liquid targets, Gaseous targets, and targets with both solid/liquid and gaseous elements +* Low energy (< 25 keV/nucleon) electronic stopping modes including: + * local (Oen-Robinson), + * nonlocal (Lindhard-Scharff), + * and equipartition +* Biersack-Varelas interpolation is also included for electronic stopping up to ~1 GeV/nucleon. Note that high energy physics beyond electronic stopping are not included, and that Biersack-Varelas may not be as accurate as other methods. +* Biersack-Haggmark treatment of high-energy free-flight paths between collisions can be included to greatly speed up high-energy simulations (i.e., by neglecting very small angle scattering). +* A wide range of interaction potentials are provided, including: + * the Kr-C, ZBL, Lenz-Jensen, and Moliere universal, screened-Coulomb potentials. + * the Lennard-Jones 12-6 and Morse attractive-repulsive potentials. +* Solving the distance-of-closest-approach problem is achieved using: + * the Newton-Raphson method for purely repulsive potentials, + * or, for attractive-repulsive potentials, an Adaptive Chebyshev Proxy Rootfinder with Automatic Subdivision algorithm and a polynomial root-finding algorithm are provided through [rcpr]. +* Multiple interaction potentials can be used in a single simulation for any number of potentials/species. + * For example, the He-W interaction can be specified using a Lennard-Jones 12-6 potential, while the W-W interaction can be defined using a Kr-C potential. +* The scattering integral can be calculated using: + * Gauss-Mehler quadrature, + * Gauss-Legendre quadrature, + * Mendenall-Weller quadrature, + * or the MAGIC algorithm (for certain screened Coulomb potentials only). +* Input files use the [TOML] format, making them both human-readable and easily parsable. +* RustBCA generates user-friendly, context-providing error messages, which help pinpoint the cause of errors and provide suggested fixes to the user. +* The simulation results are comma-delimited (`csv` format) and include: + * the energies and directions of emitted particles (reflected ions and sputtered atoms), + * the final positions of implanted ions, + * displacements, + * full trajectory tracking for both the incident ions and target atoms, + * and many other parameters such as position of origin of sputtered particles and energy loss along trajectories. +* Optionally, the code can produce energy-angle and implantation distributions when built with the `--features distributions` flag and disable space-intensive particle list output with `--features no_list_output`. +* Library functions for modeling ion reflection, implantation, and sputtering in C++/C, Python, and Fortran codes. + +## Installation + +Without optional features, `RustBCA` should compile with `cargo` alone on +Windows, MacOS, and Linux systems: + +``` +cargo build --release +``` + +will add an executable at `target/release/`. + +[HDF5] for particle list input has been tested on Windows, but version 1.10.6 must be used. + +#### Manual Dependences + +* [rustup], the [Rust] toolchain (includes `cargo`, the [Rust] package manager, `rustc`, the [Rust] compiler, and more). + +#### Automatic Dependencies + +* see [Cargo.toml](https://github.com/lcpp-org/RustBCA/blob/master/Cargo.toml) for a complete list of required and optional dependencies managed by `cargo`. + +#### Optional Dependencies + +* [HDF5] libraries +* For manipulating input files and running associated scripts, the following are suggested: + * [Python] 3.6+ + * [Python] libraries: `numpy`, `matplotlib`, `toml`, `shapely`, and `scipy`. + +### Detailed instructions for Ubuntu 18.04 LTS + +1. (Optional) Install Python 3.6+ (this comes natively in Ubuntu 18.04) +2. Install `curl`: +```bash +sudo apt-get install curl +``` +3. Install [rustup], the Rust toolchain (includes rustc, the compiler, and cargo, the package manager) from https://rustup.rs/ by running the following command and following on-screen instructions: +```bash +curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y +``` +4. (Optional) Install `pip` for [Python]: +```bash +sudo apt-get install python3-pip +``` +5. (Optional) Install [Python] libraries for making input files: +```bash +python3 -m pip install numpy matplotlib shapely scipy toml +``` +7. (Optional - should come with rustup) Install `cargo`: +```bash +sudo apt-get install cargo +``` +8. Build `RustBCA`: +```bash +git clone https://github.com/lcpp-org/RustBCA +cd RustBCA +cargo build --release +``` +9. (Optional) Build `RustBCA` with optional dependencies, `hdf5` and/or `rcpr`: +```bash +cargo build --release --features cpr_rootfinder,hdf5_input + ``` +10. `input.toml` is the input file - see [Usage](https://github.com/lcpp-org/RustBCA/wiki/Usage,-Input-File,-and-Output-Files) for more information +11. Run the required tests using: +```bash +cargo test +``` +12. (Optional) Run the required and optional tests for the desired backend(s): +```bash +cargo test --features cpr_rootfinder +``` + +### Detailed instructions for Fedora 33 + +Most of the ingredients for building `RustBCA` and running the [Python] helper +scripts are available natively in the Fedora software repository, so the setup +is relatively painless. + +The [Rust] toolchain can be aquired using: + +```bash +sudo dnf install rust rust-src rust-std-static rust-analysis rust-gdb rust-lldb rustfmt +``` + +The (optional) [Python] packages can be obtained using: + +```bash +sudo dnf install python3-numpy python3-scipy python3-matplotlib python3-toml python3-shapely +``` + +or, alternatively, using `pip3`. + +Building `RustBCA` is straightforward, and can be done using: + +```bash +git clone https://github.com/lcpp-org/rustBCA +cd RustBCA +cargo build --release +``` + +with all of the explicit dependencies listed in `Cargo.toml` handled +automatically during the build. + +## Usage + +To use `RustBCA`, modify an `input.toml` file, which is used to configure each +simulation. +To run a simulation, execute: + +```bash +./RustBCA +``` + +with `input.toml` in the same directory as `RustBCA`. + +Alternatively, `RustBCA` accepts the name of a`.toml` input file as a single +command line argument. + +```bash +./RustBCA /path/to/input.toml +``` + +Additionally, `RustBCA` accepts an input file type (one of: `0D`, `1D`, `2D`, `TRIMESH`, `SPHERE`, `HOMOGENEOUS2D` - see the wiki for more details): +```bash +./RustBCA 0D /path/to/input.toml +``` +**Warning: RustBCA defaults to the 2D triangular mesh input mode.** For more details, see [Input Files](https://github.com/lcpp-org/RustBCA/wiki/Standalone-Code:-Input-File). +Also have a look at the examples on the [Wiki] to see some examples of RustBCA input files. + +[BCA]: https://en.wikipedia.org/wiki/Binary_collision_approximation +[HDF5]: https://en.wikipedia.org/wiki/Hierarchical_Data_Format +[IPython]: https://en.wikipedia.org/wiki/IPython +[Python]: https://en.wikipedia.org/wiki/Python_(programming_language) +[rcpr]: https://github.com/drobnyjt/rcpr +[rustup]: https://rustup.rs +[Rust]: https://en.wikipedia.org/wiki/Rust_(programming_language) +[TOML]: https://en.wikipedia.org/wiki/TOML +[Wiki]: https://github.com/lcpp-org/RustBCA/wiki diff --git a/examples/layered_geometry.toml b/examples/layered_geometry.toml index aaaa434..c070292 100644 --- a/examples/layered_geometry.toml +++ b/examples/layered_geometry.toml @@ -1,7 +1,7 @@ #Use with 2D geometry option only [options] name = "2000.0eV_0.0001deg_He_TiO2_Al_Si" -track_trajectories = false +track_trajectories = true track_recoils = true track_recoil_trajectories = false write_buffer_size = 8000 @@ -54,7 +54,8 @@ points = [[0.0, -0.5], [0.01, -0.5], [0.04, -0.5], [0.5, -0.5], [0.5, 0.5], [0.0 triangles = [[0, 1, 6], [0, 6, 7], [1, 2, 5], [1, 5, 6], [2, 3, 4], [2, 4, 5]] densities = [ [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 0.0, 4.996e+10,], [ 0.0, 0.0, 0.0, 4.996e+10,],] boundary = [0, 3, 4, 7] -simulation_boundary_points = [ [ 0.6, -0.6,], [ -0.1, -0.6,], [ -0.1, 0.6,], [ 0.6, 0.6,], [ 0.6, -0.6,],] +#simulation_boundary_points = [ [ 0.6, -0.6,], [ -0.1, -0.6,], [ -0.1, 0.6,], [ 0.6, 0.6,],] +simulation_boundary_points = [[0.6, 0.6], [-0.1, 0.6], [-0.1, -0.6], [0.6, -0.6]] electronic_stopping_correction_factors = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0] energy_barrier_thickness = 2.2E-4 diff --git a/examples/multiple_interaction_potentials.toml b/examples/multiple_interaction_potentials.toml index ba1aea1..f129abe 100644 --- a/examples/multiple_interaction_potentials.toml +++ b/examples/multiple_interaction_potentials.toml @@ -48,7 +48,7 @@ length_unit = "MICRON" triangles = [ [ 0.0, 0.01, 0.0, 0.5, 0.5, -0.5,], [ 0.0, 0.01, 0.01, -0.5, 0.5, -0.5,], [ 0.01, 0.01, 0.04, -0.5, 0.5, -0.5,], [ 0.01, 0.04, 0.04, 0.5, 0.5, -0.5,], [ 0.04, 0.5, 0.04, 0.5, 0.5, -0.5,], [ 0.04, 0.5, 0.5, -0.5, 0.5, -0.5,],] densities = [ [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 0.0, 4.996e+10,], [ 0.0, 0.0, 0.0, 4.996e+10,],] material_boundary_points = [ [ 0.0, 0.5,], [ 0.5, 0.5,], [ 0.5, -0.5,], [ 0.0, -0.5,],] -simulation_boundary_points = [ [ 0.6, -0.6,], [ -0.1, -0.6,], [ -0.1, 0.6,], [ 0.6, 0.6,], [ 0.6, -0.6,],] +simulation_boundary_points = [ [ 0.6, -0.6,], [ -0.1, -0.6,], [ -0.1, 0.6,], [ 0.6, 0.6,],] electronic_stopping_correction_factors = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0] [material_parameters] diff --git a/src/enums.rs b/src/enums.rs index 908c040..d8582b6 100644 --- a/src/enums.rs +++ b/src/enums.rs @@ -3,13 +3,13 @@ use super::*; pub enum MaterialType { MESH0D(material::Material), MESH1D(material::Material), - MESH2D(material::Material), + MESH2D(material::Material), SPHERE(material::Material), #[cfg(feature = "parry3d")] BALL(material::Material), #[cfg(feature = "parry3d")] TRIMESH(material::Material), - HOMOGENEOUS2D(material::Material) + HOMOGENEOUS2D(material::Material) } #[derive(Deserialize, PartialEq, Clone, Copy)] diff --git a/src/geometry.rs b/src/geometry.rs index a915629..09f26dd 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -1,8 +1,11 @@ use super::*; -use geo::algorithm::contains::Contains; -use geo::{Polygon, LineString, Point, point, Closest}; -use geo::algorithm::closest_point::ClosestPoint; +use itertools::Itertools; + +use parry2d_f64::shape::{Polyline, TriMesh, FeatureId::*}; +use parry2d_f64::query::PointQuery; +use parry2d_f64::math::Point as Point2d; +use parry2d_f64::math::Isometry; ///Trait for a Geometry object - all forms of geometry must implement these traits to be used pub trait Geometry { @@ -18,6 +21,7 @@ pub trait Geometry { fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool; fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool; fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64); + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64); } pub trait GeometryElement: Clone { @@ -109,6 +113,10 @@ impl Geometry for Mesh0D { fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { (0., y, z) } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + (-1., 0., 0.) + } } #[derive(Deserialize, Clone)] @@ -265,6 +273,175 @@ impl Geometry for Mesh1D { (self.top, y, z) } } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + (-1., 0., 0.) + } +} + +#[derive(Deserialize, Clone)] +pub struct ParryHomogeneousMesh2DInput { + pub length_unit: String, + pub points: Vec<(f64, f64)>, + pub simulation_boundary_points: Vec<(f64, f64)>, + pub densities: Vec, + pub electronic_stopping_correction_factor: f64 +} + +#[derive(Clone)] +pub struct ParryHomogeneousMesh2D { + pub boundary: Polyline, + pub simulation_boundary: Polyline, + pub energy_barrier_thickness: f64, + pub densities: Vec, + pub electronic_stopping_correction_factor: f64, + pub concentrations: Vec +} + + +impl Geometry for ParryHomogeneousMesh2D { + + type InputFileFormat = InputHomogeneous2D; + + fn new(input: &<::InputFileFormat as GeometryInput>::GeometryInput) -> Self { + let length_unit: f64 = match input.length_unit.as_str() { + "MICRON" => MICRON, + "CM" => CM, + "MM" => MM, + "ANGSTROM" => ANGSTROM, + "NM" => NM, + "M" => 1., + _ => input.length_unit.parse() + .expect(format!( + "Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M", + &input.length_unit.as_str() + ).as_str()), + }; + + let mut boundary_points_converted: Vec> = input.points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect(); + let number_boundary_points = boundary_points_converted.len() as u32; + + for p in boundary_points_converted.iter().combinations(2) { + assert!(p[0] != p[1], "Input error: duplicate vertices in boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).") + } + + let mut simulation_boundary_points_converted: Vec> = input.simulation_boundary_points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect(); + let number_simulation_boundary_points = simulation_boundary_points_converted.len() as u32; + + let test_simulation_boundary_ccw = (0..number_simulation_boundary_points as usize) + .map(|i| (simulation_boundary_points_converted[(i + 1) % number_simulation_boundary_points as usize].x - simulation_boundary_points_converted[i].x)*(simulation_boundary_points_converted[i].y + simulation_boundary_points_converted[(i + 1) % number_simulation_boundary_points as usize].y)) + .sum::(); + + if test_simulation_boundary_ccw > 0.0 { + simulation_boundary_points_converted = simulation_boundary_points_converted.clone().into_iter().rev().collect::>>(); + } + + for p in simulation_boundary_points_converted.iter().combinations(2) { + assert!(p[0] != p[1], "Input error: duplicate vertices in simulation boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).") + } + + let test_ccw = (0..number_boundary_points as usize) + .map(|i| (boundary_points_converted[(i + 1) % number_boundary_points as usize].x - boundary_points_converted[i].x)*(boundary_points_converted[i].y + boundary_points_converted[(i + 1) % number_boundary_points as usize].y)) + .sum::(); + + if test_ccw > 0.0 { + boundary_points_converted = boundary_points_converted.clone().into_iter().rev().collect::>>(); + } + + let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor; + + let densities: Vec = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect(); + + let total_density: f64 = densities.iter().sum(); + + let energy_barrier_thickness = total_density.powf(-1./3.)/SQRTPI*2.; + + let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); + + let mut linked_boundary_points = (0..number_boundary_points).zip(1..number_boundary_points).map(|(x, y)| [x, y]).collect::>(); + linked_boundary_points.push([number_boundary_points - 1, 0]); + let boundary2 = Polyline::new(boundary_points_converted, Some(linked_boundary_points)); + + let mut linked_simulation_boundary_points = (0..number_simulation_boundary_points).zip(1..number_simulation_boundary_points).map(|(x, y)| [x, y]).collect::>(); + linked_simulation_boundary_points.push([number_simulation_boundary_points - 1, 0]); + let simulation_boundary2 = Polyline::new(simulation_boundary_points_converted, Some(linked_simulation_boundary_points)); + + ParryHomogeneousMesh2D { + densities, + simulation_boundary: simulation_boundary2, + boundary: boundary2, + electronic_stopping_correction_factor, + energy_barrier_thickness, + concentrations + } + } + + fn get_densities(&self, x: f64, y: f64, z: f64) -> &Vec { + &self.densities + } + + fn get_ck(&self, x: f64, y: f64, z: f64) -> f64 { + self.electronic_stopping_correction_factor + } + + fn get_total_density(&self, x: f64, y: f64, z: f64) -> f64 { + self.densities.iter().sum::() + } + + fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec { + &self.concentrations + } + + fn inside(&self, x: f64, y: f64, z: f64) -> bool { + let p = Point2d::new(x, y); + if self.boundary.aabb(&Isometry::identity()).contains_local_point(&p) { + let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); + point_projection.is_inside + } else { + false + } + } + + fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool { + let p = Point2d::new(x, y); + + let (point_projection, (_, _)) = self.simulation_boundary.project_local_point_assuming_solid_interior_ccw(p); + point_projection.is_inside + + } + fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { + if self.inside(x, y, z) { + true + } else { + let p = Point2d::new(x, y); + (self.boundary.distance_to_local_point(&p, true) as f64) < self.energy_barrier_thickness + } + } + + fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let p = Point2d::new(x, y); + let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); + let (x_, y_) = (point_projection.point.x, point_projection.point.y); + (x_ as f64, y_ as f64, z) + } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let p = Point2d::new(x, y); + + let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); + let (x_intersect, y_intersect, z_intersect) = self.closest_point(x, y, z); + + let dx = x_intersect - x; + let dy = y_intersect - y; + let dz = z_intersect - z; + let mag = (dx*dx + dy*dy + dz*dz).sqrt(); + + if point_projection.is_inside { + (dx/mag, dy/mag, dz/mag) + } else { + (-dx/mag, -dy/mag, -dz/mag) + } + } } #[derive(Deserialize, Clone)] @@ -276,6 +453,7 @@ pub struct HomogeneousMesh2DInput { pub electronic_stopping_correction_factor: f64 } +/* #[derive(Clone)] pub struct HomogeneousMesh2D { pub boundary: Polygon, @@ -379,7 +557,12 @@ impl Geometry for HomogeneousMesh2D { panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z); } } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + panic!("Not implemented.") + } } +*/ /// Object that contains raw mesh input data. #[derive(Deserialize, Clone)] @@ -394,6 +577,213 @@ pub struct Mesh2DInput { pub energy_barrier_thickness: f64, } +#[derive(Clone)] +pub struct ParryMesh2D { + trimesh: TriMesh, + densities: Vec>, + electronic_stopping_correction_factors: Vec, + concentrations: Vec>, + pub boundary: Polyline, + pub simulation_boundary: Polyline, + pub energy_barrier_thickness: f64 +} + +impl ParryMesh2D { + fn get_id(&self, x: f64, y: f64, z: f64) -> usize { + let p = Point2d::new(x, y); + let (point_projection, feature_id) = self.trimesh.project_local_point_and_get_feature(&p); + match feature_id { + Vertex(vertex_id) => { + for (triangle_id, triangle) in self.trimesh.triangles().enumerate() { + if triangle.contains_local_point(&self.trimesh.vertices()[vertex_id as usize]) { + return triangle_id as usize + } + } + panic!("Geometry error: point ({}, {}) located on a vertex that is not associated with any triangle. Check mesh input.", x, y); + }, + Face(triangle_id) => { + triangle_id as usize + } + Unknown => { + panic!("Geometry error: unknown feature detected near ({}, {}). Check mesh input.", x, y); + } + } + } +} + +impl Geometry for ParryMesh2D { + type InputFileFormat = Input2D; + + fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { + if self.inside(x, y, z) { + true + } else { + let p = Point2d::new(x, y); + if self.trimesh.distance_to_local_point(&p, true) < self.energy_barrier_thickness { + true + } else { + false + } + } + } + + fn inside_simulation_boundary (&self, x: f64, y: f64, z: f64) -> bool { + let p = Point2d::new(x, y); + + let (point_projection, (_, _)) = self.simulation_boundary.project_local_point_assuming_solid_interior_ccw(p); + + point_projection.is_inside + } + + fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let p = Point2d::new(x, y); + let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); + let (x_, y_) = (point_projection.point.x, point_projection.point.y); + (x_ as f64, y_ as f64, z) + } + + fn get_densities(&self, x: f64, y: f64, z: f64) -> &Vec { + &self.densities[self.get_id(x, y, z)] + } + + fn get_ck(&self, x: f64, y: f64, z: f64) -> f64 { + self.electronic_stopping_correction_factors[self.get_id(x, y, z)] + } + + /// Determine the total number density of the triangle that contains or is nearest to (x, y). + fn get_total_density(&self, x: f64, y: f64, z: f64) -> f64 { + self.get_densities(x, y, z).iter().sum::() + } + + /// Find the concentrations of the triangle that contains or is nearest to (x, y). + fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec { + &self.concentrations[self.get_id(x, y, z)] + } + + /// Determines whether the point (x, y) is inside the mesh. + fn inside(&self, x: f64, y: f64, z: f64) -> bool { + let p = Point2d::new(x, y); + if self.boundary.aabb(&Isometry::identity()).contains_local_point(&p) { + let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); + point_projection.is_inside + } else { + false + } + } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let p = Point2d::new(x, y); + + let (point_projection, (_, _)) = self.boundary.project_local_point_assuming_solid_interior_ccw(p); + let (x_intersect, y_intersect, z_intersect) = self.closest_point(x, y, z); + + let dx = x_intersect - x; + let dy = y_intersect - y; + let dz = z_intersect - z; + let mag = (dx*dx + dy*dy + dz*dz).sqrt(); + + if point_projection.is_inside { + (dx/mag, dy/mag, dz/mag) + } else { + (-dx/mag, -dy/mag, -dz/mag) + } + } + + fn new(geometry_input: &<::InputFileFormat as GeometryInput>::GeometryInput) -> ParryMesh2D { + let length_unit: f64 = match geometry_input.length_unit.as_str() { + "MICRON" => MICRON, + "CM" => CM, + "MM" => MM, + "ANGSTROM" => ANGSTROM, + "NM" => NM, + "M" => 1., + _ => geometry_input.length_unit.parse() + .expect(format!( + "Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M", + &geometry_input.length_unit.as_str() + ).as_str()), + }; + + let mut boundary_points_converted: Vec> = geometry_input.points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect(); + let number_boundary_points = boundary_points_converted.len() as u32; + + dbg!(&geometry_input.simulation_boundary_points); + + let mut simulation_boundary_points_converted: Vec> = geometry_input.simulation_boundary_points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect(); + let number_simulation_boundary_points = simulation_boundary_points_converted.len() as u32; + + dbg!(number_simulation_boundary_points); + + let test_simulation_boundary_ccw = (0..number_simulation_boundary_points as usize) + .map(|i| (simulation_boundary_points_converted[(i + 1) % number_simulation_boundary_points as usize].x - simulation_boundary_points_converted[i].x)*(simulation_boundary_points_converted[i].y + simulation_boundary_points_converted[(i + 1) % number_simulation_boundary_points as usize].y)) + .sum::(); + + dbg!(test_simulation_boundary_ccw > 0.0); + + dbg!(&simulation_boundary_points_converted); + + if test_simulation_boundary_ccw > 0.0 { + simulation_boundary_points_converted = simulation_boundary_points_converted.clone().into_iter().rev().collect::>>(); + } + + dbg!(&simulation_boundary_points_converted); + + for p in simulation_boundary_points_converted.iter().combinations(2) { + assert!(p[0] != p[1], "Input error: duplicate vertices {}, {} in simulation boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).", p[0], p[1]) + } + + let test_ccw = (0..number_boundary_points as usize) + .map(|i| (boundary_points_converted[(i + 1) % number_boundary_points as usize].x - boundary_points_converted[i].x)*(boundary_points_converted[i].y + boundary_points_converted[(i + 1) % number_boundary_points as usize].y)) + .sum::(); + + if test_ccw > 0.0 { + boundary_points_converted = boundary_points_converted.clone().into_iter().rev().collect::>>(); + } + + for p in boundary_points_converted.iter().combinations(2) { + assert!(p[0] != p[1], "Input error: duplicate vertices in boundary geometry input. Boundary must be defined ccw with no duplicate points (terminal segment will span from final to initial point).") + } + + let triangles = geometry_input.triangles.iter().map(|(i, j, k)| [*i as u32, *j as u32, *k as u32]).collect::>(); + let points_converted = geometry_input.points.iter().map(|(x, y)| Point2d::new(x*length_unit, y*length_unit)).collect::>>(); + + let trimesh = TriMesh::new(points_converted, triangles); + + let electronic_stopping_correction_factors = geometry_input.electronic_stopping_correction_factors.clone(); + + let densities: Vec> = geometry_input.densities.iter().map(|density_list| density_list.iter().map(|density| density / (length_unit).powi(3)).collect::>()).collect::>>(); + + + let energy_barrier_thickness = geometry_input.energy_barrier_thickness*length_unit; + + let concentrations: Vec> = densities.iter().map(|density_list| density_list.iter().map(|density| density / density_list.iter().sum::() ).collect::>()).collect::>>(); + + for concentration in &concentrations { + assert!((concentration.iter().sum::() - 1.0).abs() < 1e-6); + } + + let mut linked_boundary_points = (0..number_boundary_points).zip(1..number_boundary_points).map(|(x, y)| [x, y]).collect::>(); + linked_boundary_points.push([number_boundary_points - 1, 0]); + let boundary2 = Polyline::new(boundary_points_converted, Some(linked_boundary_points)); + + let number_simulation_boundary_points = simulation_boundary_points_converted.clone().len() as u32; + let mut linked_simulation_boundary_points = (0..number_simulation_boundary_points).zip(1..number_simulation_boundary_points).map(|(x, y)| [x, y]).collect::>(); + linked_simulation_boundary_points.push([number_simulation_boundary_points - 1, 0]); + let simulation_boundary2 = Polyline::new(simulation_boundary_points_converted, Some(linked_simulation_boundary_points)); + + ParryMesh2D { + trimesh: trimesh.expect("Input Error: Trimesh failed to build. Check mesh vertices and points."), + densities, + simulation_boundary: simulation_boundary2, + boundary: boundary2, + electronic_stopping_correction_factors, + energy_barrier_thickness, + concentrations + } + } +} + +/* /// Triangular mesh for rustbca. #[derive(Clone)] pub struct Mesh2D { @@ -491,26 +881,6 @@ impl Geometry for Mesh2D { cells.push(Cell2D::new(coordinate_set_converted, densities, concentrations, ck)); } - /* - for ((coordinate_set, densities), ck) in triangles.iter().zip(densities).zip(electronic_stopping_correction_factors) { - let coordinate_set_converted = ( - coordinate_set.0*length_unit, - coordinate_set.1*length_unit, - coordinate_set.2*length_unit, - coordinate_set.3*length_unit, - coordinate_set.4*length_unit, - coordinate_set.5*length_unit, - ); - - let total_density: f64 = densities.iter().sum(); - let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); - - cells.push(Cell2D::new(coordinate_set_converted, densities, concentrations, ck)); - } - */ - - - let mut boundary_points_converted = Vec::with_capacity(material_boundary_point_indices.len()); for index in material_boundary_point_indices.iter() { boundary_points_converted.push((points[*index].0*length_unit, points[*index].1*length_unit)); @@ -549,8 +919,18 @@ impl Geometry for Mesh2D { fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) { + if p.x() == x && p.y() == y { + dbg!(&self.boundary); + println!("single point"); + panic!("({} {} {}) ({} {} {})", p.x(), p.y(), z, x, y, z); + } (p.x(), p.y(), z) } else if let Closest::Intersection(p) = self.boundary.closest_point(&point!(x: x, y: y)) { + if p.x() == x && p.y() == y { + dbg!(&self.boundary); + println!("intersection"); + panic!("({} {} {}) ({} {} {})", p.x(), p.y(), z, x, y, z); + } (p.x(), p.y(), z) } else { panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z); @@ -617,6 +997,10 @@ impl Geometry for Mesh2D { fn inside(&self, x: f64, y: f64, z: f64) -> bool { self.boundary.contains(&Point::new(x, y)) } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + panic!("Not implemented."); + } } /// A mesh cell that contains a triangle and the local number densities and concentrations. @@ -662,6 +1046,7 @@ impl GeometryElement for Cell2D { self.electronic_stopping_correction_factor } } +*/ #[derive(Clone)] pub struct Layer1D { @@ -711,6 +1096,7 @@ impl GeometryElement for Layer1D { } } +/* /// A triangle in 2D, with points (x1, y1), (x2, y2), (x3, y3), and the three line segments bewtween them. #[derive(Clone)] pub struct Triangle2D { @@ -801,3 +1187,4 @@ impl Triangle2D { ((x - centroid.0)*(x - centroid.0) + (y - centroid.1)*(y - centroid.1)).sqrt() } } +*/ \ No newline at end of file diff --git a/src/lib.rs b/src/lib.rs index 41ea254..1808d12 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -56,6 +56,7 @@ use pyo3::types::*; //Load internal modules pub mod material; pub mod particle; +#[cfg(test)] pub mod tests; pub mod interactions; pub mod bca; @@ -75,7 +76,7 @@ pub use crate::consts::*; pub use crate::structs::*; pub use crate::input::{Input2D, InputHomogeneous2D, Input1D, Input0D, Options, InputFile, GeometryInput}; pub use crate::output::{OutputUnits}; -pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D}; +pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, ParryMesh2D}; pub use crate::sphere::{Sphere, SphereInput, InputSphere}; #[cfg(feature = "parry3d")] diff --git a/src/main.rs b/src/main.rs index 7da24c6..a66f4d7 100644 --- a/src/main.rs +++ b/src/main.rs @@ -58,7 +58,7 @@ pub use crate::consts::*; pub use crate::structs::*; pub use crate::input::{Input2D, InputHomogeneous2D, Input1D, Input0D, Options, InputFile, GeometryInput}; pub use crate::output::{OutputUnits}; -pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D, HomogeneousMesh2D}; +pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, ParryMesh2D, ParryHomogeneousMesh2D}; pub use crate::sphere::{Sphere, SphereInput, InputSphere}; pub use crate::physics::{physics_loop}; @@ -84,7 +84,7 @@ fn main() { "HOMOGENEOUS2D" => GeometryType::HOMOGENEOUS2D, _ => panic!("Unimplemented geometry {}.", args[1].clone()) }), - _ => panic!("Too many command line arguments. RustBCA accepts 0 (use 'input.toml') 1 () or 2 ( )"), + _ => panic!("Too many command line arguments. RustBCA accepts 0 (defaulting to 'input.toml' and 2D mesh), 1 ( defaulting to 2D mesh), or 2 ( )"), }; match geometry_type { @@ -97,8 +97,8 @@ fn main() { physics_loop::(particle_input_array, material, options, output_units); }, GeometryType::MESH2D => { - let (particle_input_array, material, options, output_units) = input::input::(input_file); - physics_loop::(particle_input_array, material, options, output_units); + let (particle_input_array, material, options, output_units) = input::input::(input_file); + physics_loop::(particle_input_array, material, options, output_units); }, GeometryType::SPHERE => { let (particle_input_array, material, options, output_units) = input::input::(input_file); @@ -115,8 +115,8 @@ fn main() { physics_loop::(particle_input_array, material, options, output_units); } GeometryType::HOMOGENEOUS2D => { - let (particle_input_array, material, options, output_units) = input::input::(input_file); - physics_loop::(particle_input_array, material, options, output_units); + let (particle_input_array, material, options, output_units) = input::input::(input_file); + physics_loop::(particle_input_array, material, options, output_units); } } } diff --git a/src/material.rs b/src/material.rs index e6931fa..1367f56 100644 --- a/src/material.rs +++ b/src/material.rs @@ -300,7 +300,7 @@ impl Material { let stopping_power = match electronic_stopping_mode { //Biersack-Varelas Interpolation - ElectronicStoppingMode::INTERPOLATED => 1./(1./S_high + 1./S_low*ck), + ElectronicStoppingMode::INTERPOLATED => ck/(1./S_high + 1./S_low), //Oen-Robinson ElectronicStoppingMode::LOW_ENERGY_LOCAL => S_low*ck, //Lindhard-Scharff @@ -339,14 +339,28 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, let leaving = !inside_now & inside_old; let entering = inside_now & !inside_old; + assert!(!x.is_nan(), "Particle x is NaN before surface binding energy calculation."); + assert!(!y.is_nan(), "Particle y is NaN before surface binding energy calculation."); + assert!(!z.is_nan(), "Particle z is NaN before surface binding energy calculation."); + + assert!(!cosx.is_nan(), "Particle dirx is NaN before surface binding energy calculation."); + assert!(!cosy.is_nan(), "Particle diry is NaN before surface binding energy calculation."); + assert!(!cosz.is_nan(), "Particle dirz is NaN before surface binding energy calculation."); + if entering { if particle_1.backreflected { particle_1.backreflected = false; } else { - let (x2, y2, z2) = material.closest_point(x, y, z); - let dx = x2 - x; - let dy = y2 - y; - let dz = z2 - z; + let (x2, y2, z2) = material.closest_point(x_old, y_old, z_old); + assert!(!x2.is_nan(), "Closest x is NaN in normal calculation."); + assert!(!y2.is_nan(), "Closest y is NaN in normal calculation."); + assert!(!z2.is_nan(), "Closest z is NaN in normal calculation."); + let dx = x2 - x_old; + let dy = y2 - y_old; + let dz = z2 - z_old; + assert!(!dx.is_nan(), "dx is NaN in normal calculation."); + assert!(!dy.is_nan(), "dy is NaN in normal calculation."); + assert!(!dz.is_nan(), "dz is NaN in normal calculation."); let mag = (dx*dx + dy*dy + dz*dz).sqrt(); let costheta = dx*cosx/mag + dy*cosy/mag + dz*cosz/mag; @@ -360,10 +374,17 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, if leaving { let (x2, y2, z2) = material.closest_point(x, y, z); + assert!(!x2.is_nan(), "Closest x is NaN in normal calculation."); + assert!(!y2.is_nan(), "Closest y is NaN in normal calculation."); + assert!(!z2.is_nan(), "Closest z is NaN in normal calculation."); let dx = x2 - x; let dy = y2 - y; let dz = z2 - z; + assert!(!dx.is_nan(), "dx is NaN in normal calculation."); + assert!(!dy.is_nan(), "dy is NaN in normal calculation."); + assert!(!dz.is_nan(), "dz is NaN in normal calculation."); let mag = (dx*dx + dy*dy + dz*dz).sqrt(); + let costheta = dx*cosx/mag + dy*cosy/mag + dz*cosz/mag; let leaving_energy = match material.surface_binding_model { @@ -414,4 +435,6 @@ pub fn boundary_condition_planar(particle_1: &mut particle::Particl if !particle_1.stopped & !particle_1.left { surface_binding_energy(particle_1, material); } + + assert!(!particle_1.pos.x.is_nan(), "Particle position is NaN: ({}, {}, {}) old: ({}, {}, {}), left? {} stopped? {}", x, y, z, particle_1.pos_old.x, particle_1.pos_old.y, particle_1.pos_old.z, particle_1.left, particle_1.stopped); } diff --git a/src/parry.rs b/src/parry.rs index 8686683..e7ca7bb 100644 --- a/src/parry.rs +++ b/src/parry.rs @@ -2,7 +2,7 @@ use super::*; use parry3d_f64::shape::{Ball, TriMesh}; use parry3d_f64::query::{PointQuery, Ray, RayCast}; use parry3d_f64::math::{Isometry, Point, Vector}; -use parry3d_f64::bounding_volume::AABB; +use parry3d_f64::bounding_volume::Aabb; #[derive(Deserialize, Clone)] pub struct InputParryBall { @@ -124,6 +124,15 @@ impl Geometry for ParryBall { let (x_, y_, z_) = (point_projection.point.x, point_projection.point.y, point_projection.point.z); (x_ as f64, y_ as f64, z_ as f64) } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let r = (x.powi(2) + y.powi(2) + z.powi(2)).sqrt(); + let R = self.radius; + let ux = x/r; + let uy = y/r; + let uz = z/r; + (ux, uy, uz) + } } @@ -171,7 +180,7 @@ pub struct ParryTriMesh { pub electronic_stopping_correction_factor: f64, pub energy_barrier_thickness: f64, pub trimesh: TriMesh, - pub boundary: AABB, + pub boundary: Aabb, } impl GeometryInput for InputParryTriMesh { @@ -204,7 +213,7 @@ impl Geometry for ParryTriMesh { let energy_barrier_thickness = total_density.powf(-1./3.)/SQRTPI*2.; let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); let points = input.vertices.iter().map(|p| Point::new(p[0]*length_unit , p[1]*length_unit , p[2]*length_unit)).collect(); - let trimesh = TriMesh::new(points, input.indices.clone()); + let trimesh = TriMesh::new(points, input.indices.clone()).expect("Input error: failed to build Trimesh. Check vertices, indices"); let boundary = trimesh.aabb(&Isometry::identity()); ParryTriMesh { @@ -263,6 +272,22 @@ impl Geometry for ParryTriMesh { let (x_, y_, z_) = (point_projection.point.x, point_projection.point.y, point_projection.point.z); (x_ as f64, y_ as f64, z_ as f64) } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let p = Point::new(x, y, z); + let point_projection = self.trimesh.project_local_point(&p, false); + + let dx = point_projection.point.x - x; + let dy = point_projection.point.y - y; + let dz = point_projection.point.z - z; + let mag = (dx*dx + dy*dy + dz*dz).sqrt(); + + if point_projection.is_inside { + (dx/mag, dy/mag, dz/mag) + } else { + (-dx/mag, -dy/mag, -dz/mag) + } + } } fn inside_trimesh(trimesh: &TriMesh, x: f64, y: f64, z: f64) -> bool { diff --git a/src/particle.rs b/src/particle.rs index 79cc09d..f4eef54 100644 --- a/src/particle.rs +++ b/src/particle.rs @@ -65,7 +65,7 @@ pub struct ParticleInput { } /// Particle object. Particles in rustbca include incident ions and material atoms. -#[derive(Clone)] +#[derive(Clone, Debug)] pub struct Particle { pub m: f64, pub Z: f64, @@ -285,6 +285,9 @@ impl Particle { //In order to keep average denisty constant, must add back previous asymptotic deflection let distance_traveled = mfp + self.asymptotic_deflection - asymptotic_deflection; + assert!(!distance_traveled.is_nan(), "Travel distance is NaN, mfp: {}, asymp. defl: {}, old asymp. defl: {}", mfp, self.asymptotic_deflection, asymptotic_deflection); + assert!(!self.dir_old.x.is_nan(), "Particle direction is NaN, dir: ({}, {}, {})", self.dir_old.x, self.dir_old.y, self.dir_old.z); + //dir has been updated, so use previous direction to advance in space self.pos.x += self.dir_old.x*distance_traveled; self.pos.y += self.dir_old.y*distance_traveled; @@ -308,10 +311,20 @@ pub fn surface_refraction(particle: &mut Particle, normal: Vector, Es: f64) { let u1x = (E/(E + Es)).sqrt()*particle.dir.x + ((-(E).sqrt()*costheta + (E*costheta.powi(2) + Es).sqrt())/(E + Es).sqrt())*normal.x; let u1y = (E/(E + Es)).sqrt()*particle.dir.y + ((-(E).sqrt()*costheta + (E*costheta.powi(2) + Es).sqrt())/(E + Es).sqrt())*normal.y; let u1z = (E/(E + Es)).sqrt()*particle.dir.z + ((-(E).sqrt()*costheta + (E*costheta.powi(2) + Es).sqrt())/(E + Es).sqrt())*normal.z; + + particle.E += Es; + + if u1x.is_nan() { + println!("normal: ({}, {}, {})", normal.x, normal.y, normal.z); + dbg!(&particle); + + } + particle.dir.x = u1x; particle.dir.y = u1y; particle.dir.z = u1z; - particle.E += Es; + + assert!(!u1x.is_nan(), "Particle direction after surface refraction is NaN. dir: ({}, {}, {}), old dir: ({}, {}, {}), normal: ({}, {}, {}), costheta: {}", u1x, u1y, u1z, particle.dir_old.x, particle.dir_old.y, particle.dir_old.z, normal.x, normal.y, normal.z, costheta); } /// Calcualte the refraction angle based on the surface binding energy of the material. diff --git a/src/sphere.rs b/src/sphere.rs index 92553b6..6f4055a 100644 --- a/src/sphere.rs +++ b/src/sphere.rs @@ -120,4 +120,13 @@ impl Geometry for Sphere { let uz = z/r; (ux*R, uy*R, uz*R) } + + fn nearest_normal_vector(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + let r = (x.powi(2) + y.powi(2) + z.powi(2)).sqrt(); + let R = self.radius; + let ux = x/r; + let uy = y/r; + let uz = z/r; + (ux, uy, uz) + } } diff --git a/src/structs.rs b/src/structs.rs index 47f20a0..fa5993d 100644 --- a/src/structs.rs +++ b/src/structs.rs @@ -1,5 +1,5 @@ /// 3D vector. -#[derive(Clone, Copy)] +#[derive(Clone, Copy, Debug)] pub struct Vector { pub x: f64, pub y: f64, @@ -45,7 +45,7 @@ impl Vector { } /// TrajectoryElement is a trajectory-tracking object that includes x, y, z, and the current energy. -#[derive(Clone)] +#[derive(Clone, Debug)] pub struct TrajectoryElement { pub E: f64, pub x: f64, @@ -65,7 +65,7 @@ impl TrajectoryElement { } /// Energy loss is an output tracker that tracks the separate nuclear and electronic energy losses. -#[derive(Clone)] +#[derive(Clone, Debug)] pub struct EnergyLoss { pub En: f64, pub Ee: f64, diff --git a/src/tests.rs b/src/tests.rs index b926984..8e2d088 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -3,6 +3,101 @@ use super::*; #[cfg(test)] use float_cmp::*; +#[cfg(test)] +fn test_parry2d() { + let points: Vec> = vec![Point2d::new(-1.0, -1.0), Point2d::new(-1.0, 1.0), Point2d::new(1.0, 1.0), Point2d::new(1.0, -1.0)]; + + let indices: Vec<[u32; 2]> = vec![[0, 3], [3, 2], [2, 1], [1, 0]]; + //indices.reverse(); + + //(0.0000006261797114005236, 0.000006009591179670447) + let query_point = Point2d::new(0.0, 0.0); + + let polyline = Polyline::new(points, Some(indices)); + + let (point_projection, (_, _)) = polyline.project_local_point_assuming_solid_interior_ccw( + query_point + ); + + assert!(point_projection.is_inside); + assert!(polyline.aabb(&Isometry::identity()).contains_local_point(&query_point)); + + let test_geometry = vec![(0.25, 0.2), (0.75, 0.2), (0.8, 0.8), (0.5, 0.5), (0.2, 0.8)]; + let num_segments = 5; + + let geometry_input_homogeneous_2D = geometry::HomogeneousMesh2DInput { + length_unit: "M".to_string(), + points: test_geometry.clone(), + densities: vec![0.03, 0.03], + simulation_boundary_points: vec![(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)], + electronic_stopping_correction_factor: 1.0, + }; + + let material_parameters = material::MaterialParameters{ + energy_unit: "EV".to_string(), + mass_unit: "AMU".to_string(), + Eb: vec![0.0, 0.0], + Es: vec![2.0, 4.0], + Ec: vec![1.0, 1.0], + Ed: vec![0.0, 0.0], + Z: vec![29., 1.], + m: vec![63.54, 1.0008], + interaction_index: vec![0, 0], + surface_binding_model: SurfaceBindingModel::PLANAR{calculation: SurfaceBindingCalculation::TARGET}, + bulk_binding_model: BulkBindingModel::INDIVIDUAL, + }; + + let material_homogeneous_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_homogeneous_2D); + + let query_points_outside = vec![[0.5, 0.1, 0.0], [0.9, 0.5, 0.0], [0.6, 0.7, 0.0], [0.4, 0.7, 0.0], [0.1, 0.5, 0.0]]; + let query_points_inside = vec![[0.5, 0.21, 0.0], [0.75, 0.21, 0.0], [0.51, 0.5, 0.0], [0.49, 0.5, 0.0], [0.25, 0.21, 0.0]]; + + for query_point in query_points_outside.clone() { + assert!(!material_homogeneous_2D.inside(query_point[0], query_point[1], query_point[2]), "Point ({}, {}, {}) failed outsidedness check.", query_point[0], query_point[1], query_point[2]); + } + + for query_point in query_points_inside.clone() { + assert!(material_homogeneous_2D.inside(query_point[0], query_point[1], query_point[2]), "Point ({}, {}, {}) failed outsidedness check.", query_point[0], query_point[1], query_point[2]); + } + + + let normals_outside = query_points_outside + .iter() + .map( |query_point| { + material_homogeneous_2D.geometry.nearest_normal_vector(query_point[0], query_point[1], query_point[2]) + } + ).collect::>(); + + let normals_inside = query_points_inside + .iter() + .map( |query_point| { + material_homogeneous_2D.geometry.nearest_normal_vector(query_point[0], query_point[1], query_point[2]) + } + ).collect::>(); + + let test_normals = (0..num_segments) + .zip(1..num_segments + 1) + .map(|(i, j)| { + let dx = test_geometry[i % num_segments].0 - test_geometry[j % num_segments].0; + let dy = test_geometry[i % num_segments].1 - test_geometry[j % num_segments].1; + let mag = (dx*dx + dy*dy).sqrt(); + (-dy/mag, dx/mag, 0.0) + }).collect::>(); + + for (test_normal, normal) in test_normals.iter().zip(normals_outside) { + assert!(approx_eq!(f64, normal.0, test_normal.0, epsilon=1E-9), "Normal vector check failed. Calculated Normal: ({}, {}, {}) Test: ({}, {}, {})", normal.0, normal.1, normal.2, test_normal.0, test_normal.1, test_normal.2); + assert!(approx_eq!(f64, normal.1, test_normal.1, epsilon=1E-9)); + assert!(approx_eq!(f64, normal.2, test_normal.2, epsilon=1E-9)); + } + + for (test_normal, normal) in test_normals.iter().zip(normals_inside) { + assert!(approx_eq!(f64, normal.0, test_normal.0, epsilon=1E-9), "Normal vector check failed. Calculated Normal: ({}, {}, {}) Test: ({}, {}, {})", normal.0, normal.1, normal.2, test_normal.0, test_normal.1, test_normal.2); + assert!(approx_eq!(f64, normal.1, test_normal.1, epsilon=1E-9)); + assert!(approx_eq!(f64, normal.2, test_normal.2, epsilon=1E-9)); + } + +} + #[test] #[cfg(feature = "cpr_rootfinder")] @@ -439,6 +534,61 @@ fn test_spherical_geometry() { } +#[test] +fn extended_test_2D_geometry() { + let mass = 1.; + let Z = 1.; + let E = 10.*EV; + let Ec = 1.*EV; + let Es = 5.76*EV; + let x = 1.5e-6; + let y = 10.0e-6; + let z = 0.; + let cosx = 0.0; + let cosy = -1.0; + let cosz = 0.; + + let mut particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, 0.0, x, y, z, cosx, cosy, cosz, false, false, 0); + + let material_parameters = material::MaterialParameters{ + energy_unit: "EV".to_string(), + mass_unit: "AMU".to_string(), + Eb: vec![0.0, 0.0], + Es: vec![2.0, 4.0], + Ec: vec![1.0, 1.0], + Ed: vec![0.0, 0.0], + Z: vec![29., 1.], + m: vec![63.54, 1.0008], + interaction_index: vec![0, 0], + surface_binding_model: SurfaceBindingModel::PLANAR{calculation: SurfaceBindingCalculation::TARGET}, + bulk_binding_model: BulkBindingModel::INDIVIDUAL, + }; + + let thickness: f64 = 1000.; + let depth: f64 = 1000.; + + let geometry_input_2D = geometry::Mesh2DInput { + length_unit: "MICRON".to_string(), + points: vec![(0.0, 0.0), (1.0, 0.0), (2.0, 0.0), (3.0, 0.0), (4.0, 0.0), (4.0, 1.0), (3.0, 1.0), (3.0, 10.0), (2.0, 10.0), (2.0, 1.0), (1.0, 1.0), (1.0, 10.0), (0.0, 10.0), (0.0, 1.0)], + triangles: vec![(13, 12, 11), (13, 11, 10), (0, 13, 10), (0, 10, 1), (1, 10, 9), (1, 9, 2), (9, 8, 7), (9, 7, 6), (2, 9, 6), (2, 6, 3), (3, 6, 5), (3, 5, 4)], + densities: vec![vec![3e10, 3e10]; 12], + boundary: vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13], + simulation_boundary_points: vec![(-0.1, -0.1), (4.1, -0.1), (4.1, 10.1), (-0.1, 10.1)], + energy_barrier_thickness: 4e-4, + electronic_stopping_correction_factors: vec![1.0; 12], + }; + + let material_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_2D); + + + assert!(material_2D.inside(1.5e-6, 0.999e-6, 0.0)); + assert!(!material_2D.inside(1.5e-6, 1.0001e-6, 0.0)); + + material::surface_binding_energy(&mut particle_1, &material_2D); + + material::boundary_condition_planar(&mut particle_1, &material_2D); +} + #[test] fn test_geometry() { let mass = 1.; @@ -477,7 +627,7 @@ fn test_geometry() { triangles: vec![(0, 1, 2), (0, 2, 3)], densities: vec![vec![0.03, 0.03], vec![0.03, 0.03]], boundary: vec![0, 1, 2, 3], - simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.), (0., 1.1*thickness/2.)], + simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.)], energy_barrier_thickness: 10., electronic_stopping_correction_factors: vec![1.0, 1.0], }; @@ -486,7 +636,7 @@ fn test_geometry() { length_unit: "ANGSTROM".to_string(), points: vec![(0., -thickness/2.), (depth, -thickness/2.), (depth, thickness/2.), (0., thickness/2.)], densities: vec![0.03, 0.03], - simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.), (0., 1.1*thickness/2.)], + simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.)], electronic_stopping_correction_factor: 1.0, }; @@ -504,8 +654,8 @@ fn test_geometry() { electronic_stopping_correction_factor: 1.0 }; - let material_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_2D); - let material_homogeneous_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_homogeneous_2D); + let material_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_2D); + let material_homogeneous_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_homogeneous_2D); let material_1D: material::Material = material::Material::::new(&material_parameters, &geometry_input_1D); let mut material_0D: material::Material = material::Material::::new(&material_parameters, &geometry_input_0D); material_0D.geometry.energy_barrier_thickness = 10.*ANGSTROM; @@ -529,7 +679,12 @@ fn test_geometry() { assert_eq!(material_2D.geometry.get_densities(0., 0., 0.), material_0D.geometry.get_densities(0., 0., 0.)); assert_eq!(material_2D.geometry.get_total_density(0., 0., 0.), material_0D.geometry.get_total_density(0., 0., 0.)); assert_eq!(material_2D.geometry.get_concentrations(0., 0., 0.), material_0D.geometry.get_concentrations(0., 0., 0.)); - assert_eq!(material_2D.geometry.closest_point(-10., 0., 5.), material_0D.geometry.closest_point(-10., 0., 5.)); + let (x2d, y2d, z2d) = material_2D.geometry.closest_point(0., 0., 0.); + let (x0d, y0d, z0d) = material_0D.geometry.closest_point(0., 0., 0.); + assert!(approx_eq!(f64, x2d, x0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d); + assert!(approx_eq!(f64, y2d, y0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d); + assert!(approx_eq!(f64, z2d, z0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d); + assert_eq!(material_2D.geometry.get_densities(-10., 0., 5.), material_0D.geometry.get_densities(-10., 0., 5.)); assert_eq!(material_2D.geometry.get_ck(-10., 0., 5.), material_0D.geometry.get_ck(-10., 0., 5.)); @@ -537,7 +692,14 @@ fn test_geometry() { assert_eq!(material_2D.geometry.get_densities(0., 0., 0.), material_homogeneous_2D.geometry.get_densities(0., 0., 0.)); assert_eq!(material_2D.geometry.get_total_density(0., 0., 0.), material_homogeneous_2D.geometry.get_total_density(0., 0., 0.)); assert_eq!(material_2D.geometry.get_concentrations(0., 0., 0.), material_homogeneous_2D.geometry.get_concentrations(0., 0., 0.)); - assert_eq!(material_2D.geometry.closest_point(-10., 0., 5.), material_homogeneous_2D.geometry.closest_point(-10., 0., 5.)); + + let (x2d, y2d, z2d) = material_2D.geometry.closest_point(-10., 0., 5.); + let (x0d, y0d, z0d) = material_0D.geometry.closest_point(-10., 0., 5.); + assert!(approx_eq!(f64, x2d, x0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d); + assert!(approx_eq!(f64, y2d, y0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d); + assert!(approx_eq!(f64, z2d, z0d, epsilon=1e-6), "2D: ({}, {}, {}) 0D: ({}, {}, {})", x2d, y2d, z2d, x0d, y0d, z0d); + + assert_eq!(material_2D.geometry.get_densities(-10., 0., 5.), material_homogeneous_2D.geometry.get_densities(-10., 0., 5.)); assert_eq!(material_2D.geometry.get_ck(-10., 0., 5.), material_homogeneous_2D.geometry.get_ck(-10., 0., 5.)); } @@ -580,7 +742,7 @@ fn test_surface_binding_energy_barrier() { triangles: vec![(0, 1, 2), (0, 2, 3)], densities: vec![vec![0.03, 0.03], vec![0.03, 0.03]], boundary: vec![0, 1, 2, 3], - simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.), (0., 1.1*thickness/2.)], + simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.)], energy_barrier_thickness: 10., electronic_stopping_correction_factors: vec![1.0, 1.0], }; @@ -591,7 +753,7 @@ fn test_surface_binding_energy_barrier() { electronic_stopping_correction_factor: 1.0 }; - let material_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_2D); + let material_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_2D); let mut material_0D: material::Material = material::Material::::new(&material_parameters, &geometry_input_0D); material_0D.geometry.energy_barrier_thickness = 10.*ANGSTROM; @@ -662,6 +824,8 @@ fn test_surface_binding_energy_barrier() { assert!(material_0D.inside_energy_barrier(1015.*ANGSTROM, 0., 0.)); } +//Deprecated with geo +/* #[test] fn test_triangle_contains() { let triangle_1 = geometry::Triangle2D::new((0., 2., 0., 2., 0., 0.)); @@ -685,6 +849,7 @@ fn test_triangle_distance_to() { assert!(approx_eq!(f64, triangle_1.distance_to(2., 0.), 0., epsilon=1E-12), "{}", triangle_1.distance_to(2., 0.)); assert!(approx_eq!(f64, triangle_1.distance_to(0., 2.), 0., epsilon=1E-12), "{}", triangle_1.distance_to(0., 2.)); } +*/ #[test] fn test_surface_refraction() { @@ -778,7 +943,7 @@ fn test_momentum_conservation() { InteractionPotential::LENZ_JENSEN, InteractionPotential::MORSE{D: 5.4971E-20, r0: 2.782E-10, alpha: 1.4198E10} ]; - + let mut rootfinders = vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-3}; 4]; //[[{"CPR"={n0=2, nmax=100, epsilon=1E-9, complex_threshold=1E-3, truncation_threshold=1E-9, far_from_zero=1E9, interval_limit=1E-12, derivative_free=true}}]] @@ -796,6 +961,13 @@ fn test_momentum_conservation() { } ); + //Aluminum + let m2 = 6.941; + let Z2 = 3.; + let Ec2 = 1.; + let Es2 = 1.; + + for direction in vec![(1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0), (-1.0, 0.0, 0.0), (0.0, -1.0, 0.0), (0.0, 0.0, -1.0)] { for energy_eV in vec![1., 10., 100., 1000., 1000.] { //Aluminum @@ -818,6 +990,19 @@ fn test_momentum_conservation() { let cosy = direction.1; let cosz = direction.2; + let thickness: f64 = 1000.; + let depth: f64 = 1000.; + let geometry_input = geometry::Mesh2DInput { + length_unit: "ANGSTROM".to_string(), + triangles: vec![(0, 1, 2), (0, 2, 3)], + points: vec![(0., -thickness/2.), (depth, -thickness/2.), (depth, thickness/2.), (0., thickness/2.)], + densities: vec![vec![0.06306], vec![0.06306]], + boundary: vec![0, 1, 2, 3], + simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.)], + electronic_stopping_correction_factors: vec![0.0, 0.0], + energy_barrier_thickness: 0., + }; + let material_parameters = material::MaterialParameters{ energy_unit: "EV".to_string(), mass_unit: "AMU".to_string(), @@ -832,20 +1017,7 @@ fn test_momentum_conservation() { bulk_binding_model: BulkBindingModel::INDIVIDUAL, }; - let thickness: f64 = 1000.; - let depth: f64 = 1000.; - let geometry_input = geometry::Mesh2DInput { - length_unit: "ANGSTROM".to_string(), - triangles: vec![(0, 1, 2), (0, 2, 3)], - points: vec![(0., -thickness/2.), (depth, -thickness/2.), (depth, thickness/2.), (0., thickness/2.)], - densities: vec![vec![0.06306], vec![0.06306]], - boundary: vec![0, 1, 2, 3], - simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.), (0., 1.1*thickness/2.)], - electronic_stopping_correction_factors: vec![0.0, 0.0], - energy_barrier_thickness: 0., - }; - - let material_1: material::Material = material::Material::::new(&material_parameters, &geometry_input); + let material_1: material::Material = material::Material::::new(&material_parameters, &geometry_input); for high_energy_free_flight_paths in vec![true, false] { for (potential, root_finder) in potentials.iter().zip(rootfinders.clone()) { @@ -857,7 +1029,7 @@ fn test_momentum_conservation() { _ => {} } - println!("Case: {} {} {} {}", energy_eV, high_energy_free_flight_paths, potential, scattering_integral); + //println!("Case: {} {} {} {}", energy_eV, high_energy_free_flight_paths, potential, scattering_integral); let mut particle_1 = particle::Particle::new(m1, Z1, E1, Ec1, Es1, 0.0, x1, y1, z1, cosx, cosy, cosz, false, false, 0); @@ -922,9 +1094,11 @@ fn test_momentum_conservation() { let binary_collision_geometries = bca::determine_mfp_phi_impact_parameter(&mut particle_1, &material_1, &options); + /* println!("Phi: {} rad p: {} Angstrom mfp: {} Angstrom", binary_collision_geometries[0].phi_azimuthal, binary_collision_geometries[0].impact_parameter/ANGSTROM, binary_collision_geometries[0].mfp/ANGSTROM); + */ let (species_index, mut particle_2) = bca::choose_collision_partner(&mut particle_1, &material_1, &binary_collision_geometries[0], &options); @@ -937,12 +1111,13 @@ fn test_momentum_conservation() { let binary_collision_result = bca::calculate_binary_collision(&particle_1, &particle_2, &binary_collision_geometries[0], &options).unwrap(); + /* println!("E_recoil: {} eV Psi: {} rad Psi_recoil: {} rad", binary_collision_result.recoil_energy/EV, binary_collision_result.psi, binary_collision_result.psi_recoil); println!("Initial Energies: {} eV {} eV", particle_1.E/EV, particle_2.E/EV); - + */ //Energy transfer to recoil particle_2.E = binary_collision_result.recoil_energy - material_1.average_bulk_binding_energy(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); @@ -963,11 +1138,13 @@ fn test_momentum_conservation() { let final_momentum = mom1_1.add(&mom2_1); + /* println!("Final Energies: {} eV {} eV", particle_1.E/EV, particle_2.E/EV); println!("X: {} {} {}% Error", initial_momentum.x/ANGSTROM/AMU, final_momentum.x/ANGSTROM/AMU, 100.*(final_momentum.x - initial_momentum.x)/initial_momentum.magnitude()); println!("Y: {} {} {}% Error", initial_momentum.y/ANGSTROM/AMU, final_momentum.y/ANGSTROM/AMU, 100.*(final_momentum.y - initial_momentum.y)/initial_momentum.magnitude()); println!("Z: {} {} {}% Error", initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, 100.*(final_momentum.z - initial_momentum.z)/initial_momentum.magnitude()); println!(); + */ //These values are in [angstrom amu / second], so very large. assert!(approx_eq!(f64, initial_momentum.x/ANGSTROM/AMU, final_momentum.x/ANGSTROM/AMU, epsilon = 1000.));