Multigrain 2-fluids simulation 'error reading dumpfile'

[StephaneMichoulier]

When trying to run a multigrain 2-fluids simulation with the dustydisc setup with more than 2 different grain sizes, an error appears. Apparently, phantom complains about third, fourth and so on dust types (type 9, 10, 11 ...).

To reproduce the error: Compile Phantom with the dustydisc setup. Use the 2-fluids method for dust. Add 3 or more different grain sizes. I tested the grain size distribution set with both 'log-space' and 'manually' option. Both produce the error with 3 or more grain sizes. Everything is fine with 1 or 2 grain sizes.

This is the error:

reading setup from file: multigraintest_00000.tmp on unit 9 FT:Phantom:2022.0.1: (hydro+dust): 21/02/2023 11:35:50.3 n(gas) = 100000 n(dust1) = 10000 n(dust2) = 10000 n(dust3) = 10000 npart = 130000 npart(total) = 130000 setting isothermal sound speed^2 (polyk) = 2.4999999999999996E-003 gamma = 1.0000000000000000
qfacdisc = 0.25000000000000000
time = 0.0000000000000000

ALLOCATING ALL ARRAYS Total memory allocated to arrays: 187.304 MB n = 130000 reading particles 1: 130000, from block 1 lims= 1- 130000 ... got 18 sink properties from 1 sink particles

ID| Mass | Racc | Macc | hsoft |

1| 1.0 | 1.0 | 0.0 | 1.0 |

particle i= 120001 is type 9 ... particle i= 130000 is type 9 ---npartoftypetot = 100000 0 0 0 0 0 10000 10000 10000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 npartoftypetotact = 100000 0 0 0 0 0 10000 10000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

ERROR! read_dump: particle type counts do not match header <<< finished reading (hydro) file

FATAL ERROR! initial: error reading dumpfile

Here, npartoftypetotact is not correctly initialize with respect to npartoftypetot.

[danieljprice]

I tracked this down to a few lines added in the count_particle_types routine:

subroutine count_particle_types(npartoftype)
 use part, only:iphase,iamtype,npart
 integer, intent(out) :: npartoftype(:)
 integer :: i, itype

 npartoftype(:) = 0
 do i = 1, npart
    itype = iamtype(iphase(i))
    if (itype  >  8) then
       print *, "particle i=", i,"is type", itype
       cycle
    endif
    npartoftype(itype) = npartoftype(itype) + 1
 enddo

end subroutine count_particle_types

a quick "git blame" shows these lines were added only recently, I think they are just some debugging lines left from a pull request to fix a few other things by @alisonkyoung1:

d08efab04 src/main/readwrite_dumps_fortran.F90 (Alison Young       2023-01-09 14:10:38 +0000 1941)     if (itype  >  8) then
24e8f1a6b src/main/readwrite_dumps_fortran.F90 (Alison Young       2022-06-27 10:29:37 +0100 1942)        print *, "particle i=", i,"is type", itype
24e8f1a6b src/main/readwrite_dumps_fortran.F90 (Alison Young       2022-06-27 10:29:37 +0100 1943)        cycle
d08efab04 src/main/readwrite_dumps_fortran.F90 (Alison Young       2023-01-09 14:10:38 +0000 1944)     endif

needless to say, we should delete these lines as they block simulations of multigrain dust...

[danieljprice]

fixed by #385