mod_nonlinear_phys.t

Go to the documentation of this file.

!> Module containing the physics routines for scalar nonlinear equation
module mod_nonlinear_phys
 
  implicit none
  private
 
  !> index of the single scalar unknown
  integer, protected, public :: rho_       = 1
 
  !> switch between burgers (i.e. rho**2) 
  !> or nonconvex flux (i.e. rho**3)
  integer, protected, public :: nonlinear_flux_type = 1
 
  !> whether the KdV source term is added
  logical, protected, public :: kdv_source_term = .false.
 
  !> Whether particles module is added
  logical, public, protected              :: nonlinear_particles = .false.
 
  ! Public methods
  public :: nonlinear_phys_init
  public :: nonlinear_get_v
 
contains
 
  !> Read this module's parameters from a file
  subroutine nonlinear_params_read(files)
    use mod_global_parameters, only: unitpar
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /nonlinear_list/ nonlinear_flux_type, kdv_source_term, nonlinear_particles
 
    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status='old')
       read(unitpar, nonlinear_list, end=111)
111    close(unitpar)
    end do
 
  end subroutine nonlinear_params_read
 
  !> Write this module's parameters to a snapshot
  subroutine nonlinear_write_info(fh)
  ! for nonlinear scalar equation, nothing to write
  ! note: this is info only stored at end of dat files, 
  !       is never read/used for restarts, only expects
  !       an integer (number of parameters) and 
  !       corresponding double values and character names
  !       and is meant for use in the python tools
    use mod_global_parameters
    integer, intent(in)                 :: fh
    integer, dimension(MPI_STATUS_SIZE) :: st
    integer                             :: er
 
    ! Write zero parameters
    call mpi_file_write(fh, 0, 1, mpi_integer, st, er)
 
  end subroutine nonlinear_write_info
 
  subroutine nonlinear_phys_init()
    use mod_global_parameters
    use mod_physics
    use mod_kdv, only: kdv_init
    use mod_particles, only: particles_init
 
    call nonlinear_params_read(par_files)
 
    physics_type = "nonlinear"
    phys_energy  = .false.
    use_particles = nonlinear_particles
 
    allocate(start_indices(number_species),stop_indices(number_species))
    ! set the index of the first flux variable for species 1
    start_indices(1)=1
    rho_ = var_set_rho()
 
    ! set number of variables which need update ghostcells
    nwgc=nwflux
 
    ! set the index of the last flux variable for species 1
    stop_indices(1)=nwflux
  
    !  Number of variables need reconstruction in w
    nw_recon=nwflux
 
    ! Check whether custom flux types have been defined
    if (.not. allocated(flux_type)) then
       allocate(flux_type(ndir, nw))
       flux_type = flux_default
    else if (any(shape(flux_type) /= [ndir, nw])) then
       call mpistop("phys_check error: flux_type has wrong shape")
    end if
 
    phys_get_cmax        => nonlinear_get_cmax
    phys_get_cbounds     => nonlinear_get_cbounds
    phys_get_flux        => nonlinear_get_flux
    phys_add_source_geom => nonlinear_add_source_geom
    phys_add_source      => nonlinear_add_source
    phys_to_conserved    => nonlinear_to_conserved
    phys_to_primitive    => nonlinear_to_primitive
    phys_get_dt          => nonlinear_get_dt
    phys_write_info      => nonlinear_write_info
 
    if (kdv_source_term) call kdv_init()
 
    ! Initialize particles module
    if (nonlinear_particles) then
       call particles_init()
    end if
 
  end subroutine nonlinear_phys_init
 
  subroutine nonlinear_to_conserved(ixI^L, ixO^L, w, x)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:^nd)
 
    ! Do nothing (primitive and conservative are equal for nonlinear module)
  end subroutine nonlinear_to_conserved
 
  subroutine nonlinear_to_primitive(ixI^L, ixO^L, w, x)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(inout) :: w(ixi^s, nw)
    double precision, intent(in)    :: x(ixi^s, 1:^nd)
 
    ! Do nothing (primitive and conservative are equal for nonlinear module)
  end subroutine nonlinear_to_primitive
 
  subroutine nonlinear_get_v(w, x, ixI^L, ixO^L, idim, v)
    use mod_global_parameters
    integer, intent(in)           :: ixi^l, ixo^l, idim
    double precision, intent(in)  :: w(ixi^s, nw), x(ixi^s, 1:^nd)
    double precision, intent(out) :: v(ixi^s)
 
   select case(nonlinear_flux_type)
    case(1)
       v(ixo^s)=w(ixo^s,rho_)
    case(2)
       v(ixo^s)=3.0d0*w(ixo^s,rho_)**2
    case default
       call mpistop('Undefined fluxtype: set nonlinear_flux_type to 1 or 2')
    end select
 
  end subroutine nonlinear_get_v
 
  subroutine nonlinear_get_cmax(w, x, ixI^L, ixO^L, idim, cmax)
    use mod_global_parameters
    integer, intent(in)                       :: ixi^l, ixo^l, idim
    double precision, intent(in)              :: w(ixi^s, nw), x(ixi^s, 1:^nd)
    double precision, intent(inout)           :: cmax(ixi^s)
 
    call nonlinear_get_v(w, x, ixi^l, ixo^l, idim, cmax)
 
    cmax(ixo^s) = abs(cmax(ixo^s))
 
  end subroutine nonlinear_get_cmax
 
  subroutine nonlinear_get_cbounds(wLC, wRC, wLp, wRp, x, ixI^L, ixO^L, idim,Hspeed, cmax, cmin)
    use mod_global_parameters
    use mod_variables
    integer, intent(in)             :: ixi^l, ixo^l, idim
    double precision, intent(in)    :: wlc(ixi^s, nw), wrc(ixi^s,nw)
    double precision, intent(in)    :: wlp(ixi^s, nw), wrp(ixi^s,nw)
    double precision, intent(in)    :: x(ixi^s, 1:^nd)
    double precision, intent(inout) :: cmax(ixi^s,1:number_species)
    double precision, intent(inout), optional :: cmin(ixi^s,1:number_species)
    double precision, intent(in)    :: hspeed(ixi^s,1:number_species)
    double precision :: wmean(ixi^s,nw)
 
    ! since get_v depends on w, the first argument should be some average over the
    ! left and right state
    wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
    call nonlinear_get_v(wmean, x, ixi^l, ixo^l, idim, cmax(ixi^s,1))
 
    if (present(cmin)) then
       cmin(ixo^s,1) = min(cmax(ixo^s,1), zero)
       cmax(ixo^s,1) = max(cmax(ixo^s,1), zero)
    else
       cmax(ixo^s,1) = maxval(abs(cmax(ixo^s,1)))
    end if
 
  end subroutine nonlinear_get_cbounds
 
  subroutine nonlinear_get_dt(w, ixI^L, ixO^L, dtnew, dx^D, x)
    use mod_global_parameters
    use mod_kdv, only: kdv_get_dt
 
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: dx^d, x(ixi^s, 1:^nd)
    double precision, intent(in)    :: w(ixi^s, 1:nw)
    double precision, intent(inout) :: dtnew
 
    dtnew = bigdouble
 
    if(kdv_source_term) then
      call kdv_get_dt(w, ixi^l, ixo^l, dtnew, dx^d, x)
    endif
  end subroutine nonlinear_get_dt
 
  ! here we select the flux according to the nonlinear_flux_type parameter
  subroutine nonlinear_get_flux(wC, w, x, ixI^L, ixO^L, idim, f)
    use mod_global_parameters
    integer, intent(in)             :: ixi^l, ixo^l, idim
    double precision, intent(in)    :: wc(ixi^s, 1:nw)
    double precision, intent(in)    :: w(ixi^s, 1:nw)
    double precision, intent(in)    :: x(ixi^s, 1:^nd)
    double precision, intent(out)   :: f(ixi^s, nwflux)
 
    select case(nonlinear_flux_type)
    case(1)
       f(ixo^s,rho_)=half*w(ixo^s,rho_)**2
    case(2)
       f(ixo^s,rho_)=w(ixo^s,rho_)**3
    case default
       call mpistop('Undefined fluxtype: set nonlinear_flux_type to 1 or 2')
    end select
 
  end subroutine nonlinear_get_flux
 
  subroutine nonlinear_add_source_geom(qdt, dtfactor, ixI^L, ixO^L, wCT,wprim, w, x)
 
    ! Add geometrical source terms to w
    ! There are no geometrical source terms 
 
    use mod_global_parameters
 
    integer, intent(in) :: ixi^l, ixo^l
    double precision, intent(in) :: qdt, dtfactor, x(ixi^s, 1:^nd)
    double precision, intent(inout) :: wct(ixi^s, 1:nw), wprim(ixi^s,1:nw),w(ixi^s, 1:nw)
 
  end subroutine nonlinear_add_source_geom
 
  subroutine nonlinear_add_source(qdt, dtfactor, ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
    use mod_global_parameters
    use mod_kdv, only: kdv_add_source
 
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: qdt,dtfactor
    double precision, intent(in)    :: wct(ixi^s, 1:nw),wctprim(ixi^s,1:nw), x(ixi^s, 1:ndim)
    double precision, intent(inout) :: w(ixi^s, 1:nw)
    logical, intent(in)             :: qsourcesplit
    logical, intent(inout)          :: active
 
    if(kdv_source_term) then
      call kdv_add_source(qdt,ixi^l,ixo^l,wct,w,x,qsourcesplit,active)
    end if
 
  end subroutine nonlinear_add_source
 
end module mod_nonlinear_phys