gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into reconnection

Browse Source
This commit is contained in:
Grzegorz Kowal 2020-11-10 12:09:48 -03:00
commit 79c4348a4f
5 changed files with 45 additions and 80 deletions

2
python/amunpy.py

@ -157,7 +157,7 @@ class AmunXML:
and all(v in variables for v in ['dens','velo']):
variables.append('ekin')
if self.attributes['eqsys'] in [ 'mhd', 'srmhd' ] \
and 'emag' in variables:
and 'magn' in variables:
variables.append('emag')
if 'velo' in variables and 'magn' in variables:
variables.append('elex')

2
sources/boundaries.F90

@ -5070,7 +5070,7 @@ module boundaries
! prepare indices for the boundaries
!
if (side(1) == 1) then
if (side(2) == 1) then
jl = 1
ju = nh - 1
else

17
sources/driver.F90

@ -86,7 +86,6 @@ program amun
use timers , only : start_timer, stop_timer, set_timer, get_timer
use timers , only : get_timer_total, timer_enabled, timer_description
use timers , only : get_count, ntimers
use user_problem , only : initialize_user_problem, finalize_user_problem
! module variables are not implicit by default
!
@ -401,15 +400,7 @@ program amun
end if
go to 3100
end if
call initialize_user_problem(problem, master, status)
if (check_status(status /= 0)) then
if (master) then
write(error_unit,"('[AMUN::program]: ', a)") &
"Problem initializing USER_PROBLEM module!"
end if
go to 3000
end if
call initialize_problems(problem, status)
call initialize_problems(problem, master, status)
if (check_status(status /= 0)) then
if (master) then
write(error_unit,"('[AMUN::program]: ', a)") &
@ -824,12 +815,6 @@ program amun
write(error_unit,"('[AMUN::program]: ', a)") &
"Problem finalizing PROBLEMS module!"
end if
3000 continue
call finalize_user_problem(status)
if (check_status(status /= 0) .and. master) then
write(error_unit,"('[AMUN::program]: ', a)") &
"Problem finalizing USER_PROBLEM module!"
end if
3100 continue
call finalize_sources(status)
if (check_status(status /= 0) .and. master) then

58
sources/problems.F90

@ -94,54 +94,37 @@ module problems
! Arguments:
!
! problem - the problem name
! verbose - a logical flag turning the information printing;
! status - an integer flag for error return value;
!
!===============================================================================
!
subroutine initialize_problems(problem, status)
subroutine initialize_problems(problem, verbose, status)
! include external procedures and variables
!
use parameters , only : get_parameter
use user_problem, only : setup_problem_user
use user_problem, only : initialize_user_problem, setup_problem_user
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
character(len=64), intent(in) :: problem
logical , intent(in) :: verbose
integer , intent(out) :: status
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! set timer descriptions
!
call set_timer('problems:: initialize', imi)
call set_timer('problems:: update' , imu)
! start accounting time for module initialization/finalization
!
call start_timer(imi)
#endif /* PROFILE */
! reset the status flag
!
status = 0
! set problem name
!
problem_name = problem
! associate the setup_problem pointer with the respective problem setup
! subroutine
!
select case(trim(problem))
! general test problems
!
case("riemann")
setup_problem => setup_problem_riemann
@ -169,11 +152,11 @@ module problems
case default
setup_problem => setup_problem_user
call initialize_user_problem(problem, verbose, status)
end select
#ifdef PROFILE
! stop accounting time for module initialization/finalization
!
call stop_timer(imi)
#endif /* PROFILE */
@ -196,33 +179,25 @@ module problems
!
subroutine finalize_problems(status)
! local variables are not implicit by default
!
use user_problem, only : finalize_user_problem
implicit none
! subroutine arguments
!
integer, intent(out) :: status
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for module initialization/finalization
!
call start_timer(imi)
#endif /* PROFILE */
! reset the status flag
!
status = 0
! nullify procedure pointers
!
nullify(setup_problem)
call finalize_user_problem(status)
#ifdef PROFILE
! stop accounting time for module initialization/finalization
!
call stop_timer(imi)
#endif /* PROFILE */
@ -2140,6 +2115,8 @@ module problems
real(kind=8), save :: eta = 1.00d-08
real(kind=8), save :: zeta = 0.00d+00
real(kind=8), save :: dens = 1.00d+00
real(kind=8), save :: pres = 5.00d-02
real(kind=8), save :: pmag = 1.00d+00
real(kind=8), save :: bamp = 1.00d+00
real(kind=8), save :: bnor = 0.00d+00
real(kind=8), save :: bgui = 0.00d+00
@ -2155,7 +2132,7 @@ module problems
! local variables
!
integer :: i, j, k = 1
real(kind=8) :: dx, dy, t, btot2 = 0.0d+00
real(kind=8) :: dx, dy, t
! local arrays
!
@ -2188,7 +2165,8 @@ module problems
ss = abs(bamp) / (sqrt(dens) * max(1.0d-16, eta))
delta = ss**(-1.0d+00/3.0d+00)
if (abs(bper) > 0.0d+00) vper = 0.0d+00
btot2 = bamp**2 + bnor**2 + bgui**2
pres = 0.5d+00 * beta
pmag = 0.5d+00 * (bamp**2 + bnor**2 + bgui**2)
first = .false.
end if
@ -2222,7 +2200,7 @@ module problems
q(ivx,:) = 0.0d+00
q(ivy,:) = 0.0d+00
q(ivz,:) = 0.0d+00
if (ipr > 0) q(ipr,:) = 0.5d+00 * beta
if (ipr > 0) q(ipr,:) = pres
! if magnetic field is present, set its antiparallel configuration
!
@ -2243,9 +2221,9 @@ module problems
! correct the gas pressure due to the variance of magnetic pressure
!
if (ipr > 0) then
q(ipr,:) = q(ipr,:) + 0.5d+00 * (btot2 - sum(q(ibx:ibz,:)**2,1))
q(ipr,:) = pres + pmag - 0.5d+00 * sum(q(ibx:ibz,:)**2,1)
else
q(idn,:) = 0.5d+00 * (beta + btot2 - sum(q(ibx:ibz,:)**2,1)) / csnd2
q(idn,:) = (pres + pmag - 0.5d+00 * sum(q(ibx:ibz,:)**2,1)) / csnd2
end if
end if

46
sources/schemes.F90

@ -1091,7 +1091,7 @@ module schemes
! calculate the fluxes for the intermediate state
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if ! sl < 0 < sr
@ -1191,7 +1191,7 @@ module schemes
! calculate the pressure of the intermediate state
!
pr = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn
pr = wl(idn) * wr(imx) / dn - wr(idn) * wl(imx) / dn
! separate intermediate states depending on the sign of the advection speed
!
@ -1346,7 +1346,7 @@ module schemes
! the total pressure, constant across the contact discontinuity
!
pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2
pt = wl(idn) * wr(imx) / dn - wr(idn) * wl(imx) / dn + b2
! separate the cases when Bₓ = 0 and Bₓ ≠ 0
!
@ -1415,7 +1415,7 @@ module schemes
! when Sₘ = 0 all variables are continuous, therefore the flux reduces
! to the HLL one
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if ! sm = 0
@ -1472,7 +1472,7 @@ module schemes
! when Sₘ = 0 all variables are continuous, therefore the flux reduces
! to the HLL one
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if ! sm = 0
@ -2056,7 +2056,7 @@ module schemes
! local variables
!
integer :: i
real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm
real(kind=8) :: sl, sr, sm, sml, smr, srml, slmm, srmm, sdif
real(kind=8) :: bx, b2, dn, dvl, dvr, ca
! local arrays to store the states
@ -2203,13 +2203,13 @@ module schemes
! the flux corresponding to the middle intermediate state
!
f(:,i) = (sml * wcr(:) - smr * wcl(:)) / (smr - sml)
f(:,i) = 0.5d+00 * (sml * wcr(:) / ca - smr * wcl(:) / ca)
else ! Bx = 0 -> Sₘ = 0
! no Alfvén wave and Sₘ = 0, so revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if
@ -2252,13 +2252,14 @@ module schemes
! the flux corresponding to the middle intermediate state
!
f(:,i) = (sml * wr(:) - sr * wcl(:)) / (sr - sml)
sdif = sr - sml
f(:,i) = sml * wr(:) / sdif - sr * wcl(:) / sdif
else ! Bx = 0 -> Sₘ = 0
! no Alfvén wave and Sₘ = 0, so revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if
@ -2297,13 +2298,14 @@ module schemes
! the flux corresponding to the middle intermediate state
!
f(:,i) = (sl * wcr(:) - smr * wl(:)) / (smr - sl)
sdif = smr - sl
f(:,i) = sl * wcr(:) / sdif - smr * wl(:) / sdif
else ! Bx = 0 -> Sₘ = 0
! no Alfvén wave and Sₘ = 0, so revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if
@ -2313,7 +2315,7 @@ module schemes
! both outer states are degenerate so revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if ! sl* < sl & sr* > sr
@ -2432,7 +2434,7 @@ module schemes
! the total pressure, constant across the contact discontinuity and Alfvén waves
!
pt = (wl(idn) * wr(imx) - wr(idn) * wl(imx)) / dn + b2
pt = wl(idn) * wr(imx) / dn - wr(idn) * wl(imx) / dn + b2
! left and right Alfvén speeds
!
@ -2560,7 +2562,7 @@ module schemes
! prepare constant primitive variables of the intermediate states
!
dv = car * dnr + cal * dnl
dv = abs(bx) * (sqrt(dnr) + sqrt(dnl))
vy = (wcr(imy) - wcl(imy)) / dv
vz = (wcr(imz) - wcl(imz)) / dv
dv = car + cal
@ -2612,7 +2614,7 @@ module schemes
! in the case when Sₘ = 0 and Bₓ² > 0, all variables are continuous, therefore
! the flux can be averaged from the Alfvén waves using a simple HLL formula;
!
f(:,i) = (sml * wcr(:) - smr * wcl(:)) / (smr - sml)
f(:,i) = sml * wcr(:) / dv - smr * wcl(:) / dv
end if ! sm = 0
@ -2674,7 +2676,7 @@ module schemes
! since Bx = 0 and Sm = 0, then revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if
@ -2765,7 +2767,7 @@ module schemes
! calculate the average flux over the right inner intermediate state
!
f(:,i) = (sm * wr(:) - sr * wcr(:)) / srmm
f(:,i) = sm * wr(:) / srmm - sr * wcr(:) / srmm
end if ! sm < 0
@ -2773,7 +2775,7 @@ module schemes
! no Alfvén wave, so revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if
@ -2858,7 +2860,7 @@ module schemes
! calculate the average flux over the left inner intermediate state
!
f(:,i) = (sm * wl(:) - sl * wcl(:)) / slmm
f(:,i) = sm * wl(:) / slmm - sl * wcl(:) / slmm
end if ! sm > 0
@ -2866,7 +2868,7 @@ module schemes
! no Alfvén wave, so revert to the HLL flux
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if
@ -2876,7 +2878,7 @@ module schemes
! so far we revert to HLL flux in the case of degeneracies
!
f(:,i) = (sl * wr(:) - sr * wl(:)) / srml
f(:,i) = sl * wr(:) / srml - sr * wl(:) / srml
end if ! sl* < sl & sr* > sr