Merge branch 'master' into reconnection

This commit is contained in:
Grzegorz Kowal 2015-05-11 12:32:51 -03:00
commit 0778c00239
6 changed files with 156 additions and 156 deletions

View File

@ -5791,7 +5791,7 @@ module boundaries
use coordinates , only : im , jm , km
use coordinates , only : faces_dp
use equations , only : nv
use interpolations , only : limiter
use interpolations , only : limiter_prol
! local variables are not implicit by default
!
@ -5854,15 +5854,15 @@ module boundaries
!
dql = qn(p,i ,j,k) - qn(p,im1,j,k)
dqr = qn(p,ip1,j,k) - qn(p,i ,j,k)
dqx = limiter(0.25d+00, dql, dqr)
dqx = limiter_prol(0.25d+00, dql, dqr)
dql = qn(p,i,j ,k) - qn(p,i,jm1,k)
dqr = qn(p,i,jp1,k) - qn(p,i,j ,k)
dqy = limiter(0.25d+00, dql, dqr)
dqy = limiter_prol(0.25d+00, dql, dqr)
dql = qn(p,i,j,k ) - qn(p,i,j,km1)
dqr = qn(p,i,j,kp1) - qn(p,i,j,k )
dqz = limiter(0.25d+00, dql, dqr)
dqz = limiter_prol(0.25d+00, dql, dqr)
! calculate the derivative terms
!
@ -6078,7 +6078,7 @@ module boundaries
use coordinates , only : im , jm , km
use coordinates , only : edges_dp
use equations , only : nv
use interpolations , only : limiter
use interpolations , only : limiter_prol
! local variables are not implicit by default
!
@ -6154,16 +6154,16 @@ module boundaries
!
dql = qn(p,i ,j,k) - qn(p,im1,j,k)
dqr = qn(p,ip1,j,k) - qn(p,i ,j,k)
dqx = limiter(0.25d+00, dql, dqr)
dqx = limiter_prol(0.25d+00, dql, dqr)
dql = qn(p,i,j ,k) - qn(p,i,jm1,k)
dqr = qn(p,i,jp1,k) - qn(p,i,j ,k)
dqy = limiter(0.25d+00, dql, dqr)
dqy = limiter_prol(0.25d+00, dql, dqr)
#if NDIMS == 3
dql = qn(p,i,j,k ) - qn(p,i,j,km1)
dqr = qn(p,i,j,kp1) - qn(p,i,j,k )
dqz = limiter(0.25d+00, dql, dqr)
dqz = limiter_prol(0.25d+00, dql, dqr)
#endif /* NDIMS == 3 */
#if NDIMS == 2
@ -6334,7 +6334,7 @@ module boundaries
use coordinates , only : im , jm , km
use coordinates , only : corners_dp
use equations , only : nv
use interpolations , only : limiter
use interpolations , only : limiter_prol
! local variables are not implicit by default
!
@ -6415,16 +6415,16 @@ module boundaries
!
dql = qn(p,i ,j,k) - qn(p,im1,j,k)
dqr = qn(p,ip1,j,k) - qn(p,i ,j,k)
dqx = limiter(0.25d+00, dql, dqr)
dqx = limiter_prol(0.25d+00, dql, dqr)
dql = qn(p,i,j ,k) - qn(p,i,jm1,k)
dqr = qn(p,i,jp1,k) - qn(p,i,j ,k)
dqy = limiter(0.25d+00, dql, dqr)
dqy = limiter_prol(0.25d+00, dql, dqr)
#if NDIMS == 3
dql = qn(p,i,j,k ) - qn(p,i,j,km1)
dqr = qn(p,i,j,kp1) - qn(p,i,j,k )
dqz = limiter(0.25d+00, dql, dqr)
dqz = limiter_prol(0.25d+00, dql, dqr)
#endif /* NDIMS == 3 */
#if NDIMS == 2

View File

@ -287,7 +287,6 @@ module evolution
! include external procedures
!
use blocks , only : set_blocks_update
use boundaries , only : boundary_variables
use mesh , only : update_mesh
! include external variables
@ -334,10 +333,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
#ifdef DEBUG
! check variables for NaNs
!
@ -516,7 +511,6 @@ module evolution
! include external procedures
!
use boundaries , only : boundary_variables
use sources , only : update_sources
! include external variables
@ -589,10 +583,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
#ifdef PROFILE
! stop accounting time for one step update
!
@ -623,7 +613,6 @@ module evolution
! include external procedures
!
use boundaries , only : boundary_variables
use sources , only : update_sources
! include external variables
@ -691,10 +680,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
! update fluxes from the intermediate stage
!
call update_fluxes()
@ -740,10 +725,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
#ifdef PROFILE
! stop accounting time for one step update
!
@ -777,7 +758,6 @@ module evolution
! include external procedures
!
use boundaries , only : boundary_variables
use sources , only : update_sources
! include external variables
@ -901,10 +881,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
end do ! n = 1, stages - 1
!= the final step: U(n+1) = 1/m U(0) + (m-1)/m [1 + dt/(m-1) L] U(m-1)
@ -954,10 +930,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
#ifdef PROFILE
! stop accounting time for one step update
!
@ -989,7 +961,6 @@ module evolution
! include external procedures
!
use boundaries , only : boundary_variables
use sources , only : update_sources
! include external variables
@ -1072,10 +1043,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
!! 2nd substep of integration
!!
! prepare the fractional time step
@ -1118,10 +1085,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
!! 3rd substep of integration
!!
! prepare the fractional time step
@ -1173,10 +1136,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
#ifdef PROFILE
! stop accounting time for one step update
!
@ -1209,7 +1168,6 @@ module evolution
! include external procedures
!
use boundaries , only : boundary_variables
use sources , only : update_sources
! include external variables
@ -1292,10 +1250,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
!= 2nd step: U(2) = U(1) + 1/2 dt F[U(1)]
!
! update fluxes
@ -1333,10 +1287,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
!= 3rd step: U(3) = 2/3 U(n) + 1/3 U(2) + 1/6 dt F[U(2)]
!
! calculate the fractional time step
@ -1379,10 +1329,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
!= the final step: U(n+1) = U(3) + 1/2 dt F[U(3)]
!
! calculate fractional time step
@ -1433,10 +1379,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
#ifdef PROFILE
! stop accounting time for one step update
!
@ -1469,7 +1411,6 @@ module evolution
! include external procedures
!
use boundaries , only : boundary_variables
use sources , only : update_sources
! include external variables
@ -1560,10 +1501,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
!= 2nd step: U(2) = U(1) + b1 dt F[U(1)]
!
! update fluxes
@ -1601,10 +1538,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
!= 3rd step: U(3) = a31 U(n) + a33 U(2) + b3 dt F[U(2)]
!
! calculate the fractional time step
@ -1647,10 +1580,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
!= 4th step: U(4) = a41 U(n) + a44 U(3) + b4 dt F[U(3)]
!
! calculate the fractional time step
@ -1697,10 +1626,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
!= the final step: U(n+1) = a53 U(2) + a55 U(4) + b5 dt F[U(4)]
!
! calculate the fractional time step
@ -1748,10 +1673,6 @@ module evolution
!
call update_variables()
! update boundaries
!
call boundary_variables()
#ifdef PROFILE
! stop accounting time for one step update
!
@ -1961,6 +1882,7 @@ module evolution
! include external procedures
!
use boundaries , only : boundary_variables
use equations , only : update_primitive_variables
use shapes , only : update_shapes
@ -1986,32 +1908,24 @@ module evolution
call start_timer(imv)
#endif /* PROFILE */
! associate the pointer with the first block on the data block list
! update primitive variables and shapes if necessary
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! associate pmeta with the corresponding meta block
!
pmeta => pdata%meta
! convert conserved variables to primitive ones for the current block and
! update shapes if necessary
!
if (pmeta%update) then
call update_primitive_variables(pdata%u, pdata%q)
call update_shapes(pdata)
end if
! assign pointer to the next block
!
pdata => pdata%next
end do
! update boundaries
!
call boundary_variables()
#ifdef PROFILE
! stop accounting time for variable update
!

View File

@ -50,7 +50,9 @@ module interpolations
! pointers to the reconstruction and limiter procedures
!
procedure(reconstruct) , pointer, save :: reconstruct_states => null()
procedure(limiter_zero) , pointer, save :: limiter => null()
procedure(limiter_zero) , pointer, save :: limiter_tvd => null()
procedure(limiter_zero) , pointer, save :: limiter_prol => null()
procedure(limiter_zero) , pointer, save :: limiter_clip => null()
! module parameters
!
@ -73,7 +75,7 @@ module interpolations
! declare public subroutines
!
public :: initialize_interpolations, finalize_interpolations
public :: reconstruct, limiter
public :: reconstruct, limiter_prol
public :: fix_positivity
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
@ -111,11 +113,15 @@ module interpolations
! local variables
!
character(len=255) :: sreconstruction = "tvd"
character(len=255) :: slimiter = "mm"
character(len=255) :: tlimiter = "mm"
character(len=255) :: plimiter = "mm"
character(len=255) :: climiter = "mm"
character(len=255) :: positivity_fix = "off"
character(len=255) :: clip_extrema = "off"
character(len=255) :: name_rec = ""
character(len=255) :: name_lim = ""
character(len=255) :: name_tlim = ""
character(len=255) :: name_plim = ""
character(len=255) :: name_clim = ""
!
!-------------------------------------------------------------------------------
!
@ -135,9 +141,11 @@ module interpolations
! obtain the user defined interpolation methods and coefficients
!
call get_parameter_string ("reconstruction" , sreconstruction)
call get_parameter_string ("limiter" , slimiter )
call get_parameter_string ("limiter" , tlimiter )
call get_parameter_string ("fix_positivity" , positivity_fix )
call get_parameter_string ("clip_extrema" , clip_extrema )
call get_parameter_string ("extrema_limiter" , climiter )
call get_parameter_string ("prolongation_limiter", plimiter )
call get_parameter_integer("nghosts" , ng )
call get_parameter_real ("eps" , eps )
call get_parameter_real ("limo3_rad" , rad )
@ -207,27 +215,67 @@ module interpolations
end if
end select
! select the limiter
! select the TVD limiter
!
select case(trim(slimiter))
select case(trim(tlimiter))
case ("mm", "minmod")
name_lim = "minmod"
limiter => limiter_minmod
name_tlim = "minmod"
limiter_tvd => limiter_minmod
case ("mc", "monotonized_central")
name_lim = "monotonized central"
limiter => limiter_monotonized_central
name_tlim = "monotonized central"
limiter_tvd => limiter_monotonized_central
case ("sb", "superbee")
name_lim = "superbee"
limiter => limiter_superbee
name_tlim = "superbee"
limiter_tvd => limiter_superbee
case ("vl", "vanleer")
name_lim = "van Leer"
limiter => limiter_vanleer
name_tlim = "van Leer"
limiter_tvd => limiter_vanleer
case ("va", "vanalbada")
name_lim = "van Albada"
limiter => limiter_vanalbada
name_tlim = "van Albada"
limiter_tvd => limiter_vanalbada
case default
name_lim = "zero derivative"
limiter => limiter_zero
name_tlim = "zero derivative"
limiter_tvd => limiter_zero
end select
! select the prolongation limiter
!
select case(trim(plimiter))
case ("mm", "minmod")
name_plim = "minmod"
limiter_prol => limiter_minmod
case ("mc", "monotonized_central")
name_plim = "monotonized central"
limiter_prol => limiter_monotonized_central
case ("sb", "superbee")
name_plim = "superbee"
limiter_prol => limiter_superbee
case ("vl", "vanleer")
name_plim = "van Leer"
limiter_prol => limiter_vanleer
case default
name_plim = "zero derivative"
limiter_prol => limiter_zero
end select
! select the clipping limiter
!
select case(trim(climiter))
case ("mm", "minmod")
name_clim = "minmod"
limiter_clip => limiter_minmod
case ("mc", "monotonized_central")
name_clim = "monotonized central"
limiter_clip => limiter_monotonized_central
case ("sb", "superbee")
name_clim = "superbee"
limiter_clip => limiter_superbee
case ("vl", "vanleer")
name_clim = "van Leer"
limiter_clip => limiter_vanleer
case default
name_clim = "zero derivative"
limiter_clip => limiter_zero
end select
! check additional reconstruction limiting
@ -249,10 +297,14 @@ module interpolations
!
if (verbose) then
write (*,"(4x,a15,8x,'=',1x,a)") "reconstruction ", trim(name_rec)
write (*,"(4x,a15,8x,'=',1x,a)") "limiter ", trim(name_lim)
write (*,"(4x,a15,8x,'=',1x,a)") "fix positivity ", trim(positivity_fix)
write (*,"(4x,a15,8x,'=',1x,a)") "clip extrema ", trim(clip_extrema)
write (*,"(4x,a14, 9x,'=',1x,a)") "reconstruction" , trim(name_rec)
write (*,"(4x,a11,12x,'=',1x,a)") "TVD limiter" , trim(name_tlim)
write (*,"(4x,a20, 3x,'=',1x,a)") "prolongation limiter", trim(name_plim)
write (*,"(4x,a14, 9x,'=',1x,a)") "fix positivity" , trim(positivity_fix)
write (*,"(4x,a12,11x,'=',1x,a)") "clip extrema" , trim(clip_extrema)
if (clip) then
write (*,"(4x,a15,8x,'=',1x,a)") "extrema limiter", trim(name_clim)
end if
end if
@ -298,7 +350,9 @@ module interpolations
! release the procedure pointers
!
nullify(reconstruct_states)
nullify(limiter)
nullify(limiter_tvd)
nullify(limiter_prol)
nullify(limiter_clip)
#ifdef PROFILE
! stop accounting time for module initialization/finalization
@ -419,7 +473,7 @@ module interpolations
! obtain the TVD limited derivative
!
df = limiter(0.5d+00, dfl, dfr)
df = limiter_tvd(0.5d+00, dfl, dfr)
! update the left and right-side interpolation states
!
@ -2483,10 +2537,15 @@ module interpolations
!
real(kind=8), intent(in) :: x, a, b
real(kind=8) :: c
! local variables
!
real(kind=8) :: y
!
!-------------------------------------------------------------------------------
!
c = (sign(x, a) + sign(x, b)) * min(abs(a), abs(b), 2.5d-01 * abs(a + b))
y = x - eps
c = (sign(y, a) + sign(y, b)) * min(abs(a), abs(b), 2.5d-01 * abs(a + b))
!-------------------------------------------------------------------------------
!
@ -2517,11 +2576,16 @@ module interpolations
!
real(kind=8), intent(in) :: x, a, b
real(kind=8) :: c
! local variables
!
real(kind=8) :: y
!
!-------------------------------------------------------------------------------
!
c = 0.5d+00 * (sign(x, a) + sign(x, b)) &
* max(min(2.0d+00 * abs(a), abs(b)), min(abs(a), 2.0d+00 * abs(b)))
y = x - eps
c = (sign(y, a) + sign(y, b)) &
* max(min(abs(a), 0.5d+00 * abs(b)), min(0.5d+00 * abs(a), abs(b)))
!-------------------------------------------------------------------------------
!
@ -2557,7 +2621,7 @@ module interpolations
!
c = a * b
if (c > 0.0d+00) then
c = 2.0d+00 * x * c / (a + b)
c = 2.0d+00 * (x - eps) * c / (a + b)
else
c = 0.0d+00
end if
@ -2719,8 +2783,9 @@ module interpolations
! local variables
!
integer :: i, im1, ip1
integer :: i, im1, ip1, ip2
real(kind=8) :: fmn, fmx
real(kind=8) :: dfl, dfr, df
!
!------------------------------------------------------------------------------
!
@ -2748,10 +2813,19 @@ module interpolations
!
if (fl(i) < fmn .or. fl(i) > fmx) then
! calculate the left and right derivatives
!
dfl = f(i ) - f(im1)
dfr = f(ip1) - f(i )
! get the limited slope
!
df = limiter_clip(0.5d+00, dfl, dfr)
! calculate new states
!
fl(i ) = f(i )
fr(im1) = f(i )
fl(i ) = f(i ) + df
fr(im1) = f(i ) - df
end if
@ -2759,10 +2833,23 @@ module interpolations
!
if (fr(i) < fmn .or. fr(i) > fmx) then
! calculate the missing index
!
ip2 = min(n, i + 2)
! calculate the left and right derivatives
!
dfl = f(ip1) - f(i )
dfr = f(ip2) - f(ip1)
! get the limited slope
!
df = limiter_clip(0.5d+00, dfl, dfr)
! calculate new states
!
fl(ip1) = f(ip1)
fr(i ) = f(ip1)
fl(ip1) = f(ip1) + df
fr(i ) = f(ip1) - df
end if

View File

@ -1055,7 +1055,7 @@ module mesh
use coordinates , only : ng, nh, in, jn, kn, im, jm, km
use coordinates , only : ib, ie, jb, je, kb, ke
use equations , only : nv
use interpolations , only : limiter
use interpolations , only : limiter_prol
! local variables are not implicit by default
!
@ -1141,16 +1141,16 @@ module mesh
dul = pdata%u(p,i ,j,k) - pdata%u(p,i-1,j,k)
dur = pdata%u(p,i+1,j,k) - pdata%u(p,i ,j,k)
dux = limiter(0.25d+00, dul, dur)
dux = limiter_prol(0.25d+00, dul, dur)
dul = pdata%u(p,i,j ,k) - pdata%u(p,i,j-1,k)
dur = pdata%u(p,i,j+1,k) - pdata%u(p,i,j ,k)
duy = limiter(0.25d+00, dul, dur)
duy = limiter_prol(0.25d+00, dul, dur)
#if NDIMS == 3
dul = pdata%u(p,i,j,k ) - pdata%u(p,i,j,k-1)
dur = pdata%u(p,i,j,k+1) - pdata%u(p,i,j,k )
duz = limiter(0.25d+00, dul, dur)
duz = limiter_prol(0.25d+00, dul, dur)
#endif /* NDIMS == 3 */
#if NDIMS == 2

View File

@ -114,7 +114,7 @@ module schemes
character(len=64) :: solver = "HLL"
character(len=64) :: statev = "primitive"
character(len=255) :: name_sol = ""
character(len=255) :: name_sts = ""
character(len=255) :: name_sts = "primitive"
!
!-------------------------------------------------------------------------------
!

View File

@ -601,7 +601,7 @@ module shapes
! local arrays
!
real(kind=8), dimension(nv,im) :: q, u
real(kind=8), dimension(nv) :: qj
real(kind=8), dimension(nv) :: qj, uj
real(kind=8), dimension(im) :: x
real(kind=8), dimension(jm) :: y
real(kind=8), dimension(km) :: z
@ -650,6 +650,7 @@ module shapes
qj(ibz) = bjet
qj(ibp) = 0.0d+00
end if ! ibx > 0
call prim2cons(1, qj(1:nv), uj(1:nv))
! prepare block coordinates
!
@ -674,6 +675,7 @@ module shapes
! copy the primitive variable vector
!
q(1:nv,1:im) = pdata%q(1:nv,1:im,j,k)
u(1:nv,1:im) = pdata%u(1:nv,1:im,j,k)
! calculate radius
!
@ -683,22 +685,19 @@ module shapes
do i = 1, im
if (x(i) <= max(dx, ljet)) then
q(1:nv,i) = qj(1:nv)
u(1:nv,i) = uj(1:nv)
end if
end do ! i = 1, im
end if ! R < Rjet
! convert the primitive variables to conservative ones
! copy the primitive variables to the current block
!
call prim2cons(im, q(1:nv,1:im), u(1:nv,1:im))
pdata%q(1:nv,1:im,j,k) = q(1:nv,1:im)
! copy the conserved variables to the current block
!
pdata%u(1:nv,1:im,j,k) = u(1:nv,1:im)
! copy the primitive variables to the current block
!
pdata%q(1:nv,1:im,j,k) = q(1:nv,1:im)
end do ! j = 1, jm
end do ! k = 1, km