!> 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_srhd 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 /srhd_list/ srhd_eos,srhd_gamma,srhd_n_tracer, & He_abundance, srhd_source_usr, & srhd_small_pressure, srhd_small_density do n = 1, size(files) open(unitpar, file=trim(files(n)), status="old") read(unitpar, srhd_list, end=111) 111 close(unitpar) end do if (srhd_small_pressure < 0.0d0) call mpistop(& "srhd_small_pressure should be positive.") if (srhd_small_density < 0.0d0) call mpistop(& "srhd_small_density should be positive.") end subroutine read_params !> Initialize the module subroutine phys_init() use mod_global_parameters integer :: idir call read_params(par_files) call phys_units() phys_energy = .true. phys_total_energy = .true. phys_gamma = srhd_gamma gamma_1=srhd_gamma-1.0d0 inv_gamma_1=1.0d0/gamma_1 gamma_to_gamma_1=srhd_gamma/gamma_1 phys_internal_e=.false. phys_partial_ionization=.false. need_global_cmax=.false. ! Whether diagonal ghost cells are required for the physics phys_req_diagonal = .false. use_particles = srhd_particles ! Determine flux variables rho_ = var_set_rho() d_=rho_ allocate(mom(ndir)) mom(:) = var_set_momentum(ndir) ! Set index of energy variable e_ = var_set_energy() p_ = e_ ! Register tracer fields tracer(1) = var_set_fluxvar("trc", "trp", 1, need_bc=.false.) ! Set index for auxiliary variables ! MUST be after the possible tracers (which have fluxes) xi_ = var_set_auxvar('xi','xi') lfac_= var_set_auxvar('lfac','lfac') ! set number of variables which need update ghostcells nwgc=nwflux+nwaux ! 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 nvector = 1 ! No. vector vars allocate(iw_vector(nvector)) iw_vector(1) = mom(1) - 1 ! use cycle, needs to be dealt with: ! ! Initialize particles module ! if (srhd_particles) then ! call particles_init() ! phys_req_diagonal = .true. ! end if end subroutine phys_init subroutine phys_units() use mod_global_parameters double precision :: mp,kB !> here no SI_UNIT used by default, to be implemented mp = mp_cgs kB = kB_cgs unit_velocity=const_c ! we assume user sets: unit_numberdensity, unit_length, He_abundance ! then together with light speed c, all units fixed unit_density=(1.0d0+4.0d0*He_abundance)*mp*unit_numberdensity unit_pressure=unit_density*unit_velocity**2 unit_temperature=unit_pressure/((2.0d0+3.0d0*He_abundance)*unit_numberdensity*kB) unit_time=unit_length/unit_velocity unit_mass=unit_density*unit_length**3 end subroutine phys_units pure subroutine to_primitive(u) use mod_con2prim real(dp), intent(inout) :: u(nw_phys) real(dp) :: rho,rhoh,pth,E real(dp) :: ssqr !! begin: call srhd_get_auxiliary(u) ssqr=(u(iw_mom(1))**2+u(iw_mom(2))**2+u(iw_mom(3))**2) if(srhd_eos)then call con2prim_eos(u(lfac_),u(xi_),u(iw_rho),ssqr,u(iw_e)) else call con2prim(u(lfac_),u(xi_),u(iw_rho),ssqr,u(iw_e)) endif !! end: call srhd_get_auxiliary(u) rho=u(iw_rho)/u(lfac_) rhoh=u(xi_)/u(lfac_)**2 !! begin: call srhd_get_pressure_eos(rho,rhoh,pth,E) if(srhd_eos)then E = (rhoh+dsqrt(rhoh**2+(srhd_gamma**2-1.0_dp)*rho**2)) & /(srhd_gamma+1.0_dp) pth = 0.5d0*gamma_1* (E-rho*(rho/E)) else pth =(rhoh-rho)/gamma_to_gamma_1 endif !! end: call srhd_get_pressure_eos(rho,rhoh,pth,E) u(iw_rho)=rho u(iw_mom(1))=u(lfac_)*u(iw_mom(1))/u(xi_) u(iw_mom(2))=u(lfac_)*u(iw_mom(2))/u(xi_) u(iw_mom(3))=u(lfac_)*u(iw_mom(3))/u(xi_) u(iw_e)=pth end subroutine to_primitive pure subroutine to_conservative(u) real(dp), intent(inout) :: u(nw_phys) real(dp) :: rho,rhoh,pth real(dp) :: E,E_th rhoh=(u(iw_mom(1))**2+u(iw_mom(2))**2+u(iw_mom(3))**2) u(lfac_)=dsqrt(1.0d0+rhoh) rho=u(iw_rho) pth=u(iw_e) !! begin: call srhd_get_enthalpy_eos(rho,pth,rhoh) if(srhd_eos) then E_th = pth*inv_gamma_1 E = E_th+dsqrt(E_th**2+rho**2) ! writing rho/E on purpose, for numerics rhoh = 0.5_dp*((srhd_gamma+1.0_dp)*E & - gamma_1*rho*(rho/E)) else rhoh = rho+gamma_to_gamma_1*pth end if !! end: call srhd_get_enthalpy_eos(rho,pth,rhoh) rhoh=rhoh*u(lfac_) u(xi_)=u(lfac_)*rhoh u(iw_rho)=u(lfac_)*rho u(iw_mom(1))=rhoh*u(iw_mom(1)) u(iw_mom(2))=rhoh*u(iw_mom(2)) u(iw_mom(3))=rhoh*u(iw_mom(3)) u(iw_e)=u(xi_)-pth-u(iw_rho) 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_srhd