diff --git a/parameters.cpp b/parameters.cpp index 792710e76d..0b64be5533 100644 --- a/parameters.cpp +++ b/parameters.cpp @@ -149,7 +149,7 @@ int P::writeRestartAsFloat = false; string P::loadBalanceAlgorithm = string(""); std::map P::loadBalanceOptions; uint P::rebalanceInterval = numeric_limits::max(); - +int P::mapOrderShift=0; vector P::outputVariableList; vector P::diagnosticVariableList; @@ -214,6 +214,7 @@ vector P::blurPassString; vector P::numPasses; bool P::artificialPADiff; +int P::seed; Realf P::PADcoefficient; Realf P::PADCFL; int P::PADvbins; @@ -245,6 +246,8 @@ bool P::addParameters() { typedef Readparameters RP; // the other default parameters we read through the add/get interface RP::add("io.diagnostic_write_interval", "Write diagnostic output every arg time steps", numeric_limits::max()); + RP::add("map_order_shift","shift the map_order seed",0); + RP::add("map_order.seed","Scalar multiplier for map_order in vlasovmover",579450); RP::addComposing( "io.system_write_t_interval", @@ -608,6 +611,8 @@ bool P::addParameters() { void Parameters::getParameters() { typedef Readparameters RP; // get numerical values of the parameters + RP::get("map_order_shift",P::mapOrderShift); + RP::get("map_order.seed",P::seed); RP::get("io.diagnostic_write_interval", P::diagnosticInterval); RP::get("io.diagnostic_write_all_data_reducers", P::diagnosticWriteAllDROs); RP::get("io.system_write_t_interval", P::systemWriteTimeInterval); @@ -640,6 +645,9 @@ void Parameters::getParameters() { // Checks for validity of io and restart parameters int myRank; MPI_Comm_rank(MPI_COMM_WORLD, &myRank); + if (myRank==MASTER_RANK) { + std::cout << "mapOrderShift=" << P::mapOrderShift << std::endl; + } const string prefix = string("./"); if (access(&(P::restartWritePath[0]), W_OK) != 0) { if (myRank == MASTER_RANK) { diff --git a/parameters.h b/parameters.h index 8f0c59cca5..7f358659e1 100644 --- a/parameters.h +++ b/parameters.h @@ -193,6 +193,7 @@ struct Parameters { static bool refineOnRestart; static bool forceRefinement; static bool shouldFilter; + static int seed; static bool useAlpha1; static Real alpha1RefineThreshold; static Real alpha1CoarsenThreshold; @@ -231,6 +232,7 @@ struct Parameters { static std::vector amrBoxCenterY; static std::vector amrBoxCenterZ; static std::vector amrBoxMaxLevel; + static int mapOrderShift; static bool amrTransShortPencils; /*!< Use short or long pencils in AMR translation.*/ static std::vector blurPassString; static std::vector numPasses; diff --git a/testpackage/run_tests.sh b/testpackage/run_tests.sh index d652d9bbbe..3079e86d58 100755 --- a/testpackage/run_tests.sh +++ b/testpackage/run_tests.sh @@ -5,16 +5,16 @@ # define tab interval sequence so that we have aligned output # this is now covering at least up proton/vg_ptensor_nonthermal_offdiagonal_0 and numbers printed at setprecision(3) with negative mantissa and exponent tabseq="1,46,62,78,94,110" -tabs $tabseq &> /dev/null # suppress special character output, list matches expand below - +tabs $tabseq &>/dev/null # suppress special character output, list matches expand below +run_dir="run" ## add absolute paths to folder names, filenames -reference_dir=$( readlink -f $reference_dir ) -run_dir=$( readlink -f $run_dir )_$( date +%Y.%m.%d_%H.%M.%S) -reference_revision_parsed=$( readlink -f $reference_dir/$reference_revision ) -bin=$( readlink -f $bin ) -diffbin=$( readlink -f $diffbin ) -test_dir=$( readlink -f $test_dir) +reference_dir=$(readlink -f $reference_dir) +run_dir=$(readlink -f $run_dir)_$(date +%H%M%S)_$$ +reference_revision_parsed=$(readlink -f $reference_dir/$reference_revision) +# bin=$( readlink -f "$bin" ) +diffbin=$(readlink -f $diffbin) +test_dir=$(readlink -f $test_dir) # for run in ${run_tests[*]} # do @@ -23,189 +23,169 @@ test_dir=$( readlink -f $test_dir) # done if [[ ! $small_run_command ]]; then - echo "No small_run_command provided in machine config, please update it!" - exit + echo "No small_run_command provided in machine config, please update it!" + exit fi -flags=$( $run_command $bin --version |grep CXXFLAGS) -solveropts=$(echo $flags|sed 's/[-+]//g' | gawk '{for(i = 1;i<=NF;i++) { if( $i=="DDP" || $i=="DFP" || index($i,"PF")|| index($i,"DVEC") || index($i,"SEMILAG") ) printf "__%s", $(i) }}') -revision=$( $run_command $bin --version |gawk '{if(flag==1) {print $1;flag=0}if ($3=="log") flag=1;}' ) - -if [ $create_verification_files == 1 ] -then - #if we create the references, then lets simply run in the reference dir and turn off tests below. Revision is - #automatically obtained from the --version output (this overwrites the reverence_revision variable from a launch script) - reference_revision=${revision}${solveropts} - echo "Computing reference results into ${reference_dir}/${reference_revision}" - if [[ ! -z $GITHUB_ACTIONS ]] - then - echo "REFERENCE_REVISION=${reference_dir}/${reference_revision}" >> "$GITHUB_ENV" - fi +flags=$($run_command $bin --version | grep CXXFLAGS) +solveropts=$(echo $flags | sed 's/[-+]//g' | gawk '{for(i = 1;i<=NF;i++) { if( $i=="DDP" || $i=="DFP" || index($i,"PF")|| index($i,"DVEC") || index($i,"SEMILAG") ) printf "__%s", $(i) }}') +revision=$($run_command $bin --version | gawk '{if(flag==1) {print $1;flag=0}if ($3=="log") flag=1;}') + +if [ $create_verification_files == 1 ]; then + #if we create the references, then lets simply run in the reference dir and turn off tests below. Revision is + #automatically obtained from the --version output (this overwrites the reverence_revision variable from a launch script) + reference_revision=${revision}${solveropts} + echo "Computing reference results into ${reference_dir}/${reference_revision}" + if [[ ! -z $GITHUB_ACTIONS ]]; then + echo "REFERENCE_REVISION=${reference_dir}/${reference_revision}" >>"$GITHUB_ENV" + fi else - echo "----------" - echo "This will be verifying ${revision}_$solveropts against $reference_revision_parsed" - echo "----------" + echo "----------" + echo "This will be verifying ${revision}_$solveropts against $reference_revision_parsed" + echo "----------" fi -if [ -d $run_dir ] -then - echo $run_dir exists? - exit +if [ -d $run_dir ]; then + echo $run_dir exists? + exit fi mkdir -p $run_dir # loop over different test cases -for run in ${run_tests[*]} -do - echo -e "\n" - echo "----------" - echo "running ${test_name[$run]} " -# directory for test results - vlsv_dir=${run_dir}/${test_name[$run]} - cfg_dir=${test_dir}/${test_name[$run]} - -# Check if folder for new run exists, if not create them, otherwise delete old results - if [ ! -d ${vlsv_dir} ]; then - mkdir -p ${vlsv_dir} - else - rm -f ${vlsv_dir}/* +for run in ${run_tests[*]}; do + echo -e "\n" + echo "----------" + echo "running ${test_name[$run]} " + # directory for test results + vlsv_dir=${run_dir}/${test_name[$run]} + cfg_dir=${test_dir}/${test_name[$run]} + + # Check if folder for new run exists, if not create them, otherwise delete old results + if [ ! -d ${vlsv_dir} ]; then + mkdir -p ${vlsv_dir} + else + rm -f ${vlsv_dir}/* + fi + + # change to run directory of the test case, e.g. test_Fluctuations + cd ${vlsv_dir} + cp ${cfg_dir}/* . + + # export OMP_NUM_THREADS=$t + export MPICH_MAX_THREAD_SAFETY=funneled + + # Run prerequisite script, if it exists + test -e test_prelude.sh && ./test_prelude.sh + + # Run the actual simulation + if [[ ${single_cell[$run]} ]]; then + $small_run_command $bin --version >VERSION.txt + $small_run_command $bin --run_config=${test_name[$run]}.cfg + else + $run_command $bin --version >VERSION.txt + $run_command $bin --run_config=${test_name[$run]}.cfg + fi + + # Run postprocessing script, if it exists + test -e test_postproc.sh && ./test_postproc.sh + + ###copy new reference data to correct folder + if [ $create_verification_files == 1 ]; then + reference_result_dir=${reference_dir}/${reference_revision}/${test_name[$run]} + if [ -e $reference_result_dir ]; then + echo "Removing previous reference results" + rm -rf $reference_result_dir fi -# change to run directory of the test case, e.g. test_Fluctuations - cd ${vlsv_dir} - cp ${cfg_dir}/* . + mkdir -p $reference_result_dir + cp * $reference_result_dir + fi - export OMP_NUM_THREADS=$t - export MPICH_MAX_THREAD_SAFETY=funneled + cd $base_dir - # Run prerequisite script, if it exists - test -e test_prelude.sh && ./test_prelude.sh + ### TESTS ##### + if [ ! $create_verification_files == 1 ]; then + ##Compare test case with right solutions + echo "--------------------------------------------------------------------------------------------" + echo "${test_name[$run]} - Verifying ${revision}_$solveropts against $reference_revision" + echo "--------------------------------------------------------------------------------------------" + reference_result_dir=${reference_dir}/${reference_revision}/${test_name[$run]} - # Run the actual simulation - if [[ ${single_cell[$run]} ]]; then - $small_run_command $bin --version > VERSION.txt - $small_run_command $bin --run_config=${test_name[$run]}.cfg + #print header + + echo "------------------------------------------------------------" + echo " ref-time | new-time | speedup |" + echo "------------------------------------------------------------" + if [ -e ${reference_result_dir}/${comparison_phiprof[$run]} ]; then + refPerf=$(grep "Propagate " ${reference_result_dir}/${comparison_phiprof[$run]} | gawk '(NR==1){print $11}') else - $run_command $bin --version > VERSION.txt - $run_command $bin --run_config=${test_name[$run]}.cfg + refPerf="NA" fi - - # Run postprocessing script, if it exists - test -e test_postproc.sh && ./test_postproc.sh - - - ###copy new reference data to correct folder - if [ $create_verification_files == 1 ] - then - reference_result_dir=${reference_dir}/${reference_revision}/${test_name[$run]} - if [ -e $reference_result_dir ] - then - echo "Removing previous reference results" - rm -rf $reference_result_dir - fi - - mkdir -p $reference_result_dir - cp * $reference_result_dir + if [ -e ${vlsv_dir}/${comparison_phiprof[$run]} ]; then + newPerf=$(grep "Propagate " ${vlsv_dir}/${comparison_phiprof[$run]} | gawk '(NR==1){print $11}') + else + newPerf="NA" fi - - cd $base_dir - - - -### TESTS ##### - if [ ! $create_verification_files == 1 ] - then -##Compare test case with right solutions - echo "--------------------------------------------------------------------------------------------" - echo "${test_name[$run]} - Verifying ${revision}_$solveropts against $reference_revision" + #print speedup if both refPerf and newPerf are numerical values + speedup=$(echo $refPerf $newPerf | gawk '{if($2 == $2 + 0 && $1 == $1 + 0 ) print $1/$2; else print "NA"}') + + tabs 1,14,33,59 &>/dev/null # match next line + echo -e " $refPerf\t| $newPerf\t| $speedup\t|" | expand -t 1,14,33,59 # match previous line + echo "------------------------------------------------------------" + tabs $tabseq &>/dev/null # reset for other printouts + echo -e " variable\t| absolute diff\t| relative diff |" | expand -t $tabseq # list matches tabs above + echo "------------------------------------------------------------" + + variables=(${variable_names[$run]// / }) + indices=(${variable_components[$run]// / }) + for vlsv in ${comparison_vlsv[$run]}; do + if [ ! -f "${vlsv_dir}/${vlsv}" ]; then + echo "Output file ${vlsv_dir}/${vlsv} not found!" echo "--------------------------------------------------------------------------------------------" - reference_result_dir=${reference_dir}/${reference_revision}/${test_name[$run]} - - #print header - - - echo "------------------------------------------------------------" - echo " ref-time | new-time | speedup |" - echo "------------------------------------------------------------" - if [ -e ${reference_result_dir}/${comparison_phiprof[$run]} ] - then - refPerf=$(grep "Propagate " ${reference_result_dir}/${comparison_phiprof[$run]} |gawk '(NR==1){print $11}') - else - refPerf="NA" - fi - if [ -e ${vlsv_dir}/${comparison_phiprof[$run]} ] - then - newPerf=$(grep "Propagate " ${vlsv_dir}/${comparison_phiprof[$run]} |gawk '(NR==1){print $11}') - else - newPerf="NA" - fi - #print speedup if both refPerf and newPerf are numerical values - speedup=$( echo $refPerf $newPerf |gawk '{if($2 == $2 + 0 && $1 == $1 + 0 ) print $1/$2; else print "NA"}') - - tabs 1,14,33,59 &> /dev/null # match next line - echo -e " $refPerf\t| $newPerf\t| $speedup\t|" | expand -t 1,14,33,59 # match previous line - echo "------------------------------------------------------------" - tabs $tabseq &> /dev/null # reset for other printouts - echo -e " variable\t| absolute diff\t| relative diff |" | expand -t $tabseq # list matches tabs above - echo "------------------------------------------------------------" - - variables=(${variable_names[$run]// / }) - indices=(${variable_components[$run]// / }) - for vlsv in ${comparison_vlsv[$run]} - do - if [ ! -f "${vlsv_dir}/${vlsv}" ]; then - echo "Output file ${vlsv_dir}/${vlsv} not found!" - echo "--------------------------------------------------------------------------------------------" - continue - fi - if [ ! -f "${reference_result_dir}/${vlsv}" ]; then - echo "Reference file ${reference_result_dir}/${vlsv} not found!" - echo "--------------------------------------------------------------------------------------------" - continue - fi - echo "Comparing file ${vlsv_dir}/${vlsv} against reference" - for i in ${!variables[*]} - do - if [[ "${variables[$i]}" == "fg_"* ]] - then - A=$( $run_command_tools $diffbin --meshname=fsgrid ${reference_result_dir}/${vlsv} ${vlsv_dir}/${vlsv} ${variables[$i]} ${indices[$i]} ) - relativeValue=$(grep "The relative 0-distance between both datasets" <<< $A |gawk '{print $8}' ) - absoluteValue=$(grep "The absolute 0-distance between both datasets" <<< $A |gawk '{print $8}' ) - #print the results - echo -e " ${variables[$i]}_${indices[$i]}\t ${absoluteValue}\t ${relativeValue}" | expand -t $tabseq #list matches tabs above - elif [[ "${variables[$i]}" == "ig_"* ]] - then - B=$( $run_command_tools $diffbin --meshname=ionosphere ${reference_result_dir}/${vlsv} ${vlsv_dir}/${vlsv} ${variables[$i]} ${indices[$i]} ) - relativeValue=$(grep "The relative 0-distance between both datasets" <<< $B |gawk '{print $8}' ) - absoluteValue=$(grep "The absolute 0-distance between both datasets" <<< $B |gawk '{print $8}' ) - #print the results - echo -e " ${variables[$i]}_${indices[$i]}\t ${absoluteValue}\t ${relativeValue}" | expand -t $tabseq # list matches tabs above - elif [ ! "${variables[$i]}" == "proton" ] - then # Regular vg_ variable - C=$( $run_command_tools $diffbin ${reference_result_dir}/${vlsv} ${vlsv_dir}/${vlsv} ${variables[$i]} ${indices[$i]} ) - relativeValue=$(grep "The relative 0-distance between both datasets" <<< $C |gawk '{print $8}' ) - absoluteValue=$(grep "The absolute 0-distance between both datasets" <<< $C |gawk '{print $8}' ) - #print the results - echo -e " ${variables[$i]}_${indices[$i]}\t ${absoluteValue}\t ${relativeValue}" | expand -t $tabseq # list matches tabs above - elif [ "${variables[$i]}" == "proton" ] - then - echo "--------------------------------------------------------------------------------------------" - echo " Distribution function diff " - echo "--------------------------------------------------------------------------------------------" - $run_command_tools $diffbin ${reference_result_dir}/${vlsv} ${vlsv_dir}/${vlsv} proton 0 - fi - done # loop over variables - - # Print also time difference, if it is not zero - timeDiff=$(grep "delta t" <<< $C |gawk '{print $8}' ) - if (( $(awk 'BEGIN{print ('$timeDiff'!= 0.0)?1:0}') )) - then - echo "WARNING! VLSV file timestamps differ by ${timeDiff}s." - else - echo "VLSV file timestamps match." - fi - echo "--------------------------------------------------------------------------------------------" - done # loop over vlsv files to compare + continue + fi + if [ ! -f "${reference_result_dir}/${vlsv}" ]; then + echo "Reference file ${reference_result_dir}/${vlsv} not found!" echo "--------------------------------------------------------------------------------------------" - fi + continue + fi + echo "Comparing file ${vlsv_dir}/${vlsv} against reference" + for i in ${!variables[*]}; do + if [[ "${variables[$i]}" == "fg_"* ]]; then + A=$($run_command_tools $diffbin --meshname=fsgrid ${reference_result_dir}/${vlsv} ${vlsv_dir}/${vlsv} ${variables[$i]} ${indices[$i]}) + relativeValue=$(grep "The relative 0-distance between both datasets" <<<$A | gawk '{print $8}') + absoluteValue=$(grep "The absolute 0-distance between both datasets" <<<$A | gawk '{print $8}') + #print the results + echo -e " ${variables[$i]}_${indices[$i]}\t ${absoluteValue}\t ${relativeValue}" | expand -t $tabseq #list matches tabs above + elif [[ "${variables[$i]}" == "ig_"* ]]; then + B=$($run_command_tools $diffbin --meshname=ionosphere ${reference_result_dir}/${vlsv} ${vlsv_dir}/${vlsv} ${variables[$i]} ${indices[$i]}) + relativeValue=$(grep "The relative 0-distance between both datasets" <<<$B | gawk '{print $8}') + absoluteValue=$(grep "The absolute 0-distance between both datasets" <<<$B | gawk '{print $8}') + #print the results + echo -e " ${variables[$i]}_${indices[$i]}\t ${absoluteValue}\t ${relativeValue}" | expand -t $tabseq # list matches tabs above + elif [ ! "${variables[$i]}" == "proton" ]; then # Regular vg_ variable + C=$($run_command_tools $diffbin ${reference_result_dir}/${vlsv} ${vlsv_dir}/${vlsv} ${variables[$i]} ${indices[$i]}) + relativeValue=$(grep "The relative 0-distance between both datasets" <<<$C | gawk '{print $8}') + absoluteValue=$(grep "The absolute 0-distance between both datasets" <<<$C | gawk '{print $8}') + #print the results + echo -e " ${variables[$i]}_${indices[$i]}\t ${absoluteValue}\t ${relativeValue}" | expand -t $tabseq # list matches tabs above + elif [ "${variables[$i]}" == "proton" ]; then + echo "--------------------------------------------------------------------------------------------" + echo " Distribution function diff " + echo "--------------------------------------------------------------------------------------------" + $run_command_tools $diffbin ${reference_result_dir}/${vlsv} ${vlsv_dir}/${vlsv} proton 0 + fi + done # loop over variables + + # Print also time difference, if it is not zero + timeDiff=$(grep "delta t" <<<$C | gawk '{print $8}') + if (($(awk 'BEGIN{print ('$timeDiff'!= 0.0)?1:0}'))); then + echo "WARNING! VLSV file timestamps differ by ${timeDiff}s." + else + echo "VLSV file timestamps match." + fi + echo "--------------------------------------------------------------------------------------------" + done # loop over vlsv files to compare + echo "--------------------------------------------------------------------------------------------" + fi done # loop over tests diff --git a/testpackage/small_test_carrington.sh b/testpackage/small_test_carrington.sh index e5ee25e85a..93561918fa 100755 --- a/testpackage/small_test_carrington.sh +++ b/testpackage/small_test_carrington.sh @@ -1,34 +1,41 @@ #!/bin/bash -#SBATCH -t 01:30:00 # Run time (hh:mm:ss) -#SBATCH --job-name=ctestpackage +#SBATCH -t 00:30:00 # Run time (hh:mm:ss) +#SBATCH --job-name=repeat_test_acc ##SBATCH -A spacephysics -#SBATCH --constraint="carrington" +#SBATCH --constraint="amd" # test short medium 20min1d 3d #SBATCH -p short #SBATCH --exclusive -#SBATCH --nodes=1 -#SBATCH -c 4 # CPU cores per task -#SBATCH -n 16 # number of tasks -#SBATCH --mem-per-cpu=5G -#SBATCH --hint=multithread +##SBATCH --nodes=1 +#SBATCH -c 1 # CPU cores per task +#SBATCH -n 1 # number of tasks +#SBATCH --array=0-9 +#SBATCH --mem=10G +#SBATCH --no-requeue +##SBATCH --hint=multithread # If 1, the reference vlsv files are generated # if 0 then we check the v1 against reference files create_verification_files=0 # folder for all reference data -reference_dir="/turso/group/spacephysics/vlasiator/testpackage/" +reference_dir="/home/siclasse/" +mainfolder="/home/siclasse/vlasiator_hile/vlasiator/testpackage/" cd $SLURM_SUBMIT_DIR -bin="/proj/USERNAME/BINARYNAME" -diffbin="/turso/group/spacephysics/vlasiator/testpackage/vlsvdiff_DP_carrington" +diffbin="../vlsvdiff_DP" #compare agains which revision #reference_revision="CI_reference" -reference_revision="current" +reference_revision="rng_ref2" -source ../modules/carrington_gcc_openmpi.sh -# export UCX_NET_DEVICES=eth5,mlx5_0:1 # This is important for multi-node performance! +module purge +module load GCC/13.2.0 +module load OpenMPI/4.1.6-GCC-13.2.0 +module load PMIx/4.2.6-GCCcore-13.2.0 +module load PAPI/7.1.0-GCCcore-13.2.0 +module load Boost/1.83.0-GCC-13.2.0 +module load xthi export UCX_TLS=dc_mlx5 export UCX_NET_DEVICES=mlx5_0:1 @@ -37,17 +44,16 @@ export OMPI_MCA_pml='ucx' export OMPI_MCA_mtl='^ofi' #Carrington has 2 x 16 cores per node, plus hyperthreading ht=2 -t=$SLURM_CPUS_PER_TASK -export OMP_NUM_THREADS=$t - +# t=$SLURM_CPUS_PER_TASK +# export OMP_NUM_THREADS=$t #command for running stuff -run_command="srun --mpi=pmix_v3 -n $SLURM_NTASKS " -small_run_command="srun --mpi=pmix_v3 -n 1" +run_command="mpirun -n 1 " +small_run_command="mpirun -n 1 " run_command_tools="mpirun -np 1 " umask 007 # Launch the OpenMP job to the allocated compute node -echo "Running $exec on $SLURM_NTASKS mpi tasks, with $t threads per task on $SLURM_NNODES nodes ($ht threads per physical core)" +# echo "Running $exec on $SLURM_NTASKS mpi tasks, with $t threads per task on $SLURM_NNODES nodes ($ht threads per physical core)" # Optional debug printouts # srun -np 1 /appl/bin/hostinfo @@ -58,5 +64,19 @@ echo "Running $exec on $SLURM_NTASKS mpi tasks, with $t threads per task on $SLU source test_definitions_small.sh wait # Run tests -source run_tests.sh -wait +jobcount=$(($SLURM_ARRAY_TASK_MAX - $SLURM_ARRAY_TASK_MIN + 1)) +index=$(($SLURM_ARRAY_TASK_ID - $SLURM_ARRAY_TASK_MIN)) + +# echo $index $jobcount $SLURM_ARRAY_TASK_COUNT $SLURM_ARRAY_TASK_ID +# echo $index +for n in $(seq $(($index * 50)) $((($index + 1) * 50 - 1))); do + echo "N=$n" + bin="/home/siclasse/vlasiator_hile/vlasiator/vlasiator --map_order_shift $n" + source run_tests.sh + + wait + # mv logfile.txt logfile_$n.txt + cd $mainfolder + # cd $SLURM_SUBMIT_DIR + wait +done diff --git a/testpackage/test_definitions_small.sh b/testpackage/test_definitions_small.sh index 7cb6cdd1b7..40d0ed6025 100644 --- a/testpackage/test_definitions_small.sh +++ b/testpackage/test_definitions_small.sh @@ -26,7 +26,7 @@ index=1 # 1 basic multipeak acceleration test (fixed timestep) test_name[${index}]="acctest_1_maxw_500k_100k_20kms_10deg" -comparison_vlsv[${index}]="fullf.0000000.vlsv fullf.0000001.vlsv fullf.0000002.vlsv fullf.0000020.vlsv" +comparison_vlsv[${index}]="fullf.0000020.vlsv" # fullf.0000001.vlsv fullf.0000002.vlsv fullf.0000020.vlsv" comparison_phiprof[${index}]="phiprof_0.txt" variable_names[${index}]="proton/vg_rho proton/vg_v proton/vg_v proton/vg_v proton" variable_components[${index}]="0 0 1 2" @@ -268,4 +268,4 @@ do done # Alternatively, set tests manually, e.g. -# run_tests=( 1 6 9 ) +run_tests=( 1 ) diff --git a/testpackage/tests/acctest_1_maxw_500k_100k_20kms_10deg/acctest_1_maxw_500k_100k_20kms_10deg.cfg b/testpackage/tests/acctest_1_maxw_500k_100k_20kms_10deg/acctest_1_maxw_500k_100k_20kms_10deg.cfg index a038833162..ebf4bb7f60 100644 --- a/testpackage/tests/acctest_1_maxw_500k_100k_20kms_10deg/acctest_1_maxw_500k_100k_20kms_10deg.cfg +++ b/testpackage/tests/acctest_1_maxw_500k_100k_20kms_10deg/acctest_1_maxw_500k_100k_20kms_10deg.cfg @@ -10,7 +10,7 @@ ParticlePopulations = proton diagnostic_write_interval = 1 write_initial_state = 0 -system_write_t_interval = 180 +system_write_t_interval = 3600 system_write_file_name = fullf system_write_distribution_stride = 1 system_write_distribution_xline_stride = 0 @@ -43,7 +43,7 @@ y_min = 0.0 y_max = 1.0e6 z_min = 0 z_max = 1.0e6 -t_max = 3600 +t_max = 3601 dt = 10.0 [proton_properties] @@ -58,9 +58,9 @@ vy_min = -2.0e6 vy_max = +2.0e6 vz_min = -2.0e6 vz_max = +2.0e6 -vx_length = 50 -vy_length = 50 -vz_length = 50 +vx_length = 100 +vy_length = 100 +vz_length = 100 [proton_sparse] minValue = 1.0e-16 @@ -71,9 +71,12 @@ periodic_z = yes [MultiPeak] #magnitude of 1.82206867e-10 gives a period of 360s, useful for testing... -Bx = 1.2e-10 -By = 0.8e-10 -Bz = 1.1135233442526334e-10 +# Bx = 1.2e-10 +# By = 0.8e-10 +# Bz = 1.1135233442526334e-10 +Bx = 0 +By = 0 +Bz = 1.82206867e-10 magXPertAbsAmp = 0 magYPertAbsAmp = 0 magZPertAbsAmp = 0 diff --git a/testpackage/vdf_statisticator.py b/testpackage/vdf_statisticator.py new file mode 100644 index 0000000000..025c48aa46 --- /dev/null +++ b/testpackage/vdf_statisticator.py @@ -0,0 +1,243 @@ + + +import numpy as np +import scipy.ndimage as ndimage +import analysator as pt +import matplotlib.pyplot as plt +import os +import scipy.stats + +#text file of the run names +#I used run_[timestamp]_[processid], so i made the list with `ls -l | grep -Po "run_\d+_\d+` > runlist.txt +#within those run folders are the test folders +runs = np.loadtxt("/home/siclasse/vlasiator_hile/vlasiator/testpackage/acctest500/runlist.txt", object) + +run_title="acctest" +do_slices=True #plot mean,std,snr vdf in slices not using analysator +do_slices_analysator=False #this should maybe be done with ohters disabled, uses analysator to plot the slices of the vdf instead +do_integrated=True #use analysator to get the integrated grid (note that i modified analysator to return the mesh and bins earlier in plot_vdf()) +do_peak=True #get the peak indices of the peak value for the vdf from the first run and make historgram of the vdf value in that cell +do_statistics=True #histogram of bulkv(vg_v),proton/numrho(proton/vg_rho) and bulkrho (vg_rhom) + +i = 0 +dit = {} +grid_init = [] +all_a = [] +max_i = len(runs) #for manually setting how many runs to plot +# max_i = 6 +peak_id = -1 +peak_vals = [] +peak_coord_og = np.array([-1, -1, -1]) + +#I love eating RAM, nom nom +grid_size=[400,400,400] +if do_integrated: + bins_runs=np.full((max_i,grid_size[0],grid_size[1]),0,dtype=np.float32) +if do_slices: + fullgrid = np.full((max_i, *grid_size), 0, dtype=np.float32) + +bulkv=[] +bulkrho=np.array([]) +numrho=np.array([]) + +filename="fullf" +for d,run in enumerate(runs): + ref = f"/home/siclasse/vlasiator_hile/vlasiator/testpackage/acctest500/{run}/acctest_1_maxw_500k_100k_20kms_10deg" + print(ref) + if not os.path.isfile(ref+f"/{filename}.0000001.vlsv"): + print(f"File {ref}/{filename} not found") + continue + file = pt.vlsvfile.VlsvReader(ref + f"/{filename}.0000001.vlsv",indexer="dict") + R_EARTH = 6.371e6 + x = 15 + y = 0 + z = 0 + + cellid=1 + # cellid = file.get_cellid_with_vdf([x*R_EARTH, y*R_EARTH, z*R_EARTH]) + if do_statistics: + bulkv.append(np.array(file.read_variable("vg_v",cellids=cellid))) + bulkrho=np.append(bulkrho,file.read_variable("vg_rhom",cellids=cellid)) + numrho=np.append(numrho,file.read_variable("proton/vg_rho",cellids=cellid)) + + print(run,cellid) + pop = "proton" + velocity_cell_map = file.read_velocity_cells(cellid, pop) + maps = list(zip(*velocity_cell_map.items())) + velocity_cell_ids = np.array(maps[0], dtype=np.int64) + velocity_cell_values = np.array(maps[1], dtype=np.float32) + + velocity_cell_coordinates = file.get_velocity_cell_coordinates( + velocity_cell_ids, pop + ) + velocity_cell_indices = file.get_velocity_cell_indices( + velocity_cell_coordinates, pop + ) + + if i == 0: + dv = (velocity_cell_coordinates[1] - velocity_cell_coordinates[0])[0] + peak_id = np.argmax(velocity_cell_values) + peak_coord = velocity_cell_coordinates[peak_id] + peak_highest_corner = np.max(velocity_cell_indices, axis=0) + peak_lowest_corner = np.min(velocity_cell_indices, axis=0) + peak_ind = velocity_cell_indices[peak_id] + peak_id = velocity_cell_ids[peak_id] + slicevecs=[[0,0,peak_coord[2]+z*dv] for z in range(-10,10)] + + #note that the analysator plot_vdf function was modified to return these values after they are set the final time + if do_integrated: + xmesh,ymesh,bins=pt.plot.plot_vdf(vlsvobj=file,cellids=[cellid],outputfile=f'./plot_vdf_{d}.png',xy=1,slicethick=1,fmin=1e-15,fmax=4e-9,box=np.array([-500,500,-1000,1000])*10**3) + bins_runs[i]=bins + plt.clf() + + if do_slices_analysator: + for slicevec in slicevecs: + os.system(f"mkdir slices_{slicevec[2]}_analysator_plot") + pt.plot.plot_vdf(vlsvobj=file,cellids=[1],outputfile=f'./slices_{slicevec[2]}_analysator_plot/plot_vdf_{slicevec[2]}_{d}_{run_title}.png',center=slicevec,xy=1,slicethick=0.5,box=np.array([-500,500,-1000,1000])*10**3,reducer="average",fmin=1e-15,fmax=4e-9) + + if do_peak: + peak_vals.append( + velocity_cell_values[np.argwhere(velocity_cell_ids == peak_id)][0][0] + ) + + if do_slices: + for id, cellindices in enumerate(velocity_cell_indices): + ind = np.asarray(cellindices) + fullgrid[i, ind[0],ind[1],ind[2]] = velocity_cell_values[id] + + i += 1 + if i == max_i: + break + + + +if do_statistics: + bulkv=np.array(bulkv) + for name,val in [("numrho",numrho),("bulkv_x",bulkv[:,0]),("bulkv_y",bulkv[:,1]),("bulkv_z",bulkv[:,2]),("bulkrho",bulkrho*10**25)]: + plt.title(f"{run_title} {name} hist") + + statistic,pval=scipy.stats.shapiro(val) + plt.hist(val, bins="auto", density=False,label=f"shapiro-wilk={statistic:.3f},pval={pval:.3f},mean={np.mean(val):.3f},std={np.std(val,ddof=1):.3f}") + #if chosen alpha level (the level at which we are comfortable falsely rejecting null(data is normally distributed) is above p-value, we may not reject the possibility that it is normally distributed) + #note that this is not the same as to say it is normally distributed or not. + + #had to scale with e-25 for some so the shapiro-wilk didn't complain SHOULD BE CHANGED IF USING DIFFERENT DATA + + plt.gca().set_xlabel("bins e-25" if name=="bulkrho" else "bins") + plt.gca().set_ylabel("counts") + plt.gcf().set_size_inches(8,6) + plt.legend() + plt.savefig(f"{name}_hist_{run_title}.png") + plt.clf() + +if do_peak: + plt.hist(peak_vals, bins="auto", density=False) + plt.savefig(f"./peak_hist_{run_title}_{i}_peakid_{peak_id}.png") + plt.clf() + + +if do_slices: + std = np.nanstd(fullgrid, axis=0, ddof=1) #should not have nans but whatever +# std = np.nan_to_num(std) + + mean = np.nanmean(fullgrid, axis=0) + print("std maximum ",np.max(std)) + tick_n = 10 + xtick = ( + np.arange(0, grid_size[0], 1) * dv + file.get_velocity_mesh_extent()[0:3][0] + dv + ) + xtick = [f"{tick * 10**-3:.0f}" for tick in xtick] + snr = np.full(std.shape, np.nan) + np.divide(mean, std, where=std > 0, out=snr) + +# snr_peak_coords=[] +# peak_snr=np.array([-1,-1,-1]) + +if do_integrated: + mean_inte=np.mean(bins_runs,axis=0) + std_inte=np.std(bins_runs,axis=0) + + snr_inte = np.full(std_inte.shape,0,dtype=np.float32) + np.divide(mean_inte, std_inte, where=std_inte > 0, out=snr_inte) + +for scaling in ["linear", "log"]: + if do_integrated: + for name, data in [("std", std_inte), ("mean", mean_inte), ("snr", snr_inte)]: + plt.gcf().set_size_inches(8,8) + plt.gca().set_xlim(-800,800) + plt.gca().set_ylim(-800,800) + plt.pcolormesh(xmesh,ymesh,data,norm=scaling) + + plt.colorbar(orientation="horizontal") + plt.title( + f"{name} of VDF over {i} runs, scaling {scaling}, run {run_title}" + ) + + fig = plt.gcf() + ax = plt.gca() + + plt.xlabel("v_x [m/s]") + plt.xticks(rotation=45) + plt.ylabel("v_y [m/s]") + plt.savefig(f"./integrated_{name}_{i}_{scaling}_{run_title}.png", dpi=300) + ax.clear() + fig.clear() + + if do_slices: + for slice_i in range(peak_ind[2]-10, peak_ind[2]+10): + for name, data in [("std", std), ("mean", mean), ("snr", snr)]: + if np.max(data[:, ::-1, slice_i].T) != np.min(data[:, ::-1, slice_i].T): + cutoff = 1 * (mean[:, :, slice_i] > 10 ** (-16)) + + x, y = np.meshgrid( + np.arange(0, grid_size[0]), np.arange(0, grid_size[1]) + ) + plt.contour( + x, + y, + cutoff.T, + 0, + colors="black", + linewidths=0.5, + linestyles="dashed", + ) + + #following was for tracking the SNR blob + # if name=="snr" and scaling=="linear": + # peak_snr=np.argwhere(data[:,:,slice_i].T==np.nanmax(data[:,:,slice_i].T))[0] + # snr_peak_coords.append([peak_snr[0],peak_snr[1],slice_i]) + # continue + # else: + # continue + # + + plt.imshow(data[:, :, slice_i].T, norm=scaling) + + plt.colorbar(orientation="horizontal") + plt.title( + f"{name} of VDF over {i} runs, scaling {scaling} slice z={slice_i}, run {run_title}" + ) + + fig = plt.gcf() + ax = plt.gca() + ax.set_xticks( + np.linspace(0.5, grid_size[0] - 0.5, grid_size[0])[::tick_n] + ) + ax.set_yticks( + np.linspace(-0.5, grid_size[1] - 1.5, grid_size[1])[::-tick_n] + ) + ax.set_xticklabels(xtick[::tick_n]) + ax.set_yticklabels(xtick[::-tick_n]) + fig.set_size_inches(6, 6) + + fig.tight_layout() + + plt.xlim(peak_lowest_corner[0] - 10, peak_highest_corner[0] + 10) + plt.ylim(peak_lowest_corner[1], peak_highest_corner[1]) + plt.xlabel("v_x [km/s]") + plt.xticks(rotation=45) + plt.ylabel("v_y [km/s]") + plt.savefig(f"{name}_{slice_i}_{i}_{scaling}_{run_title}.png", dpi=300) + ax.clear() + fig.clear() + diff --git a/vlasovsolver/vlasovmover.cpp b/vlasovsolver/vlasovmover.cpp index 8ac4d5f607..15b3ecd54a 100644 --- a/vlasovsolver/vlasovmover.cpp +++ b/vlasovsolver/vlasovmover.cpp @@ -20,8 +20,11 @@ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */ +#include +#include #include #include +#include #include #include @@ -421,10 +424,12 @@ void calculateAcceleration(const uint popID,const uint globalMaxSubcycles,const // set seed, initialise generator and get value. The order is the same // for all cells, but varies with timestep. - std::default_random_engine rndState; - rndState.seed(P::tstep); + std::knuth_b rndState; + uint64_t seed=(P::seed*(uint64_t(P::tstep)+1)<<32)+uint64_t(P::mapOrderShift); + rndState.seed(seed); uint map_order = std::uniform_int_distribution<>(0,2)(rndState); + // assert((P::tstep != P::tstep+P::mapOrderShift)); // Calculate length of step for each cell #pragma omp parallel for for (size_t c=0; c