MHD boundary conditions in Dedalus v2

[Bindesh Tripathi]

Hi all,

While solving the magnetic induction equation in a vector potential (A) formulation, (a) ∂_tA − η∇²A + ∇Ψ = u × B, together with the Coulomb gauge condition (b) ∇.A=0, the numerical output for time>0 does not respect the following boundary conditions (BCs) that I impose in a Dedalus script:

(c) A_x(z= ± L_z /2) = 0,
(d) A_y(z= ± L_z /2) = 0,
(e) Ψ  (z= ± L_z /2) = 0.

Could I ask if I am missing anything trivial here, e.g., writing a separate BC for nx=0 and ny=0 mode? (The z-axis uses Chebyshev polynomials.)

Thanks,

Bindesh

[Jeffrey S. Oishi]

Hi Bindesh,

Apologies if you already know this, but Dedalus uses the roots grid and thus does not have the boundary values on the grid. If you mean something other than Ax[0] != 0, for example, please elaborate on the problems you're seeing.

Jeff

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/acedfaad-381f-44be-9bb4-2272d0bbf569n%40googlegroups.com.

[Bindesh Tripathi]

Hi Jeff,

Thank you for your quick response and the suggestion.  I apologize for being too brief earlier.

Aware of Dedalus using the roots grid, I printed the values of Ax[‘g’][:, :, 0] and Ax[‘g’] [:, :, -1], and similarly for Ay.  Both are zeros, ~10**(-17), at t=0.  Then, as the simulation proceeds, they both quickly become appreciably larger than zero, ~10**(-3), causing the magnetic fields, curl of A, to be larger, and the current density, curl of curl of A, to become enormous that almost all the resistive dissipation in the system happens at the boundaries; there is no issue with the viscous dissipation, though, and the hydro part of the code has always behaved well.  Addition of the magnetic field is the issue here, and likely the boundary conditions are the culprits.  Hoping to get away with this issue, I have run the code with higher resolution---512^3---as well, but to no avail.

The concerning portion of the code is shown below, on the appended form(, which begins on line number 372 in the full script, https://drive.google.com/file/d/1v3iTIjZWvkYXZvj9oE2lEKy-L3XDASMP/view?usp=sharing).  My guess is that the boundary conditions on the scalar gauge function Psi are likely the source of the issue; the pressure gauge is okay.  Your notes on the MHD boundary conditions that I came to find was very helpful in building this script.  Thank you for that.

Bindesh

---

A few nomenclature below: Ax->Afx, bx->bfx,  jx->jfx

Main equations:

problem.substitutions['L_fourier(thing)'] = "d(thing, x=2) + d(thing, y=2)"

problem.substitutions['uz_z'] = "- dx(ux) - dy(uy)"  #Use divergenceless condition while using z-derv. of last component, dz(uz) = - dx(ux) - dy(uy)

problem.substitutions['Afz_z']= "- dx(Afx)- dy(Afy)"  #Use divergenceless condition while using z-derv. of last component, dz(Az) = - dx(Ax) - dy(Ay)

 

problem.add_equation("dt(ux) + dx(p) - 1/Re*L_fourier(ux) - 1/Re*dz(ux_z) + jfz*sinth/MA2              = -u_dotgrad(ux) + jf_cross_bf_x/MA2 + Fx")

problem.add_equation("dt(uy) + dy(p) - 1/Re*L_fourier(uy) - 1/Re*dz(uy_z) - jfz*costh/MA2              = -u_dotgrad(uy) + jf_cross_bf_y/MA2     ")

problem.add_equation("dt(uz) + dz(p) - 1/Re*L_fourier(uz) - 1/Re*dz(uz_z) - (jfx*sinth-jfy*costh)/MA2  = -u_dotgrad(uz) + jf_cross_bf_z/MA2     ")

problem.add_equation("dx(ux) + dy(uy) + dz(uz) = 0")

problem.add_equation("dz(ux) - ux_z = 0")

problem.add_equation("dz(uy) - uy_z = 0")

 

problem.add_equation("dt(Afx) + dx(Psi) - 1/Rm*L_fourier(Afx) - 1/Rm*dz(Afx_z) + uz*sinth            = u_cross_bf_x")

problem.add_equation("dt(Afy) + dy(Psi) - 1/Rm*L_fourier(Afy) - 1/Rm*dz(Afy_z) - uz*costh            = u_cross_bf_y")

problem.add_equation("dt(Afz) + dz(Psi) - 1/Rm*L_fourier(Afz) - 1/Rm*dz(Afz_z) - (ux*sinth-uy*costh) = u_cross_bf_z")

problem.add_equation("dx(Afx) + dy(Afy) + dz(Afz) = 0")

problem.add_equation("dz(Afx) - Afx_z = 0")

problem.add_equation("dz(Afy) - Afy_z = 0")

 

BCs:

problem.add_bc("left(ux)        = -1")

problem.add_bc("right(ux)       =  1")

problem.add_bc("left(uy)        =  0")

problem.add_bc("right(uy)       =  0")

problem.add_bc("left(uz)        =  0")

problem.add_bc("right(uz)       =  0", condition="(nx != 0) or  (ny != 0)")

problem.add_bc("right(p)        =  0", condition="(nx == 0) and (ny == 0)")

 

problem.add_bc("left(Afx)       =  0")

problem.add_bc("right(Afx)      =  0")

problem.add_bc("left(Afy)       =  0")

problem.add_bc("right(Afy)      =  0")

problem.add_bc("left(Psi)       =  0")

problem.add_bc("right(Psi)      =  0")

 

---

[Daniel Lecoanet]

Hi Bindesh,

I think Jeff’s point is that z[0] and z[-1] are not the end points of the simulation, so if you print Ax[‘g’][:,:,0] and Ax[‘g’][:,:,-1], that doesn’t give you the values at the boundaries. Instead it might be more illuminating to print interp(Ax, z=z_bot) and interp(Ax, z=z_top).

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/0178db04-fa33-4539-ac2f-e3caafc8a975n%40googlegroups.com.

[BINDESH TRIPATHI]

Hi Daniel,

 

Many thanks for offering this note and the idea.  I had missed Jeff’s point, so I appreciate your comment.

 

Interpolating, or precisely, extrapolating, Ax to find its value at z=z_bot gives a number not very different from Ax[‘g’][:,:,0].  For example, a typical value of the former is -0.0008253 and of the latter is -0.0008265.

 

Using Paraview for 3D visualization, where I use data interpolated on a uniform grid, this issue of large current density at the top and bottom of the simulation domain can be clearly observed, please see the attached image.  Left plot is at t=0, and right plot is at the next time step, t=1.  The behavior of the current density evolution at around z=0 is exactly as anticipated, though---shearing by the mean flow.  Except those tiny regions near the top and bottom, everything else is okay.  Could I ask if you may have any further suggestions to investigate or change the BCs?

 

Thank you for taking time to offer ideas.

 

Bindesh

--
You received this message because you are subscribed to a topic in the Google Groups "Dedalus Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dedalus-users/etMq1wqJjV4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/C61B0113-A8E2-4A6C-B6AF-2B43E79DC8F9%40northwestern.edu.

[Bindesh Tripathi]

Hi all,

I'm now almost sure that the boundary conditions (BCs) on the vector potential A are the trouble maker.  Computing ∇.A shows that, for time>0, ∇.A is non-zero and larger especially near the boundaries.  This suggests that the BCs on A are not correctly getting imposed.  I have tried implementing several combinations of the BCs for perfectly conducting boundaries, including the following one, but the issue persists. (Likely problematic equations are highlighted in red.)

problem.add_bc("left(Afx)        =  0")
problem.add_bc("right(Afx)      =  0")
problem.add_bc("left(Afy)         =  0")
problem.add_bc("right(Afy)       =  0")

problem.add_bc("left(Psi)          =  0")

problem.add_bc("right(Psi)        =  0", condition="(nx != 0) or  (ny != 0)")
problem.add_bc("right(Afz)        =  0", condition="(nx == 0) and (ny == 0)")

The equations for the magnetic variables, based on ∂_tA + ∇Ψ − 1/Rm* ∇²A = u × B and ∇.A=0, are:

problem.add_equation("dt(Afx) + dx(Psi) - 1/Rm*L_fourier(Afx) - 1/Rm*dz(Afx_z) + uz*sinth                   = u_cross_bf_x")
problem.add_equation("dt(Afy) + dy(Psi) - 1/Rm*L_fourier(Afy) - 1/Rm*dz(Afy_z)  - uz*costh                   = u_cross_bf_y")
problem.add_equation("dt(Afz) + dz(Psi) - 1/Rm*L_fourier(Afz) - 1/Rm*dz(Afz_z) - (ux*sinth-uy*costh) = u_cross_bf_z")
problem.add_equation("dx(Afx) + dy(Afy) + dz(Afz) = 0")
problem.add_equation("dz(Afx) - Afx_z = 0")
problem.add_equation("dz(Afy) - Afy_z = 0")

Note that Ax, Ay, and Az are written here as Afx, Afy, and Afz to denote the fluctuating components of A; the background components of A give rise to a uniform background magnetic field, B_0 = costheta \hat{x} + sintheta \hat{y} . Such a magnetic field then breeds, via u x B_0 in the equation for ∂_tA, the terms with "sinth" and "costh" above.

Since the issue of imposing the BCs on the vector potential A appears to be so common, I wonder if any one here is able to share their thoughts.  Thanks!

Bindesh

[Ben Brown]

Bindesh,

     I think you want:

problem.add_bc("right(Afz)        =  0", condition="(nx != 0) or  (ny != 0)")
problem.add_bc("right(Psi)        =  0", condition="(nx == 0) and (ny == 0)")

I think you’ve got your conditions and fields the wrong way around in your script; you only want to impose the gauge choice for Psi on a single mode, and you had it imposed on all but one.  It follows the same form as the p gauge choice. 

Separately, check your A bcs with some care; I often get them wrong myself, but when I do, I go back to Jeff Oishi’s note to relearn how to do them right (I don’t know that these are wrong, I just know it’s a spot where my intuition tends to go awry)

—Ben

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/0558a76c-9953-46aa-a371-3d2bbbe1e236n%40googlegroups.com.

[Bindesh Tripathi]

Thank you, Ben, very much for suggesting the changes in the boundary conditions.  I now see the analogy between the pressure gauge and the Psi gauge here.

Jeff Oishi's note, unless it has been updated, suggests that for perfectly electrically conducting boundaries, Bz = jx = jy = 0, or equivalently, Ax=Ay=Psi=constant(=0) at the boundaries.  I convinced myself of the same for the problem at hand, as well, and thus imposed 

problem.add_bc("left(Psi)          =  0")

problem.add_bc("right(Psi)        =  0", condition="(nx != 0) or  (ny != 0)")

problem.add_bc("right(Afz)        =  0", condition="(nx == 0) and (ny == 0)").

But your point, Ben, also makes sense.  So, I am confused what the right boundary conditions should be for perfectly conducting boundaries.  I have a background magnetic field that has no Bz component, but Bx and By exist and are uniform throughout the volume.

Could I ask if there are further thoughts on this?

Bindesh

[Bindesh Tripathi]

Dear all,

The previous suggestions were very helpful.  Dedalus now produces nice results.  There is, however, one point that is worrisome: in violation of the imposed boundary conditions, there is an enormous resistive dissipation in the simulation at the top and the bottom layers (exactly at z=z[0] and z=z[-1] of the "Dedalus domain").  When the data at these two grid points are removed by hand while post-processing, the maximum dissipation at a point in the resulting slightly-reduced simulation domain drops by six (6) orders of magnitude, and gives results that are completely anticipated; please see the attached image for a visualization (where the plotted-data have been interpolated in a uniform grid such that the first and the last z-axis points of the new uniform grid are the same as that of the non-uniform Chebyshev grid).  

In the Chebyshev grid as well, the data fall off from the top and the bottom by 4 orders of magnitude as I go to z=z[1] from z[0], or to z=z[-2] from z[-1].  Then, walking down a few additional Chebyshev grids (less than 10 in number), the data fall off by another 2 orders of magnitude; this is consistent with the finding of the 6-orders-of-magnitude fall-off in the uniform grid, as the Chebyshev grid points are clustered near the boundaries.  Is there any suggestion on how to interpret this finding in Dedalus?  I wish to get rid of the large values at z=z[0] and z=z[-1].

Thanks,

Bindesh