diff options
Diffstat (limited to 'ML_BSSN/src/ML_BSSN_lapse_evol.cc')
| -rw-r--r-- | ML_BSSN/src/ML_BSSN_lapse_evol.cc | 41 |
1 files changed, 33 insertions, 8 deletions
diff --git a/ML_BSSN/src/ML_BSSN_lapse_evol.cc b/ML_BSSN/src/ML_BSSN_lapse_evol.cc index e688b5a..0388510 100644 --- a/ML_BSSN/src/ML_BSSN_lapse_evol.cc +++ b/ML_BSSN/src/ML_BSSN_lapse_evol.cc @@ -211,22 +211,47 @@ static void ML_BSSN_lapse_evol_Body(const cGH* restrict const cctkGH, const int break; } - PDupwindNthAnti2alpha = 0.0; - PDupwindNthSymm2alpha = 0.0; + PDupwindNthAnti2alpha = ToReal(0.0); + PDupwindNthSymm2alpha = ToReal(0.0); + + CCTK_REAL_VEC alpharhsL; + +#if 1 if (lapse_gh) { - const double e4phi = conformalMethod == 1 ? 1.0 / (phiL * phiL) : exp(4.0 * phiL); - const double det_gamma = pow(e4phi, 3.0); + CCTK_REAL_VEC e4phi = kifthen(conformalMethod == 1, ToReal(1.0) / (phiL * phiL), kexp(ToReal(4.0) * phiL)); + CCTK_REAL_VEC det_gamma = kpow(e4phi, 3.0); - alpharhs[index] = -(alphaL * alphaL) * trKL + gh_eta_l_bar * pow(alphaL, gh_q) * (alphaL * alphaL) * log(pow(det_gamma, 0.5 * gh_p) / alphaL) - + LapseAdvectionCoeff * (beta1L * PDupwindNthAnti1alpha + beta3L * PDupwindNthAnti3alpha + - fabs(beta1L) * PDupwindNthSymm1alpha + fabs(beta3L) * PDupwindNthSymm3alpha); + alpharhsL = -(alphaL * alphaL) * trKL + ToReal(gh_eta_l_bar) * kpow(alphaL, gh_q) * (alphaL * alphaL) * klog(kpow(det_gamma, 0.5 * gh_p) / alphaL) + + ToReal(LapseAdvectionCoeff) * (beta1L * PDupwindNthAnti1alpha + beta3L * PDupwindNthAnti3alpha + + kfabs(beta1L) * PDupwindNthSymm1alpha + kfabs(beta3L) * PDupwindNthSymm3alpha); } else { - alpharhs[index] = -harmonicF * pow(alphaL, harmonicN) * trKL + WFactor * WL + alpharhsL = -ToReal(harmonicF) * kpow(alphaL, harmonicN) * trKL + ToReal(WFactor) * WL + + ToReal(LapseAdvectionCoeff) * (beta1L * PDupwindNthAnti1alpha + beta3L * PDupwindNthAnti3alpha + + kfabs(beta1L) * PDupwindNthSymm1alpha + kfabs(beta3L) * PDupwindNthSymm3alpha); + } +#else + { + const double t = cctkGH->cctk_time; + const double fact_harmonic = alphaL * alphaL; + const double fact_1plog = 2.0 * alphaL; + double fact; + + if (t < 5.0) + fact = fact_harmonic; + else if (t > 5.5) + fact = fact_1plog; + else + fact = fact_harmonic * exp(-36.0 * pow((t - 5.0) / 0.5, 4)) + fact_1plog * (1.0 - exp(-36.0 * pow((t - 5.0) / 0.5, 4))); + + alpharhs[index] = -fact * trKL + LapseAdvectionCoeff * (beta1L * PDupwindNthAnti1alpha + beta3L * PDupwindNthAnti3alpha + fabs(beta1L) * PDupwindNthSymm1alpha + fabs(beta3L) * PDupwindNthSymm3alpha); } +#endif + vec_store_partial_prepare(i,vecimin,vecimax); + vec_store_nta_partial(alpharhs[index],alpharhsL); } CCTK_ENDLOOP3STR(ML_BSSN_lapse_evol); } |