Issue with the production of vertical magnetic fields lines in spherical polar coordinates

[shilpasarkar30]

I am regenerating the Zhu and Stone 2018 paper where the vector potential A_phi for R>rmin is of the form aR^(m+1)+b/R (where R is in cylindrical coordinate system, r is in spherical coordinate system and a, b, m are constants). I am using spherical polar coordinates and disk.cpp (hydro) program present in the Athena++ directory. For div B = 0, I have to use this expression of A_phi. I have magnetic field components B_r and B_theta which are defined as geom_coeff_r*(aphi(k,j,i+1) – aphi(k,j,i))/pcoord->dx1f(i) and geom_coeff_theta*(aphi(k,j+1,i) – aphi(k,j,i))/pcoord->dx2f(j) (where x1 and x2 are r and theta respectively). I have few questions:

(1) B_r=1/(r sinθ) d/d θ [ sinθ A_phi], and B_ θ=-1/r d/dr (r A_phi) (curl of a vector quantity in spherical coordinate system)
When I transform them to B_R and B_z (cylindrical) analytically, we expect to get only B_z (vertical magnetic fields) component and the B_R should be 0 or negligible to computational precision. But while doing simulations I always get some residual B_R left. How can I do away with that? My magnetic fields look curved because of the presence of B_R. I have to give a cap of B_R<1e-5 as 0 to make the field lines look vertical.

The reason I found is the way we are finding the B’s. The lesser the dtheta or dr (or higher the resolution) the better we would match the analytical value (which comes from the basic definition of differentiation). But there is an upper limit to increasing resolution as well. I have attached here a .par file which I am using for generating the solution. The plots obtained with and without the cap are named as streamlines_BR<1e-7 and streamlines_all resp. The plots are generated using a program I wrote.

(2) From the program given in Athena++ directory (plot_spherical.py) the field lines look further curved. I am unable to identify the problem. The plot is attached and named as athena_program_plot.png

**I have always used visit for plotting and it always showed me straight vertical lines when I used only one processor and the output file was vtk. Now as I switched to parallel and used the hdf5 format the magnetic field lines were broken (attached visit_broken.jpg). This hdf5 file has a low resolution but is generated through the same .cpp as the other figs attached. I have loaded the .xml expressions when plotting these curves. Then I tried yt project but as far as I understood it does not support spherical polar coordinate system.

It would be very nice if you can kindly help me.

Regards,
Shilpa

fig_colormesh_index.pdf
streamlines_BR<1e-7.pdf

streamlines_all.pdf


disk_mhd_sph_mpi_hdf5_mesh.txt

[zhuzh1983]

Hi Shilpa, I don’t have the problem generator that is compatible with the latest Athena++. But I do have a copy of the old code back then. You could try my old code to see if it works. If it works, you could see how the fields were set up in the problem generator. Here I have attached the code and the plotting tools. You can download all the files from my website: http://www.physics.unlv.edu/~zhzhu/YuanFeng.tgz First, you untar the file Then, you go to the Yuanfeng directory and compile the code with: python configure.py --prob=disk_sph_CPD_newf_Photo_block --coord=spherical_polar -mpi -b -hdf5 If you use icc compiler you also need to add --cxx=icc Then, you can run the code after going to bin directory mpirun -np xx athena -i athinput.disk_sph You can change the resolution, grid structure as you want. Currently I have enabled mesh refinement and only study a 90 degree phi wedge. After you run the code, you can use my python script (athenappclass15newd.py and appnewdsimp.py) to plot the data. Since you may need to run the code for quite a while, I have also attached my file at 200 output. You can play with this output first. Then when your run reaches this point, you can see if the results are similar (Hopefully they will be similar). If you want to use my 200 data, you can try

ipython import athenappclass15newd athenappclass15newd.athdf('disk.out1',200,level=1) (this will generate python data files for plotting) import appnewdsimp appnewdsimp.twodrhovelB() (which will give you plots similar to Fig 22 in my paper) appnewdsimp.onedrfig3twoplot() appnewdsimp.onedradialfig2simple() appnewdsimp.twodstreamconfig9avgmultiwind()

These plotting tools will generate several ps or eps files. I have tried them. They work ok for me and those files are in the directory too for you to compare with. I think you only need to run one step to see the initial condition. I cannot guarantee that the python script works since it used python 2. But I think vtk and hdf5 can be read into Visit software. Thanks, Zhaohuan

…

----------------------------- Zhaohuan Zhu Associate Professor, Physics & Astronomy University of Nevada - Las Vegas Email: ***@***.*** Homepage: www.physics.unlv.edu/~zhzhu

On Jan 16, 2024, at 9:59 PM, shilpasarkar30 ***@***.***> wrote: I am regenerating the Zhu and Stone 2018 paper where the vector potential A_phi for R>rmin is of the form aR^(m+1)+b/R (where R is in cylindrical coordinate system, r is in spherical coordinate system and a, b, m are constants). I am using spherical polar coordinates and disk.cpp (hydro) program present in the Athena++ directory. For div B = 0, I have to use this expression of A_phi. I have magnetic field components B_r and B_theta which are defined as geom_coeff_r*(aphi(k,j,i+1) – aphi(k,j,i))/pcoord->dx1f(i) and geom_coeff_theta*(aphi(k,j+1,i) – aphi(k,j,i))/pcoord->dx2f(j) (where x1 and x2 are r and theta respectively). I have few questions: (1) B_r=1/(r sinθ) d/d θ [ sinθ A_phi], and B_ θ=-1/r d/dr (r A_phi) (curl of a vector quantity in spherical coordinate system) When I transform them to B_R and B_z (cylindrical) analytically, we expect to get only B_z (vertical magnetic fields) component and the B_R should be 0 or negligible to computational precision. But while doing simulations I always get some residual B_R left. How can I do away with that? My magnetic fields look curved because of the presence of B_R. I have to give a cap of B_R<1e-5 as 0 to make the field lines look vertical. The reason I found is the way we are finding the B’s. The lesser the dtheta or dr (or higher the resolution) the better we would match the analytical value (which comes from the basic definition of differentiation). But there is an upper limit to increasing resolution as well. I have attached here a .par file which I am using for generating the solution. The plots obtained with and without the cap are named as streamlines_BR<1e-7 and streamlines_all resp. The plots are generated using a program I wrote. (2) From the program given in Athena++ directory (plot_spherical.py) the field lines look further curved. I am unable to identify the problem. The plot is attached and named as athena_program_plot.png It would be very nice if you can kindly help me. fig_colormesh_index.pdf <https://github.com/PrincetonUniversity/athena/files/13959185/fig_colormesh_index.pdf> Regards, Shilpa streamlines_BR<1e-7.pdf <https://github.com/PrincetonUniversity/athena/files/13959189/streamlines_BR.1e-7.pdf> streamlines_BR.1e-7_zoom.png (view on web) <https://github.com/PrincetonUniversity/athena/assets/59567763/051b0f85-9fa7-4615-af6e-a24fcfda8f30> streamlines_all.pdf <https://github.com/PrincetonUniversity/athena/files/13959194/streamlines_all.pdf> streamlines_all_zoom.png (view on web) <https://github.com/PrincetonUniversity/athena/assets/59567763/31c1543a-8739-4ae3-8f7b-80545f1c9915> athena_program_plot.png (view on web) <https://github.com/PrincetonUniversity/athena/assets/59567763/9f6c61d5-750a-409f-a6ef-128315c6be0a> disk_mhd_sph_mpi_hdf5_mesh.txt <https://github.com/PrincetonUniversity/athena/files/13959195/disk_mhd_sph_mpi_hdf5_mesh.txt> — Reply to this email directly, view it on GitHub <#555>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACB3RPTXZIYMI3YFR6N4X7LYO5SDBAVCNFSM6AAAAABB56UNMSVHI2DSMVQWIX3LMV43ERDJONRXK43TNFXW4OZWGA4DIOBYHE>. You are receiving this because you are subscribed to this thread.

[shilpasarkar30]

Thanx a lot Prof.Zhu. I would run these files and let you know. Thank you again!