diff --git a/CMake/compile_time_options.cmake b/CMake/compile_time_options.cmake index 6c09df7b..2326f07c 100644 --- a/CMake/compile_time_options.cmake +++ b/CMake/compile_time_options.cmake @@ -6,7 +6,7 @@ set( CMAKE_BUILD_TYPE "Release" ) # preprocessor definitions in source code are defined below: set( description "(no description added yet)" ) -set( GITR_USE_CUDA 0 CACHE STRING "${description}" FORCE ) +set( GITR_USE_CUDA 1 CACHE STRING "${description}" FORCE ) set( GITR_USE_OPENMP 0 CACHE STRING "${description}" FORCE ) set( GITR_USE_MPI 0 CACHE STRING "${description}" FORCE ) set( GITR_USE_DOUBLE 1 CACHE STRING "${description}" FORCE ) diff --git a/CMake/crosslink_components.cmake b/CMake/crosslink_components.cmake index cb3b6cfe..cf15a068 100644 --- a/CMake/crosslink_components.cmake +++ b/CMake/crosslink_components.cmake @@ -76,6 +76,7 @@ target_link_libraries( GITR spectroscopy libconfig_cxx libconfig_c + particle_diagnostics utils boris surface_model @@ -158,15 +159,15 @@ target_link_libraries( coulomb_tests Catch2::Catch2WithMain ) -target_link_libraries( file_io_tests - libconfig_cxx - libconfig_c - flags - utils - boris - geometry_check - Catch2::Catch2WithMain - ) +#target_link_libraries( file_io_tests +# libconfig_cxx +# libconfig_c +# flags +# utils +# boris +# geometry_check +# Catch2::Catch2WithMain +# ) target_link_libraries( interpolator_tests Catch2::Catch2WithMain ) diff --git a/CMake/define_source_components.cmake b/CMake/define_source_components.cmake index 245fbaca..e09fcaa4 100644 --- a/CMake/define_source_components.cmake +++ b/CMake/define_source_components.cmake @@ -31,6 +31,7 @@ set( non_gpu_targets # conditionally compile as GPU targets set( gpu_targets + particle_diagnostics surface_model interp2d interpRateCoeff diff --git a/CMake/define_test_components.cmake b/CMake/define_test_components.cmake index 99280ed7..606345af 100644 --- a/CMake/define_test_components.cmake +++ b/CMake/define_test_components.cmake @@ -15,8 +15,8 @@ set( gpu_test_targets cross_field_diffusion_tests atomic_tests coulomb_tests - file_io_tests ) + #file_io_tests foreach( component IN LISTS cpu_test_targets ) diff --git a/CMake/enforce_build_order.cmake b/CMake/enforce_build_order.cmake index c3320896..aed13554 100644 --- a/CMake/enforce_build_order.cmake +++ b/CMake/enforce_build_order.cmake @@ -12,6 +12,7 @@ if( dependencies ) endif() +add_dependencies( boris particle_diagnostics ) # ensure that all source targets are built before GITR add_dependencies( GITR ${non_gpu_targets} ${gpu_targets} ) diff --git a/build.sh b/build.sh index 7d78441f..2d9c7bd0 100644 --- a/build.sh +++ b/build.sh @@ -1,12 +1,12 @@ -cmake -S $1 -B $2 -G Ninja \ --DNETCDF_C_HEADERS_DIR="/home/5n4/build_right/external/netcdf-c-install/include" \ --DNETCDF_CXX_HEADERS_DIR="/home/5n4/build_right/external/netcdf-cxx4-install/include" \ --DNETCDF_C_SHARED_LIB_DIR="/home/5n4/build_right/external/netcdf-c-install/lib" \ --DNETCDF_CXX_SHARED_LIB_DIR="/home/5n4/build_right/external/netcdf-cxx4-install/lib" \ --DLIBCONFIG_C_HEADERS_DIR="/home/5n4/build_right/external/libconfig_install/include" \ --DLIBCONFIG_CXX_HEADERS_DIR="/home/5n4/build_right/external/libconfig_install/include" \ --DLIBCONFIG_C_LIB_DIR="/home/5n4/build_right/external/libconfig_install/lib" \ --DLIBCONFIG_CXX_LIB_DIR="/home/5n4/build_right/external/libconfig_install/lib" \ +cmake -S $1 -B $2 \ +-DNETCDF_C_HEADERS_DIR="/opt/cray/pe/netcdf/4.9.0.9/nvidia/23.3/include" \ +-DNETCDF_CXX_HEADERS_DIR="/opt/cray/pe/netcdf/4.9.0.9/nvidia/23.3/include" \ +-DNETCDF_C_SHARED_LIB_DIR="/opt/cray/pe/netcdf/4.9.0.9/nvidia/23.3/lib" \ +-DNETCDF_CXX_SHARED_LIB_DIR="/opt/cray/pe/netcdf/4.9.0.9/nvidia/23.3/lib" \ +-DLIBCONFIG_C_HEADERS_DIR="/global/homes/t/tyounkin/code/libconfigBuild/gnu/include" \ +-DLIBCONFIG_CXX_HEADERS_DIR="/global/homes/t/tyounkin/code/libconfigBuild/gnu/include" \ +-DLIBCONFIG_C_LIB_DIR="/global/homes/t/tyounkin/code/libconfigBuild/gnu/lib" \ +-DLIBCONFIG_CXX_LIB_DIR="/global/homes/t/tyounkin/code/libconfigBuild/gnu/lib" \ &> $2/cmake_output.txt ; -cmake --build $2 -- -j 0 &> $2/build_output.txt +cmake --build $2 -- -j &> $2/build_output.txt diff --git a/build/make_gitr.sh b/build/make_gitr.sh new file mode 100755 index 00000000..bb60863a --- /dev/null +++ b/build/make_gitr.sh @@ -0,0 +1,29 @@ +#!/bin/bash +source ../env.default.sh + +cmake -DTHRUST_INCLUDE_DIR=$THRUST_DIR \ + -DNETCDF_CXX_INCLUDE_DIR=$NETCDFCXX4INCLUDE \ + -DNETCDF_CXX_LIBRARY=$NETCDFCXX4LIB \ + -DLIBCONFIG_INCLUDE_DIR=$LIBCONFIG_INCLUDE_DIR \ + -DLIBCONFIGPP_LIBRARY=$LIBCONFIGLIB \ + -DLIBCONFIGPP_INCLUDE_DIR=$LIBCONFIG_INCLUDE_DIR \ + -DNETCDF_DIR=$NETCDFDIR \ + -DNETCDF_INCLUDE_DIR=$NETCDFINCLUDE \ + -DNETCDF_LIBRARY=$NETCDFLIB \ + -DUSE_CUDA=1 \ + -DUSE_MPI=0 \ + -DUSE_DOUBLE=1 \ + .. +# -DNETCDF_CXX_INCLUDE_DIR=$NETCDFCXX4INCLUDE \ +# -DNETCDF_CXX_LIBRARY=$NETCDFLIB_CPP \ +# -DNETCDF_DIR=$NETCDFDIR \ +# -DNETCDF_INCLUDE_DIR=$NETCDFINCLUDE \ +# -DNETCDF_LIBRARY=$NETCDFLIB \ +# -DLIBCONFIGPP_LIBRARY=$LIBCONFIGLIB \ +# -DLIBCONFIGPP_INCLUDE_DIR=$LIBCONFIG_INCLUDE_DIR \ +# -DCPPFLAGS=-I${H5DIR}/include \ +# -DLDFLAGS=${H5DIR}/lib/libhdf5.dylib \ +# -DHDF5_LIBRARIES=${H5DIR}/lib/libhdf5.a \ +# -DHDF5_INCLUDE_DIRS=${H5DIR}/include \ +# -DHDF5_C_COMPILER_EXECUTABLE=${H5DIR}/bin/h5pcc \ +# -DHDF5_HL_LIBRARIES=${H5DIR}/lib/libhdf5_hl.dylib \ diff --git a/containers/run_podman_noninteractive.sh b/containers/run_podman_noninteractive.sh new file mode 100644 index 00000000..97f4b4bd --- /dev/null +++ b/containers/run_podman_noninteractive.sh @@ -0,0 +1,5 @@ +podman-hpc run -i -t --rm \ +--name gitr_container \ +--gpu \ +-v $(pwd):/host \ +docker.io/stonecoldhughes/gitr:noninteractive diff --git a/dependencies/build_all_dependencies.sh b/dependencies/build_all_dependencies.sh index bfee2a22..72a0f2e8 100755 --- a/dependencies/build_all_dependencies.sh +++ b/dependencies/build_all_dependencies.sh @@ -1,4 +1,4 @@ -export GITR_TOP_LEVEL=~/trash +export GITR_TOP_LEVEL=$PWD ./build_zlib.sh ./build_m4.sh ./build_hdf5.sh diff --git a/env.default.sh b/env.default.sh deleted file mode 100644 index ca5859b6..00000000 --- a/env.default.sh +++ /dev/null @@ -1,27 +0,0 @@ -export GITR_TOP_LEVEL=~/trash -#export TQD_HOME=/home/tqd -#export PATH=/home/tqd/Code/cmakebuild/bin:$PATH:$TQD_HOME/Code/pythonbuild/bin:$TQD_HOME/Code/openmpiBuild/bin:/home/tqd/Code/llvm-project/build/bin -#export CMAKE_SOURCE_DIR=../ -#export NETCDF=$TQD_HOME/Code/netcdfcbuild -export NETCDFCXX4=$GITR_TOP_LEVEL/netcdfcxxbuild -#export NETCDFDIR=$NETCDF/lib -#export NETCDFCXX4DIR=$NETCDFCXX4/lib -#export NETCDFINCLUDE=$NETCDF/include -export NETCDFCXX4INCLUDE=$NETCDFCXX4/include -export NETCDFCXX4LIB=$NETCDFCXX4/lib/libnetcdf_c++4.so -#export NETCDFLIB=$NETCDFDIR/libnetcdf.so -#export NETCDFLIB_CPP=$NETCDFCXX4DIR/libnetcdf_c++4.so -#export PATH=$PATH:/$NETCDF/bin -# -export LIBCONFIGDIR=$GITR_TOP_LEVEL/libconfigbuild -export LIBCONFIGLIB=$LIBCONFIGDIR/lib/libconfig++.so -#export LIBCONFIGPP_LIBRARIES=lconfig++ -#export LIBCONFIGPP_LIBRARY=lconfig++ -#export LIBCONFIGPP_STATIC_LIBRARY= -export LIBCONFIG_INCLUDE_DIR=$LIBCONFIGDIR/include -#export LIBCONFIGPP_INCLUDE_DIR=$LIBCONFIGDIR/include -#export LIBCONFIG_LIBRARY=lconfig -#export THRUST_INCLUDE_DIRS=$TQD_HOME/Code/thrust -#export THRUST_INCLUDE_DIR=$TQD_HOME/Code/thrust -export THRUST_DIR=$GITR_TOP_LEVEL/thrust -#export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$NETCDF/lib:$NETCDFCXX4/lib:$TQD_HOME/Code/openmpiBuild/lib diff --git a/examples/circle/input/gitrGeometry.cfg b/examples/circle/input/gitrGeometry.cfg new file mode 100644 index 00000000..72eb5a65 --- /dev/null +++ b/examples/circle/input/gitrGeometry.cfg @@ -0,0 +1,56 @@ +geom = +{ + x1 = + [ + -10000.0 + ]; + z1 = + [ + 0.0 + ]; + x2 = + [ + 10000.0 + ]; + z2 = + [ + 0.0 + ]; + slope = + [ + 0.0 + ]; + intercept = + [ + 0.0 + ]; + length = + [ + 20000.0 + ]; + Z = + [ + 1.0, + 1.0 + ]; + surface = + [ + 1 + ]; + inDir = + [ + -1 + ]; + potential = + [ + 25 + ]; + y1 = -10000.0; + y2 = 10000.0; + periodic = 0; + theta0 = 0.0; + theta1 = 0.0; + periodic_bc_x0 = -10000.0; + periodic_bc_x1 = 100000.0; + periodic_bc_x = 0; +}; diff --git a/examples/circle/input/gitrInput.cfg b/examples/circle/input/gitrInput.cfg new file mode 100755 index 00000000..24ba80ae --- /dev/null +++ b/examples/circle/input/gitrInput.cfg @@ -0,0 +1,345 @@ +geometry = +{ + fileString = "gitrGeometry.cfg"; +}; +particleSource = +{ + fileString = "particleSource.cfg"; + ncFileString = "particleSource.nc"; +}; +backgroundPlasmaProfiles = +{ + Z = 1; + amu = 2; + biasPotential = 0.0; + sheath_factor = 1.0; + Bfield = + { + r = 1.9997; + z = 0.0349; + y = 0.0; + fileString = "bFromGfile.nc"; + gridNrString = "nR"; + gridNyString = "nY"; + gridNzString = "nZ"; + gridRString = "r"; + gridYString = "y"; + gridZString = "z"; + rString = "br"; + zString = "bz"; + yString = "bt"; + }; + Efield = + { + Er = 0.0; + Ez = 0.0; + Et = 0.0; + fileString = "LcS.nc"; + gridNrString = "nR"; + gridNyString = "nY"; + gridNzString = "nZ"; + gridRString = "gridR"; + gridYString = "gridY"; + gridZString = "gridZ"; + radialComponentString = "PSEr"; + axialComponentString = "PSEz"; + toroidalComponentString = "PSEt"; + }; + Temperature = + { + ti = 18.0; + te = 30.0; + fileString = "profiles.nc"; + gridNrString = "nR"; + gridNzString = "nZ"; + gridRString = "gridR"; + gridZString = "gridZ"; + IonTempString = "ti"; + ElectronTempString = "te"; + }; + Density = + { + ni = 1.0e+19; + ne = 1.0e+19; + fileString = "profiles.nc"; + gridNrString = "nR"; + gridNzString = "nZ"; + gridRString = "gridR"; + gridZString = "gridZ"; + IonDensityString = "ne"; + ElectronDensityString = "ne"; + }; + Diffusion = + { + Dperp = 0.1; + fileString = "profiles.nc"; + gridNrString = "n_x"; + gridNzString = "n_z"; + gridRString = "gridx"; + gridZString = "gridz"; + variableString = "ni"; + }; + FlowVelocity = + { + interpolatorNumber = 0; + flowVr = -19997.0; + flowVy = 0.0; + flowVz = -349.0; + fileString = "LcS.nc"; + gridNrString = "nR"; + gridNyString = "nY"; + gridNzString = "nZ"; + gridRString = "gridR"; + gridYString = "gridY"; + gridZString = "gridZ"; + flowVrString = "flowVr"; + flowVzString = "flowVz"; + flowVtString = "flowVt"; + }; + gradT = + { + gradTeR = -3.0; + gradTeY = 0.0; + gradTeZ = 0.0; + gradTiR = -30.0; + gradTiY = 0.0; + gradTiZ = 0.0; + fileString = "profiles.nc"; + gridNrString = "nX_gradTi"; + gridNzString = "nZ_gradTi"; + gridRString = "gridx_gradTi"; + gridZString = "gridz_gradTi"; + gradTiRString = "gradTiR"; + gradTiZString = "gradTiZ"; + gradTeRString = "gradTeR"; + gradTeZString = "gradTeZ"; + }; +}; +impurityParticleSource = +{ + nP = 1000; + sourceStrength = 1e+19; + Z = 74.0; + source_material_Z = 74; + source_material_SurfaceBindingEnergy = 11.75; + initialConditions = + { + x_start = 1.4; + y_start = 0.0; + z_start = -0.8; + impurity_amu = 184; + impurity_Z = 74.0; + charge = 0; + energy_eV = 10.0; + phi = 90.0; + theta = 90.0; + }; + ionization = + { + fileString = "ADAS_Rates_W.nc"; + TempGridString = "n_Temperatures_Ionize"; + DensGridString = "n_Densities_Ionize"; + nChargeStateString = "n_ChargeStates_Ionize"; + TempGridVarName = "gridTemperature_Ionization"; + DensGridVarName = "gridDensity_Ionization"; + CoeffVarName = "IonizationRateCoeff"; + }; + recombination = + { + fileString = "ADAS_Rates_W.nc"; + TempGridString = "n_Temperatures_Recombine"; + DensGridString = "n_Densities_Recombine"; + nChargeStateString = "n_ChargeStates_Recombine"; + TempGridVarName = "gridTemperature_Recombination"; + DensGridVarName = "gridDensity_Recombination"; + CoeffVarName = "RecombinationRateCoeff"; + }; +}; +timeStep = +{ + dt = 1e-09; + nPtsPerGyroOrbit = 10000.0; + ionization_nDtPerApply = 1; + collision_nDtPerApply = 5; + nT = 5000; +}; +surfaces = +{ + useMaterialSurfaces = 1; + flux = + { + nE = 10; + E0 = 0.0; + E = 1000.0; + nA = 9; + A0 = 0.0; + A = 90.0; + }; +}; +diagnostics = +{ + trackSubSampleFactor = 1; + netx0 = 0.95; + netx1 = 2.95; + nX = 150; + nety0 = -0.03; + leakZ = 0.0; + nety1 = 0.03; + nY = 120; + netz0 = -1.5; + netz1 = 1.5; + nZ = 300; + densityChargeBins = 6; + +}; + +particle_diagnostics = +{ + times_logarithmic = true; + bin_edge_0_time = -9; + bin_edge_1_time = -5; # if bin_edge_0_time not defined, defaults to log10(dt), and bin_edge_1_time is maximum edge + n_bins_time = 100; + bin_edge_0_angle = 1.5708 + bin_edge_1_angle = 62.8319 + n_bins_angle = 100 +}; + +forceEvaluation = +{ + X0 = -0.03; + X1 = 0.03; + Y0 = -0.03; + Y1 = 0.03; + Z0 = -0.015; + Z1 = 0.03; + nR = 176; + nY = 0; + nZ = 372; + particleEnergy = 10.0; +}; +surfaceModel = +{ + fileString = "simpleSurfaceModel8ev.nc"; + nEsputtRefCoeffString = "nE"; + nAsputtRefCoeffString = "nA"; + nEsputtRefDistInString = "nE"; + nAsputtRefDistInString = "nA"; + nEsputtRefDistOutString = "nEdistBins"; + nEsputtRefDistOutStringRef = "nEdistBinsRef"; + nAsputtRefDistOutString = "nAdistBins"; + E_sputtRefCoeff = "E"; + A_sputtRefCoeff = "A"; + E_sputtRefDistIn = "E"; + A_sputtRefDistIn = "A"; + E_sputtRefDistOut = "eDistEgrid"; + E_sputtRefDistOutRef = "eDistEgridRef"; + Aphi_sputtRefDistOut = "phiGrid"; + Atheta_sputtRefDistOut = "thetaGrid"; + sputtYldString = "spyld"; + reflYldString = "rfyld"; + EDist_Y = "energyDist"; + AphiDist_Y = "cosXDist"; + AthetaDist_Y = "cosYDist"; + EDist_R = "energyDistRef"; + AphiDist_R = "cosXDistRef"; + AthetaDist_R = "cosYDistRef"; +}; +geometry_hash = +{ + nHashes = 1; + hashX0 = -0.06; + hashX1 = 0.06; + hashY0 = -0.06; + hashY1 = 0.06; + hashZ0 = -0.01; + hashZ1 = 0.11; + nR_closeGeom = 1000; + nY_closeGeom = 960; + nZ_closeGeom = 40; + n_closeGeomElements = 10; + fileString = "geomHash0.nc"; + gridNrString = "nR"; + gridNyString = "nY"; + gridNzString = "nZ"; + nearestNelementsString = "n"; + gridRString = "gridR"; + gridYString = "gridY"; + gridZString = "gridZ"; + closeGeomString = "hash"; +}; +geometry_sheath = +{ + nHashes = 1; + hashX0 = -0.06; + hashX1 = 0.06; + hashY0 = -0.06; + hashY1 = 0.06; + hashZ0 = -0.01; + hashZ1 = 0.11; + nR_closeGeom = 1000; + nY_closeGeom = 960; + nZ_closeGeom = 40; + n_closeGeomElements = 10; + fileString = "geomHash_sheath.nc"; + gridNrString = "nR"; + gridNyString = "nY"; + gridNzString = "nZ"; + nearestNelementsString = "n"; + gridRString = "gridR"; + gridYString = "gridY"; + gridZString = "gridZ"; + closeGeomString = "hash"; +}; +flags = +{ + USE_CUDA = 1; + USE_MPI = 0; + USE_OPENMP = 0; + USE_IONIZATION = 1; + USEIONIZATION = 1; + USERECOMBINATION = 1; + USEPERPDIFFUSION = 0; + USEPARDIFFUSION = 0; + USECOULOMBCOLLISIONS = 1; + USEFRICTION = 1; + USEANGLESCATTERING = 1; + USEHEATING = 1; + USETHERMALFORCE = 0; + USESURFACEMODEL = 0; + USESHEATHEFIELD = 1; + BIASED_SURFACE = 0; + USEPRESHEATHEFIELD = 0; + BFIELD_INTERP = 0; + LC_INTERP = 0; + GENERATE_LC = 0; + EFIELD_INTERP = 0; + PRESHEATH_INTERP = 0; + DENSITY_INTERP = 0; + TEMP_INTERP = 0; + FLOWV_INTERP = 0; + GRADT_INTERP = 0; + ODEINT = 0; + FIXED_SEEDS = 1; + FIXEDSEEDS = 1; + PARTICLESEEDS = 1; + GEOM_TRACE = 0; + GEOM_HASH = 0; + GEOM_HASH_SHEATH = 0; + PARTICLE_TRACKS = 1; + PARTICLE_SOURCE_SPACE = 0; + PARTICLE_SOURCE_ENERGY = 0; + PARTICLE_SOURCE_ANGLE = 0; + PARTICLE_SOURCE_FILE = 1; + SPECTROSCOPY = 3; + USE3DTETGEOM = 0; + USECYLSYMM = 0; + USEFIELDALIGNEDVALUES = 0; + FLUX_EA = 1; + FORCE_EVAL = 0; + USE_SORT = 0; + USE_ADAPTIVE_DT = 0; + CHECK_COMPATIBILITY = 1; + USE_SURFACE_POTENTIAL=0; + USE_SHEATH_DENSITY=0; + USE_PARTICLE_DIAGNOSTICS=1; +}; diff --git a/examples/line_segments_graded_2d/input/gitrInput.cfg b/examples/line_segments_graded_2d/input/gitrInput.cfg new file mode 100755 index 00000000..a56d0a12 --- /dev/null +++ b/examples/line_segments_graded_2d/input/gitrInput.cfg @@ -0,0 +1,342 @@ +geometry = +{ + fileString = "gitrGeometry.cfg"; +}; +particleSource = +{ + fileString = "particleSource.cfg"; + ncFileString = "particleSource.nc"; +}; +backgroundPlasmaProfiles = +{ + Z = 1; + amu = 2; + biasPotential = 0.0; + sheath_factor = 1.0; + Bfield = + { + r = 1.9997; + z = 0.0349; + y = 0.0; + fileString = "bFromGfile.nc"; + gridNrString = "nR"; + gridNyString = "nY"; + gridNzString = "nZ"; + gridRString = "r"; + gridYString = "y"; + gridZString = "z"; + rString = "br"; + zString = "bz"; + yString = "bt"; + }; + Efield = + { + Er = 0.0; + Ez = 0.0; + Et = 0.0; + fileString = "LcS.nc"; + gridNrString = "nR"; + gridNyString = "nY"; + gridNzString = "nZ"; + gridRString = "gridR"; + gridYString = "gridY"; + gridZString = "gridZ"; + radialComponentString = "PSEr"; + axialComponentString = "PSEz"; + toroidalComponentString = "PSEt"; + }; + Temperature = + { + ti = 18.0; + te = 30.0; + fileString = "profiles.nc"; + gridNrString = "nR"; + gridNzString = "nZ"; + gridRString = "gridR"; + gridZString = "gridZ"; + IonTempString = "ti"; + ElectronTempString = "te"; + }; + Density = + { + ni = 1.0e+19; + ne = 1.0e+19; + fileString = "profiles.nc"; + gridNrString = "nR"; + gridNzString = "nZ"; + gridRString = "gridR"; + gridZString = "gridZ"; + IonDensityString = "ne"; + ElectronDensityString = "ne"; + }; + Diffusion = + { + Dperp = 0.1; + fileString = "profiles.nc"; + gridNrString = "n_x"; + gridNzString = "n_z"; + gridRString = "gridx"; + gridZString = "gridz"; + variableString = "ni"; + }; + FlowVelocity = + { + interpolatorNumber = 0; + flowVr = -19997.0; + flowVy = 0.0; + flowVz = -349.0; + fileString = "LcS.nc"; + gridNrString = "nR"; + gridNyString = "nY"; + gridNzString = "nZ"; + gridRString = "gridR"; + gridYString = "gridY"; + gridZString = "gridZ"; + flowVrString = "flowVr"; + flowVzString = "flowVz"; + flowVtString = "flowVt"; + }; + gradT = + { + gradTeR = -3.0; + gradTeY = 0.0; + gradTeZ = 0.0; + gradTiR = -30.0; + gradTiY = 0.0; + gradTiZ = 0.0; + fileString = "profiles.nc"; + gridNrString = "nX_gradTi"; + gridNzString = "nZ_gradTi"; + gridRString = "gridx_gradTi"; + gridZString = "gridz_gradTi"; + gradTiRString = "gradTiR"; + gradTiZString = "gradTiZ"; + gradTeRString = "gradTeR"; + gradTeZString = "gradTeZ"; + }; +}; +impurityParticleSource = +{ + nP = 100000; + sourceStrength = 1e+19; + Z = 74.0; + source_material_Z = 74; + source_material_SurfaceBindingEnergy = 11.75; + initialConditions = + { + x_start = 1.4; + y_start = 0.0; + z_start = -0.8; + impurity_amu = 184; + impurity_Z = 74.0; + charge = 0; + energy_eV = 10.0; + phi = 90.0; + theta = 90.0; + }; + ionization = + { + fileString = "ADAS_Rates_W.nc"; + TempGridString = "n_Temperatures_Ionize"; + DensGridString = "n_Densities_Ionize"; + nChargeStateString = "n_ChargeStates_Ionize"; + TempGridVarName = "gridTemperature_Ionization"; + DensGridVarName = "gridDensity_Ionization"; + CoeffVarName = "IonizationRateCoeff"; + }; + recombination = + { + fileString = "ADAS_Rates_W.nc"; + TempGridString = "n_Temperatures_Recombine"; + DensGridString = "n_Densities_Recombine"; + nChargeStateString = "n_ChargeStates_Recombine"; + TempGridVarName = "gridTemperature_Recombination"; + DensGridVarName = "gridDensity_Recombination"; + CoeffVarName = "RecombinationRateCoeff"; + }; +}; +timeStep = +{ + dt = 1e-09; + nPtsPerGyroOrbit = 10000.0; + ionization_nDtPerApply = 1; + collision_nDtPerApply = 5; + nT = 5000; +}; +surfaces = +{ + useMaterialSurfaces = 1; + flux = + { + nE = 10; + E0 = 0.0; + E = 1000.0; + nA = 9; + A0 = 0.0; + A = 90.0; + }; +}; +diagnostics = +{ + trackSubSampleFactor = 1; + netx0 = 0.95; + netx1 = 2.95; + nX = 150; + nety0 = -0.03; + leakZ = 0.0; + nety1 = 0.03; + nY = 120; + netz0 = -1.5; + netz1 = 1.5; + nZ = 300; + densityChargeBins = 6; + +}; + +particle_diagnostics = +{ + times_logarithmic = true; + bin_edge_0_time = -9; + bin_edge_1_time = -5; # if bin_edge_0_time not defined, defaults to log10(dt), and bin_edge_1_time is maximum edge + n_bins_time = 100; +}; + +forceEvaluation = +{ + X0 = -0.03; + X1 = 0.03; + Y0 = -0.03; + Y1 = 0.03; + Z0 = -0.015; + Z1 = 0.03; + nR = 176; + nY = 0; + nZ = 372; + particleEnergy = 10.0; +}; +surfaceModel = +{ + fileString = "simpleSurfaceModel8ev.nc"; + nEsputtRefCoeffString = "nE"; + nAsputtRefCoeffString = "nA"; + nEsputtRefDistInString = "nE"; + nAsputtRefDistInString = "nA"; + nEsputtRefDistOutString = "nEdistBins"; + nEsputtRefDistOutStringRef = "nEdistBinsRef"; + nAsputtRefDistOutString = "nAdistBins"; + E_sputtRefCoeff = "E"; + A_sputtRefCoeff = "A"; + E_sputtRefDistIn = "E"; + A_sputtRefDistIn = "A"; + E_sputtRefDistOut = "eDistEgrid"; + E_sputtRefDistOutRef = "eDistEgridRef"; + Aphi_sputtRefDistOut = "phiGrid"; + Atheta_sputtRefDistOut = "thetaGrid"; + sputtYldString = "spyld"; + reflYldString = "rfyld"; + EDist_Y = "energyDist"; + AphiDist_Y = "cosXDist"; + AthetaDist_Y = "cosYDist"; + EDist_R = "energyDistRef"; + AphiDist_R = "cosXDistRef"; + AthetaDist_R = "cosYDistRef"; +}; +geometry_hash = +{ + nHashes = 1; + hashX0 = -0.06; + hashX1 = 0.06; + hashY0 = -0.06; + hashY1 = 0.06; + hashZ0 = -0.01; + hashZ1 = 0.11; + nR_closeGeom = 1000; + nY_closeGeom = 960; + nZ_closeGeom = 40; + n_closeGeomElements = 10; + fileString = "geomHash0.nc"; + gridNrString = "nR"; + gridNyString = "nY"; + gridNzString = "nZ"; + nearestNelementsString = "n"; + gridRString = "gridR"; + gridYString = "gridY"; + gridZString = "gridZ"; + closeGeomString = "hash"; +}; +geometry_sheath = +{ + nHashes = 1; + hashX0 = -0.06; + hashX1 = 0.06; + hashY0 = -0.06; + hashY1 = 0.06; + hashZ0 = -0.01; + hashZ1 = 0.11; + nR_closeGeom = 1000; + nY_closeGeom = 960; + nZ_closeGeom = 40; + n_closeGeomElements = 10; + fileString = "geomHash_sheath.nc"; + gridNrString = "nR"; + gridNyString = "nY"; + gridNzString = "nZ"; + nearestNelementsString = "n"; + gridRString = "gridR"; + gridYString = "gridY"; + gridZString = "gridZ"; + closeGeomString = "hash"; +}; +flags = +{ + USE_CUDA = 1; + USE_MPI = 0; + USE_OPENMP = 0; + USE_IONIZATION = 1; + USEIONIZATION = 1; + USERECOMBINATION = 1; + USEPERPDIFFUSION = 0; + USEPARDIFFUSION = 0; + USECOULOMBCOLLISIONS = 1; + USEFRICTION = 1; + USEANGLESCATTERING = 1; + USEHEATING = 1; + USETHERMALFORCE = 0; + USESURFACEMODEL = 0; + USESHEATHEFIELD = 1; + BIASED_SURFACE = 0; + USEPRESHEATHEFIELD = 0; + BFIELD_INTERP = 0; + LC_INTERP = 0; + GENERATE_LC = 0; + EFIELD_INTERP = 0; + PRESHEATH_INTERP = 0; + DENSITY_INTERP = 0; + TEMP_INTERP = 0; + FLOWV_INTERP = 0; + GRADT_INTERP = 0; + ODEINT = 0; + FIXED_SEEDS = 1; + FIXEDSEEDS = 1; + PARTICLESEEDS = 1; + GEOM_TRACE = 0; + GEOM_HASH = 0; + GEOM_HASH_SHEATH = 0; + PARTICLE_TRACKS = 0; + PARTICLE_SOURCE_SPACE = 0; + PARTICLE_SOURCE_ENERGY = 0; + PARTICLE_SOURCE_ANGLE = 0; + PARTICLE_SOURCE_FILE = 1; + SPECTROSCOPY = 3; + USE3DTETGEOM = 0; + USECYLSYMM = 0; + USEFIELDALIGNEDVALUES = 0; + FLUX_EA = 1; + FORCE_EVAL = 0; + USE_SORT = 0; + USE_ADAPTIVE_DT = 0; + CHECK_COMPATIBILITY = 1; + USE_SURFACE_POTENTIAL=0; + USE_SHEATH_DENSITY=0; + USE_PARTICLE_DIAGNOSTICS=1; +}; diff --git a/include/Particles.h b/include/Particles.h index c517a092..e4f8acb1 100644 --- a/include/Particles.h +++ b/include/Particles.h @@ -66,6 +66,7 @@ class Particles : public ManagedAllocation { sim::Array surfaceHit; sim::Array firstCollision; sim::Array transitTime; + sim::Array transitAngle; sim::Array distTraveled; sim::Array wallIndex; sim::Array perpDistanceToSurface; @@ -176,6 +177,7 @@ class Particles : public ManagedAllocation { int wHT = this->surfaceHit[indx]; int fcT = this->firstCollision[indx]; gitr_precision ttT = this->transitTime[indx]; + gitr_precision ttA = this->transitAngle[indx]; gitr_precision dtT = this->distTraveled[indx]; gitr_precision firstIonizationZT = this->firstIonizationZ[indx]; gitr_precision firstIonizationTT = this->firstIonizationT[indx]; @@ -212,6 +214,7 @@ class Particles : public ManagedAllocation { this->surfaceHit[indx] = this->surfaceHit[n]; this->firstCollision[indx] = this->firstCollision[n]; this->transitTime[indx] = this->transitTime[n]; + this->transitAngle[indx] = this->transitAngle[n]; this->distTraveled[indx] = this->distTraveled[n]; this->firstIonizationZ[indx] = this->firstIonizationZ[n]; this->firstIonizationT[indx] = this->firstIonizationT[n]; @@ -246,6 +249,7 @@ class Particles : public ManagedAllocation { this->surfaceHit[n] = wHT; this->firstCollision[n] = fcT; this->transitTime[n] = ttT; + this->transitAngle[n] = ttA; this->distTraveled[n] = dtT; this->firstIonizationZ[n] = firstIonizationZT; this->firstIonizationT[n] = firstIonizationTT; @@ -303,6 +307,7 @@ class Particles : public ManagedAllocation { surfaceHit{nParticles, -1}, firstCollision{nParticles, 1}, transitTime{nParticles, 0.0}, + transitAngle{nParticles, 0.0}, distTraveled{nParticles, 0.0}, wallIndex{nParticles,0}, perpDistanceToSurface{nParticles,0.0}, diff --git a/include/flags.hpp b/include/flags.hpp index f309b5a5..5a046381 100644 --- a/include/flags.hpp +++ b/include/flags.hpp @@ -17,12 +17,14 @@ class Flags : public ManagedAllocation const bool USE_IONIZATION; const bool FIXED_SEEDS; const bool USE_ADAPTIVE_DT; + const bool USE_PARTICLE_DIAGNOSTICS; const bool USE_SHEATH_DENSITY; CUDA_CALLABLE_MEMBER Flags(libconfig::Config &cfg) : USE_IONIZATION{initialize_flags(cfg,"USE_IONIZATION",1)}, FIXED_SEEDS{initialize_flags(cfg,"FIXED_SEEDS",1)}, USE_ADAPTIVE_DT{initialize_flags(cfg,"USE_ADAPTIVE_DT",0)}, + USE_PARTICLE_DIAGNOSTICS{initialize_flags(cfg,"USE_PARTICLE_DIAGNOSTICS",0)}, USE_SHEATH_DENSITY{initialize_flags(cfg,"USE_SHEATH_DENSITY",0)} {}; bool initialize_flags(libconfig::Config &cfg, std::string s, int default_value); }; diff --git a/include/history.h b/include/history.h index c75daf02..89d5815d 100644 --- a/include/history.h +++ b/include/history.h @@ -31,11 +31,12 @@ struct history { gitr_precision *histvz; gitr_precision *histcharge; gitr_precision *histweight; + gitr_precision *histperpDistanceToSurface; - history(Particles *_particlesPointer, int* _tt, int _nT,int _subSampleFac, int _nP, gitr_precision *_histX,gitr_precision *_histY,gitr_precision *_histZ, - gitr_precision *_histv,gitr_precision *_histvx,gitr_precision *_histvy,gitr_precision *_histvz, gitr_precision * _histcharge, gitr_precision * _histweight) : - particlesPointer(_particlesPointer), tt(_tt),nT(_nT),subSampleFac(_subSampleFac), nP(_nP), - histX(_histX),histY(_histY),histZ(_histZ),histv(_histv),histvx(_histvx),histvy(_histvy),histvz(_histvz), histcharge(_histcharge), histweight(_histweight) {} + history(Particles *_particlesPointer, int* _tt, int _nT,int _subSampleFac, int _nP, gitr_precision *_histX,gitr_precision *_histY,gitr_precision *_histZ, gitr_precision *_histv, + gitr_precision *_histvx,gitr_precision *_histvy,gitr_precision *_histvz, gitr_precision * _histcharge, gitr_precision * _histweight, gitr_precision * _histperpDistanceToSurface) : + particlesPointer(_particlesPointer), tt(_tt),nT(_nT),subSampleFac(_subSampleFac), nP(_nP), histX(_histX),histY(_histY),histZ(_histZ), + histv(_histv),histvx(_histvx),histvy(_histvy),histvz(_histvz), histcharge(_histcharge), histweight(_histweight), histperpDistanceToSurface(_histperpDistanceToSurface) {} CUDA_CALLABLE_MEMBER_DEVICE void operator()(std::size_t indx) const @@ -60,6 +61,7 @@ void operator()(std::size_t indx) const histvz[histInd] = particlesPointer->vz[indexP]; histcharge[histInd] = particlesPointer->charge[indexP]; histweight[histInd] = particlesPointer->weight[indexP]; + histperpDistanceToSurface[histInd] = particlesPointer->perpDistanceToSurface[indexP]; } //else //{ diff --git a/include/particle_diagnostics.h b/include/particle_diagnostics.h new file mode 100644 index 00000000..9704ba98 --- /dev/null +++ b/include/particle_diagnostics.h @@ -0,0 +1,57 @@ +#pragma once + +#ifdef __CUDACC__ +#define CUDA_CALLABLE_MEMBER_DEVICE __device__ +#define CUDA_CALLABLE_MEMBER_HOST __host__ +#else +#define CUDA_CALLABLE_MEMBER_DEVICE +#define CUDA_CALLABLE_MEMBER_HOST +#endif + +#include "Particles.h" +#include "Fields.h" +#include "flags.hpp" +#ifdef __CUDACC__ +#include +#include +#endif + +#ifdef __GNUC__ +#include +#include +#endif + +#if USE_DOUBLE +typedef double gitr_precision; +#else +typedef float gitr_precision; +#endif + + +struct particle_diagnostics { + Flags *flags; + Particles *particlesPointer; + Boundary *boundaryVector; + bool times_logarithmic; + gitr_precision bin_edge_0_time; + gitr_precision bin_edge_1_time; + gitr_precision bin_edge_dt; + int n_bins_time; + gitr_precision *particle_time_histogram; + bool angle_logarithmic; + gitr_precision bin_edge_0_angle; + gitr_precision bin_edge_1_angle; + gitr_precision bin_edge_dtheta; + int n_bins_angle; + gitr_precision *particle_angle_histogram; + int nSurfaces; + + particle_diagnostics(Flags *_flags, Particles *_particlesPointer, Boundary *_boundaryVector, bool _time_logarithmic, + gitr_precision _bin_edge_0_time, gitr_precision _bin_edge_1_time, gitr_precision _bin_edge_dt, + int _n_bins_time, gitr_precision *_particle_time_histogram, bool _angle_logarithmic, + gitr_precision _bin_edge_0_angle, gitr_precision _bin_edge_1_angle, gitr_precision _bin_edge_dtheta, + int _n_bins_angle, gitr_precision *_particle_angle_histogram, int _nSurfaces); + + CUDA_CALLABLE_MEMBER_DEVICE + void operator()(std::size_t indx); +}; diff --git a/src/boris.cpp b/src/boris.cpp index a6d9d332..6621a04b 100644 --- a/src/boris.cpp +++ b/src/boris.cpp @@ -838,6 +838,7 @@ void move_boris::operator()(std::size_t indx) gitr_precision q_prime = 0.0; gitr_precision coeff = 0.0; gitr_precision minDist = 0.0; + gitr_precision perp_dist = 0.0; int closestBoundaryIndex; #if ODEINT == 0 @@ -1086,7 +1087,8 @@ void move_boris::operator()(std::size_t indx) particlesPointer->z[indx] = position[2]; particlesPointer->vx[indx] = v[0]; particlesPointer->vy[indx] = v[1]; - particlesPointer->vz[indx] = v[2]; + particlesPointer->vz[indx] = v[2]; + particlesPointer->perpDistanceToSurface[indx] = minDist; particlesPointer->time[indx] = particlesPointer->time[indx]+dt; particlesPointer->dt[indx] = new_dt; particlesPointer->advance[indx] = new_advance; @@ -1107,8 +1109,10 @@ void move_boris::operator()(std::size_t indx) particlesPointer->z[indx] = position[2] + v[2] * dt; particlesPointer->vx[indx] = v[0]; particlesPointer->vy[indx] = v[1]; - particlesPointer->vz[indx] = v[2]; + particlesPointer->vz[indx] = v[2]; + particlesPointer->perpDistanceToSurface[indx] = minDist; particlesPointer->angle[indx] = particlesPointer->angle[indx] + radians_per_dt; + particlesPointer->time[indx] = particlesPointer->time[indx]+dt; if(particlesPointer->z[indx] > particlesPointer->max_z[indx]) { particlesPointer->max_z[indx] = particlesPointer->z[indx]; @@ -1144,11 +1148,11 @@ void move_boris::operator()(std::size_t indx) gitr_precision half = 0.5; v[0] = particlesPointer->vx[indx]; v[1] = particlesPointer->vy[indx]; - v[2] = particlesPointer->vz[indx]; + v[2] = particlesPointer->vz[indx]; r[0] = particlesPointer->xprevious[indx]; r[1] = particlesPointer->yprevious[indx]; - r[2] = particlesPointer->zprevious[indx]; + r[2] = particlesPointer->zprevious[indx]; #ifdef __CUDACC__ #else #endif diff --git a/src/gitr.cpp b/src/gitr.cpp index 2fa4c089..fbbf9754 100644 --- a/src/gitr.cpp +++ b/src/gitr.cpp @@ -6,6 +6,7 @@ #include "array.h" #include "boris.h" #include "boundaryInit.h" +#include "constants.h" #include "coulombCollisions.h" #include "crossFieldDiffusion.h" #include "curandInitialize.h" @@ -23,6 +24,7 @@ #include "recombine.h" #include "spectroscopy.h" #include "surfaceModel.h" +#include "particle_diagnostics.h" //#include "testRoutineCuda.h" #include "thermalForce.h" #include "utils.h" @@ -3882,6 +3884,102 @@ if( efield_interp == 1 ) } */ + ////////////////////////////////////////////////// + ///// PARTICLE DIAGNOSTICS /////////// + ///////////////////////////////////////////////// + // This section will check the input file for the required fields for particle + // diagnostics, then read and set up the required vectors to collect the data + + + // user-set variables + bool times_logarithmic; + gitr_precision bin_edge_0_time; + gitr_precision bin_edge_1_time; + gitr_precision bin_edge_dt; + int n_bins_time; + bool angle_logarithmic; + gitr_precision bin_edge_0_angle; + gitr_precision bin_edge_1_angle; + gitr_precision bin_edge_dtheta; + int n_bins_angle; + + if (gitr_flags->USE_PARTICLE_DIAGNOSTICS) + { + std::cout << "Starting particle diagnostics initialization" << std::endl; + + // Set remainder of times parameters + if (cfg.lookupValue("particle_diagnostics.times_logarithmic",times_logarithmic)) {} + else + { + times_logarithmic = 1; + std::cout << "WARNING: could not get times_logarithmic from input file, defaults to "<< times_logarithmic << std::endl;} + + if (cfg.lookupValue("particle_diagnostics.bin_edge_0_time",bin_edge_0_time)) {} + else + { + bin_edge_0_time = std::log10(dt); + std::cout << "WARNING: could not get bin_edge_0_time from input file, defaults to "<< bin_edge_0_time << std::endl;} + + if (cfg.lookupValue("particle_diagnostics.bin_edge_1_time",bin_edge_1_time)) {} + else { std::cout << "ERROR: could not get bin_edge_1_time from input file " << std::endl;} + + if (cfg.lookupValue("particle_diagnostics.n_bins_time",n_bins_time)) {} + else { std::cout << "ERROR: could not get n_bins_time from input file " << std::endl;} + + // Set remainder of angle parameters + if (cfg.lookupValue("particle_diagnostics.angle_logarithmic",angle_logarithmic)) {} + else + { + angle_logarithmic = 0; + std::cout << "WARNING: could not get angle_logarithmic from input file, defaults to "<< angle_logarithmic << std::endl;} + + if (cfg.lookupValue("particle_diagnostics.bin_edge_0_angle",bin_edge_0_angle)) {} + else + { + bin_edge_0_angle = gitr_constants::pi / 2; + std::cout << "WARNING: could not get bin_edge_0_angle from input file, defaults to "<< bin_edge_0_angle << std::endl;} + + if (cfg.lookupValue("particle_diagnostics.bin_edge_1_angle",bin_edge_1_angle)) {} + else { std::cout << "ERROR: could not get bin_edge_1_angle from input file " << std::endl;} + + if (cfg.lookupValue("particle_diagnostics.n_bins_angle",n_bins_angle)) {} + else { std::cout << "ERROR: could not get n_bins_angle from input file " << std::endl;} + } + else + { + std::cout << "WARNING: No particle diagnostics information was found in input file" + << std::endl; + } + sim::Array bin_edges_time( n_bins_time + 1 , 0.0 ); + sim::Array bin_edges_angle( n_bins_angle + 1 , 0.0 ); + + if ( gitr_flags->USE_PARTICLE_DIAGNOSTICS ) + { + //creating histogram vector for all surfaces and bins + bin_edge_dt = (bin_edge_1_time - bin_edge_0_time)/n_bins_time; + bin_edge_dtheta = (bin_edge_1_angle - bin_edge_0_angle)/n_bins_angle; + + for (int i=0; i <= n_bins_time; i++) + { + bin_edges_time[i] = bin_edge_0_time + i*bin_edge_dt; + } + + for (int i=0; i <= n_bins_angle; i++) + { + bin_edges_angle[i] = bin_edge_0_angle + i*bin_edge_dtheta; + } + + } + + sim::Array histogram_particle_time(nSurfaces*n_bins_time,0.0); + sim::Array histogram_particle_angle(nSurfaces*n_bins_angle,0.0); + + particle_diagnostics particle_diagnostics0( + gitr_flags, particleArray, &boundaries[0], times_logarithmic, bin_edge_0_time, + bin_edge_1_time, bin_edge_dt, n_bins_time, &histogram_particle_time.front(), + angle_logarithmic, bin_edge_0_angle, bin_edge_1_angle, bin_edge_dtheta, + n_bins_angle, &histogram_particle_angle.front(), nSurfaces); + int subSampleFac = 1; if (world_rank == 0) { if (cfg.lookupValue("diagnostics.trackSubSampleFactor", subSampleFac)) { @@ -3920,6 +4018,8 @@ if( efield_interp == 1 ) sim::Array chargeHistoryGather(nHistories); sim::Array weightHistory(nHistories); sim::Array weightHistoryGather(nHistories); + sim::Array perpDistanceToSurfaceHistory(nHistories); + sim::Array perpDistanceToSurfaceGather(nHistories); if( particle_tracks > 0 ) { @@ -4131,7 +4231,7 @@ if( efield_interp == 1 ) &positionHistoryZ.front(), &velocityHistory.front(), &velocityHistoryX.front(), &velocityHistoryY.front(), &velocityHistoryZ.front(), &chargeHistory.front(), - &weightHistory.front()); + &weightHistory.front(), &perpDistanceToSurfaceHistory.front()); if( force_eval > 0 ) { @@ -4463,6 +4563,11 @@ std::cout << "here 2" << std::endl; thermalForce0); } + + if (gitr_flags->USE_PARTICLE_DIAGNOSTICS) + { + thrust::for_each(thrust::device, particleBegin, particleEnd, particle_diagnostics0); + } if( surface_model > 0 ) { thrust::for_each(thrust::device, particleBegin, particleEnd, reflection0); @@ -4603,6 +4708,9 @@ for(int i=0; icharge[world_rank * nP / world_size], nP / world_size, MPI_FLOAT, &chargeGather[0], nP / world_size, MPI_FLOAT, 0, MPI_COMM_WORLD); + MPI_Gather(&particleArray->perpDistanceToSurface[world_rank * nP / world_size], + nP / world_size, MPI_FLOAT, &perpDistanceToSurfaceGather[0], nP / world_size, + MPI_FLOAT, 0, MPI_COMM_WORLD); MPI_Gather(&particleArray->hasLeaked[world_rank * nP / world_size], nP / world_size, MPI_INT, &hasLeakedGather[0], nP / world_size, MPI_INT, 0, MPI_COMM_WORLD); @@ -4692,6 +4800,10 @@ for(int i=0; isurfaceHit[0]); nc_weight0.putVar(&particleArray->weight[0]); nc_charge0.putVar(&particleArray->charge[0]); + nc_perpDistanceToSurface0.putVar(&particleArray->perpDistanceToSurface[0]); nc_leak0.putVar(&particleArray->hasLeaked[0]); nc_dist0.putVar(&particleArray->distTraveled[0]); nc_time0.putVar(&particleArray->time[0]); @@ -5001,6 +5116,35 @@ std::cout << "bound 255 " << boundaries[255].impacts << std::endl; // std::cout << "weights " << " " << weights1[i] << " " << // weightThreshold[0] << std::endl; //} + if (gitr_flags->USE_PARTICLE_DIAGNOSTICS) +{ + netCDF::NcFile ncFile_particle_hist("output/particle_histograms.nc", netCDF::NcFile::replace); + netCDF::NcDim nc_ph_nSurfs = ncFile_particle_hist.addDim("nSurfaces", nSurfaces); + netCDF::NcDim nc_ph_n_bins_time = ncFile_particle_hist.addDim("n_bins_time", n_bins_time); + netCDF::NcDim nc_ph_n_edges_time = ncFile_particle_hist.addDim("n_edges_time", n_bins_time+1); + netCDF::NcDim nc_ph_n_bins_angle = ncFile_particle_hist.addDim("n_bins_angle", n_bins_angle); + netCDF::NcDim nc_ph_n_edges_angle = ncFile_particle_hist.addDim("n_edges_angle", n_bins_angle+1); + vector dims_part_hist_time; + vector dims_part_hist_angle; + dims_part_hist_time.push_back(nc_ph_nSurfs); + dims_part_hist_time.push_back(nc_ph_n_bins_time); + dims_part_hist_angle.push_back(nc_ph_nSurfs); + dims_part_hist_angle.push_back(nc_ph_n_bins_angle); + + netCDF::NcVar nc_ph_bin_edges_time = ncFile_particle_hist.addVar("bin_edges_time",netcdf_precision, nc_ph_n_edges_time); + nc_ph_bin_edges_time.putVar(&bin_edges_time[0]); + + netCDF::NcVar nc_ph_histogram_particle_time = ncFile_particle_hist.addVar("histogram_particle_time",netcdf_precision,dims_part_hist_time); + nc_ph_histogram_particle_time.putVar(&histogram_particle_time[0]); + + netCDF::NcVar nc_ph_bin_edges_angle = ncFile_particle_hist.addVar("bin_edges_angle",netcdf_precision, nc_ph_n_edges_angle); + nc_ph_bin_edges_angle.putVar(&bin_edges_angle[0]); + + netCDF::NcVar nc_ph_histogram_particle_angle = ncFile_particle_hist.addVar("histogram_particle_angle",netcdf_precision,dims_part_hist_angle); + nc_ph_histogram_particle_angle.putVar(&histogram_particle_angle[0]); + ncFile_particle_hist.close(); +} + if( surface_model > 0 || flux_ea > 0 ) { #if USE_MPI > 0 @@ -5149,6 +5293,8 @@ std::cout << "bound 255 " << boundaries[255].impacts << std::endl; netCDF::NcVar nc_charge = ncFile_hist.addVar("charge", netCDF::ncDouble, dims_hist); netCDF::NcVar nc_weight = ncFile_hist.addVar("weight", netCDF::ncDouble, dims_hist); + netCDF::NcVar nc_perpDistanceToSurface = ncFile_hist.addVar("perpDistanceToSurface", netCDF::ncDouble, dims_hist); + #if USE_MPI > 0 // if(world_rank ==0) //{ @@ -5169,6 +5315,7 @@ std::cout << "bound 255 " << boundaries[255].impacts << std::endl; nc_charge.putVar(&chargeHistoryGather[0]); nc_weight.putVar(&weightHistoryGather[0]); + nc_perpDistanceToSurface.putVar(&perpDistanceToSurfaceHistoryGather[0]); #else nc_x.putVar(&positionHistoryX[0]); nc_y.putVar(&positionHistoryY[0]); @@ -5179,6 +5326,8 @@ std::cout << "bound 255 " << boundaries[255].impacts << std::endl; nc_vz.putVar(&velocityHistoryZ[0]); nc_charge.putVar(&chargeHistory[0]); + nc_weight.putVar(&weightHistory[0]); + nc_perpDistanceToSurface.putVar(&perpDistanceToSurfaceHistory[0]); #endif ncFile_hist.close(); } diff --git a/src/particle_diagnostics.cpp b/src/particle_diagnostics.cpp new file mode 100644 index 00000000..bdecfb64 --- /dev/null +++ b/src/particle_diagnostics.cpp @@ -0,0 +1,90 @@ +#include "particle_diagnostics.h" +#include "spectroscopy.h" + + +particle_diagnostics::particle_diagnostics(Flags *_flags, + Particles *_particlesPointer, + Boundary *_boundaryVector, + bool _times_logarithmic, + gitr_precision _bin_edge_0_time, + gitr_precision _bin_edge_1_time, + gitr_precision _bin_edge_dt, + int _n_bins_time, + gitr_precision *_particle_time_histogram, + bool _angle_logarithmic, + gitr_precision _bin_edge_0_angle, + gitr_precision _bin_edge_1_angle, + gitr_precision _bin_edge_dtheta, + int _n_bins_angle, + gitr_precision *_particle_angle_histogram, + int _nSurfaces) + + : + + flags(_flags), particlesPointer(_particlesPointer), boundaryVector(_boundaryVector), times_logarithmic(_times_logarithmic), + bin_edge_0_time(_bin_edge_0_time), bin_edge_1_time(_bin_edge_1_time), bin_edge_dt(_bin_edge_dt), + n_bins_time(_n_bins_time), particle_time_histogram(_particle_time_histogram), angle_logarithmic(_angle_logarithmic), + bin_edge_0_angle(_bin_edge_0_angle), bin_edge_1_angle(_bin_edge_1_angle), bin_edge_dtheta(_bin_edge_dtheta), + n_bins_angle(_n_bins_angle), particle_angle_histogram(_particle_angle_histogram), nSurfaces(_nSurfaces) +{ } + +CUDA_CALLABLE_MEMBER_DEVICE +void particle_diagnostics::operator()(std::size_t indx) +{ + if (particlesPointer->hitWall[indx] == 1.0 && particlesPointer->weight[indx]!=0) + { + + int wallHit = particlesPointer->surfaceHit[indx]; + int surfaceHit = boundaryVector[wallHit].surfaceNumber; + + gitr_precision p_time = particlesPointer->time[indx] - particlesPointer->transitTime[indx]; + particlesPointer->transitTime[indx] = particlesPointer->time[indx]; + + int ind_time = std::floor((std::log10(p_time) - bin_edge_0_time)/bin_edge_dt); + if (times_logarithmic != 1) + { + ind_time = std::floor((p_time - bin_edge_0_time)/bin_edge_dt); + } + + if (ind_time >= n_bins_time) ind_time = n_bins_time - 1; + + int ind_2d_time = surfaceHit*n_bins_time + ind_time; + + if (ind_2d_time >=0 && ind_2d_time < nSurfaces*n_bins_time) + { + #if USE_CUDA > 0 + atomicAdd1(&particle_time_histogram[ind_2d_time],particlesPointer->weight[indx]); + #else + #pragma omp atomic + particle_time_histogram[ind_2d_time] = particle_time_histogram[ind_2d_time] + particlesPointer->weight[indx]; + #endif + } + + gitr_precision p_angle = particlesPointer->angle[indx] - particlesPointer->transitAngle[indx]; + //if (p_angle < 0.001) printf("DEBUG P_ANGLE: pangle %e || angle %e tangle %e || charge %e \n", p_angle, particlesPointer->angle[indx], particlesPointer->transitAngle[indx], particlesPointer->charge[indx]); + + particlesPointer->transitAngle[indx] = particlesPointer->angle[indx]; + + int ind_angle = std::floor((p_angle - bin_edge_0_angle)/bin_edge_dtheta); + if (angle_logarithmic == 1) + { + ind_angle = std::floor((std::log10(p_angle) - bin_edge_0_angle)/bin_edge_dtheta); + } + //if (p_angle < 0.001) printf("P_ANGLE: %e || IND_ANGLE: %d \n", p_angle, ind_angle); + + int ind_2d_angle = surfaceHit*n_bins_angle + ind_angle; + + if (ind_2d_angle >=0 && ind_2d_angle < nSurfaces*n_bins_angle) + { + //printf("IND_2D_ANGLE: %d || WEIGHT: %e \n",ind_2d_angle, particlesPointer->weight[indx]); + #if USE_CUDA > 0 + atomicAdd1(&particle_angle_histogram[ind_2d_angle],particlesPointer->weight[indx]); + #else + #pragma omp atomic + particle_angle_histogram[ind_2d_angle] = particle_angle_histogram[ind_2d_angle] + particlesPointer->weight[indx]; + #endif + } + particlesPointer->weight[indx] = 0; + } +} + diff --git a/src/utils.cpp b/src/utils.cpp index f10f7f39..da161e0f 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -577,6 +577,17 @@ int importGeometry(libconfig::Config &cfg_geom, sim::Array &boundaries else { + std::vector fields_required{"x1", "z1","x2","z2","slope","intercept","inDir","length","surface","Z"}; + for (int j = 0; j < fields_required.size(); j++) + { + + if (geom.exists(fields_required[j])) {} + else + { + std::cout << "ERROR: could not get field from geometry file "<< fields_required[j] << std::endl; + exit(0); + } + } //int nMaterials = geom["nMaterials"]; //std::cout << "nmat " << nMaterials << std::endl; for(int i=0 ; i