diff --git a/sources/equations.F90 b/sources/equations.F90 index c6dbb3e..1d438fd 100644 --- a/sources/equations.F90 +++ b/sources/equations.F90 @@ -53,13 +53,19 @@ module equations integer(kind=4), save :: nf = 0 integer(kind=4), save :: ns = 0 +! subroutine interfaces +! + interface prim2cons + module procedure prim2cons_cell + module procedure prim2cons_stencil + end interface + ! interfaces for procedure pointers ! interface - subroutine prim2cons_iface(q, u, s) - real(kind=8), dimension(:,:), intent(in) :: q - real(kind=8), dimension(:,:), intent(out) :: u - logical , optional , intent(in) :: s + subroutine prim2cons_iface(q, u) + real(kind=8), dimension(:), intent(in) :: q + real(kind=8), dimension(:), intent(out) :: u end subroutine subroutine cons2prim_iface(u, q, s) real(kind=8), dimension(:,:), intent(in) :: u @@ -84,12 +90,12 @@ module equations ! pointers to the conversion procedures ! - procedure(prim2cons_iface) , pointer, save :: prim2cons => null() - procedure(cons2prim_iface) , pointer, save :: cons2prim => null() + procedure(prim2cons_iface) , pointer, save :: prim2cons_ptr => null() + procedure(cons2prim_iface) , pointer, save :: cons2prim => null() ! pointer to the flux procedure ! - procedure(fluxspeed_iface) , pointer, save :: fluxspeed => null() + procedure(fluxspeed_iface) , pointer, save :: fluxspeed => null() procedure(get_maximum_speeds_iface), pointer, save :: & get_maximum_speeds => null() @@ -205,13 +211,7 @@ module equations ! declare public variables and subroutines ! public :: initialize_equations, finalize_equations, print_equations - public :: cons2prim, prim2cons, fluxspeed - public :: prim2cons_hd_iso, fluxspeed_hd_iso - public :: prim2cons_hd_adi, fluxspeed_hd_adi - public :: prim2cons_mhd_iso, fluxspeed_mhd_iso - public :: prim2cons_mhd_adi, fluxspeed_mhd_adi - public :: prim2cons_srhd_adi, fluxspeed_srhd_adi - public :: prim2cons_srmhd_adi, fluxspeed_srmhd_adi + public :: prim2cons, cons2prim, fluxspeed public :: get_maximum_speeds public :: update_primitive_variables public :: fix_unphysical_cells, correct_unphysical_states @@ -350,7 +350,7 @@ module equations ! set pointers to subroutines ! - prim2cons => prim2cons_hd_iso + prim2cons_ptr => prim2cons_hd_iso cons2prim => cons2prim_hd_iso fluxspeed => fluxspeed_hd_iso get_maximum_speeds => get_maximum_speeds_hd_iso @@ -389,7 +389,7 @@ module equations ! set pointers to subroutines ! - prim2cons => prim2cons_hd_adi + prim2cons_ptr => prim2cons_hd_adi cons2prim => cons2prim_hd_adi fluxspeed => fluxspeed_hd_adi get_maximum_speeds => get_maximum_speeds_hd_adi @@ -514,7 +514,7 @@ module equations ! set pointers to the subroutines ! - prim2cons => prim2cons_mhd_iso + prim2cons_ptr => prim2cons_mhd_iso cons2prim => cons2prim_mhd_iso fluxspeed => fluxspeed_mhd_iso get_maximum_speeds => get_maximum_speeds_mhd_iso @@ -560,7 +560,7 @@ module equations ! set pointers to subroutines ! - prim2cons => prim2cons_mhd_adi + prim2cons_ptr => prim2cons_mhd_adi cons2prim => cons2prim_mhd_adi fluxspeed => fluxspeed_mhd_adi get_maximum_speeds => get_maximum_speeds_mhd_adi @@ -688,7 +688,7 @@ module equations ! set pointers to subroutines ! - prim2cons => prim2cons_srhd_adi + prim2cons_ptr => prim2cons_srhd_adi cons2prim => cons2prim_srhd_adi fluxspeed => fluxspeed_srhd_adi get_maximum_speeds => get_maximum_speeds_srhd_adi @@ -869,7 +869,7 @@ module equations ! set pointers to subroutines ! - prim2cons => prim2cons_srmhd_adi + prim2cons_ptr => prim2cons_srmhd_adi cons2prim => cons2prim_srmhd_adi fluxspeed => fluxspeed_srmhd_adi get_maximum_speeds => get_maximum_speeds_srmhd_adi @@ -1166,9 +1166,7 @@ module equations ! status = 0 -! release the procedure pointers -! - nullify(prim2cons) + nullify(prim2cons_ptr) nullify(cons2prim) nullify(fluxspeed) @@ -1483,6 +1481,7 @@ module equations n = n + 1 +#ifdef DEBUG #if NDIMS == 3 sfmt = '("n: ",i0,", [i,j,k]: ",3(1x,i0),", [x,y,z]: ",3es12.4)' write(msg,sfmt) n, i, j, k, x(i), y(j), z(k) @@ -1491,6 +1490,7 @@ module equations write(msg,sfmt) n, i, j, x(i), y(j) #endif /* NDIMS == 3 */ call print_message(msg) +#endif /* DEBUG */ idx(:,n) = [ i, j, k ] @@ -1560,7 +1560,7 @@ module equations end do #endif /* NDIMS == 3 */ - call prim2cons(q(:,1:nc), u(:,1:nc), .true.) + call prim2cons_stencil(q(:,1:nc), u(:,1:nc), .true.) do n = 1, nc i = idx(1,n) @@ -1590,6 +1590,79 @@ module equations !! !=============================================================================== ! +!=============================================================================== +! +! subroutine PRIM2CONS_CELL: +! ------------------------- +! +! Subroutine converts primitive to conservative variables. +! +! Arguments: +! +! q - primitive variables (input); +! u - conservative variables (output); +! +!=============================================================================== +! + subroutine prim2cons_cell(q, u) + + implicit none + + real(kind=8), dimension(:), intent(in) :: q + real(kind=8), dimension(:), intent(out) :: u + +!------------------------------------------------------------------------------- +! + call prim2cons_ptr(q, u) + +!------------------------------------------------------------------------------- +! + end subroutine prim2cons_cell +! +!=============================================================================== +! +! subroutine PRIM2CONS_STENCIL: +! ---------------------------- +! +! Subroutine converts primitive to conservative variables along a stencil. +! +! Arguments: +! +! q - primitive variables (input); +! u - conservative variables (output); +! s - an optional flag indicating that passive scalars have +! to be calculated too; +! +!=============================================================================== +! + subroutine prim2cons_stencil(q, u, s) + + implicit none + + real(kind=8), dimension(:,:), intent(in) :: q + real(kind=8), dimension(:,:), intent(out) :: u + logical , optional , intent(in) :: s + + integer :: i, p + +!------------------------------------------------------------------------------- +! + do i = 1, size(q,2) + call prim2cons_ptr(q(:,i), u(:,i)) + end do + + if (ns > 0 .and. present(s)) then + if (s) then + do p = isl, isu + u(p,:) = q(p,:) * u(idn,:) + end do + end if + end if + +!------------------------------------------------------------------------------- +! + end subroutine prim2cons_stencil +! !******************************************************************************* ! ! ISOTHERMAL HYDRODYNAMIC EQUATIONS @@ -1601,56 +1674,28 @@ module equations ! subroutine PRIM2CONS_HD_ISO: ! --------------------------- ! -! Subroutine converts primitive variables to their corresponding -! conservative representation. +! Subroutine converts primitive to conservative variables. ! ! Arguments: ! -! q - the input array of primitive variables; -! u - the output array of conservative variables; -! s - an optional flag indicating that passive scalars have -! to be calculated too; +! q - primitive variables (input); +! u - conservative variables (output); ! !=============================================================================== ! - subroutine prim2cons_hd_iso(q, u, s) + subroutine prim2cons_hd_iso(q, u) -! local variables are not implicit by default -! implicit none -! input/output arguments -! - real(kind=8), dimension(:,:), intent(in) :: q - real(kind=8), dimension(:,:), intent(out) :: u - logical , optional , intent(in) :: s - -! local variables -! - integer :: i, p + real(kind=8), dimension(:), intent(in) :: q + real(kind=8), dimension(:), intent(out) :: u ! !------------------------------------------------------------------------------- ! -! iterate over all positions -! - do i = 1, size(q,2) - - u(idn,i) = q(idn,i) - u(imx,i) = q(idn,i) * q(ivx,i) - u(imy,i) = q(idn,i) * q(ivy,i) - u(imz,i) = q(idn,i) * q(ivz,i) - - end do - -! update primitive passive scalars -! - if (ns > 0 .and. present(s)) then - if (s) then - do p = isl, isu - u(p,:) = q(p,:) * u(idn,:) - end do - end if - end if + u(idn) = q(idn) + u(imx) = q(idn) * q(ivx) + u(imy) = q(idn) * q(ivy) + u(imz) = q(idn) * q(ivz) !------------------------------------------------------------------------------- ! @@ -1850,61 +1895,33 @@ module equations ! subroutine PRIM2CONS_HD_ADI: ! --------------------------- ! -! Subroutine converts primitive variables to their corresponding -! conservative representation. +! Subroutine converts primitive to conservative variables. ! ! Arguments: ! -! q - the input array of primitive variables; -! u - the output array of conservative variables; -! s - an optional flag indicating that passive scalars have -! to be calculated too; +! q - primitive variables (input); +! u - conservative variables (output); ! !=============================================================================== ! - subroutine prim2cons_hd_adi(q, u, s) + subroutine prim2cons_hd_adi(q, u) -! local variables are not implicit by default -! implicit none -! input/output arguments -! - real(kind=8), dimension(:,:), intent(in) :: q - real(kind=8), dimension(:,:), intent(out) :: u - logical , optional , intent(in) :: s + real(kind=8), dimension(:), intent(in) :: q + real(kind=8), dimension(:), intent(out) :: u -! local variables -! - integer :: i, p real(kind=8) :: ek, ei ! !------------------------------------------------------------------------------- ! -! iterate over all positions -! - do i = 1, size(q,2) - - u(idn,i) = q(idn,i) - u(imx,i) = q(idn,i) * q(ivx,i) - u(imy,i) = q(idn,i) * q(ivy,i) - u(imz,i) = q(idn,i) * q(ivz,i) - ek = 0.5d+00 * (u(imx,i) * q(ivx,i) + u(imy,i) * q(ivy,i) & - + u(imz,i) * q(ivz,i)) - ei = gammam1i * q(ipr,i) - u(ien,i) = ei + ek - - end do - -! update primitive passive scalars -! - if (ns > 0 .and. present(s)) then - if (s) then - do p = isl, isu - u(p,:) = q(p,:) * u(idn,:) - end do - end if - end if + u(idn) = q(idn) + u(imx) = q(idn) * q(ivx) + u(imy) = q(idn) * q(ivy) + u(imz) = q(idn) * q(ivz) + ek = 5.0d-01 * (u(imx) * q(ivx) + u(imy) * q(ivy) + u(imz) * q(ivz)) + ei = gammam1i * q(ipr) + u(ien) = ei + ek !------------------------------------------------------------------------------- ! @@ -2125,60 +2142,32 @@ module equations ! subroutine PRIM2CONS_MHD_ISO: ! ---------------------------- ! -! Subroutine converts primitive variables to their corresponding -! conservative representation. +! Subroutine converts primitive to conservative variables. ! ! Arguments: ! -! q - the input array of primitive variables; -! u - the output array of conservative variables; -! s - an optional flag indicating that passive scalars have -! to be calculated too; +! q - primitive variables (input); +! u - conservative variables (output); ! !=============================================================================== ! - subroutine prim2cons_mhd_iso(q, u, s) + subroutine prim2cons_mhd_iso(q, u) -! local variables are not implicit by default -! implicit none -! input/output arguments -! - real(kind=8), dimension(:,:), intent(in) :: q - real(kind=8), dimension(:,:), intent(out) :: u - logical , optional , intent(in) :: s - -! local variables -! - integer :: i, p + real(kind=8), dimension(:), intent(in) :: q + real(kind=8), dimension(:), intent(out) :: u ! !------------------------------------------------------------------------------- ! -! iterate over all positions -! - do i = 1, size(q,2) - - u(idn,i) = q(idn,i) - u(imx,i) = q(idn,i) * q(ivx,i) - u(imy,i) = q(idn,i) * q(ivy,i) - u(imz,i) = q(idn,i) * q(ivz,i) - u(ibx,i) = q(ibx,i) - u(iby,i) = q(iby,i) - u(ibz,i) = q(ibz,i) - u(ibp,i) = q(ibp,i) - - end do - -! update primitive passive scalars -! - if (ns > 0 .and. present(s)) then - if (s) then - do p = isl, isu - u(p,:) = q(p,:) * u(idn,:) - end do - end if - end if + u(idn) = q(idn) + u(imx) = q(idn) * q(ivx) + u(imy) = q(idn) * q(ivy) + u(imz) = q(idn) * q(ivz) + u(ibx) = q(ibx) + u(iby) = q(iby) + u(ibz) = q(ibz) + u(ibp) = q(ibp) !------------------------------------------------------------------------------- ! @@ -2410,56 +2399,39 @@ module equations ! subroutine PRIM2CONS_MHD_ADI: ! ---------------------------- ! -! Subroutine converts primitive variables to their corresponding -! conservative representation. +! Subroutine converts primitive to conservative variables. ! ! Arguments: ! -! q - the input array of primitive variables; -! u - the output array of conservative variables; -! s - an optional flag indicating that passive scalars have -! to be calculated too; +! q - primitive variables (input); +! u - conservative variables (output); ! !=============================================================================== ! - subroutine prim2cons_mhd_adi(q, u, s) + subroutine prim2cons_mhd_adi(q, u) implicit none - real(kind=8), dimension(:,:), intent(in) :: q - real(kind=8), dimension(:,:), intent(out) :: u - logical , optional , intent(in) :: s + real(kind=8), dimension(:), intent(in) :: q + real(kind=8), dimension(:), intent(out) :: u - integer :: i, p real(kind=8) :: ei, ek, em, ep !------------------------------------------------------------------------------- ! - do i = 1, size(q,2) - - u(idn,i) = q(idn,i) - u(imx,i) = q(idn,i) * q(ivx,i) - u(imy,i) = q(idn,i) * q(ivy,i) - u(imz,i) = q(idn,i) * q(ivz,i) - u(ibx,i) = q(ibx,i) - u(iby,i) = q(iby,i) - u(ibz,i) = q(ibz,i) - u(ibp,i) = q(ibp,i) - ei = gammam1i * q(ipr,i) - ek = sum(u(imx:imz,i) * q(ivx:ivz,i)) - em = sum(q(ibx:ibz,i) * q(ibx:ibz,i)) - ep = q(ibp,i) * q(ibp,i) - u(ien,i) = ei + 5.0d-01 * (ek + em + ep) - - end do - - if (ns > 0 .and. present(s)) then - if (s) then - do p = isl, isu - u(p,:) = q(p,:) * u(idn,:) - end do - end if - end if + u(idn) = q(idn) + u(imx) = q(idn) * q(ivx) + u(imy) = q(idn) * q(ivy) + u(imz) = q(idn) * q(ivz) + u(ibx) = q(ibx) + u(iby) = q(iby) + u(ibz) = q(ibz) + u(ibp) = q(ibp) + ei = gammam1i * q(ipr) + ek = sum(u(imx:imz) * q(ivx:ivz)) + em = sum(q(ibx:ibz) * q(ibx:ibz)) + ep = q(ibp) * q(ibp) + u(ien) = ei + 5.0d-01 * (ek + em + ep) !------------------------------------------------------------------------------- ! @@ -2712,67 +2684,36 @@ module equations ! subroutine PRIM2CONS_SRHD_ADI: ! ----------------------------- ! -! Subroutine converts primitive variables to their corresponding -! conservative representation. +! Subroutine converts primitive to conservative variables. ! ! Arguments: ! -! q - the input array of primitive variables; -! u - the output array of conservative variables; -! s - an optional flag indicating that passive scalars have -! to be calculated too; +! q - primitive variables (input); +! u - conservative variables (output); ! !=============================================================================== ! - subroutine prim2cons_srhd_adi(q, u, s) + subroutine prim2cons_srhd_adi(q, u) -! local variables are not implicit by default -! implicit none -! input/output arguments -! - real(kind=8), dimension(:,:), intent(in) :: q - real(kind=8), dimension(:,:), intent(out) :: u - logical , optional , intent(in) :: s + real(kind=8), dimension(:), intent(in) :: q + real(kind=8), dimension(:), intent(out) :: u -! local variables -! - integer :: i, p real(kind=8) :: vv, vm, vs, ww -! + !------------------------------------------------------------------------------- ! -! iterate over all positions -! - do i = 1, size(q,2) + vv = sum(q(ivx:ivz) * q(ivx:ivz)) + vm = 1.0d+00 - vv + vs = sqrt(vm) + ww = (q(idn) + q(ipr) / gammaxi) / vm -! calculate the square of velocity, the Lorentz factor and specific enthalpy -! - vv = sum(q(ivx:ivz,i) * q(ivx:ivz,i)) - vm = 1.0d+00 - vv - vs = sqrt(vm) - ww = (q(idn,i) + q(ipr,i) / gammaxi) / vm - -! calculate conservative variables -! - u(idn,i) = q(idn,i) / vs - u(imx,i) = ww * q(ivx,i) - u(imy,i) = ww * q(ivy,i) - u(imz,i) = ww * q(ivz,i) - u(ien,i) = ww - q(ipr,i) - u(idn,i) - - end do - -! update primitive passive scalars -! - if (ns > 0 .and. present(s)) then - if (s) then - do p = isl, isu - u(p,:) = q(p,:) * u(idn,:) - end do - end if - end if + u(idn) = q(idn) / vs + u(imx) = ww * q(ivx) + u(imy) = ww * q(ivy) + u(imz) = ww * q(ivz) + u(ien) = ww - q(ipr) - u(idn) !------------------------------------------------------------------------------- ! @@ -3710,76 +3651,43 @@ module equations ! subroutine PRIM2CONS_SRMHD_ADI: ! ------------------------------ ! -! Subroutine converts primitive variables to their corresponding -! conservative representation. +! Subroutine converts primitive to conservative variables. ! ! Arguments: ! -! q - the input array of primitive variables; -! u - the output array of conservative variables; -! s - an optional flag indicating that passive scalars have -! to be calculated too; +! q - primitive variables (input); +! u - conservative variables (output); ! !=============================================================================== ! - subroutine prim2cons_srmhd_adi(q, u, s) + subroutine prim2cons_srmhd_adi(q, u) -! local variables are not implicit by default -! implicit none -! input/output arguments -! - real(kind=8), dimension(:,:), intent(in) :: q - real(kind=8), dimension(:,:), intent(out) :: u - logical , optional , intent(in) :: s + real(kind=8), dimension(:), intent(in) :: q + real(kind=8), dimension(:), intent(out) :: u -! local variables -! - integer :: i, p real(kind=8) :: vv, bb, vb, vm, vs, ww, wt -! + !------------------------------------------------------------------------------- ! -! iterate over all positions -! - do i = 1, size(q,2) + vv = sum(q(ivx:ivz) * q(ivx:ivz)) + bb = sum(q(ibx:ibz) * q(ibx:ibz)) + vb = sum(q(ivx:ivz) * q(ibx:ibz)) + vm = 1.0d+00 - vv + vs = sqrt(vm) + ww = (q(idn) + q(ipr) / gammaxi) / vm + wt = ww + bb -! calculate the square of velocity, the quare of magnetic field, the scalar -! product of velocity and magnetic field, the Lorentz factor, specific and -! total enthalpies -! - vv = sum(q(ivx:ivz,i) * q(ivx:ivz,i)) - bb = sum(q(ibx:ibz,i) * q(ibx:ibz,i)) - vb = sum(q(ivx:ivz,i) * q(ibx:ibz,i)) - vm = 1.0d+00 - vv - vs = sqrt(vm) - ww = (q(idn,i) + q(ipr,i) / gammaxi) / vm - wt = ww + bb - -! calculate conservative variables -! - u(idn,i) = q(idn,i) / vs - u(imx,i) = wt * q(ivx,i) - vb * q(ibx,i) - u(imy,i) = wt * q(ivy,i) - vb * q(iby,i) - u(imz,i) = wt * q(ivz,i) - vb * q(ibz,i) - u(ibx,i) = q(ibx,i) - u(iby,i) = q(iby,i) - u(ibz,i) = q(ibz,i) - u(ibp,i) = q(ibp,i) - u(ien,i) = wt - q(ipr,i) - u(idn,i) - 0.5d+00 * (vm * bb + vb * vb) - - end do - -! update primitive passive scalars -! - if (ns > 0 .and. present(s)) then - if (s) then - do p = isl, isu - u(p,:) = q(p,:) * u(idn,:) - end do - end if - end if + u(idn) = q(idn) / vs + u(imx) = wt * q(ivx) - vb * q(ibx) + u(imy) = wt * q(ivy) - vb * q(iby) + u(imz) = wt * q(ivz) - vb * q(ibz) + u(ibx) = q(ibx) + u(iby) = q(iby) + u(ibz) = q(ibz) + u(ibp) = q(ibp) + u(ien) = wt - q(ipr) - u(idn) - 5.0d-01 * (vm * bb + vb * vb) !------------------------------------------------------------------------------- ! diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 04a1398..d45cc7d 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -4105,6 +4105,15 @@ module evolution n = get_dblocks() +!$omp parallel do default(shared) private(pdata,pmeta) + do l = 1, n + pdata => data_blocks(l)%ptr + pmeta => pdata%meta + + if (pmeta%update) call update_shapes(pdata, tm, dtm) + end do +!$omp end parallel do + call boundary_variables(tm, dtm, status) #ifdef MPI @@ -4112,6 +4121,15 @@ module evolution #endif /* MPI */ if (status /= 0) go to 100 +!$omp parallel do default(shared) private(pdata,pmeta) + do l = 1, n + pdata => data_blocks(l)%ptr + pmeta => pdata%meta + + if (pmeta%boundary) call update_shapes(pdata, tm, dtm) + end do +!$omp end parallel do + !$omp parallel do default(shared) private(pdata,pmeta,s) do l = 1, n pdata => data_blocks(l)%ptr @@ -4138,6 +4156,8 @@ module evolution !$omp critical if (s /= 0) status = 1 !$omp end critical + + call update_shapes(pdata, tm, dtm) end if end do @@ -4154,15 +4174,6 @@ module evolution #endif /* MPI */ if (status /= 0) go to 100 -!$omp parallel do default(shared) private(pdata,pmeta) - do l = 1, n - pdata => data_blocks(l)%ptr - pmeta => pdata%meta - - if (pmeta%update .or. pmeta%boundary) call update_shapes(pdata, tm, dtm) - end do -!$omp end parallel do - 100 continue pmeta => list_meta diff --git a/sources/mesh.F90 b/sources/mesh.F90 index 42904db..cb184f3 100644 --- a/sources/mesh.F90 +++ b/sources/mesh.F90 @@ -1752,14 +1752,6 @@ module mesh end if end if ! pmeta belongs to the current process - -! redistribute blocks equally among all processors -! - call redistribute_blocks(status) - if (status /= 0) then - call print_message(loc, "Cannot redistribute blocks!") - go to 100 - end if #endif /* MPI */ end if ! leaf at current level selected for refinement