Cylindrical configuration

[brandortt]

Hi @mrkondratyev. Thank you for this amazing work, it's very simple and yet makes possible to understand the basics of MHD FV methods, I wasn't able to find anything like this for months, very useful and well written work.

I wanted to ask you, if you knew how this code could be adapted to simulate 2D (r,z) cylindrical axis-symmetric configurations (variations on azimuthal angle theta are neglected) . Is it enough to use some jacobian in the onestep scripts or some modifications are needed in the Riemann solver too?

Thank you again for sharing this material and for your attention.

[mrkondratyev]

Greetings, @brandortt!

I am really glad, that you have found this project useful! To handle axisymmetric problems you don't need any modifications in the Riemann solver. However, you have to use Jacobian transformation to the cylindrical geometry from the Cartesian one. In the FV framework, several additions should be done, namely (first three are necessary even for the first order approximation, while the fourth one is essential only for higher-order methods):

1) You need to calculate cell volumes and cell face surface areas for the cylindrical grid. For instance, the volume is V(cell) = 0.5*(r(i+1/2)^2 - r_(i−1/2)^2) * dz ∗ 2pi, and so on. I guess, that the simplest way to include these features into the code is to write a function similar to "uniCartGrid" in "grid_setup.py". Since the code is written with a usage of the areas and volumes, it is pretty straightforward to incorporate various coordinate systems.

2) Further, the source terms should be inserted into "fluxcalc" routines of "onestep.py" files, because of different definitions of divergence/gradient operators in the cylindrical case, and a presence of centrifugal force/hoop stress there. For instance, one should add a pressure term and a centrifugal force "(P + rho*v_theta^2)/r" into the RHS of radial momentum equation for hydrodynamics, and so one. They could be added in the end of the "flux_calc_hydro" function.

3) If you need to solve the problems with rotation or toroidal magnetic fields, the additional boundary condition should be prescribed. For hydrodynamics it is zero radial and azimuthal velocities at the axis. It is very simple to include new boundary conditions, you can add them into "boundaries.py" with some new values of boundary marker.

4) This item is optional and could be added only for high-order methods, the PCM and the second-order PLM minmod methods will work OK without further modifications as well. In principle, you still can use other provided high-order routines without any correction, but they may possibly produce some artifacts near the axis. The monotonocity and accuracy of the quasi-monotonic high-order methods could be violated near geometry singularities, such as axis, and a high-order piecewise polynomial reconstruction should be done in volumetric coordinates. One should also construct arrays for the volumetric coordinates in the grid routines. The extension of the FV methods (PLM, PPM and WENO) to curvilinear geometries can be found in the paper by [A. Mignone, J. Comp. Phys., 270, 1, 784 (2014)]. The expressions are quite cumbersome for a general case, while the exact formulas for reconstruction process are provided for uniform cylindrical grids there, which could be helpful.

I hope I have addressed your questions. Feel free to ask, if something is not clear.

With kind regards, Ilya Kondratyev

[brandortt]

Hi Ilya @mrkondratyev , thank you very much for your reply and suggestions. They were very useful to solve my issue. I followed the points 1 and 2 to make a basic version of your code in cylindrical coordinates. You were right in the fact that only grid_setup.py (I described the new surfaces and volumes) and one_step_mhd.py needed to be changed. I developed the Euler equations in cylindrical coordinates and found out that basically every variable had some source terms. I added them at the end of a new routine I called "flux_calc_cyl_mhd". I added them to the Residual calculations, so far everything seemed to work ok. But then I realized that something from your original structure was in conflict with the cylindrical configuration. In fact, referring to the 2nd dimension which I decided to be the radial one, in "flux_calc_mhd" when you need to assign the variable residuals you multiply each flux by the surface and then divide the difference by the cell volume (ResVar=ResVar + ( FVar[:,1:] Sx2[:,1:] - FVar[:,:-1] Sx2[:,:-1] ) / Vols ). (Actually, i honestly didn't get why you needed to multiply the fluxes by the surface and divide them by the volume instead of multiplying them only by the dx1/2, this would have saved me from finding the following problem:). Now, suppose you have a uniform field, then FVar[:,1:] = FVar[:,-1] and in a rectangular configuration Sx2[:,1:] = Sx2[:,:-1], as expected, this will make the contribution to ResVar 0. In a cylindrical configuration, given uniform fields, we would still expect no fluxes between the cells, but your line of code gives us finite fluxes since Sx2[:,1:] != Sx2[:,:-1] in a cylindrical geometry. This led to some numerical errors which would make the simulation very unstable especially near the axis. I solved this by rewriting the Residuals lines in the following way in the new subroutine "flux_calc_cyl_mhd": ResVar=ResVar + ( FVar[:,1:] - FVar[:,:-1] ) ( Sx2[:,1:] + Sx2[:,:-1] ) / 2 / Vols : instead of premultiplying the fluxes by the surface I move them out and make an average of them. I don't know if this is completely correct, but it works, solving numerical instability and errors also at r=0. In fact, if the field is uniform the the fluxes are equal and cancel out no matter the size of the surfaces they come from, giving 0 contribution to the residuals. I also thought that, since the grid (x,r) is still rectangular, it would be also correct to write the above equation as follows: ResVar=ResVar + ( FVar[:,1:] - FVar[:,:-1] ) dx2. I didn't try it, but it should work. Furthermore, both this approaches will work in an original rectangular configurations (ofc, the surfaces are all the same in that configuration). What do you think about this problem? I don't know if I have been clear in explaining my doubts about it and about the solutions I found. Thank you again for your interest and attention. Kind regards, Brando Rettino

[mrkondratyev]

Hello Brando @brandortt !

I use the volumes and surfaces, rather than dx1 and dx2, in the residual calculations, because an introduction of different grid geometries is more simple and natural in this case, and I plan to add later the curvilinear coordinates to the code, as well as test problems, connected to 2D astrophysical fluid models. If we simply multiply or divide the fluxes on dx1 and dx2, the conservation laws will be violated on the cylindrical grid, while the form "ResVar=ResVar + ( FVar[:,1:] Sx2[:,1:] - FVar[:,:-1] Sx2[:,:-1] ) / Vols" conserves mass and total energy up to machine accuracy. More specifically, the flux divergence operator should be discretized as follows according to a geometry-independent DIV definition:

DIV(F)_cell = 1.0 / Vol_cell SUM_OVER_ALL_FACE_NEIBS_TO_THE_CELL(n_face S_face * F_face),

where n_face is a unit normal vector to the face surface, S_face is a face surface area, F_face is a Riemann flux, while V_cell is a cell volume, and the summation is taken over the faces, which are adjacent to the given cell. Of course you can write it in the form with dx1 and dx2, but in this case you have to weight the fluxes properly with radius and so on due to the different definition of DIV operator in cylindrical geometry. I guess, that a simple averaging of the surfaces may also break down the conservation property or even approximation of the numerical scheme.

If I undestand your issue right, here is my possible solution for it. First of all we need to calculate the surface areas at the faces and the volumes. The cell volume is

cVol{ij} = (R{i+1/2}2 - R_{i-1/2}2))pidZ_{j},

where R_{i-1/2} -- is a radius of the face on the inner side of the cell. The surface area in radial direction is

SR{i-1/2,j} =R{i-1/2}dZ_{j}2*pi,

while the one in Z-direction is

SZ{i,j-1/2} = (R{i+1/2}2 - R_{i-1/2}2))*pi.

The surfaces "live" of the faces of the cell in R or Z directions with half-integer indexes. I am not sure, that these definitions will be in conflict with other routines, because I describe all the surfaces and volumes as arrays rather than constant numbers even for the Cartesian case. However, I have introduced them for the whole grid including ghost cells, which could make the flux_calc function a bit more messy in indexing. By my experience, some misprints in indexing could also be introduced in curvilinear geometries at the stage of grid volumes/areas declaration, since they have different values at different locations, as you have mentioned.

After prescribing the volumes and faces, the following FV approximation for the conservation laws, governed by the Euler equations, can be written (for instance, I introduce density and R-momentum equations, and hereafter no index corresponds to cell (i,j), while (i+-) corresponds to (i+-1/2, j) and (j+-) corresponds to (i, j+-1/2) ):

Denst + ( Fdens{i+}SR{i+} - Fdens{i-}SR{i-} )/cVol + ( Fdens{j+}SZ{j+} - Fdens{j-}SZ_{j-} )/cVol = 0,

MomRt + ( Fmomr{i+}Sr{i+} - Fmomr{i-}Sr{i-} )/cVol + ( Fmomr{j+}SZ{j+} - Fmomr{j-}SZ{j-} )/cVol = (Pres + dens*velocity\theta**2 )/ cR,

where cR an average radius -- cR{i} = (R{i+1/2} + R{i-1/2})/2 = cVol/(SR{i+} - SR_{i-}).

If you consider a simple uniform flow without motion, Dens_t will be exactly zero, because Fdens = DensVelocity, as expected, while the MomR_t will also be zero, because the source term will balance the fluxes in R-direction. Since the flux for the stationary fluid is Fmomr = DensVelR**2 + Pres = Pres = const for all faces\cells, we will get:

MomRt = -Fmomr * (SR{i+} - SR_{i-})/cVol + Pres/cR = 0 exactly (up to round-off presicion).

If you utilize such form of the discrete equations, you will only have non-zero source terms for radial and azimuthal momentums, while the proper DIV discretization of mass, mom_z and energy fluxes will be taken into account by the definiotns.

With kind regards, Ilya Kondratyev