divvector – AGILE

public subroutine divvector(qvec, ixImin1, ixImin2, ixImin3, ixImax1, ixImax2, ixImax3, ixOmin1, ixOmin2, ixOmin3, ixOmax1, ixOmax2, ixOmax3, divq, fourthorder, sixthorder)

Uses

- mod_global_parameters
Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.

Calculate divergence of a vector qvec within ixL

Arguments

Type	Intent	Optional		Name
double precision,	intent(in)		::	qvec(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3,1:ndir)
integer,	intent(in)		::	ixImin1
integer,	intent(in)		::	ixImin2
integer,	intent(in)		::	ixImin3
integer,	intent(in)		::	ixImax1
integer,	intent(in)		::	ixImax2
integer,	intent(in)		::	ixImax3
integer,	intent(in)		::	ixOmin1
integer,	intent(in)		::	ixOmin2
integer,	intent(in)		::	ixOmin3
integer,	intent(in)		::	ixOmax1
integer,	intent(in)		::	ixOmax2
integer,	intent(in)		::	ixOmax3
double precision,	intent(inout)		::	divq(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3)
logical,	intent(in),	optional	::	fourthorder	Default: false
logical,	intent(in),	optional	::	sixthorder	Default: false

Calls

Help

Called by

Help

Source Code

  subroutine divvector(qvec,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,&
     ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,divq,fourthorder,&
     sixthorder)
    use mod_global_parameters
    integer, intent(in)             :: ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
       ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3
    double precision, intent(in)    :: qvec(ixImin1:ixImax1,ixImin2:ixImax2,&
       ixImin3:ixImax3,1:ndir)
    double precision, intent(inout) :: divq(ixImin1:ixImax1,ixImin2:ixImax2,&
       ixImin3:ixImax3)
    logical, intent(in), optional   :: fourthorder
    !< Default: false
    logical, intent(in), optional   :: sixthorder
    !< Default: false
    logical                         :: use_4th_order
    logical                         :: use_6th_order
    double precision                :: qC(ixImin1:ixImax1,ixImin2:ixImax2,&
       ixImin3:ixImax3), invdx(1:ndim)
    integer                         :: jxOmin1,jxOmin2,jxOmin3,jxOmax1,jxOmax2,&
       jxOmax3, hxOmin1,hxOmin2,hxOmin3,hxOmax1,hxOmax2,hxOmax3, ixCmin1,&
       ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3, jxCmin1,jxCmin2,jxCmin3,&
       jxCmax1,jxCmax2,jxCmax3
    integer                         :: idims, ixmin1,ixmin2,ixmin3,ixmax1,&
       ixmax2,ixmax3, gxOmin1,gxOmin2,gxOmin3,gxOmax1,gxOmax2,gxOmax3, kxOmin1,&
       kxOmin2,kxOmin3,kxOmax1,kxOmax2,kxOmax3
    integer                         :: lxOmin1,lxOmin2,lxOmin3,lxOmax1,lxOmax2,&
       lxOmax3, fxOmin1,fxOmin2,fxOmin3,fxOmax1,fxOmax2,fxOmax3

    use_4th_order = .false.
    use_6th_order = .false.
    if (present(fourthorder)) use_4th_order = fourthorder
    if (present(sixthorder))  use_6th_order = sixthorder
    if(use_4th_order .and. use_6th_order) call &
       mpistop("divvector: using 4th and 6th order at the same time")

    if(use_4th_order) then
      if (.not. slab_uniform) call mpistop(&
         "divvector: 4th order only supported for slab geometry")
      ! Fourth order, stencil width is two
      ixmin1=ixOmin1-2;ixmin2=ixOmin2-2;ixmin3=ixOmin3-2;ixmax1=ixOmax1+2
      ixmax2=ixOmax2+2;ixmax3=ixOmax3+2;
    else if(use_6th_order) then
      ! Sixth order, stencil width is three
      if (.not. slab_uniform) call mpistop(&
         "divvector: 6th order only supported for slab geometry")
      ixmin1=ixOmin1-3;ixmin2=ixOmin2-3;ixmin3=ixOmin3-3;ixmax1=ixOmax1+3
      ixmax2=ixOmax2+3;ixmax3=ixOmax3+3;
    else
      ! Second order, stencil width is one
      ixmin1=ixOmin1-1;ixmin2=ixOmin2-1;ixmin3=ixOmin3-1;ixmax1=ixOmax1+1
      ixmax2=ixOmax2+1;ixmax3=ixOmax3+1;
    end if

    if (ixImin1>ixmin1.or.ixImax1<ixmax1.or.ixImin2>ixmin2.or.ixImax2<ixmax2.or.&
       ixImin3>ixmin3.or.ixImax3<ixmax3) call &
       mpistop("Error in divvector: Non-conforming input limits")

    invdx=1.d0/dxlevel
    divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3)=0.0d0

    if (slab_uniform) then
      do idims=1,ndim
        if(use_6th_order) then
          lxOmin1=ixOmin1+3*kr(idims,1);lxOmin2=ixOmin2+3*kr(idims,2)
          lxOmin3=ixOmin3+3*kr(idims,3);lxOmax1=ixOmax1+3*kr(idims,1)
          lxOmax2=ixOmax2+3*kr(idims,2);lxOmax3=ixOmax3+3*kr(idims,3);
          kxOmin1=ixOmin1+2*kr(idims,1);kxOmin2=ixOmin2+2*kr(idims,2)
          kxOmin3=ixOmin3+2*kr(idims,3);kxOmax1=ixOmax1+2*kr(idims,1)
          kxOmax2=ixOmax2+2*kr(idims,2);kxOmax3=ixOmax3+2*kr(idims,3);
          jxOmin1=ixOmin1+kr(idims,1);jxOmin2=ixOmin2+kr(idims,2)
          jxOmin3=ixOmin3+kr(idims,3);jxOmax1=ixOmax1+kr(idims,1)
          jxOmax2=ixOmax2+kr(idims,2);jxOmax3=ixOmax3+kr(idims,3);
          hxOmin1=ixOmin1-kr(idims,1);hxOmin2=ixOmin2-kr(idims,2)
          hxOmin3=ixOmin3-kr(idims,3);hxOmax1=ixOmax1-kr(idims,1)
          hxOmax2=ixOmax2-kr(idims,2);hxOmax3=ixOmax3-kr(idims,3);
          gxOmin1=ixOmin1-2*kr(idims,1);gxOmin2=ixOmin2-2*kr(idims,2)
          gxOmin3=ixOmin3-2*kr(idims,3);gxOmax1=ixOmax1-2*kr(idims,1)
          gxOmax2=ixOmax2-2*kr(idims,2);gxOmax3=ixOmax3-2*kr(idims,3);
          fxOmin1=ixOmin1-3*kr(idims,1);fxOmin2=ixOmin2-3*kr(idims,2)
          fxOmin3=ixOmin3-3*kr(idims,3);fxOmax1=ixOmax1-3*kr(idims,1)
          fxOmax2=ixOmax2-3*kr(idims,2);fxOmax3=ixOmax3-3*kr(idims,3);
          divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
             ixOmin3:ixOmax3)=divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
             ixOmin3:ixOmax3)+(-qvec(fxOmin1:fxOmax1,fxOmin2:fxOmax2,&
             fxOmin3:fxOmax3,idims)+9.d0*qvec(gxOmin1:gxOmax1,gxOmin2:gxOmax2,&
             gxOmin3:gxOmax3,idims)-45.d0*qvec(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
             hxOmin3:hxOmax3,idims)+qvec(lxOmin1:lxOmax1,lxOmin2:lxOmax2,&
             lxOmin3:lxOmax3,idims)-9.d0*qvec(kxOmin1:kxOmax1,kxOmin2:kxOmax2,&
             kxOmin3:kxOmax3,idims)+45.d0*qvec(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
             jxOmin3:jxOmax3,idims))/(60.d0*dxlevel(idims))
        else if(use_4th_order) then
          ! Use fourth order scheme
          kxOmin1=ixOmin1+2*kr(idims,1);kxOmin2=ixOmin2+2*kr(idims,2)
          kxOmin3=ixOmin3+2*kr(idims,3);kxOmax1=ixOmax1+2*kr(idims,1)
          kxOmax2=ixOmax2+2*kr(idims,2);kxOmax3=ixOmax3+2*kr(idims,3);
          jxOmin1=ixOmin1+kr(idims,1);jxOmin2=ixOmin2+kr(idims,2)
          jxOmin3=ixOmin3+kr(idims,3);jxOmax1=ixOmax1+kr(idims,1)
          jxOmax2=ixOmax2+kr(idims,2);jxOmax3=ixOmax3+kr(idims,3);
          hxOmin1=ixOmin1-kr(idims,1);hxOmin2=ixOmin2-kr(idims,2)
          hxOmin3=ixOmin3-kr(idims,3);hxOmax1=ixOmax1-kr(idims,1)
          hxOmax2=ixOmax2-kr(idims,2);hxOmax3=ixOmax3-kr(idims,3);
          gxOmin1=ixOmin1-2*kr(idims,1);gxOmin2=ixOmin2-2*kr(idims,2)
          gxOmin3=ixOmin3-2*kr(idims,3);gxOmax1=ixOmax1-2*kr(idims,1)
          gxOmax2=ixOmax2-2*kr(idims,2);gxOmax3=ixOmax3-2*kr(idims,3);
          divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
             ixOmin3:ixOmax3)=divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
             ixOmin3:ixOmax3)+(-qvec(kxOmin1:kxOmax1,kxOmin2:kxOmax2,&
             kxOmin3:kxOmax3,idims)+8.d0*qvec(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
             jxOmin3:jxOmax3,idims)-8.d0*qvec(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
             hxOmin3:hxOmax3,idims)+qvec(gxOmin1:gxOmax1,gxOmin2:gxOmax2,&
             gxOmin3:gxOmax3,idims))/(12.d0*dxlevel(idims))
        else
          ! Use second order scheme
          jxOmin1=ixOmin1+kr(idims,1);jxOmin2=ixOmin2+kr(idims,2)
          jxOmin3=ixOmin3+kr(idims,3);jxOmax1=ixOmax1+kr(idims,1)
          jxOmax2=ixOmax2+kr(idims,2);jxOmax3=ixOmax3+kr(idims,3);
          hxOmin1=ixOmin1-kr(idims,1);hxOmin2=ixOmin2-kr(idims,2)
          hxOmin3=ixOmin3-kr(idims,3);hxOmax1=ixOmax1-kr(idims,1)
          hxOmax2=ixOmax2-kr(idims,2);hxOmax3=ixOmax3-kr(idims,3);
          divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
             ixOmin3:ixOmax3)=divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
             ixOmin3:ixOmax3)+half*(qvec(jxOmin1:jxOmax1,jxOmin2:jxOmax2,&
             jxOmin3:jxOmax3,idims) - qvec(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
             hxOmin3:hxOmax3,idims))*invdx(idims)
        end if
      end do
    else
      do idims=1,ndim
        hxOmin1=ixOmin1-kr(idims,1);hxOmin2=ixOmin2-kr(idims,2)
        hxOmin3=ixOmin3-kr(idims,3);hxOmax1=ixOmax1-kr(idims,1)
        hxOmax2=ixOmax2-kr(idims,2);hxOmax3=ixOmax3-kr(idims,3);
        ixCmin1=hxOmin1;ixCmin2=hxOmin2;ixCmin3=hxOmin3;ixCmax1=ixOmax1
        ixCmax2=ixOmax2;ixCmax3=ixOmax3;
        jxCmin1=ixCmin1+kr(idims,1);jxCmin2=ixCmin2+kr(idims,2)
        jxCmin3=ixCmin3+kr(idims,3);jxCmax1=ixCmax1+kr(idims,1)
        jxCmax2=ixCmax2+kr(idims,2);jxCmax3=ixCmax3+kr(idims,3);
        if(stretched_dim(idims) .and. stretch_uncentered) then
          ! linear interpolation at cell interface along stretched dimension
          qC(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
             ixCmin3:ixCmax3)=block%surfaceC(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
             ixCmin3:ixCmax3,idims)*(qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
             ixCmin3:ixCmax3,idims)+0.5d0*block%dx(ixCmin1:ixCmax1,&
             ixCmin2:ixCmax2,ixCmin3:ixCmax3,idims)*(qvec(jxCmin1:jxCmax1,&
             jxCmin2:jxCmax2,jxCmin3:jxCmax3,idims)-qvec(ixCmin1:ixCmax1,&
             ixCmin2:ixCmax2,ixCmin3:ixCmax3,idims))/(block%x(jxCmin1:jxCmax1,&
             jxCmin2:jxCmax2,jxCmin3:jxCmax3,idims)-block%x(ixCmin1:ixCmax1,&
             ixCmin2:ixCmax2,ixCmin3:ixCmax3,idims)))
        else
          qC(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
             ixCmin3:ixCmax3)=block%surfaceC(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
             ixCmin3:ixCmax3,idims)*half*(qvec(ixCmin1:ixCmax1,ixCmin2:ixCmax2,&
             ixCmin3:ixCmax3,idims)+qvec(jxCmin1:jxCmax1,jxCmin2:jxCmax2,&
             jxCmin3:jxCmax3,idims))
        end if
        divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
           ixOmin3:ixOmax3)=divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
           ixOmin3:ixOmax3)+qC(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
           ixOmin3:ixOmax3)-qC(hxOmin1:hxOmax1,hxOmin2:hxOmax2,&
           hxOmin3:hxOmax3)
      end do
      divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
         ixOmin3:ixOmax3)=divq(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
         ixOmin3:ixOmax3)/block%dvolume(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
         ixOmin3:ixOmax3)
    end if
  end subroutine divvector

divvector Subroutine

Contents

Source Code

public subroutine divvector(qvec, ixImin1, ixImin2, ixImin3, ixImax1, ixImax2, ixImax3, ixOmin1, ixOmin2, ixOmin3, ixOmax1, ixOmax2, ixOmax3, divq, fourthorder, sixthorder)

Uses

Arguments

Calls

Called by

Source Code