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

PROBLEMS: Fix uniform field component calculation.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2014-08-04 16:10:11 -03:00
parent e9468961aa
commit 87b18a2796
1 changed files with 13 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
src/problems.F90
View File

@ -261,12 +261,14 @@ module problems
logical , save :: first = .true. logical , save :: first = .true.
real(kind=8), save :: dn_amb, dn_ovr real(kind=8), save :: dn_amb, dn_ovr
real(kind=8), save :: pr_amb, pr_ovr real(kind=8), save :: pr_amb, pr_ovr
real(kind=8), save :: bx, by
real(kind=8), save :: r2 real(kind=8), save :: r2
! local variables ! local variables
! !
integer :: i, j, k, ic, jc, kc integer :: i, j, k, ic, jc, kc
real(kind=8) :: xl, yl, zl, xu, yu, zu, rl, ru real(kind=8) :: xl, yl, zl, xu, yu, zu, rl, ru
real(kind=8) :: sn
#if NDIMS == 3 #if NDIMS == 3
real(kind=8) :: xb, yb, zb real(kind=8) :: xb, yb, zb
real(kind=8) :: xt, yt, zt real(kind=8) :: xt, yt, zt
@ -275,7 +277,7 @@ module problems
real(kind=8) :: rlu, rul real(kind=8) :: rlu, rul
real(kind=8) :: xb, yb real(kind=8) :: xb, yb
real(kind=8) :: xt, yt real(kind=8) :: xt, yt
real(kind=8) :: sn, ph real(kind=8) :: ph
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
real(kind=8) :: dx, dy, dz, dxh, dyh, dzh, dvol real(kind=8) :: dx, dy, dz, dxh, dyh, dzh, dvol
real(kind=8) :: fc_amb, fc_ovr real(kind=8) :: fc_amb, fc_ovr
@ -340,6 +342,14 @@ module problems
dn_ovr = dn_amb * ratio dn_ovr = dn_amb * ratio
end if end if
! calculate initial uniform field components
!
if (ibx > 0) then
sn = sin(d2r * angle)
bx = buni * sqrt(1.0d+00 - sn * sn)
by = buni * sn
end if
! calculate the square of radius ! calculate the square of radius
! !
r2 = radius * radius r2 = radius * radius
@ -416,8 +426,8 @@ module problems
! !
if (ibx > 0) then if (ibx > 0) then
q(ibx,:) = buni * cos(d2r * angle) q(ibx,:) = bx
q(iby,:) = buni * sin(d2r * angle) q(iby,:) = by
q(ibz,:) = 0.0d+00 q(ibz,:) = 0.0d+00
q(ibp,:) = 0.0d+00 q(ibp,:) = 0.0d+00