IENPJPUOBDWVTAFOKQTLF4GMEUABOUKCGHND73YRC3JZTWDZMDHQC
"# Needed functions and #include's:\n",
"USES INCLUDE: Symmetry.h\n",
"USES INCLUDE: Boundary.h\n",
"\n",
"# Needed Boundary Conditions function\n",
"CCTK_INT FUNCTION GetBoundarySpecification(CCTK_INT IN size, CCTK_INT OUT ARRAY nboundaryzones, CCTK_INT OUT ARRAY is_internal, CCTK_INT OUT ARRAY is_staggered, CCTK_INT OUT ARRAY shiftout)\n",
"USES FUNCTION GetBoundarySpecification\n",
"\n",
"CCTK_INT FUNCTION SymmetryTableHandleForGrid(CCTK_POINTER_TO_CONST IN cctkGH)\n",
"USES FUNCTION SymmetryTableHandleForGrid\n",
"\n",
"CCTK_INT FUNCTION Boundary_SelectVarForBC(CCTK_POINTER_TO_CONST IN GH, CCTK_INT IN faces, CCTK_INT IN boundary_width, CCTK_INT IN table_handle, CCTK_STRING IN var_name, CCTK_STRING IN bc_name)\n",
"USES FUNCTION Boundary_SelectVarForBC\n",
"\n",
"schedule BaikalETK_BoundaryConditions_evolved_gfs in MoL_PostStep\n",
"{\n",
" LANG: C\n",
" OPTIONS: LEVEL\n",
" SYNC: evol_variables\n",
"} \"Apply boundary conditions and perform AMR+interprocessor synchronization\"\n",
"\n",
"schedule GROUP ApplyBCs as BaikalETK_ApplyBCs in MoL_PostStep after BaikalETK_BoundaryConditions_evolved_gfs\n",
"{\n",
"} \"Group for applying boundary conditions\"\n",
"\n",
"\n",
"\n",
"schedule BaikalETK_BoundaryConditions_aux_gfs in MoL_PseudoEvolution after BaikalETK_BSSN_constraints\n",
"{\n",
" LANG: C\n",
" OPTIONS: LOCAL # Needed so that cctk_nghostzones[0] (the number of boundary points) is defined.\n",
" # In other words, don't use LEVEL mode here, or the number of boundary points\n",
" # filled may not match the actual number of ghost zones. Weird, huh?\n",
" SYNC: aux_variables\n",
"} \"Enforce symmetry BCs in constraint computation\"\n",
"\"\"\")\n",
" if add_stress_energy_source_terms == True:\n",
" file.write(\"\"\"\n",
"schedule BaikalETK_BSSN_to_ADM in MoL_PseudoEvolution after BaikalETK_BoundaryConditions_aux_gfs\n",
"{\n",
" LANG: C\n",
" OPTIONS: Local\n",
"} \"Perform BSSN-to-ADM conversion in MoL_PseudoEvolution. Needed for proper HydroBase integration.\"\n",
"\"\"\")\n",
" file.write(\"\"\"\n",
"schedule GROUP ApplyBCs as BaikalETK_auxgfs_ApplyBCs in MoL_PseudoEvolution after BaikalETK_BoundaryConditions_aux_gfs\n",
"{\n",
"} \"Apply boundary conditions\"\n",
" file.write(\"}\\n\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Next register symmetries\n",
"with open(append_to_make_code_defn_list(\"Symmetry_registration_oldCartGrid3D.c\"),\"w\") as file:\n",
" file.write(\"\"\"\n",
"#include \"cctk.h\"\n",
"#include \"cctk_Arguments.h\"\n",
"#include \"cctk_Parameters.h\"\n",
"#include \"Symmetry.h\"\n",
"\n",
"void BaikalETK_Symmetry_registration(CCTK_ARGUMENTS)\n",
"{\n",
" DECLARE_CCTK_ARGUMENTS;\n",
" DECLARE_CCTK_PARAMETERS;\n",
" \n",
" // Stores gridfunction parity across x=0, y=0, and z=0 planes, respectively\n",
" int sym[3];\n",
" \n",
" // Next register parities for each gridfunction based on its name \n",
" // (to ensure this algorithm is robust, gridfunctions with integers\n",
" // in their base names are forbidden in NRPy+).\n",
"\"\"\")\n",
" full_gfs_list = []\n",
" full_gfs_list.extend(evol_gfs_list)\n",
" full_gfs_list.extend(auxevol_gfs_list)\n",
" full_gfs_list.extend(aux_gfs_list)\n",
" for gf in full_gfs_list:\n",
" file.write(\"\"\"\n",
" // Default to scalar symmetry:\n",
" sym[0] = 1; sym[1] = 1; sym[2] = 1;\n",
" // Now modify sym[0], sym[1], and/or sym[2] as needed \n",
" // to account for gridfunction parity across \n",
" // x=0, y=0, and/or z=0 planes, respectively\n",
"\"\"\")\n",
" # If gridfunction name does not end in a digit, by NRPy+ syntax, it must be a scalar\n",
" if gf[len(gf) - 1].isdigit() == False:\n",
" pass # Scalar = default\n",
" elif len(gf) > 2:\n",
" # Rank-1 indexed expression (e.g., vector)\n",
" if gf[len(gf) - 2].isdigit() == False:\n",
" if int(gf[-1]) > 2:\n",
" print(\"Error: Found invalid gridfunction name: \"+gf)\n",
" sys.exit(1)\n",
" symidx = gf[-1]\n",
" file.write(\" sym[\"+symidx+\"] = -1;\\n\")\n",
" # Rank-2 indexed expression\n",
" elif gf[len(gf) - 2].isdigit() == True:\n",
" if len(gf) > 3 and gf[len(gf) - 3].isdigit() == True:\n",
" print(\"Error: Found a Rank-3 or above gridfunction: \"+gf+\", which is at the moment unsupported.\")\n",
" print(\"It should be easy to support this if desired.\")\n",
" sys.exit(1)\n",
" symidx0 = gf[-2]\n",
" file.write(\" sym[\"+symidx0+\"] *= -1;\\n\")\n",
" symidx1 = gf[-1]\n",
" file.write(\" sym[\"+symidx1+\"] *= -1;\\n\")\n",
" else:\n",
" print(\"Don't know how you got this far with a gridfunction named \"+gf+\", but I'll take no more of this nonsense.\")\n",
" print(\" Please follow best-practices and rename your gridfunction to be more descriptive\")\n",
" sys.exit(1)\n",
" file.write(\" SetCartSymVN(cctkGH, sym, \\\"BaikalETK::\"+gf+\"\\\");\\n\")\n",
" file.write(\"}\\n\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Next register with the boundary conditions thorns.\n",
"# PART 1: Set BC type to \"none\" for all variables\n",
"# Since we choose NewRad boundary conditions, we must register all \n",
"# gridfunctions to have boundary type \"none\". This is because\n",
"# NewRad is seen by the rest of the Toolkit as a modification to the\n",
"# RHSs.\n",
"\n",
"# This code is based on Kranc's McLachlan/ML_BSSN/src/Boundaries.cc code.\n",
"with open(append_to_make_code_defn_list(\"BoundaryConditions.c\"),\"w\") as file:\n",
" file.write(\"\"\"\n",
"#include \"cctk.h\"\n",
"#include \"cctk_Arguments.h\"\n",
"#include \"cctk_Parameters.h\"\n",
"#include \"cctk_Faces.h\"\n",
"#include \"util_Table.h\"\n",
"#include \"Symmetry.h\"\n",
"\n",
"// Set `none` boundary conditions on BSSN RHSs, as these are set via NewRad.\n",
"void BaikalETK_BoundaryConditions_evolved_gfs(CCTK_ARGUMENTS)\n",
"{\n",
" DECLARE_CCTK_ARGUMENTS;\n",
" DECLARE_CCTK_PARAMETERS;\n",
" \n",
" CCTK_INT ierr CCTK_ATTRIBUTE_UNUSED = 0;\n",
"\"\"\")\n",
" for gf in evol_gfs_list:\n",
" file.write(\"\"\"\n",
" ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, \"BaikalETK::\"\"\"+gf+\"\"\"\", \"none\");\n",
" if (ierr < 0) CCTK_ERROR(\"Failed to register BC for BaikalETK::\"\"\"+gf+\"\"\"!\");\n",
"\"\"\")\n",
" file.write(\"\"\"\n",
"}\n",
"\n",
"// Set `flat` boundary conditions on BSSN constraints, similar to what Lean does.\n",
"void BaikalETK_BoundaryConditions_aux_gfs(CCTK_ARGUMENTS) {\n",
" DECLARE_CCTK_ARGUMENTS;\n",
" DECLARE_CCTK_PARAMETERS;\n",
" \n",
" CCTK_INT ierr CCTK_ATTRIBUTE_UNUSED = 0;\n",
"\n",
"\"\"\")\n",
" for gf in aux_gfs_list:\n",
" file.write(\"\"\"\n",
" ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, cctk_nghostzones[0], -1, \"BaikalETK::\"\"\"+gf+\"\"\"\", \"flat\");\n",
" if (ierr < 0) CCTK_ERROR(\"Failed to register BC for BaikalETK::\"\"\"+gf+\"\"\"!\");\n",
"\"\"\")\n",
" file.write(\"}\\n\")\n",
" \n",
" \n",
"# PART 2: Set C code for calling NewRad BCs\n",
"# As explained in lean_public/LeanBSSNMoL/src/calc_bssn_rhs.F90,\n",
"# the function NewRad_Apply takes the following arguments:\n",
"# NewRad_Apply(cctkGH, var, rhs, var0, v0, radpower),\n",
"# which implement the boundary condition:\n",
"# var = var_at_infinite_r + u(r-var_char_speed*t)/r^var_radpower\n",
"# Obviously for var_radpower>0, var_at_infinite_r is the value of\n",
"# the variable at r->infinity. var_char_speed is the propagation\n",
"# speed at the outer boundary, and var_radpower is the radial \n",
"# falloff rate.\n",
"\n",
"with open(append_to_make_code_defn_list(\"BoundaryCondition_NewRad.c\"),\"w\") as file:\n",
" file.write(\"\"\"\n",
"#include <math.h>\n",
"\n",
"#include \"cctk.h\"\n",
"#include \"cctk_Arguments.h\"\n",
"#include \"cctk_Parameters.h\"\n",
"\n",
"void BaikalETK_NewRad(CCTK_ARGUMENTS) {\n",
" DECLARE_CCTK_ARGUMENTS;\n",
" DECLARE_CCTK_PARAMETERS;\n",
" \n",
"\"\"\")\n",
" for gf in evol_gfs_list:\n",
" var_at_infinite_r = \"0.0\"\n",
" var_char_speed = \"1.0\"\n",
" var_radpower = \"1.0\"\n",
" \n",
" if gf == \"alpha\":\n",
" var_at_infinite_r = \"1.0\"\n",
" if LapseCondition == \"OnePlusLog\":\n",
" var_char_speed = \"sqrt(2.0)\"\n",
" else:\n",
" pass # 1.0 (default) is fine\n",
" if \"aDD\" in gf or \"trK\" in gf: # consistent with Lean code.\n",
" var_radpower = \"2.0\"\n",
"\n",
" file.write(\" NewRad_Apply(cctkGH, \"+gf+\", \"+gf.replace(\"GF\",\"\")+\"_rhsGF, \"+var_at_infinite_r+\", \"+\n",
" var_char_speed+\", \"+var_radpower+\");\\n\")\n",