diff --git a/Jenkinsfile b/Jenkinsfile index a65e51f6..ea5d40aa 100644 --- a/Jenkinsfile +++ b/Jenkinsfile @@ -44,7 +44,7 @@ pipeline { stage('Fast Flavor'){ steps{ dir('Exec'){ - sh 'python ../Scripts/initial_conditions/st2_2beam_fast_flavor.py' + sh 'python ../Scripts/initial_conditions/st2_2beam_fast_flavor.py' sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_fast_flavor' sh 'python ../Scripts/tests/fast_flavor_test.py' sh 'rm -rf plt*' @@ -54,7 +54,7 @@ pipeline { stage('Fast Flavor k'){ steps{ dir('Exec'){ - sh 'python ../Scripts/initial_conditions/st3_2beam_fast_flavor_nonzerok.py' + sh 'python ../Scripts/initial_conditions/st3_2beam_fast_flavor_nonzerok.py' sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_fast_flavor_nonzerok' sh 'python ../Scripts/tests/fast_flavor_k_test.py' sh 'rm -rf plt*' @@ -80,26 +80,13 @@ pipeline { } } - stage('Fiducial 3F CPU HDF5'){ steps{ - dir('Exec'){ - sh 'cp ../makefiles/GNUmakefile_jenkins_HDF5 GNUmakefile' - sh 'make realclean; make generate; make -j' - sh 'python ../Scripts/initial_conditions/st4_linear_moment_ffi_3F.py' - sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.ex ../sample_inputs/inputs_1d_fiducial' - /*sh 'python3 ../Scripts/babysitting/avgfee_HDF5.py'*/ - sh 'rm -rf plt*' - } - } - } - stage('Collisions flavor instability'){ steps{ dir('Exec'){ - sh 'cp ../makefiles/GNUmakefile_jenkins_HDF5_CUDA GNUmakefile' - sh 'make realclean; make generate NUM_FLAVORS=2; make -j NUM_FLAVORS=2' sh 'python ../Scripts/initial_conditions/st8_coll_inst_test.py' sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_collisional_instability_test' sh 'python ../Scripts/data_reduction/reduce_data.py' sh 'python ../Scripts/tests/coll_inst_test.py' + archiveArtifacts artifacts: '*.pdf' sh 'rm -rf plt* *pdf' } } @@ -109,10 +96,24 @@ pipeline { sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_bc_periodic_init' sh 'python ../Scripts/collisions/writeparticleinfohdf5.py' sh 'python ../Scripts/tests/bc_empty_init_test.py' + archiveArtifacts artifacts: '*.pdf' sh 'rm -rf plt* *pdf' } } } + + stage('Fiducial 3F CPU HDF5'){ steps{ + dir('Exec'){ + sh 'cp ../makefiles/GNUmakefile_jenkins_HDF5 GNUmakefile' + sh 'make realclean; make generate; make -j' + sh 'python ../Scripts/initial_conditions/st4_linear_moment_ffi_3F.py' + sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.ex ../sample_inputs/inputs_1d_fiducial' + /*sh 'python3 ../Scripts/babysitting/avgfee_HDF5.py'*/ + sh 'rm -rf plt*' + } + } + } + stage('Collisions to equilibrium'){ steps{ dir('Exec'){ sh 'cp ../makefiles/GNUmakefile_jenkins_HDF5_CUDA GNUmakefile' @@ -120,7 +121,7 @@ pipeline { sh 'python ../Scripts/initial_conditions/st7_empty_particles.py' sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_coll_equi_test' sh 'python ../Scripts/tests/coll_equi_test.py' - sh 'rm -rf plt* *pdf' + sh 'rm -rf plt*' } } } @@ -128,7 +129,7 @@ pipeline { stage('Fermi-Dirac test'){ steps{ dir('Exec'){ sh 'python ../Scripts/initial_conditions/st9_empty_particles_multi_energy.py' - sh 'python ../Scripts/collisions/nsm_constant_background_rho_Ye_T__writer.py' + sh 'python ../Scripts/collisions/nsm_constant_background_rho_Ye_T_writer.py' sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_fermi_dirac_test' sh 'python ../Scripts/collisions/writeparticleinfohdf5.py' sh 'python ../Scripts/tests/fermi_dirac_test.py' diff --git a/Scripts/collisions/compute_dE_coll_inst_test.py b/Scripts/collisions/compute_dE_coll_inst_test.py index dec22d2c..4f4ae1c7 100644 --- a/Scripts/collisions/compute_dE_coll_inst_test.py +++ b/Scripts/collisions/compute_dE_coll_inst_test.py @@ -18,16 +18,19 @@ hc = h * c # erg cm # Simulation parameters -V = 1 # Volume of a cell ( ccm ) +V = (1.0e3)**3 # Volume of a cell ( ccm ) Ndir = 92 # Number of momentum beams isotropically distributed per cell E = 20.0 # Neutrinos and antineutrinos energy bin center ( Mev ) T = 7.0 # Background matter temperature ( Mev ) -N_eq_electron_neutrino = 3.260869565e+31 # Number of electron neutrinos at equilibrium +N_eq_electron_neutrino_density = 3.0e33 # Density of electron neutrinos at equilibrium (ccm) +N_eq_electron_neutrino = N_eq_electron_neutrino_density * V / Ndir # Number of electron neutrinos at equilibrium per beam +print(f'Number of electron neutrinos at equilibrium per beam = {N_eq_electron_neutrino}') u_electron_neutrino = 20.0 # Electron neutrino chemical potential ( Mev ) # Fermi-dirac distribution factor for electron neutrinos f_eq_electron_neutrinos = 1 / ( 1 + np.exp( ( E - u_electron_neutrino ) / T ) ) # adimentional +print(f'Fermi-Dirac distribution factor for electron neutrinos = {f_eq_electron_neutrinos}') # We know : # dE^3 = 3 * Neq * ( hc )^ 3 / ( dV * dOmega * feq ) @@ -47,9 +50,13 @@ dE=np.real(complex_deltaE) # Electron neutrino flavor -N_eq_electron_antineutrino = 2.717391304e+31 # Number of electron antineutrinos at equilibrium +N_eq_electron_antineutrino_density = 2.5e33 # Density of electron antineutrinos at equilibrium (ccm) +N_eq_electron_antineutrino = N_eq_electron_antineutrino_density * V / Ndir # Number of electron antineutrinos at equilibrium per beam +print(f'Number of electron antineutrinos at equilibrium per beam = {N_eq_electron_antineutrino}') +Vphase = V * ( 4 * np.pi / Ndir ) * delta_E_cubic / 3 +print(f'Phase space volume in ccm MeV^3 = {Vphase}') # Computing electron antineutrino chemical potential f_eq_electron_antineutrino = 3 * N_eq_electron_antineutrino * ( hc )**3 / ( V * ( 4 * np.pi / Ndir ) * delta_E_cubic ) u_electron_antineutrino = E - T * np.log( 1 / f_eq_electron_antineutrino - 1 ) -print(f'Electron neutrino chemical potential in MeV = {u_electron_antineutrino}') \ No newline at end of file +print(f'Electron antineutrino chemical potential in MeV = {u_electron_antineutrino}') \ No newline at end of file diff --git a/Scripts/collisions/convert_cellh5_inidat.py b/Scripts/collisions/convert_cellh5_inidat.py new file mode 100644 index 00000000..c0dc7695 --- /dev/null +++ b/Scripts/collisions/convert_cellh5_inidat.py @@ -0,0 +1,79 @@ +''' +Created by Erick Urquilla, Department of Physics and Astronomy, University of Tennessee, Knoxville. +This script convert the cell hdf5 files into a particle_input.dat file that can be used as an initial condition +for the Emu code. +The cell hdf5 files are generated by the script write_particles_per_cell.py +The particle data is stored in hdf5 files with the following labels: cell_i_j_k.h5 +where i, j, k are the cell indexes in the x, y, z directions respectively. +The ddf5 files has to be in the directory this script is run. +''' + +import glob +import h5py +import numpy as np +import sys +import os +importpath = os.path.dirname(os.path.realpath(__file__)) +sys.path.append(importpath) +sys.path.append(importpath+"/../initial_conditions") +from initial_condition_tools import uniform_sphere, write_particles +import amrex_plot_tools as amrex + +def read_hdf5_file(file_name): + """ + Reads an HDF5 file and returns its datasets as a dictionary. + + Parameters: + file_name (str): Path to the HDF5 file. + + Returns: + dict: A dictionary where keys are dataset names and values are NumPy arrays. + """ + with h5py.File(file_name, 'r') as file: + return {dataset: np.array(file[dataset]) for dataset in file.keys()} + +# Find all files matching the pattern +file_pattern = 'cell_*_*_*.h5' +files = glob.glob(file_pattern) + +data = read_hdf5_file(files[0]) + +labels_nf2=['pos_x','pos_y','pos_z', 'time', 'x', 'y', 'z', 'pupx', 'pupy', 'pupz', 'pupt', 'N00_Re', 'N01_Re', 'N01_Im', 'N11_Re', 'N00_Rebar', 'N01_Rebar', 'N01_Imbar', 'N11_Rebar', 'TrHN', 'Vphase'] +labels_nf3=['pos_x','pos_y','pos_z','time','x', 'y', 'z', 'pupx', 'pupy', 'pupz', 'pupt', 'N00_Re', 'N01_Re', 'N01_Im', 'N02_Re', 'N02_Im', 'N11_Re', 'N12_Re', 'N12_Im' ,'N22_Re', 'N00_Rebar', 'N01_Rebar', 'N01_Imbar', 'N02_Rebar', 'N02_Imbar', 'N11_Rebar', 'N12_Rebar' ,'N12_Imbar', 'N22_Rebar', 'TrHN', 'Vphase'] + +NF = 0 +if len(data.keys()) == len(labels_nf2): + labels = labels_nf2 + NF = 2 +elif len(data.keys()) == len(labels_nf3): + labels = labels_nf3 + NF = 3 +else: + print(f'len(data.keys()) = {len(data.keys())} does not match any of the known labels') + exit(1) + +# Get variable keys +rkey, ikey = amrex.get_particle_keys(NF, ignore_pos=True) + +# Generate the number of variables that describe each particle +n_variables = len(rkey) + +# Generate the number of particles +n_particles = len(data[labels[0]]) + +# Initialize a NumPy array to store all particles +particles = np.zeros((n_particles, n_variables)) + +# Volume of the simulation box +volume = 1e5**3 # ccm + +for label in rkey: + + print(f'Writing {label} to particles') + particles[:, rkey[label]] = data[label] + + if (label.startswith('N') | label.startswith('Vphase')): + particles[:, rkey[label]] /= volume + +# Write particles initial condition file +write_particles(np.array(particles), NF, "particle_input.dat") \ No newline at end of file diff --git a/Scripts/collisions/nsm_constant_background_rho_Ye_T__writer.py b/Scripts/collisions/nsm_constant_background_rho_Ye_T_writer.py similarity index 100% rename from Scripts/collisions/nsm_constant_background_rho_Ye_T__writer.py rename to Scripts/collisions/nsm_constant_background_rho_Ye_T_writer.py diff --git a/Scripts/collisions/single_particle_plot.py b/Scripts/collisions/single_particle_plot.py new file mode 100644 index 00000000..7aec75dd --- /dev/null +++ b/Scripts/collisions/single_particle_plot.py @@ -0,0 +1,100 @@ +''' +Created by Erick Urquilla. University of Tennessee Knoxville, USA. +This script is used to plot the single particle polarization vector in the x, y, and z directions. +This script works exclusively for the two flavor case. +The script reads the plt files generated with the writeparticleinfohdf5.py script. +''' + +import numpy as np +import glob +import h5py +import matplotlib.pyplot as plt + +file_names = np.array(glob.glob('plt*.h5')) +file_names = [file_name for file_name in file_names if '_reduced_data' not in file_name] +file_names.sort(key=lambda x: int(x.split('plt')[1].split('.h5')[0])) + +particle_momentum = [0,0,0,0] +with h5py.File(file_names[0], 'r') as hdf: + data = {} + for key in hdf.keys(): + data[key] = np.array(hdf.get(key)) + particle_momentum[0] = data['pupx'][0] + particle_momentum[1] = data['pupy'][0] + particle_momentum[2] = data['pupz'][0] + particle_momentum[3] = data['pupt'][0] + +time = np.zeros(len(file_names)) +N00_Re = np.zeros(len(file_names)) +N00_Rebar = np.zeros(len(file_names)) +N01_Im = np.zeros(len(file_names)) +N01_Imbar = np.zeros(len(file_names)) +N01_Re = np.zeros(len(file_names)) +N01_Rebar = np.zeros(len(file_names)) +N11_Re = np.zeros(len(file_names)) +N11_Rebar = np.zeros(len(file_names)) + +for i, file in enumerate(file_names): + with h5py.File(file, 'r') as hdf: + + # ['N00_Re', 'N00_Rebar', 'N01_Im', 'N01_Imbar', 'N01_Re', 'N01_Rebar', 'N11_Re', 'N11_Rebar', 'TrHN', 'Vphase', 'pos_x', 'pos_y', 'pos_z', 'pupt', 'pupx', 'pupy', 'pupz', 'time', 'x', 'y', 'z'] + + data = {} + for key in hdf.keys(): + data[key] = np.array(hdf.get(key)) + + for j in range(len(data['pupx'])): + if (data['pupx'][j] == particle_momentum[0] and + data['pupy'][j] == particle_momentum[1] and + data['pupz'][j] == particle_momentum[2] and + data['pupt'][j] == particle_momentum[3]): + time[i] = data['time'][j] + N00_Re[i] = data['N00_Re'][j] + N00_Rebar[i] = data['N00_Rebar'][j] + N01_Im[i] = data['N01_Im'][j] + N01_Imbar[i] = data['N01_Imbar'][j] + N01_Re[i] = data['N01_Re'][j] + N01_Rebar[i] = data['N01_Rebar'][j] + N11_Re[i] = data['N11_Re'][j] + N11_Rebar[i] = data['N11_Rebar'][j] + break + +P_x = 0.5 * N01_Re +P_y = 0.5 * N01_Im +P_z = 0.5 * (N00_Re - N11_Re) + +# Plot P_x vs P_y +plt.plot(P_x, P_y, color='gray', alpha=0.5) +sc = plt.scatter(P_x, P_y, c=time, cmap='viridis') +plt.colorbar(sc, label=r'$t \, (s)$') +plt.xlabel(r'$P_x$') +plt.ylabel(r'$P_y$') +plt.savefig('particle_momentum_plot_Px_Py.pdf') +plt.close() + +# Plot log(P_x) vs log(P_y) +plt.plot(np.log10(np.abs(P_x)), np.log10(np.abs(P_y)), color='gray', alpha=0.5) +sc = plt.scatter(np.log10(np.abs(P_x)), np.log10(np.abs(P_y)), c=time, cmap='viridis') +plt.colorbar(sc, label=r'$t \, (s)$') +plt.xlabel(r'$\log_{10}(|P_x|)$') +plt.ylabel(r'$\log_{10}(|P_y|)$') +plt.savefig('particle_momentum_plot_log_Px_Py.pdf') +plt.close() + +# Plot P_x vs P_z +plt.plot(P_x, P_z, color='gray', alpha=0.5) +sc = plt.scatter(P_x, P_z, c=time, cmap='viridis') +plt.colorbar(sc, label=r'$t \, (s)$') +plt.xlabel(r'$P_x$') +plt.ylabel(r'$P_z$') +plt.savefig('particle_momentum_plot_Px_Pz.pdf') +plt.close() + +# Plot P_y vs P_z +plt.plot(P_y, P_z, color='gray', alpha=0.5) +sc = plt.scatter(P_y, P_z, c=time, cmap='viridis') +plt.colorbar(sc, label=r'$t \, (s)$') +plt.xlabel(r'$P_y$') +plt.ylabel(r'$P_z$') +plt.savefig('particle_momentum_plot_Py_Pz.pdf') +plt.close() \ No newline at end of file diff --git a/Scripts/collisions/write_particles_per_cell.py b/Scripts/collisions/write_particles_per_cell.py new file mode 100755 index 00000000..2e8512d1 --- /dev/null +++ b/Scripts/collisions/write_particles_per_cell.py @@ -0,0 +1,175 @@ +''' +Created by Erick Urquilla, Department of Physics and Astronomy, University of Tennessee, Knoxville. +This script was based on the script reduced_data.py. It writes all the particle information +clasifyin it per cell. The data can be found in the directories plt*_particles +The particle data is stored in hdf5 files with the following labels: cell_i_j_k.h5 +where i, j, k are the cell indexes in the x, y, z directions respectively. +i=0, j=0, k=0 is the cell with the minimum x, y, z coordinates. + +Run this script with the command: +python write_particles_per_cell.py --fi -1 --bi 4 --bj 4 --bk 4 --nbi 96 --nbj 96 --nbk 64 +''' + +import os +os.environ['HDF5_USE_FILE_LOCKING'] = 'FALSE' +import sys +sys.path.append(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))+'/data_reduction') +import numpy as np +import glob +import multiprocessing as mp +import h5py +import amrex_plot_tools as amrex +import emu_yt_module as emu +from multiprocessing import Pool +import argparse +import time + +class GridData(object): + def __init__(self, ad): + x = ad['index','x'].d + y = ad['index','y'].d + z = ad['index','z'].d + dx = ad['index','dx'].d + dy = ad['index','dy'].d + dz = ad['index','dz'].d + self.ad = ad + self.dx = dx[0] + self.dy = dy[0] + self.dz = dz[0] + self.xmin = np.min(x-dx/2.) + self.ymin = np.min(y-dy/2.) + self.zmin = np.min(z-dz/2.) + self.xmax = np.max(x+dx/2.) + self.ymax = np.max(y+dy/2.) + self.zmax = np.max(z+dz/2.) + self.nx = int((self.xmax - self.xmin) / self.dx + 0.5) + self.ny = int((self.ymax - self.ymin) / self.dy + 0.5) + self.nz = int((self.zmax - self.zmin) / self.dz + 0.5) + print(self.nx, self.ny, self.nz) + + def get_particle_cell_ids(self,rdata): + # get coordinates + x = rdata[:,rkey["x"]] + y = rdata[:,rkey["y"]] + z = rdata[:,rkey["z"]] + ix = (x/self.dx).astype(int) + iy = (y/self.dy).astype(int) + iz = (z/self.dz).astype(int) + + # HACK - get this grid's bounds using particle locations + ix -= np.min(ix) + iy -= np.min(iy) + iz -= np.min(iz) + nx = np.max(ix)+1 + ny = np.max(iy)+1 + nz = np.max(iz)+1 + idlist = (iz + nz*iy + nz*ny*ix).astype(int) + + return idlist + +def writehdf5files(dire): + + if not os.path.exists(dire + "_particles"): + os.makedirs(dire + "_particles") + + eds = emu.EmuDataset(dire) + t = eds.ds.current_time + ad = eds.ds.all_data() + + header = amrex.AMReXParticleHeader(dire+"/neutrinos/Header") + grid_data = GridData(ad) + nlevels = len(header.grids) + assert nlevels==1 + level = 0 + + if (n_box_x < 0) and (n_box_y < 0) and (n_box_z < 0): + ngrids = len(header.grids[level]) + print(f"Number of grids in level {level}: {ngrids}") + boxes = np.arange(len(header.grids[level])) + print(f"Boxes to process: {boxes}") + else: + boxes = [box_i + n_box_x * box_j + n_box_x * n_box_y * box_k] + + print(f"Boxes to process: {boxes}") + + # loop over all boxes + for gridID in boxes: + + print(f"Processing grid {gridID} in directory {dire}") + + # read particle data on a single grid + idata, rdata = amrex.read_particle_data(dire, ptype="neutrinos", level_gridID=(level,gridID)) + + # get cell index for each particle + cell_x_index = (rdata[:,rkey['pos_x']] / grid_data.dx).astype(int) + cell_y_index = (rdata[:,rkey['pos_y']] / grid_data.dy).astype(int) + cell_z_index = (rdata[:,rkey['pos_z']] / grid_data.dz).astype(int) + cell_index = np.stack((cell_x_index, cell_y_index, cell_z_index), axis=1) + # get unique cell index + unique_cell_index = np.unique(cell_index, axis=0) + + for i in range(unique_cell_index.shape[0]): + cell_filename = f"{dire}_particles/cell_{unique_cell_index[i,0]}_{unique_cell_index[i,1]}_{unique_cell_index[i,2]}.h5" + if os.path.exists(cell_filename): + print(f"file {cell_filename} already exists") + else: + mask_cell_group = np.all(cell_index == unique_cell_index[i], axis=1) + cell_group = rdata[mask_cell_group] + with h5py.File(cell_filename, 'w') as cell_hf: + for label in labels: + cell_hf.create_dataset(label, data=cell_group[:, rkey[label]], maxshape=(None,), chunks=True) + +directories = sorted(glob.glob("plt*/neutrinos")) +directories = [directories[i].split('/')[0] for i in range(len(directories))] # remove "neutrinos" +directories.sort(key=lambda x: int(x.split('plt')[-1])) + +if not directories: + print("No new directories to process.") + sys.exit(0) + +parser = argparse.ArgumentParser(description='Process some directories.') +parser.add_argument('--fi', type=int, default=-2, help='index of the plt file to be analyzed') +parser.add_argument('--bi', type=int, default=-2, help='chosen box index in the x direction') +parser.add_argument('--bj', type=int, default=-2, help='chosen box index in the y direction') +parser.add_argument('--bk', type=int, default=-2, help='chosen box index in the z direction') +parser.add_argument('--nbi', type=int, default=-2, help='number of boxes in the x direction') +parser.add_argument('--nbj', type=int, default=-2, help='number of boxes in the y direction') +parser.add_argument('--nbk', type=int, default=-2, help='number of boxes in the z direction') +args = parser.parse_args() + +if args.fi >= -1: + if args.fi == -1: + directories = [directories[-1]] + elif args.fi < len(directories): + directories = [directories[args.fi]] + else: + print(f"Index {args.fi} is out of range. Exiting.") + sys.exit(1) + +box_i = args.bi +box_j = args.bj +box_k = args.bk +print(f"Box indices: i={box_i}, j={box_j}, k={box_k}") +n_box_x = args.nbi +n_box_y = args.nbj +n_box_z = args.nbk +print(f"Number of boxes: nbi={n_box_x}, nbj={n_box_y}, nbk={n_box_z}") + +print(f"Directories to process: {directories}") + +# get NF +eds = emu.EmuDataset(directories[0]) +NF = eds.get_num_flavors() +if NF==2: + rkey, ikey = amrex.get_particle_keys(NF) + labels=['pos_x','pos_y','pos_z', 'time', 'x', 'y', 'z', 'pupx', 'pupy', 'pupz', 'pupt', 'N00_Re', 'N01_Re', 'N01_Im', 'N11_Re', 'N00_Rebar', 'N01_Rebar', 'N01_Imbar', 'N11_Rebar', 'TrHN', 'Vphase'] +if NF==3: + rkey, ikey = amrex.get_particle_keys(NF) + labels=['pos_x','pos_y','pos_z','time','x', 'y', 'z', 'pupx', 'pupy', 'pupz', 'pupt', 'N00_Re', 'N01_Re', 'N01_Im', 'N02_Re', 'N02_Im', 'N11_Re', 'N12_Re', 'N12_Im' ,'N22_Re', 'N00_Rebar', 'N01_Rebar', 'N01_Imbar', 'N02_Rebar', 'N02_Imbar', 'N11_Rebar', 'N12_Rebar' ,'N12_Imbar', 'N22_Rebar', 'TrHN', 'Vphase'] + +start_time_gpu = time.time() +nproc = mp.cpu_count() +pool = Pool(nproc) +pool.map(writehdf5files, directories) +end_time_gpu = time.time() +print(f"Processing time: {end_time_gpu - start_time_gpu} seconds") \ No newline at end of file diff --git a/Scripts/data_reduction/convertToHDF5.py b/Scripts/data_reduction/convertToHDF5.py index 6489aaa8..1226b3a4 100755 --- a/Scripts/data_reduction/convertToHDF5.py +++ b/Scripts/data_reduction/convertToHDF5.py @@ -28,13 +28,26 @@ # misc # ######## fluid_vars = ["N","Fx","Fy","Fz"] +# fluid_vars = ["N","Fx","Fy","Fz","Pxx","Pxy","Pxz","Pyy","Pyz","Pzz"] # Use this line to write the pressure moment in the H5 file. nunubar = ["","bar"] def convert_to_HDF5(sim_directory, DELETE_ALL_BUT_LAST_RESTART=False): ######################### # loop over directories # ######################### - fluid_directories = sorted(glob.glob(sim_directory+"/plt?????")) + + # fluid_directories = sorted(glob.glob(sim_directory+"/plt*/neutrinos")) + # fluid_directories = [d[:-10] if d.endswith('/neutrinos') else d for d in fluid_directories] + fluid_directories = sorted(set(glob.glob(sim_directory+"/plt*"))) + + def get_plt_number(x): + try: + return int(x.split('plt')[-1]) + except ValueError: + return None + fluid_directories = [d for d in fluid_directories if get_plt_number(d) is not None] + # Remove duplicates and sort by plt number + fluid_directories = sorted(set(fluid_directories), key=lambda x: int(x.split('plt')[-1])) nfluid = len(fluid_directories) print(fluid_directories) @@ -49,14 +62,19 @@ def convert_to_HDF5(sim_directory, DELETE_ALL_BUT_LAST_RESTART=False): if d==fluid_directories[0]: NF = eds.get_num_flavors() allData = h5py.File(sim_directory+"/allData.h5","w") + allData["dx(cm)"] = eds.dx + allData["dy(cm)"] = eds.dy allData["dz(cm)"] = eds.dz + allData["Nx"] = eds.Nx + allData["Ny"] = eds.Ny + allData["Nz"] = eds.Nz allData.create_dataset("t(s)", data=np.zeros(0), maxshape=(None,), dtype=datatype) allData.create_dataset("it", data=np.zeros(0), maxshape=(None,), dtype=int) # create fluid data sets - maxshape = (None, eds.Nz) - chunkshape = (1, eds.Nz) - zeros = np.zeros((0,eds.Nz)) + maxshape = (None,eds.Nx,eds.Ny,eds.Nz) + chunkshape = (1,eds.Nx,eds.Ny,eds.Nz) + zeros = np.zeros((0,eds.Nx,eds.Ny,eds.Nz)) varlist = [] for v in fluid_vars: for f1 in range(NF): diff --git a/Scripts/data_reduction/convertToHDF5_split.py b/Scripts/data_reduction/convertToHDF5_split.py new file mode 100755 index 00000000..82ffebab --- /dev/null +++ b/Scripts/data_reduction/convertToHDF5_split.py @@ -0,0 +1,101 @@ +# used to make plots but now just generates a hdf5 file with domain-averaged data. +# Run in the directory of the simulation the data should be generated for. +# Still has functionality for per-snapshot plots, but the line is commented out. +# This version averages the magnitudes of off-diagonal components rather than the real/imaginary parts +# also normalizes fluxes by sumtrace of N rather than F. +# This data is used for the growth plot. +# Note - also tried a version using maxima rather than averages, and it did not make the growth plot look any better. + +import os +os.environ['HDF5_USE_FILE_LOCKING'] = 'FALSE' +import sys +sys.path.append(os.path.dirname(os.path.abspath(__file__))) +import numpy as np +import matplotlib.pyplot as plt +import yt +import glob +import multiprocessing as mp +import h5py +import amrex_plot_tools as amrex +import emu_yt_module as emu +from multiprocessing import Pool +import scipy.special +import shutil + +# NOTE - assumes particle output is done at a multiple of fluid output + +######## +# misc # +######## +fluid_vars = ["N","Fx","Fy","Fz"] +# fluid_vars = ["N","Fx","Fy","Fz","Pxx","Pxy","Pxz","Pyy","Pyz","Pzz"] # Use this line to write the pressure moment in the H5 file. +nunubar = ["","bar"] + +def convert_to_HDF5(sim_directory, DELETE_ALL_BUT_LAST_RESTART=False): + ######################### + # loop over directories # + ######################### + + # fluid_directories = sorted(glob.glob(sim_directory+"/plt*/neutrinos")) + # fluid_directories = [d[:-10] if d.endswith('/neutrinos') else d for d in fluid_directories] + fluid_directories = sorted(set(glob.glob(sim_directory+"/plt*"))) + + def get_plt_number(x): + try: + return int(x.split('plt')[-1]) + except ValueError: + return None + fluid_directories = [d for d in fluid_directories if get_plt_number(d) is not None] + # Remove duplicates and sort by plt number + fluid_directories = sorted(set(fluid_directories), key=lambda x: int(x.split('plt')[-1])) + nfluid = len(fluid_directories) + + print(fluid_directories) + + for d in fluid_directories: + print(d) + eds = emu.EmuDataset(d) + t = eds.ds.current_time + ad = eds.ds.all_data() + datatype = ad['boxlib',"N00_Re"].d.dtype + + out_file = sim_directory + "/mesh_" + d.split('/')[-1] + ".h5" + if os.path.exists(out_file): + print(f"Skipping existing file: {out_file}") + continue + + NF = eds.get_num_flavors() + allData = h5py.File(out_file,"w") + allData["dx(cm)"] = eds.dx + allData["dy(cm)"] = eds.dy + allData["dz(cm)"] = eds.dz + allData["Nx"] = eds.Nx + allData["Ny"] = eds.Ny + allData["Nz"] = eds.Nz + allData["t(s)"] = eds.ds.current_time + allData["it"] = int(d.split('plt')[-1]) + + varlist = [] + for v in fluid_vars: + for f1 in range(NF): + for f2 in range(f1,NF): + for b in nunubar: + varlist.append(v+str(f1)+str(f2)+"_Re"+b+"(1|ccm)") + if f2!=f1: + varlist.append(v+str(f1)+str(f2)+"_Im"+b+"(1|ccm)") + + for v in varlist: + data = eds.cg[v[:-7]] + allData[v] = data + + allData.close() + + if DELETE_ALL_BUT_LAST_RESTART: + particle_directories = [d[:-10] for d in sorted(glob.glob(sim_directory+"/plt*/neutrinos"))] + last_particle_directory = particle_directories[-1] + for d in fluid_directories: + if d != last_particle_directory: + shutil.rmtree(d) + +if __name__ == "__main__": + convert_to_HDF5(".") diff --git a/Scripts/initial_conditions/st10_coll_inst_genr_cond.py b/Scripts/initial_conditions/st10_coll_inst_genr_cond.py new file mode 100644 index 00000000..4e051fa6 --- /dev/null +++ b/Scripts/initial_conditions/st10_coll_inst_genr_cond.py @@ -0,0 +1,85 @@ +''' +Created by Erick Urquilla, Department of Physics and Astronomy, University of Tennessee, Knoxville. +This script is used to create the particle initial conditions that are isotropic in momentum space. +However it use energy dependent number densities as initial conditions. +This script will be used to generate the initial conditions for the CFI simulations in the gw170817 paper. +''' +import numpy as np +import sys +import os +importpath = os.path.dirname(os.path.realpath(__file__)) +sys.path.append(importpath) +sys.path.append(importpath+"/../data_reduction") +from initial_condition_tools import uniform_sphere, moment_interpolate_particles, minerbo_interpolate, write_particles +import amrex_plot_tools as amrex + +# These initial conditions are intended to replicate the collisional instability outputs in +# "Collisional Flavor Instabilities of Supernova Neutrinos" by L. Johns [2104.11369]. + +NF = 3 # Number of flavors +nphi_equator = 16 # number of direction in equator ---> theta = pi/2 + +nnue = [9.361931627879481e+28, 1.966560674916465e+30, 1.2478046652293569e+31, 4.213408072316177e+31, 8.668898381315493e+31, 1.2345408071889402e+32, 1.2650011635960715e+32, 9.075213507474793e+31, 4.363393832402667e+31, 1.3323289626121315e+31, 2.4272162118330767e+30, 2.4262494851051388e+29, 1.1746613663424375e+28] # 1/ccm +nnua = [1.4021316040267574e+28, 2.541581132928583e+29, 3.117477188768948e+30, 1.3580108877004713e+31, 3.4640494620358637e+31, 5.957577726384348e+31, 7.223385086336037e+31, 6.179574473449261e+31, 3.619607827490873e+31, 1.370308674959416e+31, 3.0878517493962613e+30, 3.72763909837056e+29, 2.117667306969499e+28] # 1/ccm +nnux = [3.362254311220344e+27, 3.328228400009195e+28, 1.3450298899247732e+29, 3.439338355099913e+29, 6.35649925739673e+29, 8.842179118437819e+29, 9.304787828388212e+29, 7.182754529428967e+29, 4.007279212720665e+29, 1.5434940087807304e+29, 3.795479165923593e+28, 5.714921116371749e+27, 4.824408960556585e+26] # 1/ccm + +# Neutrino flux factors +fnue = np.array([0.0 , 0.0 , 0.0]) +fnua = np.array([0.0 , 0.0 , 0.0]) +fnux = np.array([0.0 , 0.0 , 0.0]) + +# Energy bins from NuLib table (Cutted) +# Energy bin centers extracted from NuLib table +energies_center_Mev = [1, 3, 5.2382, 8.0097, 11.442, 15.691, 20.953, 27.468, 35.536, 45.525, 57.895, 73.212, 92.178] # Energy in Mev +# Energy bin bottom extracted from NuLib table +energies_bottom_Mev = [0, 2, 4, 6.4765, 9.543, 13.34, 18.042, 23.864, 31.073, 39.999, 51.052, 64.738, 81.685] +# Energy bin top extracted from NuLib table +energies_top_Mev = [2, 4, 6.4765, 9.543, 13.34, 18.042, 23.864, 31.073, 39.999, 51.052, 64.738, 81.685, 102.67] + +# Energies in ergs +energies_center_erg = np.array(energies_center_Mev) * 1e6*amrex.eV # Energy in ergs +energies_bottom_erg = np.array(energies_bottom_Mev) * 1e6*amrex.eV # Energy in ergs +energies_top_erg = np.array(energies_top_Mev ) * 1e6*amrex.eV # Energy in ergs + +# Get variable keys +rkey, ikey = amrex.get_particle_keys(NF, ignore_pos=True) + +# Generate the number of variables that describe each particle +n_variables = len(rkey) + +# Get the momentum distribution of the particles +phat = uniform_sphere(nphi_equator) + +particles = np.zeros((len(energies_center_erg), len(phat), n_variables)) # [particle, direction, variable] + +for i in range(len(energies_center_erg)): + + # Matrix to save the neutrino number densities + nnu = np.zeros((2,NF)) + nnu[0,0] = nnue[i] + nnu[1,0] = nnua[i] + nnu[:,1:] = nnux[i] + + # Matrix to save the neutrino number densities fluxes + fnu = np.zeros((2,NF,3)) + fnu[0,0,:] = nnu[0,0] * fnue + fnu[1,0,:] = nnu[1,0] * fnua + fnu[:,1:,:] = nnu[:,1:,np.newaxis] * fnux[np.newaxis,np.newaxis,:] + + # Generate particles + particles[i] = moment_interpolate_particles(nphi_equator, nnu, fnu, energies_center_erg[i], uniform_sphere, minerbo_interpolate) # [particle, variable] + + # Generate the number of directions + n_directions = len(particles[i]) + + # Compute the phase space volume dOmega * dE^3 / 3 + Vphase = ( 4.0 * np.pi / n_directions ) * ( ( energies_top_erg[i] ** 3 - energies_bottom_erg[i] ** 3 ) / 3.0 ) + + # Save V_phase in the last column of the particle array + particles[i,:,-1] = Vphase + +# Reshape the particles array +particles = particles.reshape(len(energies_center_erg) * len(phat), n_variables) + +# Write particles initial condition file +write_particles(np.array(particles), NF, "particle_input.dat") \ No newline at end of file diff --git a/Scripts/initial_conditions/st9_empty_particles_multi_energy.py b/Scripts/initial_conditions/st9_empty_particles_multi_energy.py index 4ec43d26..a2e76e8d 100644 --- a/Scripts/initial_conditions/st9_empty_particles_multi_energy.py +++ b/Scripts/initial_conditions/st9_empty_particles_multi_energy.py @@ -16,12 +16,28 @@ nphi_equator = 16 # number of direction in equator NF = 3 # number of flavors +''' +# Energy bins from NuLib table (No cutted) +--------- NuLib table energy bins --------- # Energy bin centers extracted from NuLib table energies_center_Mev = [1, 3, 5.23824, 8.00974, 11.4415, 15.6909, 20.9527, 27.4681, 35.5357, 45.5254, 57.8951, 73.2117, 92.1775, 115.662, 144.741, 180.748, 225.334, 280.542] # Energy in Mev # Energy bin bottom extracted from NuLib table energies_bottom_Mev = [0, 2, 4, 6.47649, 9.54299, 13.3401, 18.0418, 23.8636, 31.0725, 39.9989, 51.0519, 64.7382, 81.6853, 102.67, 128.654, 160.828, 200.668, 250] # Energy bin top extracted from NuLib table energies_top_Mev = [2, 4, 6.47649, 9.54299, 13.3401, 18.0418, 23.8636, 31.0725, 39.9989, 51.0519, 64.7382, 81.6853, 102.67, 128.654, 160.828, 200.668, 250, 311.085] +''' + +# The following energy bins are an extraction of the NuLib table +# I have delete the first five energy bin to make the fermi-dirac test facter +# If I simulation want to run with all energy bins the user must change the energy bins manually + +# Energy bins from NuLib table (Cutted) +# Energy bin centers extracted from NuLib table +energies_center_Mev = [15.6909, 20.9527, 27.4681, 35.5357, 45.5254, 57.8951, 73.2117, 92.1775, 115.662, 144.741, 180.748, 225.334, 280.542] # Energy in Mev +# Energy bin bottom extracted from NuLib table +energies_bottom_Mev = [13.3401, 18.0418, 23.8636, 31.0725, 39.9989, 51.0519, 64.7382, 81.6853, 102.67, 128.654, 160.828, 200.668, 250] +# Energy bin top extracted from NuLib table +energies_top_Mev = [18.0418, 23.8636, 31.0725, 39.9989, 51.0519, 64.7382, 81.6853, 102.67, 128.654, 160.828, 200.668, 250, 311.085] # Energies in ergs energies_center_erg = np.array(energies_center_Mev) * 1e6*amrex.eV # Energy in ergs diff --git a/Scripts/symbolic_hermitians/generate_code.py b/Scripts/symbolic_hermitians/generate_code.py index 2e223e95..6153643b 100755 --- a/Scripts/symbolic_hermitians/generate_code.py +++ b/Scripts/symbolic_hermitians/generate_code.py @@ -325,6 +325,8 @@ def delete_generated_files(): sgn = 1 line = "Real "+Vlist[icomp]+" = "+str(sgn)+"*("+M2list[icomp] + ")*PhysConst::c4/(2.*p.rdata(PIdx::pupt));" code.append(line) + line = Vlist[icomp]+" *= parms->attenuation_vacuum;" + code.append(line) write_code(code, os.path.join(args.emu_home,"Source/generated_files","Evolve.cpp_Vvac_fill")) #============================# @@ -404,11 +406,11 @@ def sgn(t,var): line = line + ");" code.append(line) code.append("") - line = "inside_parentheses = SI_partial + SI_partialbar" + line = "inside_parentheses = parms->attenuation_self_interaction * ( SI_partial + SI_partialbar )" # matter potential if("V00" in Vlist[icomp]): - line = line + " + " + rhoye + line = line + " + parms->attenuation_matter * " + rhoye line = line + ";" code.append(line) diff --git a/Scripts/tests/coll_inst_test.py b/Scripts/tests/coll_inst_test.py index 4d22f404..a99ea222 100644 --- a/Scripts/tests/coll_inst_test.py +++ b/Scripts/tests/coll_inst_test.py @@ -7,12 +7,10 @@ import numpy as np import argparse import glob -import EmuReader import sys import os importpath = os.path.dirname(os.path.realpath(__file__))+"/../data_reduction/" sys.path.append(importpath) -import amrex_plot_tools as amrex import numpy as np import h5py import glob @@ -21,6 +19,63 @@ parser.add_argument("-na", "--no_assert", action="store_true", help="If --no_assert is supplied, do not raise assertion errors if the test error > tolerance.") args = parser.parse_args() +############################################################ +############################################################ +# PLOT SETTINGS +import matplotlib as mpl +import matplotlib.pyplot as plt +from matplotlib.ticker import AutoLocator, AutoMinorLocator, LogLocator +import os +import glob + +# Font settings +mpl.rcParams['font.size'] = 22 +mpl.rcParams['font.family'] = 'serif' +mpl.rc('text', usetex=False) + +# Tick settings +mpl.rcParams['xtick.major.width'] = 2 +mpl.rcParams['xtick.major.pad'] = 8 +mpl.rcParams['xtick.minor.size'] = 4 + +mpl.rcParams['xtick.minor.width'] = 2 +mpl.rcParams['ytick.major.size'] = 7 +mpl.rcParams['ytick.major.width'] = 2 +mpl.rcParams['ytick.minor.size'] = 4 +mpl.rcParams['ytick.minor.width'] = 2 + +# Axis linewidth +mpl.rcParams['axes.linewidth'] = 2 + +# Tick direction and enabling ticks on all sides +mpl.rcParams['xtick.direction'] = 'in' +mpl.rcParams['ytick.direction'] = 'in' +mpl.rcParams['xtick.top'] = True +mpl.rcParams['ytick.right'] = True + +# Function to apply custom tick locators and other settings to an Axes object +def apply_custom_settings(ax, leg, log_scale_y=False): + + if log_scale_y: + # Use LogLocator for the y-axis if it's in log scale + ax.set_yscale('log') + ax.yaxis.set_major_locator(LogLocator(base=10.0)) + ax.yaxis.set_minor_locator(LogLocator(base=10.0, subs='auto', numticks=100)) + else: + # Use AutoLocator for regular scales + ax.yaxis.set_major_locator(AutoLocator()) + ax.yaxis.set_minor_locator(AutoMinorLocator()) + + # Apply the AutoLocator for the x-axis + ax.xaxis.set_major_locator(AutoLocator()) + ax.xaxis.set_minor_locator(AutoMinorLocator()) + + # Legend settings + leg.get_frame().set_edgecolor('w') + leg.get_frame().set_linewidth(0.0) +############################################################ +############################################################ + if __name__ == "__main__": # Create a list of data files to read @@ -43,8 +98,8 @@ t[i] = np.array(hf['t(s)'][:][0]) # Fit the exponential function ( y = a e ^ ( b x ) ) to the data - l1 = 10 # initial item for fit - l2 = 40 # last item for fit + l1 = 30 # initial item for fit + l2 = 150 # last item for fit coefficients = np.polyfit(t[l1:l2], np.log(N_avg_mag[:,0,1][l1:l2]), 1) coefficients_bar = np.polyfit(t[l1:l2], np.log(Nbar_avg_mag[:,0,1][l1:l2]), 1) a = np.exp(coefficients[1]) @@ -54,50 +109,52 @@ print(f'{b} ---> EMU : Im Omega') print(f'{bbar} ---> EMU : Im Omegabar') - # Plots - import matplotlib.pyplot as plt - - # Plots N and Nbar - plt.plot(t, N_avg_mag[:,0,0], label = r'$N_{ee}$') - plt.plot(t, N_avg_mag[:,0,1], label = r'$N_{eu}$') - plt.plot(t, N_avg_mag[:,1,1], label = r'$N_{uu}$') - plt.plot(t, Nbar_avg_mag[:,0,0], label = r'$\bar{N}_{ee}$') - plt.plot(t, Nbar_avg_mag[:,0,1], label = r'$\bar{N}_{eu}$') - plt.plot(t, Nbar_avg_mag[:,1,1], label = r'$\bar{N}_{uu}$') - plt.legend() - plt.xlabel(r'$t$ (s)') - plt.ylabel(r'$N$ and $\bar{N}$') - plt.savefig('N_and_Nbar.pdf') + # Plots + fig, ax = plt.subplots() + ax.plot(t, N_avg_mag[:,0,0], label = r'$n_{ee}$') + ax.plot(t, N_avg_mag[:,0,1], label = r'$n_{eu}$') + ax.plot(t, N_avg_mag[:,1,1], label = r'$n_{uu}$') + ax.plot(t, Nbar_avg_mag[:,0,0], label = r'$\bar{n}_{ee}$') + ax.plot(t, Nbar_avg_mag[:,0,1], label = r'$\bar{n}_{eu}$') + ax.plot(t, Nbar_avg_mag[:,1,1], label = r'$\bar{n}_{uu}$') + ax.set_xlabel(r'$t$ (s)') + ax.set_ylabel(r'$n$ and $\bar{n}$ (cm$^{-3}$)') + leg = ax.legend(framealpha=0.0, ncol=1, fontsize=15) + apply_custom_settings(ax, leg, log_scale_y=False) + plt.tight_layout() + fig.savefig('N_and_Nbar.pdf', bbox_inches='tight') plt.clf() - # Plots N and F - plt.plot(t, N_avg_mag[:,0,0], label = r'$N_{ee}$') - plt.plot(t, N_avg_mag[:,0,1], label = r'$N_{eu}$') - plt.plot(t, N_avg_mag[:,1,1], label = r'$N_{uu}$') - plt.plot(t[l1:l2], N_avg_mag[:,0,1][l1:l2], label = f'Im Omega = {b}') - plt.plot(t, F_avg_mag[:,0,0,1], label = r'$F^x_{eu}$') - plt.plot(t, F_avg_mag[:,1,0,1], label = r'$F^y_{eu}$') - plt.plot(t, F_avg_mag[:,2,0,1], label = r'$F^z_{eu}$') - plt.legend() - plt.xlabel(r'$t$ (s)') - plt.ylabel(r'$N$ and $\vec{F}$') - plt.yscale('log') - plt.savefig('N_and_F.pdf') + fig, ax = plt.subplots() + ax.plot(t, N_avg_mag[:,0,0], label = r'$n_{ee}$') + ax.plot(t, N_avg_mag[:,0,1], label = r'$n_{eu}$') + ax.plot(t, N_avg_mag[:,1,1], label = r'$n_{uu}$') + ax.plot(t[l1:l2], N_avg_mag[:,0,1][l1:l2], label = f'Im Omega = {b:.2e} s$^{{-1}}$') + ax.plot(t, F_avg_mag[:,0,0,1], label = r'$f^x_{eu}$') + ax.plot(t, F_avg_mag[:,1,0,1], label = r'$f^y_{eu}$') + ax.plot(t, F_avg_mag[:,2,0,1], label = r'$f^z_{eu}$') + ax.set_xlabel(r'$t$ (s)') + ax.set_ylabel(r'$n$ and $\vec{f}$ (cm$^{-3}$)') + leg = ax.legend(framealpha=0.0, ncol=1, fontsize=12) + apply_custom_settings(ax, leg, log_scale_y=True) + plt.tight_layout() + fig.savefig('N_and_F.pdf', bbox_inches='tight') plt.clf() - # Plots Nbar and Fbar - plt.plot(t, Nbar_avg_mag[:,0,0], label = r'$\bar{N}_{ee}$') - plt.plot(t, Nbar_avg_mag[:,0,1], label = r'$\bar{N}_{eu}$') - plt.plot(t, Nbar_avg_mag[:,1,1], label = r'$\bar{N}_{uu}$') - plt.plot(t[l1:l2], Nbar_avg_mag[:,0,1][l1:l2], label = f'Im Omega = {bbar}') - plt.plot(t, Fbar_avg_mag[:,0,0,1], label = r'$\bar{F}^x_{eu}$') - plt.plot(t, Fbar_avg_mag[:,1,0,1], label = r'$\bar{F}^y_{eu}$') - plt.plot(t, Fbar_avg_mag[:,2,0,1], label = r'$\bar{F}^z_{eu}$') - plt.legend() - plt.xlabel(r'$t$ (s)') - plt.ylabel(r'$\bar{N}$ and $\vec{\bar{F}}$') - plt.yscale('log') - plt.savefig('Nbar_and_Fbar.pdf') + fig, ax = plt.subplots() + ax.plot(t, Nbar_avg_mag[:,0,0], label = r'$\bar{n}_{ee}$') + ax.plot(t, Nbar_avg_mag[:,0,1], label = r'$\bar{n}_{eu}$') + ax.plot(t, Nbar_avg_mag[:,1,1], label = r'$\bar{n}_{uu}$') + ax.plot(t[l1:l2], Nbar_avg_mag[:,0,1][l1:l2], label = f'Im Omega = {bbar:.2e} s$^{{-1}}$') + ax.plot(t, Fbar_avg_mag[:,0,0,1], label = r'$\bar{f}^x_{eu}$') + ax.plot(t, Fbar_avg_mag[:,1,0,1], label = r'$\bar{f}^y_{eu}$') + ax.plot(t, Fbar_avg_mag[:,2,0,1], label = r'$\bar{f}^z_{eu}$') + ax.set_xlabel(r'$t$ (s)') + ax.set_ylabel(r'$\bar{n}$ and $\vec{\bar{f}}$ (cm$^{-3}$)') + leg = ax.legend(framealpha=0.0, ncol=1, fontsize=12) + apply_custom_settings(ax, leg, log_scale_y=True) + plt.tight_layout() + fig.savefig('Nbar_and_Fbar.pdf', bbox_inches='tight') plt.clf() ###################################################################################### @@ -150,12 +207,6 @@ def myassert(condition): rel_error_bar = np.abs( bbar - b_lsa ) / np.abs( ( bbar + b_lsa ) / 2 ) rel_error_max = 0.05 - print(f"{rel_error} ---> relative error in ImOmega : EMU") - print(f"{rel_error_bar} ---> relative error in ImOmegabar : EMU") - - myassert( rel_error < rel_error_max ) - myassert( rel_error_bar < rel_error_max ) - ###################################################################################### ###################################################################################### @@ -262,15 +313,25 @@ def QKE(t,y): print(f"{rel_error_j} ---> relative erroror in Julien script") # Plotting Julien, EMU and Lucas LSA data - plt.plot(t, N_avg_mag[:,0,1], label = r'$N_{eu}$ EMU') - plt.plot(t[l1:l2], N_avg_mag[:,0,1][l1:l2], label = f'Im Omega EMU = {b}', linestyle = 'dashed') - plt.plot(time,N_eu_julien, label = r'$N_{eu}$ Julien script') - plt.plot(time[p1:p2],N_eu_julien[p1:p2], label = f'Im Omega Julien = {bj}', linestyle = 'dashed') - plt.plot(time[p1:p2], 1e23*np.exp(ImOmega_Lucas_LSA*time[p1:p2]), label = f'Im Omega Lucas LSA = {ImOmega_Lucas_LSA}', linestyle = 'dashed') - plt.xlabel(r'$t \ (\mathrm{s})$') - plt.ylabel(r'$N_{eu}$') - plt.title(f"Collisional flavor instability test") - plt.yscale('log') - plt.legend() - plt.savefig('EMU_Julien_LucasLSA_Neu.pdf') + fig, ax = plt.subplots() + ax.plot(t, N_avg_mag[:,0,1], label = r'$n_{eu}$ EMU') + ax.plot(t[l1:l2], N_avg_mag[:,0,1][l1:l2], label = f'Im$(\Omega_{{eu}})$ EMU : {b:.2e} s$^{{-1}}$', linestyle = 'dashed') + ax.plot(time, N_eu_julien, label = r'$n_{eu}$ Julien script') + ax.plot(time[p1:p2], N_eu_julien[p1:p2], label = f'Im$(\Omega_{{eu}})$ Julien F. : {bj:.2e} s$^{{-1}}$', linestyle = 'dashed') + ax.plot(time[p1:p2], 1e23*np.exp(ImOmega_Lucas_LSA*time[p1:p2]), label = f"Im$(\Omega_{{eu}})$ Luke J. (LSA) : {ImOmega_Lucas_LSA:.2e} s$^{{-1}}$", linestyle = 'dashed') + ax.set_xlabel(r'$t \ (\mathrm{s})$') + ax.set_ylabel(r'$n_{eu}$ (cm$^{-3}$)') + ax.set_title(f"Collisional flavor instability test") + leg = ax.legend(framealpha=0.0, ncol=1, fontsize=12) + apply_custom_settings(ax, leg, log_scale_y=True) + plt.tight_layout() + fig.savefig('EMU_Julien_LucasLSA_Neu.pdf', bbox_inches='tight') plt.close() + + # Asserts + + print(f"{rel_error} ---> relative error in ImOmega : EMU") + print(f"{rel_error_bar} ---> relative error in ImOmegabar : EMU") + + myassert( rel_error < rel_error_max ) + myassert( rel_error_bar < rel_error_max ) \ No newline at end of file diff --git a/Scripts/tests/fermi_dirac_test.py b/Scripts/tests/fermi_dirac_test.py index ad7a5206..1459000d 100644 --- a/Scripts/tests/fermi_dirac_test.py +++ b/Scripts/tests/fermi_dirac_test.py @@ -18,12 +18,16 @@ # Sort the data file names by time step number directories = sorted(directories, key=lambda x: int(x.split("plt")[1].split(".")[0])) +# +# Energy bins from NuLib table (Cutted) +# The following energy bins should be the same as the st9_empty_particles_multi_energy initial condition script +# # Energy bin centers extracted from NuLib table -energies_center_Mev = np.array([1, 3, 5.23824, 8.00974, 11.4415, 15.6909, 20.9527, 27.4681, 35.5357, 45.5254, 57.8951, 73.2117, 92.1775, 115.662, 144.741, 180.748, 225.334, 280.542]) # Energy in Mev +energies_center_Mev = np.array([15.6909, 20.9527, 27.4681, 35.5357, 45.5254, 57.8951, 73.2117, 92.1775, 115.662, 144.741, 180.748, 225.334, 280.542]) # Energy in Mev # Energy bin bottom extracted from NuLib table -energies_bottom_Mev = np.array([0, 2, 4, 6.47649, 9.54299, 13.3401, 18.0418, 23.8636, 31.0725, 39.9989, 51.0519, 64.7382, 81.6853, 102.67, 128.654, 160.828, 200.668, 250]) +energies_bottom_Mev = np.array([13.3401, 18.0418, 23.8636, 31.0725, 39.9989, 51.0519, 64.7382, 81.6853, 102.67, 128.654, 160.828, 200.668, 250]) # Energy bin top extracted from NuLib table -energies_top_Mev = np.array([2, 4, 6.47649, 9.54299, 13.3401, 18.0418, 23.8636, 31.0725, 39.9989, 51.0519, 64.7382, 81.6853, 102.67, 128.654, 160.828, 200.668, 250, 311.085]) +energies_top_Mev = np.array([18.0418, 23.8636, 31.0725, 39.9989, 51.0519, 64.7382, 81.6853, 102.67, 128.654, 160.828, 200.668, 250, 311.085]) # Energies in ergs energies_center_erg = np.array(energies_center_Mev) * 1e6 * eV # Energy in ergs diff --git a/Source/DataReducer.cpp b/Source/DataReducer.cpp index b90799e1..fa5cb148 100644 --- a/Source/DataReducer.cpp +++ b/Source/DataReducer.cpp @@ -141,7 +141,7 @@ DataReducer::WriteReducedData0D(const amrex::Geometry& geom, Real TrN = 0; Real Vphase = p.rdata(PIdx::Vphase); #include "generated_files/DataReducer.cpp_fill_particles" - return GpuTuple{TrN,TrHN, Vphase}; + return amrex::GpuTuple{TrN, TrHN, Vphase}; }, reduce_ops); Real TrN = amrex::get<0>(particleResult); Real TrHN = amrex::get<1>(particleResult); @@ -150,7 +150,7 @@ DataReducer::WriteReducedData0D(const amrex::Geometry& geom, ParallelDescriptor::ReduceRealSum(TrHN); ParallelDescriptor::ReduceRealSum(Vphase); - printf("TrN=%g, TrHN=%g, Vphase=%g\n", TrN, TrHN, Vphase); + // printf("TrN=%g, TrHN=%g, Vphase=%g\n", TrN, TrHN, Vphase); //=============================// // Do reductions over the grid // diff --git a/Source/Evolve.H b/Source/Evolve.H index 61e5ada3..24ff2610 100644 --- a/Source/Evolve.H +++ b/Source/Evolve.H @@ -26,27 +26,14 @@ namespace GIdx void Initialize(); } -amrex::Real compute_dt(const amrex::Geometry& geom, const MultiFab& state, const FlavoredNeutrinoContainer& neutrinos, const TestParams* parms); +MultiFab compute_max_IMFP(const Geometry& geom, const MultiFab& state, const TestParams* parms); + +amrex::Real compute_dt(const amrex::Geometry& geom, MultiFab& state, FlavoredNeutrinoContainer& neutrinos, FlavoredNeutrinoContainer& neutrinos_dt, const TestParams* parms, const MultiFab& maximum_IMFP_abs); void deposit_to_mesh(const FlavoredNeutrinoContainer& neutrinos, amrex::MultiFab& state, const amrex::Geometry& geom); void interpolate_rhs_from_mesh(FlavoredNeutrinoContainer& neutrinos_rhs, const amrex::MultiFab& state, const amrex::Geometry& geom, const TestParams* parms); -/** - * @brief Sets the N and Nbar to zero for particles inside the black hole or boundary cells. - * - * This function iterates over all particles in the `FlavoredNeutrinoContainer` and sets N and Nbar to zero if pariticles are inside the black hole or within the boundary cells of the simulation domain. - * - * @param neutrinos Reference to the container holding the flavored neutrinos. - * @param parms Pointer to the structure containing test parameters, including black hole properties and domain dimensions. - * - * The function performs the following steps: - * - Iterates over all particles in the container. - * - Computes the distance of each particle from the black hole center. - * - Sets N and Nbar to zero if the particle is inside the black hole radius. - * - Sets N and Nbar to zero if the particle is within the boundary cells of the simulation domain. - * - */ void empty_particles_at_boundary_cells(FlavoredNeutrinoContainer& neutrinos, const TestParams* parms); #endif diff --git a/Source/Evolve.cpp b/Source/Evolve.cpp index eb939812..c82f2ff0 100644 --- a/Source/Evolve.cpp +++ b/Source/Evolve.cpp @@ -26,86 +26,477 @@ namespace GIdx } } -Real compute_dt(const Geometry& geom, const MultiFab& state, const FlavoredNeutrinoContainer& /* neutrinos */, const TestParams* parms) -{ - AMREX_ASSERT(parms->cfl_factor > 0.0 || parms->flavor_cfl_factor > 0.0 || parms->collision_cfl_factor > 0.0); +/** + * @brief Computes the maximum of three values. + * + * This function takes three arguments of the same type and returns the maximum value among them. + * It uses the standard library's `max` function to perform the comparison. + * + * @tparam T The type of the input values. Must support comparison via the `max` function. + * @param a The first value to compare. + * @param b The second value to compare. + * @param c The third value to compare. + * @return The maximum value among the three input values. + */ +template static inline AMREX_GPU_DEVICE T max3(T a, T b, T c) { + using std::max; + return max(max(a, b), c); +} + +/** + * @brief Computes the maximum of four values. + * + * This function takes four values of the same type and returns the maximum + * value among them. It uses the standard library's `max` function to compare + * the values. + * + * @tparam T The type of the values to be compared. This type must support + * comparison using the `max` function. + * @param a The first value to compare. + * @param b The second value to compare. + * @param c The third value to compare. + * @param d The fourth value to compare. + * @return The maximum value among the four input values. + */ +template static inline AMREX_GPU_DEVICE T max4(T a, T b, T c, T d) { + using std::max; + return max(max(a, b), max(c, d)); +} + +/** + * @brief Computes the maximum value among six given values. + * + * This function takes six input values of type T and returns the maximum + * value among them. It utilizes the `max3` function to find the maximum + * of three values and then compares the results of two such calls to + * determine the overall maximum. + * + * @tparam T The type of the input values. It should support comparison + * operations. + * @param a The first value to compare. + * @param b The second value to compare. + * @param c The third value to compare. + * @param d The fourth value to compare. + * @param e The fifth value to compare. + * @param f The sixth value to compare. + * @return The maximum value among the six input values. + */ +template static inline AMREX_GPU_DEVICE T max6(T a, T b, T c, T d, T e, T f) { + using std::max; + return max(max3(a, b, c), max3(d, e, f)); +} + +/** + * @brief Computes the maximum Inverse Mean Free Path (IMFP) for absorption and emission processes. + * + * This function calculates the maximum IMFP for absorption processes on every cell center based on the specified method in the parameters. + * It supports two methods: + * 1. Using IMFP values from the input file. + * 2. Using a NuLib table to interpolate IMFP values. + * + * @param geom The geometry of the computational domain. + * @param state The MultiFab containing the state variables. + * @param parms Pointer to the TestParams structure containing the parameters for the computation. + * @return A MultiFab containing the maximum IMFP values for each cell center. + */ +MultiFab compute_max_IMFP(const Geometry& geom, const MultiFab& state, const TestParams* parms){ + + const amrex::IntVect ngrow(0, 0, 0); // We do not need ghost cells for IMFP + MultiFab mf_IMFP(state.boxarray, state.DistributionMap(), 1, ngrow); //ncomp=1 + mf_IMFP.setVal(0.0); + + // If IMFP_method is 1, use the IMFPs from the input file and find the maximum absorption IMFP + if (parms->IMFP_method == 1) { + + // Use the IMFPs from the input file and find the maximum absorption IMFP + double max_IMFP_abs = std::numeric_limits::lowest(); + for (int i = 0; i < 2; ++i) { + for (int j = 0; j < NUM_FLAVORS; ++j) { + max_IMFP_abs = std::max(max_IMFP_abs, parms->IMFP_abs[i][j]); + } + } + + mf_IMFP.setVal(max_IMFP_abs); + + // If IMFP_method is 2, use the NuLib table to find the maximum absorption IMFP + } else if (parms->IMFP_method == 2) + { + //Create NuLib table object + using namespace nulib_private; + NuLib_tabulated NuLib_tabulated_obj(alltables_nulib, logrho_nulib, logtemp_nulib, + yes_nulib, helperVarsReal_nulib, helperVarsInt_nulib); + + int start_comp = GIdx::rho; + int num_comps = 3; //We only want to get GIdx::rho, GIdx::T and GIdx::Ye + MultiFab rho_T_ye_state(state, amrex::make_alias, start_comp, num_comps); + + for(amrex::MFIter mfi(rho_T_ye_state); mfi.isValid(); ++mfi){ + const amrex::Box& bx = mfi.validbox(); + const amrex::Array4& mf_array = rho_T_ye_state.array(mfi); + const amrex::Array4& mf_IMFP_array = mf_IMFP.array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k){ + + //Get the values from the input arrays + const Real rho = mf_array(i, j, k, GIdx::rho - start_comp); // g/ccm + const Real T_grid = mf_array(i, j, k, GIdx::T - start_comp); //erg + const Real Ye = mf_array(i, j, k, GIdx::Ye - start_comp); + + const Real temperature = T_grid/ (1e6*CGSUnitsConst::eV); //convert temperature from erg to MeV. + + double max_IMFP_abs_local = std::numeric_limits::lowest(); // Initialize to the lowest possible value + + Real IMFP_abs[NUM_FLAVORS][NUM_FLAVORS]; // Neutrino inverse mean free path matrix for nucleon absortion: diag( k_e , k_u , k_t ) + Real IMFP_absbar[NUM_FLAVORS][NUM_FLAVORS]; // Antineutrino inverse mean free path matrix for nucleon absortion: diag( kbar_e , kbar_u , kbar_t ) + + //--------------------- Values from NuLib table --------------------------- + int keyerr, anyerr; + double *helperVarsReal_nulib = NuLib_tabulated_obj.get_helperVarsReal_nulib(); + int idx_group = NULIBVAR(idx_group); + //FIXME: specify neutrino energy using the following: + // double neutrino_energy = p.rdata(PIdx::pupt); locate energy bin using this. + + //idx_species = {0 for electron neutrino, 1 for electron antineutrino and 2 for all other heavier ones} + //electron neutrino: [0, 0] + int idx_species = 0; + double absorption_opacity, scattering_opacity; + NuLib_tabulated_obj.get_opacities(rho, temperature, Ye, absorption_opacity, scattering_opacity, + keyerr, anyerr, idx_species, idx_group); + if (anyerr) assert(0); + + #ifdef DEBUG_INTERPOLATION_TABLES + printf("(Evolve.cpp) absorption_opacity[e] interpolated = %17.6g\n", absorption_opacity); + printf("(Evolve.cpp) scattering_opacity[e] interpolated = %17.6g\n", scattering_opacity); + #endif + + IMFP_abs[0][0] = absorption_opacity; + + //electron antineutrino: [1, 0] + idx_species = 1; + NuLib_tabulated_obj.get_opacities(rho, temperature, Ye, absorption_opacity, scattering_opacity, + keyerr, anyerr, idx_species, idx_group); + if (anyerr) assert(0); + + #ifdef DEBUG_INTERPOLATION_TABLES + printf("(Evolve.cpp) absorption_opacity[a] interpolated = %17.6g\n", absorption_opacity); + printf("(Evolve.cpp) scattering_opacity[a] interpolated = %17.6g\n", scattering_opacity); + #endif + + IMFP_absbar[0][0] = absorption_opacity; + + //heavier ones: muon neutrino[0,1], muon antineutruino[1,1], tau neutrino[0,2], tau antineutrino[1,2] + idx_species = 2; + NuLib_tabulated_obj.get_opacities(rho, temperature, Ye, absorption_opacity, scattering_opacity, + keyerr, anyerr, idx_species, idx_group); + if (anyerr) assert(0); + + #ifdef DEBUG_INTERPOLATION_TABLES + printf("(Evolve.cpp) absorption_opacity[x] interpolated = %17.6g\n", absorption_opacity); + printf("(Evolve.cpp) scattering_opacity[x] interpolated = %17.6g\n", scattering_opacity); + #endif + + for (int i=1; ineutrino or 1->antineutrino + IMFP_abs[i][i] = absorption_opacity ; + IMFP_absbar[i][i] = absorption_opacity ; + } - // translation part of timestep limit + //Calculate max of all IMFP_abs and IMFP_absbar. + if (NUM_FLAVORS == 2) { + max_IMFP_abs_local = max4(IMFP_abs[0][0], IMFP_abs[1][1], + IMFP_absbar[0][0], IMFP_absbar[1][1]); + } else if (NUM_FLAVORS == 3) { + max_IMFP_abs_local = max6(IMFP_abs[0][0], IMFP_abs[1][1], IMFP_abs[2][2], + IMFP_absbar[0][0], IMFP_absbar[1][1], IMFP_absbar[2][2]); + } + //----------------------------------------------------------------------- + mf_IMFP_array(i, j, k) = max_IMFP_abs_local; + }); + } + } + // Return the MultiFab with the maximum IMFP values + return mf_IMFP; +} + +/** + * @brief Computes the time step size for the simulation based on various CFL factors and conditions. + * + * This function calculates the time step size considering translation, flavor, and collision CFL factors. + * It ensures that the time step is limited by the smallest of these factors to maintain stability. + * + * @param geom The geometry of the simulation domain. + * @param state The state MultiFab containing the simulation data. + * @param neutrinos The container for flavored neutrinos. + * @param parms Pointer to the structure containing simulation parameters. + * @param maximum_IMFP_abs The MultiFab containing the maximum inverse mean free path for absorption. + * + * @return The computed time step size. + * + * @note At least one of cfl_factor, flavor_cfl_factor, or collision_cfl_factor must be greater than 0.0. + * @note The function handles periodic boundary conditions and black hole regions if specified in the parameters. + * @note The function performs a reduction operation to find the minimum time step across all cells and MPI ranks. + */ +Real compute_dt( + const Geometry& geom, + MultiFab& state, + FlavoredNeutrinoContainer& neutrinos, + FlavoredNeutrinoContainer& neutrinos_dt, + const TestParams* parms, + const MultiFab& maximum_IMFP_abs) +{ + // Initialize the maximum real value + const Real max_real = std::numeric_limits::max(); + + // Get the cell size array const auto dxi = geom.CellSizeArray(); Real dt_translation = 0.0; if (parms->cfl_factor > 0.0) { - dt_translation = std::min(std::min(dxi[0],dxi[1]), dxi[2]) / PhysConst::c * parms->cfl_factor; + // Calculate the time step size based on the translation CFL factor + // dt = (min(dx,dy,dz)/c) * cfl_factor + dt_translation = std::min({dxi[0], dxi[1], dxi[2]}) / PhysConst::c * parms->cfl_factor; } Real dt_flavor = 0.0; - if (parms->flavor_cfl_factor > 0.0 && parms->collision_cfl_factor > 0.0) { - // define the reduction operator to get the max contribution to - // the potential from matter and neutrinos - // compute "effective" potential (ergs) that produces characteristic timescale - // when multiplied by hbar - ReduceOps reduce_op; - ReduceData reduce_data(reduce_op); - using ReduceTuple = typename decltype(reduce_data)::Type; - for (MFIter mfi(state); mfi.isValid(); ++mfi) { - const Box& bx = mfi.fabbox(); - auto const& fab = state.array(mfi); - reduce_op.eval(bx, reduce_data, - [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple - { - Real V_adaptive=0, V_adaptive2=0, V_stupid=0; - #include "generated_files/Evolve.cpp_compute_dt_fill" - return {V_adaptive, V_stupid}; - }); - } - // extract the reduced values from the combined reduced data structure - auto rv = reduce_data.value(); - Real Vmax_adaptive = amrex::get<0>(rv) + FlavoredNeutrinoContainer::Vvac_max; - Real Vmax_stupid = amrex::get<1>(rv) + FlavoredNeutrinoContainer::Vvac_max; + if ( parms->time_step_method == 0 ){ + + AMREX_ASSERT_WITH_MESSAGE(parms->cfl_factor > 0.0 || parms->flavor_cfl_factor > 0.0 || parms->collision_cfl_factor > 0.0, + "Error: At least one of cfl_factor, flavor_cfl_factor, or collision_cfl_factor must be greater than 0.0."); + + if (parms->flavor_cfl_factor > 0.0 && parms->collision_cfl_factor > 0.0) { + + ReduceOps< ReduceOpMin, ReduceOpMin, ReduceOpMin > reduce_op; + ReduceData< Real , Real, Real > reduce_data(reduce_op); + using ReduceTuple = typename decltype(reduce_data)::Type; + for (MFIter mfi(state); mfi.isValid(); ++mfi) { + const Box& bx = mfi.validbox(); + auto const& fab = state.array(mfi); + auto const& multifab_IMFP = maximum_IMFP_abs.array(mfi); + Real V_vac_max = FlavoredNeutrinoContainer::Vvac_max; + reduce_op.eval(bx, reduce_data, + [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple + { + // Check if the cell is at the boundary or inside the black hole + if (parms->do_periodic_empty_bc == 1) { + // Check if the cell is at the boundary + if (i == 0 || i == parms->ncell[0] - 1 || + j == 0 || j == parms->ncell[1] - 1 || + k == 0 || k == parms->ncell[2] - 1) { + return {max_real,max_real,max_real}; + } + }else if (parms->do_blackhole == 1) { + // Check if the cell is inside the black hole + + // Calculate the cell size + double cell_size_x = parms->Lx / parms->ncell[0]; + double cell_size_y = parms->Ly / parms->ncell[1]; + double cell_size_z = parms->Lz / parms->ncell[2]; + // Calculate the cell center coordinates + double x_cell_center = (i + 0.5) * cell_size_x; + double y_cell_center = (j + 0.5) * cell_size_y; + double z_cell_center = (k + 0.5) * cell_size_z; + // Calculate the distance from the black hole center + Real distance_from_bh = sqrt(pow(x_cell_center - parms->bh_center_x, 2) + pow(y_cell_center - parms->bh_center_y, 2) + pow(z_cell_center - parms->bh_center_z, 2)); + + // Check if the cell is inside the black hole + if (distance_from_bh < parms->bh_radius) { + return {max_real,max_real,max_real}; + } + } + + // V_stupid = Max(N_ab,Nbar_ab, Ye*rho/Mp)*4.0*sqrt(2)*GF + // V_adaptive id the magnitud of vector the following vector + // |vec{H}| = | sqrt(2)*GF * ( (N_ab - Nbar_ab) + (F - Fbar) + Ye*rho/Mp ) | + + Real V_adaptive=0, V_adaptive2=0, V_stupid=0; + #include "generated_files/Evolve.cpp_compute_dt_fill" + + V_adaptive += V_vac_max; + V_stupid += V_vac_max; + + V_adaptive *= parms->attenuation_hamiltonians; + V_stupid *= parms->attenuation_hamiltonians; + + Real dt_adaptive = max_real; + Real dt_stupid = max_real; + Real dt_absorption = max_real; + + // Ensure that the minimum trace is not zero + if (std::abs(V_adaptive) > 0.0){ + dt_adaptive = parms->flavor_cfl_factor * ( PhysConst::hbar / std::abs(V_adaptive) ) ; + dt_stupid = parms->flavor_cfl_factor * ( PhysConst::hbar / std::abs(V_stupid ) ) ; + } + + // Calculate the absorption time step + if ( multifab_IMFP(i, j, k) > 0.0 ) { + dt_absorption = parms->collision_cfl_factor * ( 1.0 / ( PhysConst::c * multifab_IMFP(i, j, k) ) ); + } + + return {dt_adaptive, dt_stupid, dt_absorption}; + }); + } - // reduce across MPI ranks - ParallelDescriptor::ReduceRealMax(Vmax_adaptive); - ParallelDescriptor::ReduceRealMax(Vmax_stupid ); + // extract the reduced values from the combined reduced data structure + auto rv = reduce_data.value(); + Real min_dt_adaptive = amrex::get<0>(rv); + Real min_dt_stupid = amrex::get<1>(rv); + Real min_dt_absorption = amrex::get<2>(rv); - // define the dt associated with each method - Real dt_flavor_adaptive = std::numeric_limits::max(); - Real dt_flavor_stupid = std::numeric_limits::max(); - Real dt_flavor_absorption = std::numeric_limits::max(); // Initialize with infinity + // reduce across MPI ranks + ParallelDescriptor::ReduceRealMin(min_dt_adaptive); + ParallelDescriptor::ReduceRealMin(min_dt_stupid ); + ParallelDescriptor::ReduceRealMin(min_dt_absorption); - if (parms->attenuation_hamiltonians != 0) { - dt_flavor_adaptive = PhysConst::hbar / Vmax_adaptive * parms->flavor_cfl_factor / parms->attenuation_hamiltonians; - dt_flavor_stupid = PhysConst::hbar / Vmax_stupid * parms->flavor_cfl_factor / parms->attenuation_hamiltonians; - } + // define the dt associated with each method + Real dt_flavor_adaptive = max_real; + Real dt_flavor_stupid = max_real; + Real dt_flavor_absorption = max_real; // Initialize with infinity - if (parms->IMFP_method == 1) { - // Use the IMFPs from the input file and find the maximum absorption IMFP - double max_IMFP_abs = std::numeric_limits::lowest(); // Initialize max to lowest possible value - for (int i = 0; i < 2; ++i) { - for (int j = 0; j < NUM_FLAVORS; ++j) { - max_IMFP_abs = std::max(max_IMFP_abs, parms->IMFP_abs[i][j]); - } + if (parms->IMFP_method == 1 || parms->IMFP_method == 2) { + dt_flavor_absorption = min_dt_absorption; + } + if (parms->attenuation_hamiltonians != 0) { + dt_flavor_adaptive = min_dt_adaptive; + dt_flavor_stupid = min_dt_stupid; + } + + // pick the appropriate timestep + dt_flavor = min(dt_flavor_stupid, dt_flavor_adaptive, dt_flavor_absorption); + if(parms->max_adaptive_speedup > 1) { + dt_flavor = min(dt_flavor_stupid*parms->max_adaptive_speedup, dt_flavor_adaptive, dt_flavor_absorption); } - // Calculate dt_flavor_absorption - dt_flavor_absorption = (1 / (PhysConst::c * max_IMFP_abs)) * parms->collision_cfl_factor; } + } else if ( parms->time_step_method == 1 ){ + + AMREX_ASSERT_WITH_MESSAGE(parms->cfl_factor > 0.0 || parms->flavor_cfl_factor > 0.0, + "Error: At least one of cfl_factor, flavor_cfl_factor, or collision_cfl_factor must be greater than 0.0."); + + deposit_to_mesh(neutrinos, state, geom); + state.FillBoundary(geom.periodicity()); + neutrinos_dt.copyParticles(neutrinos, true); + interpolate_rhs_from_mesh(neutrinos_dt, state, geom, parms); - // pick the appropriate timestep - dt_flavor = min(dt_flavor_stupid, dt_flavor_adaptive, dt_flavor_absorption); - if(parms->max_adaptive_speedup>1) { - dt_flavor = min(dt_flavor_stupid*parms->max_adaptive_speedup, dt_flavor_adaptive, dt_flavor_absorption); + ReduceOps reduce_op; + ReduceData reduce_data(reduce_op); + + const int lev = 0; + FNParIter pti_time_derivative(neutrinos_dt, lev); + for (FNParIter pti(neutrinos, lev); pti.isValid(); ++pti) + { + + const int np = pti.numParticles(); + const int np_dt = pti_time_derivative.numParticles(); + + FlavoredNeutrinoContainer::ParticleType* pstruct = &(pti.GetArrayOfStructs()[0]); + FlavoredNeutrinoContainer::ParticleType* pstruct_dt = &(pti_time_derivative.GetArrayOfStructs()[0]); + + reduce_op.eval(np, reduce_data, + [=] AMREX_GPU_DEVICE (int i) -> Real + { + + FlavoredNeutrinoContainer::ParticleType& p = pstruct[i]; + FlavoredNeutrinoContainer::ParticleType& p_dt = pstruct_dt[i]; + + auto update_dt = [&] (int idx) { + if (std::isnan(p.rdata(idx))) { + amrex::Abort("Error: NaN value detected in N or Nbar."); + p_dt.rdata(idx) = -1.0; + return; + } + if (std::abs(p_dt.rdata(idx)) > 0) { + if ( std::abs(p.rdata(idx)) > 1.0){ + p_dt.rdata(idx) = parms->flavor_cfl_factor * std::abs( p.rdata(idx) / p_dt.rdata(idx)); + } else{ + p_dt.rdata(idx) = parms->flavor_cfl_factor * std::abs( 1.0 / p_dt.rdata(idx)); + } + }else{ + p_dt.rdata(idx) = max_real; + } + }; + + update_dt(PIdx::N00_Rebar); + update_dt(PIdx::N01_Rebar); + update_dt(PIdx::N01_Imbar); + update_dt(PIdx::N11_Rebar); + update_dt(PIdx::N00_Re); + update_dt(PIdx::N01_Re); + update_dt(PIdx::N01_Im); + update_dt(PIdx::N11_Re); + + #if NUM_FLAVORS == 3 + update_dt(PIdx::N02_Rebar); + update_dt(PIdx::N02_Imbar); + update_dt(PIdx::N12_Rebar); + update_dt(PIdx::N12_Imbar); + update_dt(PIdx::N22_Rebar); + update_dt(PIdx::N02_Re); + update_dt(PIdx::N02_Im); + update_dt(PIdx::N12_Re); + update_dt(PIdx::N12_Im); + update_dt(PIdx::N22_Re); + #endif + + Real min_dt = std::min({ + p_dt.rdata(PIdx::N00_Rebar), + p_dt.rdata(PIdx::N01_Rebar), + p_dt.rdata(PIdx::N01_Imbar), + p_dt.rdata(PIdx::N11_Rebar), + p_dt.rdata(PIdx::N00_Re), + p_dt.rdata(PIdx::N01_Re), + p_dt.rdata(PIdx::N01_Im), + p_dt.rdata(PIdx::N11_Re) + }); + + #if NUM_FLAVORS == 3 + min_dt = std::min({ + min_dt, + p_dt.rdata(PIdx::N02_Rebar), + p_dt.rdata(PIdx::N02_Imbar), + p_dt.rdata(PIdx::N12_Rebar), + p_dt.rdata(PIdx::N12_Imbar), + p_dt.rdata(PIdx::N22_Rebar), + p_dt.rdata(PIdx::N02_Re), + p_dt.rdata(PIdx::N02_Im), + p_dt.rdata(PIdx::N12_Re), + p_dt.rdata(PIdx::N12_Im), + p_dt.rdata(PIdx::N22_Re) + }); + #endif + return min_dt; + + }); + + ++pti_time_derivative; } + + // Extract the reduced value + Real qke_dt = amrex::get<0>(reduce_data.value()); + + // Reduce across MPI ranks + ParallelDescriptor::ReduceRealMin(qke_dt); + + dt_flavor = qke_dt; + } + + if (dt_flavor < 0.0) { + printf("Error: NaN value detected in N or Nbar. Aborting...\n"); + std::abort(); } Real dt = 0.0; if (dt_translation != 0.0 && dt_flavor != 0.0) { dt = std::min(dt_translation, dt_flavor); - } else if (dt_translation != 0.0) { - dt = dt_translation; - } else if (dt_flavor != 0.0) { - dt = dt_flavor; } else { - amrex::Error("Timestep selection failed, try using both cfl_factor and flavor_cfl_factor"); + if (dt_translation != 0.0) { + dt = dt_translation; + } else if (dt_flavor != 0.0) { + dt = dt_flavor; + } else { + amrex::Error("Timestep selection failed, both dt_translation and dt_flavor are zero. Try using both cfl_factor and flavor_cfl_factor."); + } } + + if (dtminimum_time_step) dt = parms->minimum_time_step; + // printf("dt = %g, dt_flavor = %g, dt_translation = %g\n", dt, dt_flavor, dt_translation); return dt; } @@ -167,7 +558,7 @@ void interpolate_rhs_from_mesh(FlavoredNeutrinoContainer& neutrinos_rhs, const M NuLib_tabulated NuLib_tabulated_obj(alltables_nulib, logrho_nulib, logtemp_nulib, yes_nulib, helperVarsReal_nulib, helperVarsInt_nulib); - + NuLib_energies NuLib_energies_obj(energy_bottom, energy_top); //The following commented loop can be used to print information about each particle in case debug is needed in future. /* @@ -344,9 +735,25 @@ void interpolate_rhs_from_mesh(FlavoredNeutrinoContainer& neutrinos_rhs, const M //--------------------- Values from NuLib table --------------------------- double *helperVarsReal_nulib = NuLib_tabulated_obj.get_helperVarsReal_nulib(); - int idx_group = NULIBVAR(idx_group); - //FIXME: specify neutrino energy using the following: - // double neutrino_energy = p.rdata(PIdx::pupt); locate energy bin using this. + int *helperVarsInt_nulib = NuLib_tabulated_obj.get_helperVarsInt_nulib(); + + double *energy_bottom = NuLib_energies_obj.get_energy_bottom_nulib(); + double *energy_top = NuLib_energies_obj.get_energy_top_nulib(); + + double neutrino_energy_erg = p.rdata(PIdx::pupt); //locate energy bin using this. + double neutrino_energy_MeV = neutrino_energy_erg / (1e6*CGSUnitsConst::eV); + + //Decide which energy bin to use (i.e. determine 'idx_group') + int idx_group = -1; + for (int i=0; i= energy_bottom[i] && neutrino_energy_MeV <= energy_top[i]){ + idx_group = i; + break; + } + } + + if(idx_group == -1) assert(0); //abort if energy bin cannot be found. + //printf("Given neutrino energy = %f, selected bin index = %d\n", neutrino_energy_MeV, idx_group); //idx_species = {0 for electron neutrino, 1 for electron antineutrino and 2 for all other heavier ones} //electron neutrino: [0, 0] @@ -401,19 +808,59 @@ void interpolate_rhs_from_mesh(FlavoredNeutrinoContainer& neutrinos_rhs, const M } else AMREX_ASSERT_WITH_MESSAGE(false, "only available opacity_method is 0, 1 or 2"); - // Compute equilibrium distribution functions and include Pauli blocking term if requested - if(parms->IMFP_method==1 || parms->IMFP_method==2){ + double f_nu_e[13] = { + 5.3245656748033138e-03, 1.5978206640139703e-02, 2.7340626736756864e-02, 3.2089846556718205e-02, 2.6210613940471663e-02, 1.6053556351649751e-02, 7.4581054456212296e-03, 2.5160179998209029e-03, 5.8403718278747475e-04, 8.7777097910152665e-05, 7.9876592427094722e-06, 4.0330773860229255e-07, 9.9487753023125056e-09 + }; - for (int i=0; iDo_Pauli_blocking == 1){ - IMFP_abs[i][i] = IMFP_abs[i][i] / ( 1 - f_eq[i][i] ) ; // Multiply the absortion inverse mean free path by the Pauli blocking term 1 / (1 - f_eq). - IMFP_absbar[i][i] = IMFP_absbar[i][i] / ( 1 - f_eqbar[i][i] ) ; // Multiply the absortion inverse mean free path by the Pauli blocking term 1 / (1 - f_eq). + double energybinsbottom_ergs[13] = { + 0.0000000000000000e+00, 3.2043599999999999e-06, 6.4087199999999999e-06, 1.0376518769999999e-05, 1.5289603739999999e-05, 2.1373081199999998e-05, 2.8906531560000000e-05, 3.8234423520000002e-05, 4.9784539140000002e-05, 6.4085597819999999e-05, 8.1794493359999999e-05, 1.0372192884000000e-04, 1.3087407330000000e-04 + }; + + double energybinstopMeV_ergs[13] = { + 3.2043599999999999e-06, 6.4087199999999999e-06, 1.0376518769999999e-05, 1.5289603739999999e-05, 2.1373081199999998e-05, 2.8906531560000000e-05, 3.8234423520000002e-05, 4.9784539140000002e-05, 6.4085597819999999e-05, 8.1794493359999999e-05, 1.0372192884000000e-04, 1.3087407330000000e-04, 1.6449582060000000e-04 + }; + + double energybinsMeV_ergs[13] = { + 1.6021800000000000e-06, 4.8065399999999999e-06, 8.3925392760000003e-06, 1.2832981146000000e-05, 1.8332143560000002e-05, 2.5139806379999999e-05, 3.3570477540000001e-05, 4.4008680239999999e-05, 5.6935068480000001e-05, 7.2939244500000000e-05, 9.2758211099999997e-05, 1.1729880216000001e-04, 1.4768574804000000e-04 + }; + + // Compute equilibrium distribution functions and include Pauli blocking term if requested + if(parms->IMFP_method==1 || parms->IMFP_method==2){ + + if (parms-> set_equilibrium_distribution == 1){ + for (int bin = 0; bin < 13; ++bin) { + if (p.rdata(PIdx::pupt) >= energybinsbottom_ergs[bin] && p.rdata(PIdx::pupt) < energybinstopMeV_ergs[bin]) { + f_eq[0][0] = f_nu_e[bin]; + f_eqbar[0][0] = f_nubar_e[bin]; + f_eq[1][1] = f_nu_x[bin]; + f_eqbar[1][1] = f_nu_x[bin]; + #if NUM_FLAVORS == 3 + f_eq[2][2] = f_nu_x[bin]; + f_eqbar[2][2] = f_nu_x[bin]; + #endif + break; + } + } + } else { + for (int i=0; iDo_Pauli_blocking == 1){ + IMFP_abs[i][i] = IMFP_abs[i][i] / ( 1 - f_eq[i][i] ) ; // Multiply the absortion inverse mean free path by the Pauli blocking term 1 / (1 - f_eq). + IMFP_absbar[i][i] = IMFP_absbar[i][i] / ( 1 - f_eqbar[i][i] ) ; // Multiply the absortion inverse mean free path by the Pauli blocking term 1 / (1 - f_eq). + } } } } diff --git a/Source/IO.cpp b/Source/IO.cpp index 2a839ef8..28b6e939 100644 --- a/Source/IO.cpp +++ b/Source/IO.cpp @@ -38,13 +38,18 @@ WritePlotFile (const amrex::MultiFab& state, amrex::WriteSingleLevelPlotfile(plotfilename, state, GIdx::names, geom, time, step); #endif +FlavoredNeutrinoContainer neutrinos_for_IO(geom, state.DistributionMap(), state.boxArray()); +neutrinos_for_IO.copyParticles(neutrinos, true); + if (write_plot_particles == 1) { auto neutrino_varnames = neutrinos.get_attribute_names(); #ifdef AMREX_USE_HDF5 neutrinos.CheckpointHDF5(plotfilename, "neutrinos", true, neutrino_varnames); + // neutrinos_for_IO.WritePlotFile(plotfilename, "neutrinos", neutrino_varnames); #else neutrinos.Checkpoint(plotfilename, "neutrinos", true, neutrino_varnames); + // neutrinos_for_IO.WritePlotFile(plotfilename, "neutrinos", neutrino_varnames); #endif } diff --git a/Source/NuLibTable.H b/Source/NuLibTable.H index 262460ed..bbdb71b2 100644 --- a/Source/NuLibTable.H +++ b/Source/NuLibTable.H @@ -12,7 +12,7 @@ void ReadNuLibTable(const std::string nulib_table_name); namespace nulib_private { -// table data + // table data extern double *alltables_nulib; //extern double *epstable; extern double *logrho_nulib; @@ -20,6 +20,8 @@ namespace nulib_private { extern double *yes_nulib; extern double *species_nulib; extern double *group_nulib; + extern double *energy_bottom; + extern double *energy_top; #define NULIBVAR(VAR) helperVarsReal_nulib[helperVarsEnumReal_nulib::VAR] //WARNING: If the number of variables is increased here, then memory allocation for helperVarsReal/helperVarsInt pointer in GRHydroX_ReadEOSTable.cxx will also need to be increased. diff --git a/Source/NuLibTableFunctions.H b/Source/NuLibTableFunctions.H index 95f2c198..3ccd5bd7 100644 --- a/Source/NuLibTableFunctions.H +++ b/Source/NuLibTableFunctions.H @@ -19,6 +19,11 @@ struct NuLib_tabulated { AMREX_GPU_DEVICE AMREX_GPU_HOST double *get_helperVarsReal_nulib() const { return helperVarsReal_nulib; };//get_helperVarsReal_nulib + + //--------------- get helperVarsInt pointer --------------------------------------------- + AMREX_GPU_DEVICE AMREX_GPU_HOST int *get_helperVarsInt_nulib() const { + return helperVarsInt_nulib; + };//get_helperVarsInt_nulib //--------------- get opacaties --------------------------------------------- @@ -41,4 +46,26 @@ struct NuLib_tabulated { }; //struct NuLib_tabulated + +struct NuLib_energies { + double *energy_bottom, *energy_top; + + //constructor for NuLib_energies + AMREX_GPU_DEVICE AMREX_GPU_HOST NuLib_energies() = default;//Need to keep it + AMREX_GPU_DEVICE AMREX_GPU_HOST NuLib_energies(double *energy_bottom, double *energy_top): + energy_bottom(energy_bottom), energy_top(energy_top) {} + + //--------------- get energy_bottom pointer --------------------------------------------- + AMREX_GPU_DEVICE AMREX_GPU_HOST double *get_energy_bottom_nulib() const { + return energy_bottom; + };//get_energy_bottom_nulib + + //--------------- get energy_top pointer --------------------------------------------- + AMREX_GPU_DEVICE AMREX_GPU_HOST double *get_energy_top_nulib() const { + return energy_top; + };//get_energy_top_nulib + +}; //struct NuLib_energies + + #endif // NULIBTABLEFUCNTIONS_HXX \ No newline at end of file diff --git a/Source/Parameters.H b/Source/Parameters.H index c7ff098a..8b4c946f 100644 --- a/Source/Parameters.H +++ b/Source/Parameters.H @@ -17,7 +17,9 @@ struct TestParams : public amrex::Gpu::Managed IntVect ncell; // Number of cells in the domain IntVect nppc; // Number of points of particle emission per cell in each dimension Real Lx, Ly, Lz; // Length of the domain in each dimension - int max_grid_size; // Maximum grid size subdivisions of the domain for parallelization + int max_grid_size_x; // Maximum grid size subdivisions in x of the domain for parallelization + int max_grid_size_y; // Maximum grid size subdivisions in y of the domain for parallelization + int max_grid_size_z; // Maximum grid size subdivisions in z of the domain for parallelization int nsteps; // Number of time steps Real end_time; // End time of the simulation int write_plot_every; // Write plot files every n steps @@ -30,6 +32,8 @@ struct TestParams : public amrex::Gpu::Managed bool do_restart; // Flag to restart from a previous simulation std::string restart_dir; // Directory path to restart from in case do_restart is true Real maxError; + Real minimum_time_step; // Minimum time step allowed in the simulation + int time_step_method; // Time step method flag // File name with particles' initial conditions // This includes energy, momentum, N, Nbar, phase space volume, and others. @@ -56,6 +60,9 @@ struct TestParams : public amrex::Gpu::Managed // Attenuation to Hamiltonians Real attenuation_hamiltonians; // Multiplication factor to [N,H] + Real attenuation_vacuum; + Real attenuation_matter; + Real attenuation_self_interaction; // Black hole properties int do_blackhole; // Flag to include a black hole in the simulation @@ -74,6 +81,8 @@ struct TestParams : public amrex::Gpu::Managed std::string background_rho_Ye_T_table_name; // Background matter table file name int boundary_condition_type; // Boundary condition type flag + int set_equilibrium_distribution; // Flag to set the equilibrium distribution function for neutrinos and antineutrinos + void Initialize(){ ParmParse pp; @@ -82,7 +91,9 @@ struct TestParams : public amrex::Gpu::Managed pp.get("Ly", Ly); pp.get("Lz", Lz); pp.get("nppc", nppc); - pp.get("max_grid_size", max_grid_size); + pp.get("max_grid_size_x", max_grid_size_x); + pp.get("max_grid_size_y", max_grid_size_y); + pp.get("max_grid_size_z", max_grid_size_z); pp.get("nsteps", nsteps); pp.get("end_time", end_time); pp.get("cfl_factor", cfl_factor); @@ -94,6 +105,8 @@ struct TestParams : public amrex::Gpu::Managed pp.get("do_restart", do_restart); pp.get("restart_dir", restart_dir); pp.get("maxError", maxError); + pp.get("minimum_time_step", minimum_time_step); + pp.get("time_step_method", time_step_method); // File name with particles' initial conditions // This includes energy, momentum, N, Nbar, phase space volume, and others. @@ -129,7 +142,10 @@ struct TestParams : public amrex::Gpu::Managed // attenuation to hamiltonians pp.get("attenuation_hamiltonians", attenuation_hamiltonians); - + pp.get("attenuation_vacuum", attenuation_vacuum); + pp.get("attenuation_matter", attenuation_matter); + pp.get("attenuation_self_interaction", attenuation_self_interaction); + // Boundary condition type pp.get("boundary_condition_type", boundary_condition_type); pp.get("do_periodic_empty_bc", do_periodic_empty_bc); @@ -161,6 +177,11 @@ struct TestParams : public amrex::Gpu::Managed // Flag for collision term treatment pp.get("IMFP_method", IMFP_method); + if(IMFP_method>0){ + // Set equilibrium distribution function for neutrinos and antineutrinos + pp.get("set_equilibrium_distribution", set_equilibrium_distribution); + } + if(IMFP_method==0) { // Do nothing } else if(IMFP_method==1){ diff --git a/Source/ReadNuLibTable.cpp b/Source/ReadNuLibTable.cpp index bbed8577..f55e3c3e 100644 --- a/Source/ReadNuLibTable.cpp +++ b/Source/ReadNuLibTable.cpp @@ -46,6 +46,8 @@ namespace nulib_private { double *yes_nulib; double *species_nulib; //TODO: Get rid of this? double *group_nulib; + double *energy_bottom; + double *energy_top; double *helperVarsReal_nulib; int *helperVarsInt_nulib; } @@ -135,8 +137,6 @@ void ReadNuLibTable(const std::string nulib_table_name) { } //Allocate memory for energy bin determination. - double *energy_bottom; - double *energy_top; if (!(energy_bottom = myManagedArena.allocate(ngroup_) )) { printf("(ReadNuLibTable.cpp) Cannot allocate memory for NuLib table"); assert(0); @@ -198,26 +198,7 @@ void ReadNuLibTable(const std::string nulib_table_name) { for(int i=0;i= energy_bottom[i] && given_energy <= energy_top[i]){ - idx_group_ = i; - break; - } - } - - printf("Given neutrino energy = %f, selected bin index = %d\n", given_energy, idx_group); - myManagedArena.deallocate(energy_bottom, ngroup_); - myManagedArena.deallocate(energy_top, ngroup_); - //---------------------------------------------------------------------------------------------- + //allocate memory for helperVars helperVarsReal_nulib = myManagedArena.allocate(24); @@ -226,14 +207,43 @@ void ReadNuLibTable(const std::string nulib_table_name) { const double temp0_ = exp(logtemp_nulib[0]); const double temp1_ = exp(logtemp_nulib[1]); - NULIBVAR(idx_group) = idx_group_; - NULIBVAR_INT(nrho) = nrho_; NULIBVAR_INT(ntemp) = ntemp_; NULIBVAR_INT(nye) = nye_; NULIBVAR_INT(nspecies) = nspecies_; NULIBVAR_INT(ngroup) = ngroup_; + //Ensure that temp is uniformly spaced + double dtemp_old, dtemp; + dtemp = logtemp_nulib[1] - logtemp_nulib[0]; + for (int i=1; imax_grid_size); + const IntVect mgs(parms->max_grid_size_x, parms->max_grid_size_y, parms->max_grid_size_z); + ba.maxSize(mgs); + amrex::Print() << "Number of boxes created: " << ba.size() << std::endl; // This defines the physical box, [0,1] in each dimension RealBox real_box({AMREX_D_DECL( 0.0, 0.0, 0.0)}, @@ -138,6 +140,13 @@ void evolve_flavor(const TestParams* parms) amrex::Print() << "Reading NuLib table... " << std::endl; ReadNuLibTable(parms->nulib_table_name); } + + // Compute the maximum Inverse Mean Free Path (IMFP) in the domain + MultiFab max_IMFP_absortion; + max_IMFP_absortion.setVal(0.0); + if (parms->time_step_method == 0) { + max_IMFP_absortion = compute_max_IMFP(geom, state, parms); + } // Initialize particles on the domain amrex::Print() << "Initializing particles... " << std::endl; @@ -145,6 +154,7 @@ void evolve_flavor(const TestParams* parms) // We store old-time and new-time data FlavoredNeutrinoContainer neutrinos_old(geom, dm, ba); FlavoredNeutrinoContainer neutrinos_new(geom, dm, ba); + FlavoredNeutrinoContainer neutrinos_dt(geom, dm, ba); // Track the Figure of Merit for the simulation // defined as number of particles advanced per microsecond of walltime @@ -251,15 +261,17 @@ void evolve_flavor(const TestParams* parms) const int step = integrator.get_step_number(); const Real time = integrator.get_time(); - printf("Writing reduced data to file... \n"); - rd.WriteReducedData0D(geom, state, neutrinos, time, step+1); - printf("Done. \n"); + // Write the reduced data to file + rd.WriteReducedData0D(geom, state, neutrinos, time, step+1); run_fom += neutrinos.TotalNumberOfParticles(); + // Print the current step and time + printf("step = %d, t = %g s\n", step+1, time); + // Write the Mesh Data to Plotfile if required - bool write_plotfile = parms->write_plot_every > 0 && (step+1) % parms->write_plot_every == 0; - bool write_plot_particles = parms->write_plot_particles_every > 0 && (step+1) % parms->write_plot_particles_every == 0; + bool write_plotfile = parms->write_plot_every > 0 && (step+1) % parms->write_plot_every == 0; + bool write_plot_particles = parms->write_plot_particles_every > 0 && (step+1) % parms->write_plot_particles_every == 0; if (write_plotfile || write_plot_particles) { // Only include the Particle Data if write_plot_particles_every is satisfied int write_plot_particles = parms->write_plot_particles_every > 0 && @@ -271,11 +283,11 @@ void evolve_flavor(const TestParams* parms) // Note: this won't be the same as the new-time grid data // because the last deposit_to_mesh call was at either the old time (forward Euler) // or the final RK stage, if using Runge-Kutta. - printf("Setting next timestep... \n"); - const Real dt = compute_dt(geom, state, neutrinos, parms); + // printf("Setting next timestep... \n"); + const Real dt = compute_dt(geom, state, neutrinos, neutrinos_dt, parms, max_IMFP_absortion); integrator.set_timestep(dt); //printf("current_dt = %g, dt = %g \n", current_dt, dt); - printf("Done. \n"); + // printf("Done. \n"); }; // Attach our RHS and post timestep hooks to the integrator @@ -283,8 +295,8 @@ void evolve_flavor(const TestParams* parms) integrator.set_post_timestep(post_timestep_fun); // Get a starting timestep - const Real starting_dt = compute_dt(geom, state, neutrinos_old, parms); - + const Real starting_dt = compute_dt(geom, state, neutrinos_old, neutrinos_dt, parms, max_IMFP_absortion); + // Do all the science! amrex::Print() << "Starting timestepping loop... " << std::endl; diff --git a/sample_inputs/inputs_1d_fiducial b/sample_inputs/inputs_1d_fiducial index a8148b92..926bf50a 100644 --- a/sample_inputs/inputs_1d_fiducial +++ b/sample_inputs/inputs_1d_fiducial @@ -3,6 +3,9 @@ perturbation_amplitude = 1e-6 # attenuation parameters to time derivative of N due to hamiltonians attenuation_hamiltonians = 1 +attenuation_vacuum = 1.0 +attenuation_matter = 1.0 +attenuation_self_interaction = 1.0 collision_cfl_factor = 1e-3 cfl_factor = 0.5 @@ -24,11 +27,19 @@ nppc = (1, 1, 1) particle_data_filename = "particle_input.dat" # Maximum size of each grid in the domain -max_grid_size = 16 +max_grid_size_x = 1 +max_grid_size_y = 1 +max_grid_size_z = 64 # Number of steps to run nsteps = 1000 +# Minimum timestep (s) +minimum_time_step = 1e-15 + +# Time step method +time_step_method = 0 + # Simulation end time end_time = 5.0e-9 diff --git a/sample_inputs/inputs_bc_periodic_init b/sample_inputs/inputs_bc_periodic_init index bbc2ad0b..1b43043a 100644 --- a/sample_inputs/inputs_bc_periodic_init +++ b/sample_inputs/inputs_bc_periodic_init @@ -3,10 +3,13 @@ perturbation_amplitude = 0.0 # attenuation parameters to time derivative of N due to hamiltonians attenuation_hamiltonians = 0 +attenuation_vacuum = 1.0 +attenuation_matter = 1.0 +attenuation_self_interaction = 1.0 collision_cfl_factor = 1e-1 -cfl_factor = 2 -flavor_cfl_factor = 2 +cfl_factor = 1e-2 +flavor_cfl_factor = 1e-1 max_adaptive_speedup = 0 maxError = 1e-6 @@ -24,10 +27,18 @@ nppc = (1, 1, 1) particle_data_filename = "particle_input.dat" # Maximum size of each grid in the domain -max_grid_size = 16 +max_grid_size_x = 5 +max_grid_size_y = 5 +max_grid_size_z = 5 # Number of steps to run -nsteps = 30000 +nsteps = 5000 + +# Minimum timestep (s) +minimum_time_step = 1e-15 + +# Time step method +time_step_method = 0 # Simulation end time end_time = 5.0e19 @@ -41,10 +52,10 @@ T_MeV = 7.0 Ye = 0 # Write plotfiles -write_plot_every = 3000 +write_plot_every = 500 # Write particle data in plotfiles -write_plot_particles_every = 3000 +write_plot_particles_every = 500 # checkpointing do_restart = 0 @@ -86,6 +97,7 @@ deltaCP_degrees = 0 # opacity stuff # ################# IMFP_method = 1 +set_equilibrium_distribution= 0 Do_Pauli_blocking = 0 # If 1, it will multiply the inverse mean free path by 1 / (1 - f_eq); if 0, do nothing. diff --git a/sample_inputs/inputs_bipolar_test b/sample_inputs/inputs_bipolar_test index 9781389a..e86efd52 100644 --- a/sample_inputs/inputs_bipolar_test +++ b/sample_inputs/inputs_bipolar_test @@ -6,6 +6,9 @@ maxError = 1e-6 # attenuation parameters to time derivative of N due to hamiltonians attenuation_hamiltonians = 1 +attenuation_vacuum = 1.0 +attenuation_matter = 1.0 +attenuation_self_interaction = 1.0 perturbation_type = 0 perturbation_amplitude = 0 @@ -24,11 +27,19 @@ nppc = (1, 1, 1) particle_data_filename = "particle_input.dat" # Maximum size of each grid in the domain -max_grid_size = 64 +max_grid_size_x = 1 +max_grid_size_y = 1 +max_grid_size_z = 1 # Number of steps to run nsteps = 200 +# Minimum timestep (s) +minimum_time_step = 1e-15 + +# Time step method +time_step_method = 0 + # Simulation end time end_time = 1.0 diff --git a/sample_inputs/inputs_coll_equi_test b/sample_inputs/inputs_coll_equi_test index b2d603d3..212c4740 100644 --- a/sample_inputs/inputs_coll_equi_test +++ b/sample_inputs/inputs_coll_equi_test @@ -3,9 +3,12 @@ perturbation_amplitude = 1e-6 # attenuation parameters to time derivative of N due to hamiltonians attenuation_hamiltonians = 0 +attenuation_vacuum = 1.0 +attenuation_matter = 1.0 +attenuation_self_interaction = 1.0 -collision_cfl_factor = 0.04 -cfl_factor = 1e10 +collision_cfl_factor = 2.0 +cfl_factor = 2.0 flavor_cfl_factor = 0.5 max_adaptive_speedup = 0 maxError = 1e-6 @@ -24,11 +27,19 @@ nppc = (1, 1, 1) particle_data_filename = "particle_input.dat" # Maximum size of each grid in the domain -max_grid_size = 16 +max_grid_size_x = 2 +max_grid_size_y = 2 +max_grid_size_z = 2 # Number of steps to run nsteps = 1000 +# Minimum timestep (s) +minimum_time_step = 1e-15 + +# Time step method +time_step_method = 0 + # Simulation end time end_time = 5.0e9 @@ -86,6 +97,7 @@ deltaCP_degrees = 222 # opacity stuff # ################# IMFP_method = 1 +set_equilibrium_distribution= 0 Do_Pauli_blocking = 0 # If 1, it will multiply the inverse mean free path by 1 / (1 - f_eq); if 0, do nothing. diff --git a/sample_inputs/inputs_collisional_instability_test b/sample_inputs/inputs_collisional_instability_test index cc56120e..96cc616b 100644 --- a/sample_inputs/inputs_collisional_instability_test +++ b/sample_inputs/inputs_collisional_instability_test @@ -3,10 +3,13 @@ perturbation_amplitude = 0.0 # attenuation parameters to time derivative of N due to hamiltonians attenuation_hamiltonians = 1.0 +attenuation_vacuum = 1.0 +attenuation_matter = 1.0 +attenuation_self_interaction = 1.0 -collision_cfl_factor = 1e-3 -cfl_factor = 0.5 -flavor_cfl_factor = 0.5 +collision_cfl_factor = 1e-1 +cfl_factor = 1e-1 +flavor_cfl_factor = 1e0 max_adaptive_speedup = 0 maxError = 1e-6 @@ -15,19 +18,27 @@ integration.rk.type = 4 # Domain size in 3D index space ncell = (1, 1, 1) -Lx = 1 # cm -Ly = 1 # cm -Lz = 1 # cm +Lx = 1e4 # cm +Ly = 1e4 # cm +Lz = 1e4 # cm # Number of particles per cell nppc = (1, 1, 1) particle_data_filename = "particle_input.dat" # Maximum size of each grid in the domain -max_grid_size = 16 +max_grid_size_x = 1 +max_grid_size_y = 1 +max_grid_size_z = 1 # Number of steps to run -nsteps = 40000 +nsteps = 3000 + +# Minimum timestep (s) +minimum_time_step = 1e-15 + +# Time step method +time_step_method = 1 # Simulation end time end_time = 5.0e19 @@ -41,15 +52,17 @@ T_MeV = 7.0 Ye = 0 # Write plotfiles -write_plot_every = 500 +write_plot_every = 10 # Write particle data in plotfiles -write_plot_particles_every = 500 +write_plot_particles_every = 10 # checkpointing do_restart = 0 restart_dir = "" +set_equilibrium_distribution = 0 + ############################### # NEUTRINO PHYSICS PARAMETERS # ############################### diff --git a/sample_inputs/inputs_fast_flavor b/sample_inputs/inputs_fast_flavor index fdcdebe7..53de4f6c 100644 --- a/sample_inputs/inputs_fast_flavor +++ b/sample_inputs/inputs_fast_flavor @@ -6,6 +6,9 @@ maxError = 1e-6 # attenuation parameters to time derivative of N due to hamiltonians attenuation_hamiltonians = 1 +attenuation_vacuum = 1.0 +attenuation_matter = 1.0 +attenuation_self_interaction = 1.0 perturbation_type = 0 perturbation_amplitude = 0 @@ -24,11 +27,19 @@ nppc = (1, 1, 1) particle_data_filename = "particle_input.dat" # Maximum size of each grid in the domain -max_grid_size = 64 +max_grid_size_x = 1 +max_grid_size_y = 1 +max_grid_size_z = 1 # Number of steps to run nsteps = 100 +# Minimum timestep (s) +minimum_time_step = 1e-15 + +# Time step method +time_step_method = 0 + # Simulation end time end_time = 1.0 diff --git a/sample_inputs/inputs_fast_flavor_nonzerok b/sample_inputs/inputs_fast_flavor_nonzerok index a20c6f00..0df19d68 100644 --- a/sample_inputs/inputs_fast_flavor_nonzerok +++ b/sample_inputs/inputs_fast_flavor_nonzerok @@ -4,6 +4,9 @@ perturbation_wavelength_cm = 1 # attenuation parameters to time derivative of N due to hamiltonians attenuation_hamiltonians = 1 +attenuation_vacuum = 1.0 +attenuation_matter = 1.0 +attenuation_self_interaction = 1.0 collision_cfl_factor = 1e-3 cfl_factor = 0.5 @@ -25,11 +28,19 @@ nppc = (1, 1, 1) particle_data_filename = "particle_input.dat" # Maximum size of each grid in the domain -max_grid_size = 10 +max_grid_size_x = 1 +max_grid_size_y = 1 +max_grid_size_z = 50 # Number of steps to run nsteps = 5000 +# Minimum timestep (s) +minimum_time_step = 1e-15 + +# Time step method +time_step_method = 0 + # Simulation end time end_time = 1.0e-10 diff --git a/sample_inputs/inputs_fermi_dirac_test b/sample_inputs/inputs_fermi_dirac_test index f980f5f6..7fe4a83c 100644 --- a/sample_inputs/inputs_fermi_dirac_test +++ b/sample_inputs/inputs_fermi_dirac_test @@ -3,10 +3,13 @@ perturbation_amplitude = 1e-6 # attenuation parameters to time derivative of N due to hamiltonians attenuation_hamiltonians = 0 +attenuation_vacuum = 1.0 +attenuation_matter = 1.0 +attenuation_self_interaction = 1.0 -collision_cfl_factor = 100 -cfl_factor = 100 -flavor_cfl_factor = 100 +collision_cfl_factor = 55 +cfl_factor = 55 +flavor_cfl_factor = 0 max_adaptive_speedup = 0 maxError = 1e-6 @@ -24,10 +27,18 @@ nppc = (1, 1, 1) particle_data_filename = "particle_input.dat" # Maximum size of each grid in the domain -max_grid_size = 16 +max_grid_size_x = 2 +max_grid_size_y = 2 +max_grid_size_z = 2 # Number of steps to run -nsteps = 2000 +nsteps = 6000 + +# Minimum timestep (s) +minimum_time_step = 1e-15 + +# Time step method +time_step_method = 0 # Simulation end time end_time = 5.0e9 @@ -36,10 +47,10 @@ end_time = 5.0e9 amrex.fpe_trap_invalid=1 # Write plotfiles -write_plot_every = 2000 +write_plot_every = 6000 # Write particle data in plotfiles -write_plot_particles_every = 2000 +write_plot_particles_every = 6000 # checkpointing do_restart = 0 @@ -83,6 +94,7 @@ deltaCP_degrees = 222 IMFP_method = 2 Do_Pauli_blocking = 0 # If 1, it will multiply the inverse mean free path by 1 / (1 - f_eq); if 0, do nothing. boundary_condition_type = 0 +set_equilibrium_distribution = 0 ################# # Background Ye, T and rho diff --git a/sample_inputs/inputs_msw_test b/sample_inputs/inputs_msw_test index 8e837961..c0757439 100644 --- a/sample_inputs/inputs_msw_test +++ b/sample_inputs/inputs_msw_test @@ -6,6 +6,9 @@ maxError = 1e-6 # attenuation parameters to time derivative of N due to hamiltonians attenuation_hamiltonians = 1 +attenuation_vacuum = 1.0 +attenuation_matter = 1.0 +attenuation_self_interaction = 1.0 perturbation_type = 0 perturbation_amplitude = 0 @@ -24,11 +27,19 @@ nppc = (1, 1, 1) particle_data_filename = "particle_input.dat" # Maximum size of each grid in the domain -max_grid_size = 64 +max_grid_size_x = 1 +max_grid_size_y = 1 +max_grid_size_z = 1 # Number of steps to run nsteps = 50 +# Minimum timestep (s) +minimum_time_step = 1e-15 + +# Time step method +time_step_method = 0 + # Simulation end time end_time = 1.0