diff --git a/AsterSeeds/configuration.ccl b/AsterSeeds/configuration.ccl index 7cd85446..39bb8da4 100644 --- a/AsterSeeds/configuration.ccl +++ b/AsterSeeds/configuration.ccl @@ -1,3 +1,3 @@ # Configuration definitions for thorn AsterSeeds -REQUIRES Loop EOSX +REQUIRES Loop EOSX AsterUtils diff --git a/AsterSeeds/interface.ccl b/AsterSeeds/interface.ccl index c6f782e3..f86ac869 100644 --- a/AsterSeeds/interface.ccl +++ b/AsterSeeds/interface.ccl @@ -5,6 +5,7 @@ INHERITS: HydroBaseX AsterX USES INCLUDE HEADER: loop_device.hxx USES INCLUDE HEADER: eos.hxx eos_idealgas.hxx +USES INCLUDE HEADER: aster_fd.hxx aster_interp.hxx aster_utils.hxx CCTK_REAL Avec_cent TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"' { diff --git a/AsterSeeds/param.ccl b/AsterSeeds/param.ccl index 803b9f45..14e1219a 100644 --- a/AsterSeeds/param.ccl +++ b/AsterSeeds/param.ccl @@ -26,6 +26,7 @@ KEYWORD test_case "Name of the testcase" STEERABLE=never "spherical shock" :: "" "magTOV" :: "" + "magBNS" :: "" } "Balsara1" private: @@ -114,17 +115,19 @@ REAL r0 "Characteristic radial offset" STEERABLE=ALWAYS *:* :: "Anything" } 0.0 -REAL dipole_x "x-coordinate of dipole center" STEERABLE=ALWAYS +# coordinates of dipole center, to be based on the location of maximum of rho + +REAL dipole_x[2] "x-coordinate of dipole center" STEERABLE=ALWAYS { *:* :: "Anything" } 0.0 -REAL dipole_y "y-coordinate of the dipole center" STEERABLE=ALWAYS +REAL dipole_y[2] "y-coordinate of the dipole center" STEERABLE=ALWAYS { *:* :: "Anything" } 0.0 -REAL dipole_z "z-coordinate of the dipole center" STEERABLE=ALWAYS +REAL dipole_z[2] "z-coordinate of the dipole center" STEERABLE=ALWAYS { *:* :: "Anything" } 0.0 diff --git a/AsterSeeds/schedule.ccl b/AsterSeeds/schedule.ccl index d5c93547..95e377b5 100644 --- a/AsterSeeds/schedule.ccl +++ b/AsterSeeds/schedule.ccl @@ -53,14 +53,14 @@ if (CCTK_Equals(test_type, "3DTest")) { if (CCTK_Equals(test_case, "magTOV")) { STORAGE: Avec_cent - SCHEDULE AsterSeeds_InitializeCenteredAvec AT initial AFTER HydroBaseX_PostInitial BEFORE AsterX_InitialGroup { + SCHEDULE AsterSeeds_InitializeCenteredAvec_TOV AT initial AFTER HydroBaseX_PostInitial BEFORE AsterX_InitialGroup { LANG: C READS: HydroBaseX::press(everywhere) WRITES: Avec_x_cent(everywhere) Avec_y_cent(everywhere) Avec_z_cent(everywhere) } "Set up initial conditions for the cell-centered vector potential" - SCHEDULE AsterSeeds_InitializeStagAvec AT initial AFTER AsterSeeds_InitializeCenteredAvec BEFORE AsterX_InitialGroup { + SCHEDULE AsterSeeds_InitializeStagAvec_TOV AT initial AFTER AsterSeeds_InitializeCenteredAvec_TOV BEFORE AsterX_InitialGroup { LANG: C READS: Avec_x_cent(everywhere) Avec_y_cent(everywhere) Avec_z_cent(everywhere) WRITES: AsterX::Avec_x(interior) AsterX::Avec_y(interior) AsterX::Avec_z(interior) @@ -68,4 +68,26 @@ if (CCTK_Equals(test_type, "3DTest")) { } "Set up initial conditions for the vector potential" } + + #Initial conditions for magnetic field for BNS + + if (CCTK_Equals(test_case, "magBNS")) { + + STORAGE: Avec_cent + SCHEDULE AsterSeeds_InitializeCenteredAvec_BNS AT initial AFTER HydroBaseX_PostInitial BEFORE AsterX_InitialGroup { + LANG: C + READS: HydroBaseX::press(everywhere) + WRITES: Avec_x_cent(everywhere) Avec_y_cent(everywhere) Avec_z_cent(everywhere) + } + "Set up initial conditions for the cell-centered vector potential" + + SCHEDULE AsterSeeds_InitializeStagAvec_BNS AT initial AFTER AsterSeeds_InitializeCenteredAvec_BNS BEFORE AsterX_InitialGroup { + LANG: C + READS: Avec_x_cent(everywhere) Avec_y_cent(everywhere) Avec_z_cent(everywhere) + WRITES: AsterX::Avec_x(interior) AsterX::Avec_y(interior) AsterX::Avec_z(interior) + SYNC: AsterX::Avec_x AsterX::Avec_y AsterX::Avec_z + } + "Set up initial conditions for the vector potential" + } + } diff --git a/AsterSeeds/src/1D_tests.cxx b/AsterSeeds/src/1D_tests.cxx index 4fee21c0..46e97642 100644 --- a/AsterSeeds/src/1D_tests.cxx +++ b/AsterSeeds/src/1D_tests.cxx @@ -5,15 +5,16 @@ #include #include -#include -#include -#include +#include "eos.hxx" +#include "eos_idealgas.hxx" +#include "seeds_utils.hxx" namespace AsterSeeds { using namespace std; using namespace Loop; using namespace EOSX; +using namespace AsterUtils; extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTSX_Tests1D_Initialize; @@ -29,9 +30,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { if (CCTK_EQUALS(test_case, "equilibrium")) { - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { rho(p.I) = 1.0; velx(p.I) = 0.0; vely(p.I) = 0.0; @@ -41,25 +42,25 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = 0.0; }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = 0.0; }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = 0.0; }); } else if (CCTK_EQUALS(test_case, "sound wave")) { - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { rho(p.I) = 1.0; velx(p.I) = 0.0 + amplitude * sin(M_PI * p.x); vely(p.I) = 0.0; @@ -69,19 +70,19 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = 0.0; }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = 0.0; }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = 0.0; }); } else if (CCTK_EQUALS(test_case, "Alfven wave")) { @@ -89,9 +90,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { const CCTK_REAL va = 0.5; const CCTK_REAL k = 2*M_PI; - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { rho(p.I) = 1.0; velx(p.I) = 0.0; vely(p.I) = -va * A0 * cos(k*p.x); @@ -101,27 +102,27 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = p.z*cos(k*p.x) - p.y*sin(k*p.x); }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = -p.z/2.0; }); //CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = p.x*sin(k*p.x) - p.z; }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = p.y/2.0; }); //CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = p.y - p.x*cos(k*p.x); }); } else if (CCTK_EQUALS(test_case, "shock tube")) { - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.x <= 0.0) { rho(p.I) = 2.0; velx(p.I) = 0.0; @@ -139,19 +140,19 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = 0.0; }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = 0.0; }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = 0.0; }); } else if (CCTK_EQUALS(test_case, "Balsara1")) { @@ -176,9 +177,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { if (CCTK_EQUALS(shock_dir, "x")) { - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.x <= 0.0) { rho(p.I) = rhol; velx(p.I) = vxl; @@ -197,9 +198,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.x <= 0.0) { Avec_x(p.I) = Byl * (p.z) - Bzl * (p.y); @@ -207,14 +208,14 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Avec_x(p.I) = Byr * (p.z) - Bzr * (p.y); } }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = 0.0; }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = Bxr * (p.y); }); @@ -228,9 +229,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Bzr = -1.0; Bxr = 0.0; - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.y <= 0.0) { rho(p.I) = rhol; vely(p.I) = vxl; @@ -249,23 +250,23 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.y <= 0.0) { Avec_y(p.I) = Bzl * (p.x) - Bxl * (p.z); } else { Avec_y(p.I) = Bzr * (p.x) - Bxr * (p.z); } }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = 0.0; }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = Byr * (p.z); }); @@ -279,9 +280,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Bxr = -1.0; Byr = 0.0; - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.z <= 0.0) { rho(p.I) = rhol; velz(p.I) = vxl; @@ -300,9 +301,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.z <= 0.0) { Avec_z(p.I) = Bxl * (p.y) - Byl * (p.x); @@ -310,14 +311,14 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Avec_z(p.I) = Bxr * (p.y) - Byr * (p.x); } }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = 0.0; }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = Bzr * (p.x); }); @@ -345,9 +346,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { if (CCTK_EQUALS(shock_dir, "x")) { - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.x <= 0.0) { rho(p.I) = rhol; velx(p.I) = vxl; @@ -366,9 +367,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.x <= 0.0) { Avec_x(p.I) = Byl * (p.z) - Bzl * (p.y); @@ -376,14 +377,14 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Avec_x(p.I) = Byr * (p.z) - Bzr * (p.y); } }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = 0.0; }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = Bxr * (p.y); }); @@ -397,9 +398,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Bzr = 0.7; Bxr = 0.7; - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.y <= 0.0) { rho(p.I) = rhol; vely(p.I) = vxl; @@ -418,23 +419,23 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.y <= 0.0) { Avec_y(p.I) = Bzl * (p.x) - Bxl * (p.z); } else { Avec_y(p.I) = Bzr * (p.x) - Bxr * (p.z); } }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = 0.0; }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = Byr * (p.z); }); @@ -448,9 +449,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Bxr = 0.7; Byr = 0.7; - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.z <= 0.0) { rho(p.I) = rhol; velz(p.I) = vxl; @@ -469,9 +470,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.z <= 0.0) { Avec_z(p.I) = Bxl * (p.y) - Byl * (p.x); @@ -479,14 +480,14 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Avec_z(p.I) = Bxr * (p.y) - Byr * (p.x); } }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = 0.0; }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = Bzr * (p.x); }); @@ -514,9 +515,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { if (CCTK_EQUALS(shock_dir, "x")) { - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.x <= 0.0) { rho(p.I) = rhol; velx(p.I) = vxl; @@ -535,9 +536,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.x <= 0.0) { Avec_x(p.I) = Byl * (p.z) - Bzl * (p.y); @@ -545,14 +546,14 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Avec_x(p.I) = Byr * (p.z) - Bzr * (p.y); } }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = 0.0; }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = Bxr * (p.y); }); @@ -566,9 +567,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Bzr = 0.7; Bxr = 0.7; - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.y <= 0.0) { rho(p.I) = rhol; vely(p.I) = vxl; @@ -587,23 +588,23 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.y <= 0.0) { Avec_y(p.I) = Bzl * (p.x) - Bxl * (p.z); } else { Avec_y(p.I) = Bzr * (p.x) - Bxr * (p.z); } }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = 0.0; }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = Byr * (p.z); }); @@ -617,9 +618,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Bxr = 0.7; Byr = 0.7; - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.z <= 0.0) { rho(p.I) = rhol; velz(p.I) = vxl; @@ -638,9 +639,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.z <= 0.0) { Avec_z(p.I) = Bxl * (p.y) - Byl * (p.x); @@ -648,14 +649,14 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Avec_z(p.I) = Bxr * (p.y) - Byr * (p.x); } }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = 0.0; }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = Bzr * (p.x); }); @@ -683,9 +684,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { if (CCTK_EQUALS(shock_dir, "x")) { - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.x <= 0.0) { rho(p.I) = rhol; velx(p.I) = vxl; @@ -704,9 +705,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.x <= 0.0) { Avec_x(p.I) = Byl * (p.z) - Bzl * (p.y); @@ -714,14 +715,14 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Avec_x(p.I) = Byr * (p.z) - Bzr * (p.y); } }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = 0.0; }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = Bxr * (p.y); }); @@ -735,9 +736,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Bzr = -7.0; Bxr = -7.0; - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.y <= 0.0) { rho(p.I) = rhol; vely(p.I) = vxl; @@ -756,23 +757,23 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.y <= 0.0) { Avec_y(p.I) = Bzl * (p.x) - Bxl * (p.z); } else { Avec_y(p.I) = Bzr * (p.x) - Bxr * (p.z); } }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = 0.0; }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = Byr * (p.z); }); @@ -786,9 +787,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Bxr = -7.0; Byr = -7.0; - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.z <= 0.0) { rho(p.I) = rhol; velz(p.I) = vxl; @@ -807,9 +808,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.z <= 0.0) { Avec_z(p.I) = Bxl * (p.y) - Byl * (p.x); @@ -817,14 +818,14 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Avec_z(p.I) = Bxr * (p.y) - Byr * (p.x); } }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = 0.0; }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = Bzr * (p.x); }); @@ -852,9 +853,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { if (CCTK_EQUALS(shock_dir, "x")) { - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.x <= 0.0) { rho(p.I) = rhol; velx(p.I) = vxl; @@ -873,9 +874,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.x <= 0.0) { Avec_x(p.I) = Byl * (p.z) - Bzl * (p.y); @@ -883,14 +884,14 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Avec_x(p.I) = Byr * (p.z) - Bzr * (p.y); } }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = 0.0; }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = Bxr * (p.y); }); @@ -904,9 +905,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Bzr = -0.7; Bxr = 0.5; - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.y <= 0.0) { rho(p.I) = rhol; vely(p.I) = vxl; @@ -925,23 +926,23 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.y <= 0.0) { Avec_y(p.I) = Bzl * (p.x) - Bxl * (p.z); } else { Avec_y(p.I) = Bzr * (p.x) - Bxr * (p.z); } }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = 0.0; }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = Byr * (p.z); }); @@ -955,9 +956,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Bxr = -0.7; Byr = 0.5; - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.z <= 0.0) { rho(p.I) = rhol; velz(p.I) = vxl; @@ -976,9 +977,9 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { if (p.z <= 0.0) { Avec_z(p.I) = Bxl * (p.y) - Byl * (p.x); @@ -986,14 +987,14 @@ extern "C" void Tests1D_Initialize(CCTK_ARGUMENTS) { Avec_z(p.I) = Bxr * (p.y) - Byr * (p.x); } }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = 0.0; }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = Bzr * (p.x); }); diff --git a/AsterSeeds/src/2D_tests.cxx b/AsterSeeds/src/2D_tests.cxx index 01679ac7..d1c844ef 100644 --- a/AsterSeeds/src/2D_tests.cxx +++ b/AsterSeeds/src/2D_tests.cxx @@ -5,15 +5,16 @@ #include #include -#include -#include -#include +#include "eos.hxx" +#include "eos_idealgas.hxx" +#include "seeds_utils.hxx" namespace AsterSeeds { using namespace std; using namespace Loop; using namespace EOSX; +using namespace AsterUtils; extern "C" void Tests2D_Initialize(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTSX_Tests2D_Initialize; @@ -29,9 +30,9 @@ extern "C" void Tests2D_Initialize(CCTK_ARGUMENTS) { // See Cipolletta et al (2020) and Del Zanna, Bucciantini, Londrillo (2003) if (CCTK_EQUALS(test_case, "magnetic rotor")) { - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { const CCTK_REAL r_cyl = sqrt(pow2(p.x) + pow2(p.y)); if (r_cyl < 0.1) { @@ -54,19 +55,19 @@ extern "C" void Tests2D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = 0.; }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = 0.; }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = p.y; }); } @@ -93,9 +94,9 @@ extern "C" void Tests2D_Initialize(CCTK_ARGUMENTS) { CCTK_VERROR("Invalid loop advection type"); } - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { rho(p.I) = 1.; press(p.I) = 3.; eps(p.I) = eos_th.eps_from_valid_rho_press_ye(rho(p.I), press(p.I), @@ -105,19 +106,19 @@ extern "C" void Tests2D_Initialize(CCTK_ARGUMENTS) { velz(p.I) = axial_vel; }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = 0.; }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = 0.; }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { const CCTK_REAL r_cyl = sqrt(pow2(p.x) + pow2(p.y)); const CCTK_REAL R_loop = 0.3; const CCTK_REAL A_loop = 0.001; @@ -131,9 +132,9 @@ extern "C" void Tests2D_Initialize(CCTK_ARGUMENTS) { // See Cipolletta et al (2020) else if (CCTK_EQUALS(test_case, "cylindrical blast")) { - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { const CCTK_REAL r_cyl = sqrt(pow2(p.x) + pow2(p.y)); const CCTK_REAL r_in = 0.8; const CCTK_REAL r_out = 1.0; @@ -171,19 +172,19 @@ extern "C" void Tests2D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = 0.0; }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = 0.0; }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = 0.1 * p.y; }); } diff --git a/AsterSeeds/src/3D_tests.cxx b/AsterSeeds/src/3D_tests.cxx index a31cd5a2..089007ec 100644 --- a/AsterSeeds/src/3D_tests.cxx +++ b/AsterSeeds/src/3D_tests.cxx @@ -5,15 +5,16 @@ #include #include -#include -#include -#include +#include "eos.hxx" +#include "eos_idealgas.hxx" +#include "seeds_utils.hxx" namespace AsterSeeds { using namespace std; using namespace Loop; using namespace EOSX; +using namespace AsterUtils; extern "C" void Tests3D_Initialize(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTSX_Tests3D_Initialize; @@ -29,9 +30,9 @@ extern "C" void Tests3D_Initialize(CCTK_ARGUMENTS) { if (CCTK_EQUALS(test_case, "spherical shock")) { - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { const CCTK_REAL r2 = pow2(p.x) + pow2(p.y) + pow2(p.z); if (r2 <= pow2(shock_radius)) { rho(p.I) = 2.0; @@ -50,19 +51,19 @@ extern "C" void Tests3D_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = 0.0; }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = 0.0; }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = 0.0; }); } else { diff --git a/AsterSeeds/src/atmosphere.cxx b/AsterSeeds/src/atmosphere.cxx index 45b754e0..cbff3691 100644 --- a/AsterSeeds/src/atmosphere.cxx +++ b/AsterSeeds/src/atmosphere.cxx @@ -5,15 +5,16 @@ #include #include -#include -#include -#include +#include "eos.hxx" +#include "eos_idealgas.hxx" +#include "seeds_utils.hxx" namespace AsterSeeds { using namespace std; using namespace Loop; using namespace EOSX; +using namespace AsterUtils; extern "C" void Atmosphere_Initialize(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTSX_Atmosphere_Initialize; @@ -27,9 +28,9 @@ extern "C" void Atmosphere_Initialize(CCTK_ARGUMENTS) { const eos_idealgas eos_th(gl_gamma, particle_mass, rgeps, rgrho, rgye); const CCTK_REAL dummy_ye = 0.5; - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { rho(p.I) = rho_atmosphere; velx(p.I) = 0.; @@ -40,19 +41,19 @@ extern "C" void Atmosphere_Initialize(CCTK_ARGUMENTS) { dummy_ye); }); - grid.loop_all_device<1, 0, 0>( + grid.loop_all<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = 0.; }); - grid.loop_all_device<0, 1, 0>( + grid.loop_all<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = 0.; }); - grid.loop_all_device<0, 0, 1>( + grid.loop_all<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const PointDesc &p) + [=] CCTK_HOST(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = 0.; }); } diff --git a/AsterSeeds/src/magBNS.cxx b/AsterSeeds/src/magBNS.cxx new file mode 100644 index 00000000..74e8c4d1 --- /dev/null +++ b/AsterSeeds/src/magBNS.cxx @@ -0,0 +1,118 @@ +#include + +#include +#include +#include + +#include + +#include "seeds_utils.hxx" + +namespace AsterSeeds { +using namespace std; +using namespace Loop; +using namespace AsterUtils; + +extern "C" void AsterSeeds_InitializeCenteredAvec_BNS(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_AsterSeeds_InitializeCenteredAvec_BNS; + DECLARE_CCTK_PARAMETERS; + + if (CCTK_EQUALS(Afield_config, "internal dipole")) { + + /* computing cell centered vector potential components */ + grid.loop_all<1, 1, 1>( + grid.nghostzones, + [=] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + // dummy initialization + CCTK_REAL x_local = 0.0; + CCTK_REAL y_local = 0.0; + + // For star 1 at minus side + if (p.x < 0) { + x_local = p.x - dipole_x[0]; + y_local = p.y - dipole_y[0]; + } + // For star 2 at minus side + else { + x_local = p.x - dipole_x[1]; + y_local = p.y - dipole_y[1]; + } + + CCTK_REAL Pcut = press_max * press_cut; + CCTK_REAL Pdiff = std::max(press(p.I) - Pcut, 0.0); + CCTK_REAL Aphi_local = Ab * pow(Pdiff, Avec_kappa); + Avec_x_cent(p.I) = -y_local * Aphi_local; + Avec_y_cent(p.I) = x_local * Aphi_local; + Avec_z_cent(p.I) = 0.0; + }); + + } else if (CCTK_EQUALS(Afield_config, "external dipole")) { + + /* computing cell centered vector potential components */ + grid.loop_all<1, 1, 1>( + grid.nghostzones, + [=] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + // For star 1 at minus side + CCTK_REAL x_local_s1 = p.x - dipole_x[0]; + CCTK_REAL y_local_s1 = p.y - dipole_y[0]; + CCTK_REAL z_local_s1 = p.z - dipole_z[0]; + CCTK_REAL cylrad2_s1 = + x_local_s1 * x_local_s1 + y_local_s1 * y_local_s1; + CCTK_REAL rsph_s1 = + sqrt(x_local_s1 * x_local_s1 + y_local_s1 * y_local_s1 + + z_local_s1 * z_local_s1); + + CCTK_REAL Aphi_local_s1 = + B0 * (pow(r0, 3.0) / (pow(r0, 3.0) + pow(rsph_s1, 3.0))) / + sqrt(cylrad2_s1 + 1.0e-16); + + // For star 2 at minus side + CCTK_REAL x_local_s2 = p.x - dipole_x[1]; + CCTK_REAL y_local_s2 = p.y - dipole_y[1]; + CCTK_REAL z_local_s2 = p.z - dipole_z[1]; + CCTK_REAL cylrad2_s2 = + x_local_s2 * x_local_s2 + y_local_s2 * y_local_s2; + CCTK_REAL rsph_s2 = + sqrt(x_local_s2 * x_local_s2 + y_local_s2 * y_local_s2 + + z_local_s2 * z_local_s2); + + CCTK_REAL Aphi_local_s2 = + B0 * (pow(r0, 3.0) / (pow(r0, 3.0) + pow(rsph_s2, 3.0))) / + sqrt(cylrad2_s2 + 1.0e-16); + + Avec_x_cent(p.I) = + -(y_local_s1 * Aphi_local_s1 + y_local_s2 * Aphi_local_s2); + Avec_y_cent(p.I) = + x_local_s1 * Aphi_local_s1 + x_local_s2 * Aphi_local_s2; + Avec_z_cent(p.I) = 0.0; + }); + + } else { + CCTK_ERROR("Vector potential configuration not defined"); + } +} + +extern "C" void AsterSeeds_InitializeStagAvec_BNS(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_AsterSeeds_InitializeStagAvec_BNS; + DECLARE_CCTK_PARAMETERS; + + grid.loop_int<1, 0, 0>( + grid.nghostzones, + [=] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + Avec_x(p.I) = calc_avg_c2e(Avec_x_cent, p, 0); + }); + + grid.loop_int<0, 1, 0>( + grid.nghostzones, + [=] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + Avec_y(p.I) = calc_avg_c2e(Avec_y_cent, p, 1); + }); + + grid.loop_int<0, 0, 1>( + grid.nghostzones, + [=] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + Avec_z(p.I) = calc_avg_c2e(Avec_z_cent, p, 2); + }); +} + +} // namespace AsterSeeds diff --git a/AsterSeeds/src/magtov.cxx b/AsterSeeds/src/magTOV.cxx similarity index 50% rename from AsterSeeds/src/magtov.cxx rename to AsterSeeds/src/magTOV.cxx index a859537d..d67d7a1f 100644 --- a/AsterSeeds/src/magtov.cxx +++ b/AsterSeeds/src/magTOV.cxx @@ -5,44 +5,52 @@ #include #include + #include "seeds_utils.hxx" namespace AsterSeeds { using namespace std; using namespace Loop; +using namespace AsterUtils; -extern "C" void AsterSeeds_InitializeCenteredAvec(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTSX_AsterSeeds_InitializeCenteredAvec; +extern "C" void AsterSeeds_InitializeCenteredAvec_TOV(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_AsterSeeds_InitializeCenteredAvec_TOV; DECLARE_CCTK_PARAMETERS; if (CCTK_EQUALS(Afield_config, "internal dipole")) { /* computing cell centered vector potential components */ - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + CCTK_REAL x_local = p.x - dipole_x[0]; + CCTK_REAL y_local = p.y - dipole_y[0]; CCTK_REAL Pcut = press_max * press_cut; CCTK_REAL Pdiff = std::max(press(p.I) - Pcut, 0.0); CCTK_REAL Aphi_local = Ab * pow(Pdiff, Avec_kappa); - Avec_x_cent(p.I) = -p.y * Aphi_local; - Avec_y_cent(p.I) = p.x * Aphi_local; + Avec_x_cent(p.I) = -y_local * Aphi_local; + Avec_y_cent(p.I) = x_local * Aphi_local; Avec_z_cent(p.I) = 0.0; }); } else if (CCTK_EQUALS(Afield_config, "external dipole")) { /* computing cell centered vector potential components */ - grid.loop_all_device<1, 1, 1>( + grid.loop_all<1, 1, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { - CCTK_REAL cylrad2 = p.x * p.x + p.y * p.y; - CCTK_REAL rsph = sqrt(p.x * p.x + p.y * p.y + p.z * p.z); + [=] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + CCTK_REAL x_local = p.x - dipole_x[0]; + CCTK_REAL y_local = p.y - dipole_y[0]; + CCTK_REAL z_local = p.z - dipole_z[0]; + CCTK_REAL cylrad2 = x_local * x_local + y_local * y_local; + CCTK_REAL rsph = + sqrt(x_local * x_local + y_local * y_local + z_local * z_local); CCTK_REAL rsph3 = pow(rsph, 3.0); CCTK_REAL r03 = pow(r0, 3.0); CCTK_REAL Aphi_local = B0 * (r03 / (r03 + rsph3)) / sqrt(cylrad2 + 1.0e-16); - Avec_x_cent(p.I) = -p.y * Aphi_local; - Avec_y_cent(p.I) = p.x * Aphi_local; + Avec_x_cent(p.I) = -y_local * Aphi_local; + Avec_y_cent(p.I) = x_local * Aphi_local; Avec_z_cent(p.I) = 0.0; }); @@ -51,25 +59,25 @@ extern "C" void AsterSeeds_InitializeCenteredAvec(CCTK_ARGUMENTS) { } } -extern "C" void AsterSeeds_InitializeStagAvec(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTSX_AsterSeeds_InitializeStagAvec; +extern "C" void AsterSeeds_InitializeStagAvec_TOV(CCTK_ARGUMENTS) { + DECLARE_CCTK_ARGUMENTSX_AsterSeeds_InitializeStagAvec_TOV; DECLARE_CCTK_PARAMETERS; - grid.loop_int_device<1, 0, 0>( + grid.loop_int<1, 0, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_x(p.I) = calc_avg_c2e(Avec_x_cent, p, 0); }); - grid.loop_int_device<0, 1, 0>( + grid.loop_int<0, 1, 0>( grid.nghostzones, - [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_y(p.I) = calc_avg_c2e(Avec_y_cent, p, 1); }); - grid.loop_int_device<0, 0, 1>( + grid.loop_int<0, 0, 1>( grid.nghostzones, - [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { + [=] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE { Avec_z(p.I) = calc_avg_c2e(Avec_z_cent, p, 2); }); } diff --git a/AsterSeeds/src/make.code.defn b/AsterSeeds/src/make.code.defn index 91fb98a4..864c3f03 100644 --- a/AsterSeeds/src/make.code.defn +++ b/AsterSeeds/src/make.code.defn @@ -6,7 +6,8 @@ SRCS = \ 2D_tests.cxx \ 3D_tests.cxx \ atmosphere.cxx \ - magtov.cxx + magTOV.cxx \ + magBNS.cxx # Subdirectories containing source files diff --git a/AsterSeeds/src/seeds_utils.hxx b/AsterSeeds/src/seeds_utils.hxx index aefafcff..2aa4cc24 100644 --- a/AsterSeeds/src/seeds_utils.hxx +++ b/AsterSeeds/src/seeds_utils.hxx @@ -9,30 +9,13 @@ #include #include +#include "aster_utils.hxx" + namespace AsterSeeds { using namespace std; using namespace Loop; +using namespace AsterUtils; -template -inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST T pow2(T x) { - return x * x; -} - -// Second-order average of cell-centered grid functions to edge center -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_avg_c2e(const GF3D2 &gf, const PointDesc &p, const int dir) { - T gf_avg = 0.0; - - for (int dk = 0; dk < (dir == 2 ? 1 : 2); ++dk) { - for (int dj = 0; dj < (dir == 1 ? 1 : 2); ++dj) { - for (int di = 0; di < (dir == 0 ? 1 : 2); ++di) { - gf_avg += gf(p.I - p.DI[0] * di - p.DI[1] * dj - p.DI[2] * dk); - } - } - } - return gf_avg / 4.0; -} } // namespace AsterSeeds diff --git a/AsterUtils/README b/AsterUtils/README new file mode 100644 index 00000000..f83f8a85 --- /dev/null +++ b/AsterUtils/README @@ -0,0 +1,9 @@ +Cactus Code Thorn AsterUtils +Author(s) : Jay Kalinani , Liwei Ji +Maintainer(s): Jay Kalinani , Liwei Ji +Licence : LGPL +-------------------------------------------------------------------------- + +1. Purpose + +Provides utility functions for interpolations, finite differencing and more diff --git a/AsterUtils/configuration.ccl b/AsterUtils/configuration.ccl new file mode 100644 index 00000000..9f78b563 --- /dev/null +++ b/AsterUtils/configuration.ccl @@ -0,0 +1,6 @@ +# Configuration definitions for thorn AsterUtils + +REQUIRES Loop +PROVIDES AsterUtils { +} +REQUIRES AsterUtils diff --git a/AsterUtils/interface.ccl b/AsterUtils/interface.ccl new file mode 100644 index 00000000..e602630f --- /dev/null +++ b/AsterUtils/interface.ccl @@ -0,0 +1,8 @@ +# Interface definition for thorn AsterUtils + +IMPLEMENTS: AsterUtils + +USES INCLUDE HEADER: loop_device.hxx +INCLUDES HEADER: aster_fd.hxx IN aster_fd.hxx +INCLUDES HEADER: aster_interp.hxx IN aster_interp.hxx +INCLUDES HEADER: aster_utils.hxx IN aster_utils.hxx diff --git a/AsterUtils/param.ccl b/AsterUtils/param.ccl new file mode 100644 index 00000000..5a2d0eed --- /dev/null +++ b/AsterUtils/param.ccl @@ -0,0 +1,3 @@ +# Parameter definitions for thorn AsterUtils + + diff --git a/AsterUtils/schedule.ccl b/AsterUtils/schedule.ccl new file mode 100644 index 00000000..f61e2887 --- /dev/null +++ b/AsterUtils/schedule.ccl @@ -0,0 +1,2 @@ +#Schedule definitions for thorn AsterUtils + diff --git a/AsterX/src/fd.hxx b/AsterUtils/src/aster_fd.hxx similarity index 97% rename from AsterX/src/fd.hxx rename to AsterUtils/src/aster_fd.hxx index b83164f8..cef1031a 100644 --- a/AsterX/src/fd.hxx +++ b/AsterUtils/src/aster_fd.hxx @@ -1,5 +1,5 @@ -#ifndef ASTERX_FD_HXX -#define ASTERX_FD_HXX +#ifndef ASTER_FD_HXX +#define ASTER_FD_HXX #include #include @@ -15,7 +15,7 @@ #include #include -namespace AsterX { +namespace AsterUtils { using namespace std; using namespace Loop; using namespace Arith; @@ -109,6 +109,6 @@ CCTK_DEVICE CCTK_HOST return 0.25 * (dgf1 + dgf2 + dgf3 + dgf4); } -} // namespace AsterX +} // namespace AsterUtils -#endif // ASTERX_FD_HXX +#endif // ASTER_FD_HXX diff --git a/AsterX/src/interp.hxx b/AsterUtils/src/aster_interp.hxx similarity index 97% rename from AsterX/src/interp.hxx rename to AsterUtils/src/aster_interp.hxx index d0860d97..801266ae 100644 --- a/AsterX/src/interp.hxx +++ b/AsterUtils/src/aster_interp.hxx @@ -1,5 +1,5 @@ -#ifndef ASTERX_INTERP_HXX -#define ASTERX_INTERP_HXX +#ifndef ASTER_INTERP_HXX +#define ASTER_INTERP_HXX #include #include @@ -14,7 +14,7 @@ #include #include -namespace AsterX { +namespace AsterUtils { using namespace std; using namespace Loop; using namespace Arith; @@ -147,6 +147,6 @@ CCTK_DEVICE CCTK_HOST return gf_avg / 3.0; } -} // namespace AsterX +} // namespace AsterUtils -#endif // ASTERX_INTERP_HXX +#endif // ASTER_INTERP_HXX diff --git a/AsterX/src/utils.hxx b/AsterUtils/src/aster_utils.hxx similarity index 95% rename from AsterX/src/utils.hxx rename to AsterUtils/src/aster_utils.hxx index 598415a4..20ca4c0e 100644 --- a/AsterX/src/utils.hxx +++ b/AsterUtils/src/aster_utils.hxx @@ -1,5 +1,5 @@ -#ifndef ASTERX_UTILS_HXX -#define ASTERX_UTILS_HXX +#ifndef ASTER_UTILS_HXX +#define ASTER_UTILS_HXX #include #include @@ -14,10 +14,10 @@ #include #include -#include -#include +#include "aster_fd.hxx" +#include "aster_interp.hxx" -namespace AsterX { +namespace AsterUtils { using namespace std; using namespace Loop; using namespace Arith; @@ -141,11 +141,6 @@ calc_wlorentz_zvec(const vec &zvec_up, const vec &zvec_dn) { return sqrt(1.0 + calc_contraction(zvec_up, zvec_dn)); } -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T pow2(T x) { - return x * x; -} - template CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec pow2(vec x) { @@ -169,6 +164,6 @@ calc_avg_neighbors(const vec flag, const vec u_nbs, CCTK_REAL(D); } -} // namespace AsterX +} // namespace AsterUtils -#endif // ASTERX_UTILS_HXX +#endif // ASTER_UTILS_HXX diff --git a/AsterUtils/src/make.code.defn b/AsterUtils/src/make.code.defn new file mode 100644 index 00000000..3d6ed9a0 --- /dev/null +++ b/AsterUtils/src/make.code.defn @@ -0,0 +1,8 @@ +# Main make.code.defn file for thorn AsterUtils +# $Header:$ + +# Source files in this directory +SRCS = + +# Subdirectories containing source files +SUBDIRS = diff --git a/AsterX/configuration.ccl b/AsterX/configuration.ccl index 81cefe27..8320a5e3 100644 --- a/AsterX/configuration.ccl +++ b/AsterX/configuration.ccl @@ -1,3 +1,3 @@ # Configuration definitions for thorn AsterX -REQUIRES Loop EOSX Con2PrimFactory ReconX +REQUIRES Loop EOSX Con2PrimFactory ReconX AsterUtils diff --git a/AsterX/interface.ccl b/AsterX/interface.ccl index 3d797c5b..32b69e4e 100644 --- a/AsterX/interface.ccl +++ b/AsterX/interface.ccl @@ -11,6 +11,7 @@ USES INCLUDE HEADER: eos.hxx eos_idealgas.hxx USES INCLUDE HEADER: c2p.hxx c2p_2DNoble.hxx c2p_1DPalenzuela.hxx USES INCLUDE HEADER: reconstruct.hxx USES INCLUDE HEADER: reconx_utils.hxx +USES INCLUDE HEADER: aster_fd.hxx aster_interp.hxx aster_utils.hxx PUBLIC: diff --git a/AsterX/param.ccl b/AsterX/param.ccl index 42f3264e..b3cad544 100644 --- a/AsterX/param.ccl +++ b/AsterX/param.ccl @@ -116,6 +116,10 @@ USES CCTK_REAL r_atmo USES CCTK_REAL n_rho_atmo USES CCTK_REAL n_press_atmo USES BOOLEAN thermal_eos_atmo +USES CCTK_REAL alp_thresh +USES CCTK_REAL rho_BH +USES CCTK_REAL eps_BH +USES CCTK_REAL vwlim_BH SHARES: ReconX USES KEYWORD reconstruction_method diff --git a/AsterX/schedule.ccl b/AsterX/schedule.ccl index 0b872473..85e319a1 100644 --- a/AsterX/schedule.ccl +++ b/AsterX/schedule.ccl @@ -123,6 +123,7 @@ SCHEDULE AsterX_Con2Prim IN AsterX_Con2PrimGroup AFTER AsterX_ComputedBFromdBsta { LANG: C READS: ADMBaseX::metric(interior) + READS: ADMBaseX::lapse(everywhere) READS: dens(interior) tau(interior) mom(interior) dB(interior) READS: saved_prims(interior) READS: Avec_x(interior) Avec_y(interior) Avec_z(interior) diff --git a/AsterX/src/computeBfromA.cxx b/AsterX/src/computeBfromA.cxx index adf30842..53938d23 100644 --- a/AsterX/src/computeBfromA.cxx +++ b/AsterX/src/computeBfromA.cxx @@ -4,9 +4,10 @@ #include #include -#include "utils.hxx" #include +#include "aster_utils.hxx" + struct metric { CCTK_REAL gxx, gxy, gxz, gyy, gyz, gzz; }; @@ -14,6 +15,7 @@ struct metric { namespace AsterX { using namespace Loop; using namespace Arith; +using namespace AsterUtils; template void ComputeStaggeredB(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTSX_AsterX_ComputedBstagFromA; diff --git a/AsterX/src/con2prim.cxx b/AsterX/src/con2prim.cxx index 2bcb016d..62d6ddfc 100644 --- a/AsterX/src/con2prim.cxx +++ b/AsterX/src/con2prim.cxx @@ -1,28 +1,27 @@ +#include + #include #include #include - -#include - #include #include "c2p.hxx" #include "c2p_1DPalenzuela.hxx" #include "c2p_2DNoble.hxx" -#include -#include - -#include -#include +#include "eos_1p.hxx" +#include "eos_polytropic.hxx" +#include "eos.hxx" +#include "eos_idealgas.hxx" -#include "utils.hxx" +#include "aster_utils.hxx" namespace AsterX { using namespace std; using namespace Loop; using namespace EOSX; using namespace Con2PrimFactory; +using namespace AsterUtils; enum class eos_t { IdealGas, Hybrid, Tabulated }; enum class c2p_first_t { Noble, Palenzuela }; @@ -130,6 +129,32 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, atmo.set(pv_seeds); } + // Modifying primitive seeds within BH interiors before C2Ps are called + // NOTE: By default, alp_thresh=0 so the if condition below is never + // triggered. One must be very careful when using this functionality and + // must correctly set alp_thresh, rho_BH, eps_BH and vwlim_BH in the parfile + + if (alp(p.I) < alp_thresh) { + if ((pv_seeds.rho > rho_BH) || (pv_seeds.eps > eps_BH)) { + pv_seeds.rho = rho_BH; // typically set to 0.01% to 1% of rho_max of + // initial NS or disk + pv_seeds.eps = eps_BH; + pv_seeds.Ye = Ye_atmo; + pv_seeds.press = + eos_th.press_from_valid_rho_eps_ye(rho_BH, eps_BH, Ye_atmo); + // check on velocities + CCTK_REAL wlim_BH = sqrt(1.0 + vwlim_BH * vwlim_BH); + CCTK_REAL vlim_BH = vwlim_BH / wlim_BH; + CCTK_REAL sol_v = + sqrt((pv_seeds.w_lor * pv_seeds.w_lor - 1.0)) / pv_seeds.w_lor; + if (sol_v > vlim_BH) { + pv_seeds.vel *= vlim_BH / sol_v; + pv_seeds.w_lor = wlim_BH; + } + cv.from_prim(pv_seeds, glo); + } + } + // Construct error report object: c2p_report rep_first; c2p_report rep_second; @@ -170,6 +195,33 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, if (rep_first.failed() && rep_second.failed()) { printf("Second C2P failed too :( :( \n"); rep_second.debug_message(); + + // Treatment for BH interiors after C2P failures + // NOTE: By default, alp_thresh=0 so the if condition below is never + // triggered. One must be very careful when using this functionality and + // must correctly set alp_thresh, rho_BH, eps_BH and vwlim_BH in the + // parfile + if (alp(p.I) < alp_thresh) { + if ((pv_seeds.rho > rho_BH) || (pv_seeds.eps > eps_BH)) { + pv.rho = rho_BH; // typically set to 0.01% to 1% of rho_max of initial + // NS or disk + pv.eps = eps_BH; + pv.Ye = Ye_atmo; + pv.press = + eos_th.press_from_valid_rho_eps_ye(rho_BH, eps_BH, Ye_atmo); + // check on velocities + CCTK_REAL wlim_BH = sqrt(1.0 + vwlim_BH * vwlim_BH); + CCTK_REAL vlim_BH = vwlim_BH / wlim_BH; + CCTK_REAL sol_v = sqrt((pv.w_lor * pv.w_lor - 1.0)) / pv.w_lor; + if (sol_v > vlim_BH) { + pv.vel *= vlim_BH / sol_v; + pv.w_lor = wlim_BH; + } + cv.from_prim(pv, glo); + rep_first.set_atmo = 0; + rep_second.set_atmo = 0; + } + } con2prim_flag(p.I) = 0; } @@ -219,13 +271,16 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold, pv.scatter(rho(p.I), eps(p.I), dummy_Ye, press(p.I), velx(p.I), vely(p.I), velz(p.I), wlor, Bvecx(p.I), Bvecy(p.I), Bvecz(p.I), Ex, Ey, Ez); - zvec_x(p.I) = wlor * pv.vel(0); + zvec_x(p.I) = wlor * pv.vel(0); zvec_y(p.I) = wlor * pv.vel(1); zvec_z(p.I) = wlor * pv.vel(2); - svec_x(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(0); - svec_y(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(1); - svec_z(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(2); + svec_x(p.I) = + (pv.rho + pv.rho * pv.eps + pv.press) * wlor * wlor * pv.vel(0); + svec_y(p.I) = + (pv.rho + pv.rho * pv.eps + pv.press) * wlor * wlor * pv.vel(1); + svec_z(p.I) = + (pv.rho + pv.rho * pv.eps + pv.press) * wlor * wlor * pv.vel(2); // Write back cv cv.scatter(dens(p.I), momx(p.I), momy(p.I), momz(p.I), tau(p.I), dummy_Ye, diff --git a/AsterX/src/con2prim_rpa.cxx b/AsterX/src/con2prim_rpa.cxx index 660822b1..5fc50b0c 100644 --- a/AsterX/src/con2prim_rpa.cxx +++ b/AsterX/src/con2prim_rpa.cxx @@ -1,12 +1,12 @@ +#include + #include #include #include -#include - #include -#include "utils.hxx" +#include "aster_utils.hxx" #include "reprimand/eos_thermal.h" // The EOS framework #include "reprimand/eos_idealgas.h" @@ -18,6 +18,7 @@ namespace AsterX { using namespace std; using namespace Loop; using namespace EOS_Toolkit; +using namespace AsterUtils; extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTSX_AsterX_Con2Prim; @@ -58,10 +59,10 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) { press_atm = (radial_distance > r_atmo) ? (p_atmo * pow(r_atmo / radial_distance, n_press_atmo)) : p_atmo; - //TODO: eos.at_rho_press_ye(rho_atm, press_atm, Ye_atmo).eps() does not exist in RePrimAnd - //Currently computing eps from ideal gas EOS - //eps_atm = eos.at_rho_press_ye(rho_atm, press_atm, Ye_atmo).eps(); - eps_atm = press_atm/(rho_atm*(gl_gamma - 1.)); + // TODO: eos.at_rho_press_ye(rho_atm, press_atm, Ye_atmo).eps() does not + // exist in RePrimAnd Currently computing eps from ideal gas EOS eps_atm = + // eos.at_rho_press_ye(rho_atm, press_atm, Ye_atmo).eps(); + eps_atm = press_atm / (rho_atm * (gl_gamma - 1.)); } else { eps_atm = eos_id.at_rho(rho_atm).eps(); eps_atm = eos.range_eps(rho_atm, Ye_atmo).limit_to(eps_atm); @@ -73,8 +74,8 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) { CCTK_REAL dummy_dYe = 0.5; // Get a recovery function - con2prim_mhd cv2pv(eos, rho_strict, Ye_lenient, vw_lim, B_lim, atmo, c2p_tol, - max_iter); + con2prim_mhd cv2pv(eos, rho_strict, Ye_lenient, vw_lim, B_lim, atmo, + c2p_tol, max_iter); /* Get covariant metric */ const smat glo( @@ -83,6 +84,10 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) { sm_metric3 g(sm_symt3l(glo(0, 0), glo(0, 1), glo(1, 1), glo(0, 2), glo(1, 2), glo(2, 2))); + /* Get covariant metric used for BH interior fixes */ + const smat glow( + [&](int i, int j) ARITH_INLINE { return calc_avg_v2c(gf_g(i, j), p); }); + /* Calculate inverse of 3-metric */ const CCTK_REAL spatial_detg = calc_det(glo); const CCTK_REAL sqrt_detg = sqrt(spatial_detg); @@ -104,6 +109,74 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) { {momx(p.I), momy(p.I), momz(p.I)}, {dBx(p.I), dBy(p.I), dBz(p.I)}}; + // Modifying primitive seeds within BH interiors before C2Ps are called + // NOTE: By default, alp_thresh=0 so the if condition below is never + // triggered. One must be very careful when using this functionality and + // must correctly set alp_thresh, rho_BH, eps_BH and vwlim_BH in the parfile + + if (alp(p.I) < alp_thresh) { + if ((pv_seeds.rho > rho_BH) || (pv_seeds.eps > eps_BH)) { + pv_seeds.rho = rho_BH; // typically set to 0.01% to 1% of rho_max of + // initial NS or disk + pv_seeds.eps = eps_BH; + pv_seeds.ye = Ye_atmo; + pv_seeds.press = eos.at_rho_eps_ye(rho_BH, eps_BH, Ye_atmo).press(); + // check on velocities + CCTK_REAL wlim_BH = sqrt(1.0 + vwlim_BH * vwlim_BH); + CCTK_REAL vlim_BH = vwlim_BH / wlim_BH; + CCTK_REAL sol_v = + sqrt((pv_seeds.w_lor * pv_seeds.w_lor - 1.0)) / pv_seeds.w_lor; + if (sol_v > vlim_BH) { + pv_seeds.vel *= vlim_BH / sol_v; + pv_seeds.w_lor = wlim_BH; + } + + // cv.from_prim(pv_seeds, g); + // We do not save electric field, required by cv.from_prim(pv_seeds, g) + // in RPA. Thus, we recompute the CVs explicitly below + const vec &v_up = pv_seeds.vel; + const vec v_low = calc_contraction(glow, v_up); + /* Computing B_j */ + const vec &B_up = pv_seeds.B; + const vec B_low = calc_contraction(glow, B_up); + /* Computing b^t : this is b^0 * alp */ + const CCTK_REAL bst = pv_seeds.w_lor * calc_contraction(B_up, v_low); + /* Computing b_j */ + const vec b_low = B_low / pv_seeds.w_lor + bst * v_low; + /* Computing b^mu b_mu */ + const CCTK_REAL bs2 = (calc_contraction(B_up, B_low) + bst * bst) / + (pv_seeds.w_lor * pv_seeds.w_lor); + // computing conserved from primitives + cv.dens = sqrt_detg * pv_seeds.rho * pv_seeds.w_lor; + cv.scon(0) = + sqrt_detg * + (pv_seeds.w_lor * pv_seeds.w_lor * + (pv_seeds.rho * (1.0 + pv_seeds.eps) + pv_seeds.press + bs2) * + v_low(0) - + bst * b_low(0)); + cv.scon(1) = + sqrt_detg * + (pv_seeds.w_lor * pv_seeds.w_lor * + (pv_seeds.rho * (1.0 + pv_seeds.eps) + pv_seeds.press + bs2) * + v_low(1) - + bst * b_low(1)); + cv.scon(2) = + sqrt_detg * + (pv_seeds.w_lor * pv_seeds.w_lor * + (pv_seeds.rho * (1.0 + pv_seeds.eps) + pv_seeds.press + bs2) * + v_low(2) - + bst * b_low(2)); + cv.tau = sqrt_detg * (pv_seeds.w_lor * pv_seeds.w_lor * + (pv_seeds.rho * (1.0 + pv_seeds.eps) + + pv_seeds.press + bs2) - + (pv_seeds.press + 0.5 * bs2) - bst * bst) - + cv.dens; + cv.bcons = sqrt_detg * pv_seeds.B; + cv.tracer_ye = cv.dens * pv_seeds.ye; + } + } + + // Calling C2P cv2pv(pv, cv, g, rep); // Handle incorrectable errors @@ -114,8 +187,8 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) { // need to fix pv to computed values like pv.rho instead of rho(p.I) printf( "WARNING: " - "C2P failed. Printing cons and saved prims before set to " - "atmo: \n" + "C2P failed. Printing cons and saved prims before set to atmo or " + "BH interior fix: \n" "cctk_iteration = %i \n " "x, y, z = %26.16e, %26.16e, %26.16e \n " "gxx, gxy, gxz, gyy, gyz, gzz = %f, %f, %f, %f, %f, %f \n " @@ -137,12 +210,70 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) { // velz(p.I), Bvecx(p.I), Bvecy(p.I), Bvecz(p.I), Avec_x(p.I), Avec_y(p.I), Avec_z(p.I)); } - // set to atmo - cv.bcons(0) = dBx(p.I); - cv.bcons(1) = dBy(p.I); - cv.bcons(2) = dBz(p.I); - pv.B = cv.bcons / sqrt_detg; - atmo.set(pv, cv, g); + + if (alp(p.I) < alp_thresh) { + if ((pv.rho > rho_BH) || (pv.eps > eps_BH)) { + pv.rho = rho_BH; // typically set to 0.01% to 1% of rho_max of + // initial NS or disk + pv.eps = eps_BH; + pv.ye = Ye_atmo; + pv.press = eos.at_rho_eps_ye(rho_BH, eps_BH, Ye_atmo).press(); + // check on velocities + CCTK_REAL wlim_BH = sqrt(1.0 + vwlim_BH * vwlim_BH); + CCTK_REAL vlim_BH = vwlim_BH / wlim_BH; + CCTK_REAL sol_v = sqrt((pv.w_lor * pv.w_lor - 1.0)) / pv.w_lor; + if (sol_v > vlim_BH) { + pv.vel *= vlim_BH / sol_v; + pv.w_lor = wlim_BH; + } + // cv.from_prim(pv, g); + // We do not save electric field, required by cv.from_prim(pv, g) in + // RPA. Thus, we recompute the CVs explicitly below + + const vec &v_up = pv.vel; + const vec v_low = calc_contraction(glow, v_up); + /* Computing B_j */ + const vec &B_up = pv.B; + const vec B_low = calc_contraction(glow, B_up); + /* Computing b^t : this is b^0 * alp */ + const CCTK_REAL bst = pv.w_lor * calc_contraction(B_up, v_low); + /* Computing b_j */ + const vec b_low = B_low / pv.w_lor + bst * v_low; + /* Computing b^mu b_mu */ + const CCTK_REAL bs2 = (calc_contraction(B_up, B_low) + bst * bst) / + (pv.w_lor * pv.w_lor); + // computing conserved from primitives + cv.dens = sqrt_detg * pv.rho * pv.w_lor; + cv.scon(0) = + sqrt_detg * + (pv.w_lor * pv.w_lor * + (pv.rho * (1.0 + pv.eps) + pv.press + bs2) * v_low(0) - + bst * b_low(0)); + cv.scon(1) = + sqrt_detg * + (pv.w_lor * pv.w_lor * + (pv.rho * (1.0 + pv.eps) + pv.press + bs2) * v_low(1) - + bst * b_low(1)); + cv.scon(2) = + sqrt_detg * + (pv.w_lor * pv.w_lor * + (pv.rho * (1.0 + pv.eps) + pv.press + bs2) * v_low(2) - + bst * b_low(2)); + cv.tau = sqrt_detg * (pv.w_lor * pv.w_lor * + (pv.rho * (1.0 + pv.eps) + pv.press + bs2) - + (pv.press + 0.5 * bs2) - bst * bst) - + cv.dens; + cv.bcons = sqrt_detg * pv.B; + cv.tracer_ye = cv.dens * pv.ye; + } + } else { + // set to atmo + cv.bcons(0) = dBx(p.I); + cv.bcons(1) = dBy(p.I); + cv.bcons(2) = dBz(p.I); + pv.B = cv.bcons / sqrt_detg; + atmo.set(pv, cv, g); + } } /* set flag to success */ @@ -154,13 +285,16 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) { pv.scatter(rho(p.I), eps(p.I), dumye, press(p.I), velx(p.I), vely(p.I), velz(p.I), wlor, Ex, Ey, Ez, Bvecx(p.I), Bvecy(p.I), Bvecz(p.I)); - zvec_x(p.I) = wlor * pv.vel(0); + zvec_x(p.I) = wlor * pv.vel(0); zvec_y(p.I) = wlor * pv.vel(1); zvec_z(p.I) = wlor * pv.vel(2); - svec_x(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(0); - svec_y(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(1); - svec_z(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(2); + svec_x(p.I) = + (pv.rho + pv.rho * pv.eps + pv.press) * wlor * wlor * pv.vel(0); + svec_y(p.I) = + (pv.rho + pv.rho * pv.eps + pv.press) * wlor * wlor * pv.vel(1); + svec_z(p.I) = + (pv.rho + pv.rho * pv.eps + pv.press) * wlor * wlor * pv.vel(2); // Write back cv if (rep.adjust_cons) { diff --git a/AsterX/src/fluxes.cxx b/AsterX/src/fluxes.cxx index 893bf2bf..1fe0ba84 100644 --- a/AsterX/src/fluxes.cxx +++ b/AsterX/src/fluxes.cxx @@ -4,17 +4,23 @@ #include #include +#include +#include +#include +#include + #include #include #include #include -#include "utils.hxx" +#include "eos.hxx" +#include "eos_idealgas.hxx" + #include "eigenvalues.hxx" #include "fluxes.hxx" -#include -#include -#include +#include "reconstruct.hxx" +#include "aster_utils.hxx" namespace AsterX { using namespace std; @@ -22,6 +28,7 @@ using namespace Loop; using namespace Arith; using namespace EOSX; using namespace ReconX; +using namespace AsterUtils; enum class flux_t { LxF, HLLE }; enum class eos_t { IdealGas, Hybrid, Tabulated }; diff --git a/AsterX/src/prim2con.cxx b/AsterX/src/prim2con.cxx index 60c76463..0f435d76 100644 --- a/AsterX/src/prim2con.cxx +++ b/AsterX/src/prim2con.cxx @@ -1,10 +1,9 @@ -#include "prim2con.hxx" - -#include - #include #include #include +#include + +#include "prim2con.hxx" namespace AsterX { using namespace Loop; diff --git a/AsterX/src/prim2con.hxx b/AsterX/src/prim2con.hxx index a3d038b6..49c9b876 100644 --- a/AsterX/src/prim2con.hxx +++ b/AsterX/src/prim2con.hxx @@ -10,12 +10,14 @@ #include #include #include -#include "utils.hxx" + +#include "aster_utils.hxx" namespace AsterX { using namespace std; using namespace Loop; using namespace Arith; +using namespace AsterUtils; struct prim { CCTK_REAL rho; diff --git a/AsterX/src/rhs.cxx b/AsterX/src/rhs.cxx index 948603c2..ba161b32 100644 --- a/AsterX/src/rhs.cxx +++ b/AsterX/src/rhs.cxx @@ -1,22 +1,17 @@ -#include - #include #include #include - -#include -#include "utils.hxx" - -// #ifdef AMREX_USE_GPU -// #include -// #endif - +#include #include +#include "reconstruct.hxx" +#include "aster_utils.hxx" + namespace AsterX { using namespace Loop; using namespace Arith; using namespace ReconX; +using namespace AsterUtils; enum class vector_potential_gauge_t { algebraic, generalized_lorentz }; diff --git a/AsterX/src/source.cxx b/AsterX/src/source.cxx index 7e891983..d27487c5 100644 --- a/AsterX/src/source.cxx +++ b/AsterX/src/source.cxx @@ -8,12 +8,14 @@ #include #include #include -#include "utils.hxx" + +#include "aster_utils.hxx" namespace AsterX { using namespace std; using namespace Loop; using namespace Arith; +using namespace AsterUtils; template void SourceTerms(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTSX_AsterX_SourceTerms; diff --git a/AsterX/src/tmunu.cxx b/AsterX/src/tmunu.cxx index b0e96f5c..1ac54376 100644 --- a/AsterX/src/tmunu.cxx +++ b/AsterX/src/tmunu.cxx @@ -1,19 +1,20 @@ -#include - #include #include #include -#include "utils.hxx" #include #include #include #include +#include + +#include "aster_utils.hxx" namespace AsterX { using namespace std; using namespace Loop; using namespace Arith; +using namespace AsterUtils; template void Tmunu(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTSX_AsterX_Tmunu; diff --git a/AsterX/src/unit_tests/contraction_smat_upvec.cxx b/AsterX/src/unit_tests/contraction_smat_upvec.cxx index cb77402d..da9d93a8 100644 --- a/AsterX/src/unit_tests/contraction_smat_upvec.cxx +++ b/AsterX/src/unit_tests/contraction_smat_upvec.cxx @@ -4,14 +4,14 @@ #include #include "../test.hxx" -#include "../utils.hxx" +#include "aster_utils.hxx" #include void AsterXTests::test_contraction_smat_upvec(std::mt19937_64 &engine, int repetitions) { using namespace Arith; - using namespace AsterX; + using namespace AsterUtils; using std::uniform_real_distribution; uniform_real_distribution distrib{-1.0, 1.0}; @@ -61,4 +61,4 @@ void AsterXTests::test_contraction_smat_upvec(std::mt19937_64 &engine, } } } -} \ No newline at end of file +} diff --git a/AsterX/src/unit_tests/contraction_upvec_downvec.cxx b/AsterX/src/unit_tests/contraction_upvec_downvec.cxx index cd56f02b..6d31e291 100644 --- a/AsterX/src/unit_tests/contraction_upvec_downvec.cxx +++ b/AsterX/src/unit_tests/contraction_upvec_downvec.cxx @@ -4,14 +4,14 @@ #include #include "../test.hxx" -#include "../utils.hxx" +#include "aster_utils.hxx" #include void AsterXTests::test_contraction_upvec_downvec(std::mt19937_64 &engine, int repetitions) { using namespace Arith; - using namespace AsterX; + using namespace AsterUtils; using std::uniform_real_distribution; uniform_real_distribution distrib{-1.0, 1.0}; @@ -42,4 +42,4 @@ void AsterXTests::test_contraction_upvec_downvec(std::mt19937_64 &engine, } } } -} \ No newline at end of file +} diff --git a/AsterX/src/unit_tests/cross_product.cxx b/AsterX/src/unit_tests/cross_product.cxx index 07b232ec..10bb49ab 100644 --- a/AsterX/src/unit_tests/cross_product.cxx +++ b/AsterX/src/unit_tests/cross_product.cxx @@ -4,13 +4,13 @@ #include #include "../test.hxx" -#include "../utils.hxx" +#include "aster_utils.hxx" #include void AsterXTests::test_cross_product(std::mt19937_64 &engine, int repetitions) { using namespace Arith; - using namespace AsterX; + using namespace AsterUtils; using std::uniform_real_distribution; uniform_real_distribution distrib{-1.0, 1.0}; @@ -56,4 +56,4 @@ void AsterXTests::test_cross_product(std::mt19937_64 &engine, int repetitions) { } } } -} \ No newline at end of file +} diff --git a/AsterX/src/unit_tests/wlorentz.cxx b/AsterX/src/unit_tests/wlorentz.cxx index e3e2539d..1d8d4c3f 100644 --- a/AsterX/src/unit_tests/wlorentz.cxx +++ b/AsterX/src/unit_tests/wlorentz.cxx @@ -4,12 +4,12 @@ #include #include "../test.hxx" -#include "../utils.hxx" +#include "aster_utils.hxx" // Value must be bounded void AsterXTests::test_wlorentz(std::mt19937_64 &engine, int repetitions) { using namespace Arith; - using namespace AsterX; + using namespace AsterUtils; using std::sqrt; using std::uniform_real_distribution; @@ -43,4 +43,4 @@ void AsterXTests::test_wlorentz(std::mt19937_64 &engine, int repetitions) { } } } -} \ No newline at end of file +} diff --git a/Con2PrimFactory/configuration.ccl b/Con2PrimFactory/configuration.ccl index 92500653..b53e1dcc 100644 --- a/Con2PrimFactory/configuration.ccl +++ b/Con2PrimFactory/configuration.ccl @@ -1,5 +1,5 @@ #Configuration definitions for thorn Con2PrimFactory -REQUIRES Loop EOSX Algo +REQUIRES Loop EOSX Algo AsterUtils PROVIDES Con2PrimFactory { } diff --git a/Con2PrimFactory/interface.ccl b/Con2PrimFactory/interface.ccl index f6298a19..0470f717 100644 --- a/Con2PrimFactory/interface.ccl +++ b/Con2PrimFactory/interface.ccl @@ -9,3 +9,4 @@ INCLUDES HEADER: c2p_1DPalenzuela.hxx IN c2p_1DPalenzuela.hxx USES INCLUDE HEADER: eos.hxx USES INCLUDE HEADER: eos_idealgas.hxx USES INCLUDE HEADER: roots.hxx +USES INCLUDE HEADER: aster_utils.hxx diff --git a/Con2PrimFactory/param.ccl b/Con2PrimFactory/param.ccl index 6cec8390..36b2ec39 100644 --- a/Con2PrimFactory/param.ccl +++ b/Con2PrimFactory/param.ccl @@ -92,3 +92,24 @@ BOOLEAN thermal_eos_atmo "Whether to use the thermal (evolution) EOS for setting { } "no" +# Parameters for BH treatment + +CCTK_REAL alp_thresh "Lapse threshold below which primitives are modified for stability, specifically inside BH interior. Value should be set such that it is lesser than the lapse at the apparent horizon " STEERABLE=ALWAYS +{ + 0.0:* :: "" +} 0.0 + +CCTK_REAL rho_BH " Maximum density allowed in regions where alp #include #include -#include "c2p_utils.hxx" -#include "eos.hxx" -#include "eos_idealgas.hxx" - #include + #include "prims.hxx" #include "cons.hxx" #include "atmo.hxx" #include "c2p_report.hxx" +#include "eos.hxx" +#include "eos_idealgas.hxx" + +#include "c2p_utils.hxx" + namespace Con2PrimFactory { +using namespace AsterUtils; + constexpr CCTK_INT X = 0; constexpr CCTK_INT Y = 1; constexpr CCTK_INT Z = 2; diff --git a/Con2PrimFactory/src/c2p_1DPalenzuela.hxx b/Con2PrimFactory/src/c2p_1DPalenzuela.hxx index 44ba27b9..d06a866c 100644 --- a/Con2PrimFactory/src/c2p_1DPalenzuela.hxx +++ b/Con2PrimFactory/src/c2p_1DPalenzuela.hxx @@ -3,6 +3,7 @@ #include "c2p.hxx" #include "roots.hxx" + namespace Con2PrimFactory { class c2p_1DPalenzuela : public c2p { diff --git a/Con2PrimFactory/src/c2p_utils.hxx b/Con2PrimFactory/src/c2p_utils.hxx index a3e01208..bdde3085 100644 --- a/Con2PrimFactory/src/c2p_utils.hxx +++ b/Con2PrimFactory/src/c2p_utils.hxx @@ -1,6 +1,8 @@ #ifndef C2PUTILS_HXX #define C2PUTILS_HXX +#include + #include #include #include @@ -13,12 +15,10 @@ #include #include #include -#include + +#include "aster_utils.hxx" namespace Con2PrimFactory { -using namespace std; -using namespace Loop; -using namespace Arith; enum class ROOTSTAT { SUCCESS, @@ -26,348 +26,6 @@ enum class ROOTSTAT { NOT_BRACKETED }; -// Computes the contraction of smat and vec -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec -calc_contraction(const smat &g, const vec &v) { - return [&](int i) ARITH_INLINE { - return sum([&](int j) ARITH_INLINE { return g(i, j) * v(j); }); - }; -} - -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_contraction(const vec &v_up, const vec &v_dn) { - return sum([&](int i) ARITH_INLINE { return v_up(i) * v_dn(i); }); -} - -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_contraction(const smat &g_dn, const smat &T_up) { - return sum_symm([&](int i, int j) - ARITH_INLINE { return g_dn(i, j) * T_up(i, j); }); -} - -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_contraction(const smat &g, const vec &v1, - const vec &v2) { - // return calc_contraction(v2, calc_contraction(g, v1)); - return sum_symm([&](int i, int j) - ARITH_INLINE { return g(i, j) * v1(i) * v2(j); }); -} - -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec -calc_cross_product(const vec &B, const vec &v) { - return ([&](int i) ARITH_INLINE { - const int j = (i + 1) % 3, k = (i + 2) % 3; - return B(j) * v(k) - B(k) * v(j); - }); -} - -// Contraction for rc case: consider both sides of the face -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec, D> -calc_contraction(const smat &g, const vec, D> &v_rc) { - return ([&](int i) ARITH_INLINE { - return vec([&](int f) ARITH_INLINE { - return sum([&](int j) ARITH_INLINE { return g(i, j) * v_rc(j)(f); }); - }); - }); -} - -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec -calc_contraction(const vec, D> &vup_rc, - const vec, D> &vdn_rc) { - return ([&](int f) ARITH_INLINE { - return sum([&](int i) - ARITH_INLINE { return vup_rc(i)(f) * vdn_rc(i)(f); }); - }); -} - -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec, 3> -calc_cross_product(const vec, 3> &B_rc, - const vec, 3> &v_rc) { - return ([&](int i) ARITH_INLINE { - const int j = (i + 1) % 3, k = (i + 2) % 3; - return vec([&](int f) ARITH_INLINE { - return B_rc(j)(f) * v_rc(k)(f) - B_rc(k)(f) * v_rc(j)(f); - }); - }); -} - -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec, 3> -calc_cross_product(const vec &B, const vec, 3> &v_rc) { - return ([&](int i) ARITH_INLINE { - const int j = (i + 1) % 3, k = (i + 2) % 3; - return vec([&](int f) ARITH_INLINE { - return B(j) * v_rc(k)(f) - B(k) * v_rc(j)(f); - }); - }); -} - -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec -pow2(vec x) { - return ([&](int f) { return x(f) * x(f); }); -} - -// Computes the norm of vec, measured with smat -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_norm(const vec &v, const smat &g) { - return sum_symm([&](int i, int j) - ARITH_INLINE { return g(i, j) * v(i) * v(j); }); -} - -// Computes the transpose of vec -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec, F> -calc_transpose(const vec, D> &vv) { - return ([&](int f) ARITH_INLINE { - return vec([&](int d) ARITH_INLINE { return vv(d)(f); }); - }); -} - -// Computes the Lorentz factor -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_wlorentz(const vec &v_up, const vec &v_dn) { - return 1.0 / sqrt(1.0 - calc_contraction(v_up, v_dn)); -} - -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T pow2(T x) { - return x * x; -} - -// FD2: vertex centered input, vertex centered output, oneside stencil -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_fd2_v2v_oneside(const GF3D2 &gf, const PointDesc &p, - const int dir, const int sign) { - return -sign * - (gf(p.I + 2 * sign * p.DI[dir]) - 4.0 * gf(p.I + sign * p.DI[dir]) + - 3.0 * gf(p.I)) * - (0.5 / p.DX[dir]); -} - -// FD2: cell centered input, cell centered output -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_fd2_c2c(const GF3D2 &gf, const PointDesc &p, const int dir) { - return (0.5 / p.DX[dir]) * (gf(p.I + p.DI[dir]) - gf(p.I - p.DI[dir])); -} - -// FD2: vertex centered input, edge centered output -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_fd2_v2e(const GF3D2 &gf, const PointDesc &p, const int dir) { - return (gf(p.I + p.DI[dir]) - gf(p.I)) / p.DX[dir]; -} - -// FD2: vertex centered input, cell centered output -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_fd2_v2c(const GF3D2 &gf, const PointDesc &p, int dir) { - T dgf1, dgf2, dgf3, dgf4; - int dir1, dir2; - if (dir == 0) { - dir1 = 1; - dir2 = 2; - } else if (dir == 1) { - dir1 = 0; - dir2 = 2; - } else { - dir1 = 0; - dir2 = 1; - } - - dgf1 = (gf(p.I + p.DI[dir]) - gf(p.I)) / p.DX[dir]; - dgf2 = (gf(p.I + p.DI[dir1] + p.DI[dir]) - gf(p.I + p.DI[dir1])) / p.DX[dir]; - dgf3 = (gf(p.I + p.DI[dir2] + p.DI[dir]) - gf(p.I + p.DI[dir2])) / p.DX[dir]; - dgf4 = (gf(p.I + p.DI[dir1] + p.DI[dir2] + p.DI[dir]) - - gf(p.I + p.DI[dir1] + p.DI[dir2])) / - p.DX[dir]; - - return 0.25 * (dgf1 + dgf2 + dgf3 + dgf4); -} - -// FD4: cell centered input, cell centered output -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_fd4_c2c(const GF3D2 &gf, const PointDesc &p, const int dir) { - return (1.0 / (12.0 * p.DX[dir])) * - (-gf(p.I + 2 * p.DI[dir]) + 8.0 * gf(p.I + p.DI[dir]) - - 8.0 * gf(p.I - p.DI[dir]) + gf(p.I - 2 * p.DI[dir])); -} - -// FD4: vertex centered input, cell centered output -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_fd4_v2c(const GF3D2 &gf, const PointDesc &p, int dir) { - T dgf1, dgf2, dgf3, dgf4; - T dgf5, dgf6, dgf7, dgf8; - int dir1, dir2; - - if (dir == 0) { - dir1 = 1; - dir2 = 2; - } else if (dir == 1) { - dir1 = 0; - dir2 = 2; - } else { - dir1 = 0; - dir2 = 1; - } - - dgf1 = - ((1.0 / 24.0) * gf(p.I + 2 * p.DI[dir]) - (9.0 / 8.0) * gf(p.I + p.DI[dir]) + - (9.0 / 8.0) * gf(p.I) - (1.0 / 24.0) * gf(p.I - p.DI[dir])) / - p.DX[dir]; - dgf2 = ((1.0 / 24.0) * gf(p.I + p.DI[dir1] + 2 * p.DI[dir]) - - (9.0 / 8.0) * gf(p.I + p.DI[dir1] + p.DI[dir]) + - (9.0 / 8.0) * gf(p.I + p.DI[dir1]) - - (1.0 / 24.0) * gf(p.I + p.DI[dir1] - p.DI[dir])) / - p.DX[dir]; - dgf3 = ((1.0 / 24.0) * gf(p.I + p.DI[dir2] + 2 * p.DI[dir]) - - (9.0 / 8.0) * gf(p.I + p.DI[dir2] + p.DI[dir]) + - (9.0 / 8.0) * gf(p.I + p.DI[dir2]) - - (1.0 / 24.0) * gf(p.I + p.DI[dir2] - p.DI[dir])) / - p.DX[dir]; - dgf4 = ((1.0 / 24.0) * gf(p.I + p.DI[dir1] + p.DI[dir2] + 2 * p.DI[dir]) - - (9.0 / 8.0) * gf(p.I + p.DI[dir1] + p.DI[dir2] + p.DI[dir]) + - (9.0 / 8.0) * gf(p.I + p.DI[dir1] + p.DI[dir2]) - - (1.0 / 24.0) * gf(p.I + p.DI[dir1] + p.DI[dir2] - p.DI[dir])) / - p.DX[dir]; - - dgf5 = ((1.0 / 24.0) * gf(p.I + 2 * p.DI[dir1] + 2 * p.DI[dir2] + 2 * p.DI[dir]) - - (9.0 / 8.0) * gf(p.I + 2 * p.DI[dir1] + 2 * p.DI[dir2] + p.DI[dir]) + - (9.0 / 8.0) * gf(p.I + 2 * p.DI[dir1] + 2 * p.DI[dir2]) - - (1.0 / 24.0) * gf(p.I + 2 * p.DI[dir1] + 2 * p.DI[dir2] - p.DI[dir])) / - p.DX[dir]; - - dgf6 = ((1.0 / 24.0) * gf(p.I + 2 * p.DI[dir1] - p.DI[dir2] + 2 * p.DI[dir]) - - (9.0 / 8.0) * gf(p.I + 2 * p.DI[dir1] - p.DI[dir2] + p.DI[dir]) + - (9.0 / 8.0) * gf(p.I + 2 * p.DI[dir1] - p.DI[dir2]) - - (1.0 / 24.0) * gf(p.I + 2 * p.DI[dir1] - p.DI[dir2] - p.DI[dir])) / - p.DX[dir]; - - dgf7 = ((1.0 / 24.0) * gf(p.I - p.DI[dir1] + 2 * p.DI[dir2] + 2 * p.DI[dir]) - - (9.0 / 8.0) * gf(p.I - p.DI[dir1] + 2 * p.DI[dir2] + p.DI[dir]) + - (9.0 / 8.0) * gf(p.I - p.DI[dir1] + 2 * p.DI[dir2]) - - (1.0 / 24.0) * gf(p.I - p.DI[dir1] + 2 * p.DI[dir2] - p.DI[dir])) / - p.DX[dir]; - - dgf8 = ((1.0 / 24.0) * gf(p.I - p.DI[dir1] - p.DI[dir2] + 2 * p.DI[dir]) - - (9.0 / 8.0) * gf(p.I - p.DI[dir1] - p.DI[dir2] + p.DI[dir]) + - (9.0 / 8.0) * gf(p.I - p.DI[dir1] - p.DI[dir2]) - - (1.0 / 24.0) * gf(p.I - p.DI[dir1] - p.DI[dir2] - p.DI[dir])) / - p.DX[dir]; - - return (27.0 / 14.0) * (dgf1 + dgf2 + dgf3 + dgf4) - - (10.0 / 7.0) * (dgf5 + dgf6 + dgf7 + dgf8); -} - -// Second-order average of vertex-centered grid functions to cell center -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_avg_v2c(const GF3D2 &gf, const PointDesc &p) { - T gf_avg = 0.0; - - for (int dk = 0; dk < 2; ++dk) { - for (int dj = 0; dj < 2; ++dj) { - for (int di = 0; di < 2; ++di) { - gf_avg += gf(p.I + p.DI[0] * di + p.DI[1] * dj + p.DI[2] * dk); - } - } - } - return gf_avg / 8.0; -} - -// Second-order average of edge-centered grid functions to vertex-centered -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_avg_e2v(const GF3D2 &gf, const PointDesc &p, const int dir) { - T gf_avg = 0.0; - - for (int di = 0; di < 2; ++di) { - gf_avg += gf(p.I - p.DI[dir] * di); - } - return gf_avg / 2.0; -} - -// Second-order average of edge-centered grid functions (along dir) to cell -// center -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_avg_e2c(const GF3D2 &gf, const PointDesc &p, const int dir) { - T gf_avg = 0.0; - - for (int dk = 0; dk < (dir == 2 ? 1 : 2); ++dk) { - for (int dj = 0; dj < (dir == 1 ? 1 : 2); ++dj) { - for (int di = 0; di < (dir == 0 ? 1 : 2); ++di) { - gf_avg += gf(p.I + p.DI[0] * di + p.DI[1] * dj + p.DI[2] * dk); - } - } - } - return gf_avg / 4.0; -} - -// Second-order average of vertex-centered grid functionsto face -// center (perp to dir) -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_avg_v2f(const GF3D2 &gf, const PointDesc &p, const int dir) { - T gf_avg = 0.0; - - for (int dk = 0; dk < (dir == 2 ? 1 : 2); ++dk) { - for (int dj = 0; dj < (dir == 1 ? 1 : 2); ++dj) { - for (int di = 0; di < (dir == 0 ? 1 : 2); ++di) { - gf_avg += gf(p.I + p.DI[0] * di + p.DI[1] * dj + p.DI[2] * dk); - } - } - } - return gf_avg / 4.0; -} - -// Second-order average of cell-centered grid functions to edge center -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_avg_c2e(const GF3D2 &gf, const PointDesc &p, const int dir) { - T gf_avg = 0.0; - - for (int dk = 0; dk < (dir == 2 ? 1 : 2); ++dk) { - for (int dj = 0; dj < (dir == 1 ? 1 : 2); ++dj) { - for (int di = 0; di < (dir == 0 ? 1 : 2); ++di) { - gf_avg += gf(p.I - p.DI[0] * di - p.DI[1] * dj - p.DI[2] * dk); - } - } - } - return gf_avg / 4.0; -} - -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline vec -get_neighbors(const GF3D2 &gf, const PointDesc &p) { - return {gf(p.I - p.DI[0]), gf(p.I + p.DI[0]), gf(p.I - p.DI[1]), - gf(p.I + p.DI[1]), gf(p.I - p.DI[2]), gf(p.I + p.DI[2])}; -} - -template -CCTK_DEVICE CCTK_HOST CCTK_ATTRIBUTE_ALWAYS_INLINE inline T -calc_avg_neighbors(const vec flag, const vec u_nbs, - const vec u_saved_nbs) { - return sum([&](int i) ARITH_INLINE { - return (flag(i) * u_nbs(i) + (1.0 - flag(i)) * u_saved_nbs(i)); - }) / - CCTK_REAL(D); -} - } // namespace Con2PrimFactory #endif // #ifndef C2PUTILS_HXX diff --git a/Con2PrimFactory/src/cons.hxx b/Con2PrimFactory/src/cons.hxx index 43a6118f..9f0019bf 100644 --- a/Con2PrimFactory/src/cons.hxx +++ b/Con2PrimFactory/src/cons.hxx @@ -9,6 +9,8 @@ namespace Con2PrimFactory { +using namespace AsterUtils; + // Structure to represent conserved variables for mhd struct cons_vars { CCTK_REAL dens; diff --git a/Con2PrimFactory/src/prims.hxx b/Con2PrimFactory/src/prims.hxx index 6f4458b0..c39a26c7 100644 --- a/Con2PrimFactory/src/prims.hxx +++ b/Con2PrimFactory/src/prims.hxx @@ -9,6 +9,8 @@ namespace Con2PrimFactory { +using namespace AsterUtils; + struct prim_vars { CCTK_REAL rho; CCTK_REAL eps; diff --git a/Con2PrimFactory/src/test.cxx b/Con2PrimFactory/src/test.cxx index f91fd845..c80b592c 100644 --- a/Con2PrimFactory/src/test.cxx +++ b/Con2PrimFactory/src/test.cxx @@ -3,17 +3,20 @@ #include #include -#include "c2p_utils.hxx" +#include "eos.hxx" +#include "eos_idealgas.hxx" + #include "c2p.hxx" #include "c2p_1DPalenzuela.hxx" #include "c2p_2DNoble.hxx" -#include -#include + +#include "c2p_utils.hxx" namespace Con2PrimFactory { using namespace Arith; using namespace EOSX; +using namespace AsterUtils; extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; diff --git a/Docs/thornlist/asterx.th b/Docs/thornlist/asterx.th index d0066337..0603ec4a 100644 --- a/Docs/thornlist/asterx.th +++ b/Docs/thornlist/asterx.th @@ -274,6 +274,7 @@ SpacetimeX/Z4c !REPO_PATH= $2 !CHECKOUT = AsterX/AsterSeeds +AsterX/AsterUtils AsterX/AsterX AsterX/Con2PrimFactory AsterX/EOSX diff --git a/scripts/asterx.th b/scripts/asterx.th index 21a57dbf..fb70e3a4 100644 --- a/scripts/asterx.th +++ b/scripts/asterx.th @@ -743,6 +743,7 @@ SpacetimeX/Z4c !TYPE = ignore !CHECKOUT = AsterX/AsterSeeds +AsterX/AsterUtils AsterX/AsterX AsterX/Con2PrimFactory AsterX/EOSX