From e68a1899169e67d1abc0948d1bc26bcf87b891f0 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 17 May 2011 20:23:44 -0300 Subject: [PATCH 1/4] Implement Keplerian orbit for the satellite in binaries problem. - now the satellite changes positions and velocity accordind to Keplers laws; correct Keplerian orbit has been implemented in shape_binaries(); - correct initialization of the satellite in init_binaries(); --- src/problem.F90 | 180 ++++++++++++++++++++++++++++++++---------------- 1 file changed, 119 insertions(+), 61 deletions(-) diff --git a/src/problem.F90 b/src/problem.F90 index 70f4c10..ceb2987 100644 --- a/src/problem.F90 +++ b/src/problem.F90 @@ -709,7 +709,7 @@ module problem use config , only : im, jm, km use config , only : dens, pres, csnd2, gamma use config , only : dnfac, dnrat - use config , only : rstar, vstar, rsat, dsat, vsat, tsat + use config , only : rstar, vstar, rsat, dsat, vsat, tsat, esat use constants, only : dpi use coords , only : ax, ay, az, adr use scheme , only : prim2cons @@ -731,14 +731,14 @@ module problem ! local variables ! - integer :: i, j, k - real :: dr - real :: dnamb, pramb - real :: dnstar, prstar, vrstar, rc - real :: dnsat , prsat , vrsat , rs, xs, ys - real :: om = 0.0d0, sn = 0.0d0, cs = 1.0d0 - real :: asat, bsat, xsat - real :: xsh, ysh, xvl, yvl + integer :: i, j, k + real :: dr, rc, rs, xs, ys + real :: om = 0.0d0, sn = 0.0d0, cs = 1.0d0 + real, save :: dnamb, pramb + real, save :: dnstar, prstar, vrstar + real, save :: dnsat , prsat , vrsat + real, save :: asat, bsat, xsat = - 1.0d0 + real, save :: xsh, ysh, xvl, yvl ! local arrays ! @@ -749,32 +749,50 @@ module problem ! !------------------------------------------------------------------------------- ! +! calculate the interior and orbit parameters only once +! + if (xsat .le. 0.0d0) then + ! calculate pressure from sound speed ! #ifdef ISO - pres = csnd2 * dens + pres = csnd2 * dens #endif /* ISO */ #ifdef ADI - pres = csnd2 * dens / gamma + pres = csnd2 * dens / gamma #endif /* ADI */ -! calculate parameters +! calculate star interior parameters ! - dnamb = dens - pramb = pres - dnstar = dnamb * dnfac * dnrat - dnsat = dnamb * dnfac - prstar = pramb * dnfac - prsat = pramb * dnfac / dnrat - vrstar = vstar / rstar - vrsat = vsat / rsat - if (tsat .gt. 0.0d0) then - om = dpi / tsat - sn = 0.0d0 - cs = 1.0d0 + dnamb = dens + pramb = pres + dnstar = dnamb * dnfac * dnrat + dnsat = dnamb * dnfac + prstar = pramb * dnfac + prsat = pramb * dnfac / dnrat + vrstar = vstar / rstar + vrsat = vsat / rsat + +! calculate orbit parameters +! + asat = dsat / (1.0d0 - esat) + xsat = asat * esat + bsat = sqrt(asat * asat - xsat * xsat) + +! angular speed and triginometric functions +! + om = dpi / tsat + sn = 0.0d0 + cs = 1.0d0 + +! calculate satellite initial position and speed +! + xsh = asat * cs - xsat + ysh = bsat * sn + xvl = - asat * sn * om + yvl = bsat * cs * om + end if - xvl = - asat * om * sn - yvl = bsat * om * cs ! obtain the diagonal size of a cell at the current level ! @@ -1385,20 +1403,23 @@ module problem ! local variables ! - integer :: i, j, k - real :: dr - real :: dnamb, pramb - real :: dnstar, vrstar, prstar, rc - real :: dnsat , vrsat , prsat , rs, xs, ys + integer :: i, j, k + real :: dr, rs, xs, ys, rc + real :: om = 0.0d0, sn = 0.0d0, cs = 1.0d0 + real :: rhs, res + real, save :: dnamb, pramb + real, save :: dnstar, vrstar, prstar + real, save :: dnsat , vrsat , prsat #ifdef ADI - real :: ekin, ekstar, eksat + real, save :: ekstar, eksat + real :: ekin #ifdef MHD - real :: emag + real :: emag #endif /* MHD */ #endif /* ADI */ - real :: om = 0.0d0, ang = 0.0d0, sn = 0.0d0, cs = 1.0d0 - real :: asat, bsat, xsat - real :: xsh, ysh, xvl, yvl + real, save :: asat, bsat, xsat = -1.0d0 + real, save :: xsh, ysh, xvl, yvl + real, save :: tpre = -1.0d0, ang = 0.0d0 ! local arrays ! @@ -1408,39 +1429,76 @@ module problem ! !------------------------------------------------------------------------------- ! -! calculate parameters +! calculate the star and satellite orbit parameters only once +! + if (xsat .lt. 0.0d0) then + +! calculate parameters for the star interiors ! #ifdef ISO - pres = csnd2 * dens + pres = csnd2 * dens #endif /* ISO */ #ifdef ADI - pres = csnd2 * dens / gamma + pres = csnd2 * dens / gamma #endif /* ADI */ - dnamb = dens - pramb = pres - dnstar = dnamb * dnfac * dnrat - dnsat = dnamb * dnfac - prstar = pramb * dnfac - prsat = pramb * dnfac / dnrat - vrstar = vstar / rstar - vrsat = vsat / rsat + dnamb = dens + pramb = pres + dnstar = dnamb * dnfac * dnrat + dnsat = dnamb * dnfac + prstar = pramb * dnfac + prsat = pramb * dnfac / dnrat + vrstar = vstar / rstar + vrsat = vsat / rsat #ifdef ADI - ekstar = 0.5d0 * dnstar * vstar * vstar - eksat = 0.5d0 * dnsat * vsat * vsat + ekstar = 0.5d0 * dnstar * vstar * vstar + eksat = 0.5d0 * dnsat * vsat * vsat #endif /* ADI */ - if (tsat .gt. 0.0d0) then - om = dpi / tsat - ang = om * t - sn = sin(ang) - cs = cos(ang) + +! calculate orbit parameters +! + asat = dsat / (1.0d0 - esat) + xsat = asat * esat + bsat = sqrt(asat * asat - xsat * xsat) + + end if + +! recalculate the angles and coordinates +! + if (t .gt. tpre) then + +! calculate the new value of eccentric anomaly +! + rhs = dpi * t / tsat + res = 1.0d0 + i = 1 + do while(res .ge. 1.0e-14 .and. i .le. 1000) + sn = esat * sin(ang) + ang = rhs + sn + res = abs((ang - sn) - rhs) + i = i + 1 + end do + +! calculate trigonometric functions +! + sn = sin(ang) + cs = cos(ang) + +! calculate the angular velocity +! + om = dpi / tsat / (1.0d0 - esat * cs) + +! calculate position coordinates and velocity components +! + xsh = asat * cs - xsat + ysh = bsat * sn + xvl = - asat * sn * om + yvl = bsat * cs * om + +! substitute the previous time +! + tpre = t + end if - asat = dsat / (1.0d0 - esat) - xsat = asat * esat - bsat = sqrt(asat * asat - xsat * xsat) - xsh = asat * cs - xsat - ysh = bsat * sn - xvl = - asat * om * sn - yvl = bsat * om * cs ! obtain the diagonal size of a cell at the current level ! From 3e586c94d336749136296688eec3af21686e0319 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 17 May 2011 20:26:33 -0300 Subject: [PATCH 2/4] Fix finding of dxmin in find_new_timestep(). --- src/evolution.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/evolution.F90 b/src/evolution.F90 index bcd9922..88e1bfb 100644 --- a/src/evolution.F90 +++ b/src/evolution.F90 @@ -249,17 +249,17 @@ module evolution end do ! meta blocks + end if ! maxlev > 1 + ! find the smallest spacial step ! #if NDIMS == 2 - dxmin = min(adx(lev), ady(lev)) + dxmin = min(adx(lev), ady(lev)) #endif /* NDIMS == 2 */ #if NDIMS == 3 - dxmin = min(adx(lev), ady(lev), adz(lev)) + dxmin = min(adx(lev), ady(lev), adz(lev)) #endif /* NDIMS == 3 */ - end if ! maxlev > 1 - ! iterate over all data blocks in order to find the maximum speed among them ! pdata => list_data From acbce0c1258023b256786b3f0f6da2214d966567 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 17 May 2011 20:29:41 -0300 Subject: [PATCH 3/4] Fix initial angular velocity in init_binaries(). --- src/problem.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/problem.F90 b/src/problem.F90 index ceb2987..3d55479 100644 --- a/src/problem.F90 +++ b/src/problem.F90 @@ -781,7 +781,7 @@ module problem ! angular speed and triginometric functions ! - om = dpi / tsat + om = dpi / tsat / (1.0d0 - esat) sn = 0.0d0 cs = 1.0d0 From 6c68231c6aef2bd870e9ea15ff249ca1323a102b Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 17 May 2011 21:01:20 -0300 Subject: [PATCH 4/4] Reduce number of local variables in binaries problem. --- src/problem.F90 | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/src/problem.F90 b/src/problem.F90 index 3d55479..d13c4e2 100644 --- a/src/problem.F90 +++ b/src/problem.F90 @@ -737,7 +737,7 @@ module problem real, save :: dnamb, pramb real, save :: dnstar, prstar, vrstar real, save :: dnsat , prsat , vrsat - real, save :: asat, bsat, xsat = - 1.0d0 + real, save :: asat, bsat = - 1.0d0 real, save :: xsh, ysh, xvl, yvl ! local arrays @@ -751,7 +751,7 @@ module problem ! ! calculate the interior and orbit parameters only once ! - if (xsat .le. 0.0d0) then + if (bsat .le. 0.0d0) then ! calculate pressure from sound speed ! @@ -776,18 +776,17 @@ module problem ! calculate orbit parameters ! asat = dsat / (1.0d0 - esat) - xsat = asat * esat - bsat = sqrt(asat * asat - xsat * xsat) + bsat = asat * sqrt(1.0d0 - esat * esat) ! angular speed and triginometric functions ! - om = dpi / tsat / (1.0d0 - esat) sn = 0.0d0 cs = 1.0d0 + om = dpi / tsat / (1.0d0 - esat * cs) ! calculate satellite initial position and speed ! - xsh = asat * cs - xsat + xsh = asat * (cs - esat) ysh = bsat * sn xvl = - asat * sn * om yvl = bsat * cs * om @@ -1417,7 +1416,7 @@ module problem real :: emag #endif /* MHD */ #endif /* ADI */ - real, save :: asat, bsat, xsat = -1.0d0 + real, save :: asat, bsat = -1.0d0 real, save :: xsh, ysh, xvl, yvl real, save :: tpre = -1.0d0, ang = 0.0d0 @@ -1431,7 +1430,7 @@ module problem ! ! calculate the star and satellite orbit parameters only once ! - if (xsat .lt. 0.0d0) then + if (bsat .lt. 0.0d0) then ! calculate parameters for the star interiors ! @@ -1457,8 +1456,7 @@ module problem ! calculate orbit parameters ! asat = dsat / (1.0d0 - esat) - xsat = asat * esat - bsat = sqrt(asat * asat - xsat * xsat) + bsat = asat * sqrt(1.0d0 - esat * esat) end if @@ -1471,7 +1469,7 @@ module problem rhs = dpi * t / tsat res = 1.0d0 i = 1 - do while(res .ge. 1.0e-14 .and. i .le. 1000) + do while(res .ge. 1.0e-15 .and. i .le. 100) sn = esat * sin(ang) ang = rhs + sn res = abs((ang - sn) - rhs) @@ -1489,7 +1487,7 @@ module problem ! calculate position coordinates and velocity components ! - xsh = asat * cs - xsat + xsh = asat * (cs - esat) ysh = bsat * sn xvl = - asat * sn * om yvl = bsat * cs * om