 |
Bimetric 3+1 toolkit for spherical symmetry
|
|
Bimetric 3+1 toolkit for spherical symmetry
|
Data Structures | |
| class | MoL |
| Time evolution using the Method of Lines (MoL). More... | |
| class | BandLUDecomposition |
BandLUDecomposition implements a method for solving linear equations A * x = b for a band-diagonal matrix A using LU decomposition. More... | |
| class | CubicSpline |
| CubicSpline encapsulates the natural cubic spline of degree 3 with continuity C2. More... | |
| class | Vector< T > |
| A vector (or an array) of elements. More... | |
| class | Matrix< T > |
| A matrix (a two-dimensional array) of elements. More... | |
| class | IntegFace |
| Interface used for MoL integration that prepares and finalizes integration steps. More... | |
| struct | MoLDescriptor |
| A method of lines (MoL) descriptor. More... | |
Defines | |
| #define | GF_r(f) ( 0.5 * ( f(m,n+1) - f(m,n-1) ) * inv_delta_r ) |
| #define | GF_rr(f) ( ( f(m,n+1) - 2 * f(m,n) + f(m,n-1) ) * inv_delta_rr ) |
| #define | extrapolate_R(f, m, n) |
| extrapolate_R is an optimized version of 4th order in accuracy extrapolation using the 4th order Taylor expansion. | |
| #define | extrapolate_lin(f, m, n) |
| extrapolate_lin is a linear extrapolation of the 4th order in accuracy (used for derivatives). | |
| #define | KreissOligerTerm(f, dt) ( - ( GF(f,m,n-1) - 2* GF(f,m,n) + GF(f,m,n+1) ) * dissip_delta_r * (dt) ) |
KreissOligerTerm is a macro that gives a Kreiss-Oliger dissipation term. | |
| #define | assertBounds(v, e) |
Typedefs | |
| typedef const Vector< Real > | VecReal_I |
| typedef Vector< Real > | VecReal |
| typedef Vector< Real > | VecReal_O |
| typedef Vector< Real > | VecReal_IO |
| typedef Vector< Int > | VecInt |
| typedef const Matrix< Real > | MatReal_I |
| typedef Matrix< Real > | MatReal |
| typedef Matrix< Real > | MatReal_O |
| typedef Matrix< Real > | MatReal_IO |
Functions | |
| void | CubicSpline_sanityCheck () |
| void | MoL::integStep_Euler (Int mNext, Int m) |
| Evolves the state variables using the FT scheme. | |
| void | MoL::integStep_AvgICN (Int mNext, Int mMid, Int m) |
| Evolves the sate variables using the ICN averaging. | |
| void | MoL::integrate_AvgICN (int ICN_steps) |
| Time evolution using the iterated Crank-Nicolson (ICN). | |
| void | MoL::integStep_MoL (Int m1, Int m, Real alpha) |
| The intermediate method of lines (MoL) steps. | |
| void | MoL::integStep_MoLInit (Int m1, Int m, Real beta_delta_t) |
| Prepares the state variables for the intermediate MoL steps. | |
| void | MoL::integrate_MoL (const MoLDescriptor &MoL) |
| Time evolution using the method of lines (MoL). | |
| BandLUDecomposition::BandLUDecomposition (MatReal_I &A, const Int mm1, const Int mm2) | |
| The constructor which does the LU decomposition. | |
| void | BandLUDecomposition::solve (VecReal_I &b, VecReal_O &x) |
Given a right-hand side vector b[0..n-1], solves the band-diagonal linear equation A * x = b. | |
| static void | BandLUDecomposition::sanityCheck () |
| Performs a sanity check test with sample data. | |
Variables | |
| MoLDescriptor | ICN2 |
| Iterative Crank-Nicolson, 2-iterations. | |
| MoLDescriptor | ICN3 |
| Iterative Crank-Nicolson, 3-iterations. | |
| MoLDescriptor | RK1 |
| Runge-Kutta, 1st order (the Euler method) | |
| MoLDescriptor | RK2 |
| Runge-Kutta, 2nd order. | |
| MoLDescriptor | RK3 |
| Runge-Kutta, 3rd order. | |
| MoLDescriptor | RK4 |
| Runge-Kutta, 4th order. | |
| static std::vector< MoL * > | MoL::knownIntegrators |
| Keep the list of all created integrator objects (used by signalHandler). | |
| static std::map< std::string, int > | MoL::knownMethods |
| A list of the implemented integration methods. | |
| #define assertBounds | ( | v, | |
| e | |||
| ) |
| #define extrapolate_lin | ( | f, | |
| m, | |||
| n | |||
| ) |
GF(f,m,n) = \
1./4. * GF(f,m,n-5) - 4./3. * GF(f,m,n-4) + 3 * GF(f,m,n-3) \
- 4 * GF(f,m,n-2) + 37./12. * GF(f,m,n-1)
extrapolate_lin is a linear extrapolation of the 4th order in accuracy (used for derivatives).
Definition at line 72 of file finiteDifferences.h.
| #define extrapolate_R | ( | f, | |
| m, | |||
| n | |||
| ) |
GF(f,m,n) = \
- 7./48. * GF(f,m,n-8) + 209./144. * GF(f,m,n-7) - 103./16. * GF(f,m,n-6) \
+ 793./48. * GF(f,m,n-5) - 3835./144. * GF(f,m,n-4) + 439./16. * GF(f,m,n-3) \
- 281./16. * GF(f,m,n-2) + 917./144. * GF(f,m,n-1)
extrapolate_R is an optimized version of 4th order in accuracy extrapolation using the 4th order Taylor expansion.
Definition at line 64 of file finiteDifferences.h.
| #define GF_r | ( | f | ) | ( 0.5 * ( f(m,n+1) - f(m,n-1) ) * inv_delta_r ) |
Definition at line 33 of file finiteDifferences.h.
| #define GF_rr | ( | f | ) | ( ( f(m,n+1) - 2 * f(m,n) + f(m,n-1) ) * inv_delta_rr ) |
Definition at line 35 of file finiteDifferences.h.
| #define KreissOligerTerm | ( | f, | |
| dt | |||
| ) | ( - ( GF(f,m,n-1) - 2* GF(f,m,n) + GF(f,m,n+1) ) * dissip_delta_r * (dt) ) |
KreissOligerTerm is a macro that gives a Kreiss-Oliger dissipation term.
Coefficients d_k at gridpoints for the centred finite-difference discretisation of the derivative partial_x^{2N} (of the dissipative Kreiss-Oliger operator).
---------------------------------------------
N -4 -3 -2 -1 0 +1 +2 +3 +4
---------------------------------------------
1 +1 -2 +1
2 -1 +4 -6 +4 -1
3 +1 -6 +15 -20 +15 -6 +1
4 -1 +8 -28 +56 -70 +56 -28 +8 -1
--------------------------------------------- Here, N >= 1 is an integer and D_h^{2N} is a centered difference operator approximating a partial spatial derivative of order 2N, e.g., a second- (N = 1) or fourth-order (N = 2) spatial derivative.
To be subtracted from the evolution equation:
`epsilon (-1)^N delta_t {delta_x}^{2N-1} D_h^{2N}`where:
`D_h^{2N} = (2 delta_x)^{-2N} * ( sum_k d_k f[n+k] )` For N = 2, subtract from the evolution equation u_t:
`dissip(f,dt) = epsilon 2^{-2N} * dt/delta_x * ( sum_k d_k f[n+k] )`Definition at line 108 of file finiteDifferences.h.
| typedef Matrix<Real> MatReal_IO |
| typedef Vector<Real> VecReal_IO |
| BandLUDecomposition::BandLUDecomposition | ( | MatReal_I & | A, |
| const Int | mm1, | ||
| const Int | mm2 | ||
| ) |
The constructor which does the LU decomposition.
Given an n*n band-diagonal matrix A with m1 subdiagonal rows and m2 superdiagonal rows, compactly stored in the array A[0..n-1][0..m1+m2], an LU decomposition of a rowwise permutation of A is constructed.
It takes as arguments the compactly stored matrix A, and the integers m1 >= 0 (subdiagonal size) and m2 >= 0 (superdiagonal size). The upper and lower triangular matrices are stored in U and L, respectively.
The upper and lower triangular matrices are stored in U and L, respectively. The stored vector indx[0..n-1] records the row permutation effected by the partial pivoting; d is +/-1 depending on whether the number of row interchanges was even or odd, respectively.
Some pivoting is possible within the storage limitations, and the routine does take advantage of the opportunity. Also, when TINY is returned as a diagonal element of U, then the original matrix (perhaps as modified by roundoff error) is in fact singular.
| A | Compactly stored diagonal-band matrix |
| mm1 | Size of subdiagonal |
| mm2 | Size of superdiagonal |
Definition at line 102 of file bandSol.h.
: nn(A.nrows()), m1(mm1), m2(mm2), L(nn,m1), U(A), indx(nn) { const Real TINY = 1e-40; Int mm = m1 + m2 + 1; Int l = m1; for( Int i = 0; i < m1; ++i ) { for( Int j = m1 - i; j < mm; ++j ) { U[i][j-l] = U[i][j]; } l--; for( Int j = mm - l - 1; j < mm; ++j ) { U[i][j] = 0; } } Real d = 1; l = m1; for( Int k = 0; k < nn; ++k ) { Real dum = U[k][0]; Int i = k; if( l < nn ) { l++; } for( Int j = k + 1; j < l; ++j ) { if( abs( U[j][0] ) > abs( dum ) ) { dum = U[j][0]; i = j; } } indx[k] = i + 1; if( dum == 0 ) { U[k][0]=TINY; } if (i != k) { d = -d; for( Int j = 0; j < mm; ++j ) { SWAP( U[k][j], U[i][j] ); } } for( i = k + 1; i < l; ++i ) { dum = U[i][0] / U[k][0]; L[k][i-k-1] = dum; for( Int j = 1; j < mm; ++j ) { U[i][j-1] = U[i][j] - dum * U[k][j]; } U[i][mm-1] = 0; } } }
| void CubicSpline_sanityCheck | ( | ) |
Definition at line 117 of file cubicSpline.h.
{
Real pts[5][2] =
{
{0.005, 4.892186615295462}, {0.045, 44.02967953765916},
{0.105, 102.7359189212047 }, {0.145, 40.74442799816186},
{0.195, 10.719645882407347}
};
CubicSpline ncs( 5 );
ncs.initialize( pts );
ncs.dump ();
for( Real x = 0; x < 0.2001; x += 0.01 ) {
std::cout << x << "\t" << ncs( x ) << std::endl;
}
}
| void MoL::integrate_AvgICN | ( | int | ICN_steps | ) | [private] |
Time evolution using the iterated Crank-Nicolson (ICN).
Integrates the evolution equations using the iterated Crank-Nicolson (ICN) method.
If ICN_steps == 0, then the method effectively becomes FTCS. See [9], [3]
Definition at line 418 of file integrator.h.
{
slog << "Employing Method of Lines, Iterative Crank-Nicolson (avg)"
<< ", N = " << ICN_steps << std::endl << std::endl;
cur_t = t_0; // The initial `t`
running = true; // Enable integration
for( Int n = 0; n < nTotal; ++n ) {
GF( fld::t, 0, n ) = t_0;
}
Int mStep = 0; // The step counter (used to filter checkpoints)
while( running && abs(cur_t) < abs(t_1) )
{
for( Int m = 0; running && m < mLen && abs(cur_t) < abs(t_1); ++m )
{
Int mNext = m + 1 >= mLen ? 0 : m + 1;
/////////////////////////////////////////////////////////////////////////////
// Predictor. Use FTCS for the first step.
// Note that this will be the only step, if ICN_steps == 0.
//
Int m_1 = ICN_steps >= 1 ? mLen + 1 : mNext;
integStep_Begin ( m ); if( ! running ) break;
integStep_Euler ( m_1, m );
integStep_End ( m_1, m );
/////////////////////////////////////////////////////////////////////////////
// Corrector iterations. There will be no iterations, if ICN_steps == 0.
//
for( int i = 0; running && i < ICN_steps; ++i )
{
Int m_2 = ( i + 1 >= ICN_steps ? mNext : m_1 + 1 );
integStep_Begin ( m_1 ); if( ! running ) break;
integStep_AvgICN ( m_2, m_1, m );
integStep_End ( m_2, m_1 );
m_1 = m_2;
}
if( ( mStep++ % output->get_mSkip () ) == 0 ) {
integStep_Checkpoint( m );
}
cur_t = GF( fld::t, m, nGhost ); // The current time step
}
}
}
| void MoL::integrate_MoL | ( | const MoLDescriptor & | MoL | ) | [private] |
Time evolution using the method of lines (MoL).
The method of lines (MoL) separates the time integration from the rest of an evolution scheme.
A N-step MoL method evolves the equation:
partial_t Y = L[Y], using the following algorithm:
Y^{(0)} = Y^{m},
Y^{(i+1)} = Sum_{j=0}^{i} ( alpha_{ij} Y^{(j)} ) + <- integStep_MoL
Delta t beta_i L[Y^{(i)}], for i = 0, ..., N-1, <- integStep_MoLInit
Y^{m+1} = Y^{(N)}. The variables Y^{(i)}, i = 0, ..., N, are intermediate.
The method is completely specified by N, and arrays alpha and beta in the structure MoLDescriptor.
Definition at line 522 of file integrator.h.
{
slog << "Employing Method of Lines, " << MoL.name
<< ", N = " << MoL.N << std::endl << std::endl;
cur_t = t_0; // The initial `t`
running = true; // Enable integration
for( Int n = 0; n < nTotal; ++n ) {
GF( fld::t, 0, n ) = t_0;
}
// The array mk[] contains redirections to grid rows (m's):
// mk[0] points to Y^{m} == Y^{(0)},
// mk[1], ..., mk[Mol.N-1] point to the intermediate Y^{(1)}, ..., Y^{(N-1)},
// mk[MoL.N] points to Y^{m+1} == Y^{(N)}.
//
int mk[ 10 ]; // int mk[ MoL.N + 1 ]; // Note that N is at most 8
for( int i = 0; i <= MoL.N; ++i ) {
mk[i] = mLen + i;
}
Int mStep = 0; // The step counter
while( running && abs(cur_t) < abs(t_1) )
{
for( Int m = 0; running && m < mLen && abs(cur_t) < abs(t_1); ++m )
{
mk[0] = m; // Y^{(m)}
mk[MoL.N] = m + 1 >= mLen ? 0 : m + 1; // Y^{(m+1)}
for( int i = 0; running && i < MoL.N; ++i )
{
integStep_Begin( mk[i] ); if( ! running ) break;
integStep_MoLInit( mk[i+1], mk[i], MoL.beta[i] * delta_t );
for( int j = 0; j <= i; j++ ) {
integStep_MoL( mk[i+1], mk[j], MoL.alpha[i][j] );
}
integStep_End( mk[i+1], mk[i] );
}
if( ( mStep++ % output->get_mSkip () ) == 0 ) {
integStep_Checkpoint( m );
}
cur_t = GF( fld::t, m, nGhost ); // The current time step
}
}
}
| void MoL::integStep_AvgICN | ( | Int | mNext, |
| Int | mMid, | ||
| Int | m | ||
| ) | [private] |
Evolves the sate variables using the ICN averaging.
Definition at line 388 of file integrator.h.
{
const Real half_dt = delta_t / 2;
OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n )
{
GF( fld::t, m2, n ) = GF( fld::t, m, n ) + delta_t; // Evolve time
for( auto f: constantGF )
{
GF( f, m2, n ) = GF( f, m, n );
}
for( auto e: evolvedGF )
{
GF( e.f, m2, n ) = GF( e.f, m, n )
+ half_dt * ( GF( e.f_t, m1, n) + GF( e.f_t, m, n ) )
- KreissOligerTerm( e.f, delta_t );
}
}
}
| void MoL::integStep_Euler | ( | Int | mNext, |
| Int | m | ||
| ) | [private] |
Evolves the state variables using the FT scheme.
Definition at line 368 of file integrator.h.
{
OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n )
{
GF( fld::t, m1, n ) = GF( fld::t, m, n ) + delta_t; // Evolve time
for( auto f: constantGF )
{
GF( f, m1, n ) = GF( f, m, n );
}
for( auto e: evolvedGF )
{
GF( e.f, m1, n ) = GF( e.f, m, n ) + delta_t * GF( e.f_t, m, n )
- KreissOligerTerm( e.f, delta_t );
}
}
}
| void MoL::integStep_MoL | ( | Int | m1, |
| Int | m, | ||
| Real | alpha | ||
| ) | [private] |
| void MoL::integStep_MoLInit | ( | Int | m1, |
| Int | m, | ||
| Real | beta_delta_t | ||
| ) | [private] |
Prepares the state variables for the intermediate MoL steps.
Definition at line 483 of file integrator.h.
{
OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n )
{
// There is no explicit dependence on t, and t_t(m,n) = 1.
// Also, r is constant along the `line`, i.e., r_t(m,n) = 0.
//
GF( fld::t, m1, n ) = bc_delta_t; // Evolve time
for( auto f: constantGF )
{
GF( f, m1, n ) = GF( f, m, n );
}
for( auto e: evolvedGF )
{
GF( e.f, m1, n ) = bc_delta_t * GF( e.f_t, m, n )
- KreissOligerTerm( e.f, bc_delta_t );
}
}
}
| void BandLUDecomposition::sanityCheck | ( | ) | [static] |
Performs a sanity check test with sample data.
Sanity check for BandLUDecomposition.
Definition at line 206 of file bandSol.h.
{
const Real A_values[] = {
// -2 -1 0 +1 <--- the matrix diagonal lies at [0]
0, 0, 3, 1,
0, 4, 1, 5,
9, 2, 6, 5,
3, 5, 8, 9,
7, 9, 3, 2,
3, 8, 4, 6,
2, 4, 4, 0
};
const Real b_values[] = {
7, 6, 5, 4, 3, 2, 1
};
MatReal_I A( 7, 4, A_values );
VecReal_I b( 7, b_values );
VecReal_O x( 7 );
BandLUDecomposition( A, 2, 1 ).solve( b, x );
for( Int i = 0; i < 7; ++i ) {
std::cout << x[i] << ", ";
}
std::cout << std::endl << std::endl << "Compare to:" << std::endl << std::endl
<< "4.66912, -7.00737, -1.13382, -3.24088, 6.29092, 10.616, -13.5114"
<< std::endl << std::endl;
}
| void BandLUDecomposition::solve | ( | VecReal_I & | b, |
| VecReal_O & | x | ||
| ) |
Given a right-hand side vector b[0..n-1], solves the band-diagonal linear equation A * x = b.
Definition at line 165 of file bandSol.h.
{
Int mm = m1 + m2 + 1;
Int l = m1;
for( Int k = 0; k < nn; ++k ) {
x[k] = b[k];
}
for( Int k = 0; k < nn; ++k )
{
Int j = indx[k] - 1;
if( j != k ) {
SWAP( x[k], x[j] );
}
if( l < nn ) {
l++;
}
for( j = k + 1; j < l; ++j ) {
x[j] -= L[k][j-k-1] * x[k];
}
}
l = 1;
for( Int i = nn - 1; i >= 0; --i )
{
Real dum = x[i];
for( Int k = 1; k < l; ++k ) {
dum -= U[i][k]*x[k+i];
}
x[i] = dum / U[i][0];
if( l < mm ) {
l++;
}
}
}
{
"Iterative Crank-Nicolson", 2,
{ { 1., 0. },
{ 0., 1. } },
{ 1./2., 1. }
}
Iterative Crank-Nicolson, 2-iterations.
Definition at line 76 of file methodOfLines.h.
{
"Iterative Crank-Nicolson", 3,
{ { 1., 0., 0. },
{ 0., 1., 0. },
{ 0., 0., 1. } },
{ 1./2., 1./2., 1. }
}
Iterative Crank-Nicolson, 3-iterations.
Definition at line 84 of file methodOfLines.h.
std::vector< MoL * > MoL::knownIntegrators [static, private] |
Keep the list of all created integrator objects (used by signalHandler).
Definition at line 55 of file integrator.h.
std::map< std::string, int > MoL::knownMethods [static, private] |
{
{ "Euler", 0 }, { "avgICN2", 1 },
{ "avgICN3", 2 }, { "avgICN4", 3 },
{ "ICN2", 4 }, { "ICN3", 5 },
{ "RK1", 6 }, { "RK2", 7 },
{ "RK3", 8 }, { "RK4", 9 }
}
A list of the implemented integration methods.
Definition at line 51 of file integrator.h.
{
"Runge-Kutta", 1,
{ { 1. } },
{ 1. }
}
Runge-Kutta, 1st order (the Euler method)
Definition at line 93 of file methodOfLines.h.
{
"Runge-Kutta", 2,
{ { 1., 0. },
{ 1./2., 1./2. } },
{ 1., 1./2. }
}
Runge-Kutta, 2nd order.
Definition at line 100 of file methodOfLines.h.
{
"Runge-Kutta", 3,
{ { 1., 0., 0. },
{ 3./4., 1./4., 0. },
{ 1./3., 0., 2./3. } },
{ 1., 1./4., 2./3. }
}
Runge-Kutta, 3rd order.
Definition at line 108 of file methodOfLines.h.
{
"Runge-Kutta", 4,
{ { 1., 0., 0., 0. },
{ 1., 0., 0., 0. },
{ 1., 0., 0., 0. },
{ -1./3., 1./3., 2./3., 1./3. } },
{ 1./2., 1./2., 1., 1./6. }
}
Runge-Kutta, 4th order.
Definition at line 117 of file methodOfLines.h.
1.8.0