mod_hyperdiffusivity.t

Go to the documentation of this file.

!> Subroutines for Roe-type Riemann solver for HD
module mod_hyperdiffusivity
 
 
  implicit none
 
  !all public
  !private
  !public :: hyperdiffusivity_init
 
 
contains
 
  subroutine hyperdiffusivity_init()
    use mod_global_parameters
    use mod_geometry
 
    !print*, slab_uniform !!THIS IS FALSE, why??
 
    !if(coordinate .ne. Cartesian .or. .not. slab_uniform) then
    if(coordinate .ne. cartesian) then
      call mpistop("Hyperdiffusivity only implemented for Cartesian uniform grid")
    endif
 
    nghostcells = max(nghostcells,3)
 
  end subroutine hyperdiffusivity_init
 
 
  subroutine hyp_coeff(ixI^L, ixO^L, var, idimm, nu_hyp)
    use mod_global_parameters
    integer, intent(in)                :: ixI^L, idimm
    integer, intent(out)                :: ixO^L
    double precision, intent(in)       :: var(ixI^S)
    double precision, intent(out)       :: nu_hyp(ixI^S)
 
    double precision        :: tmp(ixI^S), tmp2(ixI^S)
    integer :: hx^L, hx1f^L, hx1b^L, hx2b^L
 
    double precision, parameter :: eps=1e-8
 
 
   hxmin^d=iximin^d+2*kr(idimm,^d); 
   hxmax^d=iximax^d-kr(idimm,^d); 
 
   hx1f^l=hx^l+kr(idimm,^d);  
   hx1b^l=hx^l-kr(idimm,^d);  
   hx2b^l=hx^l-2*kr(idimm,^d);  
 
  !store d3 in tmp
  tmp(hx^s)= abs(3d0 * (var(hx^s) - var(hx1b^s)) - (var(hx1f^s)-var(hx2b^s)))
  !store d1 in tmp2
  tmp2(hx^s)= abs((var(hx^s) - var(hx1b^s)))
 
  ixomin^d=hxmin^d+kr(idimm,^d);
  ixomax^d=hxmax^d-kr(idimm,^d);
  hx1f^l=ixo^l+kr(idimm,^d);  
  hx1b^l=ixo^l-kr(idimm,^d);  
 
  nu_hyp(ixo^s) = max(tmp(hx1b^s),tmp(ixo^s),tmp(hx1f^s))/max(tmp2(hx1b^s),tmp2(ixo^s),tmp2(hx1f^s),eps)  
 
  !print*, "HYP IXO ", ixO^L 
 
  end subroutine hyp_coeff
 
 
 
  subroutine div_vel_coeff(ixI^L, ixO^L, vel, idimm, nu_vel)
    use mod_global_parameters
    integer, intent(in)                :: ixI^L, idimm
    integer, intent(out)                :: ixO^L
    double precision, intent(in)       :: vel(ixI^S,1:ndir)
    double precision, intent(out)       :: nu_vel(ixI^S)
    integer :: hx1f^L, hx1b^L
 
    integer :: ii
 
   ixomin^d=iximin^d+1; 
   ixomax^d=iximax^d-1; 
 
   ixomin^d=ixomin^d+kr(idimm,^d); 
 
   nu_vel(ixo^s) = 0d0 
 
  ! ndim should always be smaller or equal to ndir !
   do ii = 1, ndim
    hx1b^l=ixo^l-kr(ii,^d);
    if(ii==idimm) then
      nu_vel(ixo^s) = nu_vel(ixo^s) + (vel(ixo^s,ii) - vel(hx1b^s,ii))/dxlevel(ii) 
    else
      hx1f^l=ixo^l+kr(ii,^d);  
      nu_vel(ixo^s) = nu_vel(ixo^s) + (vel(hx1f^s,ii) - vel(hx1b^s,ii))/(2d0*dxlevel(ii)) 
    endif
 
    where(nu_vel(ixo^s) < 0d0)
      nu_vel(ixo^s) = -nu_vel(ixo^s)
    elsewhere
      nu_vel(ixo^s) = 0d0
    endwhere
 
   enddo 
 
  !print*, "DIV_VEL IXO ", ixO^L 
 
  end subroutine div_vel_coeff
 
 
 
  !var has cell centered values
  !nu_hyper is defined at the interfaces
  subroutine second_same_deriv(ixI^L, ixO^L, nu_hyper, var, idimm, res)
    use mod_global_parameters
    integer, intent(in)                :: ixI^L, idimm
    integer, intent(out)                :: ixO^L
    double precision, intent(in)       :: nu_hyper(ixI^S),var(ixI^S)
    double precision, intent(out)       :: res(ixI^S)
 
    integer :: hxf^L, hxb^L
 
    ixomin^d=iximin^d+3;
    ixomax^d=iximax^d-3;
 
    hxf^l=ixo^l+kr(idimm,^d);
    hxb^l=ixo^l-kr(idimm,^d);
 
    res(ixo^s) = 1d0/(dxlevel(idimm)**2)*&
          (nu_hyper(hxf^s) * (var(hxf^s)-var(ixo^s))-&
          nu_hyper(ixo^s) * (var(ixo^s)-var(hxb^s)))
   !print*, "SECOND SAME DERIV IXO ", ixO^L 
 
  end subroutine second_same_deriv
 
  !var has cell centered values
  !var2 has cell centered values
  !nu_hyper is defined at the interfaces
  subroutine second_same_deriv2(ixI^L, ixO^L, nu_hyper, var2, var, idimm, res)
    use mod_global_parameters
    integer, intent(in)                :: ixI^L, idimm
    integer, intent(out)                :: ixO^L
    double precision, intent(in)       :: nu_hyper(ixI^S),var(ixI^S), var2(ixI^S)
    double precision, intent(out)       :: res(ixI^S)
 
    integer :: hxf^L, hxb^L
 
    ixomin^d=iximin^d+3;
    ixomax^d=iximax^d-3;
 
    hxf^l=ixo^l+kr(idimm,^d);
    hxb^l=ixo^l-kr(idimm,^d);
 
    res(ixo^s) = 1d0/(2d0*dxlevel(idimm)**2)*&
          (nu_hyper(hxf^s) * (var2(hxf^s)+var2(ixo^s)) *  (var(hxf^s)-var(ixo^s))-&
          nu_hyper(ixo^s) * (var2(hxb^s)+var2(ixo^s)) * (var(ixo^s)-var(hxb^s)))
    !print*, "SECOND SAME DERIV2 IXO ", ixO^L 
 
  end subroutine second_same_deriv2
 
  !idimm inner derivative, idimm2 outer
  ! deriv_idimm2(nu * deriv_idimm (u) )
  !var has cell centered values
  !nu_hyper is defined at the interfaces
  subroutine second_cross_deriv(ixI^L, ixO^L, nu_hyper, var, idimm, idimm2, res)
    use mod_global_parameters
    integer, intent(in)                :: ixI^L, idimm,idimm2
    integer, intent(out)                :: ixO^L
    double precision, intent(in)       :: nu_hyper(ixI^S),var(ixI^S)
    double precision, intent(out)       :: res(ixI^S)
 
    integer :: hxfi^L, hxbi^L, hxfifj^L, hxbibj^L, hxfibj^L, hxbifj^L
 
    ixomin^d=iximin^d+3;
    ixomax^d=iximax^d-3;
 
    hxfi^l=ixo^l+kr(idimm2,^d);
    hxbi^l=ixo^l-kr(idimm2,^d);
 
    hxfifj^l=hxfi^l+kr(idimm,^d);
    hxfibj^l=hxfi^l-kr(idimm,^d);
 
    hxbifj^l=hxbi^l+kr(idimm,^d);
    hxbibj^l=hxbi^l-kr(idimm,^d);
 
    res(ixo^s) = 1d0/(8d0*dxlevel(idimm) * dxlevel(idimm2))*&
          ((nu_hyper(hxfifj^s) + nu_hyper(hxfi^s)) * (var(hxfifj^s)-var(hxfibj^s))-&
          (nu_hyper(hxbifj^s) + nu_hyper(hxbi^s)) * (var(hxbifj^s)-var(hxbibj^s)))
 
    !print*, "SECOND CROSS DERIV IXO ", ixO^L 
 
  end subroutine second_cross_deriv
 
 
  !idimm inner derivative, idimm2 outer
  ! deriv_idimm2(nu * var2 * deriv_idimm (u) )
  !var has cell centered values
  !var2 has cell centered values
  !nu_hyper is defined at the interfaces (with numbering index left center): center_i-1 interface_i center_i
  subroutine second_cross_deriv2(ixI^L, ixO^L, nu_hyper, var2, var, idimm, idimm2, res)
    use mod_global_parameters
    integer, intent(in)                :: ixI^L, idimm,idimm2
    integer, intent(out)                :: ixO^L
    double precision, intent(in)       :: nu_hyper(ixI^S),var(ixI^S),var2(ixI^S)
    double precision, intent(out)       :: res(ixI^S)
 
    integer :: hxfi^L, hxbi^L, hxfifj^L, hxbibj^L, hxfibj^L, hxbifj^L
 
    ixomin^d=iximin^d+3;
    ixomax^d=iximax^d-3;
 
    hxfi^l=ixo^l+kr(idimm2,^d);
    hxbi^l=ixo^l-kr(idimm2,^d);
 
    hxfifj^l=hxfi^l+kr(idimm,^d);
    hxfibj^l=hxfi^l-kr(idimm,^d);
 
    hxbifj^l=hxbi^l+kr(idimm,^d);
    hxbibj^l=hxbi^l-kr(idimm,^d);
 
    res(ixo^s) = 1d0/(8d0*dxlevel(idimm) * dxlevel(idimm2))*&
          ((nu_hyper(hxfifj^s) + nu_hyper(hxfi^s)) * var2(hxfi^s) * (var(hxfifj^s)-var(hxfibj^s))-&
          (nu_hyper(hxbifj^s) + nu_hyper(hxbi^s)) * var2(hxbi^s) * (var(hxbifj^s)-var(hxbibj^s)))
 
    !print*, "SECOND CROSS DERIV2 IXO ", ixO^L 
 
  end subroutine second_cross_deriv2
 
end module mod_hyperdiffusivity