danieljprice / phantom

Phantom Smoothed Particle Hydrodynamics and Magnetohydrodynamics code
https://phantomsph.github.io
Other
103 stars 223 forks source link

Error appeared when evoluting star.in using star or grstar setup(GR=yes) #402

Closed Shigaharuki3012 closed 1 year ago

Shigaharuki3012 commented 1 year ago

Steps followed:

/phantom/scripts/writemake.sh star > Makefile
make SETUP=star GR=yes
make setup
./phantomsetup star.setup

I choose Polytrope , 524000 particles and automatically relaxing the star. After normally relaxing the star: ./phantom star.in Then errors:

<<< finished reading (hydro) file 

 ERROR: units are incorrect for GR, need G = c = 1
  ...but we have G =    1.00000000000000       and c =    686.518075865455     

 Centre of mass is at (x,y,z) = ( 4.657E-18, -4.573E-18,  4.490E-18)

 FATAL ERROR! initial: errors in particle data from file : errors = 1

Same error appeared when using grstar setup.

danieljprice commented 1 year ago

If you override compile time options on the command line you have to use GR=yes for both phantom and phantomsetup. Or better, just use grstar as the setup in the original writemake command…

Shigaharuki3012 commented 1 year ago

In fact if I use grstar as the setup in the original writemake command like :

/phantom/scripts/writemake.sh grstar > Makefile
make 
make setup
./phantomsetup star

After last step of the above:

--> ALLOCATING PART ARRAYS

 FATAL ERROR! memory: metrics allocation error

will appear.

Shigaharuki3012 commented 1 year ago

To add,same error above would not appear when using gfortran compiler.But when relaxing star:

rstar  =    1.0000000085730454       mstar =   0.99935806369265401       tdyn =    1.1110774265470595     
 WARNING: Accrete particles: but Metric = Minkowski
 Centre of mass is at (x,y,z) = (-2.254E-18,  2.859E-20, -1.885E-20)
 Particle setup OK

 RELAX-A-STAR-O-MATIC: Etherm:  0.500     Epot: -0.856     R*:   1.00    
       WILL stop WHEN: dens error <   1.00% AND Ekin/Epot <   1.000E-07 OR Iter=1000

-------->   TIME =    0.000    : full dump written to file relax_00000   <--------

 ERROR! get_u0 in rho2dens: 1/sqrt(-v_mu v^mu) ---> non-negative: v_mu v^mu

 ERROR! get_u0 in prim2cons: 1/sqrt(-v_mu v^mu) ---> non-negative: v_mu v^mu

 ERROR! get_u0 in rho2dens: 1/sqrt(-v_mu v^mu) ---> non-negative: v_mu v^mu

 ERROR! get_u0 in prim2cons: 1/sqrt(-v_mu v^mu) ---> non-negative: v_mu v^mu

 ERROR! get_u0 in rho2dens: 1/sqrt(-v_mu v^mu) ---> non-negative: v_mu v^mu

 ERROR! get_u0 in prim2cons: 1/sqrt(-v_mu v^mu) ---> non-negative: v_mu v^mu

 ERROR! get_u0 in rho2dens: 1/sqrt(-v_mu v^mu) ---> non-negative: v_mu v^mu

 ERROR! get_u0 in prim2cons: 1/sqrt(-v_mu v^mu) ---> non-negative: v_mu v^mu

FATAL ERROR! cons2prim: could not solve rootfinding on particle 6277

would appear.

This is a different error compared to #333. I don't know whether the coder in #333 can follow all the simulation steps normally after using gfortran compiler.

Shigaharuki3012 commented 1 year ago

Just to update,errors would not appear when relaxing star if using uniform density profile instead of polytrope for the star(For grstar case.That is,when GR=yes).

Update:for uniform distribution of density the error will be :(Using ifort compiler)

>>>>>>  s  t  r  e   t    c     h       m     a    p   p  i  n  g  <<<<<<
 stretching to match tabulated density profile in r direction
 density at r =   2.000000000000000E-004  is   0.238732414637067
 total mass      =   0.999999999988753
 >>>>>> done

 rstar  =    1.00000000000000       mstar =   0.999999999991998       tdyn =    1.11072073454764
  0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000  0.480869792729328      -0.280291103893986      -0.768026436902166       0.0000000000000.768026436902166       0.000000000000000E+000 -0.280291103893986       -5.06864549880007        2.37948028304181       -163138.76802643690228304181       0.5893.0846313844251      -0.768026436902166        2.37948028304181       0.582984190562906        4.77523376614390       0.00000E+    4.7752337                   N0000000000000E+000  -13.0846313844251        4.77523376614390       -26.9673372499710                          NaN
forrtl: severe (24): end-of-file during read, unit -4, file /proc/38447/fd/0
Image              PC                Routine            Line        Source
libifcoremt.so.5   00002B2107FEA97F  for__io_return        Unknown  Unknown
libifcoremt.so.5   00002B2108032C6D  for_read_seq_lis      Unknown  Unknown
phantomsetup       000000000052868A  Unknown               Unknown  Unknown
phantomsetup       0000000000524933  Unknown               Unknown  Unknown
phantomsetup       0000000000519717  Unknown               Unknown  Unknown
phantomsetup       0000000000520701  Unknown               Unknown  Unknown
phantomsetup       0000000000531C43  Unknown               Unknown  Unknown
phantomsetup       0000000000403A9E  Unknown               Unknown  Unknown
libc-2.17.so       00002B210A296555  __libc_start_main     Unknown  Unknown
phantomsetup       00000000004039A9  Unknown               Unknown  Unknown                                                                                     
~                                                                                

which seems like an error of reading arrays in some file.

danieljprice commented 1 year ago

I am closing this as from trying to reproduce this it seems like the issue is you are setting up a star with M=R=1 in relativity which implies a star with radius R = 1GM/c^2, which is smaller than it's Schwarzschild radius. Hence it is totally wrong to assume that the Newtonian self-gravity is a small perturbation to GR in this case

What is missing is that the code should at least give an error if you try to do this...