diff --git a/.github/workflows/krome.yml b/.github/workflows/krome.yml index 0e3c3c715..2331e83c7 100644 --- a/.github/workflows/krome.yml +++ b/.github/workflows/krome.yml @@ -7,7 +7,7 @@ env: PREFIX: /usr/local/ PHANTOM_DIR: ${{ github.workspace }} KROMEPATH: ${{ github.workspace }}/krome - FC: gfortran # default if not overwritten by matrix + PHANTOM_FC: gfortran # default if not overwritten by matrix # A workflow run is made up of one or more jobs that can run sequentially or in parallel jobs: @@ -31,7 +31,15 @@ jobs: - name: "Install LAPACK dependencies (Ubuntu)" if: matrix.os == 'ubuntu-latest' run: sudo apt-get update && sudo apt-get install -y liblapack-dev - + + - name: "Install HDF5 dependencies (Ubuntu)" + if: matrix.os == 'ubuntu-latest' + run: sudo apt-get update && sudo apt-get install -y libhdf5-serial-dev && echo "HDF5_DIR=/usr/lib/x86_64-linux-gnu/hdf5/serial" >> "$GITHUB_ENV" + + - name: "Install HDF5 dependencies (macOS)" + if: matrix.os == 'macos-latest' + run: brew update && brew install hdf5 && echo "HDF5_DIR=$(brew --prefix hdf5)" >> "$GITHUB_ENV" + - name: "Clone phantom" uses: actions/checkout@v4 @@ -42,13 +50,13 @@ jobs: run: echo -e "-n=networks/react_umist\n-iRHS\n-noSinkCheck\n-noRecCheck\n-noTlimits\n-unsafe\n-skipODEthermo\n-skipJacobian" > krome.setup && cat krome.setup - name: "Compile phantom and link with krome" - run: make SYSTEM=${{ env.FC }} SETUP=wind KROME=yes KROMEPATH=${KROMEPATH} PREFIX=${PREFIX} RUNDIR=${{ github.workspace }} + run: make SYSTEM=$PHANTOM_FC SETUP=wind KROME=yes RUNDIR=${{ github.workspace }} - name: "Compile phantomsetup and link with krome" - run: make SYSTEM=${{ env.FC }} SETUP=wind KROME=yes KROMEPATH=${KROMEPATH} PREFIX=${PREFIX} RUNDIR=${PREFIX} setup + run: make SYSTEM=$PHANTOM_FC SETUP=wind KROME=yes RUNDIR=${{ github.workspace }} setup - name: "Compile phantomanalysis and link with krome" - run: make SYSTEM=${{ env.FC }} SETUP=wind KROME=yes KROMEPATH=${KROMEPATH} PREFIX=${PREFIX} RUNDIR=${PREFIX} analysis + run: make SYSTEM=$PHANTOM_FC SETUP=wind KROME=yes RUNDIR=${{ github.workspace }} analysis - name: "Compile phantommoddump and link with krome" - run: make SYSTEM=${{ env.FC }} SETUP=wind KROME=yes KROMEPATH=${KROMEPATH} PREFIX=${PREFIX} RUNDIR=${PREFIX} moddump + run: make SYSTEM=$PHANTOM_FC SETUP=wind KROME=yes RUNDIR=${{ github.workspace }} moddump diff --git a/AUTHORS b/AUTHORS index 2b95c98b5..00582b408 100644 --- a/AUTHORS +++ b/AUTHORS @@ -35,7 +35,9 @@ Spencer Magnall Caitlyn Hardiman Amber Tilly Steven Rieder +Ali Pourmand Sergei Biriukov +Camille Landri Giovanni Dipierro Roberto Iaconi Siméon @@ -45,10 +47,8 @@ Stephen Nielson Sahl Rowther Martina Toscani Benedetta Veronesi -Camille Landri Simon Glover Thomas Reichardt -Ali Pourmand Chunliang Mu <86601204+chunliangmu@users.noreply.github.com> Jean-François Gonzalez Christopher Russell diff --git a/build/Makefile b/build/Makefile index c0c7b49e1..335123eec 100644 --- a/build/Makefile +++ b/build/Makefile @@ -266,6 +266,7 @@ ifeq ($(SINKTREE), yes) endif ifeq ($(KROME), yes) + HDF5 := yes FPPFLAGS += -DKROME LDFLAGS += -llapack endif @@ -404,7 +405,7 @@ ifeq (X$(HDF5LIB), X) endif FFLAGS+= -I$(HDF5INCLUDE) CCFLAGS+= -I$(HDF5INCLUDE) - LDFLAGS+= -L$(HDF5_DIR)/lib -lhdf5 -lhdf5_fortran + LDFLAGS+= -L$(HDF5LIB) -Wl,-rpath,$(HDF5LIB) -lhdf5 -lhdf5_fortran endif IDFLAGS=$(FPPFLAGS) @@ -884,7 +885,7 @@ allanalysis: analysis_structurefn.f90) \ analysisfile="utils_timing.f90 io_structurefn.f90 utils_sort.f90 random.f90 struct_part.f90 $$analysisfile"; \ ;; \ - analysis_krome.f90|analysis_mcfost.f90) \ + analysis_krome.F90|analysis_mcfost.f90) \ continue; \ ;; \ esac; \ diff --git a/docs/physics/composition.rst b/docs/physics/composition.rst index e395dde06..bf5fbd709 100644 --- a/docs/physics/composition.rst +++ b/docs/physics/composition.rst @@ -1,10 +1,6 @@ -Composition tracking in phantom +Composition tracking in Phantom ===================================== -This file documents how to track chemical abundances and mixing in phantom. - -Composition tracking with fixed abundances --------------------------------------------- Tracking the chemical composition of the gas in phantom with fixed abundances is straightforward, since phantom is a Lagrangian code and the particle identifiers are preserved throughout the simulation. @@ -14,27 +10,185 @@ always written to the dump files in the same order, so the particle id is simply the particle index in the dump file. In the MPI code, the particle id is stored in the 'iorig' array in the dump file. -Hence composition tracking can done as a post-processing step. +Hence composition tracking can be done as a post-processing step. + +Tracking chemical abundances with KROME +===================================== + +Phantom contains routines to post-process simulation outputs with KROME, a chemical kinetics package that solves time-dependent chemical evolution according to a specified network of chemical reactions. +KROME has only been robustly tested in Phantom for post-processing, and not for on-the-fly chemistry. +The publicly available routines of Phantom+KROME are meant to be used to process wind simulations, and are not guaranteed to be accurate for other types of simulations. +Below we provide a guide on how to install KROME, and run phantomanalysis to calculate the evolution of chemical species in a simulation. + +To calculate chemical abundances in Phantom dumps, you first need to install KROME and HDF5. + +Downloading KROME +~~~~~~~~~~~~~~~~~~ +You can download KROME from the publicly available repository (https://bitbucket.org/tgrassi/krome/src/master/). For an introduction on how KROME works, follow the tutorials provided in their documentation. +Make sure to store it in a directory that will be accessible to your Phantom installation. You do not need to compile KROME right now, it will be compiled automatically when you compile Phantom with KROME support. +You should set the environment variable ``KROMEPATH`` to point to the installation directory of KROME: +:: + export KROMEPATH=/path/to/krome + +Installing HDF5 +~~~~~~~~~~~~~~~~~~ +HDF5 is a library and file format for storing large amounts of data. The KROME routines within Phantom make use of the HDF5 library to read and write output files, because chemical abundance files are typically too large +to be stored in ASCII format. +There are several ways to install HDF5, you can get it directly from your package manager (on Linux or macOS), you can download precompiled binaries from the official release on the GitHub page (https://github.com/HDFGroup/hdf5/releases/tag/2.1.1), +or you can compile it from source (also using their GitHub repository). If you work on a cluster, it is likely to already be installed as a module, and you should be able to load the HDF5 module installed on the system (for questions, contact your cluster admin). + +To decide how to install HDF5, **you must check which compiler you intend to use for Phantom**. When linking libraries, it is important that they are compiled with the same compiler and version, otherwise you may encounter linking errors when compiling your code. +Thus, if you intend to compile Phantom with gfortran, you should install HDF5 with gfortran, same goes for ifort or ifx. Pre-compiled binaries available on the HDF5 GitHub releases or your own package manager are compiled with gfortran or ifx, but not ifort. +If you are using ifort, you will need to compile HDF5 from source. +Another constraint to account for is that **KROME is only parallelized with ifort and ifx**, so if you intend to use KROME in parallel, you will need to compile HDF5 with ifort or ifx as well. +You can technically compile KROME with gfortran, but this will strongly limit the performance of the post-processing routines. + +If you already have HDF5 installed, you can check which compiler was used to compile it by running the command ``h5fc -showconfig`` in the terminal. + +Below we detail the three installation methods for gfortran, ifort, and ifx, and more specific instructions for macOS users. + +gfortran +--------- +If you're planning on compiling Phantom with gfortran, you can simply install HDF5 from your package manager, for example on Ubuntu/Debian, with the following commands: +:: + # Ubuntu/Debian + sudo apt-get install libhdf5-serial-dev + +You should then set the environment variable ``HDF5_DIR`` to point to the installation directory of HDF5, which is typically ``/usr/lib/x86_64-linux-gnu/hdf5/serial``: +:: + export HDF5_DIR=/usr/lib/x86_64-linux-gnu/hdf5/serial +You can check where it is installed by running the command ``h5fc -showconfig`` in the terminal, which will show details about the installation of HDF5. + + +ifx +--------- +You can download ifx precompiled binaries from the official release on the GitHub page (https://github.com/HDFGroup/hdf5/releases/tag/2.1.1). +Make sure to download the file containing "intel" in the filename, for example ``hdf5-2.1.1-ubuntu-2404-intel.tar.gz``. +The downloaded archive can be extracted: +:: + tar -xvf hdf5-2.1.1-ubuntu-2404-intel.tar.gz + +which will create a directory containing another archive, which needs to be extracted as well: +:: + tar -xvf HDF5-2.1.1-Linux.tar.gz +The extracted directory contains a series of subdirectories, including a ``HDF5/``, which contains all the libraries and binaries needed to link HDF5 to Phantom. +You can move this directory to a permanent location that is more accessible, and set the environment variable ``HDF5_DIR`` to point to it: +:: + export HDF5_DIR=/path/to/HDF5 + +ifort +----- +Since no precompiled binaries are available for ifort, you will need to compile HDF5 from source. You can download the source code from the official release on the GitHub page (https://github.com/HDFGroup/hdf5/releases/tag/2.1.1). +Extract the downloaded archive to a permanent location and store the location in the environment variable ``HDF5_DIR``: +:: + tar -xvf hdf5-2.1.1.tar.gz + cd /path/to/hdf5-2.1.1 + export HDF5_DIR=/path/to/hdf5-2.1.1 +Then, you can set up the compilation of HDF5 with ifort using the following commands: +:: + cmake --preset ci-StdShar-Intel -B build/ifort \ + -DCMAKE_Fortran_COMPILER=ifort \ + -DCMAKE_INSTALL_PREFIX=$HDF5_DIR \ + -DHDF5_BUILD_FORTRAN=ON \ + -DHDF5_BUILD_JAVA=OFF \ + -DHDF5_ENABLE_JNI=OFF \ + -DHDF5_BUILD_CPP_LIB=OFF \ + -DHDF5_BUILD_TOOLS=OFF \ + -DHDF5_BUILD_EXAMPLES=OFF \ + -DBUILD_TESTING=OFF \ + -DHDF5_ENABLE_PLUGIN_SUPPORT=OFF +You can then compile and install HDF5 with: +:: + cmake --build build/ifort -j + cmake --install build/ifort +Finally, you need to specify two new environment variables to point to the HDF5 libraries and binaries: +:: + export HDF5INCLUDE=$HDF5_DIR/mod/shared + export HDF5LIB=$HDF5_DIR/lib + +macOS +----- +On macOS, you can only compile Phantom with gfortran, since Intel compilers are not supported anymore. +This means that you can only compile Phantom and KROME with gfortran, and **you can only use KROME in serial**. +You can install HDF5 with Homebrew: +:: + brew install hdf5 +Using the following command, you can store the path to HDF5 (likely something like ``/opt/homebrew/opt/hdf5``) and automatically store it to the environment variable ``HDF5_DIR``: +:: + export HDF5_DIR=$(brew --prefix hdf5) + + +Compiling Phantom with KROME +~~~~~~~~~~~~~~~~~~ +Once both KROME and HDF5 are installed, you can compile Phantom with KROME libraries. +You can use KROME post-processing on pre-existing models, or you can create one for the occasion. + +Once you have selected the model to process, you need to create a setup file for KROME in the directory where your phantom makefile is stored. +This file, called krome.setup, contains compilation options for KROME. For wind models, it can be created using the following command: +:: + echo -e "-n=networks/react_umist\n-iRHS\n-noSinkCheck\n-noRecCheck\n-noTlimits\n-unsafe\n-skipODEthermo\n-skipJacobian" > krome.setup + +Each flag is explained in the KROME documentation, but the most important one is ``-n``, which specifies the chemical network to use. +The chemical network files are stored in the ``networks/`` directory of the KROME directory, and you may choose any network, +although **there is no guarantee of the chemical network being suitable for your simulation**. +The network we are using here, ``react_umist`` is a detailed network for winds of AGB stars based on observations and detailed 1D chemical simulations, +but you may want to explore other networks depending on your needs. + +Once the setup file is created, you can compile Phantom with KROME by running: +:: + make KROME=yes + +The compilation will only work if the paths to HDF5 and KROME directories are set up properly, so remember to set the environment variables ``KROMEPATH`` and ``HDF5_DIR`` before compiling Phantom. +If you intend to use KROME with Phantom often, you may consider adding these variables to your ``~/.bashrc`` file. + +You can then compile the KROME routines for post-processing by running: +:: + make analysis KROME=yes ANALYSIS=analysis_krome.F90 + + +.. note:: + - You need to have compiled Phantom with KROME support before attempting to compile or use the KROME analysis routines. + - Forgetting to set KROME=yes will result in a regular compilation of Phantom **without** KROME support, and the analysis routines will not compile. + - If you need to compile other phantom utilities (phantomsetup, phantommoddump, etc.) while you are using KROME analysis, you need to compile them **with** KROME support as well, otherwise you will get linking errors when trying to run the different routines. + - Switching between KROME and non-KROME compilations of Phantom is possible, but you need to recompile Phantom with the appropriate flags each time you switch and **you need to clean the build directory** with ``make clean KROME=yes`` to avoid conflicts. + - When using ``make KROME=yes``, the compilation of Phantom takes longer than usual, since it compiles KROME as well, thus you may want to compile Phantom with KROME only when you need to use the KROME routines. + - The compilation of KROME is quiet, so you will not see much output from the compilation process itself, and since it takes a few minutes to complete it may look like it is stuck. + - If the compilation fails with missing routines or variables with "krome" in their name, KROME likely failed to compile properly. Run ``make clean KROME=yes``, check your environment variables then recompile Phantom with KROME. + + +Phantomanalysis with KROME +~~~~~~~~~~~~~~~~~~ +Your simulation is now ready to be post-processed with KROME the same way you would post-process it with any other analysis routine: +:: + ./phantomanalysis dump_????? -Writing a .comp file -~~~~~~~~~~~~~~~~~~~~~ -In practice, the composition of the gas can be tracked by writing a .comp file -either once for the entire simulation, or one per dump file. The .comp file -is a simple ascii file with one row per particle, where the columns are the abundances. -You should also give the first line a header with labels for each element:: +The routines will read the first dump file, if it finds a chemistry file corresponding to the dump (dump_?????.h5), it will read the chemical abundances from that file and use them as initial conditions for the chemical evolution, otherwise it will use default initial abundances. +The default initial abundances are set in ``analysis_krome.F90`` (located in ``src/utils/``), and correspond to AGB stellar winds abundances. - # h1, he3, he4, c12, n14, o16, ne20, mg24, si28, s32, ar36, ca40, ti44, cr48, fe52, ni56 - 7.1119142075E-01 9.3180507858E-05 2.7341075111E-01 ... - 7.1119141814E-01 9.3180815248E-05 2.7341075311E-01 ... - 7.0179434509E-01 3.6580185766E-06 2.8264795926E-01 ... - ... +In case you need to change the initial abundances, you can edit the subroutine ``set_initial_abundances`` in ``analysis_krome.F90``, recompile the analysis routines, and run ``phantomanalysis`` again. +Please copy the ``analysis_krome.F90`` file to a different name before editing it, so we can keep the original version intact as you make any changes you want to your copy. You then need to change how ``phantomanalysis`` is compiled: +:: + make analysis KROME=yes ANALYSIS=analysis_krome_mycopy.F90 -If the phantom dump is called `foo_0000`, the .comp file should be called `foo_0000.comp` -or simply `foo.comp` if the composition is constant on each particle. +If your analysis stopped before processing all the dumps, you can restart it by using a bit of shell scripting. You can store the following commands in a bash script and run it, or include them in the script you use to run ``phantomanalysis`` on your cluster: +:: + export FILE=dump_name + export DIR=/path/to/dumps -This file is automatically read by splash and used to create extra columns in the visualisation. + start=00010 #the last dump number that was processed + files=$(printf '%s\n' $DIR/"$FILE"_* \ + | awk -F_ -v s="$start" '$NF ~ /^[0-9]{5}$/ && ($NF+0)>=s' \ + | sort) -See also --------- + if [ -z "$files" ]; then + echo "No files found in $DIR matching pattern $FILE_* with number >= $start" + exit 1 + fi + ./phantomanalysis $files +This script will find all the dump files in the directory that match the pattern and have a number greater than or equal to the last processed dump (that you must indicate by hand), and run ``phantomanalysis`` on the list of files. -- :doc:`Setting up stars and tidal disruption events ` +.. WARNING:: + The KROME routines are not guaranteed to work for all types of simulations, and have only been robustly tested for post-processing wind simulations. + More importantly, no one but you can guarantee that the chemical network you are using is suitable for your simulation. + KROME developers provide a large number of chemical networks, but it is up to you to choose the one that is appropriate for your simulation, or build your own. + Please make sure to check the KROME documentation and tutorials to understand how to choose a suitable chemical network, and how to set up initial abundances and physical conditions for your simulation. diff --git a/src/main/eos.f90 b/src/main/eos.f90 index d75ad799b..01c4da0a0 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -50,8 +50,9 @@ module eos ! ! :Dependencies: dim, dump_utils, eos_HIIR, eos_barotropic, eos_gasradrec, ! eos_helmholtz, eos_idealplusrad, eos_mesa, eos_piecewise, eos_shen, -! eos_stamatellos, eos_stratified, eos_tillotson, infile_utils, io, -! ionization_mod, mesa_microphysics, part, physcon, units +! eos_stamatellos, eos_stratified, eos_tillotson, eos_zerotemp, +! infile_utils, io, ionization_mod, mesa_microphysics, part, physcon, +! units ! use part, only:ien_etotal,ien_entropy,ien_type use dim, only:gr,do_radiation @@ -529,7 +530,7 @@ subroutine equationofstate(eos_type,ponrhoi,spsoundi,rhoi,xi,yi,zi,tempi,eni,gam ponrhoi = presi/rhoi gammai = 1.d0 + presi/(eni*rhoi) spsoundi = sqrt(gammai*ponrhoi) - case (25) ! zero temperature EOS + case (25) ! zero temperature EOS cgsrhoi = rhoi * unit_density call get_zerotemp_pressure(cgsrhoi,cgspresi) @@ -538,7 +539,7 @@ subroutine equationofstate(eos_type,ponrhoi,spsoundi,rhoi,xi,yi,zi,tempi,eni,gam presi = cgspresi/unit_pressure ponrhoi = presi/rhoi spsoundi = cgsspsoundi / unit_velocity - tempi = 0. + tempi = 0. case default spsoundi = 0. ! avoids compiler warnings ponrhoi = 0. @@ -647,7 +648,6 @@ subroutine init_eos(eos_type,ierr) if (ierr_ra > 0) call warning('init_eos','Failed to read EOS file') call init_coolra - case(25) ! ! zero temperature @@ -989,7 +989,7 @@ subroutine calc_temp_and_ene(eos_type,rho,pres,ene,temp,ierr,guesseint,mu_local, case(24) ! Stamatellos temp = pres /(rho * Rg) * mu call getintenerg_opdep(temp, rho, ene) - case(25) ! zero temp eos + case(25) ! zero temp eos call get_zerotemp_u(rho,ene) temp = 0 case default @@ -1035,7 +1035,7 @@ subroutine calc_rho_from_PT(eos_type,pres,temp,rho,ierr,mu_local,X_local,Z_local rho = pres / (temp * Rg) * mu case(12) ! Ideal gas + radiation call get_idealplusrad_rhofrompresT(pres,temp,mu,rho) - case(25) ! zero temperature eos + case(25) ! zero temperature eos call get_zerotemp_rhofrompres(pres,rho,ierr) case default ierr = 1 @@ -1511,7 +1511,7 @@ logical function eos_allows_shock_and_work(ieos) integer, intent(in) :: ieos select case(ieos) - case(2,5,10,12,15,16,21,22,24,25) + case(2,5,10,12,15,16,21,22,24,25) eos_allows_shock_and_work = .true. case default eos_allows_shock_and_work = .false. @@ -1538,7 +1538,7 @@ end function eos_requires_polyk ! a non-zero pressure even if no thermal energy is set !+ !----------------------------------------------------------------------- -logical function eos_has_pressure_without_u(ieos) +logical function eos_has_pressure_without_u(ieos) integer, intent(in) :: ieos eos_has_pressure_without_u = eos_requires_isothermal(ieos) .or. & diff --git a/src/main/eos_zerotemp.f90 b/src/main/eos_zerotemp.f90 index 4c59d9cda..326ad7b19 100644 --- a/src/main/eos_zerotemp.f90 +++ b/src/main/eos_zerotemp.f90 @@ -6,18 +6,27 @@ !--------------------------------------------------------------------------! module eos_zerotemp ! -! Implements zero temperature equation of state, e.g. for white dwarfs. Meant to +! Implements zero temperature equation of state, e.g. for white dwarfs. Meant to ! ! :References: Kippenhahn & Weigert, Stellar Structure and Evolution, section 15.2 ! -! :Dependencies: physcon +! :Owner: Ali Pourmand +! +! :Runtime parameters: +! - xc : *Carbon mass fraction* +! - xh : *Hydrogen mass fraction* +! - xhe : *Helium mass fraction* +! - xmg : *Magnesium mass fraction* +! - xne : *Neon mass fraction* +! - xo : *Oxygen mass fraction* +! +! :Dependencies: infile_utils, io, physcon, units ! use units, only:unit_density,unit_velocity use physcon, only:pi,atomic_mass_unit,mass_electron_cgs,planckh,c implicit none real :: mu_e ! mean molecular weight per free electron - public :: eos_zerotemp_init,read_options_eos_zerotemp,write_options_eos_zerotemp,& f_chandra,get_zerotemp_pressure,get_zerotemp_u,get_zerotemp_spsoundi,& get_zerotemp_rhofrompres,eos_zerotemp_calc_mu_e,eos_zerotemp_eosinfo @@ -31,8 +40,8 @@ module eos_zerotemp real :: xc = 0.5 real :: xo = 0.5 real :: xne = 0.0 - real :: xmg = 0.0 - + real :: xmg = 0.0 + integer, parameter :: speciesmax = 6 character(len=10) :: speciesname(speciesmax) real :: xmass(speciesmax) ! mass fraction of species @@ -40,7 +49,6 @@ module eos_zerotemp real :: Zion(speciesmax) ! number of protons contains - !---------------------------------------------------------------- !+ ! initialise zero-temperature EOS composition @@ -48,7 +56,7 @@ module eos_zerotemp !---------------------------------------------------------------- subroutine eos_zerotemp_init(ierr) - use io, only: warning, fatal + use io, only:warning, fatal integer, intent(out) :: ierr ierr = 0 @@ -98,7 +106,7 @@ end subroutine eos_zerotemp_init !---------------------------------------------------------------- subroutine read_options_eos_zerotemp(db,nerr) - use infile_utils, only: inopts, read_inopt + use infile_utils, only:inopts, read_inopt type(inopts), intent(inout) :: db(:) integer, intent(inout) :: nerr @@ -118,7 +126,7 @@ end subroutine read_options_eos_zerotemp !---------------------------------------------------------------- subroutine write_options_eos_zerotemp(iunit) - use infile_utils, only: write_inopt + use infile_utils, only:write_inopt integer, intent(in) :: iunit call write_inopt(xh ,'xh' ,'Hydrogen mass fraction',iunit) @@ -171,8 +179,6 @@ real function f_chandra(x) result(fx) fx = x*(2*x**2 - 3)*sqrt(1+x**2) + 3*log(x + sqrt(1+x**2)) end function f_chandra - - ! ---------------------------------------------------------------- !+ ! Calculates the zero temperature pressure for a fully degenerate electron gas, i.e. a @@ -191,7 +197,6 @@ subroutine get_zerotemp_pressure(rhoi,presi) ! This is a good approximation for white dwarfs where the electrons are highly degenerate but the ions are not. ! Note also that this assumes a fully ionised gas, so mu is the mean molecular weight per free electron - ne = rhoi / (mu_e * atomic_mass_unit) x = ( (3.0 * ne * planckh**3) / & @@ -201,7 +206,6 @@ subroutine get_zerotemp_pressure(rhoi,presi) end subroutine get_zerotemp_pressure - !----------------------------------------------------------------------- !+! Inputs and outputs in cgs units !+! get internal energy from density for zero temperature EOS @@ -221,7 +225,6 @@ subroutine get_zerotemp_u(rhoi,u) u = u/rhoi !previous equation is erg/cm^3 so here I convert it end subroutine get_zerotemp_u - !---------------------------------------------------------------- !+ ! Calculates sound speed from density (derivative of f(x), analytically found), input output in cgs @@ -273,7 +276,7 @@ subroutine get_zerotemp_rhofrompres(presi,densi,ierr) do while (f_chandra(xhi) < fx) xhi = 2.0*xhi - end do + enddo ierr = 0 @@ -285,11 +288,11 @@ subroutine get_zerotemp_rhofrompres(presi,densi,ierr) xhi = xmid else xlo = xmid - end if + endif if (abs(xhi-xlo)/(xmid+1.e-300) < tolerance) exit - end do + enddo if (iter == iter_max) ierr = 1 @@ -302,5 +305,4 @@ subroutine get_zerotemp_rhofrompres(presi,densi,ierr) end subroutine get_zerotemp_rhofrompres - end module eos_zerotemp \ No newline at end of file diff --git a/src/main/ptmass_radiation.f90 b/src/main/ptmass_radiation.f90 index ef7cf2177..c31b09fcb 100644 --- a/src/main/ptmass_radiation.f90 +++ b/src/main/ptmass_radiation.f90 @@ -284,16 +284,20 @@ subroutine get_dust_temperature(npart,xyzh,eos_vars,nptmass,xyzmh_ptmass,dust_te if (iget_tdust == 4) then ! update tau_Lucy if (idust_opacity == 2) then - call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, nucleation(:,ikappa), iray_resolution, tau_lucy) + call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, nucleation(:,ikappa),& + iray_resolution, .false., tau_lucy) else - call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, calc_kappa_bowen(dust_temp(1:npart)), iray_resolution, tau_lucy) + call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, calc_kappa_bowen(dust_temp(1:npart)),& + iray_resolution, .false., tau_lucy) endif elseif (itau_alloc == 1) then ! update tau if (idust_opacity == 2) then - call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, nucleation(:,ikappa), iray_resolution, tau) + call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, nucleation(:,ikappa),& + iray_resolution, .true., tau) else - call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, calc_kappa_bowen(dust_temp(1:npart)), iray_resolution, tau) + call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, calc_kappa_bowen(dust_temp(1:npart)),& + iray_resolution, .true., tau) endif endif ! diff --git a/src/main/utils_raytracer.f90 b/src/main/utils_raytracer.f90 index 1a14eddf8..0f1e3a45e 100644 --- a/src/main/utils_raytracer.f90 +++ b/src/main/utils_raytracer.f90 @@ -20,12 +20,12 @@ module raytracer ! ! :Runtime parameters: None ! -! :Dependencies: dim, healpix, kernel, neighkdtree, part, units +! :Dependencies: healpix, kernel, neighkdtree, part, units ! use healpix implicit none - public :: get_all_tau + public :: get_all_tau, get_all_tau_single private @@ -43,24 +43,26 @@ module raytracer ! IN: xyzh: The array containing the particles position+smooting lenght ! IN: kappa_cgs: The array containing the opacities of all SPH particles ! IN: order: The healpix order which is used for the uniform ray sampling + ! IN: outwards: The direction of the ray, true for outwards false for inwards !+ ! OUT: tau: The array of optical depths for each SPH particle !+ !------------------------------------------------------------------------------------ -subroutine get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, kappa_cgs, order, tau) +subroutine get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, kappa_cgs, order, outwards, tau) use part, only: iReff - integer, intent(in) :: npart, order, nptmass - real, intent(in) :: kappa_cgs(:), xyzh(:,:), xyzmh_ptmass(:,:) - real, intent(out) :: tau(:) + integer, intent(in) :: npart, order, nptmass + real, intent(in) :: kappa_cgs(:), xyzh(:,:), xyzmh_ptmass(:,:) + logical, intent(in) :: outwards + real, intent(out) :: tau(:) real :: Rinject Rinject = xyzmh_ptmass(iReff,1) if (nptmass == 2 ) then call get_all_tau_companion(npart, xyzmh_ptmass(1:3,1), xyzmh_ptmass(iReff,1), xyzh, kappa_cgs, & - Rinject, xyzmh_ptmass(1:3,2), xyzmh_ptmass(iReff,2), order, tau) + Rinject, xyzmh_ptmass(1:3,2), xyzmh_ptmass(iReff,2), order, outwards, tau) else call get_all_tau_single(npart, xyzmh_ptmass(1:3,1), xyzmh_ptmass(iReff,1), xyzh,& - kappa_cgs, Rinject, order, tau) + kappa_cgs, Rinject, order, outwards, tau) endif end subroutine get_all_tau @@ -78,20 +80,22 @@ end subroutine get_all_tau ! IN: Rstar: The radius of the primary star ! IN: Rinject: The particles injection radius ! IN: order: The healpix order which is used for the uniform ray sampling + ! IN: outwards: The direction of the ray, true for outwards false for inwards !+ ! OUT: taus: The array of optical depths to each SPH particle !+ !--------------------------------------------------------------------------------- -subroutine get_all_tau_single(npart, primary, Rstar, xyzh, kappa, Rinject, order, tau) +subroutine get_all_tau_single(npart, primary, Rstar, xyzh, kappa, Rinject, order, outwards, tau) use part, only:isdead_or_accreted - integer, intent(in) :: npart,order - real, intent(in) :: primary(3), kappa(:), Rstar, Rinject, xyzh(:,:) - real, intent(out) :: tau(:) + integer, intent(in) :: npart,order + real, intent(in) :: primary(3), Rstar, xyzh(:,:), kappa(:), Rinject + logical, intent(in) :: outwards + real, intent(out) :: tau(:) integer :: i, nrays, nsides real :: ray_dir(3),part_dir(3) - real, allocatable :: rays_dist(:,:), rays_tau(:,:) - integer, allocatable :: rays_dim(:) + real, dimension(:,:), allocatable :: rays_dist, rays_tau + integer, dimension(:), allocatable :: rays_dim integer, parameter :: ndim = 200 ! maximum number of points along the ray where tau is calculated nrays = 12*4**order ! The number of rays traced given the healpix order @@ -108,13 +112,13 @@ subroutine get_all_tau_single(npart, primary, Rstar, xyzh, kappa, Rinject, order !$omp parallel default(none) & !$omp private(ray_dir) & - !$omp shared(nrays,nsides,primary,kappa,xyzh,Rstar,Rinject,rays_dist,rays_tau,rays_dim) + !$omp shared(nrays,nsides,primary,kappa,xyzh,Rstar,Rinject,rays_dist,rays_tau,rays_dim,outwards) !$omp do do i = 1, nrays !returns ray_dir, the unit vector identifying a ray (index i-1 because healpix starts counting from index 0) call pix2vec_nest(nsides, i-1, ray_dir) !calculate the properties along the ray (tau, distance, number of points) - call ray_tracer(primary,ray_dir,xyzh,kappa,Rstar,Rinject,rays_tau(:,i),rays_dist(:,i),rays_dim(i)) + call ray_tracer(primary,ray_dir,xyzh,kappa,Rstar,Rinject,outwards,rays_tau(:,i),rays_dist(:,i),rays_dim(i)) enddo !$omp enddo !$omp end parallel @@ -157,21 +161,23 @@ end subroutine get_all_tau_single ! IN: companion: The xyz coordinates of the companion ! IN: Rcomp: The radius of the companion ! IN: order: The healpix order which is used for the uniform ray sampling + ! IN: outwards: The direction of the ray, true for outwards false for inwards !+ ! OUT: tau: The array of optical depths for each SPH particle !+ !-------------------------------------------------------------------------- -subroutine get_all_tau_companion(npart, primary, Rstar, xyzh, kappa, Rinject, companion, Rcomp, order, tau) +subroutine get_all_tau_companion(npart, primary, Rstar, xyzh, kappa, Rinject, companion, Rcomp, order, outwards, tau) use part, only:isdead_or_accreted - integer, intent(in) :: npart, order - real, intent(in) :: primary(3), companion(3), kappa(:), Rstar, Rinject, xyzh(:,:), Rcomp - real, intent(out) :: tau(:) + integer, intent(in) :: npart, order + real, intent(in) :: primary(3), companion(3), kappa(:), Rstar, Rinject, xyzh(:,:), Rcomp + logical, intent(in) :: outwards + real, intent(out) :: tau(:) integer :: i, nrays, nsides real :: normCompanion,theta0,phi,cosphi,sinphi,theta,sep,root real :: ray_dir(3),part_dir(3),uvecCompanion(3) - real, allocatable :: rays_dist(:,:), rays_tau(:,:) - integer, allocatable :: rays_dim(:) + real, dimension(:,:), allocatable :: rays_dist, rays_tau + integer, dimension(:), allocatable :: rays_dim integer, parameter :: ndim = 200 ! maximum number of points along the ray where tau is calculated nrays = 12*4**order ! The number of rays traced given the healpix order @@ -197,7 +203,7 @@ subroutine get_all_tau_companion(npart, primary, Rstar, xyzh, kappa, Rinject, co !$omp parallel default(none) & !$omp private(ray_dir,theta,root,sep) & !$omp shared(nrays,nsides,primary,kappa,xyzh,Rstar,Rinject,Rcomp,rays_dist,rays_tau,rays_dim) & - !$omp shared(uvecCompanion,normCompanion,cosphi,sinphi,theta0) + !$omp shared(uvecCompanion,normCompanion,cosphi,sinphi,theta0,outwards) !$omp do do i = 1, nrays !returns ray_dir, the unit vector identifying a ray (index i-1 because healpix starts counting from index 0) @@ -210,9 +216,9 @@ subroutine get_all_tau_companion(npart, primary, Rstar, xyzh, kappa, Rinject, co if (theta < theta0) then root = sqrt(Rcomp**2-normCompanion**2*sin(theta)**2) sep = normCompanion*cos(theta)-root - call ray_tracer(primary,ray_dir,xyzh,kappa,Rstar,Rinject,rays_tau(:,i),rays_dist(:,i),rays_dim(i), sep) + call ray_tracer(primary,ray_dir,xyzh,kappa,Rstar,Rinject,outwards,rays_tau(:,i),rays_dist(:,i),rays_dim(i), sep) else - call ray_tracer(primary,ray_dir,xyzh,kappa,Rstar,Rinject,rays_tau(:,i),rays_dist(:,i),rays_dim(i)) + call ray_tracer(primary,ray_dir,xyzh,kappa,Rstar,Rinject,outwards,rays_tau(:,i),rays_dist(:,i),rays_dim(i)) endif enddo !$omp enddo @@ -262,9 +268,9 @@ end subroutine get_all_tau_companion !+ !-------------------------------------------------------------------------- subroutine interpolate_tau(nsides, vec, rays_tau, rays_dist, rays_dim, tau) - integer, intent(in) :: nsides, rays_dim(:) - real, intent(in) :: vec(:), rays_dist(:,:), rays_tau(:,:) - real, intent(out) :: tau + integer, intent(in) :: nsides, rays_dim(:) + real, intent(in) :: vec(:), rays_dist(:,:), rays_tau(:,:) + real, intent(out) :: tau integer :: rayIndex, neighbours(8), nneigh, i, k real :: tautemp, ray(3), vectemp(3), weight, tempdist(8), distRay_sq, vec_norm2 @@ -327,9 +333,9 @@ end subroutine interpolate_tau !+ !-------------------------------------------------------------------------- subroutine get_tau_on_ray(distance, tau_along_ray, dist_along_ray, len, tau) - real, intent(in) :: distance, tau_along_ray(:), dist_along_ray(:) - integer, intent(in) :: len - real, intent(out) :: tau + real, intent(in) :: distance, tau_along_ray(:), dist_along_ray(:) + integer, intent(in) :: len + real, intent(out) :: tau integer :: L, R, m ! left, right and middle index for binary search @@ -366,6 +372,7 @@ end subroutine get_tau_on_ray ! IN: kappa: The array containing the particles opacity ! IN: Rstar: The radius of the primary star ! IN: Rinject: The particles injection radius + ! IN: outwards: The direction of the ray, true for outwards false for inwards !+ ! OUT: tau_along_ray: The vector of cumulative optical depth along the ray ! OUT: dist_along_ray: The vector of distances from the primary along the ray @@ -374,13 +381,14 @@ end subroutine get_tau_on_ray ! OPT: maxDistance: The maximal distance the ray needs to be traced !+ !-------------------------------------------------------------------------- -subroutine ray_tracer(primary, ray, xyzh, kappa, Rstar, Rinject, tau_along_ray, dist_along_ray, len, maxDistance) +subroutine ray_tracer(primary, ray, xyzh, kappa, Rstar, Rinject, outwards, tau_along_ray, dist_along_ray, len, maxDistance) use units, only:unit_opacity use part, only:itauL_alloc - real, intent(in) :: primary(3), ray(3), Rstar, Rinject, xyzh(:,:), kappa(:) - real, intent(out) :: dist_along_ray(:), tau_along_ray(:) - integer, intent(out) :: len + real, intent(in) :: primary(3), ray(3), Rstar, Rinject, xyzh(:,:), kappa(:) + logical, optional :: outwards real, optional :: maxDistance + real, intent(out) :: dist_along_ray(:), tau_along_ray(:) + integer, intent(out) :: len real, parameter :: tau_max = 99. real :: dr, next_dr, h, dtaudr, previousdtaudr, nextdtaudr, distance @@ -418,9 +426,11 @@ subroutine ray_tracer(primary, ray, xyzh, kappa, Rstar, Rinject, tau_along_ray, endif len = i - if (itauL_alloc > 0) then - !reverse integration start from zero inward + if (.not. outwards) then + !reverse integration start from zero outward and integrate inwards tau_along_ray(1:len) = tau_along_ray(len) - tau_along_ray(1:len) + endif + if (itauL_alloc > 0) then !find the first point where tau_lucy < 2/3 if (tau_along_ray(1) > 2./3.) then L = 1 @@ -442,7 +452,7 @@ end subroutine ray_tracer logical function hasNext(inext, tau, distance, maxDistance) integer, intent(in) :: inext - real, intent(in) :: tau, distance + real, intent(in) :: tau, distance real, optional :: maxDistance real :: tau_max = 99. if (present(maxDistance)) then @@ -474,17 +484,16 @@ end function hasNext !-------------------------------------------------------------------------- subroutine find_next(inpoint, h, ray, xyzh, kappa, dtaudr, distance, inext) use neighkdtree, only:getneigh_pos,leaf_is_active,listneigh - use kernel, only:radkern,cnormk,wkern - use part, only:hfact,rhoh,massoftype,igas - use dim, only:maxpsph + use kernel, only:radkern,cnormk,wkern + use part, only:hfact,rhoh,massoftype,igas real, intent(in) :: xyzh(:,:), kappa(:), inpoint(:), ray(:), h integer, intent(inout) :: inext real, intent(out) :: distance, dtaudr integer, parameter :: nmaxcache = 0 - real :: xyzcache(0,nmaxcache) + real :: xyzcache(3, max(1,nmaxcache)) - integer :: nneigh, i, prev,j + integer :: nneigh, i, prev real :: dmin, vec(3), dr, raydistance, q, norm_sq prev = inext @@ -498,23 +507,21 @@ subroutine find_next(inpoint, h, ray, xyzh, kappa, dtaudr, distance, inext) dmin = huge(0.) !loop over all neighbours do i=1,nneigh - j = listneigh(i) - if (j > maxpsph) cycle - vec = xyzh(1:3,j) - inpoint + vec = xyzh(1:3,listneigh(i)) - inpoint norm_sq = dot_product(vec,vec) - q = sqrt(norm_sq)/xyzh(4,j) + q = sqrt(norm_sq)/xyzh(4,listneigh(i)) !add optical depth contribution from each particle - dtaudr = dtaudr+wkern(q*q,q)*kappa(j)*rhoh(xyzh(4,j), massoftype(igas)) + dtaudr = dtaudr+wkern(q*q,q)*kappa(listneigh(i))*rhoh(xyzh(4,listneigh(i)), massoftype(igas)) ! find the next particle : among the neighbours find the particle located the closest to the ray - if (j /= prev) then + if (listneigh(i) /= prev) then dr = dot_product(vec,ray) !projected distance along the ray if (dr>0.) then !distance perpendicular to the ray direction raydistance = norm_sq - dr**2 if (raydistance < dmin) then dmin = raydistance - inext = j + inext = listneigh(i) distance = dr endif endif diff --git a/src/setup/setup_grtde.f90 b/src/setup/setup_grtde.f90 index fb5482a02..7074af903 100644 --- a/src/setup/setup_grtde.f90 +++ b/src/setup/setup_grtde.f90 @@ -38,9 +38,9 @@ module setup ! - z2 : *pos z star 2* ! ! :Dependencies: eos, externalforces, gravwaveutils, infile_utils, io, -! kernel, mpidomain, options, orbits, part, physcon, relaxstar, -! setbinary, setorbit, setstar, setunits, setup_params, systemutils, -! timestep, units, vectorutils +! kernel, metric_tools, mpidomain, options, orbits, part, physcon, +! relaxstar, setbinary, setorbit, setstar, setunits, setup_params, +! systemutils, timestep, units, vectorutils ! use setstar, only:star_t diff --git a/src/tests/test_setstar.f90 b/src/tests/test_setstar.f90 index 89fc47f69..d5cbada09 100644 --- a/src/tests/test_setstar.f90 +++ b/src/tests/test_setstar.f90 @@ -52,7 +52,6 @@ subroutine test_setstar(ntests,npass) ! test red supergiant call test_redsupergiant(ntests,npass) - ! test white dwarf call test_whitedwarf(ntests,npass) @@ -364,11 +363,9 @@ subroutine test_redsupergiant(ntests,npass) end subroutine test_redsupergiant - - !----------------------------------------------------------------------- !+ -! test that we can successfully relax a white dwarf with the +! test that we can successfully relax a white dwarf with the ! zero temperature EoS and Helmholtz eos(Ali Pourmand, similar to test_redsupergiant) !+ !----------------------------------------------------------------------- @@ -448,14 +445,13 @@ subroutine test_whitedwarf(ntests,npass) rhozero=rhozero,npart_total=ntot,mask=i_belong,ierr=ierr,& write_files=.false.,density_error=rmserr,energy_error=ekin) - call checkval(ierr,expected_error,0,nfail(1),'set_star runs with expected ierr') call update_test_scores(ntests,nfail,npass) call check_setup(nerror,nwarn,restart=.false.) call checkval(nerror+nwarn,0,0,nfail(1),'no errors or warnings') call update_test_scores(ntests,nfail,npass) - + call checkval(rmserr,0.0,0.05,nfail(1),'error in density profile') call update_test_scores(ntests,nfail,npass) @@ -478,6 +474,4 @@ subroutine test_whitedwarf(ntests,npass) end subroutine test_whitedwarf - - end module testsetstar diff --git a/src/utils/analysis_krome.F90 b/src/utils/analysis_krome.F90 new file mode 100644 index 000000000..1b48408a1 --- /dev/null +++ b/src/utils/analysis_krome.F90 @@ -0,0 +1,407 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2026 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.github.io/ ! +!--------------------------------------------------------------------------! +module analysis +! +! Analysis routine comparing time between dumps +! +! :References: None +! +! :Owner: Camille Landri +! +! :Runtime parameters: None +! +! :Dependencies: eos, hdf5, io, krome_main, krome_user, neighkdtree, +! omp_lib, part, physcon, raytracer, units +! + use krome_user, only:krome_nmols + use part, only: maxp + use raytracer, only: get_all_tau + use hdf5 + #ifdef _OPENMP + use omp_lib, only: omp_set_num_threads, omp_get_max_threads, omp_get_wtime + #endif + implicit none + character(len=20), parameter, public :: analysistype = 'krome' + public :: do_analysis + logical, allocatable :: mask(:) + + real(8), allocatable :: abundance(:,:), abundance_prev(:,:), one(:) + character(len=16) :: abundance_label(krome_nmols) + integer(8), allocatable :: iorig_old(:) + integer, allocatable :: iprev(:) + logical :: done_init = .false. + real :: AuvAv = 4.65, albedo = 0.5 + +contains + +subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) + use part, only: isdead_or_accreted, iorig, rhoh, nptmass, xyzmh_ptmass, iReff, iboundary, igas, iphase, iamtype + use neighkdtree, only:build_tree + use units, only: utime,unit_density,udist + use physcon, only: atomic_mass_unit + use eos, only: get_temperature, ieos, gamma,gmw, init_eos + use io, only: fatal, iverbose + use krome_main, only: krome_init, krome + use krome_user, only: krome_get_names,krome_set_user_Auv,krome_set_user_xi,& + krome_set_user_alb,krome_set_user_AuvAv + character(len=*), intent(in) :: dumpfile + integer, intent(in) :: num,npart,iunit + real, intent(in) :: xyzh(:,:),vxyzu(:,:), xzymh_ptmass_copy(:,:) + real, intent(in) :: particlemass,time + real, save :: tprev = 0. + integer, save :: nprev = 0 + real :: dt_cgs, rho_cgs + real :: numberdensity, T_gas, gammai, mui, AUV, xi + real :: abundance_part(krome_nmols), Y(krome_nmols), column_density(npart), xyzh_copy(4,npart) + real :: max_radius, radius, tstart + integer :: i, j, k, isize=0, ierr, completed_iterations, npart_copy = 0, hdferr, i_radius = 1 + + +#ifdef _OPENMP +#ifdef __GFORTRAN__ + print*, "Setting number of threads to 1 (KROME is not thread-safe when compiled with gfortran)" + call omp_set_num_threads(1) +#else + print*, "running with ", omp_get_max_threads(), " threads" +#endif +#else + print*, "running without OpenMP" +#endif + + if (.not.done_init) then + done_init = .true. + + ! initialize HDF5 + call h5open_f(hdferr) + if (hdferr /= 0) then + call fatal(analysistype, 'Failed to initialize HDF5') + endif + + print*, "initialising KROME" + call krome_init() + print*, "Initialised KROME" + + abundance_label(:) = krome_get_names() + allocate(abundance(krome_nmols,maxp)) + abundance = 0. + allocate(abundance_prev(krome_nmols,maxp)) + abundance_prev = 0. + allocate(one(maxp)) + one = 1. + allocate(iorig_old(maxp)) + iorig_old = 0 + allocate(iprev(maxp)) + iprev = 0 + + !check if a file with abundances already exists, if so read it in and use it as initial abundances + inquire(file=trim(dumpfile)//'.h5', size=isize) + if (isize > 0) then + print*, "Found existing abundance file, reading in abundances" + call read_chem(npart, dumpfile) + else + print*, "setting abundances" + !$omp parallel do default(none) & + !$omp shared(npart,xyzh,vxyzu,dt_cgs,nprev,iorig,iorig_old,iprev) & + !$omp shared(abundance,abundance_prev,particlemass,unit_density) & + !$omp shared(ieos,rho_cgs,T_gas,j) & + !$omp private(i,abundance_part) + do i=1, npart + if (.not.isdead_or_accreted(xyzh(4,i))) then + call chem_init(abundance_part) + abundance(:,i) = abundance_part + endif + enddo + call write_chem(npart, dumpfile) + endif + call init_eos(ieos, ierr) + if (ierr /= 0) call fatal(analysistype, "Failed to initialise EOS") + + else + dt_cgs = (time - tprev)*utime + completed_iterations = 0 + print*, "not first step data, timestep = ",dt_cgs, "npart = ",npart, "nprev = ",nprev + print*, "Building neighbour tree..." +#ifdef _OPENMP + tstart = omp_get_wtime() +#endif + xyzmh_ptmass_copy(:,:) = xyzmh_ptmass(:,:) !to avoid overwriting the original ptmass array + xyzmh_ptmass_copy(iReff,1) = 2. + npart_copy = npart + xyzh_copy = xyzh(:,:npart) + call build_tree(npart_copy,npart_copy,xyzh_copy,vxyzu,xyzmh_ptmass_copy) +#ifdef _OPENMP + print*, " - Took ", omp_get_wtime() - tstart, " seconds" +#endif + print*, "Calculating column density..." + tstart = omp_get_wtime() + call get_all_tau(npart, nptmass, xyzmh_ptmass_copy, xyzh, one, 5, .false., column_density) + max_radius = 0.0 + do i = 1, npart + if (.not.isdead_or_accreted(xyzh(4, i))) then + radius = sqrt(xyzh(1, i)**2 + xyzh(2, i)**2 + xyzh(3, i)**2) + if (radius > max_radius) then + max_radius = radius + i_radius = i + endif + endif + enddo + column_density = column_density + rhoh(xyzh(4,i_radius),particlemass)*unit_density * max_radius * udist + print*, " - Took ", omp_get_wtime() - tstart, " seconds" + + print*, "Running KROME" +#ifdef _OPENMP + tstart = omp_get_wtime() +#endif + !$omp parallel do default(none) & + !$omp shared(npart,xyzh,vxyzu,dt_cgs,nprev,iorig,iorig_old,iprev,iverbose) & + !$omp shared(abundance,abundance_prev,particlemass,unit_density,udist,iphase) & + !$omp shared(ieos,gamma,gmw,time,completed_iterations,column_density,AuvAv,albedo) & + !$omp private(i,j,abundance_part,Y,rho_cgs,numberdensity,T_gas,gammai,mui,AUV,xi) + + outer: do i=1,npart + if (.not.isdead_or_accreted(xyzh(4,i))) then + if (i <= nprev) then + if (iorig(i) == iorig_old(i)) then + iprev(i) = i + j = i + endif + endif + if (iprev(i) == 0) then + inner: do k=1,nprev + if (iorig(i) == iorig_old(k)) then + iprev(i) = k + j = k + exit inner + endif + enddo inner + endif + + if (j == iprev(i)) then + abundance_part(:) = abundance_prev(:,iprev(i)) + else + call chem_init(abundance_part) + endif + if (iamtype(iphase(i)) /= iboundary .and. i > 2460) then ! 2460 is the amount of boundary particles + !Thermodynamic quantities + rho_cgs = rhoh(xyzh(4,i),particlemass)*unit_density + gammai = gamma + mui = gmw + numberdensity = rho_cgs / (mui * atomic_mass_unit) + T_gas = get_temperature(ieos,xyzh(1:3, i),rhoh(xyzh(4,i),particlemass),vxyzu(:,i),gammai,mui) + T_gas = max(T_gas,20.0d0) + + !Radiation quantities + AUV = AuvAv * column_density(i) / (mui * atomic_mass_unit) / 1.87e21 + xi = get_xi(AUV) + + call krome_set_user_Auv(AUV) + call krome_set_user_xi(xi) + call krome_set_user_alb(ALBEDO) + call krome_set_user_AuvAv(AuvAv) + + Y = abundance_part*numberdensity + call krome(Y,T_gas,dt_cgs) + abundance_part = Y/numberdensity + abundance(:,i) = abundance_part + endif + endif + if (iverbose > 1) then + !$omp atomic + completed_iterations = completed_iterations + 1 + print*, 'Completed ', completed_iterations, ' of ', npart + endif + enddo outer +#ifdef _OPENMP + print*, " - Took ", omp_get_wtime() - tstart, " seconds" + tstart = omp_get_wtime() +#endif + call write_chem(npart, dumpfile) +#ifdef _OPENMP + print*, " - Took ", omp_get_wtime() - tstart, " seconds" +#endif + endif + + ! keep track of previous abundances for next dump + nprev = npart + tprev = time + iorig_old(1:npart) = iorig(1:npart) + abundance_prev(:,1:npart) = abundance(:,1:npart) +end subroutine do_analysis + +real function get_xi(AUV) + use physcon, only:pi + real, intent(in) :: AUV + real :: xi + real :: W(6), GA(6), ceta + integer :: i + + W(1) = 0.17132449 + W(2) = 0.36076157 + W(3) = 0.46791393 + W(4) = W(1) + W(5) = W(2) + W(6) = W(3) + GA(1) = 0.93246951 + GA(2) = 0.66120939 + GA(3) = 0.23861919 + GA(4) = -GA(1) + GA(5) = -GA(2) + GA(6) = -GA(3) + + xi = 0.0 + do i=1,6 + ceta = (pi*GA(i)+pi)/2.0 + xi=xi+(W(i)*(sin(ceta)*exp((-AUV*ceta)/sin(ceta)))) + enddo + xi = (pi/4.0)*xi + + get_xi = xi + +end function get_xi + +subroutine write_chem(npart, dumpfile) + integer, intent(in) :: npart + character(len=*), intent(in) :: dumpfile + integer :: hdferr + integer(hid_t) :: file_id, dspace_id, dset_id, group_id, type_id + integer(hsize_t) :: dims(2) + integer(hsize_t) :: dims_labels(1) + integer(size_t) :: str_len + + ! create HDF5 file + print "(1x,a)",'Writing to '//trim(dumpfile)//'.h5' + call h5fcreate_f(trim(dumpfile)//'.h5', H5F_ACC_TRUNC_F, file_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to create HDF5 file' + return + endif + + ! create group for particle data + call h5gcreate_f(file_id, 'chemistry', group_id, hdferr) + ! create dataspace for particle datasets + dims(1) = krome_nmols + dims(2) = npart + call h5screate_simple_f(2, dims, dspace_id, hdferr) + ! create one 2D dataset for all species abundances: + call h5dcreate_f(group_id, 'abundances', H5T_NATIVE_DOUBLE, dspace_id, dset_id, hdferr) + call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, abundance(:,1:npart), dims, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to write dataset abundances' + return + endif + call h5dclose_f(dset_id, hdferr) + + ! store species labels so abundances row i is always identifiable + dims_labels(1) = krome_nmols + call h5screate_simple_f(1, dims_labels, dspace_id, hdferr) + call h5tcopy_f(H5T_FORTRAN_S1, type_id, hdferr) + str_len = len(abundance_label(1)) + call h5tset_size_f(type_id, str_len, hdferr) + call h5dcreate_f(group_id, 'species_labels', type_id, dspace_id, dset_id, hdferr) + call h5dwrite_f(dset_id, type_id, abundance_label, dims_labels, hdferr) + + call h5dclose_f(dset_id, hdferr) + call h5tclose_f(type_id, hdferr) + call h5sclose_f(dspace_id, hdferr) + call h5gclose_f(group_id, hdferr) + call h5fclose_f(file_id, hdferr) +end subroutine write_chem + +subroutine read_chem(npart, dumpfile) + integer, intent(in) :: npart + character(len=*), intent(in) :: dumpfile + integer :: hdferr + integer :: i + integer(hid_t) :: file_id, dset_id, group_id, filespace_id, type_id + integer(hsize_t) :: max_dims(2), file_dims(2) + integer(hsize_t) :: max_dims_labels(1), file_dims_labels(1) + integer(size_t) :: str_len + character(len=16) :: labels_file(krome_nmols) + + ! open HDF5 file + print "(1x,a)",'Reading from '//trim(dumpfile)//'.h5' + call h5fopen_f(trim(dumpfile)//'.h5', H5F_ACC_RDONLY_F, file_id, hdferr) + if (hdferr /= 0) then + print*,'ERROR: Failed to open HDF5 file' + return + endif + + ! open group for particle data + call h5gopen_f(file_id, 'chemistry', group_id, hdferr) + + call h5dget_space_f(dset_id, filespace_id, hdferr) + call h5sget_simple_extent_dims_f(filespace_id, file_dims, max_dims, hdferr) + if (file_dims(1) /= krome_nmols .or. file_dims(2) /= npart) then + print*,'ERROR: abundances shape mismatch in HDF5 file' + return + endif + + call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, abundance(:,1:npart), file_dims, hdferr) + call h5dclose_f(dset_id, hdferr) + call h5sclose_f(filespace_id, hdferr) + + ! read and validate species labels to guarantee row-label mapping + call h5dopen_f(group_id, 'species_labels', dset_id, hdferr) + call h5dget_space_f(dset_id, filespace_id, hdferr) + call h5sget_simple_extent_dims_f(filespace_id, file_dims_labels, max_dims_labels, hdferr) + if (file_dims_labels(1) /= krome_nmols) then + call fatal(analysistype, 'species_labels length mismatch in HDF5 file') + endif + + call h5tcopy_f(H5T_FORTRAN_S1, type_id, hdferr) + str_len = len(labels_file(1)) + call h5tset_size_f(type_id, str_len, hdferr) + call h5dread_f(dset_id, type_id, labels_file, file_dims_labels, hdferr) + + do i = 1, krome_nmols + if (trim(adjustl(labels_file(i))) /= trim(adjustl(abundance_label(i)))) then + print*,'ERROR: species_labels mismatch at index ', i + call fatal(analysistype, 'species_labels mismatch in HDF5 file') + endif + enddo + + call h5dclose_f(dset_id, hdferr) + call h5sclose_f(filespace_id, hdferr) + call h5tclose_f(type_id, hdferr) + call h5gclose_f(group_id, hdferr) + call h5fclose_f(file_id, hdferr) +end subroutine read_chem + +subroutine chem_init(abundance_part) + use krome_user, only:krome_idx_H2,krome_idx_He,krome_idx_CO,krome_idx_C2H2,& + krome_idx_HCN,krome_idx_N2,krome_idx_SiC2,krome_idx_CS,& + krome_idx_SiS,krome_idx_SiO,krome_idx_CH4,krome_idx_H2O,& + krome_idx_HCl,krome_idx_C2H4,krome_idx_NH3,krome_idx_HCP,& + krome_idx_HF,krome_idx_H2S,krome_idx_e,krome_get_electrons + real, intent(out) :: abundance_part(krome_nmols) + + ! Initial abundances for the krome model taken from Agúndez et al. (2020) + ! H2, He, CO, C2H2, HCN, N2, SiC2, CS, SiS, SiO, CH4, H2O, HCl, C2H4, NH3, HCP, HF, H2S + abundance_part(:) = 0. + abundance_part(krome_idx_H2) = 0.5d0 + abundance_part(krome_idx_He) = 8.5d-2 + abundance_part(krome_idx_CO) = 4d-4 + abundance_part(krome_idx_C2H2) = 2.19d-5 + abundance_part(krome_idx_HCN) = 2.045d-5 + abundance_part(krome_idx_N2) = 2d-5 + abundance_part(krome_idx_SiC2) = 9.35d-6 + abundance_part(krome_idx_CS) = 5.3d-6 + abundance_part(krome_idx_SiS) = 2.99d-6 + abundance_part(krome_idx_SiO) = 2.51d-6 + abundance_part(krome_idx_CH4) = 1.75d-6 + abundance_part(krome_idx_H2O) = 1.275d-6 + abundance_part(krome_idx_HCl) = 1.625d-7 + abundance_part(krome_idx_C2H4) = 3.425d-8 + abundance_part(krome_idx_NH3) = 3d-8 + abundance_part(krome_idx_HCP) = 1.25d-8 + abundance_part(krome_idx_HF) = 8.5d-9 + abundance_part(krome_idx_H2S) = 2d-9 + abundance_part(krome_idx_e) = krome_get_electrons(abundance_part(:)) + +end subroutine chem_init + +end module analysis diff --git a/src/utils/analysis_krome.f90 b/src/utils/analysis_krome.f90 deleted file mode 100644 index 4d69e7f0f..000000000 --- a/src/utils/analysis_krome.f90 +++ /dev/null @@ -1,236 +0,0 @@ -!--------------------------------------------------------------------------! -! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! -! Copyright (c) 2007-2026 The Authors (see AUTHORS) ! -! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.github.io/ ! -!--------------------------------------------------------------------------! -module analysis -! -! Analysis routine comparing time between dumps -! -! :References: None -! -! :Owner: Mats Esseldeurs -! -! :Runtime parameters: None -! -! :Dependencies: eos, io, krome_main, krome_user, neighkdtree, part, -! physcon, raytracer, units -! - use krome_user, only:krome_nmols - use part, only: maxp - use raytracer, only: get_all_tau - implicit none - character(len=20), parameter, public :: analysistype = 'krome' - public :: do_analysis - logical, allocatable :: mask(:) - - real, allocatable :: abundance(:,:), abundance_prev(:,:), one(:) - character(len=16) :: abundance_label(krome_nmols) - integer(8), allocatable :: iorig_old(:) - integer, allocatable :: iprev(:) - logical :: done_init = .false. - real :: AuvAv = 4.65, albedo = 0.5 - -contains - -subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) - use part, only: isdead_or_accreted, iorig, rhoh, nptmass, xyzmh_ptmass, iReff - use neighkdtree, only:build_tree - use units, only: utime,unit_density - use eos, only: get_temperature, ieos, gamma,gmw, init_eos - use io, only: fatal - use krome_main, only: krome_init, krome - use krome_user, only: krome_get_names,krome_set_user_Auv,krome_set_user_xi,& - krome_set_user_alb,krome_set_user_AuvAv - character(len=*), intent(in) :: dumpfile - integer, intent(in) :: num,npart,iunit - real, intent(in) :: xyzh(:,:),vxyzu(:,:) - real, intent(in) :: particlemass,time - real, save :: tprev = 0. - integer, save :: nprev = 0 - real :: dt_cgs, rho_cgs, numberdensity, T_gas, gammai, mui, AUV, xi - real :: abundance_part(krome_nmols), Y(krome_nmols), column_density(npart), xyzh_copy(4,npart) - integer :: i, j, ierr, completed_iterations, npart_copy = 0 - - if (.not.done_init) then - done_init = .true. - print*, "initialising KROME" - call krome_init() - print*, "Initialised KROME" - abundance_label(:) = krome_get_names() - maxp = maxp /2 - allocate(abundance(krome_nmols,maxp)) - abundance = 0. - allocate(abundance_prev(krome_nmols,maxp)) - abundance_prev = 0. - allocate(one(maxp)) - one = 1. - allocate(iorig_old(maxp)) - iorig_old = 0 - allocate(iprev(maxp)) - iprev = 0 - print*, "setting abundances" - !$omp parallel do default(none) & - !$omp shared(npart,xyzh,vxyzu,dt_cgs,nprev,iorig,iorig_old,iprev) & - !$omp shared(abundance,abundance_prev,particlemass,unit_density) & - !$omp shared(ieos,rho_cgs,T_gas,j) & - !$omp private(i,abundance_part) - do i=1, npart - if (.not.isdead_or_accreted(xyzh(4,i))) then - call chem_init(abundance_part) - abundance(:,i) = abundance_part - endif - enddo - call init_eos(ieos, ierr) - if (ierr /= 0) call fatal(analysistype, "Failed to initialise EOS") - else - dt_cgs = (time - tprev)*utime - completed_iterations = 0 - print*, "not first step data, timestep = ",dt_cgs, "npart = ",npart, "nprev = ",nprev - xyzmh_ptmass(iReff,1) = 2. - npart_copy = npart - xyzh_copy = xyzh(:,:npart) - call build_tree(npart_copy,npart_copy,xyzh_copy,vxyzu) - call get_all_tau(npart, nptmass, xyzmh_ptmass, xyzh, one, 5, .false., column_density) - !$omp parallel do default(none) & - !$omp shared(npart,xyzh,vxyzu,dt_cgs,nprev,iorig,iorig_old,iprev) & - !$omp shared(abundance,abundance_prev,particlemass,unit_density) & - !$omp shared(ieos,gamma,gmw,time,completed_iterations,column_density,AuvAv,albedo) & - !$omp private(i,j,abundance_part,Y,rho_cgs,numberdensity,T_gas,gammai,mui,AUV,xi) - outer: do i=1,npart - if (.not.isdead_or_accreted(xyzh(4,i))) then - inner: do j=1,nprev - if (iorig(i) == iorig_old(j)) then - iprev(i) = j - exit inner - endif - enddo inner - if (j == iprev(i)) then - abundance_part(:) = abundance_prev(:,iprev(i)) - else - call chem_init(abundance_part) - endif - !Thermodynamic quantities - rho_cgs = rhoh(xyzh(4,i),particlemass)*unit_density - gammai = gamma - mui = gmw - numberdensity = rho_cgs / (mui * 1.6605E-24) - T_gas = get_temperature(ieos,xyzh(1:3, i),rhoh(xyzh(4,i),particlemass),vxyzu(:,i),gammai,mui) - T_gas = max(T_gas,20.0d0) - - !Radiation quantities - AUV = 4.65 * column_density(i) / 1.87e21 - xi = get_xi(AUV) - - call krome_set_user_Auv(AUV) - call krome_set_user_xi(xi) - call krome_set_user_alb(ALBEDO) - call krome_set_user_AuvAv(AuvAv) - - Y = abundance_part*numberdensity - call krome(Y,T_gas,dt_cgs) - abundance_part = Y/numberdensity - abundance(:,i) = abundance_part - endif - !$omp atomic - completed_iterations = completed_iterations + 1 - print*, 'Completed ', completed_iterations, ' of ', npart - enddo outer - endif - - call write_chem(npart, dumpfile) - nprev = npart - tprev = time - iorig_old(1:npart) = iorig(1:npart) - abundance_prev(:,1:npart) = abundance(:,1:npart) -end subroutine do_analysis - -real function get_xi(AUV) - use physcon, only:pi - real, intent(in) :: AUV - real :: xi - real :: W(6), GA(6), ceta - integer :: i - - W(1) = 0.17132449 - W(2) = 0.36076157 - W(3) = 0.46791393 - W(4) = W(1) - W(5) = W(2) - W(6) = W(3) - GA(1) = 0.93246951 - GA(2) = 0.66120939 - GA(3) = 0.23861919 - GA(4) = -GA(1) - GA(5) = -GA(2) - GA(6) = -GA(3) - - xi = 0.0 - do i=1,6 - ceta = (pi*GA(i)+pi)/2.0 - xi=xi+(W(i)*(sin(ceta)*exp((-AUV*ceta)/sin(ceta)))) - enddo - xi = (pi/4.0)*xi - - get_xi = xi - -end function get_xi - -subroutine write_chem(npart, dumpfile) - use krome_user, only:krome_idx_He,krome_idx_C,krome_idx_N,krome_idx_O,& - krome_idx_H,krome_idx_S,krome_idx_Fe,krome_idx_Si,krome_idx_Mg,& - krome_idx_Na,krome_idx_P,krome_idx_F,krome_idx_CO,krome_idx_C2H2,& - krome_idx_C2H,krome_idx_H2,krome_idx_SiNC,krome_idx_e - integer, intent(in) :: npart - character(len=*), intent(in) :: dumpfile - integer :: i, iu - - open(newunit=iu,file=dumpfile//'.comp',status='replace',action='write') - write(iu, *) '# H, He, C, N, O, S, Fe, Si, Mg, Na, P, F, CO, C2H2, C2H, H2, SiNC, e-' - do i=1, npart - write(iu, *) abundance(krome_idx_H, i), abundance(krome_idx_He, i), abundance(krome_idx_C, i), & - abundance(krome_idx_N, i), abundance(krome_idx_O, i), abundance(krome_idx_S, i), & - abundance(krome_idx_Fe, i), abundance(krome_idx_Si, i), abundance(krome_idx_Mg, i), & - abundance(krome_idx_Na, i), abundance(krome_idx_P, i), abundance(krome_idx_F, i), & - abundance(krome_idx_CO, i), abundance(krome_idx_C2H2, i), abundance(krome_idx_C2H, i), & - abundance(krome_idx_H2, i), abundance(krome_idx_SiNC, i), abundance(krome_idx_e, i) - enddo - close(iu) - -end subroutine write_chem - -subroutine chem_init(abundance_part) - use krome_user, only:krome_idx_H2,krome_idx_He,krome_idx_CO,krome_idx_C2H2,& - krome_idx_HCN,krome_idx_N2,krome_idx_SiC2,krome_idx_CS,& - krome_idx_SiS,krome_idx_SiO,krome_idx_CH4,krome_idx_H2O,& - krome_idx_HCl,krome_idx_C2H4,krome_idx_NH3,krome_idx_HCP,& - krome_idx_HF,krome_idx_H2S,krome_idx_e,krome_get_electrons - real, intent(out) :: abundance_part(krome_nmols) - - ! Initial abundances for the krome model taken from Agúndez et al. (2020) - ! H2, He, CO, C2H2, HCN, N2, SiC2, CS, SiS, SiO, CH4, H2O, HCl, C2H4, NH3, HCP, HF, H2S - abundance_part(:) = 0. - abundance_part(krome_idx_H2) = 0.5d0 - abundance_part(krome_idx_He) = 8.5d-2 - abundance_part(krome_idx_CO) = 4d-4 - abundance_part(krome_idx_C2H2) = 2.19d-5 - abundance_part(krome_idx_HCN) = 2.045d-5 - abundance_part(krome_idx_N2) = 2d-5 - abundance_part(krome_idx_SiC2) = 9.35d-6 - abundance_part(krome_idx_CS) = 5.3d-6 - abundance_part(krome_idx_SiS) = 2.99d-6 - abundance_part(krome_idx_SiO) = 2.51d-6 - abundance_part(krome_idx_CH4) = 1.75d-6 - abundance_part(krome_idx_H2O) = 1.275d-6 - abundance_part(krome_idx_HCl) = 1.625d-7 - abundance_part(krome_idx_C2H4) = 3.425d-8 - abundance_part(krome_idx_NH3) = 3d-8 - abundance_part(krome_idx_HCP) = 1.25d-8 - abundance_part(krome_idx_HF) = 8.5d-9 - abundance_part(krome_idx_H2S) = 2d-9 - abundance_part(krome_idx_e) = krome_get_electrons(abundance_part(:)) - -end subroutine chem_init - -end module analysis diff --git a/src/utils/analysis_raytracer.f90 b/src/utils/analysis_raytracer.f90 index 8b022f824..73fc2e611 100644 --- a/src/utils/analysis_raytracer.f90 +++ b/src/utils/analysis_raytracer.f90 @@ -351,11 +351,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) print*, 'Start calculating optical depth outward: ', trim(jstring) if (primsec(1,2) == 0. .and. primsec(2,2) == 0. .and. primsec(3,2) == 0.) then call system_clock(start) - call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, j, tau) + call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, j, .true., tau) call system_clock(finish) else call system_clock(start) - call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, j, tau) + call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, j, .true., tau) call system_clock(finish) endif timeTau = (finish-start)/1000. @@ -501,11 +501,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) call build_tree(npart2,npart2,xyzh2,vxyzu) if (primsec(1,2) == 0. .and. primsec(2,2) == 0. .and. primsec(3,2) == 0.) then call system_clock(start) - call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, order, tau) + call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, order, .true., tau) call system_clock(finish) else call system_clock(start) - call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, order, tau) + call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, order, .true., tau) call system_clock(finish) endif timeTau = (finish-start)/1000. @@ -521,11 +521,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) print*,'Start doing scaling analysis with order =',order if (primsec(1,2) == 0. .and. primsec(2,2) == 0. .and. primsec(3,2) == 0.) then call system_clock(start) - call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, order, tau) + call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, order, .true., tau) call system_clock(finish) else call system_clock(start) - call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, order, tau) + call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, order, .true., tau) call system_clock(finish) endif timeTau = (finish-start)/1000. @@ -574,11 +574,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) print*, 'Start calculating optical depth outward: ', trim(jstring) if (primsec(1,2) == 0. .and. primsec(2,2) == 0. .and. primsec(3,2) == 0.) then call system_clock(start) - call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, j, tau) + call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, j, .true., tau) call system_clock(finish) else call system_clock(start) - call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, j, tau) + call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, j, .true., tau) call system_clock(finish) endif timeTau = (finish-start)/1000. @@ -636,11 +636,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) print*,'Start calculating optical depth' if (primsec(1,2) == 0. .and. primsec(2,2) == 0. .and. primsec(3,2) == 0.) then call system_clock(start) - call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, order, tau) + call get_all_tau(npart2, 1, xyzmh_ptmass, xyzh2, kappa, order, .true., tau) call system_clock(finish) else call system_clock(start) - call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, order, tau) + call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, order, .true., tau) call system_clock(finish) endif timeTau = (finish-start)/1000. @@ -668,7 +668,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) order = 7 print*, 'Start calculating optical depth outward, order=',order call system_clock(start) - call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, order, tau) + call get_all_tau(npart2, 2, xyzmh_ptmass, xyzh2, kappa, order, .true., tau) call system_clock(finish) timeTau = (finish-start)/1000. print*,'Time = ',timeTau,' seconds.' diff --git a/src/utils/utils_raytracer_all.f90 b/src/utils/utils_raytracer_all.f90 index 966b5f98d..4d3b8f6ef 100644 --- a/src/utils/utils_raytracer_all.f90 +++ b/src/utils/utils_raytracer_all.f90 @@ -1064,7 +1064,7 @@ subroutine get_tau_inwards(point, primary, xyzh, neighbors, kappa, Rstar, tau) integer :: i, next, previous, nneigh integer, parameter :: nmaxcache = 0 - real :: xyzcache(0,nmaxcache) + real :: xyzcache(3, max(1,nmaxcache)) real :: ray(3), nextDist, previousDist, maxDist, dtaudr, previousdtaudr, nextdtaudr ray = primary - xyzh(1:3,point)