diff --git a/.github/workflows/rustbca_compile_check.yml b/.github/workflows/rustbca_compile_check.yml index 17ed2de..427930b 100644 --- a/.github/workflows/rustbca_compile_check.yml +++ b/.github/workflows/rustbca_compile_check.yml @@ -26,7 +26,7 @@ jobs: sudo apt-get install python3-pip - name: Install Python libraries run: | - pip3 + pip3 python3 -m pip install numpy shapely scipy - name: Install python TOML library from source run: | @@ -42,4 +42,6 @@ jobs: - name: Run Examples run: | cargo run --release examples/boron_nitride.toml + cargo run --release 0D examples/boron_nitride_0D.toml + cargo run --release 0D examples/titanium_dioxide_0D.toml cargo run --release --features distributions examples/layered_geometry.toml diff --git a/Cargo.toml b/Cargo.toml index c004727..53ee5ff 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -9,12 +9,13 @@ name="rustBCA" [dependencies] rand = "0.8.3" -geo = "0.17.1" + toml = "0.5.8" anyhow = "1.0.38" itertools = "0.10.0" rayon = "1.5.0" -indicatif = {version="0.15.0", features=["rayon"]} +geo = {version = "0.17.1", optional = false} +indicatif = {version = "0.15.0", features=["rayon"]} serde = { version = "1.0.123", features = ["derive"] } hdf5 = {version = "0.7.1", optional = true} openblas-src = {version = "0.9", optional = true} @@ -39,3 +40,5 @@ cpr_rootfinder_netlib = ["rcpr", "netlib-src"] cpr_rootfinder_intel_mkl = ["rcpr", "intel-mkl-src"] distributions = ["ndarray"] no_list_output = [] +mesh_0D = [] +mesh_2D = [] diff --git a/examples/boron_nitride.toml b/examples/boron_nitride.toml index 8359226..047951e 100644 --- a/examples/boron_nitride.toml +++ b/examples/boron_nitride.toml @@ -44,7 +44,7 @@ dir = [ [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,], [ 0.9999999999984769 interaction_index = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0] particle_input_filename="" -[mesh_2d_input] +[geometry_input] length_unit = "MICRON" energy_barrier_thickness = 1.7645653881793786e-4 triangles = [ [ 0.0, 0.0, 0.0012558103905862675, 0.0, 0.02, 0.019960534568565433,], [ 0.0, 0.0012558103905862675, 0.0025066646712860853, 0.0, 0.019960534568565433, 0.019842294026289557,], [ 0.0, 0.0025066646712860853, 0.003747626291714493, 0.0, 0.019842294026289557, 0.019645745014573775,], [ 0.0, 0.003747626291714493, 0.004973797743297096, 0.0, 0.019645745014573775, 0.01937166322257262,], [ 0.0, 0.004973797743297096, 0.0061803398874989484, 0.0, 0.01937166322257262, 0.019021130325903073,], [ 0.0, 0.0061803398874989484, 0.00736249105369356, 0.0, 0.019021130325903073, 0.018595529717765027,], [ 0.0, 0.00736249105369356, 0.008515585831301454, 0.0, 0.018595529717765027, 0.01809654104932039,], [ 0.0, 0.008515585831301454, 0.009635073482034308, 0.0, 0.01809654104932039, 0.01752613360087727,], [ 0.0, 0.009635073482034308, 0.010716535899579934, 0.0, 0.01752613360087727, 0.0168865585100403,], [ 0.0, 0.010716535899579934, 0.011755705045849463, 0.0, 0.0168865585100403, 0.016180339887498948,], [ 0.0, 0.011755705045849463, 0.012748479794973795, 0.0, 0.016180339887498948, 0.015410264855515783,], [ 0.0, 0.012748479794973795, 0.013690942118573775, 0.0, 0.015410264855515783, 0.014579372548428232,], [ 0.0, 0.013690942118573775, 0.014579372548428232, 0.0, 0.014579372548428232, 0.013690942118573773,], [ 0.0, 0.014579372548428232, 0.015410264855515785, 0.0, 0.013690942118573773, 0.012748479794973793,], [ 0.0, 0.015410264855515785, 0.016180339887498948, 0.0, 0.012748479794973793, 0.011755705045849461,], [ 0.0, 0.016180339887498948, 0.0168865585100403, 0.0, 0.011755705045849461, 0.01071653589957993,], [ 0.0, 0.0168865585100403, 0.017526133600877274, 0.0, 0.01071653589957993, 0.009635073482034304,], [ 0.0, 0.017526133600877274, 0.018096541049320392, 0.0, 0.009635073482034304, 0.008515585831301454,], [ 0.0, 0.018096541049320392, 0.01859552971776503, 0.0, 0.008515585831301454, 0.007362491053693557,], [ 0.0, 0.01859552971776503, 0.019021130325903073, 0.0, 0.007362491053693557, 0.006180339887498949,], [ 0.0, 0.019021130325903073, 0.01937166322257262, 0.0, 0.006180339887498949, 0.004973797743297095,], [ 0.0, 0.01937166322257262, 0.019645745014573775, 0.0, 0.004973797743297095, 0.0037476262917144902,], [ 0.0, 0.019645745014573775, 0.019842294026289557, 0.0, 0.0037476262917144902, 0.0025066646712860853,], [ 0.0, 0.019842294026289557, 0.019960534568565433, 0.0, 0.0025066646712860853, 0.001255810390586266,], [ 0.0, 0.019960534568565433, 0.02, 0.0, 0.001255810390586266, -3.2162452993532734e-18,], [ 0.0, 0.02, 0.019960534568565433, 0.0, -3.2162452993532734e-18, -0.001255810390586268,], [ 0.0, 0.019960534568565433, 0.019842294026289557, 0.0, -0.001255810390586268, -0.0025066646712860875,], [ 0.0, 0.019842294026289557, 0.019645745014573772, 0.0, -0.0025066646712860875, -0.0037476262917144967,], [ 0.0, 0.019645745014573772, 0.01937166322257262, 0.0, -0.0037476262917144967, -0.004973797743297097,], [ 0.0, 0.01937166322257262, 0.019021130325903073, 0.0, -0.004973797743297097, -0.006180339887498951,], [ 0.0, 0.019021130325903073, 0.01859552971776503, 0.0, -0.006180339887498951, -0.00736249105369356,], [ 0.0, 0.01859552971776503, 0.01809654104932039, 0.0, -0.00736249105369356, -0.008515585831301454,], [ 0.0, 0.01809654104932039, 0.01752613360087727, 0.0, -0.008515585831301454, -0.00963507348203431,], [ 0.0, 0.01752613360087727, 0.0168865585100403, 0.0, -0.00963507348203431, -0.010716535899579938,], [ 0.0, 0.0168865585100403, 0.016180339887498948, 0.0, -0.010716535899579938, -0.011755705045849461,], [ 0.0, 0.016180339887498948, 0.015410264855515785, 0.0, -0.011755705045849461, -0.012748479794973795,], [ 0.0, 0.015410264855515785, 0.014579372548428228, 0.0, -0.012748479794973795, -0.013690942118573775,], [ 0.0, 0.014579372548428228, 0.01369094211857377, 0.0, -0.013690942118573775, -0.014579372548428234,], [ 0.0, 0.01369094211857377, 0.01274847979497379, 0.0, -0.014579372548428234, -0.015410264855515788,], [ 0.0, 0.01274847979497379, 0.011755705045849465, 0.0, -0.015410264855515788, -0.016180339887498948,], [ 0.0, 0.011755705045849465, 0.010716535899579934, 0.0, -0.016180339887498948, -0.0168865585100403,], [ 0.0, 0.010716535899579934, 0.009635073482034304, 0.0, -0.0168865585100403, -0.01752613360087727,], [ 0.0, 0.009635073482034304, 0.00851558583130145, 0.0, -0.01752613360087727, -0.018096541049320392,], [ 0.0, 0.00851558583130145, 0.007362491053693555, 0.0, -0.018096541049320392, -0.01859552971776503,], [ 0.0, 0.007362491053693555, 0.006180339887498942, 0.0, -0.01859552971776503, -0.019021130325903073,], [ 0.0, 0.006180339887498942, 0.004973797743297097, 0.0, -0.019021130325903073, -0.01937166322257262,], [ 0.0, 0.004973797743297097, 0.0037476262917144915, 0.0, -0.01937166322257262, -0.019645745014573775,], [ 0.0, 0.0037476262917144915, 0.002506664671286082, 0.0, -0.019645745014573775, -0.019842294026289557,], [ 0.0, 0.002506664671286082, 0.0012558103905862628, 0.0, -0.019842294026289557, -0.019960534568565433,], [ 0.0, 0.0012558103905862628, -6.432490598706547e-18, 0.0, -0.019960534568565433, -0.02,], [ 0.0, -6.432490598706547e-18, -0.0012558103905862669, 0.0, -0.02, -0.019960534568565433,], [ 0.0, -0.0012558103905862669, -0.0025066646712860858, 0.0, -0.019960534568565433, -0.019842294026289557,], [ 0.0, -0.0025066646712860858, -0.0037476262917144954, 0.0, -0.019842294026289557, -0.019645745014573772,], [ 0.0, -0.0037476262917144954, -0.0049737977432971, 0.0, -0.019645745014573772, -0.01937166322257262,], [ 0.0, -0.0049737977432971, -0.0061803398874989545, 0.0, -0.01937166322257262, -0.019021130325903073,], [ 0.0, -0.0061803398874989545, -0.007362491053693567, 0.0, -0.019021130325903073, -0.018595529717765024,], [ 0.0, -0.007362491053693567, -0.008515585831301454, 0.0, -0.018595529717765024, -0.01809654104932039,], [ 0.0, -0.008515585831301454, -0.009635073482034308, 0.0, -0.01809654104932039, -0.01752613360087727,], [ 0.0, -0.009635073482034308, -0.010716535899579936, 0.0, -0.01752613360087727, -0.0168865585100403,], [ 0.0, -0.010716535899579936, -0.011755705045849468, 0.0, -0.0168865585100403, -0.016180339887498944,], [ 0.0, -0.011755705045849468, -0.0127484797949738, 0.0, -0.016180339887498944, -0.015410264855515781,], [ 0.0, -0.0127484797949738, -0.013690942118573775, 0.0, -0.015410264855515781, -0.014579372548428232,], [ 0.0, -0.013690942118573775, -0.014579372548428232, 0.0, -0.014579372548428232, -0.013690942118573773,], [ 0.0, -0.014579372548428232, -0.015410264855515788, 0.0, -0.013690942118573773, -0.01274847979497379,], [ 0.0, -0.015410264855515788, -0.016180339887498948, 0.0, -0.01274847979497379, -0.011755705045849465,], [ 0.0, -0.016180339887498948, -0.016886558510040308, 0.0, -0.011755705045849465, -0.010716535899579927,], [ 0.0, -0.016886558510040308, -0.01752613360087727, 0.0, -0.010716535899579927, -0.009635073482034306,], [ 0.0, -0.01752613360087727, -0.018096541049320396, 0.0, -0.009635073482034306, -0.008515585831301443,], [ 0.0, -0.018096541049320396, -0.01859552971776503, 0.0, -0.008515585831301443, -0.007362491053693556,], [ 0.0, -0.01859552971776503, -0.019021130325903073, 0.0, -0.007362491053693556, -0.006180339887498951,], [ 0.0, -0.019021130325903073, -0.019371663222572624, 0.0, -0.006180339887498951, -0.004973797743297089,], [ 0.0, -0.019371663222572624, -0.019645745014573775, 0.0, -0.004973797743297089, -0.003747626291714493,], [ 0.0, -0.019645745014573775, -0.019842294026289557, 0.0, -0.003747626291714493, -0.0025066646712860745,], [ 0.0, -0.019842294026289557, -0.019960534568565433, 0.0, -0.0025066646712860745, -0.001255810390586264,], [ 0.0, -0.019960534568565433, -0.02, 0.0, -0.001255810390586264, -3.673940397442059e-18,], [ 0.0, -0.02, -0.019960534568565433, 0.0, -3.673940397442059e-18, 0.0012558103905862745,], [ 0.0, -0.019960534568565433, -0.019842294026289557, 0.0, 0.0012558103905862745, 0.0025066646712860845,], [ 0.0, -0.019842294026289557, -0.019645745014573772, 0.0, 0.0025066646712860845, 0.003747626291714503,], [ 0.0, -0.019645745014573772, -0.01937166322257262, 0.0, 0.003747626291714503, 0.004973797743297099,], [ 0.0, -0.01937166322257262, -0.019021130325903073, 0.0, 0.004973797743297099, 0.006180339887498945,], [ 0.0, -0.019021130325903073, -0.018595529717765024, 0.0, 0.006180339887498945, 0.007362491053693565,], [ 0.0, -0.018595529717765024, -0.018096541049320392, 0.0, 0.007362491053693565, 0.008515585831301452,], [ 0.0, -0.018096541049320392, -0.017526133600877267, 0.0, 0.008515585831301452, 0.009635073482034314,], [ 0.0, -0.017526133600877267, -0.0168865585100403, 0.0, 0.009635073482034314, 0.010716535899579936,], [ 0.0, -0.0168865585100403, -0.01618033988749894, 0.0, 0.010716535899579936, 0.011755705045849473,], [ 0.0, -0.01618033988749894, -0.015410264855515781, 0.0, 0.011755705045849473, 0.0127484797949738,], [ 0.0, -0.015410264855515781, -0.014579372548428232, 0.0, 0.0127484797949738, 0.013690942118573773,], [ 0.0, -0.014579372548428232, -0.013690942118573766, 0.0, 0.013690942118573773, 0.014579372548428239,], [ 0.0, -0.013690942118573766, -0.012748479794973793, 0.0, 0.014579372548428239, 0.015410264855515788,], [ 0.0, -0.012748479794973793, -0.011755705045849453, 0.0, 0.015410264855515788, 0.016180339887498955,], [ 0.0, -0.011755705045849453, -0.010716535899579927, 0.0, 0.016180339887498955, 0.016886558510040305,], [ 0.0, -0.010716535899579927, -0.009635073482034308, 0.0, 0.016886558510040305, 0.01752613360087727,], [ 0.0, -0.009635073482034308, -0.008515585831301445, 0.0, 0.01752613360087727, 0.018096541049320396,], [ 0.0, -0.008515585831301445, -0.007362491053693557, 0.0, 0.018096541049320396, 0.01859552971776503,], [ 0.0, -0.007362491053693557, -0.006180339887498935, 0.0, 0.01859552971776503, 0.019021130325903076,], [ 0.0, -0.006180339887498935, -0.00497379774329709, 0.0, 0.019021130325903076, 0.019371663222572624,], [ 0.0, -0.00497379774329709, -0.0037476262917144937, 0.0, 0.019371663222572624, 0.019645745014573775,], [ 0.0, -0.0037476262917144937, -0.002506664671286076, 0.0, 0.019645745014573775, 0.019842294026289557,], [ 0.0, -0.002506664671286076, -0.0012558103905862654, 0.0, 0.019842294026289557, 0.019960534568565433,], [ 0.0, -0.0012558103905862654, -4.898587196589413e-18, 0.0, 0.019960534568565433, 0.02,],] diff --git a/examples/boron_nitride_0D.toml b/examples/boron_nitride_0D.toml new file mode 100644 index 0000000..4aad3ca --- /dev/null +++ b/examples/boron_nitride_0D.toml @@ -0,0 +1,50 @@ +[options] +name = "boron_dust_grain_" +track_trajectories = false +track_recoils = true +track_recoil_trajectories = false +write_buffer_size = 8000 +weak_collision_order = 0 +suppress_deep_recoils = false +high_energy_free_flight_paths = false +num_threads = 1 +num_chunks = 1 +use_hdf5 = false +electronic_stopping_mode = "LOW_ENERGY_NONLOCAL" +mean_free_path_model = "LIQUID" +interaction_potential = [["KR_C"]] +scattering_integral = [["MENDENHALL_WELLER"]] +root_finder = [[{"NEWTON"={max_iterations=100, tolerance=1E-3}}]] +track_displacements = false +track_energy_losses = false + +[material_parameters] +energy_unit = "EV" +mass_unit = "AMU" +Eb = [ 0.0, 0.0,] +Es = [ 5.76, 0.0,] +Ec = [ 1.0, 1.0,] +Z = [ 5, 7,] +m = [ 10.811, 14,] +interaction_index = [0, 0] +surface_binding_model = "AVERAGE" + +[particle_parameters] +length_unit = "MICRON" +energy_unit = "EV" +mass_unit = "AMU" +N = [ 10000 ] +m = [ 1.008 ] +Z = [ 1 ] +E = [ 1000.0 ] +Ec = [ 1.0 ] +Es = [ 10.0 ] +pos = [ [ 0.0, 0.0, 0.0,] ] +dir = [ [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,] ] +interaction_index = [ 0 ] +particle_input_filename="" + +[geometry_input] +length_unit = "ANGSTROM" +electronic_stopping_correction_factor = 1.0 +densities = [0.065, 0.065] diff --git a/examples/layered_geometry.toml b/examples/layered_geometry.toml index 5f22a46..01da35b 100644 --- a/examples/layered_geometry.toml +++ b/examples/layered_geometry.toml @@ -45,10 +45,10 @@ Ec = [ 1.0,] Es = [ 0.0,] interaction_index = [0] pos = [ [ -1.7453292519934434e-8, 0.0, 0.0,],] -dir = [ [ 0.999, 0.0001, 0.0,],] +dir = [ [ 0.999, 0.001, 0.0,],] particle_input_filename = "" -[mesh_2d_input] +[geometry_input] 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,],] diff --git a/examples/multiple_interaction_potentials.toml b/examples/multiple_interaction_potentials.toml index f5a49ea..f10fa32 100644 --- a/examples/multiple_interaction_potentials.toml +++ b/examples/multiple_interaction_potentials.toml @@ -42,7 +42,7 @@ pos = [ [ -1.7453292519934434e-8, 0.0, 0.0,], [ -1.7453292519934434e-8, 0.0, 0.0 dir = [ [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,], [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,],] particle_input_filename = -[mesh_2d_input] +[geometry_input] energy_barrier_thickness = 2.2E-4 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,],] diff --git a/examples/titanium_dioxide_0D.toml b/examples/titanium_dioxide_0D.toml new file mode 100644 index 0000000..ff1e7f3 --- /dev/null +++ b/examples/titanium_dioxide_0D.toml @@ -0,0 +1,65 @@ +[options] +name = "2000.0eV_0.0001deg_He_TiO2" +track_trajectories = false +track_recoils = true +track_recoil_trajectories = false +write_buffer_size = 8000 +weak_collision_order = 3 +suppress_deep_recoils = false +high_energy_free_flight_paths = false +electronic_stopping_mode = "LOW_ENERGY_NONLOCAL" +mean_free_path_model = "LIQUID" +interaction_potential = [["KR_C"]] +scattering_integral = [["MENDENHALL_WELLER"]] +use_hdf5 = false +root_finder = [[{"NEWTON"={max_iterations=100, tolerance=1E-6}}]] +num_threads = 4 +num_chunks = 10 +track_displacements = false +track_energy_losses = false +energy_min = 0.0 +energy_max = 2000.0 +energy_num = 20 +angle_min = 0.0 +angle_max = 90.0 +angle_num = 20 +x_min = 0.0 +y_min = -0.1 +z_min = -0.1 +x_max = 0.1 +y_max = 0.1 +z_max = 0.1 +x_num = 20 +y_num = 20 +z_num = 20 + +[particle_parameters] +length_unit = "1e-6" +energy_unit = "EV" +mass_unit = "AMU" +N = [ 100000,] +m = [ 4.002602,] +Z = [ 2,] +E = [ 2000.0,] +Ec = [ 1.0,] +Es = [ 0.0,] +interaction_index = [0] +pos = [ [ 0.0, 0.0, 0.0,],] +dir = [ [ 0.5, 0.5, 0.0,],] +particle_input_filename = "" + +[geometry_input] +length_unit = "MICRON" +densities = [ 3.0e+10, 6.0e+10, 0.0, 0.0 ] +electronic_stopping_correction_factor = 1.0 + +[material_parameters] +energy_unit = "EV" +mass_unit = "AMU" +Eb = [ 3.0, 2.58, 3.0, 0.0,] +Es = [ 4.84, 2.58, 3.39, 4.72,] +Ec = [ 3.5, 2.0, 3.0, 1.5,] +Z = [ 22, 8, 13, 14,] +m = [ 47.867, 15.9994, 26.98, 28.08553,] +interaction_index = [0, 0, 0, 0] +surface_binding_model = "TARGET" diff --git a/src/bca.rs b/src/bca.rs index 0fea7b7..d1ed73f 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -1,5 +1,4 @@ use super::*; -use geo::{Closest}; #[cfg(any(feature = "cpr_rootfinder_openblas", feature = "cpr_rootfinder_netlib", feature = "cpr_rootfinder_intel_mkl"))] use rcpr::chebyshev::*; @@ -63,7 +62,7 @@ impl fmt::Display for BinaryCollisionResult { } /// This function takes a single particle, a material, and an options object and runs a binary-collision-approximation trajectory for that particle in that material, producing a final particle list that consists of the original ion and any material particles displaced thereby. -pub fn single_ion_bca(particle: particle::Particle, material: &material::Material, options: &Options) -> Vec { +pub fn single_ion_bca(particle: particle::Particle, material: &material::Material, options: &Options) -> Vec { let mut particles: Vec = Vec::new(); particles.push(particle); @@ -77,6 +76,8 @@ pub fn single_ion_bca(particle: particle::Particle, material: &material::Materia //Remove particle from top of vector as particle_1 let mut particle_1 = particles.pop().unwrap(); + //println!("Particle start Z: {} E: {} ({}, {}, {})", particle_1.Z, particle_1.E/Q, particle_1.pos.x/ANGSTROM, particle_1.pos.y/ANGSTROM, particle_1.pos.z/ANGSTROM); + //BCA loop while !particle_1.stopped & !particle_1.left { @@ -96,7 +97,7 @@ pub fn single_ion_bca(particle: particle::Particle, material: &material::Materia &binary_collision_geometry, &options); //If recoil location is inside, proceed with binary collision loop - if material.inside(particle_2.pos.x, particle_2.pos.y) & material.inside_energy_barrier(particle_1.pos.x, particle_1.pos.y) { + if material.inside(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z) & material.inside_energy_barrier(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z) { //Determine scattering angle from binary collision let binary_collision_result = bca::calculate_binary_collision(&particle_1, @@ -105,6 +106,8 @@ pub fn single_ion_bca(particle: particle::Particle, material: &material::Materia particle_1.pos.x, particle_2.pos.x, &binary_collision_geometry)) .unwrap(); + //println!("{}", binary_collision_result); + //Only use 0th order collision for local electronic stopping if k == 0 { normalized_distance_of_closest_approach = binary_collision_result.normalized_distance_of_closest_approach; @@ -113,7 +116,7 @@ pub fn single_ion_bca(particle: particle::Particle, material: &material::Materia } //Energy transfer to recoil - particle_2.E = binary_collision_result.recoil_energy - material.average_bulk_binding_energy(particle_2.pos.x, particle_2.pos.y); + particle_2.E = binary_collision_result.recoil_energy - material.average_bulk_binding_energy(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); particle_2.energy_origin = particle_2.E; //Accumulate asymptotic deflections for primary particle @@ -128,7 +131,7 @@ pub fn single_ion_bca(particle: particle::Particle, material: &material::Materia particle::rotate_particle(&mut particle_2, -binary_collision_result.psi_recoil, binary_collision_geometry.phi_azimuthal); - + particle_2.dir_old.x = particle_2.dir.x; particle_2.dir_old.y = particle_2.dir.y; particle_2.dir_old.z = particle_2.dir.z; @@ -148,23 +151,22 @@ pub fn single_ion_bca(particle: particle::Particle, material: &material::Materia let Ma: f64 = particle_1.m; let Mb: f64 = particle_2.m; - let n = material.total_number_density(particle_2.pos.x, particle_2.pos.y); + let n = material.total_number_density(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); //We just need the lindhard screening length here, so the particular potential is not important let a: f64 = interactions::screening_length(Za, Zb, InteractionPotential::MOLIERE); let reduced_energy: f64 = LINDHARD_REDUCED_ENERGY_PREFACTOR*a*Mb/(Ma+Mb)/Za/Zb*E; let estimated_range_of_recoils = (reduced_energy.powf(0.3) + 0.1).powf(3.)/n/a/a; - if let Closest::SinglePoint(p2) = material.closest_point(particle_2.pos.x, particle_2.pos.y) { - let dx = p2.x() - particle_2.pos.x; - let dy = p2.y() - particle_2.pos.y; - let distance_to_surface = (dx*dx + dy*dy).sqrt(); + let (x2, y2, z2) = material.closest_point(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); + let dx = x2 - particle_2.pos.x; + let dy = y2 - particle_2.pos.y; + let dz = z2 - particle_2.pos.z; + let distance_to_surface = (dx*dx + dy*dy + dz*dz).sqrt(); - if (distance_to_surface < estimated_range_of_recoils) & (particle_2.E > particle_2.Ec) { - particles.push(particle_2); - } - } else { - panic!("Numerical error: geometry algorithm failed to find distance from particle to surface.") + if (distance_to_surface < estimated_range_of_recoils) & (particle_2.E > particle_2.Ec) { + particles.push(particle_2); } + //If transferred energy > cutoff energy, add recoil to particle vector } else if options.track_recoils & (particle_2.E > particle_2.Ec) { particles.push(particle_2); @@ -180,13 +182,17 @@ pub fn single_ion_bca(particle: particle::Particle, material: &material::Materia bca::update_particle_energy(&mut particle_1, &material, distance_traveled, total_energy_loss, normalized_distance_of_closest_approach, strong_collision_Z, strong_collision_index, &options); + //println!("Particle finished collision loop Z: {} E: {} ({}, {}, {})", particle_1.Z, particle_1.E/Q, particle_1.pos.x/ANGSTROM, particle_1.pos.y/ANGSTROM, particle_1.pos.z/ANGSTROM); + + //println!("{} {} {}", energy_0/EV, energy_1/EV, (energy_1 - energy_0)/EV); //Check boundary conditions on leaving and stopping - material::boundary_condition_2D_planar(&mut particle_1, &material); + material::boundary_condition_planar(&mut particle_1, &material); //Set particle index to topmost particle particle_index = particles.len(); } + //println!("Particle stopped or left Z: {} E: {} ({}, {}, {})", particle_1.Z, particle_1.E/Q, particle_1.pos.x/ANGSTROM, particle_1.pos.y/ANGSTROM, particle_1.pos.z/ANGSTROM); particle_output.push(particle_1); } particle_output @@ -194,19 +200,18 @@ pub fn single_ion_bca(particle: particle::Particle, material: &material::Materia /// For a particle in a material, determine the mean free path and choose the azimuthal angle and impact parameter. /// The mean free path can be exponentially distributed (gaseous) or constant (amorphous solid/liquid). Azimuthal angles are chosen uniformly. Impact parameters are chosen for collision partners distributed uniformly on a disk of density-dependent radius. -pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle::Particle, material: &material::Material, options: &Options) -> Vec { +pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle::Particle, material: &material::Material, options: &Options) -> Vec { let x = particle_1.pos.x; let y = particle_1.pos.y; let z = particle_1.pos.z; - let mut mfp = material.mfp(x, y); + let mut mfp = material.mfp(x, y, z); let mut phis_azimuthal = Vec::with_capacity(options.weak_collision_order + 1); let mut binary_collision_geometries = Vec::with_capacity(options.weak_collision_order + 1); //Each weak collision gets its own aziumuthal angle in annuli around collision point - //azimuthal angle randomly selected (0..2pi) for k in 0..options.weak_collision_order + 1 { phis_azimuthal.push(2.*PI*rand::random::()); } @@ -214,11 +219,11 @@ pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle::Particle, m if options.high_energy_free_flight_paths { let Ma: f64 = particle_1.m; - let Mb: f64 = material.average_mass(x, y); + let Mb: f64 = material.average_mass(x, y, z); let Za: f64 = particle_1.Z; - let Zb: f64 = material.average_Z(x, y); - let n: &Vec = material.number_densities(x, y); - let ck: f64 = material.electronic_stopping_correction_factor(x, y); + let Zb: f64 = material.average_Z(x, y, z); + let n: &Vec = material.number_densities(x, y, z); + let ck: f64 = material.electronic_stopping_correction_factor(x, y, z); let E: f64 = particle_1.E; let Ec: f64 = particle_1.Ec; //We just need the Lindhard screening length here, so the particular potential is not important @@ -232,19 +237,19 @@ pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle::Particle, m //Free flight path formulation here described in SRIM textbook chapter 7, and Eckstein 1991 7.5.2 let ep = (reduced_energy*reduced_energy_min).sqrt(); let mut pmax = a/(ep + ep.sqrt() + 0.125*ep.powf(0.1)); - let mut ffp = 1./(material.total_number_density(x, y)*pmax*pmax*PI); + let mut ffp = 1./(material.total_number_density(x, y, z)*pmax*pmax*PI); let stopping_powers = material.electronic_stopping_cross_sections(particle_1, ElectronicStoppingMode::INTERPOLATED); - let delta_energy_electronic = stopping_powers.iter().zip(material.number_densities(x, y)).map(|(&se, number_density)| se*number_density).sum::()*ffp*ck; + let delta_energy_electronic = stopping_powers.iter().zip(material.number_densities(x, y, z)).map(|(&se, number_density)| se*number_density).sum::()*ffp*ck; //If losing too much energy, scale free-flight-path down //5 percent limit set in original TRIM paper, Biersack and Haggmark 1980 if delta_energy_electronic > 0.05*E { ffp *= 0.05*E/delta_energy_electronic; - pmax = (1./(material.total_number_density(x, y)*PI*ffp)).sqrt() + pmax = (1./(material.total_number_density(x, y, z)*PI*ffp)).sqrt() } - //If free-flight-path less than the interatomic spacing, revert to solid model + //If free-flight-path less than the interatomic spacing, revert to amorphous solid model //Mentioned in Eckstein 1991, Ziegler, Biersack, and Ziegler 2008 (SRIM textbook 7-8) if ffp < mfp { @@ -325,7 +330,7 @@ pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle::Particle, m } /// For a particle in a material, and for a particular binary collision geometry, choose a species for the collision partner. -pub fn choose_collision_partner(particle_1: &particle::Particle, material: &material::Material, binary_collision_geometry: &BinaryCollisionGeometry, options: &Options) -> (usize, particle::Particle) { +pub fn choose_collision_partner(particle_1: &particle::Particle, material: &material::Material, binary_collision_geometry: &BinaryCollisionGeometry, options: &Options) -> (usize, particle::Particle) { let x = particle_1.pos.x; let y = particle_1.pos.y; let z = particle_1.pos.z; @@ -348,7 +353,7 @@ pub fn choose_collision_partner(particle_1: &particle::Particle, material: &mate let z_recoil: f64 = z + mfp*cosz + impact_parameter*(sinphi*cosy - cosphi*cosx*cosz)/sinx; //Choose recoil Z, M - let (species_index, Z_recoil, M_recoil, Ec_recoil, Es_recoil, interaction_index) = material.choose(x_recoil, y_recoil); + let (species_index, Z_recoil, M_recoil, Ec_recoil, Es_recoil, interaction_index) = material.choose(x_recoil, y_recoil, z_recoil); return (species_index, particle::Particle::new( @@ -405,7 +410,7 @@ fn distance_of_closest_approach(particle_1: &particle::Particle, particle_2: &pa } /// Update a particle's energy following nuclear and electronic interactions of a single BCA step. -pub fn update_particle_energy(particle_1: &mut particle::Particle, material: &material::Material, distance_traveled: f64, +pub fn update_particle_energy(particle_1: &mut particle::Particle, material: &material::Material, distance_traveled: f64, recoil_energy: f64, x0: f64, strong_collision_Z: f64, strong_collision_index: usize, options: &Options) { //If particle energy drops below zero before electronic stopping calcualtion, it produces NaNs @@ -417,12 +422,13 @@ pub fn update_particle_energy(particle_1: &mut particle::Particle, material: &ma let x = particle_1.pos.x; let y = particle_1.pos.y; + let z = particle_1.pos.z; - if material.inside(x, y) { + if material.inside(x, y, z) { let interaction_potential = options.interaction_potential[particle_1.interaction_index][material.interaction_index[strong_collision_index]]; let electronic_stopping_powers = material.electronic_stopping_cross_sections(particle_1, options.electronic_stopping_mode); - let n = material.number_densities(x, y); + let n = material.number_densities(x, y, z); let delta_energy = match options.electronic_stopping_mode { ElectronicStoppingMode::INTERPOLATED => electronic_stopping_powers.iter().zip(n).map(|(se, number_density)| se*number_density).collect::>().iter().sum::()*distance_traveled, diff --git a/src/enums.rs b/src/enums.rs index a1dca03..0da78e5 100644 --- a/src/enums.rs +++ b/src/enums.rs @@ -1,5 +1,16 @@ use super::*; +pub enum MaterialType { + MESH0D(material::Material), + MESH2D(material::Material) +} + +#[derive(Deserialize)] +pub enum GeometryType { + MESH0D, + MESH2D, +} + /// Mode of electronic stopping to use. #[derive(Deserialize, PartialEq, Clone, Copy)] pub enum ElectronicStoppingMode { diff --git a/src/input.rs b/src/input.rs index de4d91a..4188b1c 100644 --- a/src/input.rs +++ b/src/input.rs @@ -1,18 +1,85 @@ use super::*; +pub trait InputFile: GeometryInput { + fn new(string: &str) -> Self; + fn get_options(&self) -> &Options; + fn get_material_parameters(&self) -> &material::MaterialParameters; + fn get_particle_parameters(&self) -> &particle::ParticleParameters; + fn get_geometry_input(&self) -> &::GeometryInput; +} + +pub trait GeometryInput { + type GeometryInput; +} /// Rustbca's internal representation of an input file. -#[derive(Deserialize)] -pub struct Input { +#[derive(Deserialize, Clone)] +pub struct Input2D { pub options: Options, pub material_parameters: material::MaterialParameters, pub particle_parameters: particle::ParticleParameters, - pub mesh_2d_input: mesh::Mesh2DInput, + pub geometry_input: mesh::Mesh2DInput, +} + +impl GeometryInput for Input2D { + type GeometryInput = mesh::Mesh2DInput; +} + +impl InputFile for Input2D { + + fn new(string: &str) -> Input2D { + toml::from_str(string).unwrap() + } + + fn get_options(&self) -> &Options{ + &self.options + } + fn get_material_parameters(&self) -> &material::MaterialParameters{ + &self.material_parameters + } + fn get_particle_parameters(&self) -> &particle::ParticleParameters{ + &self.particle_parameters + } + fn get_geometry_input(&self) -> &Self::GeometryInput{ + &self.geometry_input + } +} + +#[derive(Deserialize, Clone)] +pub struct Input0D { + pub options: Options, + pub material_parameters: material::MaterialParameters, + pub particle_parameters: particle::ParticleParameters, + pub geometry_input: mesh::Mesh0DInput, +} + +impl GeometryInput for Input0D { + type GeometryInput = mesh::Mesh0DInput; +} + +impl InputFile for Input0D { + + fn new(string: &str) -> Input0D { + toml::from_str(string).expect("Could not parse TOML file.") + } + + fn get_options(&self) -> &Options{ + &self.options + } + fn get_material_parameters(&self) -> &material::MaterialParameters{ + &self.material_parameters + } + fn get_particle_parameters(&self) ->& particle::ParticleParameters{ + &self.particle_parameters + } + fn get_geometry_input(&self) -> &Self::GeometryInput{ + &self.geometry_input + } } /// Rustbca's internal representation of the simulation-level options. #[cfg(not(feature = "distributions"))] -#[derive(Deserialize)] +#[derive(Deserialize, Clone)] pub struct Options { pub name: String, pub track_trajectories: bool, @@ -35,7 +102,7 @@ pub struct Options { } #[cfg(feature = "distributions")] -#[derive(Deserialize)] +#[derive(Deserialize, Clone)] pub struct Options { pub name: String, pub track_trajectories: bool, @@ -72,14 +139,16 @@ pub struct Options { pub z_num: usize, } -pub fn input() -> (Vec, material::Material, Options, OutputUnits){ +pub fn input() -> (Vec, material::Material, Options, OutputUnits) +where ::InputFileFormat: Deserialize<'static> + 'static, ::GeometryInput: Clone { let args: Vec = env::args().collect(); - let input_file = match args.len() { - 1 => "input.toml".to_string(), - 2 => args[1].clone(), - _ => panic!("Too many command line arguments. RustBCA accepts 0 (use 'input.toml') or 1 (input file name).") + let (input_file, geometry_type) = match args.len() { + 1 => ("input.toml".to_string(), GeometryType::MESH2D), + 2 => (args[1].clone(), GeometryType::MESH2D), + 3 => (args[2].clone(), match args[1].as_str() { "0D" => GeometryType::MESH0D, "2D" => GeometryType::MESH2D, _ => panic!("Unimplemented geometry {}.", args[1].clone()) }), + _ => panic!("Too many command line arguments. RustBCA accepts 0 (use 'input.toml') 1 () or 2 ( )"), }; //Read input file, convert to string, and open with toml @@ -91,13 +160,14 @@ pub fn input() -> (Vec, material::Material, Options, Ou .open(&input_file) .expect(format!("Input errror: could not open input file {}.", &input_file).as_str()); file.read_to_string(&mut input_toml).context("Could not convert TOML file to string.").unwrap(); - let input: Input = toml::from_str(&input_toml) - .expect("Input error: failed to parse TOML file. Check that all floats have terminal 0s and that there are no missing/extra delimiters."); + + let input: ::InputFileFormat = InputFile::new(&input_toml); //Unpack toml information into structs - let material = material::Material::new(input.material_parameters, input.mesh_2d_input); - let options = input.options; - let particle_parameters = input.particle_parameters; + let options = (*input.get_options()).clone(); + let particle_parameters = (*input.get_particle_parameters()).clone(); + let material_parameters = (*input.get_material_parameters()).clone(); + let material: material::Material = material::Material::::new(&material_parameters, input.get_geometry_input()); //Ensure nonsensical threads/chunks options crash on input assert!(options.num_threads > 0, "Input error: num_threads must be greater than zero."); diff --git a/src/main.rs b/src/main.rs index abd4aba..c34d906 100644 --- a/src/main.rs +++ b/src/main.rs @@ -62,85 +62,108 @@ pub mod structs; pub use crate::enums::*; pub use crate::consts::*; pub use crate::structs::*; -pub use crate::input::{Input, Options}; +pub use crate::input::{Input2D, Input0D, Options, InputFile, GeometryInput}; pub use crate::output::{OutputUnits}; +pub use crate::mesh::{Geometry, GeometryElement, Mesh0D, Mesh2D}; -fn main() { - //Open and process input_file - let (particle_input_array, material, options, output_units) = input::input(); - - println!("Processing {} ions...", particle_input_array.len()); - - let total_count: u64 = particle_input_array.len() as u64; - assert!(total_count/options.num_chunks > 0, "Input error: chunk size == 0 - reduce num_chunks or increase particle count."); - - #[cfg(not(feature = "no_list_output"))] - let mut output_list_streams = output::open_output_lists(&options); - - let mut summary = output::SummaryPerSpecies::new(&options); - - #[cfg(feature = "distributions")] - let mut distributions = output::Distributions::new(&options); - - //Initialize threads with rayon - println!("Initializing with {} threads...", options.num_threads); - if options.num_threads > 1 {let pool = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global().unwrap();}; - - //Create and configure progress bar - let bar: ProgressBar = ProgressBar::new(total_count); - bar.set_style(ProgressStyle::default_bar() - .template("[{elapsed_precise}][{bar:40.cyan/blue}][{eta_precise}] {percent}%") - .progress_chars("#>-")); - - //Main loop - for (chunk_index, particle_input_chunk) in particle_input_array.chunks((total_count/options.num_chunks) as usize).enumerate() { - - let mut finished_particles: Vec = Vec::new(); - - if options.num_threads > 1 { - // BCA loop is implemented as parallelized extension of a per-chunk initially empty - // finished particle array via map from particle -> finished particles via BCA - finished_particles.par_extend( - particle_input_chunk.into_par_iter() - .map(|particle_input| { - bar.tick(); - bar.inc(1); - bca::single_ion_bca(particle::Particle::from_input(*particle_input, &options), &material, &options) - }).flatten() - ); - } else { - finished_particles.extend( - particle_input_chunk.iter() - .map(|particle_input| { - bar.tick(); - bar.inc(1); - bca::single_ion_bca(particle::Particle::from_input(*particle_input, &options), &material, &options) - }).flatten() - ); - } +fn physics_loop(particle_input_array: Vec, material: material::Material, options: Options, output_units: OutputUnits) { - // Process this chunk of finished particles for output - for particle in finished_particles { + println!("Processing {} ions...", particle_input_array.len()); - summary.update(&particle); - - #[cfg(feature = "distributions")] - distributions.update(&particle, &output_units); + let total_count: u64 = particle_input_array.len() as u64; + assert!(total_count/options.num_chunks > 0, "Input error: chunk size == 0 - reduce num_chunks or increase particle count."); + #[cfg(not(feature = "no_list_output"))] + let mut output_list_streams = output::open_output_lists(&options); + + let mut summary = output::SummaryPerSpecies::new(&options); + + #[cfg(feature = "distributions")] + let mut distributions = output::Distributions::new(&options); + + //Initialize threads with rayon + println!("Initializing with {} threads...", options.num_threads); + if options.num_threads > 1 {let pool = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global().unwrap();}; + + //Create and configure progress bar + let bar: ProgressBar = ProgressBar::new(total_count); + bar.set_style(ProgressStyle::default_bar() + .template("[{elapsed_precise}][{bar:40.cyan/blue}][{eta_precise}] {percent}%") + .progress_chars("#>-")); + + //Main loop + for (chunk_index, particle_input_chunk) in particle_input_array.chunks((total_count/options.num_chunks) as usize).enumerate() { + + let mut finished_particles: Vec = Vec::new(); + + if options.num_threads > 1 { + // BCA loop is implemented as parallelized extension of a per-chunk initially empty + // finished particle array via map from particle -> finished particles via BCA + finished_particles.par_extend( + particle_input_chunk.into_par_iter() + .map(|particle_input| { + bar.tick(); + bar.inc(1); + bca::single_ion_bca(particle::Particle::from_input(*particle_input, &options), &material, &options) + }).flatten() + ); + } else { + finished_particles.extend( + particle_input_chunk.iter() + .map(|particle_input| { + bar.tick(); + bar.inc(1); + bca::single_ion_bca(particle::Particle::from_input(*particle_input, &options), &material, &options) + }).flatten() + ); + } + + // Process this chunk of finished particles for output + for particle in finished_particles { + + summary.update(&particle); + + #[cfg(feature = "distributions")] + distributions.update(&particle, &output_units); + + #[cfg(not(feature = "no_list_output"))] + output::output_lists(&mut output_list_streams, particle, &options, &output_units); + + } + //Flush all file streams before dropping to ensure all data is written #[cfg(not(feature = "no_list_output"))] - output::output_lists(&mut output_list_streams, particle, &options, &output_units); - + output::output_list_flush(&mut output_list_streams); } - //Flush all file streams before dropping to ensure all data is written - #[cfg(not(feature = "no_list_output"))] - output::output_list_flush(&mut output_list_streams); - } - summary.print(&options, &output_units); + summary.print(&options, &output_units); - //Write distributions to file - #[cfg(feature = "distributions")] - distributions.print(&options); + //Write distributions to file + #[cfg(feature = "distributions")] + distributions.print(&options); + + println!("Finished!"); +} + +fn main() { + + let args: Vec = env::args().collect(); + + let (input_file, geometry_type) = match args.len() { + 1 => ("input.toml".to_string(), GeometryType::MESH2D), + 2 => (args[1].clone(), GeometryType::MESH2D), + 3 => (args[2].clone(), match args[1].as_str() { "0D" => GeometryType::MESH0D, "2D" => GeometryType::MESH2D, _ => panic!("Unimplemented geometry {}.", args[1].clone()) }), + _ => panic!("Too many command line arguments. RustBCA accepts 0 (use 'input.toml') 1 () or 2 ( )"), + }; + + match geometry_type { + GeometryType::MESH0D => { + let (particle_input_array, material, options, output_units) = input::input::(); + physics_loop::(particle_input_array, material, options, output_units); + }, + GeometryType::MESH2D => { + let (particle_input_array, material, options, output_units) = input::input::(); + physics_loop::(particle_input_array, material, options, output_units); + }, + } - println!("Finished!"); } diff --git a/src/material.rs b/src/material.rs index 86ef326..174f3e0 100644 --- a/src/material.rs +++ b/src/material.rs @@ -1,11 +1,7 @@ use super::*; -use geo::algorithm::contains::Contains; -use geo::algorithm::closest_point::ClosestPoint; -use geo::{point, Closest}; - /// Holds material input parameters from [material_params]. -#[derive(Deserialize)] +#[derive(Deserialize, Clone)] pub struct MaterialParameters { pub energy_unit: String, pub mass_unit: String, @@ -19,20 +15,19 @@ pub struct MaterialParameters { } /// Material in rustbca. Includes the material properties and the mesh that defines the material geometry. -pub struct Material { +pub struct Material { pub m: Vec, pub Z: Vec, pub Eb: Vec, pub Es: Vec, pub Ec: Vec, pub interaction_index: Vec, - pub mesh_2d: mesh::Mesh2D, + pub geometry: Box, pub surface_binding_model: SurfaceBindingModel - } -impl Material { - /// Constructs a new material object from a material parameters object and a mesh_2d_input object. - pub fn new(material_parameters: MaterialParameters, mesh_2d_input: mesh::Mesh2DInput) -> Material { +impl Material { + + pub fn new(material_parameters: &MaterialParameters, geometry_input: &<::InputFileFormat as GeometryInput>::GeometryInput) -> Material { let energy_unit: f64 = match material_parameters.energy_unit.as_str() { "EV" => EV, @@ -56,33 +51,33 @@ impl Material { Material { m: material_parameters.m.iter().map(|&i| i*mass_unit).collect(), - Z: material_parameters.Z, + Z: material_parameters.Z.clone(), Eb: material_parameters.Eb.iter().map(|&i| i*energy_unit).collect(), Es: material_parameters.Es.iter().map(|&i| i*energy_unit).collect(), Ec: material_parameters.Ec.iter().map(|&i| i*energy_unit).collect(), - interaction_index: material_parameters.interaction_index, - mesh_2d: mesh::Mesh2D::new(mesh_2d_input), + interaction_index: material_parameters.interaction_index.clone(), surface_binding_model: material_parameters.surface_binding_model, + geometry: Box::new(T::new(geometry_input)), } } /// Gets concentrations of triangle that contains or is nearest to (x, y) - pub fn get_concentrations(&self, x: f64, y: f64) -> Vec { + pub fn get_concentrations(&self, x: f64, y: f64, z: f64) -> Vec { - let total_number_density: f64 = if self.inside(x, y) { - self.mesh_2d.get_total_density(x,y) + let total_number_density: f64 = if self.inside(x, y, z) { + self.geometry.get_total_density(x, y, z) } else { - self.mesh_2d.nearest_cell_to(x, y).densities.iter().sum::() + self.geometry.get_densities_nearest_to(x, y, z).iter().sum::() }; - return self.mesh_2d.get_densities(x, y).iter().map(|&i| i / total_number_density).collect(); + return self.geometry.get_densities(x, y, z).iter().map(|&i| i / total_number_density).collect(); } /// Gets cumulative concentrations of triangle that contains or is nearest to (x, y). /// Used for weighting the random choice of one of the constituent species. - pub fn get_cumulative_concentrations(&self, x: f64, y: f64) -> Vec { + pub fn get_cumulative_concentrations(&self, x: f64, y: f64, z: f64) -> Vec { let mut sum: f64 = 0.; - let concentrations = self.mesh_2d.get_concentrations(x, y); + let concentrations = self.geometry.get_concentrations(x, y, z); let mut cumulative_concentrations = Vec::with_capacity(concentrations.len()); for concentration in concentrations { sum += concentration; @@ -92,97 +87,74 @@ impl Material { } /// Determines whether (x, y) is inside the material. - pub fn inside(&self, x: f64, y: f64) -> bool { - return self.mesh_2d.inside(x, y); + pub fn inside(&self, x: f64, y: f64, z: f64) -> bool { + return self.geometry.inside(x, y, z); } /// Gets electronic stopping correction factor for LS and OR - pub fn electronic_stopping_correction_factor(&self, x: f64, y: f64) -> f64 { - if self.inside(x, y) { - return self.mesh_2d.get_ck(x, y); + pub fn electronic_stopping_correction_factor(&self, x: f64, y: f64, z: f64) -> f64 { + if self.inside(x, y, z) { + return self.geometry.get_ck(x, y, z); } else { - return self.mesh_2d.nearest_cell_to(x, y).electronic_stopping_correction_factor; + return self.geometry.get_ck_nearest_to(x, y, z); } } /// Determines the local mean free path from the formula sum(n(x, y))^(-1/3) - pub fn mfp(&self, x: f64, y: f64) -> f64 { - return self.total_number_density(x, y).powf(-1./3.); + pub fn mfp(&self, x: f64, y: f64, z: f64) -> f64 { + return self.total_number_density(x, y, z).powf(-1./3.); } /// Total number density of triangle that contains or is nearest to (x, y). - pub fn total_number_density(&self, x: f64, y: f64) -> f64 { - if self.inside(x, y) { - return self.mesh_2d.get_densities(x, y).iter().sum::(); + pub fn total_number_density(&self, x: f64, y: f64, z: f64) -> f64 { + if self.inside(x, y, z) { + return self.geometry.get_densities(x, y, z).iter().sum::(); } else { - return self.mesh_2d.nearest_cell_to(x, y).densities.iter().sum::(); + return self.geometry.get_densities_nearest_to(x, y, z).iter().sum::(); } - //return self.n.iter().sum::(); } /// Lists number density of each species of triangle that contains or is nearest to (x, y). - pub fn number_densities(&self, x: f64, y: f64) -> &Vec { - if self.inside(x, y) { - return &self.mesh_2d.get_densities(x, y); + pub fn number_densities(&self, x: f64, y: f64, z: f64) -> &Vec { + if self.inside(x, y, z) { + return &self.geometry.get_densities(x, y, z); } else { - return &self.mesh_2d.nearest_cell_to(x, y).densities; + return &self.geometry.get_densities_nearest_to(x, y, z); } - //return self.n.iter().sum::(); } /// Determines whether a point (x, y) is inside the energy barrier of the material. - pub fn inside_energy_barrier(&self, x: f64, y: f64) -> bool { - - if self.mesh_2d.inside(x, y) { - true - } else { - let nearest_cell = self.mesh_2d.nearest_cell_to(x, y); - let distance = nearest_cell.distance_to(x, y); - distance < self.mesh_2d.energy_barrier_thickness - } - //let p = point!(x: x, y: y); - //return self.energy_surface.contains(&p); + pub fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { + self.geometry.inside_energy_barrier(x, y, z) } /// Determines whether a point (x, y) is inside the simulation boundary. - pub fn inside_simulation_boundary(&self, x:f64, y: f64) -> bool { - let p = point!(x: x, y: y); - return self.mesh_2d.simulation_boundary.contains(&p); + pub fn inside_simulation_boundary(&self, x:f64, y: f64, z: f64) -> bool { + return self.geometry.inside_simulation_boundary(x, y, z); } /// Finds the closest point on the material boundary to the point (x, y). - pub fn closest_point(&self, x: f64, y: f64) -> Closest { - let p = point!(x: x, y: y); - return self.mesh_2d.boundary.closest_point(&p); - //return self.surface.closest_point(&p) - } - - /// Finds the closest point on the simulation boundary to the point (x, y). - pub fn closest_point_on_simulation_surface(&self, x: f64, y: f64) -> Closest { - let p = point!(x: x, y: y); - return self.mesh_2d.simulation_boundary.closest_point(&p) + pub fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + self.geometry.closest_point(x, y, z) } /// Finds the average, concentration-weighted atomic number, Z_effective, of the triangle that contains or is nearest to (x, y). - pub fn average_Z(&self, x: f64, y: f64) -> f64 { - let concentrations = self.mesh_2d.get_concentrations(x, y); + pub fn average_Z(&self, x: f64, y: f64, z: f64) -> f64 { + let concentrations = self.geometry.get_concentrations(x, y, z); return self.Z.iter().zip(concentrations).map(|(charge, concentration)| charge*concentration).collect::>().iter().sum(); - //return self.Z.iter().sum::()/self.Z.len() as f64; } /// Finds the average, concentration-weighted atomic mass, m_effective, of the triangle that contains or is nearest to (x, y). - pub fn average_mass(&self, x: f64, y: f64) -> f64 { - let concentrations = self.mesh_2d.get_concentrations(x, y); + pub fn average_mass(&self, x: f64, y: f64, z: f64) -> f64 { + let concentrations = self.geometry.get_concentrations(x, y, z); return self.m.iter().zip(concentrations).map(|(mass, concentration)| mass*concentration).collect::>().iter().sum(); - //return self.m.iter().sum::()/self.m.len() as f64; } /// Finds the average, concentration-weighted bulk binding energy of the triangle that contains or is nearest to (x, y). - pub fn average_bulk_binding_energy(&self, x: f64, y: f64) -> f64 { + pub fn average_bulk_binding_energy(&self, x: f64, y: f64, z: f64) -> f64 { //returns average bulk binding energy - let concentrations = self.mesh_2d.get_concentrations(x, y); + let concentrations = self.geometry.get_concentrations(x, y, z); return self.Eb.iter().zip(concentrations).map(|(bulk_binding_energy, concentration)| bulk_binding_energy*concentration).collect::>().iter().sum(); - //return self.Eb.iter().sum::()/self.Eb.len() as f64; } /// Finds the concentration-dependent surface binding energy of the triangle that contains or is nearest to (x, y). @@ -190,8 +162,8 @@ impl Material { /// 1. INDIVIDUAL: the surface binding energies are set individually for each species, as Es. /// 2. TARGET: the surface binding energy is calculated as the local concentration-weighted average of the target surface binding energies, unless the particle has Es = 0, in which case it is 0. /// 3. AVERAGE: the surface binding energy is the average of the particle surface binding energy and the TARGET surface binding energy, unless either is 0 in which case it is 0. - pub fn actual_surface_binding_energy(&self, particle: &particle::Particle, x: f64, y: f64) -> f64 { - let concentrations = self.mesh_2d.get_concentrations(x, y); + pub fn actual_surface_binding_energy(&self, particle: &particle::Particle, x: f64, y: f64, z: f64) -> f64 { + let concentrations = self.geometry.get_concentrations(x, y, z); match self.surface_binding_model { SurfaceBindingModel::INDIVIDUAL => particle.Es, @@ -227,9 +199,9 @@ impl Material { } ///Choose the parameters of a target atom as a concentration-weighted random draw from the species in the triangle that contains or is nearest to (x, y). - pub fn choose(&self, x: f64, y: f64) -> (usize, f64, f64, f64, f64, usize) { + pub fn choose(&self, x: f64, y: f64, z: f64) -> (usize, f64, f64, f64, f64, usize) { let random_number: f64 = rand::random::(); - let cumulative_concentrations = self.get_cumulative_concentrations(x, y); + let cumulative_concentrations = self.get_cumulative_concentrations(x, y, z); for (component_index, cumulative_concentration) in cumulative_concentrations.iter().enumerate() { if random_number < *cumulative_concentration { @@ -253,15 +225,15 @@ impl Material { let x = particle_1.pos.x; let y = particle_1.pos.y; - let ck = self.electronic_stopping_correction_factor(x, y); + let z = particle_1.pos.y; + let ck = self.electronic_stopping_correction_factor(x, y, z); - for (n, Zb) in self.number_densities(x, y).iter().zip(&self.Z) { - //let n = self.number_density(particle_1.pos.x, particle_1.pos.y); - //let Zb = self.Z_eff(particle_1.pos.x, particle_1.pos.y); + for (n, Zb) in self.number_densities(x, y, z).iter().zip(&self.Z) { let beta = (1. - (1. + E/Ma/C.powf(2.)).powf(-2.)).sqrt(); let v = beta*C; + // This term is an empirical fit to the mean ionization potential let I0 = match *Zb < 13. { true => 12. + 7./Zb, false => 9.76 + 58.5*Zb.powf(-1.19), @@ -301,22 +273,25 @@ impl Material { /// Calculate the effects of the planar surface binding potential of a material on a particle. /// These effects include surface reflection and refraction of particles with non-zero surface binding energies. -pub fn surface_binding_energy(particle_1: &mut particle::Particle, material: &material::Material) { +pub fn surface_binding_energy(particle_1: &mut particle::Particle, material: &material::Material) { let x = particle_1.pos.x; let y = particle_1.pos.y; + let z = particle_1.pos.z; let x_old = particle_1.pos_old.x; let y_old = particle_1.pos_old.y; + let z_old = particle_1.pos_old.z; let cosx = particle_1.dir.x; let cosy = particle_1.dir.y; let cosz = particle_1.dir.z; let E = particle_1.E; //Actual surface binding energies - let Es = material.actual_surface_binding_energy(particle_1, x_old, y_old); + let Es = material.actual_surface_binding_energy(particle_1, x_old, y_old, z_old); + //println!("Actual Es: {}", Es); let Ec = particle_1.Ec; - let inside_now = material.inside_energy_barrier(x, y); - let inside_old = material.inside_energy_barrier(x_old, y_old); + let inside_now = material.inside_energy_barrier(x, y, z); + let inside_old = material.inside_energy_barrier(x_old, y_old, z_old); let leaving = !inside_now & inside_old; let entering = inside_now & !inside_old; @@ -324,75 +299,68 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, material: &ma if entering { if particle_1.backreflected { particle_1.backreflected = false; - } else if let Closest::SinglePoint(p2) = material.closest_point(x, y) { - let dx = p2.x() - x; - let dy = p2.y() - y; - let mag = (dx*dx + dy*dy).sqrt(); + } else { + let (x2, y2, z2) = material.closest_point(x, y, z); + let dx = x2 - x; + let dy = y2 - y; + let dz = z2 - z; + let mag = (dx*dx + dy*dy + dz*dz).sqrt(); - let costheta = dx*cosx/mag + dy*cosy/mag; + let costheta = dx*cosx/mag + dy*cosy/mag + dz*cosz/mag; particle_1.E += Es; - //Surface refraction via Snell's law + //Surface refraction through energy barrier via Snell's law let delta_theta = particle::refraction_angle(costheta, E, E + Es); particle::rotate_particle(particle_1, delta_theta, 0.); - - //particle_1.add_trajectory(); - } else { - panic!("Numerical error: surface boundary algorithm encountered an error. Check geometry."); } } if leaving { - if let Closest::SinglePoint(p2) = material.closest_point(x, y) { - - let dx = p2.x() - x; - let dy = p2.y() - y; - let mag = (dx*dx + dy*dy).sqrt(); - let costheta = dx*cosx/mag + dy*cosy/mag; - let leaving_energy = E*costheta*costheta; + let (x2, y2, z2) = material.closest_point(x, y, z); + let dx = x2 - x; + let dy = y2 - y; + let dz = z2 - z; + let mag = (dx*dx + dy*dy + dz*dz).sqrt(); + let costheta = dx*cosx/mag + dy*cosy/mag + dz*cosz/mag; + let leaving_energy = E*costheta*costheta; - if costheta < 0. { - if leaving_energy > Es { + if costheta < 0. { + if leaving_energy > Es { - //particle_1.left = true; - particle_1.E += -Es; + particle_1.E += -Es; - //Surface refraction via Snell's law - let delta_theta = particle::refraction_angle(costheta, E, E - Es); - particle::rotate_particle(particle_1, delta_theta, 0.); + //Surface refraction via Snell's law + let delta_theta = particle::refraction_angle(costheta, E, E - Es); + particle::rotate_particle(particle_1, delta_theta, 0.); - //particle_1.add_trajectory(); + } else { - } else { - - //Specular reflection at local surface normal - particle_1.dir.x = -2.*(costheta)*dx/mag + cosx; - particle_1.dir.y = -2.*(costheta)*dy/mag + cosy; + //Specular reflection at local surface normal + particle_1.dir.x = -2.*(costheta)*dx/mag + cosx; + particle_1.dir.y = -2.*(costheta)*dy/mag + cosy; + particle_1.dir.z = -2.*(costheta)*dz/mag + cosz; - particle_1.backreflected = true; - //particle_1.add_trajectory(); - } + particle_1.backreflected = true; } - } else { - panic!("Numerical error: surface boundary algorithm encountered an error. Check geometry."); } } } /// Apply the boundary conditions of a material on a particle, including stopping, leaving, and reflection/refraction by/through the surface binding potential. -pub fn boundary_condition_2D_planar(particle_1: &mut particle::Particle, material: &material::Material) { +pub fn boundary_condition_planar(particle_1: &mut particle::Particle, material: &material::Material) { let x = particle_1.pos.x; let y = particle_1.pos.y; + let z = particle_1.pos.z; let E = particle_1.E; let Ec = particle_1.Ec; - if !material.inside_simulation_boundary(x, y) { + if !material.inside_simulation_boundary(x, y, z) { particle_1.left = true; particle_1.add_trajectory(); } - if (E < Ec) & !particle_1.left & material.inside_energy_barrier(x, y) { + if (E < Ec) & !particle_1.left & material.inside_energy_barrier(x, y, z) { particle_1.stopped = true; particle_1.add_trajectory(); } diff --git a/src/mesh.rs b/src/mesh.rs index af1ded2..6985696 100644 --- a/src/mesh.rs +++ b/src/mesh.rs @@ -1,10 +1,134 @@ use super::*; + use geo::algorithm::contains::Contains; -use geo::{Polygon, LineString, Point}; +use geo::{Polygon, LineString, Point, point, Closest}; +use geo::algorithm::closest_point::ClosestPoint; + +///Trait for a Geometry object - all forms of geometry must implement these traits to be used +pub trait Geometry: GeometryInput { + + type InputFileFormat: InputFile + Clone; + + fn new(input: &<::InputFileFormat as GeometryInput>::GeometryInput) -> Self; + fn get_densities(&self, x: f64, y: f64, z: f64) -> &Vec; + fn get_ck(&self, x: f64, y: f64, z: f64) -> f64; + fn get_total_density(&self, x: f64, y: f64, z: f64) -> f64; + fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec; + fn get_densities_nearest_to(&self, x: f64, y: f64, z: f64) -> &Vec; + fn get_ck_nearest_to(&self, x: f64, y: f64, z: f64) -> f64; + fn inside(&self, x: f64, y: f64, z: f64) -> bool; + 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 get_energy_barrier_thickness(&self) -> f64; + fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64); +} + +pub trait GeometryElement: Clone { + fn contains(&self, x: f64, y: f64, z: f64) -> bool; + fn distance_to(&self, x: f64, y: f64, z: f64) -> f64; + fn get_densities(&self) -> &Vec; + fn get_concentrations(&self) -> &Vec; + fn get_electronic_stopping_correction_factor(&self) -> f64; +} + +#[derive(Deserialize, Clone)] +pub struct Mesh0DInput { + pub length_unit: String, + pub densities: Vec, + pub electronic_stopping_correction_factor: f64, +} + +#[derive(Clone)] +pub struct Mesh0D { + pub densities: Vec, + pub concentrations: Vec, + pub electronic_stopping_correction_factor: f64, + pub energy_barrier_thickness: f64, +} + +impl GeometryInput for Mesh0D { + type GeometryInput = Mesh0DInput; +} + +impl Geometry for Mesh0D { + + type InputFileFormat = Input0D; + + fn new(input: &::GeometryInput) -> Mesh0D { + + let length_unit: f64 = match input.length_unit.as_str() { + "MICRON" => MICRON, + "CM" => CM, + "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 electronic_stopping_correction_factor = input.electronic_stopping_correction_factor; + + let densities: Vec = input.densities.iter().map(|element| element/(length_unit).powf(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::>(); + + Mesh0D { + densities, + concentrations, + electronic_stopping_correction_factor, + energy_barrier_thickness + } + } + + 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 get_densities_nearest_to(&self, x: f64, y: f64, z: f64) -> &Vec { + &self.densities + } + fn get_ck_nearest_to(&self, x: f64, y: f64, z: f64) -> f64 { + self.electronic_stopping_correction_factor + } + fn inside(&self, x: f64, y: f64, z: f64) -> bool { + x > 0.0 + } + + fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool { + //println!("x: {} energy_barrier_thickness: {}", x/ANGSTROM, self.energy_barrier_thickness/ANGSTROM); + x > -10.*self.energy_barrier_thickness + } + + fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { + x > -self.energy_barrier_thickness + } + + fn get_energy_barrier_thickness(&self) -> f64 { + self.energy_barrier_thickness + } + + fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + (0., y, z) + } +} /// Object that contains raw mesh input data. -#[derive(Deserialize)] +#[derive(Deserialize, Clone)] pub struct Mesh2DInput { pub length_unit: String, pub triangles: Vec<(f64, f64, f64, f64, f64, f64)>, @@ -16,39 +140,67 @@ pub struct Mesh2DInput { } /// Triangular mesh for rustbca. +#[derive(Clone)] pub struct Mesh2D { mesh: Vec, pub boundary: Polygon, pub simulation_boundary: Polygon, pub energy_barrier_thickness: f64 } + impl Mesh2D { - /// Constructor for Mesh2D object from mesh_2d_input. - pub fn new(mesh_2d_input: Mesh2DInput) -> Mesh2D { + /// Finds the cell that is nearest to (x, y). + fn nearest_to(&self, x: f64, y: f64, z: f64) -> &Cell2D { - let triangles = mesh_2d_input.triangles; - let material_boundary_points = mesh_2d_input.material_boundary_points; - let simulation_boundary_points = mesh_2d_input.simulation_boundary_points; - let electronic_stopping_correction_factors = mesh_2d_input.electronic_stopping_correction_factors; + let mut min_distance: f64 = std::f64::MAX; + let mut index: usize = 0; + + for (cell_index, cell) in self.mesh.iter().enumerate() { + let distance_to = cell.distance_to(x, y, z); + if distance_to < min_distance { + min_distance = distance_to; + index = cell_index; + } + } + + return &self.mesh[index]; + } +} + +impl GeometryInput for Mesh2D { + type GeometryInput = Mesh2DInput; +} + +impl Geometry for Mesh2D { + + type InputFileFormat = Input2D; + + /// Constructor for Mesh2D object from geometry_input. + fn new(geometry_input: &::GeometryInput) -> Mesh2D { + + let triangles = geometry_input.triangles.clone(); + let material_boundary_points = geometry_input.material_boundary_points.clone(); + let simulation_boundary_points = geometry_input.simulation_boundary_points.clone(); + let electronic_stopping_correction_factors = geometry_input.electronic_stopping_correction_factors.clone(); let n = triangles.len(); let mut cells: Vec = Vec::with_capacity(n); //Multiply all coordinates by value of geometry unit. - let length_unit: f64 = match mesh_2d_input.length_unit.as_str() { + let length_unit: f64 = match geometry_input.length_unit.as_str() { "MICRON" => MICRON, "CM" => CM, "ANGSTROM" => ANGSTROM, "NM" => NM, "M" => 1., - _ => mesh_2d_input.length_unit.parse() + _ => 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", &mesh_2d_input.length_unit.as_str() + "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 densities: Vec> = mesh_2d_input.densities + let densities: Vec> = geometry_input.densities .iter() .map( |row| row.iter().map(|element| element/(length_unit).powf(3.)).collect() ).collect(); @@ -89,7 +241,7 @@ impl Mesh2D { } let simulation_boundary: Polygon = Polygon::new(LineString::from(simulation_boundary_points_converted), vec![]); - let energy_barrier_thickness = mesh_2d_input.energy_barrier_thickness*length_unit; + let energy_barrier_thickness = geometry_input.energy_barrier_thickness*length_unit; Mesh2D { mesh: cells, @@ -99,10 +251,45 @@ impl Mesh2D { } } + fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { + if self.inside(x, y, z) { + true + } else { + let nearest_cell = self.nearest_to(x, y, z); + let distance = nearest_cell.distance_to(x, y, z); + distance < self.get_energy_barrier_thickness() + } + } + + fn get_ck_nearest_to(&self, x: f64, y: f64, z: f64) -> f64 { + self.nearest_to(x, y, z).get_electronic_stopping_correction_factor() + } + + fn get_densities_nearest_to(&self, x: f64, y: f64, z: f64) -> &Vec { + self.nearest_to(x, y, z).get_densities() + } + + fn get_energy_barrier_thickness(&self) -> f64 { + self.energy_barrier_thickness + } + + fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool { + self.simulation_boundary.contains(&point!(x: x, y: y)) + } + + 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)) { + (p.x(), p.y(), z) + } else { + panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z); + } + + } + /// Find the number densities of the triangle that contains or is nearest to (x, y). - pub fn get_densities(&self, x: f64, y: f64) -> &Vec { + fn get_densities(&self, x: f64, y: f64, z: f64) -> &Vec { for cell in &self.mesh { - if cell.contains(x, y) { + if cell.contains(x, y, z) { return &cell.densities; } } @@ -110,9 +297,9 @@ impl Mesh2D { } /// Find the number densities of the triangle that contains or is nearest to (x, y). - pub fn get_ck(&self, x: f64, y: f64) -> f64 { + fn get_ck(&self, x: f64, y: f64, z: f64) -> f64 { for cell in &self.mesh { - if cell.contains(x, y) { + if cell.contains(x, y, z) { return cell.electronic_stopping_correction_factor; } } @@ -120,9 +307,9 @@ impl Mesh2D { } /// Determine the total number density of the triangle that contains or is nearest to (x, y). - pub fn get_total_density(&self, x: f64, y: f64) -> f64 { + fn get_total_density(&self, x: f64, y: f64, z: f64) -> f64 { for cell in &self.mesh { - if cell.contains(x, y) { + if cell.contains(x, y, z) { return cell.densities.iter().sum::(); } } @@ -130,43 +317,27 @@ impl Mesh2D { } /// Find the concentrations of the triangle that contains or is nearest to (x, y). - pub fn get_concentrations(&self, x: f64, y: f64) -> &Vec { - if self.inside(x, y) { + fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec { + if self.inside(x, y, z) { for cell in &self.mesh { - if cell.contains(x, y) { + if cell.contains(x, y, z) { return &cell.concentrations; } } panic!("Geometry error: method inside() is returning true for points outside all cells. Check boundary points.") } else { - return &self.nearest_cell_to(x, y).concentrations; + return &self.nearest_to(x, y, z).concentrations; } } /// Determines whether the point (x, y) is inside the mesh. - pub fn inside(&self, x: f64, y: f64) -> bool { + fn inside(&self, x: f64, y: f64, z: f64) -> bool { self.boundary.contains(&Point::new(x, y)) } - - /// Finds the cell that is nearest to (x, y). - pub fn nearest_cell_to(&self, x: f64, y: f64) -> &Cell2D { - - let mut min_distance: f64 = std::f64::MAX; - let mut index: usize = 0; - - for (cell_index, cell) in self.mesh.iter().enumerate() { - let distance_to = cell.distance_to(x, y); - if distance_to < min_distance { - min_distance = distance_to; - index = cell_index; - } - } - - return &self.mesh[index]; - } } /// A mesh cell that contains a triangle and the local number densities and concentrations. +#[derive(Clone)] pub struct Cell2D { triangle: Triangle2D, pub densities: Vec, @@ -183,19 +354,34 @@ impl Cell2D { electronic_stopping_correction_factor } } +} +impl GeometryElement for Cell2D { /// Determines whether this cell contains the point (x, y). - pub fn contains(&self, x: f64, y: f64) -> bool { + fn contains(&self, x: f64, y: f64, z: f64) -> bool { self.triangle.contains(x, y) } /// Computes the shortest distance between this cell and the point (x, y). - pub fn distance_to(&self, x: f64, y: f64) -> f64 { + fn distance_to(&self, x: f64, y: f64, z: f64) -> f64 { self.triangle.distance_to(x, y) } + + fn get_densities(&self) -> &Vec { + &self.densities + } + + fn get_concentrations(&self) -> &Vec { + &self.concentrations + } + + fn get_electronic_stopping_correction_factor(&self) -> f64 { + self.electronic_stopping_correction_factor + } } /// A triangle in 2D, with points (x1, y1), (x2, y2), (x3, y3), and the three line segments bewtween them. +#[derive(Clone)] pub struct Triangle2D { x1: f64, x2: f64, diff --git a/src/particle.rs b/src/particle.rs index 95ef3f4..8387417 100644 --- a/src/particle.rs +++ b/src/particle.rs @@ -1,7 +1,7 @@ use super::*; /// Rustbca's internal representation of the particle_parameters input. -#[derive(Deserialize)] +#[derive(Deserialize, Clone)] pub struct ParticleParameters { pub particle_input_filename: String, pub length_unit: String, @@ -210,7 +210,6 @@ pub fn particle_advance(particle_1: &mut particle::Particle, mfp: f64, asymptoti /// Calcualte the refraction angle based on the surface binding energy of the material. pub fn refraction_angle(costheta: f64, energy_old: f64, energy_new: f64) -> f64 { - //println!("energy_old: {} energy_new: {} costheta: {}", energy_old/EV, energy_new/EV, costheta); let costheta = if costheta.abs() > 1. {costheta.signum()} else {costheta}; let sintheta0 = (1. - costheta*costheta).sqrt(); let sintheta1 = sintheta0*(energy_old/energy_new).sqrt(); diff --git a/src/tests.rs b/src/tests.rs index 53ce3b6..d6eaebc 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -3,6 +3,72 @@ use super::*; #[cfg(test)] use float_cmp::*; +#[test] +fn test_geometry() { + let mass = 1.; + let Z = 1.; + let E = 10.*EV; + let Ec = 1.*EV; + let Es = 5.76*EV; + let x = 0.; + let y = 0.; + let z = 0.; + let cosx = 1./(2.0_f64).sqrt(); + let cosy = 1./(2.0_f64).sqrt(); + let cosz = 0.; + let mut particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, 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], + Z: vec![29., 1.], + m: vec![63.54, 1.0008], + interaction_index: vec![0, 0], + surface_binding_model: SurfaceBindingModel::TARGET + }; + + let thickness: f64 = 1000.; + let depth: f64 = 1000.; + + let geometry_input_2D = mesh::Mesh2DInput { + length_unit: "ANGSTROM".to_string(), + triangles: vec![(0., depth, 0., thickness/2., thickness/2., -thickness/2.), (depth, depth, 0., thickness/2., -thickness/2., -thickness/2.)], + densities: vec![vec![0.03, 0.03], vec![0.03, 0.03]], + material_boundary_points: vec![(0., thickness/2.), (depth, thickness/2.), (depth, -thickness/2.), (0., -thickness/2.), (0., 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.), (0., 1.1*thickness/2.)], + energy_barrier_thickness: 10., + electronic_stopping_correction_factors: vec![1.0, 1.0], + }; + + let geometry_input_0D = mesh::Mesh0DInput { + length_unit: "ANGSTROM".to_string(), + densities: vec![0.03, 0.03], + electronic_stopping_correction_factor: 1.0 + }; + + 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; + + particle_1.pos.x = 500.*ANGSTROM; + particle_1.pos.y = 0.; + + particle_1.pos_old.x = -500.*ANGSTROM; + particle_1.pos_old.y = 0.; + + assert_eq!(material_2D.geometry.get_ck(0., 0., 0.), material_0D.geometry.get_ck(0., 0., 0.)); + assert_eq!(material_2D.geometry.get_energy_barrier_thickness(), material_0D.geometry.get_energy_barrier_thickness()); + 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.)); + assert_eq!(material_2D.geometry.get_densities_nearest_to(-10., 0., 5.), material_0D.geometry.get_densities_nearest_to(-10., 0., 5.)); + assert_eq!(material_2D.geometry.get_ck_nearest_to(-10., 0., 5.), material_0D.geometry.get_ck_nearest_to(-10., 0., 5.)); +} + #[test] fn test_surface_binding_energy_barrier() { let mass = 1.; @@ -32,17 +98,26 @@ fn test_surface_binding_energy_barrier() { let thickness: f64 = 1000.; let depth: f64 = 1000.; - let mesh_2d_input = mesh::Mesh2DInput { + + let geometry_input_2D = mesh::Mesh2DInput { length_unit: "ANGSTROM".to_string(), triangles: vec![(0., depth, 0., thickness/2., thickness/2., -thickness/2.), (depth, depth, 0., thickness/2., -thickness/2., -thickness/2.)], densities: vec![vec![0.03, 0.03], vec![0.03, 0.03]], material_boundary_points: vec![(0., thickness/2.), (depth, thickness/2.), (depth, -thickness/2.), (0., -thickness/2.), (0., 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.), (0., 1.1*thickness/2.)], energy_barrier_thickness: 10., - electronic_stopping_correction_factors: vec![0.0, 0.0], + electronic_stopping_correction_factors: vec![1.0, 1.0], + }; + + let geometry_input_0D = mesh::Mesh0DInput { + length_unit: "ANGSTROM".to_string(), + densities: vec![0.03, 0.03], + electronic_stopping_correction_factor: 1.0 }; - let material_1 = material::Material::new(material_parameters, mesh_2d_input); + 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; particle_1.pos.x = 500.*ANGSTROM; particle_1.pos.y = 0.; @@ -50,35 +125,65 @@ fn test_surface_binding_energy_barrier() { particle_1.pos_old.x = -500.*ANGSTROM; particle_1.pos_old.y = 0.; - let inside = material_1.inside(particle_1.pos.x, particle_1.pos.y); - let inside_old = material_1.inside(particle_1.pos_old.x, particle_1.pos_old.y); + let inside_0D = material_0D.inside(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z); + let inside_old_0D = material_0D.inside(particle_1.pos_old.x, particle_1.pos_old.y, particle_1.pos_old.z); + + let inside_2D = material_2D.inside(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z); + let inside_old_2D = material_2D.inside(particle_1.pos_old.x, particle_1.pos_old.y, particle_1.pos_old.z); //println!("{} {}", inside, inside_old); - assert!(inside); - assert!(!inside_old); + assert!(inside_0D); + assert!(!inside_old_0D); + + assert!(inside_2D); + assert!(!inside_old_2D); //Test concentration-dependent surface binding energy - let surface_binding_energy = material_1.actual_surface_binding_energy(&particle_1, particle_1.pos.x, particle_1.pos.y); - assert!(approx_eq!(f64, surface_binding_energy/EV, (2. + 4.)/2., epsilon=1E-12)); - //println!("sbv: {}", surface_binding_energy/EV); + let surface_binding_energy_2D = material_2D.actual_surface_binding_energy(&particle_1, particle_1.pos.x, particle_1.pos.y, particle_1.pos.z); + let surface_binding_energy_0D = material_0D.actual_surface_binding_energy(&particle_1, particle_1.pos.x, particle_1.pos.y, particle_1.pos.z); + + assert!(approx_eq!(f64, surface_binding_energy_0D/EV, (2. + 4.)/2., epsilon=1E-24)); + assert!(approx_eq!(f64, surface_binding_energy_2D/EV, (2. + 4.)/2., epsilon=1E-24)); + + println!("sbv 0D: {}", surface_binding_energy_0D/EV); + println!("sbv 2D: {}", surface_binding_energy_2D/EV); //Test leftmost boundary - assert!(material_1.inside_energy_barrier(500.*ANGSTROM, 0.)); - assert!(material_1.inside_energy_barrier(-5.*ANGSTROM, 0.)); - assert!(!material_1.inside_energy_barrier(-15.*ANGSTROM, 0.)); + assert!(material_2D.inside_energy_barrier(500.*ANGSTROM, 0., 0.)); + assert!(material_2D.inside_energy_barrier(-5.*ANGSTROM, 0., 0.)); + assert!(!material_2D.inside_energy_barrier(-15.*ANGSTROM, 0., 0.)); + + assert!(material_0D.inside_energy_barrier(500.*ANGSTROM, 0., 0.)); + assert!(material_0D.inside_energy_barrier(-5.*ANGSTROM, 0., 0.)); + assert!(!material_0D.inside_energy_barrier(-15.*ANGSTROM, 0., 0.)); + + assert!(material_0D.inside_energy_barrier(500.*ANGSTROM, 0., 0.)); + assert!(!material_0D.inside_energy_barrier(-500.*ANGSTROM, 0., 0.)); + assert!(material_0D.inside_energy_barrier(-9.*ANGSTROM, 0., 0.)); + assert!(!material_0D.inside_energy_barrier(-11.*ANGSTROM, 0., 0.)); //Test top boundary - assert!(material_1.inside_energy_barrier(500.*ANGSTROM, 0.)); - assert!(material_1.inside_energy_barrier(500.*ANGSTROM, 505.*ANGSTROM)); - assert!(!material_1.inside_energy_barrier(500.*ANGSTROM, 515.*ANGSTROM)); + assert!(material_2D.inside_energy_barrier(500.*ANGSTROM, 0., 0.)); + assert!(material_2D.inside_energy_barrier(500.*ANGSTROM, 505.*ANGSTROM, 0.)); + assert!(!material_2D.inside_energy_barrier(500.*ANGSTROM, 515.*ANGSTROM, 0.)); + + assert!(material_0D.inside_energy_barrier(500.*ANGSTROM, 0., 0.)); + assert!(material_0D.inside_energy_barrier(500.*ANGSTROM, 505.*ANGSTROM, 0.)); + assert!(material_0D.inside_energy_barrier(500.*ANGSTROM, 515.*ANGSTROM, 0.)); //Test bottom boundary - assert!(material_1.inside_energy_barrier(500.*ANGSTROM, -505.*ANGSTROM)); - assert!(!material_1.inside_energy_barrier(500.*ANGSTROM, -515.*ANGSTROM)); + assert!(material_2D.inside_energy_barrier(500.*ANGSTROM, -505.*ANGSTROM, 0.)); + assert!(!material_2D.inside_energy_barrier(500.*ANGSTROM, -515.*ANGSTROM, 0.)); + + assert!(material_0D.inside_energy_barrier(500.*ANGSTROM, -505.*ANGSTROM, 0.)); + assert!(material_0D.inside_energy_barrier(500.*ANGSTROM, -515.*ANGSTROM, 0.)); //Test rightmost boundary - assert!(material_1.inside_energy_barrier(1005.*ANGSTROM, 0.)); - assert!(!material_1.inside_energy_barrier(1015.*ANGSTROM, 0.)); + assert!(material_2D.inside_energy_barrier(1005.*ANGSTROM, 0., 0.)); + assert!(!material_2D.inside_energy_barrier(1015.*ANGSTROM, 0., 0.)); + + assert!(material_0D.inside_energy_barrier(1005.*ANGSTROM, 0., 0.)); + assert!(material_0D.inside_energy_barrier(1015.*ANGSTROM, 0., 0.)); } #[test] @@ -215,7 +320,7 @@ fn test_momentum_conservation() { let thickness: f64 = 1000.; let depth: f64 = 1000.; - let mesh_2d_input = mesh::Mesh2DInput { + let geometry_input = mesh::Mesh2DInput { length_unit: "ANGSTROM".to_string(), triangles: vec![(0., depth, 0., thickness/2., thickness/2., -thickness/2.), (depth, depth, 0., thickness/2., -thickness/2., -thickness/2.)], densities: vec![vec![0.06306], vec![0.06306]], @@ -225,7 +330,7 @@ fn test_momentum_conservation() { energy_barrier_thickness: 0., }; - let material_1 = material::Material::new(material_parameters, mesh_2d_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 in vec![InteractionPotential::KR_C, InteractionPotential::MOLIERE, InteractionPotential::ZBL, InteractionPotential::LENZ_JENSEN] { @@ -281,7 +386,7 @@ fn test_momentum_conservation() { 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.E = binary_collision_result.recoil_energy - material_1.average_bulk_binding_energy(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); //Rotate particle 1, 2 by lab frame scattering angles particle::rotate_particle(&mut particle_1, binary_collision_result.psi,