Characteristic Boundary Conditions

as implemented in bc_lower_magnetic

Introduction

The boundary conditions coded in bc_lower_magnetic are based on characteristics; i.e. that the (essentially) flux conservative equations have been re-written as a set of linked advection equations

$dU/dt+\lambda dU/dz = Q - L$

(where the sources Q and sinks L do not contain any derivatives of the main primitive variables (r, px, py, pz, e, Bx, By, Bz). Once the onerous task of converting the equations into this form is done, the setting of boundary conditions is very simple. In principle. This process can be summarized by the adage that outgoing characteristics (or waves if you wish) are computed using only internal information (grid points), while incoming characteritics are set according the wishes of the simulator. In our case, usually (or in fact always), no incoming waves. When set up in this fashion the equations look as follows:

Of course in Bifrost we work with the conservative variables, so we convert the du{xyz} to dp{xyz} and the equations actually implemented at the boundary are then

There are a number of d's that appear in these equations. They essentially represent the information that either is going out of the system or coming in in, at the various wave speeds that the system can contain. In the case of the MHD equations as implemented in Bifrost there are 7 different wave speeds/waves; the entropy wave (d2), the Alfvenic wave (d3/d4), the slow mode or acoustic wave (d5/d6) and the fast or magnetosonic waves (d7/d8). These d's can be derived as follows:

where V defines the vector of primitive varibles and S the matrix, the inverse of which converts the primitive variables into the characteristic variables that move through the system at the various charateristic wave velocities (\lambda + uz) and (\lambda - uz), such that

The right had side C's can be summarized as the terms that do not contain any 'z' derivatives (since we are working on a boundary that is implemented on upper 'z'-axis boundary).

Hence, the d's can be split into two sets, outgoing and incoming. The outgoing waves should only contain information that is interior to the computional box, while the incoming should be defined by the requirements of the simulation being run.

• the outgoing d's are given by

• the incoming d's are set, in the case of Bifrost by the requirement that each characteristic has unchanging amplitude, S{-1}dV/dt = 0 which gives

(NB the numbering of the characteristics in Bifrost is not the standard found in many papers. Usually the characteristics are numbered from the in terms of wave speeds from largest negative to largets positive. We did not know this at the time we worked these equations out and have not had the heart to change the numbering scheme.)

Implementation

As noted above, in principle this is very straightforward. However, in actual practice it has been quite difficult to code up these boundaries in a way that has led to the stable non-reflecting boundaries we desire.

Points to consider

• The d's and other terms in the equations should not be dependent on the ghost zones.
• This is very important and the reason that special stagger routines ending in _lbc have been written for derivaties and interpolations in the 'z'-axis direction.
• There remains the possibility that gridzones 2 and 3 which use some ghost zones values in interpolations and derivations may be coupling hurting this requirement. It is worth testing, but will require a rewrite of some of the stagger_mpi routines.
• The small effect of the ghostzones seems smallest when using bctypelower = "mcce".

• Furthermore, the Um(:,:,:,3) and diff_r routines are dependent on ghost zones. The impact of this should be checked.

• The most sensitive d's (by far) are d3/d4, the Alfvenic characteristic, where the (p+B{2}/2) term has given a lot of trouble.

• Setting this term to zero is in principle wrong, but seems to give more stable results, and with velocities and a general evolution that is not that different than when the term is included.
• The ordering of interpolations and derivations of the various terms in this expression actaully matters. Be careful when changing!

• The term also appears in d5/d6 and d7/d8, but does not seem to play as significant a role there.

• The diffusion terms coded in diff_r, diff_pxyz and diff_bxyz are a bit tricky to code, but are important and cannot be removed without introducing much noise on the boundary.

The current implentation seems to work very well, and much better than the older versions which used ghost zone values and required recurrent smoothing and div(B) cleaning. This seems no longer necessary and reflections have almost disappeared from the boudary. Note however that most tests have been done in 2d, there may be new problems that arise in 3d.

As coded now the bc_lower_magnetic has many fewer possibilies of being "controlled" by smoothing or similar (i.e. div(B) cleaning of ghost zones), this is good, but may also be bad... (re-)introduce factors to control d's behavior to give more control back to user???

An update

Some months later and several 3d models have been run for extended periods of time. The boundary conditions seem to hold up well when using the input parameters given in the next section. However, it should be noted that occasionally, as in certain magnetic configurations perhaps there are still some problems that have not been adaquately understood to isolate the concrete issue at fault. In particular it appears that in certain situations the d3/d4 outgoing characteristics are giving trouble (these are dependent on "cross" derivatives so are not very important in 2d runs perhaps). This can be mitigated by setting the factors tau_d3_bcl/tau_d4_bcl to some large factor, say 100.

do k=izb,izb
     do j=ysb,yeb
        do i=xsb,xeb
! Incoming d3 characteristic given by
! -sz Ry C(px) + sz Rx C(py) + Ry C(Bx) - Rx C(By)(i,j,k)
           scrb3(i,j,k) = -sz(i,j,k)*ry(i,j,k)*cpx(i,j,k) &
                          +sz(i,j,k)*rx(i,j,k)*cpy(i,j,k)
! Outgoing d3 wave when uz < -cz [sqrt(Bz^2/r)]
         if((zupuz(i,j,zs) <= -cz(i,j,zs)).or.noincoming) then
           d3_h(i,j,1) = (zupuz(i,j,k)+cz(i,j,k)) * &
                      (-scrp1(i,j,1) + scrp2(i,j,1))
         else
           d3_h(i,j,1) = scrb3(i,j,k)/rsm(i,j,k)/tau_d3_bcl &
                      + sqrt(fudge(i,j,k)/rsm(i,j,k)) &
                       *(ry(i,j,k)*cbx(i,j,k)-rx(i,j,k)*cby(i,j,k))
         endif
! Outgoing d4 wave when uz < cz [sqrt(Bz^2/r)]

In 3d simulations it also seems necessary to sometimes run smoothing at the upper boundary, say every 1000 or 2000 timesteps, but this has not been verified properly (yet).

One thing to watch for is a growth of the horizontal field Bx and/or By at the top boundary. This occurs (very) occasionally when tau_{3,4}d_bcl are set to their nominal value of 1.0, but does not happen when these parameters are set to a higher value. This may indicate that there is an issue with the centering of these characteristics and may be worth pursuing for someone with extra time on their hands.

Usage in practice

The below values work well in the models that have been tested.

; ************************* From    bc_lower_magnetic *************
; $Id: bc_lower_magnetic (a6dd9e0) Mon Nov 30 20:25:50 2020 by Mikolaj $  
    bctypelower = "mcce"
       bcldebug =    0
        tau_bcl =  1.00000000E+00
     tau_ee_bcl =  1.00000000E+02
     tau_d2_bcl =  1.00000000E+02
     tau_d3_bcl =  1.00000000E+00
     tau_d4_bcl =  1.00000000E+00
             s0 =  1.66687336E+01
             e0 =  4.83814836E-01
             r0 =  1.43913195E-01
            cs0 =  7.05358803E-01
             p0 =  4.28083465E-02
    nsmooth_bcl =      1000
      naver_bcl =        -1
     nclean_bcl =      1000
     nclean_lbl =        -4
     nclean_ubl =         5