 |
Bimetric 3+1 toolkit for spherical symmetry
|
|
Bimetric 3+1 toolkit for spherical symmetry
|
00001 /** 00002 * @file bim-solver.cpp 00003 * @brief Standard 3+1 evolution for spherically symmetric bimetric spacetimes. 00004 * @author Mikica Kocic 00005 * @copyright GNU General Public License (GPLv3). 00006 */ 00007 00008 #include <iostream> 00009 #include <iomanip> 00010 #include <fstream> 00011 #include <string> 00012 #include <cmath> 00013 #include <cstdio> 00014 #include <chrono> 00015 00016 #include "numMethods.h" // The implemented numerical modules 00017 #include "sys/slog.h" // For writing both to cerr and cout simultaneously 00018 #include "sys/trackUsedTime.h" // Keep track of the elapsed time of the application 00019 #include "sys/paramsHolder.h" // Holds 'key=value' pairs got from the parameter file 00020 #include "sys/hpc.h" // High-Performance Computing (HPC) support 00021 00022 ///////////////////////////////////////////////////////////////////////////////////////// 00023 // Declare our grid functions (they must be in advance known to the grid-driver) 00024 00025 #include "grid/gridFunctions.h" 00026 00027 namespace fld 00028 { 00029 ///////////////////////////////////////////////////////////////////////////////////// 00030 /// The bimetric grid functions (indices on the grid). 00031 enum bimIndex { bimFirst = GFCNT - 1, 00032 ///////////////////////////////////////////////////////////////////////////////////// 00033 00034 gA, gB, gK, gKD, gDA, gDB, gSig, //!< Evolved variables in the `g`-sector 00035 fA, fB, fK, fKD, fDA, fDB, fSig, //!< Evolved variables in the `f`-sector 00036 pfD, pfS, pftau, //!< State variables for the PF 00037 00038 p, //!< Separation between two metrics (relative shift) 00039 q, //!< Overall (geometric mean) shift 00040 gAlp, gDAlp, fAlp, fDAlp, //!< Lapses and their log-derivatives 00041 00042 Lt, //!< Lorentz factor (derived from `p`) 00043 R, //!< Ratio `fB/gB` 00044 pfv, //!< The three-velocity of the PF 00045 00046 mpiBoundary = pfv, //!< The last grid function for which BC must be fixed. 00047 //!< Used in MPI when calling defineGhostChunk. 00048 00049 g_rho, g_JK1, g_JK2, //!< Cumulative sources in th `g`-sector 00050 f_rho, f_JK1, f_JK2, //!< Cumulative sources in the `f`-sector 00051 00052 gW, fW, cW, //!< Lapse ratio, `( gW * gAlp + fW * fAlp ) / cW = 0` 00053 00054 // Diagnostics 00055 // 00056 C_1, C_2, C_3, C_4, //!< Constraints (more precise: constraint violations) 00057 gHorz, fHorz, //!< Apparent horizon finders 00058 p_g, p_f, //!< Alt. expr. for the separation between two metrics 00059 00060 // The RHS of the evolution equations 00061 // 00062 gA_t, gB_t, gK_t, gKD_t, gDA_t, gDB_t, gSig_t, 00063 fA_t, fB_t, fK_t, fKD_t, fDA_t, fDB_t, fSig_t, 00064 pfD_t, pfS_t, pftau_t, p_t, 00065 gAlp_t, gDAlp_t, 00066 00067 // The cached values of the spatial derivatives 00068 // 00069 gA_r, gB_r, gK_r, gKD_r, gDA_r, gDB_r, 00070 fA_r, fB_r, fK_r, fKD_r, fDA_r, fDB_r, 00071 pfD_r, pfS_r, pftau_r, pfv_r, 00072 p_r, p_rr, q_r, q_rr, eq_pr_r, eq_qr_r, 00073 gAlp_r, gDAlp_r, fAlp_r, fAlp_rr, fDAlp_r, 00074 00075 // Grid functions uses for storing temporary data 00076 // 00077 tmp, gdbg, fdbg, 00078 00079 ///////////////////////////////////////////////////////////////////////////////////// 00080 bimLast }; 00081 #undef GFCNT // Clear the old number of grid functions then... 00082 #define GFCNT fld::bimLast // update the number of grid functions on a grid point. 00083 ///////////////////////////////////////////////////////////////////////////////////// 00084 00085 /** The grid functions that will be read from the initial data. 00086 */ 00087 static const std::vector<Int> bimInput = { 00088 gA, gB, gK, gKD, gDA, gDB, gSig, 00089 fA, fB, fK, fKD, fDA, fDB, fSig, 00090 q, gAlp, fAlp, gDAlp, 00091 p, 00092 pfD, pfS, pftau 00093 }; 00094 00095 /** The grid functions that are involved in time. 00096 */ 00097 static const std::vector<EvolvedBy> bimEvolvedGF = { 00098 { p, p_t }, 00099 { gA, gA_t }, { gB, gB_t }, { gK, gK_t }, { gKD, gKD_t }, 00100 { gDA, gDA_t }, { gDB, gDB_t }, { gSig, gSig_t }, 00101 { fA, fA_t }, { fB, fB_t }, { fK, fK_t }, { fKD, fKD_t }, 00102 { fDA, fDA_t }, { fDB, fDB_t }, { fSig, fSig_t }, 00103 { pfD, pfD_t }, { pfS, pfS_t }, { pftau, pftau_t } 00104 }; 00105 00106 /** The grid functions which will be written to the output. 00107 */ 00108 static const std::vector<GF_Descriptor> bimOutput = 00109 { 00110 { gA, "gA", "A" }, 00111 { gB, "gB", "B" }, 00112 { gK, "gK", "K" }, 00113 { gKD, "gKD", "K_\\Delta" }, 00114 { fA, "fA", "\\tilde A" }, 00115 { fB, "fB", "\\tilde B" }, 00116 { fK, "fK", "\\tilde K" }, 00117 { fKD, "fKD", "\\tilde K_\\Delta" }, 00118 // 00119 { gA_t, "gA_t", "\\partial_t A" }, 00120 { gB_t, "gB_t", "\\partial_t B" }, 00121 { gK_t, "gK_t", "\\partial_t K" }, 00122 { gKD_t, "gKD_t", "\\partial_t K_\\Delta" }, 00123 { fA_t, "fA_t", "\\partial_t \\tilde A" }, 00124 { fB_t, "fB_t", "\\partial_t \\tilde B" }, 00125 { fK_t, "fK_t", "\\partial_t \\tilde K" }, 00126 { fKD_t, "fKD_t", "\\partial_t \\tilde K_\\Delta" }, 00127 // 00128 { gDA, "gDA", "D_A" }, 00129 { gDB, "gDB", "D_B" }, 00130 { gSig, "gSig", "\\sigma_g" }, 00131 { gW, "gW", "W_g" }, 00132 { fDA, "fDA", "\\tilde D_A" }, 00133 { fDB, "fDB", "\\tilde D_B" }, 00134 { fSig, "fSig", "\\sigma_f" }, 00135 { fW, "fW", "W_f" }, 00136 // 00137 { gDA_t, "gDA_t", "\\partial_t D_A" }, 00138 { gDB_t, "gDB_t", "\\partial_t D_B" }, 00139 { gSig_t, "gSig_t", "\\partial_t \\sigma_g" }, 00140 { gAlp, "gAlp", "\\alpha" }, 00141 { fDA_t, "fDA_t", "\\partial_t \\tilde D_A" }, 00142 { fDB_t, "fDB_t", "\\partial_t \\tilde D_B" }, 00143 { fSig_t, "fSig_t", "\\partial_t \\sigma_f" }, 00144 { fAlp, "fAlp", "\\tilde\\alpha" }, 00145 // 00146 { p, "p", "p" }, 00147 { p_t, "p_t", "\\partial_t p" }, 00148 { q, "q", "q" }, 00149 { cW, "cW", "W_c" }, 00150 { C_1, "C_1", "C_1" }, 00151 { C_2, "C_2", "C_2" }, 00152 { C_3, "C_3", "C_3" }, 00153 { C_4, "C_4", "C_4" }, 00154 // 00155 { g_rho, "g_rho", "\\rho" }, 00156 { g_JK1, "g_JK1", "{JK}_1" }, 00157 { g_JK2, "g_JK2", "{JK}_2" }, 00158 { gHorz, "gHorz", "z_g" }, 00159 { f_rho, "f_rho", "\\tilde{\\rho}" }, 00160 { f_JK1, "f_JK1", "\\tilde{JK}_1" }, 00161 { f_JK2, "f_JK2", "\\tilde{JK}_2" }, 00162 { fHorz, "fHorz", "z_f" }, 00163 // 00164 { pfD, "pfD", "\\hat D" }, 00165 { pfS, "pfS", "\\hat S" }, 00166 { pftau, "pftau", "\\hat \\tau" }, 00167 { pfv, "pfv", "\\hat v" }, 00168 { pfD_t, "pfD_t", "\\partial_t \\hat D" }, 00169 { pfS_t, "pfS_t", "\\partial_t \\hat S" }, 00170 { pftau_t, "pftau_t", "\\partial_t \\hat\\tau" }, 00171 { gDAlp, "gDAlp", "D_\\alpha" } 00172 }; 00173 00174 /** Additional output diagnostics (e.g., with spatial derivatives). 00175 */ 00176 static const std::vector<GF_Descriptor> bimShowDiagnostics = 00177 { 00178 { p_g, "p_g", "p_g" }, 00179 { gA_r, "gA_r", "\\partial_r A" }, 00180 { gB_r, "gB_r", "\\partial_r B" }, 00181 { gK_r, "gK_r", "\\partial_r K" }, 00182 { gKD_r, "gKD_r", "\\partial_r K_\\Delta" }, 00183 { gDA_r, "gDA_r", "\\partial_r D_A" }, 00184 { gDB_r, "gDB_r", "\\partial_r D_B" }, 00185 { gdbg, "gdbg", "g\\text{dbg}" }, 00186 // 00187 { p_f, "p_f", "p_f" }, 00188 { fA_r, "fA_r", "\\partial_r \\tilde A" }, 00189 { fB_r, "fB_r", "\\partial_r \\tilde B" }, 00190 { fK_r, "fK_r", "\\partial_r \\tilde K" }, 00191 { fKD_r, "fKD_r", "\\partial_r \\tilde K_\\Delta" }, 00192 { fDA_r, "fDA_r", "\\partial_r \\tilde D_A" }, 00193 { fDB_r, "fDB_r", "\\partial_r \\tilde D_B" }, 00194 { fdbg, "fdbg", "f\\text{dbg}" }, 00195 // 00196 { pfD_r, "pfD_r", "\\partial_r \\hat D" }, 00197 { pfS_r, "pfS_r", "\\partial_r \\hat S" }, 00198 { pftau_r, "pftau_r", "\\partial_r \\hat \\tau" }, 00199 { pfv_r, "pfv_r", "\\partial_r \\hat v" }, 00200 // 00201 { p_r, "p_r", "\\partial_r p" }, 00202 { p_rr, "p_rr", "\\partial_{rr} p" }, 00203 { q_r, "q_r", "\\partial_r q" }, 00204 { q_rr, "q_rr", "\\partial_{rr} q" }, 00205 { eq_pr_r, "eq_pr_r", "\\partial_r\\text{\"(p/r)\"}" }, 00206 { eq_qr_r, "eq_qr_r", "\\partial_r\\text{\"(q/r)\"}" }, 00207 // 00208 { gAlp_r, "gAlp_r", "\\partial_r \\alpha" }, 00209 { gDAlp_r, "gDAlp_r", "\\partial_r D_\\alpha" }, 00210 { fAlp_r, "fAlp_r", "\\partial_r \\tilde\\alpha" }, 00211 { fAlp_rr, "fAlp_rr", "\\partial_{rr} \\tilde\\alpha" }, 00212 { fDAlp_r, "fDAlp_r", "\\partial_r \\tilde D_\\alpha" } 00213 }; 00214 } 00215 00216 ///////////////////////////////////////////////////////////////////////////////////////// 00217 // After defining the grid functions, we can use the grid-driver and the other modules 00218 00219 #include "grid/gridDriver.h" 00220 #include "grid/gridInitialData.h" 00221 #include "grid/gridOutput.h" 00222 #include "grid/integrator.h" 00223 00224 ///////////////////////////////////////////////////////////////////////////////////////// 00225 00226 #include "bimetricModel.h" 00227 00228 /** BimetricEvolve encapsulates a 3+1 evolution solver for bimetric spacetimes. 00229 * 00230 * @todo EvolvedBy f, f_t --> add dissipation here? 00231 * @todo Implement PIRK? 00232 */ 00233 class BimetricEvolve 00234 : BimetricModel, // Based on the bimetric model 00235 GridUser, // To access grid functions on the grid 00236 public IntegFace // Implement the integration interface 00237 { 00238 enum Slicing //!< Known slicing methods 00239 { 00240 SLICE_CONSTG = 0, // Constant slice in g (f calculated) 00241 SLICE_CONSTGF = 1, // Constant slice in both g and f 00242 SLICE_MS2 = 2, // Maximal slicing, 2nd order FD 00243 SLICE_MS4 = 3, // Maximal slicing, 4th order FD 00244 SLICE_MS2OPT = 4, // Maximal slicing, 2nd order FD, optimized algorithm 00245 SLICE_BM = 5 // Bona-Masso slicing 00246 }; 00247 00248 Int slicing; //!< Select the slicing: maximal, Bona-Masso, ... 00249 Int lin2n; //!< Left grid-zone linear smoothing (default: `nGhost`) 00250 Int cub2n; //!< Left grid-zone cubic spline smoothing (default: `5*nGhost/2+6`) 00251 Real delta_t; //!< The integration step (obtained from the integrator) 00252 Int smooth; //!< Smooth the fields (level of smoothness) 00253 00254 ///////////////////////////////////////////////////////////////////////////////////// 00255 /** @defgroup g5 Macros to access data in a grid */ 00256 ///////////////////////////////////////////////////////////////////////////////////// 00257 /*@{*/ 00258 ///////////////////////////////////////////////////////////////////////////////////// 00259 00260 emitField(gA) emitField(gA_t) emitField(fA) emitField(fA_t) 00261 emitField(gB) emitField(gB_t) emitField(fB) emitField(fB_t) 00262 emitField(gK) emitField(gK_t) emitField(fK) emitField(fK_t) 00263 emitField(gKD) emitField(gKD_t) emitField(fKD) emitField(fKD_t) 00264 emitField(gDA) emitField(gDA_t) emitField(fDA) emitField(fDA_t) 00265 emitField(gDB) emitField(gDB_t) emitField(fDB) emitField(fDB_t) 00266 emitField(gSig) emitField(gSig_t) emitField(fSig) emitField(fSig_t) 00267 00268 emitField(p) emitField(p_g) emitField(p_f) emitField(p_t) 00269 00270 emitField(q) 00271 emitField(gAlp) emitField(gAlp_t) 00272 emitField(gDAlp) emitField(gDAlp_t) 00273 emitField(fAlp) emitField(fDAlp) 00274 emitField(gW) emitField(fW) emitField(cW) 00275 00276 emitField(g_rho) emitField(f_rho) 00277 emitField(g_JK1) emitField(f_JK1) 00278 emitField(g_JK2) emitField(f_JK2) 00279 00280 emitField(pfD) emitField(pfD_t) 00281 emitField(pfS) emitField(pfS_t) 00282 emitField(pftau) emitField(pftau_t) 00283 emitField(pfv) emitField(pfv_r) 00284 00285 emitField(C_1) emitField(C_2) emitField(C_3) emitField(C_4) 00286 emitField(gHorz) emitField(fHorz) emitField(Lt) emitField(R) 00287 /*@}*/ 00288 ///////////////////////////////////////////////////////////////////////////////////// 00289 /** @defgroup g6 Functions of the prime state variables */ 00290 ///////////////////////////////////////////////////////////////////////////////////// 00291 /*@{*/ 00292 emitDerivative_r( gA ) emitDerivative_r( fA ) 00293 emitDerivative_r( gB ) emitDerivative_r( fB ) 00294 emitDerivative_r( gK ) emitDerivative_r( fK ) 00295 emitDerivative_r( gKD ) emitDerivative_r( fKD ) 00296 emitDerivative_r( gDA ) emitDerivative_r( fDA ) 00297 emitDerivative_r( gDB ) emitDerivative_r( fDB ) 00298 // emitDerivative_r( gSig ) emitDerivative_r( fSig ) 00299 00300 emitDerivative_r( p ) 00301 emitDerivative_r( q ) 00302 emitDerivative_r( gAlp ) emitDerivative_r( fAlp ) 00303 emitDerivative_r( gDAlp ) emitDerivative_r( fDAlp ) 00304 00305 emitDerivative_r( pfD ) 00306 emitDerivative_r( pfS ) 00307 emitDerivative_r( pftau ) 00308 00309 // These are the only 2nd spatial derivatives which are needed: 00310 00311 emitDerivative_rr( p ) // used in eq_gDA_t and eq_fDA_t 00312 emitDerivative_rr( q ) // used in eq_gDA_t and eq_fDA_t 00313 emitDerivative_rr( fAlp ) // used in fDAlp_r 00314 00315 // The extrinsic curvature relations 00316 00317 #if 1 00318 Real gK1 ( Int m, Int n ) { return ( gK(m,n) + 2 * gKD(m,n) ) / 3; } 00319 Real gK2 ( Int m, Int n ) { return ( gK(m,n) - gKD(m,n) ) / 3; } 00320 Real fK1 ( Int m, Int n ) { return ( fK(m,n) + 2 * fKD(m,n) ) / 3; } 00321 Real fK2 ( Int m, Int n ) { return ( fK(m,n) - fKD(m,n) ) / 3; } 00322 #else 00323 Real gK ( Int m, Int n ) { return gK1(m,n) + 2 * gK2(m,n); } 00324 Real gK_r ( Int m, Int n ) { return gK1_r(m,n) + 2 * gK2_r(m,n); } 00325 Real fK ( Int m, Int n ) { return fK1(m,n) + 2 * fK2(m,n); } 00326 Real fK_r ( Int m, Int n ) { return fK1_r(m,n) + 2 * fK2_r(m,n); } 00327 #endif 00328 00329 // The following macros point to equations and not to the grid functions 00330 /// @todo gBet, fBet to be grid functions (after calc of gAlp, fAlp) 00331 00332 #define gBet(m,n) eq_gBet(m,n) 00333 #define fBet(m,n) eq_fBet(m,n) 00334 #define gBet_r(m,n) eq_gBet_r(m,n) 00335 #define fBet_r(m,n) eq_fBet_r(m,n) 00336 00337 ///////////////////////////////////////////////////////////////////////////////////// 00338 00339 #define eq_pr(m,n) (p(m,n)/r(m,n)) 00340 #define eq_qr(m,n) (q(m,n)/r(m,n)) 00341 00342 emitDerivative_r(eq_pr) 00343 emitDerivative_r(eq_qr) 00344 00345 ///////////////////////////////////////////////////////////////////////////////////// 00346 // The evolution equations for gAlp and gDAlp 00347 // Here: Bona-Masso slicing condition with f(alp) = 2/alp 00348 00349 inline Real eq_BM_gAlp_t( Int m, Int n ) { 00350 return -2 * gAlp(m,n) * gK(m,n); 00351 } 00352 00353 inline Real eq_BM_gDAlp_t( Int m, Int n ) { 00354 return -2 * gK_r(m,n); 00355 } 00356 00357 ///////////////////////////////////////////////////////////////////////////////////// 00358 // The time derivatives (these are the evolution equations) 00359 ///////////////////////////////////////////////////////////////////////////////////// 00360 00361 Real eq_gA_t ( Int m, Int n ); Real eq_fA_t ( Int m, Int n ); 00362 Real eq_gB_t ( Int m, Int n ); Real eq_fB_t ( Int m, Int n ); 00363 Real eq_gDA_t ( Int m, Int n ); Real eq_fDA_t ( Int m, Int n ); 00364 Real eq_gDB_t ( Int m, Int n ); Real eq_fDB_t ( Int m, Int n ); 00365 Real eq_gSig_t ( Int m, Int n ); Real eq_fSig_t ( Int m, Int n ); 00366 00367 ///////////////////////////////////////////////////////////////////////////////////// 00368 // The split form of operators for the evolution of gK, gKD, fK, and fKD 00369 ///////////////////////////////////////////////////////////////////////////////////// 00370 00371 Real eq_base_gK_t( Int m, Int n ); 00372 Real eq_invr_gK_t( Int m, Int n ); 00373 00374 Real eq_gK_t( Int m, Int n ) { 00375 return eq_base_gK_t(m,n) + eq_invr_gK_t(m,n) / r(m,n); 00376 } 00377 00378 Real eq_base_gKD_t( Int m, Int n ); 00379 Real eq_invr_gKD_t( Int m, Int n ); 00380 00381 Real eq_gKD_t( Int m, Int n ) { 00382 return eq_base_gKD_t(m,n) + eq_invr_gKD_t(m,n) / r(m,n); 00383 } 00384 00385 Real eq_base_fK_t( Int m, Int n ); 00386 Real eq_invr_fK_t( Int m, Int n ); 00387 00388 Real eq_fK_t( Int m, Int n ) { 00389 return eq_base_fK_t(m,n) + eq_invr_fK_t(m,n) / r(m,n); 00390 } 00391 00392 Real eq_base_fKD_t( Int m, Int n ); 00393 Real eq_invr_fKD_t( Int m, Int n ); 00394 00395 Real eq_fKD_t( Int m, Int n ) { 00396 return eq_base_fKD_t(m,n) + eq_invr_fKD_t(m,n) / r(m,n); 00397 } 00398 00399 // Here defined for completeness (though never used): 00400 // 00401 Real eq_base_gK1_t( Int m, Int n ); 00402 Real eq_invr_gK1_t( Int m, Int n ); 00403 Real eq_base_gK2_t( Int m, Int n ); 00404 Real eq_invr_gK2_t( Int m, Int n ); 00405 Real eq_base_fK1_t( Int m, Int n ); 00406 Real eq_invr_fK1_t( Int m, Int n ); 00407 Real eq_base_fK2_t( Int m, Int n ); 00408 Real eq_invr_fK2_t( Int m, Int n ); 00409 00410 ///////////////////////////////////////////////////////////////////////////////////// 00411 // The p-equations 00412 ///////////////////////////////////////////////////////////////////////////////////// 00413 00414 Real eq_p_t ( Int m, Int n ); 00415 Real eq_p1_t ( Int m, Int n ); 00416 00417 Real eq_base_p_g ( Int m, Int n ); 00418 Real eq_invr_p_g ( Int m, Int n ); 00419 00420 Real eq_p_g( Int m, Int n ) { 00421 return eq_base_p_g(m,n) + eq_invr_p_g(m,n) / r(m,n); 00422 } 00423 00424 Real eq_base_p_f ( Int m, Int n ); 00425 Real eq_invr_p_f ( Int m, Int n ); 00426 00427 Real eq_p_f( Int m, Int n ) { 00428 return eq_base_p_f(m,n) + eq_invr_p_f(m,n) / r(m,n); 00429 } 00430 00431 ///////////////////////////////////////////////////////////////////////////////////// 00432 // The constraints (here used as the error estimators) and the equation for `p` 00433 ///////////////////////////////////////////////////////////////////////////////////// 00434 00435 Real eq_gC_rho ( Int m, Int n ); 00436 Real eq_fC_rho ( Int m, Int n ); 00437 Real eq_gC_j ( Int m, Int n ); 00438 Real eq_fC_j ( Int m, Int n ); 00439 Real eq_C_bimConsLaw ( Int m, Int n ); 00440 Real eq_j_over_p ( Int m, Int n ); 00441 00442 Real eq_C_1 ( Int m, Int n ); 00443 Real eq_C_2 ( Int m, Int n ); 00444 Real eq_C_3 ( Int m, Int n ); 00445 Real eq_C_4 ( Int m, Int n ); 00446 00447 ///////////////////////////////////////////////////////////////////////////////////// 00448 // The sources 00449 ///////////////////////////////////////////////////////////////////////////////////// 00450 00451 Real eq_g_rho ( Int m, Int n ); Real eq_f_rho ( Int m, Int n ); 00452 Real eq_g_j ( Int m, Int n ); Real eq_f_j ( Int m, Int n ); 00453 Real eq_g_JK1 ( Int m, Int n ); Real eq_f_JK1 ( Int m, Int n ); 00454 Real eq_g_JK2 ( Int m, Int n ); Real eq_f_JK2 ( Int m, Int n ); 00455 00456 ///////////////////////////////////////////////////////////////////////////////////// 00457 // The perfect fluid 00458 ///////////////////////////////////////////////////////////////////////////////////// 00459 00460 Real eq_pf_D ( Int m, Int n ); 00461 Real eq_pf_S ( Int m, Int n ); 00462 Real eq_pf_tau ( Int m, Int n ); 00463 00464 Real eq_pf_v ( Int m, Int n ); 00465 Real eq_pf_v_r ( Int m, Int n ); 00466 00467 Real eq_pf_D_t ( Int m, Int n ); 00468 Real eq_pf_S_t ( Int m, Int n ); 00469 Real eq_pf_tau_t ( Int m, Int n ); 00470 00471 Real eq_pf_rho ( Int m, Int n ); 00472 Real eq_pf_j ( Int m, Int n ); 00473 Real eq_pf_J1 ( Int m, Int n ); 00474 Real eq_pf_J2 ( Int m, Int n ); 00475 00476 ///////////////////////////////////////////////////////////////////////////////////// 00477 // The gauge grid-functions 00478 ///////////////////////////////////////////////////////////////////////////////////// 00479 00480 Real eq_q ( Int m, Int n ); 00481 Real eq_gAlp ( Int m, Int n ); 00482 Real eq_gW ( Int m, Int n ); 00483 Real eq_fW ( Int m, Int n ); 00484 Real eq_cW ( Int m, Int n ); 00485 Real eq_fAlp ( Int m, Int n ); 00486 00487 ///////////////////////////////////////////////////////////////////////////////////// 00488 // The apparent horizon finders 00489 ///////////////////////////////////////////////////////////////////////////////////// 00490 00491 Real eq_gHorz ( Int m, Int n ); 00492 Real eq_fHorz ( Int m, Int n ); 00493 00494 ///////////////////////////////////////////////////////////////////////////////////// 00495 // The shift equations 00496 ///////////////////////////////////////////////////////////////////////////////////// 00497 00498 Real eq_gBet ( Int m, Int n ); 00499 Real eq_gBet_r ( Int m, Int n ); 00500 Real eq_gBetr_r ( Int m, Int n ); 00501 Real eq_fBet ( Int m, Int n ); 00502 Real eq_fBet_r ( Int m, Int n ); 00503 Real eq_fBetr_r ( Int m, Int n ); 00504 /*@}*/ 00505 ///////////////////////////////////////////////////////////////////////////////////// 00506 /** @defgroup g7 Maximal slicing */ 00507 ///////////////////////////////////////////////////////////////////////////////////// 00508 /*@{*/ 00509 void maximalSlice_2( Int m, Int N, Real gAlp_at_N ); 00510 void maximalSlice_4( Int m, Int N, Real gAlp_at_N ); 00511 void maximalSlice_2opt( Int m, Int N, Real gAlp_at_N ); 00512 void maximalSlice_PostSteps( Int m, Int N ); 00513 /*@}*/ 00514 ///////////////////////////////////////////////////////////////////////////////////// 00515 /** @defgroup g8 Integration */ 00516 ///////////////////////////////////////////////////////////////////////////////////// 00517 /*@{*/ 00518 /** Calculate the derived variables R, Lt, and pfv. 00519 * @note TINY_Real is added to the denominator of pfv to avoid dividing by zero. 00520 */ 00521 void calculateDerivedVariables( Int m, Int n ) 00522 { 00523 R(m,n) = fB(m,n) / gB(m,n); 00524 Lt(m,n) = sqrt( 1 + p(m,n) * p(m,n) ); 00525 pfv(m,n) = pfS(m,n) == 0 ? 0 : eq_pf_v(m,n); 00526 } 00527 00528 /** Determine gAlp (e.g. by maximal slicing). Then find fAlp and also calculate 00529 * the derivatives of the lapse related grid-functions. 00530 */ 00531 void determineGaugeFunctions( Int m ); 00532 00533 /** Sommerfeld outgoing (radiative) wave condition on a variable. 00534 */ 00535 void applySommerfeldBC 00536 ( 00537 Int m, //!< The time slice 00538 Int n, //!< The spatial position 00539 Int m_old, //!< The earlier time 00540 Int gf, //!< A grid function index 00541 Real background = 0.0, //!< Background value; 1.0 for ALPHA 00542 Int fall_off = 1 //!< Fall-off power exponent; 2 for VET, BET 00543 ); 00544 00545 ///////////////////////////////////////////////////////////////////////////////////// 00546 // IntegFace implementation 00547 ///////////////////////////////////////////////////////////////////////////////////// 00548 00549 /** Calculate the dependent variables from the prime state variables 00550 * which are needed for the integration. 00551 */ 00552 virtual void integStep_Prepare( Int m ); 00553 00554 /** Prepare the right-hand side of the evolution equations. 00555 */ 00556 virtual void integStep_CalcEvolutionRHS( Int m ); 00557 00558 /** Perform the additional steps after each integration step. 00559 */ 00560 virtual void integStep_Finalize( Int mNext, Int mPrev ); 00561 00562 /** At each checkpoint calculate diagnostics variables and check for eventual NaNs. 00563 * Returns 'false' as the indicator to abort the integration. 00564 */ 00565 virtual bool integStep_Diagnostics( Int m, Int chkNaNs_nFrom, Int chkNaNs_nTo ); 00566 00567 /** Fix the state variables at the left boundary. 00568 */ 00569 virtual void applyLeftBoundaryCondition( Int m ); 00570 00571 /** Fix the state variables at the right boundary. 00572 */ 00573 virtual void applyRightBoundaryCondition( Int m ); 00574 /*@}*/ 00575 ///////////////////////////////////////////////////////////////////////////////////// 00576 /** @defgroup g9 Public methods */ 00577 ///////////////////////////////////////////////////////////////////////////////////// 00578 /*@{*/ 00579 public: 00580 00581 /** Creates and configures the bimetric solver from the given parameters. 00582 */ 00583 BimetricEvolve( Parameters& params, 00584 UniformGrid& ug, GridOutputWriter& output, MoL& integ ); 00585 }; 00586 00587 ///////////////////////////////////////////////////////////////////////////////////////// 00588 // The constructor 00589 ///////////////////////////////////////////////////////////////////////////////////////// 00590 00591 BimetricEvolve::BimetricEvolve( Parameters& params, 00592 UniformGrid& ug, GridOutputWriter& output, MoL& integ ) 00593 : BimetricModel( params ), GridUser( ug ) 00594 { 00595 static std::map<std::string,int> knownSlicings = 00596 { 00597 { "const", SLICE_CONSTG }, 00598 { "constg", SLICE_CONSTG }, { "constgf", SLICE_CONSTGF }, 00599 { "MS2OPT", SLICE_MS2OPT }, { "MS2", SLICE_MS2 }, 00600 { "MS4", SLICE_MS4 }, { "BM", SLICE_BM } 00601 }; 00602 std::string name = params.get( "slicing.method", slicing, 0, knownSlicings ); 00603 00604 params.get( "slicing.lin2n", lin2n, nGhost ); 00605 params.get( "slicing.cub2n", cub2n, 5 * nGhost / 2 + 6 ); 00606 params.get( "slicing.smooth", smooth, 0 ); 00607 00608 slog << "Bimetric Solver:" << std::endl << std::endl 00609 << " slicing = " << name << " (#" << slicing << ")" 00610 << ", lin2n = " << lin2n << ", cub2n = " << cub2n 00611 << ", smooth = " << smooth 00612 << std::endl << std::endl; 00613 00614 if ( mpiSize() > 1 && 00615 ( slicing == SLICE_MS2 || slicing == SLICE_MS2OPT || slicing == SLICE_MS4 ) ) 00616 { 00617 slog << "*** Error: Maximal slicing is not compatible with MPI." << std::endl; 00618 gridDriver->quit( -1 ); 00619 } 00620 00621 // Cached from the integrator 00622 // 00623 delta_t = integ.dt (); 00624 00625 // Sign up for the integration 00626 // 00627 integ.addToEvolution( this ); 00628 00629 // Add our grid functions to the evolution 00630 // 00631 integ.keepConstant( { fld::q } ); // GFs that are kept constant in time 00632 integ.keepEvolved( fld::bimEvolvedGF ); // GFs that are evolved by the integrator 00633 00634 if( isGR () || slicing == SLICE_CONSTGF ) { 00635 integ.keepConstant( { fld::fAlp } ); 00636 } 00637 00638 if ( slicing == SLICE_CONSTG || slicing == SLICE_CONSTGF ) { 00639 integ.keepConstant( { fld::gAlp, fld::gDAlp } ); 00640 } 00641 else if ( slicing == SLICE_BM ) 00642 { 00643 const static std::vector<fld::EvolvedBy> evolvedGaugeGF = { 00644 { fld::gAlp, fld::gAlp_t }, 00645 { fld::gDAlp, fld::gDAlp_t } 00646 }; 00647 integ.keepEvolved( evolvedGaugeGF ); 00648 } 00649 00650 // The list of the grid functions to be written to the output. 00651 // 00652 output.gridFunctions( fld::bimOutput ); 00653 00654 bool showDiagnostics = false; 00655 params.get( "out.diagnostics", showDiagnostics, false ); 00656 00657 if( showDiagnostics ) { 00658 output.gridFunctions( fld::bimShowDiagnostics ); 00659 } 00660 } 00661 00662 ///////////////////////////////////////////////////////////////////////////////////////// 00663 // Implementation of the integration interface (BimetricEvolve::IntegFace methods) 00664 ///////////////////////////////////////////////////////////////////////////////////////// 00665 00666 void BimetricEvolve::applyLeftBoundaryCondition( Int m ) 00667 { 00668 for( Int i = 0; i < nGhost; ++i ) 00669 { 00670 Int n = nGhost - i - 1; 00671 Int nR = nGhost + i; 00672 00673 /// Here we impose the parity conditions at the left boundary. 00674 00675 t (m,n) = t (m,nR); 00676 r (m,n) = -r (m,nR); 00677 00678 p (m,n) = -p (m,nR); 00679 00680 gA (m,n) = gA (m,nR); fA (m,n) = fA (m,nR); 00681 gB (m,n) = gB (m,nR); fB (m,n) = fB (m,nR); 00682 gK (m,n) = gK (m,nR); fK (m,n) = fK (m,nR); 00683 gKD (m,n) = gKD (m,nR); fKD (m,n) = fKD (m,nR); 00684 gDA (m,n) = -gDA (m,nR); fDA (m,n) = -fDA (m,nR); 00685 gDB (m,n) = -gDB (m,nR); fDB (m,n) = -fDB (m,nR); 00686 gSig (m,n) = -gSig (m,nR); fSig (m,n) = -fSig (m,nR); 00687 00688 pfD (m,n) = pfD (m,nR); 00689 pfS (m,n) = -pfS (m,nR); 00690 pftau (m,n) = pftau (m,nR); 00691 00692 q (m,n) = -q (m,nR); 00693 gAlp (m,n) = gAlp (m,nR); 00694 fAlp (m,n) = fAlp (m,nR); 00695 gDAlp (m,n) = -gDAlp (m,nR); 00696 00697 Lt (m,n) = Lt (m,nR); 00698 R (m,n) = R (m,nR); 00699 pfv (m,n) = -pfv (m,nR); 00700 } 00701 } 00702 00703 void BimetricEvolve::applyRightBoundaryCondition( Int m ) 00704 { 00705 for( Int n = nGhost + nLen; n < nTotal; ++n ) 00706 { 00707 t(m,n) = t(m,n-1); 00708 r(m,n) = r(m,n-1) + delta_r; 00709 00710 extrapolate_R( fld::p, m, n ); // Extrapolate [n] from [n-1], [n-2], ... 00711 00712 extrapolate_R( fld::gA, m, n ); extrapolate_R( fld::fA, m, n ); 00713 extrapolate_R( fld::gB, m, n ); extrapolate_R( fld::fB, m, n ); 00714 extrapolate_R( fld::gK, m, n ); extrapolate_R( fld::fK, m, n ); 00715 extrapolate_R( fld::gKD, m, n ); extrapolate_R( fld::fKD, m, n ); 00716 extrapolate_R( fld::gDA, m, n ); extrapolate_R( fld::fDA, m, n ); 00717 extrapolate_R( fld::gDB, m, n ); extrapolate_R( fld::fDB, m, n ); 00718 extrapolate_R( fld::gSig, m, n ); extrapolate_R( fld::fSig, m, n ); 00719 00720 extrapolate_R( fld::pfD, m, n ); 00721 extrapolate_R( fld::pfS, m, n ); 00722 extrapolate_R( fld::pftau, m, n ); 00723 00724 extrapolate_R( fld::q, m, n ); 00725 extrapolate_R( fld::gAlp, m, n ); 00726 extrapolate_R( fld::gDAlp, m, n ); 00727 extrapolate_R( fld::fAlp, m, n ); 00728 00729 calculateDerivedVariables( m, n ); // Calculate R, Lt, and pfv. 00730 00731 extrapolate_R( fld::pfv, m, n ); // Nevertheless, we extrapolate pfv! 00732 } 00733 } 00734 00735 void BimetricEvolve::applySommerfeldBC 00736 ( 00737 Int m1, Int n, Int m, Int gf, 00738 const Real background, 00739 const Int fall_off 00740 ) 00741 { 00742 const Real wave_speed = 1.0; 00743 const Real Courant = wave_speed * delta_t / delta_r; 00744 00745 for( Int n = nGhost + nLen + 1; n < nTotal - 1; ++n ) 00746 { 00747 const Real r0 = r( m, n-2 ); 00748 const Real r1 = r( m, n-1 ); 00749 const Real r2 = r( m, n ); 00750 const Real r_last = Courant * r1 + ( 1.0 - Courant ) * r2; 00751 00752 Real factor = r_last / r2; 00753 if( fall_off == 2 ) { 00754 factor *= factor; 00755 } 00756 00757 const Real P0 = GF( gf, m, n-2 ); 00758 const Real P1 = GF( gf, m, n-1 ); 00759 const Real P2 = GF( gf, m, n ); 00760 const Real P01 = ( ( r_last - r1 ) * P0 + ( r0 - r_last ) * P1 ) / ( r0 - r1 ); 00761 const Real P12 = ( ( r_last - r2 ) * P1 + ( r1 - r_last ) * P2 ) / ( r1 - r2 ); 00762 const Real P012 = ( ( r_last - r2 ) * P01 + ( r0 - r_last ) * P12 ) / ( r0 - r2 ); 00763 00764 GF( gf, m1, n ) = factor * ( P012 - background ) + background; 00765 } 00766 } 00767 00768 void BimetricEvolve::integStep_Prepare( Int m ) 00769 { 00770 // Move gKD and fKD to 0 at r = 0 00771 // Real dgKD = gKD(m,nGhost)/2, dfKD = fKD(m,nGhost)/2; 00772 00773 OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n ) 00774 { 00775 calculateDerivedVariables( m, n ); 00776 // gKD(m,n) -= dgKD; 00777 // fKD(m,n) -= dfKD; 00778 } 00779 } 00780 00781 void BimetricEvolve::determineGaugeFunctions( Int m ) 00782 { 00783 ///////////////////////////////////////////////////////////////////////////////////// 00784 /// - Maximal slicing (optional) -- this will calculate gAlp and gDAlp 00785 00786 Int sliceBC = smooth >= 1 ? round(40/delta_r) : nLen/2; 00787 switch( slicing ) 00788 { 00789 case SLICE_MS2OPT: maximalSlice_2opt ( m, sliceBC, 1 ); break; 00790 case SLICE_MS2: maximalSlice_2 ( m, sliceBC, 1 ); break; 00791 case SLICE_MS4: maximalSlice_4 ( m, sliceBC, 1 ); break; 00792 } 00793 00794 #if 1 00795 #define KDT 00796 #else 00797 #define KDT 0* 00798 if( true ) 00799 { 00800 OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n ) 00801 { 00802 gW ( m, n ) = eq_gW ( m, n ); 00803 fW ( m, n ) = eq_fW ( m, n ); 00804 cW ( m, n ) = eq_cW ( m, n ); 00805 fAlp ( m, n ) = eq_fAlp ( m, n ); 00806 00807 gDAlp_r( m, n ) = 0; 00808 fAlp_r ( m, n ) = 0; 00809 fAlp_rr( m, n ) = 0; 00810 fDAlp ( m, n ) = 0; 00811 fDAlp_r( m, n ) = 0; 00812 } 00813 } else 00814 #endif 00815 00816 /// @todo Go through the order of evaluation for the gauge grid functions 00817 /// @todo GF_r and GF_rr requires fixed ghost zones 00818 /// 00819 if ( isGR () || slicing == SLICE_CONSTGF ) 00820 { 00821 OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n ) 00822 { 00823 gAlp_r ( m, n ) = GF_r( gAlp, m, n ); 00824 gDAlp_r( m, n ) = GF_r( gDAlp, m, n ); 00825 fAlp_r ( m, n ) = GF_r( fAlp, m, n ); 00826 fAlp_rr( m, n ) = GF_rr( fAlp, m, n ); 00827 fDAlp ( m, n ) = fAlp_r(m,n) / fAlp(m,n); 00828 fDAlp_r( m, n ) = ( fAlp_rr(m,n) - fAlp_r(m,n) * fAlp_r(m,n) / fAlp(m,n) ) 00829 / fAlp(m,n); 00830 } 00831 } 00832 else // if not GR 00833 { 00834 // gAlp and gDAlp must the BC fixed at this point 00835 // 00836 OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n ) { 00837 gAlp_r ( m, n ) = GF_r( gAlp, m, n ); 00838 gDAlp_r( m, n ) = GF_r( gDAlp, m, n ); 00839 } 00840 00841 OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n ) { 00842 cW ( m, n ) = eq_cW ( m, n ); 00843 gW ( m, n ) = eq_gW ( m, n ); 00844 fW ( m, n ) = eq_fW ( m, n ); 00845 fAlp ( m, n ) = eq_fAlp ( m, n ); 00846 } 00847 00848 if( smooth ) { 00849 smoothenGF( m, fld::fAlp, fld::tmp, fld::fAlp, 1 ); 00850 } 00851 else { 00852 applyBoundaryConditions( m, fld::fAlp, +1 ); 00853 } 00854 00855 OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n ) { 00856 fAlp_r( m, n ) = GF_r( fAlp, m, n ); 00857 } 00858 00859 if( smooth ) { 00860 smoothenGF( m, fld::fAlp_r, fld::tmp, fld::fAlp_r, -1 ); 00861 } 00862 else { 00863 applyBoundaryConditions( m, fld::fAlp_r, -1 ); 00864 } 00865 00866 OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n ) { 00867 fAlp_rr( m, n ) = GF_r( fAlp_r, m, n ); 00868 } 00869 00870 if( smooth ) { 00871 smoothenGF( m, fld::fAlp_rr, fld::tmp, fld::fAlp_rr, 1 ); 00872 } 00873 else { 00874 applyBoundaryConditions( m, fld::fAlp_r, -1 ); 00875 } 00876 00877 OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n ) 00878 { 00879 fDAlp ( m, n ) = fAlp_r(m,n) / fAlp(m,n); 00880 fDAlp_r( m, n ) = ( fAlp_rr(m,n) - fAlp_r(m,n) * fAlp_r(m,n) / fAlp(m,n) ) 00881 / fAlp(m,n); 00882 } 00883 } 00884 } 00885 00886 void BimetricEvolve::integStep_CalcEvolutionRHS( Int m ) 00887 { 00888 ///////////////////////////////////////////////////////////////////////////////////// 00889 /// - First, calculate the values of the spatial derivatives 00890 /// @todo fixme: Severely lost accuracy in all second derivatives 00891 00892 OMP_parallel_for( Int n = CFDS_ORDER/2; n < nTotal - CFDS_ORDER/2; ++n ) 00893 { 00894 gA_r (m,n) = GF_r( gA, m,n ); fA_r (m,n) = GF_r( fA, m,n ); 00895 gB_r (m,n) = GF_r( gB, m,n ); fB_r (m,n) = GF_r( fB, m,n ); 00896 gK_r (m,n) = GF_r( gK, m,n ); fK_r (m,n) = GF_r( fK, m,n ); 00897 gKD_r (m,n) = GF_r( gKD, m,n ); fKD_r (m,n) = GF_r( fKD, m,n ); 00898 gDA_r (m,n) = GF_r( gDA, m,n ); fDA_r (m,n) = GF_r( fDA, m,n ); 00899 gDB_r (m,n) = GF_r( gDB, m,n ); fDB_r (m,n) = GF_r( fDB, m,n ); 00900 00901 p_r (m,n) = GF_r( p, m,n ); q_r (m,n) = GF_r( q, m,n ); 00902 p_rr (m,n) = GF_rr( p, m,n ); q_rr (m,n) = GF_rr( q, m,n ); 00903 eq_pr_r(m,n) = GF_r( eq_pr, m,n ); eq_qr_r(m,n) = GF_r( eq_qr, m,n ); 00904 00905 pfD_r (m,n) = GF_r( pfD, m,n ); 00906 pfS_r (m,n) = GF_r( pfS, m,n ); 00907 pftau_r(m,n) = GF_r( pftau, m,n ); 00908 00909 // pfv_r(m,n) should be last as it depends on pfD_r, pfS_r, and pftau_r 00910 // 00911 pfv_r(m,n) = pfS(m,n) == 0 ? 0 : eq_pf_v_r( m, n ); 00912 } 00913 00914 if( smooth >= 2 ) 00915 { 00916 smoothenGF( m, fld::p_r, fld::tmp, fld::p_r, +1 ); 00917 smoothenGF( m, fld::p_rr, fld::tmp, fld::p_rr, -1 ); 00918 smoothenGF( m, fld::eq_pr_r, fld::tmp, fld::eq_pr_r, +1 ); 00919 00920 smoothenGF( m, fld::gA_r, fld::tmp, fld::gA_r, -1 ); 00921 smoothenGF( m, fld::gB_r, fld::tmp, fld::gB_r, -1 ); 00922 smoothenGF( m, fld::gK_r, fld::tmp, fld::gK_r, -1 ); 00923 smoothenGF( m, fld::gKD_r, fld::tmp, fld::gKD_r, -1 ); 00924 smoothenGF( m, fld::gDA_r, fld::tmp, fld::gDA_r, +1 ); 00925 smoothenGF( m, fld::gDB_r, fld::tmp, fld::gDB_r, +1 ); 00926 00927 smoothenGF( m, fld::fA_r, fld::tmp, fld::fA_r, -1 ); 00928 smoothenGF( m, fld::fB_r, fld::tmp, fld::fB_r, -1 ); 00929 smoothenGF( m, fld::fK_r, fld::tmp, fld::fK_r, -1 ); 00930 smoothenGF( m, fld::fKD_r, fld::tmp, fld::fKD_r, -1 ); 00931 smoothenGF( m, fld::fDA_r, fld::tmp, fld::fDA_r, +1 ); 00932 smoothenGF( m, fld::fDB_r, fld::tmp, fld::fDB_r, +1 ); 00933 } 00934 00935 /// - Calculate the gauge (e.g., do maximal slicing) 00936 /// 00937 determineGaugeFunctions( m ); 00938 00939 ///////////////////////////////////////////////////////////////////////////////////// 00940 /// - Calculate the intermediate variables that do not depend on derivatives. 00941 00942 OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n ) 00943 { 00944 // The sources (dependent on both p and the primary dynamical fields) 00945 // 00946 g_rho (m,n) = eq_g_rho (m,n); f_rho (m,n) = eq_f_rho (m,n); 00947 g_JK1 (m,n) = eq_g_JK1 (m,n); f_JK1 (m,n) = eq_f_JK1 (m,n); 00948 g_JK2 (m,n) = eq_g_JK2 (m,n); f_JK2 (m,n) = eq_f_JK2 (m,n); 00949 } 00950 00951 ///////////////////////////////////////////////////////////////////////////////////// 00952 /// - Calculate the variables that depend on the spatial derivatives 00953 00954 OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n ) 00955 { 00956 // The time evolution right-hand side 00957 // 00958 p_t (m,n) = isGR() ? 0 : KDT eq_p_t(m,n); // + eq_C_bimConsLaw (m,n); 00959 00960 gA_t (m,n) = eq_gA_t (m,n); fA_t (m,n) = eq_fA_t (m,n); 00961 gB_t (m,n) = eq_gB_t (m,n); fB_t (m,n) = eq_fB_t (m,n); 00962 gK_t (m,n) = eq_gK_t (m,n); fK_t (m,n) = eq_fK_t (m,n); 00963 gKD_t (m,n) = KDT eq_gKD_t (m,n); fKD_t (m,n) = KDT eq_fKD_t (m,n); 00964 gDA_t (m,n) = KDT eq_gDA_t (m,n); fDA_t (m,n) = KDT eq_fDA_t (m,n); 00965 gDB_t (m,n) = KDT eq_gDB_t (m,n); fDB_t (m,n) = KDT eq_fDB_t (m,n); 00966 gSig_t (m,n) = KDT eq_gSig_t (m,n); fSig_t (m,n) = KDT eq_fSig_t (m,n); 00967 00968 pfD_t (m,n) = eq_pf_D_t (m,n); 00969 pfS_t (m,n) = eq_pf_S_t (m,n); 00970 pftau_t (m,n) = eq_pf_tau_t(m,n); 00971 00972 if( slicing == SLICE_BM ) 00973 { 00974 gAlp_t (m,n) = eq_BM_gAlp_t (m,n); 00975 gDAlp_t (m,n) = eq_BM_gDAlp_t (m,n); 00976 } 00977 } 00978 00979 if( smooth >= 2 ) 00980 { 00981 // smoothenGF2( m, fld::p_t, fld::tmp, fld::p_t, -1 ); 00982 // smoothenGF( m, fld::gA_t, fld::tmp, fld::gA_t, +1 ); 00983 // smoothenGF( m, fld::gB_t, fld::tmp, fld::gB_t, +1 ); 00984 smoothenGF2( m, fld::gK_t, fld::tmp, fld::gK_t, +1 ); 00985 smoothenGF2( m, fld::gKD_t, fld::tmp, fld::gKD_t, +1 ); 00986 // smoothenGF2( m, fld::gDA_t, fld::tmp, fld::gDA_t, -1 ); 00987 // smoothenGF2( m, fld::gDB_t, fld::tmp, fld::gDB_t, -1 ); 00988 // smoothenGF( m, fld::fA_t, fld::tmp, fld::fA_t, +1 ); 00989 // smoothenGF( m, fld::fB_t, fld::tmp, fld::fB_t, +1 ); 00990 smoothenGF2( m, fld::fK_t, fld::tmp, fld::fK_t, +1 ); 00991 smoothenGF2( m, fld::fKD_t, fld::tmp, fld::fKD_t, +1 ); 00992 // smoothenGF2( m, fld::fDA_t, fld::tmp, fld::fDA_t, -1 ); 00993 // smoothenGF2( m, fld::fDB_t, fld::tmp, fld::fDB_t, -1 ); 00994 } 00995 00996 ///////////////////////////////////////////////////////////////////////////////////// 00997 /// - Smoothen the time derivatives inside the grid zone near the outer boundary 00998 /// @todo Smoothen the time derivatives inside the left grid zone 00999 01000 if ( gridDriver->isLastInRank() ) 01001 { 01002 for( Int n = nGhost + nLen - CFDS_ORDER; n < nGhost + nLen; ++n ) 01003 { 01004 for( auto e: fld::bimEvolvedGF ) { 01005 extrapolate_lin( e.f_t, m, n ); 01006 } 01007 } 01008 } 01009 } 01010 01011 ///////////////////////////////////////////////////////////////////////////////////////// 01012 /// Post-evolution smoothing and calculation of diagnostics 01013 /// 01014 void BimetricEvolve::integStep_Finalize( Int m1, Int m ) 01015 { 01016 } 01017 01018 ///////////////////////////////////////////////////////////////////////////////////////// 01019 /// Find horizons, calculate constraints and similar. Also check for NaNs. 01020 /// 01021 bool BimetricEvolve::integStep_Diagnostics( Int m, Int chkNaNs_nFrom, Int chkNaNs_nTo ) 01022 { 01023 OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n ) 01024 { 01025 /// - Apparent horizon finder (dependent only on the primary dynamical fields) 01026 /// 01027 gHorz (m,n) = eq_gHorz (m,n); 01028 fHorz (m,n) = eq_fHorz (m,n); 01029 01030 /// - The constraint violations 01031 /// 01032 C_1 (m,n) = eq_C_1 (m,n); 01033 C_2 (m,n) = eq_C_2 (m,n); 01034 C_3 (m,n) = eq_C_3 (m,n); 01035 C_4 (m,n) = eq_C_bimConsLaw (m,n); 01036 01037 // C_1 (m,n) = abs( eq_C_1 (m,n) / (1e-10 + g_rho (m,n) ) ); 01038 // C_2 (m,n) = abs( eq_C_2 (m,n) / (1e-10 + f_rho (m,n) ) ); 01039 // C_3 (m,n) = eq_C_3 (m,n); 01040 // C_4 (m,n) = abs( eq_C_bimConsLaw (m,n) / (1e-10 + p_t (m,n) ) ); 01041 01042 // C_1 (m,n) = abs( eq_C_1 (m,n) ); 01043 // C_2 (m,n) = abs( eq_C_2 (m,n) ); 01044 // C_3 (m,n) = abs( eq_C_3 (m,n) ); 01045 // C_4 (m,n) = abs( eq_C_bimConsLaw (m,n) ); 01046 01047 /// - Calculate p in an alternative way using the algebraic equations 01048 /// 01049 p_g (m,n) = eq_p_g (m,n); 01050 p_f (m,n) = eq_p_f (m,n); 01051 } 01052 01053 if ( chkNaNs_nTo < 0 ) { // CheckNaNs is disabled 01054 return true; 01055 } 01056 01057 /// - Check for NaNs in gA in the given zone. 01058 /// 01059 Int nFrom = nGhost + std::min( nLen, std::max( Int(0), chkNaNs_nFrom - nOffset ) ); 01060 Int nTo = nGhost + std::min( nLen, std::max( Int(0), chkNaNs_nTo - nOffset ) ); 01061 01062 // if( false ) 01063 for( Int n = nFrom; n < nTo; ++n ) 01064 { 01065 if( std::isnan( gA( m, n ) ) ) { 01066 std::cerr << "*** Detected gA NaN at t = " << t(m,n) 01067 << ", r = " << r(m,n) << std::endl; 01068 return false; 01069 } 01070 } 01071 return true; 01072 } 01073 01074 ///////////////////////////////////////////////////////////////////////////////////////// 01075 // Equations of Motion (generated in Mathematica) 01076 ///////////////////////////////////////////////////////////////////////////////////////// 01077 01078 #include "eom-std.h" 01079 01080 ///////////////////////////////////////////////////////////////////////////////////////// 01081 // Maximal slicing (implements a BVP solver) 01082 ///////////////////////////////////////////////////////////////////////////////////////// 01083 01084 #include "maximalSlice.h" 01085 01086 ///////////////////////////////////////////////////////////////////////////////////////// 01087 /// The main entry point of `bim-solver`. 01088 /// 01089 int main( int argc, char* argv[] ) 01090 { 01091 TrackUsedTime timer; 01092 01093 /// - Read the run-time configuration parameters 01094 /// 01095 Parameters params( argc >= 2 ? argv[1] : "config.ini" ); 01096 01097 /// - Create the grid driver 01098 /// 01099 UniformGrid ugrid( params ); 01100 01101 /// - Read the initial data (populating the grid slice at t_0) 01102 /// 01103 if( ! GridInitialData( params, ugrid ).addGFs( fld::bimInput ).load () ) { 01104 return -1; 01105 } 01106 01107 /// - Create the output sink (for storing the results) 01108 /// 01109 GridOutputWriter output ( params, ugrid ); 01110 if ( ! output.open () ) { 01111 return -1; 01112 } 01113 01114 /// - Setup the MoL integrator on the grid (whose results will go to the output sink) 01115 /// 01116 MoL integrator( params, ugrid, output ); 01117 01118 /// - Create the bimetric model and define the equations of motion (to be integrated) 01119 /// 01120 BimetricEvolve bim( params, ugrid, output, integrator ); 01121 01122 /// - Evolve the equations of motion 01123 /// 01124 return integrator.evolveEquations () ? 0 : -1; 01125 }
1.8.0 
