Compare commits
213 Commits
Author | SHA1 | Date | |
---|---|---|---|
ae0c81c575 | |||
5f1902c7dd | |||
161595791d | |||
8048597303 | |||
c353dc4ae5 | |||
c44eeb31a0 | |||
8fc234269b | |||
7c59e1c0bc | |||
6f4e800932 | |||
3de0ec092d | |||
8c340364bd | |||
dde07cf6bc | |||
c89e5fe002 | |||
f79846d1bc | |||
ef9c94d006 | |||
e15fcc2f91 | |||
1ab20ef426 | |||
7376abea36 | |||
2afa59d16c | |||
78c6d1d83f | |||
96cc183c10 | |||
795b28e426 | |||
ab3667c34d | |||
364696661a | |||
0ea3ef0841 | |||
d0f5d0c888 | |||
e59470ea30 | |||
6b3f1d74ea | |||
4616a9144a | |||
a77903d51c | |||
37b9a957ee | |||
7cf63ec5d5 | |||
d8c86a0f48 | |||
54eeca1d32 | |||
493625ce47 | |||
9cf9a57933 | |||
6aa20caf88 | |||
68683937d7 | |||
aa261bd9b8 | |||
8599e0ede3 | |||
86469c2043 | |||
73ae99ff02 | |||
44fe494299 | |||
f846d6183b | |||
f0df248bfa | |||
b229885aa9 | |||
825662e466 | |||
c577680f80 | |||
8cb77828c1 | |||
69a027c8f5 | |||
c97c443f20 | |||
19e832cf68 | |||
4dfdcb66c5 | |||
fccdff630f | |||
f55e6bd773 | |||
8a3df6b216 | |||
b0963e3e61 | |||
24bfdb9591 | |||
60c674b887 | |||
15714dedef | |||
7193a96dba | |||
8be69137b1 | |||
03cdfa1a11 | |||
e5ff27373e | |||
4c9553ae69 | |||
aa34f0aa50 | |||
2c3b48f5ef | |||
78ff4f99c4 | |||
552bbeffa1 | |||
84b2284699 | |||
c2d0fe8a3a | |||
7cbb86ff98 | |||
be396d072d | |||
52ab497962 | |||
aa1daeee84 | |||
7433177351 | |||
84de161562 | |||
66cff8ab98 | |||
1518dae615 | |||
ea06ef2d9f | |||
5f85e090a9 | |||
41d0e0e48d | |||
c7c68922b7 | |||
8959948437 | |||
6b69c67e80 | |||
59d470f8e8 | |||
ac9ca20c8a | |||
53e1b4551c | |||
af19185cd7 | |||
3d9c7c6990 | |||
ea01f254ea | |||
86204879e0 | |||
c5f89d42f8 | |||
b863d56e9b | |||
a0ae13b87b | |||
1c8341ac53 | |||
6639495441 | |||
3434305791 | |||
e220bc7cbd | |||
d4494677a7 | |||
cb10f52017 | |||
dc60992354 | |||
6ed0ba096c | |||
7ec79e40aa | |||
0cc641d7d1 | |||
4aee2f3eb0 | |||
c85071406a | |||
63d6dfe0c4 | |||
01cdba7766 | |||
530b5a5945 | |||
216a76cb9a | |||
b1b6940db9 | |||
780d9c2cd1 | |||
13719a8f29 | |||
0cc32eaf80 | |||
170f589550 | |||
fc8aa7e3ef | |||
3fa051604b | |||
352369b57c | |||
5c7706ca18 | |||
e7222201ab | |||
81284a4f9b | |||
53c53a8e4b | |||
ef58bc17a3 | |||
8502126bed | |||
d87814d3f5 | |||
6c4725518f | |||
3127ddd7eb | |||
a07e85fd3b | |||
10179fb71b | |||
cb831885cb | |||
cf404e4b52 | |||
65c0a0a235 | |||
545540a132 | |||
2d02bfff65 | |||
646e2a236c | |||
246ae65000 | |||
12a3216ea4 | |||
3c42fc79a3 | |||
b95377b7c1 | |||
a8e0273969 | |||
a3b31f3141 | |||
8b82d8444c | |||
12c0dd69c3 | |||
8f7a795045 | |||
67af4550db | |||
3c81cf068f | |||
45c7706af8 | |||
7c51d3d096 | |||
24963a7a27 | |||
e226a5ebb6 | |||
dc172efdb3 | |||
68ff188148 | |||
91f8adda59 | |||
b1525655e2 | |||
48450383f7 | |||
8a23933ae2 | |||
60e4941170 | |||
80be187082 | |||
67bc11cd5a | |||
6eee6eca65 | |||
d638ffc513 | |||
943d292ae2 | |||
084df89cdd | |||
7536c3ddba | |||
fdc94d64ea | |||
09bdebd109 | |||
5d07afc2e7 | |||
ec0da99053 | |||
269f77e144 | |||
bba29cd7ed | |||
8fc8173d26 | |||
96c158ea28 | |||
3e05170ba4 | |||
d35956e4d9 | |||
aff5ca21e2 | |||
1689a41d91 | |||
3f904f0df0 | |||
c39095c9b5 | |||
6ea3428f7d | |||
6ec3101a9b | |||
1dc305bb87 | |||
15ea95f932 | |||
32cf0c51ca | |||
d111c7cf50 | |||
7bc5807c80 | |||
5cac79f637 | |||
65be8c3d5a | |||
00860cd62b | |||
c6a54ba2c9 | |||
e11d41fb02 | |||
180d8b8e36 | |||
957ed24caf | |||
c89ef4415f | |||
b8bef681a7 | |||
6bc9493674 | |||
9f08daa277 | |||
39c22fc053 | |||
07da0e7a91 | |||
06fe62faf7 | |||
f7597e66ac | |||
d75bba8773 | |||
36ef146c9c | |||
b71f592697 | |||
a925395ded | |||
ba4ad674fb | |||
a9f90fe5ac | |||
1fa217c4b1 | |||
44ac795aac | |||
03bd0ace96 | |||
944bfad52c | |||
ba1b7d4a7c | |||
a8a296a450 |
74
problems/binaries.in
Normal file
74
problems/binaries.in
Normal file
@ -0,0 +1,74 @@
|
||||
# problem name and parameters
|
||||
#
|
||||
problem = "binaries"
|
||||
gamma = 1.1d+00
|
||||
mstar = 2.0d+00
|
||||
mcomp = 1.0d+00
|
||||
rstar = 5.0d-02
|
||||
rcomp = 5.0d-02
|
||||
dstar = 1.0d+07
|
||||
dratio = 1.0d+02
|
||||
msstar = 3.0d+00
|
||||
mscomp = 1.5d+01
|
||||
tstar = 2.0d+01
|
||||
tcomp = 1.0d+01
|
||||
vstar = 5.0d-01
|
||||
vcomp = 2.5d+00
|
||||
distance = 2.0d-01
|
||||
eccentricity = 0.9d+00
|
||||
period = 1.0d+02
|
||||
buni = 1.0d-03
|
||||
|
||||
# physics
|
||||
#
|
||||
equation_system = "mhd"
|
||||
equation_of_state = "adi"
|
||||
|
||||
# numerical methods
|
||||
#
|
||||
time_advance = "ssprk(m,2)"
|
||||
stages = 4
|
||||
riemann_solver = "hlld"
|
||||
reconstruction = "tvd"
|
||||
limiter = "mc"
|
||||
fix_positivity = "off"
|
||||
clip_extrema = "off"
|
||||
|
||||
# mesh parameters
|
||||
#
|
||||
xblocks = 3
|
||||
yblocks = 2
|
||||
zblocks = 1
|
||||
xmin = -4.5d+00
|
||||
xmax = 4.5d+00
|
||||
ymin = -3.0d+00
|
||||
ymax = 3.0d+00
|
||||
zmin = -1.5d+00
|
||||
zmax = 1.5d+00
|
||||
|
||||
# refinement control
|
||||
#
|
||||
ncells = 16
|
||||
nghosts = 2
|
||||
minlev = 1
|
||||
maxlev = 5
|
||||
|
||||
# boundary conditions
|
||||
#
|
||||
xlbndry = "open"
|
||||
xubndry = "open"
|
||||
ylbndry = "open"
|
||||
yubndry = "open"
|
||||
zlbndry = "open"
|
||||
zubndry = "open"
|
||||
enable_shapes = "on"
|
||||
|
||||
# runtime control parameters
|
||||
#
|
||||
tmax = 1.2d+02
|
||||
cfl = 3.0d-01
|
||||
|
||||
# data output control
|
||||
#
|
||||
precise_snapshots = "on"
|
||||
snapshot_interval = 1.0d-01
|
@ -31,6 +31,56 @@ module user_problem
|
||||
|
||||
implicit none
|
||||
|
||||
! flag indicating the stars are orbiting
|
||||
!
|
||||
logical, save :: orbiting = .false.
|
||||
|
||||
! default problem parameter values
|
||||
!
|
||||
real(kind=8), save :: mstar = 2.00d+00
|
||||
real(kind=8), save :: mcomp = 1.00d+00
|
||||
real(kind=8), save :: rstar = 1.00d-01
|
||||
real(kind=8), save :: rcomp = 1.00d-01
|
||||
real(kind=8), save :: rbuffer = 8.00d-02
|
||||
real(kind=8), save :: dstar = 1.00d+06
|
||||
real(kind=8), save :: dratio = 1.00d+02
|
||||
real(kind=8), save :: msstar = 3.00d+00
|
||||
real(kind=8), save :: mscomp = 1.50d+01
|
||||
real(kind=8), save :: tstar = 0.00d+00
|
||||
real(kind=8), save :: tcomp = 0.00d+00
|
||||
real(kind=8), save :: omstar = 0.00d+00
|
||||
real(kind=8), save :: omcomp = 0.00d+00
|
||||
real(kind=8), save :: vstar = 5.00d-01
|
||||
real(kind=8), save :: vcomp = 2.50d+00
|
||||
real(kind=8), save :: xshock = 0.00d+00
|
||||
real(kind=8), save :: cshock = 0.00d+00
|
||||
real(kind=8), save :: xshift = 0.00d+00
|
||||
real(kind=8), save :: dist = 2.00d-01
|
||||
real(kind=8), save :: period = 1.00d+02
|
||||
real(kind=8), save :: ecc = 9.00d-01
|
||||
real(kind=8), save :: buni = 1.00d-03
|
||||
|
||||
real(kind=8), save :: tol = 1.00d-14
|
||||
integer , save :: maxit = 20
|
||||
|
||||
real(kind=8), save :: dcomp, pstar, pcomp
|
||||
real(kind=8), save :: r2stari, r2compi, r2staro, r2compo
|
||||
real(kind=8), save :: acomp , bcomp
|
||||
real(kind=8), save :: om
|
||||
|
||||
real(kind=8), save :: ms, mc
|
||||
real(kind=8), save :: asemi, omega, ean
|
||||
real(kind=8), save :: uxs, uxc
|
||||
real(kind=8), save :: tprev
|
||||
|
||||
real(kind=8), save :: xps = 1.0d+00, yps = 0.0d+00
|
||||
real(kind=8), save :: xpc = -1.0d+00, ypc = 0.0d+00
|
||||
real(kind=8), save :: uys, uyc, xrs, yrs, xrc, yrc
|
||||
|
||||
real(kind=8) :: ida = 1.0d+00, idm = 1.0d+00
|
||||
|
||||
integer, parameter :: idx = NDIMS - 1
|
||||
|
||||
private
|
||||
public :: initialize_user_problem, finalize_user_problem
|
||||
public :: user_statistics
|
||||
@ -58,6 +108,15 @@ module user_problem
|
||||
!
|
||||
subroutine initialize_user_problem(problem, rcount, verbose, status)
|
||||
|
||||
use algebra , only : quadratic
|
||||
use constants , only : pi2
|
||||
use equations , only : ipr
|
||||
use equations , only : adiabatic_index
|
||||
use helpers , only : print_section, print_parameter
|
||||
use mesh , only : setup_problem
|
||||
use parameters , only : get_parameter
|
||||
use shapes , only : update_shapes
|
||||
|
||||
implicit none
|
||||
|
||||
character(len=64), intent(in) :: problem
|
||||
@ -65,10 +124,205 @@ module user_problem
|
||||
logical , intent(in) :: verbose
|
||||
integer , intent(out) :: status
|
||||
|
||||
character(len=32) :: orb = "off"
|
||||
|
||||
integer :: n
|
||||
real(kind=8) :: as, ac
|
||||
real(kind=8), dimension(3) :: a
|
||||
real(kind=8), dimension(2) :: x
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
status = 0
|
||||
|
||||
! star masses, radia, densities and sonic Mach numbers
|
||||
!
|
||||
call get_parameter("mstar" , mstar )
|
||||
call get_parameter("mcomp" , mcomp )
|
||||
call get_parameter("rstar" , rstar )
|
||||
call get_parameter("rcomp" , rcomp )
|
||||
call get_parameter("rbuffer" , rbuffer)
|
||||
call get_parameter("dstar" , dstar )
|
||||
call get_parameter("dratio" , dratio)
|
||||
call get_parameter("msstar" , msstar)
|
||||
call get_parameter("mscomp" , mscomp)
|
||||
|
||||
! wind speeds
|
||||
!
|
||||
call get_parameter("vstar" , vstar )
|
||||
call get_parameter("vcomp" , vcomp )
|
||||
|
||||
! orbit parameters
|
||||
!
|
||||
call get_parameter("xshift" , xshift)
|
||||
call get_parameter("distance" , dist )
|
||||
call get_parameter("period" , period)
|
||||
call get_parameter("eccentricity", ecc )
|
||||
|
||||
! the initial shock curvature
|
||||
!
|
||||
call get_parameter("shock_curvature", cshock)
|
||||
|
||||
! star rotation periods
|
||||
!
|
||||
call get_parameter("tstar" , tstar )
|
||||
call get_parameter("tcomp" , tcomp )
|
||||
|
||||
! magnetic field profile parameters
|
||||
!
|
||||
call get_parameter("buni" , buni )
|
||||
|
||||
! Kepler's equation solver parameters
|
||||
!
|
||||
call get_parameter("orbiting" , orb )
|
||||
call get_parameter("tolerance" , tol )
|
||||
call get_parameter("maxit" , maxit )
|
||||
|
||||
! check if stars should be orbiting
|
||||
!
|
||||
select case(trim(orb))
|
||||
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
|
||||
orbiting = .true.
|
||||
case default
|
||||
orbiting = .false.
|
||||
end select
|
||||
|
||||
! indices
|
||||
!
|
||||
ida = idx * adiabatic_index
|
||||
idm = 5.0d-01 * NDIMS
|
||||
|
||||
! correct rbuffer
|
||||
!
|
||||
rbuffer = max(0.0d+00, rbuffer)
|
||||
|
||||
! calculate the square of radia
|
||||
!
|
||||
r2stari = rstar**2
|
||||
r2compi = rcomp**2
|
||||
r2staro = (rstar + rbuffer)**2
|
||||
r2compo = (rcomp + rbuffer)**2
|
||||
|
||||
! calculate densities and pressures
|
||||
!
|
||||
dcomp = dstar / dratio
|
||||
if (ipr > 0) then
|
||||
pstar = (vstar / msstar)**2 * dstar / adiabatic_index
|
||||
pcomp = (vcomp / mscomp)**2 * dcomp / adiabatic_index
|
||||
end if
|
||||
|
||||
! calculate the rotation speeds
|
||||
!
|
||||
if (tstar > 0.0d+00) omstar = 1.0d+00 / tstar
|
||||
if (tcomp > 0.0d+00) omcomp = 1.0d+00 / tcomp
|
||||
|
||||
! calculate orbit parameters
|
||||
!
|
||||
acomp = dist / (1.0d+00 - ecc)
|
||||
bcomp = acomp * sqrt(1.0d+00 - ecc * ecc)
|
||||
om = (pi2 / period) / (1.0d+00 - ecc)
|
||||
|
||||
if (.not. orbiting) then
|
||||
|
||||
! fix positions if the start are not orbiting
|
||||
!
|
||||
xps = xshift + 0.5d+00 * acomp
|
||||
yps = 0.0d+00
|
||||
xpc = xshift - 0.5d+00 * acomp
|
||||
ypc = 0.0d+00
|
||||
uxs = 0.0d+00
|
||||
uys = 0.0d+00
|
||||
uxc = 0.0d+00
|
||||
uyc = 0.0d+00
|
||||
|
||||
else
|
||||
|
||||
! calculate initial positions and velocities
|
||||
!
|
||||
xps = mcomp / (mstar + mcomp) * dist + xshift
|
||||
xpc = - mstar / (mstar + mcomp) * dist + xshift
|
||||
uys = xps * om
|
||||
uyc = - xpc * om
|
||||
|
||||
! calculate mass fractions
|
||||
!
|
||||
ms = mcomp / (mstar + mcomp)
|
||||
mc = mstar / (mstar + mcomp)
|
||||
|
||||
! calculate orbit parameters
|
||||
!
|
||||
asemi = dist / (1.0d+00 - ecc)
|
||||
omega = pi2 / period
|
||||
ean = 0.0d+00
|
||||
|
||||
! set the initial positions
|
||||
!
|
||||
xrs = ms * dist + xshift
|
||||
yrs = 0.0d+00
|
||||
xrc = - mc * dist + xshift
|
||||
yrc = 0.0d+00
|
||||
|
||||
end if
|
||||
|
||||
as = 1.0d+00 / (vstar * dstar * rstar**2)
|
||||
ac = 1.0d+00 / (vcomp * dcomp * rcomp**2)
|
||||
|
||||
a(1) = as * xps**2 - ac * xpc**2
|
||||
a(2) = 2.0d+00 * (ac * xpc - as * xps)
|
||||
a(3) = as - ac
|
||||
|
||||
n = quadratic(a, x)
|
||||
if (a(3) >= 0.0d+00) then
|
||||
xshock = minval(x)
|
||||
else
|
||||
xshock = maxval(x)
|
||||
end if
|
||||
|
||||
! set the previous time
|
||||
!
|
||||
tprev = 0.0d+00
|
||||
|
||||
! proceed if no errors
|
||||
!
|
||||
if (status == 0) then
|
||||
|
||||
! print information about the user problem setup
|
||||
!
|
||||
call print_section(verbose, "Parameters (* - set, otherwise calculated)")
|
||||
call print_parameter(verbose, '(*) adiabatic_index' , adiabatic_index)
|
||||
call print_parameter(verbose, '(*) dratio = dstar/dcomp', dratio)
|
||||
call print_parameter(verbose, '(*) dstar' , dstar)
|
||||
call print_parameter(verbose, '( ) dcomp' , dcomp)
|
||||
call print_parameter(verbose, '(*) vstar' , vstar)
|
||||
call print_parameter(verbose, '(*) vcomp' , vcomp)
|
||||
call print_parameter(verbose, '(*) msstar' , msstar)
|
||||
call print_parameter(verbose, '(*) mscomp' , mscomp)
|
||||
call print_parameter(verbose, '( ) pstar' , pstar)
|
||||
call print_parameter(verbose, '( ) pcomp' , pcomp)
|
||||
call print_parameter(verbose, '(*) buni' , buni)
|
||||
call print_parameter(verbose, '(*) rstar' , rstar)
|
||||
call print_parameter(verbose, '(*) rcomp' , rcomp)
|
||||
call print_parameter(verbose, '(*) rbuffer' , rbuffer)
|
||||
call print_parameter(verbose, '(*) tstar' , tstar)
|
||||
call print_parameter(verbose, '(*) tcomp' , tcomp)
|
||||
if (orbiting) then
|
||||
call print_parameter(verbose, '(*) distance' , dist)
|
||||
call print_parameter(verbose, '(*) eccentricity', ecc)
|
||||
call print_parameter(verbose, '( ) acomp' , acomp)
|
||||
call print_parameter(verbose, '( ) bcomp' , bcomp)
|
||||
else
|
||||
call print_parameter(verbose, '(*) distance' , acomp)
|
||||
end if
|
||||
call print_parameter(verbose, '( ) shock position' , xshock)
|
||||
call print_parameter(verbose, '(*) shock curvature', cshock)
|
||||
|
||||
end if ! status
|
||||
|
||||
! associate subroutines to problem setup and shape update
|
||||
!
|
||||
setup_problem => setup_user_problem
|
||||
update_shapes => update_user_shapes
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine initialize_user_problem
|
||||
@ -91,7 +345,7 @@ module user_problem
|
||||
implicit none
|
||||
|
||||
integer, intent(out) :: status
|
||||
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
status = 0
|
||||
@ -105,7 +359,8 @@ module user_problem
|
||||
! subroutine SETUP_USER_PROBLEM:
|
||||
! -----------------------------
|
||||
!
|
||||
! Subroutine sets the initial conditions for the user specific problem.
|
||||
! Subroutine sets the initial conditions for the wind interaction in
|
||||
! the binary star problem.
|
||||
!
|
||||
! Arguments:
|
||||
!
|
||||
@ -116,14 +371,180 @@ module user_problem
|
||||
!
|
||||
subroutine setup_user_problem(pdata)
|
||||
|
||||
use blocks, only : block_data
|
||||
use blocks , only : block_data
|
||||
use coordinates, only : nn => bcells
|
||||
use coordinates, only : ax, ay, az
|
||||
use equations , only : prim2cons, magnetized
|
||||
use equations , only : nv
|
||||
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
|
||||
|
||||
implicit none
|
||||
|
||||
type(block_data), pointer, intent(inout) :: pdata
|
||||
|
||||
integer :: i, j, k
|
||||
real(kind=8) :: xs, ys, zs = 0.0d+00
|
||||
real(kind=8) :: xc, yc, zc = 0.0d+00
|
||||
real(kind=8) :: rs2, rc2, rs, rc, rd, ra
|
||||
real(kind=8) :: dns, prs, vxs, vys, vzs, bzs
|
||||
real(kind=8) :: dnc, prc, vxc, vyc, vzc, bzc
|
||||
|
||||
real(kind=8), dimension(nv,nn) :: q, u
|
||||
real(kind=8), dimension(nn) :: x, y
|
||||
#if NDIMS == 3
|
||||
real(kind=8), dimension(nn) :: z
|
||||
#endif /* NDIMS == 3 */
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#if NDIMS == 2
|
||||
k = 1
|
||||
#endif /* NDIMS == 2 */
|
||||
|
||||
x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
|
||||
y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
|
||||
#if NDIMS == 3
|
||||
z(:) = pdata%meta%zmin + az(pdata%meta%level,:)
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
! set magnetic field components
|
||||
!
|
||||
if (magnetized) then
|
||||
|
||||
q(ibx,:) = 0.0d+00
|
||||
q(iby,:) = 0.0d+00
|
||||
q(ibz,:) = buni
|
||||
q(ibp,:) = 0.0d+00
|
||||
|
||||
end if
|
||||
|
||||
#if NDIMS == 3
|
||||
! iterate over all positions in the YZ plane
|
||||
!
|
||||
do k = 1, nn
|
||||
|
||||
! calculate the Z coordinates of the central and companion stars
|
||||
!
|
||||
zs = z(k)
|
||||
zc = z(k)
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
do j = 1, nn
|
||||
|
||||
! calculate the Y coordinates of the central and companion stars
|
||||
!
|
||||
ys = y(j) - yps
|
||||
yc = y(j) - ypc
|
||||
|
||||
! distance from the line connecting both stars
|
||||
!
|
||||
#if NDIMS == 3
|
||||
ra = y(j)**2 + z(k)**2
|
||||
#else /* NDIMS == 3 */
|
||||
ra = y(j)**2
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
! sweep along the X coordinate
|
||||
!
|
||||
do i = 1, nn
|
||||
|
||||
! calculate the X coordinates of the central and companion stars
|
||||
!
|
||||
xs = x(i) - xps
|
||||
xc = x(i) - xpc
|
||||
|
||||
! calculate the distances from the centers of the central and companion stars
|
||||
!
|
||||
#if NDIMS == 3
|
||||
rs2 = max(xs * xs + ys * ys + zs * zs, 1.0d-16)
|
||||
rc2 = max(xc * xc + yc * yc + zc * zc, 1.0d-16)
|
||||
#else /* NDIMS == 3 */
|
||||
rs2 = max(xs * xs + ys * ys, 1.0d-16)
|
||||
rc2 = max(xc * xc + yc * yc, 1.0d-16)
|
||||
#endif /* NDIMS == 3 */
|
||||
rs = sqrt(rs2)
|
||||
rc = sqrt(rc2)
|
||||
|
||||
! calculate profiles from the central star
|
||||
!
|
||||
rd = max(1.0d+00, rs / rstar)
|
||||
|
||||
dns = dstar / rd**idx
|
||||
if (ipr > 0) prs = pstar / rd**ida
|
||||
if (ibz > 0) bzs = buni / rd**idm
|
||||
|
||||
! set the central star wind velocity
|
||||
!
|
||||
vxs = vstar * xs / rs
|
||||
vys = vstar * ys / rs
|
||||
vzs = vstar * zs / rs
|
||||
|
||||
! calculate profiles from the companion star
|
||||
!
|
||||
rd = max(1.0d+00, rc / rcomp)
|
||||
|
||||
dnc = dcomp / rd**idx
|
||||
if (ipr > 0) prc = pcomp / rd**ida
|
||||
if (ibz > 0) bzc = buni / rd**idm
|
||||
|
||||
! set the companion star wind velocity
|
||||
!
|
||||
vxc = vcomp * xc / rc
|
||||
vyc = vcomp * yc / rc
|
||||
vzc = vcomp * zc / rc
|
||||
|
||||
! set the variables
|
||||
!
|
||||
if (rs2 <= r2stari) then
|
||||
q(idn,i) = dns
|
||||
q(ivx,i) = vxs * rs / rstar - omstar * ys
|
||||
q(ivy,i) = vys * rs / rstar + omstar * xs + uys
|
||||
q(ivz,i) = vzs * rs / rstar
|
||||
if (ipr > 0) q(ipr,i) = prs
|
||||
if (ibz > 0) q(ibz,i) = bzs
|
||||
else if (rc2 <= r2compi) then
|
||||
q(idn,i) = dnc
|
||||
q(ivx,i) = vxc * rs / rcomp - omcomp * yc
|
||||
q(ivy,i) = vyc * rs / rcomp + omcomp * xc + uyc
|
||||
q(ivz,i) = vzc * rs / rcomp
|
||||
if (ipr > 0) q(ipr,i) = prc
|
||||
if (ibz > 0) q(ibz,i) = bzc
|
||||
else
|
||||
if (x(i) > xshock - cshock * ra) then
|
||||
q(idn,i) = dns
|
||||
q(ivx,i) = vxs
|
||||
q(ivy,i) = vys
|
||||
q(ivz,i) = vzs
|
||||
if (ipr > 0) q(ipr,i) = prs
|
||||
if (ibz > 0) q(ibz,i) = bzs
|
||||
else
|
||||
q(idn,i) = dnc
|
||||
q(ivx,i) = vxc
|
||||
q(ivy,i) = vyc
|
||||
q(ivz,i) = vzc
|
||||
if (ipr > 0) q(ipr,i) = prc
|
||||
if (ibz > 0) q(ibz,i) = bzc
|
||||
end if
|
||||
end if
|
||||
|
||||
end do ! i = 1, nn
|
||||
|
||||
! convert the primitive variables to conservative ones
|
||||
!
|
||||
call prim2cons(q(:,:), u(:,:))
|
||||
|
||||
! copy the primitive variables to the current block
|
||||
!
|
||||
pdata%q(:,:,j,k) = q(:,:)
|
||||
|
||||
! copy the conserved variables to the current block
|
||||
!
|
||||
pdata%u(:,:,j,k) = u(:,:)
|
||||
|
||||
end do ! j = 1, nn
|
||||
#if NDIMS == 3
|
||||
end do ! k = 1, nn
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -177,15 +598,194 @@ module user_problem
|
||||
!
|
||||
subroutine update_user_shapes(pdata, time, dt)
|
||||
|
||||
use blocks, only : block_data
|
||||
use blocks , only : block_data
|
||||
use coordinates, only : nn => bcells
|
||||
use coordinates, only : ax, ay, az
|
||||
use equations , only : nv
|
||||
use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
|
||||
use equations , only : prim2cons
|
||||
|
||||
implicit none
|
||||
|
||||
type(block_data), pointer, intent(inout) :: pdata
|
||||
real(kind=8) , intent(in) :: time, dt
|
||||
|
||||
integer :: i, j, k
|
||||
real(kind=8) :: xs, ys, zs
|
||||
real(kind=8) :: xc, yc, zc
|
||||
real(kind=8) :: dv, ds
|
||||
real(kind=8) :: sn, cs, man, res
|
||||
real(kind=8) :: rs2, rc2, rs, rc, rd
|
||||
|
||||
real(kind=8), dimension(nv) :: qs, qc
|
||||
real(kind=8), dimension(nn) :: x, y
|
||||
#if NDIMS == 3
|
||||
real(kind=8), dimension(nn) :: z
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! if time changes, we have to recalculate the star positions
|
||||
!
|
||||
if (time > tprev .and. orbiting) then
|
||||
|
||||
! solve the Kepler's equation to obtain the new value of eccentric anomaly
|
||||
!
|
||||
man = omega * time
|
||||
res = 1.0d+00
|
||||
i = 1
|
||||
do while(abs(res) >= tol .and. i <= maxit)
|
||||
res = (ean - ecc * sin(ean) - man) / (1.0d+00 - ecc * cos(ean))
|
||||
ean = ean - res
|
||||
i = i + 1
|
||||
end do
|
||||
if (abs(res) >= tol) then
|
||||
print *, "Kepler's equations could not be solved!"
|
||||
print *, "it: ", i, ", E: ", ean, ", dE: ", res
|
||||
print *, ""
|
||||
end if
|
||||
|
||||
! calculate trigonometric functions of the true anomaly
|
||||
!
|
||||
dv = 1.0d+00 - ecc * cos(ean)
|
||||
sn = sqrt(1.0d+00 - ecc**2) * sin(ean) / dv
|
||||
cs = (cos(ean) - ecc) / dv
|
||||
|
||||
! calculate the position and velocity of the companion star
|
||||
!
|
||||
ds = asemi * dv
|
||||
rs = ms * ds
|
||||
rc = mc * ds
|
||||
xps = rs * cs + xshift
|
||||
yps = rs * sn
|
||||
xpc = - rc * cs + xshift
|
||||
ypc = - rc * sn
|
||||
ds = time - tprev
|
||||
uxs = (xps - xrs) / ds
|
||||
uys = (yps - yrs) / ds
|
||||
uxc = (xpc - xrc) / ds
|
||||
uyc = (ypc - yrc) / ds
|
||||
|
||||
! update tprev, previous positions
|
||||
!
|
||||
tprev = time
|
||||
xrs = xps
|
||||
yrs = yps
|
||||
xrc = xpc
|
||||
yrc = ypc
|
||||
|
||||
end if ! time /= tprev
|
||||
|
||||
x(:) = pdata%meta%xmin + ax(pdata%meta%level,:)
|
||||
y(:) = pdata%meta%ymin + ay(pdata%meta%level,:)
|
||||
#if NDIMS == 3
|
||||
z(:) = pdata%meta%zmin + az(pdata%meta%level,:)
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
! set density and pressure in the stars' interiors
|
||||
!
|
||||
qs(idn) = dstar
|
||||
qc(idn) = dcomp
|
||||
if (ipr > 0) then
|
||||
qs(ipr) = pstar
|
||||
qc(ipr) = pcomp
|
||||
end if
|
||||
|
||||
! set magnetic field in the stars' interiors
|
||||
!
|
||||
if (ibx > 0) then
|
||||
qs(ibx) = 0.0d+00
|
||||
qs(iby) = 0.0d+00
|
||||
qs(ibz) = buni
|
||||
qs(ibp) = 0.0d+00
|
||||
qc(ibx) = 0.0d+00
|
||||
qc(iby) = 0.0d+00
|
||||
qc(ibz) = buni
|
||||
qc(ibp) = 0.0d+00
|
||||
end if
|
||||
|
||||
#if NDIMS == 3
|
||||
do k = 1, nn
|
||||
zs = z(k)
|
||||
zc = z(k)
|
||||
#else /* NDIMS == 3 */
|
||||
do k = 1, 1
|
||||
zs = 0.0d+00
|
||||
zc = 0.0d+00
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
do j = 1, nn
|
||||
ys = y(j) - yps
|
||||
yc = y(j) - ypc
|
||||
|
||||
do i = 1, nn
|
||||
xs = x(i) - xps
|
||||
xc = x(i) - xpc
|
||||
|
||||
! calculate distances from the centers of the central and companion stars
|
||||
!
|
||||
rs2 = xs * xs + ys * ys + zs * zs
|
||||
|
||||
if (rs2 <= r2staro) then
|
||||
|
||||
rs = max(sqrt(rs2), 1.0d-16)
|
||||
rd = max(1.0d+00, rs / rstar)
|
||||
|
||||
qs(idn) = dstar / rd**idx
|
||||
if (ipr > 0) qs(ipr) = pstar / rd**ida
|
||||
if (ibz > 0) qs(ibz) = buni / rd**idm
|
||||
|
||||
if (rs2 <= r2stari) then
|
||||
qs(ivx) = vstar * xs / rstar
|
||||
qs(ivy) = vstar * ys / rstar
|
||||
qs(ivz) = vstar * zs / rstar
|
||||
else
|
||||
qs(ivx) = vstar * xs / rs
|
||||
qs(ivy) = vstar * ys / rs
|
||||
qs(ivz) = vstar * zs / rs
|
||||
end if
|
||||
|
||||
qs(ivx) = qs(ivx) - omstar * ys + uxs
|
||||
qs(ivy) = qs(ivy) + omstar * xs + uys
|
||||
|
||||
pdata%q(:,i,j,k) = qs(:)
|
||||
|
||||
call prim2cons(qs(:), pdata%u(:,i,j,k))
|
||||
|
||||
end if
|
||||
|
||||
rc2 = xc * xc + yc * yc + zc * zc
|
||||
|
||||
if (rc2 <= r2compo) then
|
||||
|
||||
rc = max(sqrt(rc2), 1.0d-16)
|
||||
rd = max(1.0d+00, rc / rcomp)
|
||||
|
||||
qc(idn) = dcomp / rd**idx
|
||||
if (ipr > 0) qc(ipr) = pcomp / rd**ida
|
||||
if (ibz > 0) qc(ibz) = buni / rd**idm
|
||||
|
||||
if (rc2 <= r2compi) then
|
||||
qc(ivx) = vcomp * xc / rcomp
|
||||
qc(ivy) = vcomp * yc / rcomp
|
||||
qc(ivz) = vcomp * zc / rcomp
|
||||
else
|
||||
qc(ivx) = vcomp * xc / rc
|
||||
qc(ivy) = vcomp * yc / rc
|
||||
qc(ivz) = vcomp * zc / rc
|
||||
end if
|
||||
|
||||
qc(ivx) = qc(ivx) - omcomp * yc + uxc
|
||||
qc(ivy) = qc(ivy) + omcomp * xc + uyc
|
||||
|
||||
pdata%q(:,i,j,k) = qc(:)
|
||||
|
||||
call prim2cons(qc(:), pdata%u(:,i,j,k))
|
||||
end if
|
||||
end do
|
||||
|
||||
end do
|
||||
end do
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
|
Loading…
x
Reference in New Issue
Block a user