All sections and subsections can be shown/hidden by clicking on the header.
Two stars orbit about some point $\mathbf{R}_\text{cm}$ with angular
velocity $\mathbf{\Omega} = \Omega\ \hat{z}$. For a frame co-rotating
about one of the stars rather than about $\mathbf{R}_\text{cm}$, a gas
element in this frame will experience a force $\mathbf{f}_r$ of
$$
\mathbf{f}_r
= \mathbf{g}
- 2 \mathbf{\Omega} \times \mathbf{u}
- \mathbf{\Omega} \times \mathbf{\Omega} \times (\mathbf{r}-\mathbf{R}_\text{cm})
$$
The primary star is placed at grid center with the secondary at $\mathbf{r}_S = R_\text{sep}\ \hat{x}$
Since the code runs in spherical coordinates, $\mathbf{f}_r$ is best calculated by
converting $\mathbf{g}$ and $\mathbf{\Omega}$ to spherical coordiantes. This is done
using the following transformation:
$$
\left( \begin{array}{c}
\hat{r} \\
\hat{\theta} \\
\hat{\phi}
\end{array} \right)
=
\left( \begin{array}{ccc}
\sin\theta\cos\phi & \sin\theta\sin\phi & \cos\theta \\
\cos\theta\cos\phi & \cos\theta\sin\phi & -\sin\theta \\
-\sin\phi & \cos\phi & 0
\end{array} \right)
\left( \begin{array}{c}
\hat{x} \\
\hat{y} \\
\hat{z}
\end{array} \right)
$$
1.1 Converting $\mathbf{g}$ to spherical
expand «
Since the grid is centered about the primary, this force $\mathbf{g}_P$ is simply
$$
\mathbf{g}_P =
\left( \begin{array}{c}
\kappa(r)\ \frac{GM_P}{r^2} \\ 0 \\ 0
\end{array} \right)
$$
For the companion,
$$
\mathbf{g}_S = - \frac{GM_S}{\left|\mathbf{r}-\mathbf{R}_\text{sep}\right|^3}
\left(\mathbf{r}-\mathbf{R}_\text{sep}\right)
$$
The stars are placed along the $x$-axis so that the separation between the two stars
$\mathbf{R}_\text{sep} = R_\text{sep}\ \hat{x}$ on the Yin grid and
$-R_\text{sep}\ \hat{x}$ on the Yang.
In spherical, these become
$$
\mathbf{R}_\text{sep} =
\pm R_\text{sep}
\left( \begin{array}{c}
\sin\theta \cos\phi \\
\cos\theta \cos\phi \\
-\sin\phi
\end{array} \right)
$$
The total gravity due to the companion is then
$$
\mathbf{g}_S = -
\frac{GM_S}{\left|\mathbf{r}-\mathbf{R}_\text{sep}\right|^3}
\left( \begin{array}{c}
r \mp R_\text{sep} \sin\theta\cos\phi \\
\mp R_\text{sep} \cos\theta\cos\phi \\
\pm R_\text{sep} \sin\phi
\end{array} \right)
$$
1.2 Converting $\mathbf{\Omega}$ to spherical
expand «
For $\mathbf{\Omega} = \Omega\ \hat{z}$, we have
$$
\left( \begin{array}{ccc}
\sin\theta\cos\phi & \sin\theta\sin\phi & \cos\theta \\
\cos\theta\cos\phi & \cos\theta\sin\phi & -\sin\theta \\
-\sin\phi & \cos\phi & 0
\end{array} \right)
\left( \begin{array}{c}
0 \\ 0 \\ \Omega
\end{array} \right)
= \Omega
\left( \begin{array}{c}
\cos\theta \\
-\sin\theta \\
0
\end{array} \right),
$$
and for $\mathbf{\Omega} = \Omega\ \hat{y}$ (Yang grid),
$$
\left( \begin{array}{ccc}
\sin\theta\cos\phi & \sin\theta\sin\phi & \cos\theta \\
\cos\theta\cos\phi & \cos\theta\sin\phi & -\sin\theta \\
-\sin\phi & \cos\phi & 0
\end{array} \right)
\left( \begin{array}{c}
0 \\ \Omega \\ 0
\end{array} \right)
= \Omega
\left( \begin{array}{c}
\sin\theta \sin\phi \\
\cos\theta \sin\phi \\
\cos\phi
\end{array} \right)
$$
1.3 Centrifugal Force
expand «
For the Yin grid
$$
\mathbf{\Omega}\times\mathbf{r} = r\Omega
\left| \begin{array}{ccc}
\hat{r} & \hat{\theta} & \hat{\phi} \\
\cos\theta & -\sin\theta & 0 \\
1 & 0 & 0
\end{array} \right| = r\Omega
\left( \begin{array}{c}
0 \\ 0 \\ \sin\theta
\end{array} \right)
$$
$$
\mathbf{\Omega}\times\mathbf{\Omega}\times\mathbf{r} = r\Omega^2
\left| \begin{array}{ccc}
\hat{r} & \hat{\theta} & \hat{\phi} \\
\cos\theta & -\sin\theta & 0 \\
0 & 0 & \sin\theta
\end{array} \right| = r\Omega^2
\left( \begin{array}{c}
-\sin^2\theta \\
-\sin\theta\cos\theta \\
0
\end{array} \right)
$$
And for the Yang
$$
\mathbf{\Omega}\times\mathbf{r} = r\Omega
\left| \begin{array}{ccc}
\hat{r} & \hat{\theta} & \hat{\phi} \\
\sin\theta \sin\phi & \cos\theta \sin\phi & \cos\phi \\
1 & 0 & 0
\end{array} \right| = r\Omega
\left( \begin{array}{c}
0 \\ \cos\phi \\ -\cos\theta \sin\phi
\end{array} \right)
$$
$$
\mathbf{\Omega}\times\mathbf{\Omega}\times\mathbf{r} = r\Omega^2
\left| \begin{array}{ccc}
\hat{r} & \hat{\theta} & \hat{\phi} \\
\sin\theta \sin\phi & \cos\theta \sin\phi & \cos\phi \\
0 & \cos\phi & -\cos\theta \sin\phi
\end{array} \right| = r\Omega^2
\left( \begin{array}{c}
-\cos^2\theta \sin^2\phi - \cos^2\phi \\
\sin\theta \cos\theta \sin^2\phi \\
\sin\theta \sin\phi \cos\phi
\end{array} \right)
$$
1.4 Translational Force
expand «
For the Yin grid
$$
\mathbf{\Omega}\times\mathbf{R}_\text{cm} = R_\text{cm}\Omega
\left| \begin{array}{ccc}
\hat{r} & \hat{\theta} & \hat{\phi} \\
\cos\theta & -\sin\theta & 0 \\
\sin\theta \cos\phi & \cos\theta \cos\phi & -\sin\phi
\end{array} \right| = R_\text{cm}\Omega
\left( \begin{array}{c}
\sin\theta \sin\phi \\
\cos\theta \sin\phi \\
\cos\phi
\end{array} \right)
$$
$$
\mathbf{\Omega}\times\mathbf{\Omega}\times\mathbf{R}_\text{cm}
= R_\text{cm}\Omega^2
\left| \begin{array}{ccc}
\hat{r} & \hat{\theta} & \hat{\phi} \\
\cos\theta & -\sin\theta & 0 \\
\sin\theta \sin\phi & \cos\theta \sin\phi & \cos\phi
\end{array} \right| = R_\text{cm}\Omega^2
\left( \begin{array}{c}
-\sin\theta \cos\phi \\
-\cos\theta \cos\phi \\
\sin\phi
\end{array} \right)
$$
and for the Yang grid,
$$
\mathbf{\Omega}\times\mathbf{\Omega}\times\mathbf{R}_\text{cm}
= -R_\text{cm}\Omega^2
\left( \begin{array}{c}
-\sin\theta \cos\phi \\
-\cos\theta \cos\phi \\
\sin\phi
\end{array} \right)
$$
1.5 Coriolis Force
expand «
Calculating $\Omega \times \mathbf{u}$ on the Yin grid
$$
\Omega\times\mathbf{u} = \Omega
\left| \begin{array}{ccc}
\hat{r} & \hat{\theta} & \hat{\phi} \\
\cos\theta & -\sin\theta & 0 \\
u & v & w
\end{array} \right| = \Omega
\left( \begin{array}{c}
-w\sin\theta \\
-w\cos\theta \\
v\cos\theta + u\sin\theta
\end{array} \right)
$$
and for the Yang grid,
$$
\Omega\times\mathbf{u} = \Omega
\left| \begin{array}{ccc}
\hat{r} & \hat{\theta} & \hat{\phi} \\
\sin\theta\sin\phi & \cos\theta\sin\phi & \cos\phi \\
u & v & w
\end{array} \right| = \Omega
\left( \begin{array}{c}
w\cos\theta\sin\phi - v\cos\phi \\
u\cos\phi - w\sin\theta\sin\phi \\
v\sin\theta\sin\phi - u\cos\theta\sin\phi
\end{array} \right)
$$
1.6 Total Forces
contract «
$\mathbf{g}$ on the Yin and Yang grids.
note: $\left|\mathbf{r}-\mathbf{R}_\text{sep}\right|^3$
is still correct on the Yang grid as $\mathbf{R}_\text{sep}$ changes
$$
^n \mathbf{g} =
\frac{GM_P}{r^2}
\left( \begin{array}{c}
\kappa(r) \\ 0 \\ 0
\end{array} \right)
- \frac{GM_S}{\left|\mathbf{r}-\mathbf{R}_\text{sep}\right|^3}
\left( \begin{array}{c}
r - R_\text{sep} \sin\theta \cos\phi \\
- R_\text{sep} \cos\theta \cos\phi \\
R_\text{sep} \sin\phi
\end{array} \right)
$$
$$
^g \mathbf{g} =
\frac{GM_P}{r^2}
\left( \begin{array}{c}
\kappa(r) \\ 0 \\ 0
\end{array} \right)
- \frac{GM_S}{\left|\mathbf{r}-\mathbf{R}_\text{sep}\right|^3}
\left( \begin{array}{c}
r + R_\text{sep} \sin\theta \cos\phi \\
R_\text{sep} \cos\theta \cos\phi \\
- R_\text{sep} \sin\phi
\end{array} \right)
$$
Pseudoforce $\mathbf{a}$ on the Yin and Yang grids
$$
^n \mathbf{a} = r\Omega^2
\left( \begin{array}{c}
\sin^2\theta \\
\sin\theta \cos\theta \\
0
\end{array} \right)
- R_\text{cm}\Omega^2
\left( \begin{array}{c}
\sin\theta \cos\phi \\
\cos\theta \cos\phi \\
-\sin\phi
\end{array} \right)
- 2 \Omega
\left( \begin{array}{c}
-w \sin\theta \\
-w \cos\theta \\
v \cos\theta + u \sin\theta
\end{array} \right)
$$
$$
^g \mathbf{a} = r\Omega^2
\left( \begin{array}{c}
\cos^2\theta \sin^2\phi + \cos^2\phi \\
- \sin\theta \cos\theta \sin^2\phi \\
- \sin\theta \sin\phi \cos\phi
\end{array} \right)
+ R_\text{cm}\Omega^2
\left( \begin{array}{c}
\sin\theta \cos\phi \\
\cos\theta \cos\phi \\
-\sin\phi
\end{array} \right)
- 2 \Omega
\left( \begin{array}{c}
w \cos\theta \sin\phi - v \cos\phi \\
u \cos\phi - w \sin\theta \sin\phi \\
v \sin\theta \sin\phi - u \cos\theta \sin\phi
\end{array} \right)
$$
1.7 x-Directional Sweep
expand «
Forces in the x-direction are
$$
\hat{i} \cdot\ ^n \mathbf{g}
= \kappa(r)\ \frac{GM_P}{r^2}
- \frac{GM_S}{\left|\mathbf{r}-\mathbf{R}_\text{sep}\right|^3}
(r - R_\text{sep}\sin\theta\cos\phi)
$$
$$
\hat{i} \cdot\ ^n \mathbf{a}
= r\Omega^2 \sin^2\theta
- R_\text{cm}\Omega^2 \sin\theta \cos\phi
+ 2 \Omega w \sin\theta
$$
and in code,
do n = nmin-4, nmax+5
dSx = xf(n)*stheta*cphi - sep
dSy = xf(n)*stheta*sphi
dSz = xf(n)*ctheta
radS = max(sqrt(dSx**2+dSy**2+dSz**2), minRad)
grav(n) = + GMP / xf(n)**2 * kap(n) &
- GMS / radS**3 * (xf(n)-sep*stheta*cphi)
fict(n) = + (w(n)*w(n)+v(n)*v(n))/xf(n) & ! grid force
+ xf(n)*omega**2*stheta**2 & ! centrifugal
- rcm*omega**2*stheta*cphi & ! translation
+ 2.0*omega*w(n)*stheta ! coriolis
end do
For the Yang grid,
$$
\hat{i} \cdot\ ^g \mathbf{g}
= \kappa(r)\ \frac{GM_P}{r^2}
- \frac{GM_S}{\left|\mathbf{r}-\mathbf{R}_\text{sep}\right|^3}
(r + R_\text{sep}\sin\theta\cos\phi)
$$
$$
\hat{i} \cdot\ ^g \mathbf{a}
= r\Omega^2(\cos^2\theta \sin^2\phi + \cos^2\phi)
+ R_\text{cm}\Omega^2\sin\theta\cos\phi
+ 2 \Omega (v \cos\phi - w \cos\theta \sin\phi)
$$
do n = nmin-4, nmax+5
dSx = xf(n)*stheta*cphi + sep
dSy = xf(n)*stheta*sphi
dSz = xf(n)*ctheta
radS = max(sqrt(dSx**2+dSy**2+dSz**2), minRad)
grav(n) = + GMP / xf(n)**2 * kap(n) &
- GMS / radS**3 * (xf(n)+sep*stheta*cphi)
fict(n) = + (w(n)*w(n)+v(n)*v(n))/xf(n) & ! grid force
+ xf(n)*omega**2*(cphi**2+(ctheta*sphi)**2) & ! centrifugal
+ rcm*omega**2*stheta*cphi & ! translation
+ 2.0*omega*(v(n)*cphi-w(n)*ctheta*sphi) ! coriolis
end do
1.8 y-Directional Sweep
expand «
Forces in the y-direction are
$$
\hat{j} \cdot\ ^n \mathbf{g}
= \frac{GM_S}{\left|\mathbf{r}-\mathbf{R}_\text{sep}\right|^3}
R_\text{sep} \cos\theta \cos\phi
$$
$$
\hat{j} \cdot\ ^n \mathbf{a}
= r\Omega^2 \sin\theta \cos\theta
- R_\text{cm}\Omega^2 \cos\theta \cos\phi
+ 2 \Omega w \cos\theta
$$
note for j sweep, $u \to w, v \to u, w \to v$
do n = nmin-4, nmax+5
dSx = radius*sin(xf(n))*cphi - sep
dSy = radius*sin(xf(n))*sphi
dSz = radius*cos(xf(n))
radS = max(sqrt(dSx**2+dSy**2+dSz**2), minRad)
grav(n) = + GMS / radS**3 * sep*cos(xf(n))*cphi
fict(n) = + v(n)*v(n)*cos(xf(n))/(radius*sin(xf(n))) & ! grid force
- u(n)*w(n)/radius & ! grid force
+ radius*omega**2*sin(xf(n))*cos(xf(n)) & ! centrifugal
- rcm*omega**2*cos(xf(n))*cphi & ! translation
+ 2.0*omega*v(n)*cos(xf(n)) ! coriolis
end do
For the Yang grid,
$$
\hat{j} \cdot\ ^g \mathbf{g}
= - \frac{GM_S}{\left|\mathbf{r}-\mathbf{R}_\text{sep}\right|^3}
R_\text{sep} \cos\theta \cos\phi
$$
$$
\hat{j} \cdot\ ^g \mathbf{a}
= - r\Omega^2 \sin\theta \cos\theta \sin^2\phi
+ R_\text{cm}\Omega^2 \cos\theta \cos\phi
- 2 \Omega (u\cos\phi - w\sin\theta\sin\phi)
$$
do n = nmin-4, nmax+5
dSx = radius*sin(xf(n))*cphi + sep
dSy = radius*sin(xf(n))*sphi
dSz = radius*cos(xf(n))
radS = max(sqrt(dSx**2+dSy**2+dSz**2), minRad)
grav(n) = + GMS / radS**3 * sep*cos(xf(n))*cphi
fict(n) = + v(n)*v(n)*cos(xf(n))/(radius*sin(xf(n))) & ! grid force
- u(n)*w(n)/radius & ! grid force
- radius*omega**2*sin(xf(n))*cos(xf(n))*sphi**2 & ! centrifugal
+ rcm*omega**2*cos(xf(n))*cphi & ! translation
+ 2.0*omega*(v(n)*sin(xf(n))*sphi-w(n)*cphi) ! coriolis
end do
1.9 z-Directional Sweep
expand «
Forces in the z-direction are
$$
\hat{k} \cdot\ ^n \mathbf{g}
= - \frac{GM_S}{\left|\mathbf{r}-\mathbf{R}_\text{sep}\right|^3}
R_\text{sep} \sin\phi
$$
$$
\hat{k} \cdot\ ^n \mathbf{a}
= R_\text{cm}\Omega^2 \sin\phi
- 2 \Omega (v\cos\theta + u\sin\theta)
$$
note for k sweep, $u \to v, v \to w, w \to u$
do n = nmin-4, nmax+5
dSx = radius*sin(xf(n))*cphi - sep
dSy = radius*sin(xf(n))*sphi
dSz = radius*cos(xf(n))
radS = max(sqrt(dSx**2+dSy**2+dSz**2), minRad)
grav(n) = - GMS / radS**3 * sep*sin(xf(n))
fict(n) = - u(n)*w(n)/radius*ctheta & ! grid force
- u(n)*v(n)/radius*stheta & ! grid force
+ 0.0 & ! centrifugal
+ rcm*omega**2*sin(xf(n)) & ! translation
- 2.0*omega*(w(n)*ctheta+v(n)*stheta) ! coriolis
end do
For the Yang grid,
$$
\hat{k} \cdot\ ^g \mathbf{g}
= \frac{GM_S}{\left|\mathbf{r}-\mathbf{R}_\text{sep}\right|^3}
R_\text{sep} \sin\phi
$$
$$
\hat{k} \cdot\ ^g \mathbf{a}
= - r\Omega^2 \sin\theta \sin\phi \cos\phi
- R_\text{cm}\Omega^2 \sin\phi
- 2 \Omega (v\sin\theta\sin\phi - u\cos\theta\sin\phi)
$$
do n = nmin-4, nmax+5
dSx = radius*sin(xf(n))*cphi + sep
dSy = radius*sin(xf(n))*sphi
dSz = radius*cos(xf(n))
radS = max(sqrt(dSx**2+dSy**2+dSz**2), minRad)
grav(n) = + GMS / radS**3 * sep*sin(xf(n))
fict(n) = - u(n)*w(n)/radius*ctheta & ! grid force
- u(n)*v(n)/radius*stheta & ! grid force
- (radius/stheta)*omega**2*stheta*sin(xf(n))*cos(xf(n)) & !centrifugal
- rcm*omega**2*sin(xf(n)) & ! translation
+ 2.0*omega*(v(n)*ctheta*sin(xf(n))-w(n)*stheta*sin(xf(n)))
end do
The boundary code has 4 steps.
- Transform $\mathbf{u}$ to global cartesian velocity in an inertial frame,
- Find location on inner grid for each boundary coordinate,
- Interpolate all data values for each coordinate,
- Transform $\mathbf{u}$ back to local spherical velocity.
The term 'global' means independent of grid. Velocities are handled in global cartesian for
ease of interpolation and for ease of converting back to spherical when considering the
time-dependent rotational vector.
2.1 Velocity Transformation
expand «
The outer grid will be cauclated in a frame that is translating but not rotating with the primary star.
Velocity is transformed during boundary value generation to avoid repetitive calculation during run-time. The
Grid Rotation effect is removed using the relation
$\mathbf{u}' = \mathbf{u} + \mathbf{\Omega} \times \mathbf{r}$.
The Grid Translation effect may be removed using the relation
$\mathbf{u}' = \mathbf{u} - \mathbf{\Omega} \times \mathbf{R}_\text{cm}$,
but we leave the term in for simplicty (further explained in 2.4).
Velocity $\mathbf{u}_t$ in the translating frame is then
$$\mathbf{u}_t = \mathbf{u}_r + \mathbf{\Omega} \times \mathbf{r}$$
Data is initially loaded into arrays nData and gData, which hold $(\rho, P, u, v, w)$
information for the Yin and Yang grids, respectively. Coordinate arrays are given a _n or _g
suffix to designate which grid they are associated with.
For the Yin grid,
do k = 1, kmax_n
sph = sin(zzc_n(k))
cph = cos(zzc_n(k))
do j = 1, jmax_n
sth = sin(zyc_n(j))
cth = cos(zyc_n(j))
do i = 1, imax_n
! global cartesian coordinates
rGlobal(1) = zxc_n(i)*sth*cph
rGlobal(2) = zxc_n(i)*sth*sph
rGlobal(3) = zxc_n(i)*cth
! local cartesian velocity : load
vLocal(1) = nData(3,i,j,k)*sth*cph + nData(4,i,j,k)*cth*cph - nData(5,i,j,k)*sph
vLocal(2) = nData(3,i,j,k)*sth*sph + nData(4,i,j,k)*cth*sph + nData(5,i,j,k)*cph
vLocal(3) = nData(3,i,j,k)*cth - nData(4,i,j,k)*sth
! local cartesian velocity -> global cartesian velocity
vGlobal(1) = vLocal(1)
vGlobal(2) = vLocal(2)
vGlobal(3) = vLocal(3)
! global cartesian velocity : remove rotation
vGlobal(1) = vGlobal(1) - omega*rGlobal(2)
vGlobal(2) = vGlobal(2) + omega*rGlobal(1)
vGlobal(3) = vGlobal(3)
! global cartesian velocity : save
nData(3,i,j,k) = vGlobal(1)
nData(4,i,j,k) = vGlobal(2)
nData(5,i,j,k) = vGlobal(3)
end do
end do
end do
and similarly for the Yang grid,
do k = 1, kmax_g
sph = sin(zzc_g(k))
cph = cos(zzc_g(k))
do j = 1, jmax_g
sth = sin(zyc_g(j))
cth = cos(zyc_g(j))
do i = 1, imax_g ! same as imax_n
! global cartesian coordinates
rGlobal(1) =-zxc_n(i)*sth*cph
rGlobal(2) = zxc_n(i)*cth
rGlobal(3) = zxc_n(i)*sth*sph
! local cartesian velocity : load
vLocal(1) = gData(3,i,j,k)*sth*cph + gData(4,i,j,k)*cth*cph - gData(5,i,j,k)*sph
vLocal(2) = gData(3,i,j,k)*sth*sph + gData(4,i,j,k)*cth*sph + gData(5,i,j,k)*cph
vLocal(3) = gData(3,i,j,k)*cth - gData(4,i,j,k)*sth
! local cartesian velocity -> global cartesian velocity
vGlobal(1) =-vLocal(1)
vGlobal(2) = vLocal(3)
vGlobal(3) = vLocal(2)
! global cartesian velocity : remove rotation
vGlobal(1) = vGlobal(1) - omega*rGlobal(2)
vGlobal(2) = vGlobal(2) + omega*rGlobal(1)
vGlobal(3) = vGlobal(3)
! global cartesian velocity : save
gData(3,i,j,k) = vGlobal(1)
gData(4,i,j,k) = vGlobal(2)
gData(5,i,j,k) = vGlobal(3)
end do
end do
end do
2.2 Determining Boundary Location
expand «
For each angular coordinate (on each grid), six radii are required to form the boundary at each time-step.
Since the inner grid is effectively steady-state, it is simplest to just save 128 'rotation' zones for each boundary
zone, and determine which rotation zones to pull from during runtime.
At any time $t$, the inner grid has rotated by $\Omega t$. To create boundary data for that time at any point
$\mathbf{r}_o$, we first translate to the inner grid (shift the $x$-coordinate by $R_\text{cm}$), and then rotate
in the $\textit{opposite}$ direction by $\Omega t$ (since that point will now be at $\mathbf{r}_o$). Note for the
'translational grid' method, the grid is not shifted since each grid is centered about the primary.
$$
\left(\begin{array}{c}
x_i \\ y_i \\ z_i
\end{array} \right) =
\left(\begin{array}{c}
\cos\left(-\Omega t\right) & - \sin\left(-\Omega t\right) & 0 \\
\sin\left(-\Omega t\right) & \cos\left(-\Omega t\right) & 0 \\
0 & 0 & 1
\end{array} \right)
\left[
\left(\begin{array}{c}
x_o \\ y_o \\ z_o
\end{array} \right) -
\left(\begin{array}{c}
R_\text{cm} \\ 0 \\ 0
\end{array} \right)
\right]
$$
There are various ways to optimize this, but I will present the most detailed. The array rad stores the
six boundary radii of the outer grid, and the array phi stores the 128 'rotation' angles.
do k = 1, kmax_n
sph = sin(zzc_n(k))
cph = cos(zzc_n(k))
do j = 1, jmax_n
sth = sin(zyc_n(j))
cth = cos(zyc_n(j))
do i = 1, nzones_rad
! LOCAL coords (cartesian)
rLocal(1) = rad(i)*sth*cph
rLocal(2) = rad(i)*sth*sph
rLocal(3) = rad(i)*cth
! GLOBAL coords (cartesian) with CM shift
rGlobal(1) = rLocal(1) - rcm
rGlobal(2) = rLocal(2)
rGlobal(3) = rLocal(3)
! rotate about the z-axis
do n = 1, nzones_phi
cart(1) = rGlobal(1)*cos(phi(n)) + rGlobal(2)*sin(phi(n))
cart(2) =-rGlobal(1)*sin(phi(n)) + rGlobal(2)*cos(phi(n))
cart(3) = rGlobal(3)
cart(4) = sqrt(cart(1)**2 + cart(2)**2)
! find GLOBAL spherical coordinates
nrad = sqrt(cart(1)**2 + cart(2)**2 + cart(3)**2)
nthe = atan2(cart(4),cart(3))
nphi = atan2(cart(2),cart(1))
end do
end do
end do
end do
Once $^n \theta$ and $^n \phi$ are determined, we can determine the grid this location is on
simply by checking $^n \theta$.
yin = .true.
if (nthe < pi/4. .or. nthe > 3*pi/4.) yin = .false.
if (yin) then
[...determine interpolation points based on YIN grid...]
else
[...determine interpolation points based on YANG grid...]
end if
2.3 Interpolating Boundary Values
expand «
Values for the boundary zones are linearly interpolated using the eight zones closest to the point found above.
A subscript/suffix b is used to designate a point before the $i$-, $j$-, and $k$-coordinate and
a is used to designate the point after. Similarly, a subscript/suffix 1 is used to
designate the weight of the point before, and 2 is used for the point after. e.g., for
interpolation along the radial,
$$
\rho = i_1\rho(i_b, j, k) + i_2\rho(i_a, j, k),
$$
where
$$
i_1 = (r_a - r)/(r_a-r_b)
$$
$$
i_2 = (r - r_b)/(r_a-r_b)
$$
Some adjustments must be made in code if a point goes off its respective grid or encounters a periodic boundary,
but the basic algorithm is as follows
iA = 2
do while (zxc_n(iA) < rad(i))
iA = iA + 1
end do
iB = iA - 1
i1 = (zxc_n(iA)-rad(i))/(zxc_n(iA)-zxc_n(iB))
i2 = (rad(i)-zxc_n(iB))/(zxc_n(iA)-zxc_n(iB))
dthe = zyc_n(2)-zyc_n(1)
jB = 1 + (nthe-zyc_n(1))/dthe
jA = 1 + jB
j1 = (zyc_n(jA)-nthe)/dthe
j2 = (nthe-zyc_n(jB))/dthe
dphi = zzc_n(2)-zzc_n(1)
kB = 1 + (nphi-zzc_n(1))/dphi
kA = 1 + kB
k1 = (zzc_n(kA)-nphi)/dphi
k2 = (nphi-zzc_n(kB))/dphi
do nv = 1, nvars
c1 = i1*nData(nv,iB,jB,kB) + i2*nData(nv,iA,jB,kB)
c2 = i1*nData(nv,iB,jB,kA) + i2*nData(nv,iA,jB,kA)
c3 = i1*nData(nv,iB,jA,kB) + i2*nData(nv,iA,jA,kB)
c4 = i1*nData(nv,iB,jA,kA) + i2*nData(nv,iA,jA,kA)
nBound(nv,n,i,j,k) = j1*(k1*c1 + k2*c2) + j2*(k1*c3 + k2*c4)
end do
To correctly interpolate points very near the edge, the following method is used
- Which grid is the rotated boundary zone on?
- Yin - Find 8 closest zones based on Yin coordinate arrays (e.g.
zyc_n)
- Yang - Find 8 closest zones based on Yin coordinate arrays
Find global coordinates of these 8 closest points
- For each, determine which grid that point is on
- If a point is on the same grid as its rotated boundary zone,
determine its c-value for each variable.
! c-value for c1 point (jB, kB)
do nv = 1, nvars
c1(nv) = i1*nData(nv,iB,jB,kB) + i2*nData(nv,iA,jB,kB)
end do
- If a point is on a different grid as its rotated boundary zone,
interpolate the point onto opposing grid using a separate set of 8 zones, then
determine its c-value for each variable.
! c-value for c3 point (jA, kB) - here, zyc_n(jA) is off-grid
athe = zyc_n(jB) + (zyc_n(2)-zyc_n(1))
aphi = zzc_n(kB)
bthe = acos (sin(athe) *sin(aphi)) ! angular coords on opposing grid
bphi = atan2(cos(athe),-sin(athe)*cos(aphi))
dthe = zyc_g(2)-zyc_g(1)
nB = 1 + (bthe-zyc_g(1))/dthe
nA = 1 + nB
n1 = (zyc_g(nA)-bthe)/dthe
n2 = (bthe-zyc_g(nB))/dthe
dphi = zzc_g(2)-zzc_g(1)
mB = 1 + (bphi-zzc_g(1))/dphi
mA = 1 + mB
m1 = (zzc_g(mA)-bphi)/dphi
m2 = (bphi-zzc_g(mB))/dphi
do nv = 1, nvars
b1 = i1*gData(nv,iB,nB,mB) + i2*gData(nv,iA,nB,mB)
b2 = i1*gData(nv,iB,nB,mA) + i2*gData(nv,iA,nB,mA)
b3 = i1*gData(nv,iB,nA,mB) + i2*gData(nv,iA,nA,mB)
b4 = i1*gData(nv,iB,nA,mA) + i2*gData(nv,iA,nA,mA)
c3(nv) = n1*(m1*b1 + m2*b2) + n2*(m1*b3 + m2*b4)
end do
- Once all c-values are collected, final values are interpolated
do nv = 1, nvars
nBound(nv,n,i,j,k) = j1*(k1*c1(nv) + k2*c2(nv)) + j2*(k1*c3(nv) + k2*c4(nv))
end do
2.4 Velocity Transformation (again)
expand «
The final step is to rotate $\mathbf{u}$ and then transform it to local spherical velocity.
Since $\mathbf{u}$ was pulled from a point that has since been rotated by $\Omega t$, $\mathbf{u}$
is rotated as follows
$$
\left( \begin{array}{c}
u_f \\ v_f \\ w_f
\end{array} \right) =
\left( \begin{array}{c}
\cos\left(\Omega t\right) & - \sin\left(\Omega t\right) & 0 \\
\sin\left(\Omega t\right) & \cos\left(\Omega t\right) & 0 \\
0 & 0 & 1
\end{array} \right)
\left( \begin{array}{c}
u \\ v \\ w
\end{array} \right)
$$
For the 'translating grid' method, we must also account for velocity due to the grid's translation. Fortunately,
the contributing term $\mathbf{\Omega} \times \mathbf{R}_\text{cm}$ rotates about $\mathbf{\Omega}$ as well,
allowing us to ignore it entirely (otherwise we would need to subtract it out in section 2.1, then rotate it by
$\Omega t$ and add it back in here).
In code,
! YIN grid
do k = 1, kmax_n
sph = sin(zzc_n(k))
cph = cos(zzc_n(k))
do j = 1, jmax_n
sth = sin(zyc_n(j))
cth = cos(zyc_n(j))
do i = 1, nzones_rad
do n = 1, nzones_phi
! global cartesian velocity : read
vGlobal(1) = nBound(3,n,i,j,k)
vGlobal(2) = nBound(4,n,i,j,k)
vGlobal(3) = nBound(5,n,i,j,k)
! global cartesian velocity : rotate about z-axis
rGlobal(1) = vGlobal(1)*cos(phi(n)) - vGlobal(2)*sin(phi(n))
rGlobal(2) = vGlobal(1)*sin(phi(n)) + vGlobal(2)*cos(phi(n))
rGlobal(3) = vGlobal(3)
! global cartesian velocity -> local cartesian velocity
vLocal(1) = rGlobal(1)
vLocal(2) = rGlobal(2)
vLocal(3) = rGlobal(3)
! local cartesian velocity -> local spherical velocity
nBound(3,n,i,j,k) = vLocal(1)*sth*cph + vLocal(2)*sth*sph + vLocal(3)*cth
nBound(4,n,i,j,k) = vLocal(1)*cth*cph + vLocal(2)*cth*sph - vLocal(3)*sth
nBound(5,n,i,j,k) =-vLocal(1)* sph + vLocal(2)* cph
end do
end do
end do
end do
! YANG grid
do k = 1, kmax_g
sph = sin(zzc_g(k))
cph = cos(zzc_g(k))
do j = 1, jmax_g
sth = sin(zyc_g(j))
cth = cos(zyc_g(j))
do i = 1, nzones_rad
do n = 1, nzones_phi
! global cartesian velocity : read
vGlobal(1) = gBound(3,n,i,j,k)
vGlobal(2) = gBound(4,n,i,j,k)
vGlobal(3) = gBound(5,n,i,j,k)
! global cartesian velocity : rotate about z-axis
rGlobal(1) = vGlobal(1)*cos(phi(n)) - vGlobal(2)*sin(phi(n))
rGlobal(2) = vGlobal(1)*sin(phi(n)) + vGlobal(2)*cos(phi(n))
rGlobal(3) = vGlobal(3)
! global cartesian velocity -> local cartesian velocity
vLocal(1) =-rGlobal(1)
vLocal(2) = rGlobal(3)
vLocal(3) = rGlobal(2)
! local cartesian velocity -> local spherical velocity
gBound(3,n,i,j,k) = vLocal(1)*sth*cph + vLocal(2)*sth*sph + vLocal(3)*cth
gBound(4,n,i,j,k) = vLocal(1)*cth*cph + vLocal(2)*cth*sph - vLocal(3)*sth
gBound(5,n,i,j,k) =-vLocal(1)* sph + vLocal(2)* cph
end do
end do
end do
end do
The inner grid is run to a quasi-steady state so that time-dependence of the outer grid's inner boundary may be simplified by
using a snapshot of the inner grid data. For each angular zone (j,k), the code loads a 'ring' of data ranging $0$ to $2\pi$,
each point of which contains six ghost zones.
The location $\alpha$ on this ring is determined by the simulation time such that $\alpha = \Omega t$, and all boundary values
are generated via interpolation of the two nearest zones.
alpha = omega * time
do while (alpha > 2.0*pi)
alpha = alpha - 2.0*pi
end do
nBelow = 1 + (alpha - angle(1))/(angle(2)-angle(1))
nAbove = 1 + nBelow
! adjust for periodic boundary
if (nBelow == 0) then
a1 = angle(nAbove)
a2 = a1 - (angle(2)-angle(1))
nBelow = amax
else if (nAbove == amax+1) then
a2 = angle(nBelow)
a1 = a2 + (angle(2)-angle(1))
nAbove = 1
else
a2 = angle(nBelow)
a1 = angle(nAbove)
end if
! interpolation coefficients
n1 = (a1-alpha)/(angle(2)-angle(1))
n2 = (alpha-a2)/(angle(2)-angle(1))
Once interpolation coefficients are calculated, the code simply loops through each angular zone and populates the boundary zones
as follows. (Here, rin(i,j,k) is the $i$-th ghost zone for the ($j$,$k$) coordinate.)
do k = 1, ks
do j = 1, js
if (ongrid(jcol*js+j,krow*ks+k) < 2) then
do nv = 1, nvars
val = n2*iBound(nv,nAbove,i,j,k) + n1*iBound(nv,nBelow,i,j,k)
if (nv == 1) rin(i,j,k) = val
if (nv == 2) pin(i,j,k) = val
if (nv == 3) uin(i,j,k) = val
if (nv == 4) vin(i,j,k) = val
if (nv == 5) win(i,j,k) = val
end do
end if
end do
end do
