search

PrincetonUniversity / athena-public-version

(MOVED) Athena++ GRMHD code and adaptive mesh refinement (AMR) framework. Active repository --->
https://github.com/PrincetonUniversity/athena
BSD 3-Clause "New" or "Revised" License
160 stars 118 forks source link

Pressure defaulting to floor at lower polar boundary in 3D #74

Closed paytonrodman closed 3 years ago

paytonrodman commented 3 years ago

Hi all,

I'm a PhD student trying to construct a 3D MHD accretion disk simulation. I've just "extended" my simulation to cover [0,pi] in the azimuthal angle using the "polar" boundary condition, but I'm having problems on the lower boundary (theta/x2=0): the pressure in the first cell initially defaults to the pressure floor (1e-8) for some reason.

The pressure (primitive): pressure

The total energy (conserved): energy

Input file (note I've tried without mesh refinement and the issue is still there): athinput.disk_base_pol.txt

Configured with: python configure.py --prob=$prob --coord=spherical_polar --cxx=g++ --eos=adiabatic -hdf5 -b -mpi --hdf5_path=/usr/local/software/spack/spack-0.11.2/opt/spack/linux-rhel7-x86_64/gcc-5.4.0/hdf5-1.10.1-57yz6quppilpxttsbitop5nk6b3rtagg

I don't have any instances of w(IPR,k,j,i) or prim(IPR,k,j,i) in my problem generator, nor is there any discontinuity in any other variables (dens/rho, Bcc, etc), and no errors are given.

Any help or advice you can give would be much appreciated!

Cheers, Payton

felker commented 3 years ago

What is $prob defined as in your environment? Did you modify the problem generator?

paytonrodman commented 3 years ago

$prob is disk_base_pol, problem generator is based on the keplerian disk but with added magnetic fields and a few extra tidbits (different gravitational potential in source term)

The full pgen is here: disk_base_pol.cpp.txt

c-white commented 3 years ago

I see that you are checking for nan values in A3(), which is good, because at the pole x2 is 0 and you end up taking the log of 0. However, it looks like when you calculate B there are zeros that are not caught, leading to nan values. In particular, Face2Area() and Edge3Length() will be 0 exactly on the pole. Thus B near the pole becomes nan, corrupting the energy. This is masked in the output, because the code does a variable inversion between initialization and output, defaulting to floors in the bad cells.

I think everything will work if you just check the B initialization at the poles to make sure it doesn't become indeterminate.

paytonrodman commented 3 years ago

Oh, I never even suspected Face2Area() and Edge3Length(), thank you so much!