mod_calculate_xw.t

Go to the documentation of this file.

!> Handles computations for coordinates and variables in output
module mod_calculate_xw
  implicit none
 
contains
 
  !> Compute both corner as well as cell-centered values for output
  subroutine calc_grid(qunit,igrid,xC,xCC,xC_TMP,xCC_TMP,wC_TMP,wCC_TMP,normconv,&
                       ixC^L,ixCC^L,first)
    ! this subroutine computes both corner as well as cell-centered values
    ! it handles how we do the center to corner averaging, as well as
    ! whether we switch to cartesian or want primitive or conservative output,
    ! handling the addition of B0 in B0+B1 cases, ...
    !
    ! the normconv is passed on to usr_aux_output for extending with
    ! possible normalization values for the nw+1:nw+nwauxio entries
 
    use mod_usr_methods, only: usr_aux_output
    use mod_global_parameters
    use mod_limiter
    use mod_physics, only: phys_energy, physics_type, phys_to_primitive
 
    integer, intent(in) :: qunit, igrid
    double precision, intent(in), dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
    double precision, intent(in), dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC
    integer :: ixC^L,ixCC^L
    logical, intent(in) :: first
 
    double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
    double precision, dimension(ixMlo^D:ixMhi^D,ndim)   :: xCC_TMP
    double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio)   :: wC_TMP
    double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio)     :: wCC_TMP
    double precision,dimension(0:nw+nwauxio),intent(out)       :: normconv
 
    double precision :: ldw(ixG^T), dwC(ixG^T)
    double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio)   :: wC
    double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio)     :: wCC
    double precision, dimension(ixG^T,1:nw+nwauxio)   :: w
    integer :: nxCC^D,idims,jxC^L,iwe
    integer :: nx^D, nxC^D, ix^D, ix, iw, level, idir
    logical, save :: subfirst=.true.
 
    ! initialize w
    w=0.d0
 
    ixcmin^d=ixmlo^d-1; ixcmax^d=ixmhi^d; ! Corner indices
    ixccmin^d=ixmlo^d; ixccmax^d=ixmhi^d; ! Center indices
 
    nx^d=ixmhi^d-ixmlo^d+1;
    level=node(plevel_,igrid)
 
    normconv(0) = length_convert_factor
    normconv(1:nw) = w_convert_factor
    block=>ps(igrid)
    w(ixg^t,1:nw)=ps(igrid)%w(ixg^t,1:nw)
 
    if (nwextra>0) then
     ! here we actually fill the ghost layers for the nwextra variables using
     ! continuous extrapolation (as these values do not exist normally in ghost cells)
     do idims=1,ndim
      select case(idims)
       {case(^d)
         jxcmin^dd=ixghi^d+1-nghostcells^d%jxCmin^dd=ixglo^dd;
         jxcmax^dd=ixghi^dd;
         do ix^d=jxcmin^d,jxcmax^d
             w(ix^d^d%jxC^s,nw-nwextra+1:nw) = w(jxcmin^d-1^d%jxC^s,nw-nwextra+1:nw)
         end do
         jxcmin^dd=ixglo^dd;
         jxcmax^dd=ixglo^d-1+nghostcells^d%jxCmax^dd=ixghi^dd;
         do ix^d=jxcmin^d,jxcmax^d
             w(ix^d^d%jxC^s,nw-nwextra+1:nw) = w(jxcmax^d+1^d%jxC^s,nw-nwextra+1:nw)
         end do \}
      end select
     end do
    end if
 
    ! next lines needed when usr_aux_output uses gradients
    ! and later on when dwlimiter2 is used
    ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
    if(nwauxio>0)then
      ! auxiliary io variables can be computed and added by user
      ! next few lines ensure correct usage of routines like divvector etc
      ! default (no) normalization for auxiliary variables
      normconv(nw+1:nw+nwauxio)=one
 
      if (.not. associated(usr_aux_output)) then
         call mpistop("usr_aux_output not defined")
      else
         call usr_aux_output(ixg^ll,ixm^ll^ladd1,w,ps(igrid)%x,normconv)
      end if
    endif
 
    ! In case primitives to be saved: use primitive subroutine
    !  extra layer around mesh only needed when storing corner values and averaging
    if(saveprim.and.first) call phys_to_primitive(ixg^ll,ixm^ll^ladd1,w(ixg^t,1:nw),ps(igrid)%x)
 
    if(b0field) then
    ! B0+B1 split handled here
      if(.not.saveprim.and.phys_energy) then
        w(ixg^t,iw_e)=w(ixg^t,iw_e)+0.5d0*sum(ps(igrid)%B0(ixg^t,:,0)**2,dim=ndim+1) &
              + sum(w(ixg^t,iw_mag(:))*ps(igrid)%B0(ixg^t,:,0),dim=ndim+1)
      end if
      w(ixg^t,iw_mag(:))=w(ixg^t,iw_mag(:))+ps(igrid)%B0(ixg^t,:,0)
    end if
    ! compute the cell-center values for w first
    ! cell center values obtained from mere copy
    wcc(ixcc^s,:)=w(ixcc^s,:)
 
    ! compute the corner values for w now by averaging
 
    if(slab_uniform) then
       ! for slab symmetry: no geometrical info required
       do iw=1,nw+nwauxio
         {do ix^db=ixcmin^db,ixcmax^db\}
            wc(ix^d,iw)=sum(w(ix^d:ix^d+1,iw))/dble(2**ndim)
         {end do\}
       end do
    else
       do iw=1,nw+nwauxio
         {do ix^db=ixcmin^db,ixcmax^db\}
           wc(ix^d,iw)=sum(w(ix^d:ix^d+1,iw)*ps(igrid)%dvolume(ix^d:ix^d+1)) &
                    /sum(ps(igrid)%dvolume(ix^d:ix^d+1))
         {end do\}
       end do
    endif
 
    if(nocartesian) then
      ! keep the coordinate and vector components
      xc_tmp(ixc^s,1:ndim)          = xc(ixc^s,1:ndim)
      wc_tmp(ixc^s,1:nw+nwauxio)    = wc(ixc^s,1:nw+nwauxio)
      xcc_tmp(ixcc^s,1:ndim)        = xcc(ixcc^s,1:ndim)
      wcc_tmp(ixcc^s,1:nw+nwauxio)  = wcc(ixcc^s,1:nw+nwauxio)
    else
      ! do all conversions to cartesian coordinates and vector components
      ! start for the corner values
      call to_cartesian(xc_tmp,wc_tmp,ixc^l,xc,wc)
      ! then cell center values
      call to_cartesian(xcc_tmp,wcc_tmp,ixcc^l,xcc,wcc)
    endif
 
    ! Warning: differentiate between idl/idlCC/tecplot/tecplotCC/vtu(B)/vtu(B)CC
    if(nwaux>0 .and. mype==0 .and. first.and.subfirst) then
      ! when corner values are computed and auxiliaries present: warn!
      if(convert_type=='idl'.or.convert_type=='tecplot' &
       .or.convert_type=='vtu'.or.convert_type=='vtuB') &
          write(*,*) 'Warning: also averaged auxiliaries within calc_grid'
      subfirst=.false.
    endif
 
  end subroutine calc_grid
 
  !> convert to cartesian coordinates and vector components
  subroutine to_cartesian(x_TMP,w_TMP,ix^L,xC,wC)
    ! conversion of coordinate and vector components from cylindrical/spherical
    ! to cartesian coordinates and components done here
    ! Also: nullifying values lower than smalldouble
    use mod_global_parameters
    use mod_geometry
 
    integer :: ix^L, ix^D, idim, iw, ivector, iw0
    integer, dimension(nw) :: vectoriw
    double precision :: x_TEC(ndim), w_TEC(nw+nwauxio)
    double precision, dimension(ndim,ndim) :: normal
 
    double precision, dimension(ix^S,ndim) :: xC
    double precision, dimension(ix^S,nw+nwauxio)   :: wC
 
    double precision, dimension(ix^S,ndim) :: x_TMP
    double precision, dimension(ix^S,nw+nwauxio)   :: w_TMP
 
    iw0=0
    vectoriw=-1
    if(nvector>0) then
      do ivector=1,nvector
         do idim=1,ndim
            vectoriw(iw_vector(ivector)+idim)=iw_vector(ivector)
         end do
      end do
    endif
    x_tec=0.d0
    {do ix^db=ixmin^db,ixmax^db\}
       select case (coordinate)
       case (cartesian,cartesian_stretched,cartesian_expansion)
          x_tec(1:ndim)=xc(ix^d,1:ndim)
          w_tec(1:nw+nwauxio)=wc(ix^d,1:nw+nwauxio)
       case (cylindrical)
          {^ifoned
          x_tec(1)=xc(ix^d,1)}
          {^iftwod
          select case (phi_)
          case (2)
             x_tec(1)=xc(ix^d,1)*cos(xc(ix^d,2))
             x_tec(2)=xc(ix^d,1)*sin(xc(ix^d,2))
          case default
             x_tec(1)=xc(ix^d,1)
             x_tec(2)=xc(ix^d,2)
          end select}
          {^ifthreed
          x_tec(1)=xc(ix^d,1)*cos(xc(ix^d,phi_))
          x_tec(2)=xc(ix^d,1)*sin(xc(ix^d,phi_))
          x_tec(3)=xc(ix^d,z_)}
 
          if (nvector>0) then
             {^ifoned normal(1,1)=one}
 
             {^iftwod
             select case (phi_)
             case (2)
                normal(1,1)=cos(xc(ix^d,2))
                normal(1,2)=-sin(xc(ix^d,2))
                normal(2,1)=sin(xc(ix^d,2))
                normal(2,2)=cos(xc(ix^d,2))
             case default
                normal(1,1)=one
                normal(1,2)=zero
                normal(2,1)=zero
                normal(2,2)=one
             end select}
 
             {^ifthreed
             normal(1,1)=cos(xc(ix^d,phi_))
             normal(1,phi_)=-sin(xc(ix^d,phi_))
             normal(1,z_)=zero
 
             normal(2,1)=sin(xc(ix^d,phi_))
             normal(2,phi_)=cos(xc(ix^d,phi_))
             normal(2,z_)=zero
 
             normal(3,1)=zero
             normal(3,phi_)=zero
             normal(3,z_)=one}
          end if
          do iw=1,nw+nwauxio
             if (iw<=nw) iw0=vectoriw(iw)
             if (iw0>=0.and.iw<=iw0+ndim.and.iw<=nw) then
                idim=iw-iw0
                w_tec(iw0+idim)={^d&wc(ix^dd,iw0+^d)*normal(idim,^d)+}
             else
                w_tec(iw)=wc(ix^d,iw)
             end if
          end do
       case (spherical)
          x_tec(1)=xc(ix^d,1){^nooned*sin(xc(ix^d,2))}{^ifthreed*cos(xc(ix^d,3))}
          {^iftwod
          x_tec(2)=xc(ix^d,1)*cos(xc(ix^d,2))}
          {^ifthreed
          x_tec(2)=xc(ix^d,1)*sin(xc(ix^d,2))*sin(xc(ix^d,3))
          x_tec(3)=xc(ix^d,1)*cos(xc(ix^d,2))}
 
          if (nvector>0) then
             {^ifoned normal(1,1)=one}
             {^nooned
             normal(1,1)=sin(xc(ix^d,2)){^ifthreed*cos(xc(ix^d,3))}
             normal(1,2)=cos(xc(ix^d,2)){^ifthreed*cos(xc(ix^d,3))
             normal(1,3)=-sin(xc(ix^d,3))}}
 
             {^iftwod
             normal(2,1)=cos(xc(ix^d,2))
             normal(2,2)=-sin(xc(ix^d,2))}
             {^ifthreed
             normal(2,1)=sin(xc(ix^d,2))*sin(xc(ix^d,3))
             normal(2,2)=cos(xc(ix^d,2))*sin(xc(ix^d,3))
             normal(2,3)=cos(xc(ix^d,3))
 
             normal(3,1)=cos(xc(ix^d,2))
             normal(3,2)=-sin(xc(ix^d,2))
             normal(3,3)=zero}
          end if
          do iw=1,nw+nwauxio
             if (iw<=nw) iw0=vectoriw(iw)
             if (iw0>=0.and.iw<=iw0+ndim.and.iw<=nw) then
                idim=iw-iw0
                w_tec(iw0+idim)={^d&wc(ix^dd,iw0+^d)*normal(idim,^d)+}
             else
                w_tec(iw)=wc(ix^d,iw)
             end if
          end do
       case default
          write(*,*) "No converter for coordinate=",coordinate
       end select
       x_tmp(ix^d,1:ndim)=x_tec(1:ndim)
       w_tmp(ix^d,1:nw+nwauxio)=w_tec(1:nw+nwauxio)
       ! Be aware that small values are nullified here!!!
       where(dabs(w_tmp(ix^d,1:nw+nwauxio))<smalldouble)
             w_tmp(ix^d,1:nw+nwauxio)=zero
       endwhere
    {end do\}
 
  end subroutine to_cartesian
 
  !> get all variables names
  subroutine getheadernames(wnamei,xandwnamei,outfilehead)
    ! this collects all variables names in the wnamei character array, getting the info from
    ! the prim_wnames/cons_wnames strings (depending on saveprim). It combines this info with names
    ! for the dimensional directions in the xandwnamei array. In the outfilehead, it collects
    ! the dimensional names, and only those names from the nw variables for output (through w_write)
    ! together with the added names for nwauxio variables
    use mod_usr_methods, only: usr_add_aux_names
    use mod_global_parameters
    use mod_geometry
 
    character(len=name_len)   :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
    character(len=1024) :: outfilehead
 
    integer::  space_position,iw,ind
    character(len=name_len)::  wname
    character(len=std_len):: aux_variable_names
    character(len=std_len)::  scanstring
 
    logical, save:: first=.true.
 
    ! in case additional variables are computed and stored for output
    if (nwauxio>0) then
       if (.not. associated(usr_add_aux_names)) then
          call mpistop("usr_add_aux_names not defined")
       else
          call usr_add_aux_names(aux_variable_names)
       end if
    end if
 
    ! --- part to provide variable names from prim_wnames/cons_wnames strings
    if(saveprim) then
       scanstring=trim(aux_variable_names)
       wnamei(1:nw)=prim_wnames(1:nw)
    else
       scanstring=trim(aux_variable_names)
       wnamei(1:nw)=cons_wnames(1:nw)
    endif
 
    space_position=index(scanstring,' ')
    do iw=nw+1,nw+nwauxio
       do while (space_position==1)
         scanstring=scanstring(2:)
         space_position=index(scanstring,' ')
       enddo
       wname=scanstring(:space_position-1)
       scanstring=scanstring(space_position+1:)
       space_position=index(scanstring,' ')
 
       ! fill all names, even those that we will not write away from the first nw
       wnamei(iw)=trim(wname)
    enddo
    ! --- end of part to provide variable names
 
    select case (coordinate)
       case( spherical )
          xandwnamei(1)="r";{^nooned xandwnamei(2)="Theta"};{^ifthreed xandwnamei(3)="Phi"}
       case( cylindrical )
          xandwnamei(1)="R";
          {^nooned
          if( phi_ == 2 )then
             xandwnamei(2)="Phi"
          else
             xandwnamei(2)="Z"
          endif}
          {^ifthreed
          if( phi_ == 2 )then
             xandwnamei(3)="Z"
          else
             xandwnamei(3)="Phi"
          endif}
       case default
          xandwnamei(1)="X";{^nooned xandwnamei(2)="Y"};{^ifthreed xandwnamei(3)="Z"}
    end select
 
    xandwnamei(ndim+1:ndim+nw+nwauxio)=wnamei(1:nw+nwauxio)
 
    ! in outfilehead, collect the dimensional names, and all output variable names
    ! first all dimensions
    write(outfilehead,'(a)') trim(xandwnamei(1))
    {^nooned
    do iw=2,ndim
       wname=xandwnamei(iw)
    write(outfilehead,'(a)')outfilehead(1:len_trim(outfilehead))//" "//trim(wname)
    enddo
    }
    ! then all nw variables, with w_write control for inclusion
    do iw=ndim+1,ndim+nw
       wname=xandwnamei(iw)
       if(w_write(iw-ndim)) then
    write(outfilehead,'(a)')outfilehead(1:len_trim(outfilehead))//" "//trim(wname)
       endif
    enddo
    ! then all nwauxio variables
    if(nwauxio>0) then
      do iw=ndim+nw+1,ndim+nw+nwauxio
         wname=xandwnamei(iw)
    write(outfilehead,'(a)')outfilehead(1:len_trim(outfilehead))//" "//trim(wname)
      enddo
    endif
 
    if(first.and.mype==0)then
      print*,'-------------------------------------------------------------------------------'
      write(unitterm,*)'Saving visual data. Coordinate directions and variable names are:'
      ind=0
      do iw=1,ndim
        ind=ind+1
        print *,ind,xandwnamei(iw)
      enddo
      do iw=1+ndim,nw+ndim
        if(w_write(iw-ndim)) then
          ind=ind+1
          write(*,*) ind,wnamei(iw-ndim)
        end if
      end do
      do iw=ndim+nw+1,ndim+nw+nwauxio
        ind=ind+1
        print *,ind,wnamei(iw-ndim)
      enddo
      write(unitterm,*)'time =', global_time
      print*,'-------------------------------------------------------------------------------'
      first=.false.
    endif
 
  end subroutine getheadernames
 
  !> computes cell corner (xC) and cell center (xCC) coordinates
  subroutine calc_x(igrid,xC,xCC)
    use mod_global_parameters
 
    integer, intent(in)               :: igrid
    double precision, intent(out)     :: xC(ixMlo^D-1:ixMhi^D,ndim)
    double precision, intent(out)     :: xCC(ixMlo^D:ixMhi^D,ndim)
 
    integer                           :: ixC^L, idims, level, ix
 
    level=node(plevel_,igrid)
 
    ! coordinates of cell centers
    xcc(ixm^t,:)=ps(igrid)%x(ixm^t,:)
 
    ! coordinates of cell corners
    ixcmin^d=ixmlo^d-1; ixcmax^d=ixmhi^d;
    if(slab_uniform)then
       do idims=1,ndim
         xc(ixc^s,idims)=ps(igrid)%x(ixc^s,idims)+0.5d0*dx(idims,level)
       end do
    else
       ! for any non-cartesian or stretched coordinate (allow multiple stretched directions)
       {do ix=ixcmin^d,ixcmax^d
         xc(ix^d%ixC^s,^d)=ps(igrid)%x(ix^d%ixC^s,^d)+0.5d0*ps(igrid)%dx(ix^d%ixC^s,^d)
       end do\}
    endif
 
  end subroutine calc_x
 
end module mod_calculate_xw