mod_physics_mhd.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_mhd
  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 /mhd_list/ mhd_energy, mhd_gamma, mhd_glm_alpha, mhd_gravity,&
      mhd_n_tracer, mhd_radiative_cooling, He_abundance, mhd_eta, mhd_source_usr, &
      mhd_resistivity, mhd_hyperbolic_thermal_conduction, &
      mhd_hyperbolic_thermal_conduction_anisotropic, &
      tc_kappa_par, tc_kappa_perp, &
      mhd_energy_only

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

    
    

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

    integer      :: idir, idum

    call read_params(par_files)
    call phys_units()

    phys_energy  = mhd_energy
    phys_total_energy  = mhd_energy
    phys_internal_e = .false.
    phys_gamma = mhd_gamma
    phys_partial_ionization=mhd_partial_ionization
    need_global_cmax=.true.
    gamma_1=mhd_gamma-1.0_dp
    inv_gamma_1=1.0_dp/gamma_1
 

    use_particles = mhd_particles

    ! Determine flux variables
    rho_ = var_set_rho()
    

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

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

    ! Set index for heat flux variable(s)
    allocate(q_(ndir))
    q_ = var_set_qvec(ndir, need_bc=.false.)
    

    allocate(mag(ndir))
    mag(:) = var_set_bfield(ndir)
    

    psi_ = var_set_fluxvar('psi', 'psi', need_bc=.false.)
    

    !> GLM MHD uses split source addition in psi:
    any_source_split = .true.
    

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

    ! Register tracer fields
        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

    ! BnormLF fix:
    flux_type(:,psi_)=flux_tvdlf
    do idir=1,ndir
       flux_type(idir,mag(idir))=flux_tvdlf
    end do
    

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


! use cycle, needs to be dealt with:
!    ! Initialize particles module
!    if (mhd_particles) then
!       call particles_init()
!       phys_req_diagonal = .true.
!    end if

    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

    

    if (tc_kappa0_par <= 0.0d0 .and. tc_kappa_par <= 0.0d0) &
      tc_kappa0_par = 8.0d-7 * unit_temperature**3.5d0 &
                    / (unit_length * unit_density * unit_velocity**3.0d0)
    
    
    htc_Cchi = 0.823d0 * (4.753567596681522d6 / 20.0d0) &
             * unit_magneticfield * unit_temperature**1.5d0 / unit_numberdensity
    
    
  end subroutine phys_units
  pure subroutine to_primitive(u)
    
    real(dp), intent(inout) :: u(nw_phys)

    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)
    u(iw_e)=gamma_1*(u(iw_e)-0.5_dp*&
      (u(iw_rho)*(u(iw_mom(1))**2+u(iw_mom(2))**2+u(iw_mom(3))**2)+&
       u(iw_mag(1))**2+u(iw_mag(2))**2+u(iw_mag(3))**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)*(u(iw_mom(1))**2+u(iw_mom(2))**2+u(iw_mom(3))**2)+&
       u(iw_mag(1))**2+u(iw_mag(2))**2+u(iw_mag(3))**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_mhd