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

SOURCES: Plug the GRAVITY module in.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2017-03-06 11:12:25 -03:00
parent 28be951d3f
commit 9256b58905
2 changed files with 14 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/makefile
View File

@ -159,8 +159,8 @@ schemes.o : schemes.F90 algebra.o coordinates.o equations.o \
interpolations.o timers.o
shapes.o : shapes.F90 blocks.o constants.o coordinates.o equations.o \
parameters.o timers.o
sources.o : sources.F90 blocks.o coordinates.o equations.o operators.o \
parameters.o timers.o
sources.o : sources.F90 blocks.o coordinates.o equations.o gravity.o \
operators.o parameters.o timers.o
timers.o : timers.F90
utils.o : utils.F90 error.o

39
src/sources.F90
View File

@ -49,10 +49,6 @@ module sources
!
integer , save :: glm_type = 0
! gravitational acceleration coefficient
!
real(kind=8), save :: gpoint = 0.0d+00
! viscosity coefficient
!
real(kind=8), save :: viscosity = 0.0d+00
@ -144,10 +140,6 @@ module sources
glm_type = 0
end select
! get acceleration coefficient
!
call get_parameter_real("gpoint" , gpoint )
! get viscosity coefficient
!
call get_parameter_real("viscosity" , viscosity)
@ -161,7 +153,6 @@ module sources
if (verbose) then
write (*,"(4x,a,1x,a) ") "glm source terms =", trim(tglm)
write (*,"(4x,a,1x,1e9.2)") "point mass constant =", gpoint
write (*,"(4x,a,1x,1e9.2)") "viscosity =", viscosity
write (*,"(4x,a,1x,1e9.2)") "resistivity =", resistivity
@ -242,6 +233,7 @@ module sources
use equations , only : nv, inx, iny, inz
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ien
use equations , only : ibx, iby, ibz, ibp
use gravity , only : gravity_enabled, gravitational_acceleration
use operators , only : divergence, gradient, laplace, curl
! local variables are not implicit by default
@ -257,7 +249,7 @@ module sources
!
integer :: i , j , k
real(kind=8) :: fc, gc
real(kind=8) :: r2, r3, gx, gy, gz
real(kind=8) :: gx, gy, gz
real(kind=8) :: dbx, dby, dbz
real(kind=8) :: dvxdx, dvxdy, dvxdz, divv
real(kind=8) :: dvydx, dvydy, dvydz
@ -265,7 +257,7 @@ module sources
! local arrays
!
real(kind=8), dimension(3) :: dh
real(kind=8), dimension(3) :: ga, dh
real(kind=8), dimension(im) :: x
real(kind=8), dimension(jm) :: y
real(kind=8), dimension(km) :: z
@ -280,9 +272,9 @@ module sources
call start_timer(imu)
#endif /* PROFILE */
! proceed only if the gravitational acceleration coefficient is not zero
! proceed only if the gravitational term is enabled
!
if (gpoint /= 0.0d+00) then
if (gravity_enabled) then
! prepare block coordinates
!
@ -298,23 +290,16 @@ module sources
do j = 1, jm
do i = 1, im
! calculate distance from the origin
! get gravitational acceleration components
!
#if NDIMS == 2
r2 = x(i) * x(i) + y(j) * y(j)
#endif /* NDIMS == 2 */
#if NDIMS == 3
r2 = x(i) * x(i) + y(j) * y(j) + z(k) * z(k)
#endif /* NDIMS == 3 */
r3 = r2 * sqrt(r2)
call gravitational_acceleration(x(i), y(j), z(k), ga(1:3))
! calculate gravitational acceleration factors
! calculate the gravitational source terms
!
gc = gpoint * pdata%q(idn,i,j,k) / max(1.0d-16, r3)
gx = gc * x(i)
gy = gc * y(j)
gx = pdata%q(idn,i,j,k) * ga(1)
gy = pdata%q(idn,i,j,k) * ga(2)
#if NDIMS == 3
gz = gc * z(k)
gz = pdata%q(idn,i,j,k) * ga(3)
#endif /* NDIMS == 3 */
! add source terms to momentum equations
@ -345,7 +330,7 @@ module sources
end do ! j = 1, jm
end do ! k = 1, km
end if ! gpoint is not zero
end if ! gravity enabled
! proceed only if the viscosity coefficient is not zero
!