This report documents the procedure used to incorporate the consistent boundary flux method (CBF) into the ConMan code (a finite element code solving thermal convection). CBF is an accurate method to compute dynamic topography and details are presented in Zhong et al. [1993]. In this report, we will show how to modify your version of the ConMan code in order to incorporate the CBF method. One test example, including input files, is provided. Necessary subroutines may be obtained from Shijie Zhong or Michael Gurnis (email address: zhong@gps.caltech.edu or gurnis@caltech.edu).

In order to use the CBF in your ConMan code, several minor modifications need to be made in the subroutines elminp, timdrv, f_vstf, f_vrres and eg2, and common.h of your original ConMan code. The actual calculation of the CBF takes place in a new subroutine, str_CBF, which itself consists of the following new subroutines: stresn, gstiff_b, and gbrhs .

In subroutine elminp, the following statements need to be added to define arrays which will be used in str_CBF.

mpnelb1 = mpoint ( ' nelb1  ', nelx , 0    ,   0   ,      1 )

mpnelb2 = mpoint ( ' nelb2  ', nelx , 0    ,   0   ,      1 )

mpbstft    = mpoint ( 'bstifft ', nelx , 36   ,   0   ,  mprec )

mpbsrhst = mpoint ( 'bsrhst  ', nelx , 8    ,   0   ,  mprec )

mpbstfb   = mpoint ( 'bstiffb ', nelx , 36   ,   0   ,  mprec )

mpbsrhsb= mpoint ( 'bsrhsb  ', nelx , 8    ,   0   ,  mprec )

In subroutine timdrv, before any subroutine is called, variables yb, yt, nelx1, and nelz1 are assigned as follows:

yb = x(2,    1)

yt = x(2,numnp)

nelx1 = nelx + 1

nelz1 = nelz + 1

and within time step loop of subroutine timdrv, a call to EGlib is needed.

if(mod(lstep, nsmprt).eq.0) then

do 6500 neg = 1, numeg

           call EGlib (negrou(neg), 'str_new ')

6500     continue

end if

This subroutine is called only for those time steps in which normal stresses are to be computed. In addition, variables yb, yt, nelx1, and nelz1, and those array names (i.e. mpnelb1, mpnelb2, mpbstft, mpbsrhst, mpbstfb, mpbsrhsb) must be added into common.h.

In subroutine f_vstf, before the loop over the element blocks (i.e. the loop for iblk = 1, nelblk), several variables and arrays must be initilized as follows:

indext = 0

indexb = 0

if(mod(lstep, nsmprt).eq.0) then

do 9 i = 1, nelx

       nelb1(i) = 0

       nelb2(i) = 0

do 9 j = 1, 36

       bstifft(i,j) = zero

       bstiffb(i,j) = zero

9            continue

end if

Also in subroutine f_vstf, a call to subroutine gstiff_b needs to be added right before assembling elemental stiffness matrices into global stiffness matrix (i.e. directly before the call to subroutine g2assm).

if(mod(lstep, nsmprt).eq.0) then

call gstiff_b(stiff, x, iel, ien, indext, indexb,

&              bstifft, nelb1, bstiffb, nelb2)

end if

Subroutine gstiff_b is called to obtain local stiffness matrices for boundary elements and is called only for those time steps in which normal stresses are to be computed.

In subroutine f_vrres, before the loop over the element blocks (i.e. the loop for iblk = 1, nelblk), two variables must be initilized as follows:

indext= 0

indexb= 0

Again in subroutine f_vrres, a call to subroutine gbrhs must be added within the loop over the element blocks (i.e. the loop for iblk = 1, nelblk) and after elemental force vectors (i.e. el_rhs) are formed.

if(mod(lstep, nsmprt).eq.0) then

call gbrhs(el_rhs,x,ien,indext,indexb,iel,bsrhst,bsrhsb)

end if

Subroutine gbrhs is to obtain elemental forces for boundary elements.

In subroutine eg2, a call to subroutine stresn must be added

if (task .eq. 'str_new ') then

call stresn(a(mpv   ), a(mpvl   ), a(mpx    ), a(mpien  ),

&            a(mpcblk), a(mpbstft), a(mpnelb1), a(mpbstfb),

&            a(mpnelb2),a(mpbsrhst),a(mpbsrhsb))

return

endif

In subroutine stresn, normal stresses on the bottom and top boundaries are computed and are output in the file associated with file unit 'imout'. The first, second, and final columns are for x coordinates, normal stresses on top and bottom boundaries, respectively.

In subroutine eg2, the statements to call the subroutines f_vstf and f_vrres must be modified in order to pass arrays which are necessary for computing normal stresses.

           call  f_vStf(a(mpvisc) , a(mpalam) , a(mpien ) , a(mpshl ) ,

&                              ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

&                              a(mpbstft) , a(mpnelb1), a(mpbstfb) , a(mpnelb2))


                            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

                   call  f_vRes(a(mpien ) , a(mpshl ) , a(mpdet ) , a(mptl  ) ,

&                              ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

&                                 a(mpalam) , a(mpx   ) , a(mpbsrhst), a(mpbsrhsb))

Accordingly, necessary modifications must be made in subroutines f_vstf and f_vrres in order for passing these arrays.

A example is provided to check whether the CBF method works properly in your ConMan code. The input files, runfile, in.bm32 and geom.bm32, and the file with the output normal stresses, mean.bm32 are in the appendix.

Reference

Zhong, S., M. Gurnis, and G. Hulbert, Accurate determination of surface normal stress in viscous flow from a consistent boundary flux method, Phys. Earth Planet. Int., 78, 1-8, 1993.