PLUTO  4.4-patch2
Functions
rotate.c File Reference

Rotate a 1D problem. More...

#include "pluto.h"
#include "rotate.h"

Functions

double RotateGet_ta (void)
 
double RotateGet_tb (void)
 
void RotateSet (int jx, int jy, int kx, int kz)
 
void RotateVector (double *v, int s)
 

Detailed Description

Given the initial condition in the 1D, non-rotated frame ( $\Sigma'$), we rotate quantities by an angle $\gamma $ around the y-axis and an angle $\alpha$ around the z-axis. The resulting transformation leaves scalar quantities invariant and produce vector rotations $\vec{q} = \tens{R}_{\gamma\alpha}\vec{q}'$ where

\[ \tens{R}_{\gamma\alpha} = \left(\begin{array}{ccc} \cos\alpha\cos\gamma & -\sin\alpha & -\cos\alpha\sin\gamma \\ \noalign{\medskip} \sin\alpha\cos\gamma & \cos\alpha & -\sin\alpha\sin\gamma \\ \noalign{\medskip} \sin\gamma & 0 & \cos\gamma \end{array}\right) \,,\qquad \tens{R}^{-1}_{\gamma\alpha} = \left(\begin{array}{ccc} \cos\alpha\cos\gamma & \sin\alpha\cos\gamma & \sin\gamma \\ \noalign{\medskip} -\sin\alpha & \cos\alpha & 0 \\ \noalign{\medskip} -\cos\alpha\sin\gamma & -\sin\alpha\sin\gamma & \cos\gamma \end{array}\right) \]

where $\tan\gamma = \cos\alpha\tan\beta$. In 2D, simply set $\gamma = 0$.

In a discrete domain, the angles $\alpha$ and $\beta$ must be chosen in such a way that an integer shift of cells satisfies, for any flow quantity q, the translational invariance expressed by $ q(\vec{x}+\vec{s}) = q (\vec{x})$ where $\vec{s} = (n_x\Delta x, n_y \Delta y, n_z \Delta z)$ is orthogonal to $\vec{n}_p$:

\[ \vec{n}_p\cdot\vec{s} = n_x\Delta x + n_y\Delta y\tan\alpha + n_z\Delta z\tan\beta = 0 \]

The mesh spacing does need to be the same in all directions. For a given pair of rotation angles, there can be more than one set of integers that satisfies this equation. For instance, taking $\tan\alpha = 1,\;\tan\beta=2 $ one may choose (if mesh spacing is the same) $\vec{s} = \pm(1,1,-1) $, $\vec{s} = \pm (-1,1,0)$ (useful at a y-boundary) or $\vec{s} = \pm (2,0,-1) $ (useful at a z-boundary). In practice, we find convenient to provide two pairs of integer numbers, $(j_x,\, j_y)$ and $ (k_x,\, k_z)$ giving the integer shifts along two directions at a time. These numbers are used to find the rotation angles:

\[ \left\{\begin{array}{lcl} \DS j_x\Delta x + j_y\Delta y \tan\alpha &=& 0 \\ \noalign{\medskip} \DS k_x\Delta x + k_z\Delta z \tan\beta &=& 0 \,,\qquad \end{array}\right. \qquad\rightarrow\qquad \tan\alpha = -\frac{j_x}{j_y}\frac{\Delta x}{\Delta y} ,\quad \tan\beta = -\frac{k_x}{k_z}\frac{\Delta x}{\Delta z} \]

This choice allows us to assign boundary conditions using only 2 shifts (rather than 3): one along the x direction and the other in the direction normal to the boundary plane.

Plane of Discontinuity. A plane of discontinuity has equation $x' = 0$ in the non-rotated frame. Using the inverse transformation we find the equation of the plane in the rotated $\Sigma$-plane:

\[ x' = \left(\tens{R}^{-1}_{\gamma} \vec{x}\right)_x = (\cos\alpha\cos\gamma)\, x + (\sin\alpha\sin\gamma)\, y + (\sin\gamma) z = 0 \quad\rightarrow\quad \boxed{x + y\tan\alpha + z\tan\beta = 0} \]

The normal to this plane is the vector $\vec{n}_p = (1, \tan\alpha, \tan\beta)$.

Note on periodic domains. In a periodic domain, an integer number of wavelengths must be contained in each direction. Hence the rotation is usuallt applied by specifying $\tan\alpha = k_y/k_x$ and $\tan\beta = k_z/k_x$ which express the ratios between the $y$- and $z$- components of the wavevector with the $x$ component. For one wavelength in each direction, $k_x = 2\pi/L_x$, $k_y = 2\pi/L_y$, $k_z = 2\pi/L_z$. The integer shifts can then be specified using

\[ \tan\alpha = -\frac{j_x}{j_y}\frac{\Delta x}{\Delta y} = \frac{L_x}{L_y} , \qquad \tan\beta = -\frac{k_x}{k_z}\frac{\Delta x}{\Delta z} = \frac{L_x}{L_z} \]

Note on the final time. Usually, the rotated domain size differs from the typical 1D configuration since solutions are compared by looking at the profiles in the x direction. The final time is then adjusted in such a way that the intersection of the rotated plane with the x axis travels the same distance covered in the 1D case. Imagine a point in 1D that travels a distance $\delta L_1 = v\delta t_1$: in 2- or 3-D the same point will be traveling along the direction given by $\DS \vec{x}_p(t) = \hvec{n}_p vt$ where $ \hvec{n}_p = \vec{n}_p/|\vec{n}_p| $. In the rotated problem the solution is constant on the plane of equation $ \hvec{n}_p\cdot(\vec{x}-\vec{x}_p(t)) = 0$ and its intersection with the x axis is ( $ y = z = 0 $) covers the distance $\delta L_1$

\[ \hat{n}_x\delta L_1 = v\delta t \qquad\rightarrow\qquad \delta t = n_x\delta t_1 = \frac{\delta t_1}{\sqrt{1 + \tan^2\alpha + \tan^2\beta}} \]

Author
A. Mignone (migno.nosp@m.ne@t.nosp@m.o.inf.nosp@m.n.it)
Date
Apr 05, 2020

Reference:

Function Documentation

◆ RotateGet_ta()

double RotateGet_ta ( void  )

Return the rotation angle tan(alpha)

◆ RotateGet_tb()

double RotateGet_tb ( void  )

Return the rotation angle tan(beta)

◆ RotateSet()

void RotateSet ( int  jx,
int  jy,
int  kx,
int  kz 
)

For a given choice of the input shifts, initialize the rotation angles rot_ta $ = \tan\alpha$ and rot_tb $ = \tan\beta $.

◆ RotateVector()

void RotateVector ( double *  v,
int  s 
)

Rotate a vector <v[0], v[1], v[2]>;

  • s = 1 –> clockwise rotation [transform vector to the 1D unrotated frame]
  • s = -1 –> counter-clockwise rotation [transform vector to the 2/3D frame]