diff --git a/CurvADM/README b/CurvADM/README new file mode 100644 index 0000000..bfe3fa5 --- /dev/null +++ b/CurvADM/README @@ -0,0 +1,9 @@ +Cactus Code Thorn CurvADM +Author(s) : Liwei Ji +Maintainer(s): Liwei Ji +Licence : LGPL +-------------------------------------------------------------------------- + +1. Purpose + +Provide storage for the ADM variables diff --git a/CurvADM/configuration.ccl b/CurvADM/configuration.ccl new file mode 100644 index 0000000..aa4c5a3 --- /dev/null +++ b/CurvADM/configuration.ccl @@ -0,0 +1,3 @@ +# Configuration definition for thorn CurvADM + +REQUIRES Loop diff --git a/CurvADM/interface.ccl b/CurvADM/interface.ccl new file mode 100644 index 0000000..96e88bd --- /dev/null +++ b/CurvADM/interface.ccl @@ -0,0 +1,29 @@ +# Interface definition for thorn CurvADM + +IMPLEMENTS: CurvADM + +USES INCLUDE HEADER: loop_device.hxx + + +PUBLIC: + +CCTK_REAL metric TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { gxx gxy gxz gyy gyz gzz } "ADM 3-metric g_ij" + +CCTK_REAL excurv TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { kxx kxy kxz kyy kyz kzz } "ADM extrinsic curvature K_ij" + +CCTK_REAL lapse TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { alp } "ADM lapse function alpha" + +CCTK_REAL dtlapse TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { dtalp } "Time derivative of ADM lapse function" + +CCTK_REAL shift TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { betax betay betaz } "ADM shift vector beta^i" + +CCTK_REAL dtshift TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { dtbetax dtbetay dtbetaz } "Time derivative of ADM shift vector" + + +# Variables necessary to calculate second time derivatives of the four-metric + +CCTK_REAL dtexcurv TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { dtkxx dtkxy dtkxz dtkyy dtkyz dtkzz } "Time derivative of ADM extrinsic curvature K_ij" + +CCTK_REAL dt2lapse TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { dt2alp } "Second time derivative of ADM lapse function" + +CCTK_REAL dt2shift TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { dt2betax dt2betay dt2betaz } "Second time derivative of ADM shift vector" diff --git a/CurvADM/param.ccl b/CurvADM/param.ccl new file mode 100644 index 0000000..4f7ae98 --- /dev/null +++ b/CurvADM/param.ccl @@ -0,0 +1,47 @@ +# Parameter definitions for thorn CurvADM + +RESTRICTED: + +KEYWORD initial_data "Initial metric and extrinsic curvature datasets" +{ + "none" :: "Do not set up initial conditions" + "Cartesian Minkowski" :: "Minkowski geometry in Cartesian coordinates" + "linear wave" :: "linear wave" +} "Cartesian Minkowski" + +KEYWORD initial_lapse "Initial lapse value" +{ + "none" :: "Do not set up initial conditions" + "one" :: "Uniform lapse" +} "one" + +KEYWORD initial_shift "Initial shift value" +{ + "none" :: "Do not set up initial conditions" + "zero" :: "Shift is zero" +} "zero" + +KEYWORD initial_dtlapse "Initial dtlapse value" +{ + "none" :: "Do not set up initial conditions" + "zero" :: "Dtlapse is zero" +} "zero" + +KEYWORD initial_dtshift "Initial dtshift value" +{ + "none" :: "Do not set up initial conditions" + "zero" :: "Dtshift is zero" +} "zero" + + +PRIVATE: + +CCTK_REAL linear_wave_amplitude "Linear wave amplitude" +{ + 0.0:* :: "" +} 1.0e-8 + +CCTK_REAL linear_wave_wavelength "Linear wave wavelength" +{ + 0.0:* :: "" +} 1.0 diff --git a/CurvADM/schedule.ccl b/CurvADM/schedule.ccl new file mode 100644 index 0000000..70d7236 --- /dev/null +++ b/CurvADM/schedule.ccl @@ -0,0 +1,105 @@ +# Schedule definitions for thorn CurvADM + +if (CCTK_IsThornActive("ODESolvers")) { + + SCHEDULE GROUP CurvADM_InitialData IN ODESolvers_Initial + { + } "Schedule group for calculating ADM initial data" + + SCHEDULE GROUP CurvADM_InitialGauge IN ODESolvers_Initial AFTER CurvADM_InitialData + { + } "Schedule group for the ADM initial gauge condition" + + SCHEDULE GROUP CurvADM_PostInitial IN ODESolvers_Initial AFTER (CurvADM_InitialData CurvADM_InitialGauge) + { + } "Schedule group for modifying the ADM initial data, such as e.g. adding noise" + + SCHEDULE GROUP CurvADM_SetADMVars IN ODESolvers_PostStep + { + } "Set ADM variables in this group" + + SCHEDULE GROUP CurvADM_SetADMRHS IN ODESolvers_PostStep + { + } "Set ADM RHS variables in this group" + +} else { + + SCHEDULE GROUP CurvADM_SetADMVars AT post_recover_variables + { + } "Set ADM variables in this group" + + SCHEDULE GROUP CurvADM_SetADMRHS AT post_recover_variables + { + } "Set ADM RHS variables in this group" + + SCHEDULE GROUP CurvADM_SetADMVars AT postregrid + { + } "Set ADM variables in this group" + + SCHEDULE GROUP CurvADM_SetADMRHS AT postregrid + { + } "Set ADM RHS variables in this group" + + SCHEDULE GROUP CurvADM_SetADMVars AT postrestrict + { + } "Set ADM variables in this group" + + SCHEDULE GROUP CurvADM_SetADMRHS AT postrestrict + { + } "Set ADM RHS variables in this group" + + SCHEDULE GROUP CurvADM_SetADMVars AT poststep + { + } "Set ADM variables in this group" + + SCHEDULE GROUP CurvADM_SetADMRHS AT poststep + { + } "Set ADM RHS variables in this group" + +} + +if (CCTK_EQUALS(initial_data, "Cartesian Minkowski")) { + SCHEDULE CurvADM_initial_data IN CurvADM_InitialData + { + LANG: C + WRITES: metric(everywhere) curv(everywhere) + } "Set up Cartesian Minkowski initial data" +} else if (CCTK_EQUALS(initial_data, "linear wave")) { + SCHEDULE CurvADM_linear_wave IN CurvADM_InitialData + { + LANG: C + WRITES: metric(everywhere) curv(everywhere) + } "Set up linear wave initial data" +} + +if (CCTK_EQUALS(initial_lapse, "one")) { + SCHEDULE CurvADM_initial_lapse IN CurvADM_InitialGauge + { + LANG: C + WRITES: lapse(everywhere) + } "Set lapse to one" +} + +if (CCTK_EQUALS(initial_dtlapse, "zero")) { + SCHEDULE CurvADM_initial_dtlapse IN CurvADM_InitialGauge + { + LANG: C + WRITES: dtlapse(everywhere) + } "Set dtlapse to zero" +} + +if (CCTK_EQUALS(initial_shift, "zero")) { + SCHEDULE CurvADM_initial_shift IN CurvADM_InitialGauge + { + LANG: C + WRITES: shift(everywhere) + } "Set shift to zero" +} + +if (CCTK_EQUALS(initial_dtshift, "zero")) { + SCHEDULE CurvADM_initial_dtshift IN CurvADM_InitialGauge + { + LANG: C + WRITES: dtshift(everywhere) + } "Set dtshift to zero" +} diff --git a/CurvADM/src/adm.cxx b/CurvADM/src/adm.cxx new file mode 100644 index 0000000..2790a36 --- /dev/null +++ b/CurvADM/src/adm.cxx @@ -0,0 +1,80 @@ +#include + +#include +#include +#include + +#include + +namespace CurvADM { +using namespace std; +using namespace Loop; + +extern "C" void CurvADM_initial_data(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_CurvADM_initial_data; + DECLARE_CCTK_PARAMETERS; + + grid.loop_all_device<1, 1, 1>(grid.nghostzones, + [=] CCTK_DEVICE(const PointDesc &p) + CCTK_ATTRIBUTE_ALWAYS_INLINE { + gxx(p.I) = 1; + gxy(p.I) = 0; + gxz(p.I) = 0; + gyy(p.I) = 1; + gyz(p.I) = 0; + gzz(p.I) = 1; + + kxx(p.I) = 0; + kxy(p.I) = 0; + kxz(p.I) = 0; + kyy(p.I) = 0; + kyz(p.I) = 0; + kzz(p.I) = 0; + }); +} + +extern "C" void CurvADM_initial_lapse(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_CurvADM_initial_lapse; + DECLARE_CCTK_PARAMETERS; + + grid.loop_all_device<1, 1, 1>( + grid.nghostzones, [=] CCTK_DEVICE(const PointDesc &p) + CCTK_ATTRIBUTE_ALWAYS_INLINE { alp(p.I) = 1; }); +} + +extern "C" void CurvADM_initial_dtlapse(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_CurvADM_initial_dtlapse; + DECLARE_CCTK_PARAMETERS; + + grid.loop_all_device<1, 1, 1>( + grid.nghostzones, [=] CCTK_DEVICE(const PointDesc &p) + CCTK_ATTRIBUTE_ALWAYS_INLINE { dtalp(p.I) = 0; }); +} + +extern "C" void CurvADM_initial_shift(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_CurvADM_initial_shift; + DECLARE_CCTK_PARAMETERS; + + grid.loop_all_device<1, 1, 1>(grid.nghostzones, + [=] CCTK_DEVICE(const PointDesc &p) + CCTK_ATTRIBUTE_ALWAYS_INLINE { + betax(p.I) = 0; + betay(p.I) = 0; + betaz(p.I) = 0; + }); +} + +extern "C" void CurvADM_initial_dtshift(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_CurvADM_initial_dtshift; + DECLARE_CCTK_PARAMETERS; + + grid.loop_all_device<1, 1, 1>(grid.nghostzones, + [=] CCTK_DEVICE(const PointDesc &p) + CCTK_ATTRIBUTE_ALWAYS_INLINE { + dtbetax(p.I) = 0; + dtbetay(p.I) = 0; + dtbetaz(p.I) = 0; + }); +} + +} // namespace CurvADM diff --git a/CurvADM/src/linear_wave.cxx b/CurvADM/src/linear_wave.cxx new file mode 100644 index 0000000..20c1b5b --- /dev/null +++ b/CurvADM/src/linear_wave.cxx @@ -0,0 +1,47 @@ +#include + +#include +#include +#include + +#include + +namespace CurvADM { +using namespace Loop; +using namespace std; + +extern "C" void CurvADM_linear_wave(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_CurvADM_linear_wave; + DECLARE_CCTK_PARAMETERS; + + const CCTK_REAL t = cctk_time; + + // See arXiv:1111.2177 [gr-qc], (74-75) + + const auto b = [&](const PointDesc &p) { + return linear_wave_amplitude * + sin(2 * CCTK_REAL(M_PI) * (p.x - t) / linear_wave_wavelength); + }; + + const auto bt = [&](const PointDesc &p) { + return -2 * CCTK_REAL(M_PI) * linear_wave_amplitude / + linear_wave_wavelength * + cos(2 * CCTK_REAL(M_PI) * (p.x - t) / linear_wave_wavelength); + }; + + loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { gxx(p.I) = 1; }); + loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { gxy(p.I) = 0; }); + loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { gxz(p.I) = 0; }); + loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { gyy(p.I) = 1 + b(p); }); + loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { gyz(p.I) = 0; }); + loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { gzz(p.I) = 1 - b(p); }); + + loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { kxx(p.I) = 0; }); + loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { kxy(p.I) = 0; }); + loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { kxz(p.I) = 0; }); + loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { kyy(p.I) = bt(p) / 2; }); + loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { kyz(p.I) = 0; }); + loop_all<1, 1, 1>(cctkGH, [&](const PointDesc &p) { kzz(p.I) = -bt(p) / 2; }); +} + +} // namespace CurvADM diff --git a/CurvADM/src/make.code.defn b/CurvADM/src/make.code.defn new file mode 100644 index 0000000..e0d1b0a --- /dev/null +++ b/CurvADM/src/make.code.defn @@ -0,0 +1,7 @@ +# Main make.code.defn file for thorn CurvADM + +# Source files in this directory +SRCS = adm.cxx linear_wave.cxx + +# Subdirectories containing source files +SUBDIRS =