diff --git a/src/equations.F90 b/src/equations.F90 index 3d78648..5a8c821 100644 --- a/src/equations.F90 +++ b/src/equations.F90 @@ -158,9 +158,9 @@ module equations integer , save :: nrmax = 100 integer , save :: nrext = 2 -! flags for reconstruction corrections +! flags for corrections ! - logical , save :: positivity = .false. + logical , save :: fix_unphysical_cells = .false. ! by default everything is private ! @@ -174,6 +174,7 @@ module equations public :: maxspeed, reset_maxspeed, get_maxspeed public :: eigensystem_roe public :: update_primitive_variables + public :: fix_unphysical_cells, correct_unphysical_states public :: gamma public :: csnd, csnd2 public :: cmax, cmax2 @@ -231,7 +232,7 @@ module equations character(len=255) :: name_eqsys = "" character(len=255) :: name_eos = "" character(len=255) :: name_c2p = "" - character(len=255) :: positivity_fix = "off" + character(len=255) :: unphysical_fix = "off" ! !------------------------------------------------------------------------------- ! @@ -854,17 +855,17 @@ module equations ! msfac = 1.0d+00 / (gamma * msmax**2) -! get the positivity fix flag +! get the state correction flags ! - call get_parameter_string("fix_positivity", positivity_fix ) + call get_parameter_string("fix_unphysical_cells", unphysical_fix ) -! check additional reconstruction limiting +! check if the correction of unphysical cells is on ! - select case(trim(positivity_fix)) + select case(trim(unphysical_fix)) case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") - positivity = .true. + fix_unphysical_cells = .true. case default - positivity = .false. + fix_unphysical_cells = .false. end select ! print information about the equation module @@ -876,6 +877,8 @@ module equations if (relativistic) then write (*,"(4x,a,1x,a)" ) "variable conversion =", trim(name_c2p) end if + write (*,"(4x,a20, 3x,'=',1x,a)") "fix unphysical cells" & + , trim(unphysical_fix) end if @@ -1043,7 +1046,6 @@ module equations ! temporary variables ! integer :: i, j, k - real(kind=8) :: pmin ! !------------------------------------------------------------------------------- ! @@ -1059,37 +1061,6 @@ module equations end do ! j = jb, je end do ! k = kb, ke -! fix negative pressure is desired -! - if (positivity .and. ipr > 0) then - -! iterate over block interior -! - do k = kb, ke - do j = jb, je - do i = ib, ie - -! fix the cells where pressure is negative -! - if (qq(ipr,i,j,k) <= 0.0d+00) then - -! calculate pressure from the sonic Mach number limit and local velocity -! - pmin = msfac * qq(idn,i,j,k) * sum(qq(ivx:ivz,i,j,k)**2) - -! update total energy and pressure -! - uu(ien,i,j,k) = uu(ien,i,j,k) + gammam1i * (pmin - qq(ipr,i,j,k)) - qq(ipr,i,j,k) = pmin - - end if - - end do ! i = ib, ie - end do ! j = jb, je - end do ! k = kb, ke - - end if - ! fill out the borders ! do i = ibl, 1, -1 @@ -1118,6 +1089,173 @@ module equations end subroutine update_primitive_variables ! !=============================================================================== +! +! subroutine CORRECT_UNPHYSICAL_STATES: +! ------------------------------------ +! +! Subroutine seeks for unphysical states (cells with negative density or +! pressure) and try to fix them by averaging their values from physical +! neighbours. +! +! Arguments: +! +! id - the block id where the states are being checked; +! qq - the output array of primitive variables; +! uu - the input array of conservative variables; +! +!=============================================================================== +! + subroutine correct_unphysical_states(id, qq, uu) + +! include external procedures and variables +! + use coordinates, only : im, jm, km + use error , only : print_warning, print_error + +! local variables are not implicit by default +! + implicit none + +! input/output arguments +! + integer(kind=4) , intent(in) :: id + real(kind=8), dimension(nv,im,jm,km), intent(inout) :: qq + real(kind=8), dimension(nv,im,jm,km), intent(inout) :: uu + +! temporary variables +! + character(len=255) :: msg + integer :: n, p, nc, np + integer :: i, il, iu + integer :: j, jl, ju + integer :: k, kl, ku + +! temporary arrays +! + logical, dimension(im,jm,km) :: physical + +! allocatable arrays +! + integer , dimension(:,:), allocatable :: idx + real(kind=8), dimension(:,:), allocatable :: q, u + +! local parameters +! + character(len=*), parameter :: loc = 'EQUATIONS::correct_unphysical_states()' +! +!------------------------------------------------------------------------------- +! +! search for negative density or pressure +! + physical(:,:,:) = qq(idn,:,:,:) > 0.0d+00 + if (ipr > 0) then + physical(:,:,:) = physical(:,:,:) .and. qq(ipr,:,:,:) > 0.0d+00 + end if + +! apply averaging for unphysical states +! + if (.not. all(physical)) then + +! count unphysical cells +! + nc = count(.not. physical) + +! inform about the encountered unphysical states +! + write(msg,'(i4,1x,a,1x,i6)') nc, "unphysical states in block", id + call print_warning(loc, trim(msg)) + +! allocate temporary vectors for unphysical states +! + allocate(q(nv,nc), u(nv,nc), idx(3,nc)) + +! iterate over block cells +! + n = 0 + do k = 1, km + do j = 1, jm + do i = 1, im + +! fix unphysical states +! + if (.not. physical(i,j,k)) then + + n = n + 1 + + idx(:,n) = (/ i, j, k /) + +! increase the region until we find at least two physical cells, but no more +! than 2 cells away +! + np = 0 + p = 1 + do while (n < 2 .and. p < 3) + il = max( 1, i - p) + iu = min(im, i + p) + jl = max( 1, j - p) + ju = min(jm, j + p) + kl = max( 1, k - p) + ku = min(km, k + p) + + np = count(physical(il:iu,jl:ju,kl:ku)) + + p = p + 1 + end do + +! average primitive variables +! + if (np >= 2) then + + do p = 1, nv + q(p,n) = sum(qq(p,il:iu,jl:ju,kl:ku), & + physical(il:iu,jl:ju,kl:ku)) & + / np + end do + + else + +! print error, since no physical cells found for averaging +! + write(msg,'(a,a)') & + "Not sufficient number of physical cells found!" & + , "Cannot correct the unphysical cell." + call print_error(loc, trim(msg)) + stop + + end if ! not sufficient number of physical cells for averaging + + end if ! not physical + + end do ! i = 1, im + end do ! j = 1, jm + end do ! k = 1, km + +! convert the vector of primitive variables to conservative ones +! + call prim2cons(nc, q(1:nc,1:nc), u(1:nv,1:nc)) + +! update block variables +! + do n = 1, nc + i = idx(1,n) + j = idx(2,n) + k = idx(3,n) + + qq(1:nv,i,j,k) = q(1:nv,n) + uu(1:nv,i,j,k) = u(1:nv,n) + end do + +! deallocate temporary vectors +! + deallocate(q, u, idx) + + end if ! there are unphysical cells + +!------------------------------------------------------------------------------- +! + end subroutine correct_unphysical_states +! +!=============================================================================== !! !!*** PRIVATE SUBROUTINES **************************************************** !! @@ -3337,6 +3475,10 @@ module equations ! subroutine cons2prim_srhd_adi(n, u, q) +! include external procedures +! + use error, only : print_warning + ! local variables are not implicit by default ! implicit none @@ -3353,6 +3495,10 @@ module equations integer :: i real(kind=8) :: mm, bb, mb, en, dn real(kind=8) :: w , vv, vm, vs + +! local parameters +! + character(len=*), parameter :: loc = 'EQUATIONS::cons2prim_srhd_adi()' ! !------------------------------------------------------------------------------- ! @@ -3393,15 +3539,17 @@ module equations q(ivz,i) = u(imz,i) / w q(ipr,i) = w - en - else ! unphysical state + else ! cannot find physical solution - write(*,"(a,1x,a)" ) "ERROR in" & - , "EQUATIONS::cons2prim_srhd_adi()" - write(*,"(a,5(1x,1e24.16e3))") "Unphysical state for U = ", u(1:nv,i) - write(*,"(a,3(1x,1e24.16e3))") " D, |m|², E = ", dn, mm, en - stop + call print_warning(loc, "Conversion to physical primitive state failed!") + write(*,"(a,5(1x,1e24.16e3))") "U = ", u(1:nv,i) + write(*,"(a,3(1x,1e24.16e3))") "D, |m|², E = ", dn, mm, en - end if ! unphysical state +! set pressure to zero so we can hopefully fix it later +! + q(ipr,i) = 0.0d+00 + + end if end do ! i = 1, n @@ -4436,6 +4584,10 @@ module equations ! subroutine cons2prim_srmhd_adi(n, u, q) +! include external procedures +! + use error, only : print_warning + ! local variables are not implicit by default ! implicit none @@ -4452,6 +4604,10 @@ module equations integer :: i real(kind=8) :: mm, mb, bb, en, dn real(kind=8) :: w, wt, vv, vm, vs, vb, fc + +! local parameters +! + character(len=*), parameter :: loc = 'EQUATIONS::cons2prim_srmhd_adi()' ! !------------------------------------------------------------------------------- ! @@ -4506,25 +4662,30 @@ module equations ! if (q(ipr,i) <= 0.0d+00) then - write(*,*) - write(*,"(a,1x,a)" ) "WARNING in" & - , "EQUATIONS::cons2prim_srmhd_adi()" - write(*,"(a,9(1x,1e24.16e3))") "Negative pressure for U = ", u(1:nv,i) + call print_warning(loc, & + "Conversion to physical primitive state" // & + " resulted in negative pressure!") + write(*,"(a,9(1x,1e24.16e3))") "U = ", u(1:nv,i) write(*,"(a,6(1x,1e24.16e3))") " D, |m|², m.B, |B|², E, W = " & , dn, mm, mb, bb, en, w - write(*,"(a,1(1x,1e24.16e3))") "Pressure corrected to ", pmin - q(ipr,i) = pmin + +! set pressure to zero so we can hopefully fix it later +! + q(ipr,i) = 0.0d+00 end if ! p <= 0 else ! unphysical state - write(*,*) - write(*,"(a,1x,a)" ) "ERROR in" & - , "EQUATIONS::cons2prim_srmhd_adi()" - write(*,"(a,9(1x,1e24.16e3))") "Unphysical state for U = ", u(1:nv,i) - write(*,"(a,5(1x,1e24.16e3))") " D, |m|², m.B, |B|², E = ", dn, mm, mb & - , bb, en + call print_warning(loc, & + "Conversion to physical primitive state failed!") + write(*,"(a,9(1x,1e24.16e3))") "U = ", u(1:nv,i) + write(*,"(a,5(1x,1e24.16e3))") "D, |m|², m.B, |B|², E = ", dn, mm, mb & + , bb, en + +! set pressure to zero so we can hopefully fix it later +! + q(ipr,i) = 0.0d+00 end if ! unphysical state diff --git a/src/evolution.F90 b/src/evolution.F90 index 9767fb8..486663e 100644 --- a/src/evolution.F90 +++ b/src/evolution.F90 @@ -1990,6 +1990,7 @@ module evolution ! use boundaries , only : boundary_variables use equations , only : update_primitive_variables + use equations , only : fix_unphysical_cells, correct_unphysical_states use shapes , only : update_shapes ! include external variables @@ -2033,6 +2034,20 @@ module evolution ! call boundary_variables(tm, dtm) +! correct unphysical states if detected +! + if (fix_unphysical_cells) then + pdata => list_data + do while (associated(pdata)) + pmeta => pdata%meta + + if (pmeta%update) call correct_unphysical_states(pmeta%id & + , pdata%q, pdata%u) + + pdata => pdata%next + end do + end if + ! apply shapes in blocks which need it ! pdata => list_data diff --git a/src/makefile b/src/makefile index 971aaa1..779a9f1 100644 --- a/src/makefile +++ b/src/makefile @@ -127,7 +127,8 @@ driver.o : driver.F90 blocks.o coordinates.o equations.o evolution.o \ gravity.o integrals.o interpolations.o io.o mesh.o \ mpitools.o operators.o parameters.o problems.o random.o \ refinement.o schemes.o shapes.o sources.o user_problem.o -equations.o : equations.F90 algebra.o coordinates.o parameters.o timers.o +equations.o : equations.F90 algebra.o coordinates.o error.o parameters.o \ + timers.o error.o : error.F90 evolution.o : evolution.F90 blocks.o boundaries.o coordinates.o \ equations.o mesh.o mpitools.o parameters.o schemes.o \