From e75089ed19ae9c5c93d0d8288b360d6639bdd026 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 15 Aug 2020 01:09:52 -0300 Subject: [PATCH 01/12] INTERPOLATIONS: Check correctness of ngp just after reading it. Signed-off-by: Grzegorz Kowal --- sources/interpolations.F90 | 72 ++++++++++++++++---------------------- 1 file changed, 31 insertions(+), 41 deletions(-) diff --git a/sources/interpolations.F90 b/sources/interpolations.F90 index f329792..0db9c30 100644 --- a/sources/interpolations.F90 +++ b/sources/interpolations.F90 @@ -231,6 +231,17 @@ module interpolations ! kappa = min(kappa, (1.0d+00 - cfl) / cfl) +! check ngp +! + if (mod(ngp,2) == 0 .or. ngp < 3) then + if (verbose) then + write(error_unit,"('[', a, ']: ', a)") trim(loc), & + "The parameter ngp has to be an odd integer >= 3. "// & + "Resetting to default value of ngp = 5." + end if + ngp = 5 + end if + ! calculate mgp ! mgp = (ngp - 1) / 2 @@ -349,22 +360,12 @@ module interpolations reconstruct_states => reconstruct_gp order = ngp -! check the parameters -! - if (mod(ngp,2) == 0) then - if (verbose) then - write(*,*) - write(*,"(1x,a)") "ERROR!" - write(*,"(1x,a)") "The parameter ngp has to be an odd integer >= 3." - end if - status = 1 - else - ng = 2 - do while(ng < (ngp + 1) / 2) - ng = ng + 2 - end do - nghosts = max(nghosts, ng) - end if + ng = 2 + do while(ng < (ngp + 1) / 2) + ng = ng + 2 + end do + nghosts = max(nghosts, ng) + case ("mgp", "MGP") write(stmp, '(f16.1)') sgp write(name_rec, & @@ -379,35 +380,24 @@ module interpolations ! if (status == 0) call prepare_mgp(status) - interfaces => interfaces_mgp + interfaces => interfaces_mgp + + ng = 2 + do while(ng < (ngp + 1) / 2) + ng = ng + 2 + end do + nghosts = max(nghosts, ng) -! check the parameters -! - if (mod(ngp,2) == 0) then - if (verbose) then - write(*,*) - write(*,"(1x,a)") "ERROR!" - write(*,"(1x,a)") "The parameter ngp has to be an odd integer >= 3." - end if - status = 1 - else - ng = 2 - do while(ng < (ngp + 1) / 2) - ng = ng + 2 - end do - nghosts = max(nghosts, ng) - end if case default if (verbose) then - write(*,*) - write(*,"(1x,a)") "ERROR!" - write(*,"(1x,a)") "The selected reconstruction method is not " // & - "implemented: " // trim(sreconstruction) - write(*,"(1x,a)") "Available methods: 'tvd' 'limo3', 'ppm'," // & - " 'weno3', 'weno5z', 'weno5yc', 'weno5ns'," // & - " 'crweno5z', 'crweno5yc', 'crweno5ns','mp5', " // & - " 'mp7', 'crmp5', 'crmp5ld', 'crmp7', 'gp', 'mgp'." + write(error_unit,"('[', a, ']: ', a)") trim(loc), & + "The selected reconstruction method is not " // & + "implemented: " // trim(sreconstruction) // "." // & + "Available methods: 'tvd' 'limo3', 'ppm'," // & + " 'weno3', 'weno5z', 'weno5yc', 'weno5ns'," // & + " 'crweno5z', 'crweno5yc', 'crweno5ns','mp5', " // & + " 'mp7', 'mp9', 'crmp5', 'crmp5ld', 'crmp7', 'gp', 'mgp'." end if status = 1 From 9e507b283cdcbc5199960ca9abceed3cbd213a38 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 15 Aug 2020 01:12:41 -0300 Subject: [PATCH 02/12] SOURCES: Remove unused variable for 2D case. Signed-off-by: Grzegorz Kowal --- sources/sources.F90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/sources/sources.F90 b/sources/sources.F90 index 9b8be9e..2109fe0 100644 --- a/sources/sources.F90 +++ b/sources/sources.F90 @@ -321,7 +321,10 @@ module sources ! use blocks , only : block_data use coordinates , only : nn => bcells - use coordinates , only : ax, ay, az, adx, ady, adz + use coordinates , only : ax, adx, ay, ady, adz +#if NDIMS == 3 + use coordinates , only : az +#endif /* NDIMS == 3 */ use equations , only : inx, iny, inz use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ien use equations , only : ibx, iby, ibz, ibp From 02b2d5d1bd231fb4d94c78dbc8e9b14776ed35d5 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 15 Aug 2020 01:14:04 -0300 Subject: [PATCH 03/12] CMAKE: Do not warn about dummy arguments. Signed-off-by: Grzegorz Kowal --- CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9212d40..f2a36c4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,8 +19,8 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) endif() if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU") - add_compile_options("$<$:-march=native;-pipe;-pedantic;-ftree-vectorize;-fno-unsafe-math-optimizations;-frounding-math;-fsignaling-nans;-finline-limit=10000;-fdiagnostics-color=always>") - add_compile_options("$<$:-Og;-pedantic;-W;-Wall>") + add_compile_options("$<$:-march=native;-pipe;-ftree-vectorize;-fno-unsafe-math-optimizations;-frounding-math;-fsignaling-nans;-finline-limit=10000;-fdiagnostics-color=always>") + add_compile_options("$<$:-Og;-pedantic;-W;-Wall;-Wno-unused-dummy-argument>") endif() if(CMAKE_Fortran_COMPILER_ID MATCHES "Intel") From d05e3062fc98f9719f1d7268dfa91ac4661c7934 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 15 Aug 2020 01:25:31 -0300 Subject: [PATCH 04/12] FORCING: Remove unused variables for 2D case. Signed-off-by: Grzegorz Kowal --- sources/forcing.F90 | 84 +++++++++++++++++++++++++++++++++++---------- 1 file changed, 65 insertions(+), 19 deletions(-) diff --git a/sources/forcing.F90 b/sources/forcing.F90 index c9d04fc..69a1509 100644 --- a/sources/forcing.F90 +++ b/sources/forcing.F90 @@ -171,7 +171,9 @@ module forcing ! import external procedures and variables ! +#if NDIMS == 3 use constants , only : pi2 +#endif /* NDIMS == 3 */ use iso_fortran_env, only : error_unit use parameters , only : get_parameter use random , only : randuni, randnorz @@ -191,9 +193,15 @@ module forcing character(len=64) :: injection = "none" character(len=64) :: profile_type = "gauss" character(len=64) :: profile_energy = "off" - integer :: i, j, k = 0, l, k2 + integer :: i, j, l, k2 +#if NDIMS == 3 + integer :: k = 0 +#endif /* NDIMS == 3 */ real(kind=8) :: kl2, ku2, kv2, kv - real(kind=8) :: fa, fi, uu, phi + real(kind=8) :: fa, fi, uu +#if NDIMS == 3 + real(kind=8) :: phi +#endif /* NDIMS == 3 */ ! local vectors ! @@ -913,7 +921,9 @@ module forcing ! local variables ! integer :: ni, n +#if NDIMS == 3 real(kind=8) :: tmp +#endif /* NDIMS == 3 */ real(kind=8), dimension(3) :: xp, ap ! !------------------------------------------------------------------------------- @@ -993,7 +1003,9 @@ module forcing ! import external procedures and variables ! +#if NDIMS == 3 use constants, only : pi2 +#endif /* NDIMS == 3 */ use random , only : randuni, randnorz ! local variables are not implicit by default @@ -1007,7 +1019,10 @@ module forcing ! local variables ! integer :: l - real(kind=8) :: acoeff, dcoeff, phi + real(kind=8) :: acoeff, dcoeff +#if NDIMS == 3 + real(kind=8) :: phi +#endif /* NDIMS == 3 */ real(kind=8) :: dinj ! local vectors @@ -1108,8 +1123,15 @@ module forcing ! local variables ! integer :: l - real(kind=8) :: th1, th2, phi, psi, ga, gb, dinj, sqdt - complex(kind=8) :: aran, bran, xi1, xi2 + real(kind=8) :: th1, dinj, sqdt +#if NDIMS == 3 + real(kind=8) :: th2, phi, psi, ga, gb +#endif /* NDIMS == 3 */ + complex(kind=8) :: aran +#if NDIMS == 3 + complex(kind=8) :: bran + complex(kind=8) :: xi1, xi2 +#endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- ! @@ -1284,8 +1306,10 @@ module forcing ! use blocks , only : block_data use coordinates, only : nm => bcells - use coordinates, only : ax, ay, az - use coordinates, only : xlen, ylen, zlen + use coordinates, only : ax, ay, xlen, ylen +#if NDIMS == 3 + use coordinates, only : az, zlen +#endif /* NDIMS == 3 */ use coordinates, only : periodic use equations , only : idn, imx, imy, imz, ien @@ -1301,7 +1325,10 @@ module forcing ! local variables ! integer :: i, j, k = 1 - real(kind=8) :: x2, y2, z2, r2 + real(kind=8) :: x2, y2, r2 +#if NDIMS == 3 + real(kind=8) :: z2 +#endif /* NDIMS == 3 */ real(kind=8) :: fx, fy, fz, fp, e1, e2 ! local arrays @@ -1503,7 +1530,10 @@ module forcing use blocks , only : block_data use constants , only : pi2 use coordinates, only : nm => bcells, nb, ne - use coordinates, only : ax, ay, az, advol + use coordinates, only : ax, ay, advol +#if NDIMS == 3 + use coordinates, only : az +#endif /* NDIMS == 3 */ use equations , only : idn, imx, imy, imz, ien ! local variables are not implicit by default @@ -1518,14 +1548,19 @@ module forcing ! local variables ! integer :: i, j, k = 1, l, n - real(kind=8) :: cs, sn, tt + real(kind=8) :: cs, sn +#if NDIMS == 3 + real(kind=8) :: tt +#endif /* NDIMS == 3 */ real(kind=8) :: dvol ! local arrays ! - real(kind=8), dimension(nm):: x, y, z, kx, ky, kz - real(kind=8), dimension(nm):: snkx, snky, snkz - real(kind=8), dimension(nm):: cskx, csky, cskz + real(kind=8), dimension(nm):: x, y, z + real(kind=8), dimension(nm):: kx, ky, snkx, snky, cskx, csky +#if NDIMS == 3 + real(kind=8), dimension(nm):: kz, snkz, cskz +#endif /* NDIMS == 3 */ #if NDIMS == 3 real(kind=8), dimension(3,nm,nm,nm) :: acc real(kind=8), dimension( nm,nm,nm) :: den @@ -1719,8 +1754,14 @@ module forcing use blocks , only : block_data use constants , only : pi2 use coordinates, only : nm => bcells, nb, ne - use coordinates, only : ax, ay, az, advol - use equations , only : ivx, ivy, ivz + use coordinates, only : ax, ay, advol +#if NDIMS == 3 + use coordinates, only : az +#endif /* NDIMS == 3 */ + use equations , only : ivx, ivy +#if NDIMS == 3 + use equations , only : ivz +#endif /* NDIMS == 3 */ ! local variables are not implicit by default ! @@ -1733,14 +1774,19 @@ module forcing ! local variables ! integer :: i, j, k = 1, l - real(kind=8) :: cs, sn, tt, dvol + real(kind=8) :: cs, sn, dvol +#if NDIMS == 3 + real(kind=8) :: tt +#endif /* NDIMS == 3 */ complex(kind=8) :: cf ! local arrays ! - real(kind=8), dimension(nm):: x, y, z, kx, ky, kz - real(kind=8), dimension(nm):: snkx, snky, snkz - real(kind=8), dimension(nm):: cskx, csky, cskz + real(kind=8), dimension(nm):: x, y, z + real(kind=8), dimension(nm):: kx, ky, snkx, snky, cskx, csky +#if NDIMS == 3 + real(kind=8), dimension(nm):: kz, snkz, cskz +#endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- ! From 17395086fbe559c2bbb23fab3568ee53c7619dd9 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 15 Aug 2020 01:27:08 -0300 Subject: [PATCH 05/12] REFINEMENT: Remove unused variables for 2D case. Signed-off-by: Grzegorz Kowal --- sources/refinement.F90 | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/sources/refinement.F90 b/sources/refinement.F90 index aa46166..64bb3ac 100644 --- a/sources/refinement.F90 +++ b/sources/refinement.F90 @@ -502,8 +502,14 @@ module refinement ! integer :: i, im1, ip1 integer :: j, jm1, jp1 - integer :: k = 1, km1, kp1 - real(kind=8) :: fl, fr, fc, fx, fy, fz + integer :: k = 1 +#if NDIMS == 3 + integer :: km1, kp1 +#endif /* NDIMS == 3 */ + real(kind=8) :: fl, fr, fc, fx, fy +#if NDIMS == 3 + real(kind=8) :: fz +#endif /* NDIMS == 3 */ ! local parameters ! From 8f33e4aeb55c345a641280b31dbe0dba955d5d32 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 15 Aug 2020 01:28:28 -0300 Subject: [PATCH 06/12] EQUATIONS: Remove unused variables for 2D case. Signed-off-by: Grzegorz Kowal --- sources/equations.F90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/sources/equations.F90 b/sources/equations.F90 index ddb888e..0449be2 100644 --- a/sources/equations.F90 +++ b/sources/equations.F90 @@ -1564,7 +1564,10 @@ module equations integer :: n, p, nc, np integer :: i = 1, il = 1, iu = 1 integer :: j = 1, jl = 1, ju = 1 - integer :: k = 1, kl = 1, ku = 1 + integer :: k = 1 +#if NDIMS == 3 + integer :: kl = 1, ku = 1 +#endif /* NDIMS == 3 */ ! temporary arrays ! From b8e71a4fb0a981d323db4981faca64592f506496 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 15 Aug 2020 01:33:43 -0300 Subject: [PATCH 07/12] INTERPOLATIONS: Remove unused variables in 2D case. Signed-off-by: Grzegorz Kowal --- sources/interpolations.F90 | 31 +++++++++++++++++++++++-------- 1 file changed, 23 insertions(+), 8 deletions(-) diff --git a/sources/interpolations.F90 b/sources/interpolations.F90 index 0db9c30..08e637f 100644 --- a/sources/interpolations.F90 +++ b/sources/interpolations.F90 @@ -651,8 +651,14 @@ module interpolations ! local variables ! logical :: flag - integer :: i, j, i1, j1, k1, i2, j2, k2 - real(kind=max_prec) :: sig, fc, fx, fy, fz, xl, xr, yl, yr, zl, zr + integer :: i, j, i1, j1, i2, j2 +#if NDIMS == 3 + integer :: k1, k2 +#endif /* NDIMS == 3 */ + real(kind=max_prec) :: sig, fc, fx, fy, xl, xr, yl, yr +#if NDIMS == 3 + real(kind=max_prec) :: fz, zl, zr +#endif /* NDIMS == 3 */ ! local arrays for derivatives ! @@ -824,7 +830,10 @@ module interpolations ! local variables ! integer :: i , j , k = 1 - integer :: im1, jm1, km1, ip1, jp1, kp1 + integer :: im1, jm1, ip1, jp1 +#if NDIMS == 3 + integer :: km1, kp1 +#endif /* NDIMS == 3 */ ! local vectors ! @@ -1067,10 +1076,13 @@ module interpolations ! local variables ! - logical :: flag - integer :: i , j , k = 1 - integer :: il , jl , kl , iu , ju , ku - integer :: im1, jm1, km1, ip1, jp1, kp1 + logical :: flag + integer :: i, il, iu, im1, ip1 + integer :: j, jl, ju, jm1, jp1 + integer :: k = 1 +#if NDIMS == 3 + integer :: kl, ku, km1, kp1 +#endif /* NDIMS == 3 */ ! local arrays for derivatives ! @@ -1259,7 +1271,10 @@ module interpolations ! integer :: i, im1, ip1 integer :: j, jm1, jp1 - integer :: k = 1, km1, kp1 + integer :: k = 1 +#if NDIMS == 3 + integer :: km1, kp1 +#endif /* NDIMS == 3 */ integer :: m #if NDIMS == 3 integer :: n, np1, np2 From 14635cb2086e0793cd56f03337e21594612406a5 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 15 Aug 2020 16:25:12 -0300 Subject: [PATCH 08/12] CMake: Update home page. Signed-off-by: Grzegorz Kowal --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f2a36c4..4ee5a1e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,7 @@ cmake_minimum_required(VERSION 3.13) project(amun VERSION 1.0 DESCRIPTION "AMUN Code - a code to perform numerical simulations in astrophysics" - HOMEPAGE_URL "https://gitlab.com/gkowal/amun-code" + HOMEPAGE_URL "https://amuncode.org" LANGUAGES Fortran ) From fda340b21b51d1ca59dac5b4880577ed5ddd8b55 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 16 Aug 2020 16:41:28 -0300 Subject: [PATCH 09/12] EQUATIONS: Rename gamma to adiabatic_index. Keyword 'gamma' is reserved in Fortran for calculating Gamma function. Signed-off-by: Grzegorz Kowal --- sources/equations.F90 | 32 ++++++++++++++++---------------- sources/integrals.F90 | 8 ++++---- sources/io.F90 | 26 +++++++++++++------------- sources/problems.F90 | 8 ++++---- sources/schemes.F90 | 2 +- sources/shapes.F90 | 4 ++-- sources/user_problem.F90 | 4 ++-- 7 files changed, 42 insertions(+), 42 deletions(-) diff --git a/sources/equations.F90 b/sources/equations.F90 index 0449be2..abe2587 100644 --- a/sources/equations.F90 +++ b/sources/equations.F90 @@ -170,7 +170,7 @@ module equations ! adiabatic heat ratio ! - real(kind=8), save :: gamma = 5.0d+00 / 3.0d+00 + real(kind=8), save :: adiabatic_index = 5.0d+00 / 3.0d+00 ! additional adiabatic parameters ! @@ -239,7 +239,7 @@ module equations public :: eigensystem_roe public :: update_primitive_variables public :: fix_unphysical_cells, correct_unphysical_states - public :: gamma, relativistic, magnetized + public :: adiabatic_index, relativistic, magnetized public :: csnd, csnd2 public :: cmax, cmax2 public :: nv, nf, ns @@ -1121,13 +1121,13 @@ module equations ! obtain the adiabatic specific heat ratio ! - call get_parameter("gamma", gamma) + call get_parameter("adiabatic_index", adiabatic_index) ! calculate additional parameters ! - gammam1 = gamma - 1.0d+00 + gammam1 = adiabatic_index - 1.0d+00 gammam1i = 1.0d+00 / gammam1 - gammaxi = gammam1 / gamma + gammaxi = gammam1 / adiabatic_index ! obtain the isothermal sound speed ! @@ -1183,7 +1183,7 @@ module equations ! calculate the sonic Mach number factor ! - msfac = 1.0d+00 / (gamma * msmax**2) + msfac = 1.0d+00 / (adiabatic_index * msmax**2) ! get the tolerance ! @@ -2369,7 +2369,7 @@ module equations do i = 1, size(q,2) - cs = sqrt(gamma * q(ipr,i) / q(idn,i)) + cs = sqrt(adiabatic_index * q(ipr,i) / q(idn,i)) c(1,i) = q(ivx,i) - cs c(2,i) = q(ivx,i) + cs @@ -2452,7 +2452,7 @@ module equations ! calculate the adiabatic speed of sound ! - c = sqrt(gamma * qq(ipr,i,j,k) / qq(idn,i,j,k)) + c = sqrt(adiabatic_index * qq(ipr,i,j,k) / qq(idn,i,j,k)) ! calculate the maximum speed ! @@ -3483,7 +3483,7 @@ module equations do i = 1, size(q,2) - fa = gamma * q(ipr,i) + fa = adiabatic_index * q(ipr,i) fb = fa + bb(i) fc = max(0.0d+00, fb * fb - 4.0d+00 * fa * bx2(i)) cf = sqrt(max(0.5d+00 * (fb + sqrt(fc)), bb(i)) / q(idn,i)) @@ -3569,7 +3569,7 @@ module equations ! calculate the fast magnetosonic speed ! - c = sqrt((gamma * qq(ipr,i,j,k) + bb) / qq(idn,i,j,k)) + c = sqrt((adiabatic_index * qq(ipr,i,j,k) + bb) / qq(idn,i,j,k)) ! calculate the maximum of speed ! @@ -3665,7 +3665,7 @@ module equations ! prepare coefficients ! - gammam2 = gamma - 2.0d+00 + gammam2 = adiabatic_index - 2.0d+00 ! reset all elements ! @@ -4187,7 +4187,7 @@ module equations do i = 1, size(q,2) ww = q(idn,i) + q(ipr,i) / gammaxi - c2 = gamma * q(ipr,i) / ww + c2 = adiabatic_index * q(ipr,i) / ww vv = sum(q(ivx:ivz,i) * q(ivx:ivz,i)) ss = c2 * (1.0d+00 - vv) / (1.0d+00 - c2) fc = 1.0d+00 + ss @@ -4274,7 +4274,7 @@ module equations ! calculate the square of the sound speed ! ww = qq(idn,i,j,k) + qq(ipr,i,j,k) / gammaxi - c2 = gamma * qq(ipr,i,j,k) / ww + c2 = adiabatic_index * qq(ipr,i,j,k) / ww ss = c2 * (1.0d+00 - vv) / (1.0d+00 - c2) fc = 1.0d+00 + ss cc = sqrt(ss * (fc - vv)) @@ -5428,7 +5428,7 @@ module equations ! prepare parameters for this case ! - c2 = gamma * q(ipr,i) / rh + c2 = adiabatic_index * q(ipr,i) / rh v1 = abs(q(ivx,i)) v2 = v1 * v1 @@ -5462,7 +5462,7 @@ module equations ! ! prepare parameters for this case ! - c2 = gamma * q(ipr,i) / rh + c2 = adiabatic_index * q(ipr,i) / rh cc = (1.0d+00 - c2) / vm(i) gn = b2(i) - c2 * vb(i) * vb(i) @@ -5488,7 +5488,7 @@ module equations rh = q(idn,i) + q(ipr,i) / gammaxi vs = sqrt(vm(i)) rt = rh + b2(i) - c2 = gamma * q(ipr,i) / rh + c2 = adiabatic_index * q(ipr,i) / rh ca = (q(ibx,i) * vs + vb(i) * q(ivx,i) / vs)**2 ! prepare polynomial coefficients diff --git a/sources/integrals.F90 b/sources/integrals.F90 index 1f805df..8dc78e0 100644 --- a/sources/integrals.F90 +++ b/sources/integrals.F90 @@ -323,7 +323,7 @@ module integrals use coordinates , only : advol, voli use equations , only : idn, ipr, ivx, ivz, ibx, iby, ibz, ibp use equations , only : ien, imx, imy, imz - use equations , only : magnetized, gamma, csnd + use equations , only : magnetized, adiabatic_index, csnd use evolution , only : step, time, dtn use forcing , only : einj, rinj, arms #ifdef MPI @@ -533,11 +533,11 @@ module integrals ! get average, minimum and maximum values of sonic Mach number ! if (ipr > 0) then - tmp(:,:,:) = sqd(:,:,:) * vel(:,:,:) & + tmp(:,:,:) = sqd(:,:,:) * vel(:,:,:) / & #if NDIMS == 3 - / sqrt(gamma * pdata%q(ipr,nb:ne,nb:ne,nb:ne)) + sqrt(adiabatic_index * pdata%q(ipr,nb:ne,nb:ne,nb:ne)) #else /* NDIMS == 3 */ - / sqrt(gamma * pdata%q(ipr,nb:ne,nb:ne, : )) + sqrt(adiabatic_index * pdata%q(ipr,nb:ne,nb:ne, : )) #endif /* NDIMS == 3 */ else tmp(:,:,:) = vel(:,:,:) / csnd diff --git a/sources/io.F90 b/sources/io.F90 index 07f8d69..ef07185 100644 --- a/sources/io.F90 +++ b/sources/io.F90 @@ -2470,7 +2470,7 @@ module io use coordinates , only : zmin, zmax #endif /* NDIMS == 3 */ use coordinates , only : bdims => domain_base_dims - use equations , only : eqsys, eos, nv, pvars, gamma, csnd + use equations , only : eqsys, eos, nv, pvars, adiabatic_index, csnd use evolution , only : step, time, dt, dtn use iso_fortran_env, only : error_unit use mpitools , only : nprocs, nproc @@ -2565,20 +2565,20 @@ module io write(lun,"(a)") '' write(lun,"(a)") '' write(lun,"(a)") '' - call write_attribute_xml(lun, "problem" , problem_name) + call write_attribute_xml(lun, "problem" , problem_name) write(lun,"(a)") '' write(lun,"(a)") '' - call write_attribute_xml(lun, "nprocs" , nprocs) - call write_attribute_xml(lun, "nproc" , nproc) + call write_attribute_xml(lun, "nprocs" , nprocs) + call write_attribute_xml(lun, "nproc" , nproc) write(lun,"(a)") '' write(lun,"(a)") '' - call write_attribute_xml(lun, "eqsys" , eqsys) - call write_attribute_xml(lun, "eos" , eos) - call write_attribute_xml(lun, "nvars" , nv) - call write_attribute_xml(lun, "gamma" , gamma) - call write_attribute_xml(lun, "csnd" , csnd) - call write_attribute_xml(lun, "viscosity" , viscosity) - call write_attribute_xml(lun, "resistivity", resistivity) + call write_attribute_xml(lun, "eqsys" , eqsys) + call write_attribute_xml(lun, "eos" , eos) + call write_attribute_xml(lun, "nvars" , nv) + call write_attribute_xml(lun, "adiabatic_index", adiabatic_index) + call write_attribute_xml(lun, "csnd" , csnd) + call write_attribute_xml(lun, "viscosity" , viscosity) + call write_attribute_xml(lun, "resistivity" , resistivity) write(lun,"(a)") '' write(lun,"(a)") '' call write_attribute_xml(lun, "ndims" , NDIMS) @@ -3931,7 +3931,7 @@ module io use coordinates , only : bdims => domain_base_dims use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax use coordinates , only : periodic - use equations , only : eqsys, eos, gamma, csnd + use equations , only : eqsys, eos, adiabatic_index, csnd use evolution , only : step, time, dt, dtn use forcing , only : nmodes, einj, fcoefs use hdf5 , only : hid_t @@ -4033,7 +4033,7 @@ module io call write_attribute(gid, 'dt' , dt ) call write_attribute(gid, 'dtn' , dtn ) if (eos == 'adi') then - call write_attribute(gid, 'gamma', gamma) + call write_attribute(gid, 'adiabatic_index', adiabatic_index) end if if (eos == 'iso') then call write_attribute(gid, 'csnd' , csnd ) diff --git a/sources/problems.F90 b/sources/problems.F90 index 37de4ce..99c2e98 100644 --- a/sources/problems.F90 +++ b/sources/problems.F90 @@ -380,7 +380,7 @@ module problems use coordinates, only : az, adz #endif /* NDIMS == 3 */ use equations , only : prim2cons - use equations , only : gamma + use equations , only : adiabatic_index use equations , only : nv use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp use parameters , only : get_parameter @@ -487,7 +487,7 @@ module problems ! calculate parallel and perpendicular pressures from sound speeds ! - pr_amb = dens * csnd * csnd / gamma + pr_amb = dens * csnd * csnd / adiabatic_index pr_ovr = pr_amb * ratio else @@ -853,7 +853,7 @@ module problems use coordinates, only : az, adz #endif /* NDIMS == 3 */ use equations , only : prim2cons - use equations , only : gamma + use equations , only : adiabatic_index use equations , only : nv use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp use parameters , only : get_parameter @@ -965,7 +965,7 @@ module problems dn_amb = dens dn_ovr = dn_amb pr_amb = pres - pr_ovr = (gamma - 1.0d+00) * eexp / dvol + pr_ovr = (adiabatic_index - 1.0d+00) * eexp / dvol ! calculate initial uniform field components ! diff --git a/sources/schemes.F90 b/sources/schemes.F90 index 8d7bcbf..ed247cb 100644 --- a/sources/schemes.F90 +++ b/sources/schemes.F90 @@ -1847,7 +1847,7 @@ module schemes vv = sum(vs(1:3) * vs( 1: 3)) vb = sum(vs(1:3) * us(ibx:ibz)) -! calculate inverse gamma +! calculate inverse of Lorentz factor ! gi = 1.0d+00 - vv diff --git a/sources/shapes.F90 b/sources/shapes.F90 index 75d55c6..ad3fa18 100644 --- a/sources/shapes.F90 +++ b/sources/shapes.F90 @@ -357,7 +357,7 @@ module shapes use coordinates , only : az, adz #endif /* NDIMS == 3 */ use equations , only : prim2cons - use equations , only : gamma + use equations , only : adiabatic_index use equations , only : nv use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp use parameters , only : get_parameter @@ -432,7 +432,7 @@ module shapes ! if (ipr > 0) then dn_ovr = dens - pr_ovr = dens * ratio * csnd * csnd / gamma + pr_ovr = dens * ratio * csnd * csnd / adiabatic_index else dn_ovr = dens * ratio end if diff --git a/sources/user_problem.F90 b/sources/user_problem.F90 index ba3558d..66ca9f6 100644 --- a/sources/user_problem.F90 +++ b/sources/user_problem.F90 @@ -80,7 +80,7 @@ module user_problem ! include external procedures and variables ! - use equations , only : eos, csnd, csnd2, gamma + use equations , only : eos, csnd, csnd2, adiabatic_index use helpers , only : print_section, print_parameter use parameters, only : get_parameter @@ -119,7 +119,7 @@ module user_problem if (eos == "iso") then pres = csnd2 * dens else - pres = csnd2 * dens / gamma + pres = csnd2 * dens / adiabatic_index end if bmag = sqrt(2.0d+00 * pres / beta) From 6bf3afcf0b681b37e0c769eecc506f520c50f46b Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 16 Aug 2020 16:46:27 -0300 Subject: [PATCH 10/12] EQUATIONS: Rename parameter 'csnd' to 'sound_speed'. Signed-off-by: Grzegorz Kowal --- sources/equations.F90 | 2 +- sources/io.F90 | 4 ++-- sources/problems.F90 | 12 ++++++------ sources/shapes.F90 | 12 ++++++------ 4 files changed, 15 insertions(+), 15 deletions(-) diff --git a/sources/equations.F90 b/sources/equations.F90 index abe2587..37d7511 100644 --- a/sources/equations.F90 +++ b/sources/equations.F90 @@ -1131,7 +1131,7 @@ module equations ! obtain the isothermal sound speed ! - call get_parameter("csnd" , csnd ) + call get_parameter("sound_speed", csnd ) ! calculate additional parameters ! diff --git a/sources/io.F90 b/sources/io.F90 index ef07185..5e20609 100644 --- a/sources/io.F90 +++ b/sources/io.F90 @@ -2576,7 +2576,7 @@ module io call write_attribute_xml(lun, "eos" , eos) call write_attribute_xml(lun, "nvars" , nv) call write_attribute_xml(lun, "adiabatic_index", adiabatic_index) - call write_attribute_xml(lun, "csnd" , csnd) + call write_attribute_xml(lun, "sound_speed" , csnd) call write_attribute_xml(lun, "viscosity" , viscosity) call write_attribute_xml(lun, "resistivity" , resistivity) write(lun,"(a)") '' @@ -4036,7 +4036,7 @@ module io call write_attribute(gid, 'adiabatic_index', adiabatic_index) end if if (eos == 'iso') then - call write_attribute(gid, 'csnd' , csnd ) + call write_attribute(gid, 'sound_speed', csnd) end if ! store the vector attributes diff --git a/sources/problems.F90 b/sources/problems.F90 index 99c2e98..96e3442 100644 --- a/sources/problems.F90 +++ b/sources/problems.F90 @@ -462,12 +462,12 @@ module problems ! get problem parameters ! - call get_parameter("dens" , dens ) - call get_parameter("ratio" , ratio ) - call get_parameter("radius", radius) - call get_parameter("csnd" , csnd ) - call get_parameter("buni" , buni ) - call get_parameter("angle" , angle ) + call get_parameter("dens" , dens ) + call get_parameter("ratio" , ratio ) + call get_parameter("radius" , radius) + call get_parameter("sound_speed", csnd ) + call get_parameter("buni" , buni ) + call get_parameter("angle" , angle ) #if NDIMS == 3 ! get the fine grid resolution diff --git a/sources/shapes.F90 b/sources/shapes.F90 index ad3fa18..98ebc83 100644 --- a/sources/shapes.F90 +++ b/sources/shapes.F90 @@ -421,12 +421,12 @@ module shapes ! get problem parameters ! - call get_parameter("dens" , dens ) - call get_parameter("ratio" , ratio ) - call get_parameter("radius", radius) - call get_parameter("csnd" , csnd ) - call get_parameter("buni" , buni ) - call get_parameter("angle" , angle ) + call get_parameter("dens" , dens ) + call get_parameter("ratio" , ratio ) + call get_parameter("radius" , radius) + call get_parameter("sound_speed", csnd ) + call get_parameter("buni" , buni ) + call get_parameter("angle" , angle ) ! set the conditions inside the radius ! From 2327e6ffa53830d2aec131a4d554c274ca6a1cb1 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 16 Aug 2020 16:49:24 -0300 Subject: [PATCH 11/12] PROBLEMS: Update parameters for sound speed and adiabatic index. Signed-off-by: Grzegorz Kowal --- problems/blast.in | 4 ++-- problems/implosion.in | 2 +- problems/rayleigh-taylor.in | 2 +- problems/sedov-taylor.in | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/problems/blast.in b/problems/blast.in index 551ef8b..3c63d15 100644 --- a/problems/blast.in +++ b/problems/blast.in @@ -6,8 +6,8 @@ ratio = 1.0000d+02 radius = 1.0000d-01 buni = 0.0000d+00 angle = 4.5000d+01 -csnd = 4.0825d-01 -gamma = 1.6667d+00 +sound_speed = 4.0825d-01 +adiabatic_index = 1.6667d+00 # physics # diff --git a/problems/implosion.in b/problems/implosion.in index 11f94e6..5da46e0 100644 --- a/problems/implosion.in +++ b/problems/implosion.in @@ -1,7 +1,7 @@ # problem name and parameters # problem = "implosion" -gamma = 1.40d+00 +adiabatic_index = 1.40d+00 shock_line = 5.00d-01 # physics diff --git a/problems/rayleigh-taylor.in b/problems/rayleigh-taylor.in index 3751bcf..b9f0b9d 100644 --- a/problems/rayleigh-taylor.in +++ b/problems/rayleigh-taylor.in @@ -1,7 +1,7 @@ # problem name and parameters # problem = "rayleigh-taylor" -gamma = 1.4d+00 +adiabatic_index = 1.4d+00 # random number generator parameters # diff --git a/problems/sedov-taylor.in b/problems/sedov-taylor.in index 990aae0..313124a 100644 --- a/problems/sedov-taylor.in +++ b/problems/sedov-taylor.in @@ -1,7 +1,7 @@ # problem name and parameters # problem = "sedov-taylor" -gamma = 1.4d+00 +adiabatic_index = 1.4d+00 # physics # From 5886c2039a2ef098a39b2a4a91dd1cf72130bae3 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 16 Aug 2020 17:00:50 -0300 Subject: [PATCH 12/12] PYTHON: Update adiabatic index name. Signed-off-by: Grzegorz Kowal --- python/amunpy.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/amunpy.py b/python/amunpy.py index b6d17e8..c93d1f1 100644 --- a/python/amunpy.py +++ b/python/amunpy.py @@ -149,7 +149,7 @@ class AmunXML: variables.append('cury') variables.append('curz') variables.append('curr') - if 'pres' in variables and 'gamma' in self.attributes: + if 'pres' in variables and 'adiabatic_index' in self.attributes: variables.append('eint') if all(v in variables for v in ['dens','pres']): variables.append('temp') @@ -513,7 +513,7 @@ class AmunXML: dset = np.sqrt(dset) elif var == 'eint': dset = self.read_binary_data(n, 'pres') - dset *= 1.0 / (self.attributes('gamma') - 1.0) + dset *= 1.0 / (self.attributes('adiabatic_index') - 1) dset = np.reshape(dset, cm) elif var == 'temp': dset = self.read_binary_data(n, 'pres') @@ -530,7 +530,7 @@ class AmunXML: tmp = self.read_binary_data(n, 'dens') dset *= tmp tmp = self.read_binary_data(n, 'pres') - dset += 2.0 / (self.attributes('gamma') - 1.0) * tmp + dset += 2.0 / (self.attributes('adiabatic_index') - 1) * tmp if 'magn' in self.variables: tmp = self.read_binary_data(n, 'magx') dset = tmp**2 @@ -925,7 +925,7 @@ def amun_dataset(fname, vname, shrink = 1, interpolation = 'rebin', progress = F nc = amun_attribute(fname, 'nprocs') nl = amun_attribute(fname, 'nleafs') if eos == 'adi': - gm = amun_attribute(fname, 'gamma') + gm = amun_attribute(fname, 'adiabatic_index') # get block dimensions and the maximum level #