Appendix / Natural Convection
|
|
Rama
physics Appendix / Natural Convection |
This section explains how we implement the numerical fluid simulation method presented in the Natural convection and wind speed section. The explanations are interspersed with the corresponding source code, using a kind of literate programming style. The source code can be divided in several sections:
We start with the definition of the physical and Rama constants that we will need:
#include <cmath> #include <iostream> #include <fstream> // ----------------------------------------------------------------------------- // Physical and Rama constants // ----------------------------------------------------------------------------- // Stefan-Bolztmann constant, in W.m^-2.K^-4. const double sigma = 5.67e-8; // Specific gas constant of air, in J.kg^-1.K^-1. const double Rs = 8.314 / 0.0289; // Heat capacity at constant volume of air, in J.kg^-1.K^-1. const double cv = 716.8; // Radius of Rama's surface, in m. const double rs = 8000.0; // Angular speed of Rama in rad.s^-1. const double omega = 2.0 * M_PI / 180.0; // Hydrostatic equilibrium temperature, in K. const double T0 = 288.0; // Hyrostatic equilibrium pressure at the ground level, in Pa. const double ps = 101325.0; // Emissivity coefficient of the atmosphere at the ground level, in m^-1. const double kE = 2.58e-4; // Emissivity density coefficient of the atmosphere, in m^2.kg^-1. const double kappaE = kE / (ps / (Rs * T0));
To solve this issue, we introduce a halfint class, which represents a half integer. By definition, halfint(i) represents i+0.5:
class halfint {
public:
explicit halfint(int v) : v(v) {}
int operator()() const { return v; }
private:
int v;
};
We can then define operators to add and subtract integers and half integers:
halfint operator+(halfint i, int j) { return halfint(i() + j); }
halfint operator+(int i, halfint j) { return halfint(i + j()); }
int operator+(halfint i, halfint j) { return i() + j() + 1; }
halfint operator-(halfint i, int j) { return halfint(i() - j); }
halfint operator-(int i, halfint j) { return halfint(i - j() - 1); }
int operator-(halfint i, halfint j) { return i() - j(); }
Finally, we introduce a "half" constant holding the value 1/2, that we will need
very often:
const halfint half(0);We can now define some classes to manipulate staggered grids. First, we define the resolution $N_r \times N_{\theta}$ we want to use for these grids, and the corresponding grid cell sizes $\Delta_r=r_s/N_r$ and $\Delta_{\theta}=2\pi/3N_{\theta}$:
const int Nr = 64; const int Nt = 128; const double DELTAr = rs / Nr; const double DELTAt = 2.0 * M_PI / 3.0 / Nt;We are now ready to define the staggered grid classes. The goal here is to have one class per type of grid (radial, tangential, or centered, e.g. for the radial velocity, the tangential velocity, and the temperature or pressure). We also want arithmetic operators to combine grids together, as well as interpolation and partial differentiation operators. In order to implement this in a way that is efficient at runtime, we use expression templates.
We start by defining the main template from which all grid expression templates will inherit. This template provides a single operator to access the elements of the grid. This operator is itself a template, taking the types of the grid indices as arguments (these types will be either int or halfint, depending on the kind of grid).
template<class E> class GridExpr {
public:
operator E&() { return static_cast<E&>(*this); }
operator E const&() const { return static_cast<const E&>(*this); }
template<typename I, typename J> double operator()(I i, J j) const {
return static_cast<E const&>(*this)(i, j);
}
};
We can now define arithmetic operators to combine grids together. For each
arithmetic operation, here +, we first define a grid expression class inheriting
from our abstract grid expression template, such that an access to an element of
this grid returns the addition of the same elements in two other grids...
template<class E1, class E2>
class GridAdd : public GridExpr< GridAdd<E1, E2> > {
public:
GridAdd(const GridExpr<E1>& g, const GridExpr<E2>& h) : g(g), h(h) {}
template<typename I, typename J> double operator()(I i, J j) const {
return g(i, j) + h(i, j);
}
private:
const E1& g;
const E2& h;
};
... and we then overload the corresponding C++ operator so that we can use it on
grid expressions too (allowing us to write a+b, where a and b are staggered
grids, as a shorthand for GridAdd(a, b)):
template<class E1, class E2>
const GridAdd<E1,E2> operator+(const GridExpr<E1>& g, const GridExpr<E2>& h) {
return GridAdd<E1,E2>(g, h);
}
We repeat this pattern with the subtraction, multiplication and division
operators:
template<class E1, class E2>
class GridSub : public GridExpr< GridSub<E1, E2> > {
public:
GridSub(const GridExpr<E1>& g, const GridExpr<E2>& h) : g(g), h(h) {}
template<typename I, typename J> double operator()(I i, J j) const {
return g(i, j) - h(i, j);
}
private:
const E1& g;
const E2& h;
};
template<class E1, class E2>
const GridSub<E1,E2> operator-(const GridExpr<E1>& g, const GridExpr<E2>& h) {
return GridSub<E1,E2>(g, h);
}
template<class E>
class GridScale : public GridExpr< GridScale<E> > {
public:
GridScale(double s, const GridExpr<E>& g) : s(s), g(g) {}
template<typename I, typename J> double operator()(I i, J j) const {
return s * g(i, j);
}
private:
const double s;
const E& g;
};
template<class E> const GridScale<E> operator*(double s, const GridExpr<E>& g) {
return GridScale<E>(s, g);
}
template<class E> const GridScale<E> operator*(const GridExpr<E>& g, double s) {
return GridScale<E>(s, g);
}
template<class E1, class E2>
class GridMul : public GridExpr< GridMul<E1, E2> > {
public:
GridMul(const GridExpr<E1>& g, const GridExpr<E2>& h) : g(g), h(h) {}
template<typename I, typename J> double operator()(I i, J j) const {
return g(i, j) * h(i, j);
}
private:
const E1& g;
const E2& h;
};
template<class E1, class E2>
const GridMul<E1,E2> operator*(const GridExpr<E1>& g, const GridExpr<E2>& h) {
return GridMul<E1,E2>(g, h);
}
template<class E1, class E2>
class GridDiv : public GridExpr< GridDiv<E1, E2> > {
public:
GridDiv(const GridExpr<E1>& g, const GridExpr<E2>& h) : g(g), h(h) {}
template<typename I, typename J> double operator()(I i, J j) const {
return g(i, j) / h(i, j);
}
private:
const E1& g;
const E2& h;
};
template<class E1, class E2>
const GridDiv<E1,E2> operator/(const GridExpr<E1>& g, const GridExpr<E2>& h) {
return GridDiv<E1,E2>(g, h);
}
The averaging operators $\overline{f}^r$ and $\overline{f}^{\theta}$ can be
defined in a very similar way, except that we don't have native C++ operators
for them. We thus use C++ functions instead, named avgR and avgT. Note the use
of half integers here, which gives a source code that is very similar to the
equations defining there operators:
$(\overline{f}^r)_{i,\,j}=(f_{i+\h,\,j}+f_{i-\h,\,j})/2$ and
$(\overline{f}^{\theta})_{i,\,j}=f_{i,\,j+\h}+f_{i,\,j-\h})/2$.
template<class E> class GridAvgR : public GridExpr< GridAvgR<E> > {
public:
GridAvgR(const GridExpr<E>& g) : g(g) {}
template<typename I, typename J> double operator()(I i, J j) const {
return (g(i + half, j) + g(i - half, j)) * 0.5;
}
private:
const E& g;
};
template<class E> const GridAvgR<E> avgR(const GridExpr<E>& g) {
return GridAvgR<E>(g);
}
template<class E> class GridAvgT : public GridExpr< GridAvgT<E> > {
public:
GridAvgT(const GridExpr<E>& g) : g(g) {}
template<typename I, typename J> double operator()(I i, J j) const {
return (g(i, j + half) + g(i, j - half)) * 0.5;
}
private:
const E& g;
};
template<class E> const GridAvgT<E> avgT(const GridExpr<E>& g) {
return GridAvgT<E>(g);
}
The partial differentiation operators $\partial_r f$ and $\partial_{\theta} f$
are implemented in a similar, straightforward way from their mathematical
definition, in functions named partialR and partialT:
template<class E> class GridPartialR : public GridExpr< GridPartialR<E> > {
public:
GridPartialR(const GridExpr<E>& g) : g(g) {}
template<typename I, typename J> double operator()(I i, J j) const {
return (g(i + half, j) - g(i - half, j)) * (1.0 / DELTAr);
}
private:
const E& g;
};
template<class E> const GridPartialR<E> partialR(const GridExpr<E>& g) {
return GridPartialR<E>(g);
}
template<class E> class GridPartialT : public GridExpr< GridPartialT<E> > {
public:
GridPartialT(const GridExpr<E>& g) : g(g) {}
template<typename I, typename J> double operator()(I i, J j) const {
return (g(i, j + half) - g(i, j - half)) * (1.0 / DELTAt);
}
private:
const E& g;
};
template<class E> const GridPartialT<E> partialT(const GridExpr<E>& g) {
return GridPartialT<E>(g);
}
We now have all the operators we need to combine staggered grids, but we don't
have any concrete class to instantiate a staggered grid. Let's first define a
class for axisymmetric grids, i.e. grids whose elements with the same radial
index all have the same value. These grids are defined for any kind of indices,
so we provide element access operators for each combination of index types. We
also provide convenient element access operators taking only a radial index as
argument. The element values are stored in two 1D arrays, one for the integer
indices, the other for the half integer indices. Note that we declare $N_r+3$
integer indices for $N_r$ cells, because the averaging and partial
differentiation operators may require one external value on each side of the
grid, which itself contains $N_r+1$ values.
class AxisymmetricGrid : public GridExpr<AxisymmetricGrid> {
public:
double operator()(int i, int j) const { return fInt[i]; }
double operator()(int i, halfint j) const { return fInt[i]; }
double operator()(halfint i, int j) const { return fHalfInt[i()]; }
double operator()(halfint i, halfint j) const { return fHalfInt[i()]; }
double& operator()(int i) { return fInt[i]; }
double& operator()(halfint i) { return fHalfInt[i()]; }
private:
double fInt[Nr + 3];
double fHalfInt[Nr + 2];
};
We continue by defining an abstract parent class for staggered grids,
encapsulating a 2D array of values initialized to 0. Note that we declare
NR+2$\times$NT+2 elements for a NR$\times$NT grid, because the averaging and
partial differentiation operators may require one external value on each side of
the grid. We also provide a setBoundaryConditions method to initialize these
outside values, using periodic boundary conditions for $\theta$ and
Neumann
(using sign=1) or Dirichlet (using
sign=-1) boundary conditions for $r$ (using sign=0 is also possible, to set
these values to 0).
template<class E, int NR, int NT> class StaggeredGrid : public GridExpr<E> {
public:
StaggeredGrid () {
for (int i = 0; i < NR + 2; ++i) {
for (int j = 0; j < NT + 2; ++j) {
values[i][j] = 0.0;
}
}
}
void setBoundaryConditions(double sign) {
for (int j = 1; j <= NT; ++j) {
values[0][j] = sign * values[1 + (NR % 2)][j];
values[NR + 1][j] = sign * values[NR - (NR % 2)][j];
}
for (int i = 0; i < NR + 2; ++i) {
values[i][0] = values[i][NT];
values[i][NT + 1] = values[i][1];
}
}
protected:
double values[NR + 2][NT + 2];
};
Finally, we define concrete classes for radial, tangential and centered grids.
Each class provides an element access operator, and only one, of the type
expected by the kind of the grid. The goal is to get compilation errors when
trying to combine grids in a way that does not make sense (e.g. trying to assign
avgR(ur) to a radial grid, assuming ur is radial grid - in which case avgR(ur)
is a centered grid).
Each class also provides a template assignment operator, so that we can store the result of any grid expression into a new grid, e.g. f = g + h, where f, g and h are grids. Note that these operators only evaluate the grid elements: they do not evaluate nor assign values to the external rows and columns on each side of the grid - this is done separately, using the inherited setBoundaryConditions method.
class RadialGrid : public StaggeredGrid<RadialGrid, Nr + 1, Nt> {
public:
template<class E> RadialGrid& operator=(const GridExpr<E>& grid) {
// Exclude the boundary itself (r=0 and r=rs, i.e. i=1 and i=Nr+1), this is
// not needed since the velocity is always 0 on this boundary.
for (int i = 2; i <= Nr; ++i) {
for (int j = 1; j <= Nt; ++j) {
values[i][j] = grid(i, j + half);
}
}
return *this;
}
double operator()(int i, halfint j) const { return values[i][j()]; }
double& operator()(int i, halfint j) { return values[i][j()]; }
};
class TangentialGrid : public StaggeredGrid<TangentialGrid, Nr, Nt> {
public:
template<class E> TangentialGrid& operator=(const GridExpr<E>& grid) {
for (int i = 1; i <= Nr; ++i) {
for (int j = 1; j <= Nt; ++j) {
values[i][j] = grid(i + half, j);
}
}
return *this;
}
double operator()(halfint i, int j) const { return values[i()][j]; }
double& operator()(halfint i, int j) { return values[i()][j]; }
};
class CenteredGrid : public StaggeredGrid<CenteredGrid, Nr, Nt> {
public:
template<class E> CenteredGrid& operator=(const GridExpr<E>& grid) {
for (int i = 1; i <= Nr; ++i) {
for (int j = 1; j <= Nt; ++j) {
values[i][j] = grid(i + half, j + half);
}
}
return *this;
}
double operator()(halfint i, halfint j) const { return values[i()][j()]; }
double& operator()(halfint i, halfint j) { return values[i()][j()]; }
};
A tridiagonal linear system can be solved in linear time. Here we simply adapt the solve_tridiagonal_in_place_reusable C++ implementation of the Thomas algorithm, using complex numbers instead of reals for the x vector. TridiagonalSolver encapsulates a tridiagonal matrix, which can be initialized with the init() method (the ai, bi, ci represent the lower, middle and upper diagonals, respectively). Once initialized, an object of this class can be used any number of time to solve many tridiagonal systems having the same matrix but not the same input vector.
class TridiagonalSolver {
public:
void init(int i, double ai, double bi, double ci) {
a[i] = ai;
b[i] = bi;
c[i] = ci;
}
void solve(double* x) {
d[0] = c[0] / b[0];
x[0] = x[0] / b[0];
x[1] = x[1] / b[0];
for (int i = 1; i < Nr; i++) {
double m = 1.0 / (b[i] - a[i] * d[i - 1]);
d[i] = c[i] * m;
x[2 * i] = (x[2 * i] - a[i] * x[2 * (i - 1)]) * m;
x[2 * i + 1] = (x[2 * i + 1] - a[i] * x[2 * (i - 1) + 1]) * m;
}
for (int i = Nr - 2; i >= 0; i--) {
x[2 * i] = x[2 * i] - d[i] * x[2 * (i + 1)];
x[2 * i + 1] = x[2 * i + 1] - d[i] * x[2 * (i + 1) + 1];
}
}
private:
double a[Nr];
double b[Nr];
double c[Nr];
static double d[Nr];
};
double TridiagonalSolver::d[Nr];
For the fast Fourier transform we simply reuse the source code from
"A Simple and Efficient FFT Implementation in C++" Vlodymyr Myrnyy, 2007,
which is based on templates:
template<int N, int S> class DanielsonLanczos {
public:
void apply(double* data) {
next.apply(data);
next.apply(data + N);
const double wpr = -2.0 * sin(M_PI / N) * sin(M_PI / N);
const double wpi = S * sin(2.0 * M_PI / N);
double wr = 1.0;
double wi = 0.0;
for (unsigned i = 0; i < N; i += 2) {
double tempr = data[i + N] * wr - data[i + N + 1] * wi;
double tempi = data[i + N] * wi + data[i + N + 1] * wr;
data[i + N] = data[i] - tempr;
data[i + N + 1] = data[i+1] - tempi;
data[i] += tempr;
data[i + 1] += tempi;
double wtemp = wr;
wr += wr * wpr - wi * wpi;
wi += wi * wpr + wtemp * wpi;
}
}
private:
DanielsonLanczos<N / 2, S> next;
};
template<int S> class DanielsonLanczos<1, S> {
public:
void apply(double* data) { }
};
template<int N> class FFT {
public:
void fft(double* data) { scramble(data); recursion.apply(data); }
void ifft(double* data) { scramble(data); irecursion.apply(data); }
private:
void scramble(double* data) {
int j = 1;
for (int i = 1; i < 2 * N; i += 2) {
if (j > i) {
std::swap(data[j - 1], data[i - 1]);
std::swap(data[j], data[i]);
}
int m = N;
while (m >= 2 && j > m) {
j -= m;
m >>= 1;
}
j += m;
}
}
DanielsonLanczos<N, -1> recursion;
DanielsonLanczos<N, 1> irecursion;
};
AxisymmetricGrid r; AxisymmetricGrid rho0; CenteredGrid fourPikEdeltaJ; RadialGrid ur; TangentialGrid ut; CenteredGrid deltapOverRho0; CenteredGrid deltaT; RadialGrid convr; TangentialGrid convt; RadialGrid buoyancyr; RadialGrid corior; TangentialGrid coriot; CenteredGrid deltaTconv; CenteredGrid deltaTsource;We now want to implement a function to compute the convection and force terms $(Conv.)_r$, $(Conv.)_{\theta}$, $(Buo.)_r$, $(Corio.)_r$, $(Corio.)_{\theta}$ from the current velocity and temperature. Note that these terms involve averaging and partial differentiation operators that may require the velocity and the temperature outside the grid itself. Thus, before computing these terms, we make sure these outside values are correctly set, by calling the setBoundaryConditions method on each grid (note that we set $u_r=0$ outside the domain, unlike what is explained in the Natural convection and wind speed section; in fact these values are not used in our implementation because the radial grid assignment operator does not evaluate values on the boundary).
After boundary conditions are set, computing the convection and force terms is straightforward, thanks to our grid expression template classes. For instance, the first term of $(Conv.)_r$, equal to $\frac{1}{r}\partial_r \overline{r\rho_0 u_r}^r \overline{u_r}^r$, is simply implemented as partialR(avgR(r * rho0 * ur) * avgR(ur)) / r:
void computeForces() {
ur.setBoundaryConditions(0.0);
ut.setBoundaryConditions(-1.0);
deltaT.setBoundaryConditions(1.0);
convr = (partialR(avgR(r * rho0 * ur) * avgR(ur)) +
partialT(avgR(rho0 * ut) * avgT(ur)) -
avgR(avgT(rho0 * ut) * avgT(ut))) / r;
convt = (partialR(avgT(r * rho0 * ur) * avgR(ut)) +
partialT(avgT(rho0 * ut) * avgT(ut)) +
avgT(avgT(rho0 * ut) * avgR(ur))) / r;
buoyancyr = rho0 * r * (omega * omega / T0) * avgR(deltaT);
corior = (-2.0 * omega) * avgR(rho0 * avgT(ut));
coriot = (2.0 * omega) * avgT(rho0 * avgR(r * ur)) / r;
deltaTconv = (partialR(r * rho0 * ur * cv * avgR(deltaT)) +
partialT(rho0 * ut * cv * avgT(deltaT))) / r;
deltaTsource =
fourPikEdeltaJ - (16.0 * kappaE * sigma * T0 * T0 * T0) * rho0 * deltaT;
}
TridiagonalSolver tridiagonalSolvers[Nt];
void initPressureSolver() {
for (int k = 0; k < Nt; ++k) {
double cosWk = cos((2.0 * M_PI * k) / Nt);
for (int i = 0; i < Nr; ++i) {
double a = r(i + 1) * rho0(i + 1) / (DELTAr * DELTAr);
double c = r(i + 2) * rho0(i + 2) / (DELTAr * DELTAr);
double d = rho0(i + 1 + half) / (r(i + 1 + half) * DELTAt * DELTAt);
double b = -(a + c + 2.0 * d) + 2.0 * d * cosWk;
if (i == 0) {
tridiagonalSolvers[k].init(i, 0.0, a + b, c);
} else if (i == Nr - 1) {
tridiagonalSolvers[k].init(i, a, b + c, 0.0);
} else {
tridiagonalSolvers[k].init(i, a, b, c);
}
}
}
}
The pressure itself is computed in the computePressure function below. This
function first computes the centered grid $g$ which is the right hand side of
the pressure's Poisson equation. Since this involves partial differentiation
operators that may need values outside the grid, we first make sure these values
are defined by using the setBoundaryConditions method. Note that we don't add
the extra term $\frac{r_s}{\Delta r}\left((Conv.)_r+(Buo.)_r+
(Corio.)_r\right)_{N_r+1,j+\frac{3}{2}}$ coming from the pressure boundary
condition because, in our implementation, it is equal to 0 by construction (see
the RadialGrid assignment operator).
We then compute the Fourier transform of each column of $g$ in f. Once this is done we solve the tri-diagonal linear systems, one per row in f. Finally, we perform an inverse Fourier transform for each column of f, yielding the pressure $\delta p$ (divided by $\rho_0$). At the end we use the setBoundaryConditions method on the pressure field so that we can later safely compute its gradient. Since the extra term mentionned above is 0, these boundary conditions are simple Neumann boundary conditions.
void computePressure(double dt) {
static FFT<Nt> fftSolver;
static double f[Nr][2 * Nt];
static double x[2 * Nr];
static CenteredGrid g;
convr.setBoundaryConditions(0.0);
convt.setBoundaryConditions(1.0);
buoyancyr.setBoundaryConditions(0.0);
corior.setBoundaryConditions(0.0);
coriot.setBoundaryConditions(1.0);
g = partialR(r * (rho0 * ur * (1.0 / dt) - convr - corior - buoyancyr)) +
partialT(rho0 * ut * (1.0 / dt) - convt - coriot);
for (int i = 0; i < Nr; ++i) {
for (int j = 0; j < Nt; ++j) {
f[i][2 * j] = g(i + 1 + half, j + 1 + half);
f[i][2 * j + 1] = 0.0;
}
fftSolver.fft(f[i]);
}
for (int j = 0; j < Nt; ++j) {
for (int i = 0; i < Nr; ++i) {
x[2 * i] = f[i][2 * j];
x[2 * i + 1] = f[i][2 * j + 1];
}
tridiagonalSolvers[j].solve(x);
for (int i = 0; i < Nr; ++i) {
f[i][2 * j] = x[2 * i];
f[i][2 * j + 1] = x[2 * i + 1];
}
}
for (int i = 0; i < Nr; ++i) {
fftSolver.ifft(f[i]);
for (int j = 0; j < Nt; ++j) {
deltapOverRho0(i + 1 + half, j + 1 + half) = f[i][2 * j] / Nt;
}
}
deltapOverRho0.setBoundaryConditions(1.0);
}
k = f(y0) y1 = y0 + B[0] * k for i = 1; i < 4 k = f(y0 + A[i] * k) y1 += B[i] * kwhere A={0,1/2,1/2,1} and B={1/6,1/3,1/3,1/6}. In our case $y$ is the fluid state given by $u_r$, $u_{\theta}$ and $\delta T$ and $f$ is the time derivative of these fields. The Runge-Kutta method for our case is therefore:
void timeStep(double dt) {
static RadialGrid ur0, dUr, ur1;
static TangentialGrid ut0, dUt, ut1;
static CenteredGrid deltaT0, dDeltaT, deltaT1;
static const double A[] = { 0.0, 1.0 / 2.0, 1.0 / 2.0, 1.0 };
static const double B[] = { 1.0 / 6.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 6.0 };
ur0 = ur;
ut0 = ut;
deltaT0 = deltaT;
computeForces();
computePressure(dt);
dUr = (convr + corior + buoyancyr + rho0 * partialR(deltapOverRho0)) / rho0;
dUt = (convt + coriot + rho0 * partialT(deltapOverRho0) / r) / rho0;
dDeltaT = (deltaTconv - deltaTsource) / (rho0 * cv);
ur1 = ur0 - dt * B[0] * dUr;
ut1 = ut0 - dt * B[0] * dUt;
deltaT1 = deltaT0 - dt * B[0] * dDeltaT;
for (int i = 1; i < 4; ++i) {
ur = ur0 - dt * A[i] * dUr;
ut = ut0 - dt * A[i] * dUt;
deltaT = deltaT0 - dt * A[i] * dDeltaT;
computeForces();
computePressure(dt);
dUr = (convr + corior + buoyancyr + rho0 * partialR(deltapOverRho0)) / rho0;
dUt = (convt + coriot + rho0 * partialT(deltapOverRho0) / r) / rho0;
dDeltaT = (deltaTconv - deltaTsource) / (rho0 * cv);
ur1 = ur1 - dt * B[i] * dUr;
ut1 = ut1 - dt * B[i] * dUt;
deltaT1 = deltaT1 - dt * B[i] * dDeltaT;
}
ur = ur1;
ut = ut1;
deltaT = deltaT1;
}
template<class E> double integral(const GridExpr<E>& g) {
double sum = 0.0;
for (int i = 1; i <= Nr; ++i) {
for (int j = 1; j <= Nt; ++j) {
sum += r(i + half) * g(i + half, j + half);
}
}
return sum * DELTAr * DELTAt;
}
void computeEnergy(double* kinetic, double* internal, double *potential) {
ur.setBoundaryConditions(0.0);
ut.setBoundaryConditions(-1.0);
*kinetic =
integral(0.5 * (avgR(r * rho0 * ur * ur) + avgT(r * rho0 * ut * ut)) / r);
*internal = integral(rho0 * cv * deltaT);
*potential = integral(rho0 * r * r * (omega * omega / (2.0 * T0)) * deltaT);
}
void init() {
for (int i = 0; i <= Nr + 2; ++i) {
const double C = omega * omega / (2.0 * Rs * T0);
double ri = (i - 1.0) * DELTAr;
r(i) = ri;
rho0(i) = ps / (Rs * T0) * exp(-C * (rs * rs - ri * ri));
if (i < Nr + 2) {
double rh = (i - 0.5) * DELTAr;
r(i + half) = rh;
rho0(i + half) = ps / (Rs * T0) * exp(-C * (rs * rs - rh * rh));
}
}
for (int i = 1; i <= Nr; ++i) {
double ri = r(i + half);
for (int j = 1; j <= Nt; ++j) {
double thetaj = (j - 0.5) * DELTAt;
// Approximate the temperature variations found using radiative transfer
// for a static atmosphere (without the spikes near the light sources).
double deltaT = -0.39 * pow(ri / rs, 3.2) * cos(3.0 * thetaj);
fourPikEdeltaJ(i + half, j + half) =
16.0 * (kappaE * sigma * T0 * T0 * T0) * rho0(i + half) * deltaT;
}
}
}
The main simulation loop is then straightforward to implement. The only subtelty
is that the time step dt should be adapted to the current maximum velocity, for
the stability of the method. We set it to 10% of the time needed to flow
through the largest grid cell at the maximum atmosphere velocity.
int main(int argc, char* argv[]) {
init();
initPressureSolver();
const double maxCellSize = std::min(rs / Nr, 2.0 * M_PI * rs / (3.0 * Nt));
double t = 0;
while (true) {
double kinetic;
double internal;
double potential;
computeEnergy(&kinetic, &internal, &potential);
printf("Energy at t=%g s: %g %g %g / Total: %g (J.m^-1)\n",
t, kinetic, internal, potential, kinetic + internal + potential);
double vMax = 0.0;
for (int i = 1; i <= Nr; ++i) {
for (int j = 1; j <= Nt; ++j) {
double vr = (ur(i, j + half) + ur(i + 1, j + half)) * 0.5;
double vt = (ut(i + half, j) + ut(i + half, j + 1)) * 0.5;
vMax = fmaxf(vMax, sqrt(vr * vr + vt * vt));
}
}
double dt = 0.1 * maxCellSize / std::max(1.0, vMax);
timeStep(dt);
t += dt;
}
return 0;
}