mod_physics_hd.f90 Source File

!> This module defines the procedures of a physics module. It contains function
!> pointers for the various supported routines. An actual physics module has to
!> set these pointers to its implementation of these routines.


module mod_physics_hd
  use mod_variables
  use mod_physics_vars
  use mod_comm_lib, only: mpistop
  
  implicit none
  public
  

  procedure(sub_check_params), pointer    :: phys_check_params           => &
     null()
  procedure(sub_boundary_adjust), pointer :: phys_boundary_adjust        => &
     null()
  procedure(sub_global_source), pointer   :: phys_global_source_after    => &
     null()
  procedure(sub_implicit_update), pointer :: phys_implicit_update        => &
     null()
  procedure(sub_evaluate_implicit),pointer:: phys_evaluate_implicit      => &
     null()
  procedure(sub_clean_divb), pointer      :: phys_clean_divb             => &
     null()
  procedure(sub_special_advance), pointer :: phys_special_advance        => &
     null()
  ! set the equilibrium variables
  procedure(sub_set_equi_vars), pointer   :: phys_set_equi_vars          => &
     null()
  ! subroutine with no parameters which creates EUV images
  procedure(sub_check_params), pointer    :: phys_te_images              => &
     null()
  procedure(sub_face_to_center), pointer  :: phys_face_to_center         => &
     null()
  procedure(sub_write_info), pointer      :: phys_write_info             => &
     null()

  interface
     
     subroutine sub_check_params
     end subroutine sub_check_params

     subroutine sub_boundary_adjust(igrid,psb)
       use mod_global_parameters
       integer, intent(in) :: igrid
       type(state), target :: psb(max_blocks)
     end subroutine sub_boundary_adjust
     
     subroutine sub_face_to_center(ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,&
        ixOmax3,s)
       use mod_global_parameters
       integer, intent(in)                :: ixOmin1,ixOmin2,ixOmin3,ixOmax1,&
          ixOmax2,ixOmax3
       type(state)                        :: s
     end subroutine sub_face_to_center

     subroutine sub_write_info(file_handle)
       integer, intent(in) :: file_handle
     end subroutine sub_write_info

     !> Add special advance in each advect step
     subroutine sub_special_advance(qt, psa)
       use mod_global_parameters
       double precision, intent(in) :: qt
       !< Current time
       type(state), target :: psa(max_blocks)
       !< Compute based on this state
     end subroutine sub_special_advance

     !> Add global source terms on complete domain (potentially implicit)
     subroutine sub_global_source(qdt, qt, active)
       use mod_global_parameters
       double precision, intent(in) :: qdt
       !< Current time step
       double precision, intent(in) :: qt
       !< Current time
       logical, intent(inout)       :: active
       !< Output if the source is active
     end subroutine sub_global_source

     subroutine sub_evaluate_implicit(qtC,psa)
       use mod_global_parameters
       type(state), target :: psa(max_blocks)
       double precision, intent(in) :: qtC
     end subroutine sub_evaluate_implicit

     subroutine sub_implicit_update(dtfactor,qdt,qtC,psa,psb)
       use mod_global_parameters
       type(state), target :: psa(max_blocks)
       type(state), target :: psb(max_blocks)
       double precision, intent(in) :: qdt
       double precision, intent(in) :: qtC
       double precision, intent(in) :: dtfactor
     end subroutine sub_implicit_update
     
     !> set equilibrium variables other than b0 (e.g. p0 and rho0)
     subroutine sub_set_equi_vars(igrid)
       integer, intent(in) :: igrid
     end subroutine sub_set_equi_vars

     !> clean initial divb
     subroutine sub_clean_divb(qdt, qt, active)
       use mod_global_parameters
       double precision, intent(in) :: qdt
       !< Current time step
       double precision, intent(in) :: qt
       !< Current time
       logical, intent(inout)       :: active
       !< Output if the source is active
     end subroutine sub_clean_divb

  end interface

contains

    !> Read this module's parameters from a file
  subroutine read_params(files)
    use mod_global_parameters
    character(len=*), intent(in) :: files(:)
    integer                      :: n

    namelist /hd_list/ hd_energy, hd_gamma, hd_adiab, hd_partial_ionization,&
        hd_force_diagonal, hd_particles, hd_gravity, hd_n_tracer, hd_radiative_cooling, He_abundance

    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status="old")
       read(unitpar, hd_list, end=111)
111    close(unitpar)
    end do

    

  end subroutine read_params
    !> Initialize the module
  subroutine phys_init()
    use mod_global_parameters
!    use mod_particles, only: particles_init
    use mod_radiative_cooling, only: rc_fl, radiative_cooling_init_params, radiative_cooling_init

    call read_params(par_files)
    call phys_units()

    phys_energy  = hd_energy
    phys_total_energy  = hd_energy
    phys_internal_e = .false.
    phys_gamma = hd_gamma
    phys_partial_ionization=hd_partial_ionization
    gamma_1 = hd_gamma - 1.0d0
    inv_gamma_1 = 1.0d0/gamma_1
    

    use_particles = hd_particles

    ! Determine flux variables
    rho_ = var_set_rho()
    

    allocate(mom(ndir))
    mom(:) = var_set_momentum(ndir)
    

    ! Set index of energy variable
    if (hd_energy) then
       e_ = var_set_energy()
       p_ = e_
    else
       e_ = -1
       p_ = -1
    end if
    

    ! Whether diagonal ghost cells are required for the physics
    phys_req_diagonal = .false.

    if (hd_force_diagonal) then
       ! ensure corners are filled, otherwise divide by zero when getting primitives
       !  --> only for debug purposes
       phys_req_diagonal = .true.
    endif

        tracer(1) = var_set_fluxvar("trc", "trp", 1, need_bc=.false.)
    

    ! set number of variables which need update ghostcells
    nwgc=nwflux
    

    ! Define custom flux types:
    if (.not. allocated(flux_type)) then
       allocate(flux_type(ndir, nw_flux))
       flux_type = flux_default
    else if (any(shape(flux_type) /= [ndir, nw_flux])) then
       call mpistop("phys_check error: flux_type has wrong shape")
    end if
    
    
! use cycle, needs to be dealt with:    
!    ! Initialize particles module
!    if (hd_particles) then
!       call particles_init()
!       phys_req_diagonal = .true.
!    end if

    nvector      = 1 ! No. vector vars
    allocate(iw_vector(nvector))
    iw_vector(1) = mom(1) - 1
    
    

    call radiative_cooling_init_params(phys_gamma,He_abundance)
    call radiative_cooling_init(rc_fl)
    
    

  end subroutine phys_init
  subroutine phys_units()
    use mod_global_parameters
    double precision :: mp, kB
    double precision :: a,b

    !> here no SI_UNIT used by default, to be implemented
    mp = mp_cgs
    kB = kB_cgs
    !> eq_state_units by default, to be implemented
    a = 1.d0+4.d0*He_abundance
    b = 2.d0+3.d0*He_abundance

    if(unit_density/=1.d0 .or. unit_numberdensity/=1.d0) then
      if(unit_density/=1.d0) then
        unit_numberdensity=unit_density/(a*mp)
      else if(unit_numberdensity/=1.d0) then
        unit_density=a*mp*unit_numberdensity
      end if
      if(unit_temperature/=1.d0) then
        unit_pressure=b*unit_numberdensity*kB*unit_temperature
        unit_velocity=dsqrt(unit_pressure/unit_density)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      else if(unit_pressure/=1.d0) then
        unit_temperature=unit_pressure/(b*unit_numberdensity*kB)
        unit_velocity=dsqrt(unit_pressure/unit_density)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      else if(unit_velocity/=1.d0) then
        unit_pressure=unit_density*unit_velocity**2
        unit_temperature=unit_pressure/(b*unit_numberdensity*kB)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      else if(unit_time/=1.d0) then
        unit_velocity=unit_length/unit_time
        unit_pressure=unit_density*unit_velocity**2
        unit_temperature=unit_pressure/(b*unit_numberdensity*kB)
      end if
    else if(unit_temperature/=1.d0) then
      ! units of temperature and velocity are dependent
      if(unit_pressure/=1.d0) then
        unit_numberdensity=unit_pressure/(b*unit_temperature*kB)
        unit_density=a*mp*unit_numberdensity
        unit_velocity=dsqrt(unit_pressure/unit_density)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      end if
    else if(unit_pressure/=1.d0) then
      if(unit_velocity/=1.d0) then
        unit_density=unit_pressure/unit_velocity**2
        unit_numberdensity=unit_density/(a*mp)
        unit_temperature=unit_pressure/(b*unit_numberdensity*kB)
        if(unit_length/=1.d0) then
          unit_time=unit_length/unit_velocity
        else if(unit_time/=1.d0) then
          unit_length=unit_velocity*unit_time
        end if
      else if(unit_time/=0.d0) then
        unit_velocity=unit_length/unit_time
        unit_density=unit_pressure/unit_velocity**2
        unit_numberdensity=unit_density/(a*mp)
        unit_temperature=unit_pressure/(b*unit_numberdensity*kB)
      end if
    end if
    unit_mass=unit_density*unit_length**3

    
  end subroutine phys_units
pure subroutine to_primitive(u)
  
  real(dp), intent(inout) :: u(nw_phys)

  ! Compute velocity from momentum
      u(iw_mom(1)) = u(iw_mom(1))/u(iw_rho)
      u(iw_mom(2)) = u(iw_mom(2))/u(iw_rho)
      u(iw_mom(3)) = u(iw_mom(3))/u(iw_rho)

  ! Compute pressure from energy
  u(iw_e) = gamma_1 * (u(iw_e) - 0.5_dp * u(iw_rho) * &
     sum(u(iw_mom(1:ndim))**2) )

end subroutine to_primitive
pure subroutine to_conservative(u)
  
  real(dp), intent(inout) :: u(nw_phys)

  ! Compute energy from pressure and kinetic energy
  u(iw_e) = u(iw_e) * inv_gamma_1 + 0.5_dp * u(iw_rho) * &
     sum(u(iw_mom(1:ndim))**2)

  ! Compute momentum from density and velocity components
      u(iw_mom(1)) = u(iw_rho) * u(iw_mom(1))
      u(iw_mom(2)) = u(iw_rho) * u(iw_mom(2))
      u(iw_mom(3)) = u(iw_rho) * u(iw_mom(3))

end subroutine to_conservative
  subroutine phys_activate()
    call phys_init()
  end subroutine phys_activate

  
  subroutine phys_check()

    use mod_comm_lib, only: mpistop

    if (physics_type == "") call mpistop("Error: no physics module loaded")

    ! Checks whether the required physics methods have been defined
    if (.not. associated(phys_check_params)) phys_check_params => &
       dummy_check_params
    
    if (.not. associated(phys_boundary_adjust)) phys_boundary_adjust => &
       dummy_boundary_adjust

    if (.not. associated(phys_evaluate_implicit)) phys_evaluate_implicit => &
       dummy_evaluate_implicit

    if (.not. associated(phys_implicit_update)) phys_implicit_update => &
       dummy_implicit_update

    if (.not. associated(phys_face_to_center)) phys_face_to_center => &
       dummy_face_to_center

    if (.not. associated(phys_write_info)) phys_write_info => dummy_write_info

  end subroutine phys_check

  !> Transform primitive variables into conservative ones
  subroutine phys_to_conserved(ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,&
      ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3, w, x)
    use mod_global_parameters
    integer, intent(in)             :: ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
       ixImax3, ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3
    double precision, intent(inout) :: w(ixImin1:ixImax1,ixImin2:ixImax2,&
         ixImin3:ixImax3, nw_phys) 
    double precision, intent(in)    :: x(ixImin1:ixImax1,ixImin2:ixImax2,&
         ixImin3:ixImax3, 1:ndim)
    integer                         :: ix1,ix2,ix3

    do ix3=ixOmin3,ixOmax3
       do ix2=ixOmin2,ixOmax2
          do ix1=ixOmin1,ixOmax1

             call to_conservative(w(ix1,ix2,ix3,:))

          end do
       end do
    end do

  end subroutine phys_to_conserved
  
  !> Transform conservative variables into primitive ones
  subroutine phys_to_primitive(ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,&
      ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3, w, x)
    use mod_global_parameters
    integer, intent(in)             :: ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
       ixImax3, ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3
    double precision, intent(inout) :: w(ixImin1:ixImax1,ixImin2:ixImax2,&
       ixImin3:ixImax3, nw_phys)
    double precision, intent(in)    :: x(ixImin1:ixImax1,ixImin2:ixImax2,&
       ixImin3:ixImax3, 1:ndim)
    integer                         :: ix1,ix2,ix3
    
    do ix3=ixOmin3,ixOmax3
       do ix2=ixOmin2,ixOmax2
          do ix1=ixOmin1,ixOmax1

             call to_primitive(w(ix1,ix2,ix3,:))

          end do
       end do
    end do

  end subroutine phys_to_primitive
  
  subroutine dummy_boundary_adjust(igrid,psb)
    use mod_global_parameters
    integer, intent(in) :: igrid
    type(state), target :: psb(max_blocks)
  end subroutine dummy_boundary_adjust
  
  subroutine dummy_evaluate_implicit(qtC,psa)
    use mod_global_parameters
    type(state), target :: psa(max_blocks)
    double precision, intent(in) :: qtC
    integer :: iigrid, igrid

    ! Just copy in nul state
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       psa(igrid)%w = 0.0d0*psa(igrid)%w
       if(stagger_grid) psa(igrid)%ws = 0.0d0*psa(igrid)%ws
    end do

  end subroutine dummy_evaluate_implicit

  subroutine dummy_implicit_update(dtfactor,qdt,qtC,psa,psb)
    use mod_global_parameters
    type(state), target :: psa(max_blocks)
    type(state), target :: psb(max_blocks)
    double precision, intent(in) :: qdt
    double precision, intent(in) :: qtC
    double precision, intent(in) :: dtfactor
    integer :: iigrid, igrid

    ! Just copy in psb state when using the scheme without implicit part
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       psa(igrid)%w = psb(igrid)%w
       if(stagger_grid) psa(igrid)%ws = psb(igrid)%ws
    end do

  end subroutine dummy_implicit_update

  subroutine dummy_check_params
  end subroutine dummy_check_params

  subroutine dummy_face_to_center(ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,&
     ixOmax3,s)
    use mod_global_parameters
    integer, intent(in)                :: ixOmin1,ixOmin2,ixOmin3,ixOmax1,&
       ixOmax2,ixOmax3
    type(state)                        :: s
  end subroutine dummy_face_to_center

  subroutine dummy_write_info(fh)
    use mod_global_parameters
    integer, intent(in)                 :: fh
    !< File handle
    integer, dimension(MPI_STATUS_SIZE) :: st
    integer                             :: er

    ! Number of physics parameters
    integer, parameter                  :: n_par = 0

    call MPI_FILE_WRITE(fh, n_par, 1, MPI_INTEGER, st, er)
  end subroutine dummy_write_info
  
end module mod_physics_hd