From 8b1cb9387844f5c6bd62427e32489bb9b2d69020 Mon Sep 17 00:00:00 2001 From: Hengrui Zhu Date: Sat, 16 Nov 2024 11:10:21 -0500 Subject: [PATCH 1/2] horizon find and cartesian grid dump --- src/CMakeLists.txt | 2 + src/athena.hpp | 5 +- src/driver/driver.cpp | 2 +- src/outputs/restart.cpp | 2 +- src/pgen/pgen.cpp | 2 +- src/tasklist/numerical_relativity.hpp | 1 + src/utils/cart_grid.cpp | 282 +++++++++++++++++++++++++ src/utils/cart_grid.hpp | 51 +++++ src/z4c/compact_object_tracker.hpp | 2 +- src/z4c/horizon_dump.cpp | 284 ++++++++++++++++++++++++++ src/z4c/horizon_dump.hpp | 62 ++++++ src/z4c/z4c.cpp | 41 +++- src/z4c/z4c.hpp | 63 +++--- src/z4c/z4c_amr.cpp | 28 +-- src/z4c/z4c_tasks.cpp | 33 ++- 15 files changed, 799 insertions(+), 61 deletions(-) create mode 100644 src/utils/cart_grid.cpp create mode 100644 src/utils/cart_grid.hpp create mode 100644 src/z4c/horizon_dump.cpp create mode 100644 src/z4c/horizon_dump.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index cace5efe..1a6113e0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -147,8 +147,10 @@ add_executable( utils/show_config.cpp utils/lagrange_interpolator.cpp utils/tr_table.cpp + utils/cart_grid.cpp z4c/compact_object_tracker.cpp + z4c/horizon_dump.cpp z4c/tmunu.cpp z4c/z4c.cpp z4c/z4c_adm.cpp diff --git a/src/athena.hpp b/src/athena.hpp index 6eca7e7c..aef84d84 100644 --- a/src/athena.hpp +++ b/src/athena.hpp @@ -130,6 +130,10 @@ template using DualArray2D = Kokkos::DualView; template using DualArray3D = Kokkos::DualView; +template +using DualArray4D = Kokkos::DualView; +template +using DualArray5D = Kokkos::DualView; // template declarations for construction of Kokkos::View in scratch memory template @@ -477,4 +481,3 @@ struct reduction_identity< array_sum::GlobalSum > { } #endif // ATHENA_HPP_ - diff --git a/src/driver/driver.cpp b/src/driver/driver.cpp index 6a29b189..3ae676de 100644 --- a/src/driver/driver.cpp +++ b/src/driver/driver.cpp @@ -258,7 +258,7 @@ Driver::Driver(ParameterInput *pin, Mesh *pmesh, Real wtlim, Kokkos::Timer* ptim } else { std::cout << "### FATAL ERROR in " << __FILE__ << " at line " << __LINE__ << std::endl << "integrator=" << integrator << " not implemented. " - << "Valid choices are [rk1,rk2,rk3,imex2,imex3]." << std::endl; + << "Valid choices are [rk1,rk2,rk3,rk4,imex2,imex3]." << std::endl; exit(EXIT_FAILURE); } } diff --git a/src/outputs/restart.cpp b/src/outputs/restart.cpp index 4914f69f..3ca3a6be 100644 --- a/src/outputs/restart.cpp +++ b/src/outputs/restart.cpp @@ -225,7 +225,7 @@ void RestartOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin) { // output puncture tracker data if (nco > 0) { for (auto & pt : pz4c->ptracker) { - resfile.Write_any_type(pt.GetPos(), 3*sizeof(Real), "byte"); + resfile.Write_any_type(pt->GetPos(), 3*sizeof(Real), "byte"); } } // turbulence driver internal RNG diff --git a/src/pgen/pgen.cpp b/src/pgen/pgen.cpp index 9dac7b0e..868e3e7c 100644 --- a/src/pgen/pgen.cpp +++ b/src/pgen/pgen.cpp @@ -202,7 +202,7 @@ ProblemGenerator::ProblemGenerator(ParameterInput *pin, Mesh *pm, IOWrapper resf #if MPI_PARALLEL_ENABLED MPI_Bcast(&pos[0], 3*sizeof(Real), MPI_CHAR, 0, MPI_COMM_WORLD); #endif - pt.SetPos(&pos[0]); + pt->SetPos(&pos[0]); } } diff --git a/src/tasklist/numerical_relativity.hpp b/src/tasklist/numerical_relativity.hpp index 18125490..4dac1266 100644 --- a/src/tasklist/numerical_relativity.hpp +++ b/src/tasklist/numerical_relativity.hpp @@ -82,6 +82,7 @@ enum TaskName { Z4c_ClearRW, Z4c_Wave, Z4c_PT, + Z4c_DumpHorizon, Z4c_NTASKS }; diff --git a/src/utils/cart_grid.cpp b/src/utils/cart_grid.cpp new file mode 100644 index 00000000..c62b6bc3 --- /dev/null +++ b/src/utils/cart_grid.cpp @@ -0,0 +1,282 @@ +//======================================================================================== +// AthenaXXX astrophysical plasma code +// Copyright(C) 2020 James M. Stone and the Athena code team +// Licensed under the 3-clause BSD License (the "LICENSE") +//======================================================================================== +//! \file cart_grid.cpp +// \brief Initializes a Cartesian grid to interpolate data onto + +// C/C++ headers +#include +#include +#include + +// AthenaK headers +#include "athena.hpp" +#include "coordinates/cell_locations.hpp" +#include "mesh/mesh.hpp" +#include "coordinates/coordinates.hpp" +#include "cart_grid.hpp" + +//---------------------------------------------------------------------------------------- +// constructor, initializes data structures and parameters + +CartesianGrid::CartesianGrid(MeshBlockPack *pmy_pack, Real center[3], + Real extend[3], int numpoints[3], bool is_cheb): + pmy_pack(pmy_pack), + interp_indcs("interp_indcs",1,1,1,1), + interp_wghts("interp_wghts",1,1,1,1,1), + interp_vals("interp_vals",1,1,1) { + // initialize parameters for the grid + // uniform grid or spectral grid + is_cheby = is_cheb; + + // grid center + center_x1 = center[0]; + center_x2 = center[1]; + center_x3 = center[2]; + + // grid center + extend_x1 = extend[0]; + extend_x2 = extend[1]; + extend_x3 = extend[2]; + + // lower bound + min_x1 = center_x1 - extend_x1; + min_x2 = center_x2 - extend_x2; + min_x3 = center_x3 - extend_x3; + + // upper bound + max_x1 = center_x1 + extend_x1; + max_x2 = center_x2 + extend_x2; + max_x3 = center_x3 + extend_x3; + + // number of points + nx1 = numpoints[0]; + nx2 = numpoints[1]; + nx3 = numpoints[2]; + + // resolution + d_x1 = (max_x1-min_x1)/(nx1-1); + d_x2 = (max_x2-min_x2)/(nx2-1); + d_x3 = (max_x3-min_x3)/(nx3-1); + + // allocate memory for interpolation coordinates, indices, and weights + int &ng = pmy_pack->pmesh->mb_indcs.ng; + Kokkos::realloc(interp_indcs,nx1,nx2,nx3,4); + Kokkos::realloc(interp_wghts,nx1,nx2,nx3,2*ng,3); + + // Call functions to prepare CartesianGrid object for interpolation + // SetInterpolationCoordinates(); + SetInterpolationIndices(); + SetInterpolationWeights(); + + return; +} + +CartesianGrid::~CartesianGrid() {} + +void CartesianGrid::ResetCenter(Real center[3]) { + // grid center + center_x1 = center[0]; + center_x2 = center[1]; + center_x3 = center[2]; + + // lower bound + min_x1 = center_x1 - extend_x1; + min_x2 = center_x2 - extend_x2; + min_x3 = center_x3 - extend_x3; + + // upper bound + max_x1 = center_x1 + extend_x1; + max_x2 = center_x2 + extend_x2; + max_x3 = center_x3 + extend_x3; + SetInterpolationIndices(); + SetInterpolationWeights(); +} + +void CartesianGrid::SetInterpolationIndices() { + auto &size = pmy_pack->pmb->mb_size; + int nmb1 = pmy_pack->nmb_thispack - 1; + + auto &iindcs = interp_indcs; + for (int nx=0; nx= x1min && x1 <= x1max) && + (x2 >= x2min && x2 <= x2max) && + (x3 >= x3min && x3 <= x3max)) { + iindcs.h_view(nx,ny,nz,0) = m; + iindcs.h_view(nx,ny,nz,1) = + static_cast(std::floor((x1-(x1min+dx1/2.0))/dx1)); + iindcs.h_view(nx,ny,nz,2) = + static_cast(std::floor((x2-(x2min+dx2/2.0))/dx2)); + iindcs.h_view(nx,ny,nz,3) = + static_cast(std::floor((x3-(x3min+dx3/2.0))/dx3)); + } + } + } + } + } + + // sync dual arrays + interp_indcs.template modify(); + interp_indcs.template sync(); + + return; +} + +void CartesianGrid::SetInterpolationWeights() { + auto &indcs = pmy_pack->pmesh->mb_indcs; + auto &size = pmy_pack->pmb->mb_size; + int &ng = indcs.ng; + + auto &iindcs = interp_indcs; + auto &iwghts = interp_wghts; + for (int nx=0; nx(); + interp_wghts.template sync(); + + return; +} + +//---------------------------------------------------------------------------------------- +//! \fn void CartesianGrid::InterpolateToGrid +//! \brief interpolate Cartesian data to cart_grid for output + +void CartesianGrid::InterpolateToGrid(int ind, DvceArray5D &val) { + // reinitialize interpolation indices and weights if AMR + //if (pmy_pack->pmesh->adaptive) { + // SetInterpolationIndices(); + // SetInterpolationWeights(); + //} + + // capturing variables for kernel + auto &indcs = pmy_pack->pmesh->mb_indcs; + int &is = indcs.is; int &js = indcs.js; int &ks = indcs.ks; + int &ng = indcs.ng; + int n_x1 = nx1 - 1; + int n_x2 = nx2 - 1; + int n_x3 = nx3 - 1; + int index = ind; + + // reallocate container + Kokkos::realloc(interp_vals,nx1,nx2,nx3); + + auto &iindcs = interp_indcs; + auto &iwghts = interp_wghts; + auto &ivals = interp_vals; + par_for("int2cart",DevExeSpace(),0,n_x1,0,n_x2,0,n_x3, + KOKKOS_LAMBDA(int nx, int ny, int nz) { + int ii0 = iindcs.d_view(nx,ny,nz,0); + int ii1 = iindcs.d_view(nx,ny,nz,1); + int ii2 = iindcs.d_view(nx,ny,nz,2); + int ii3 = iindcs.d_view(nx,ny,nz,3); + if (ii0==-1) { // point not on this rank + ivals.d_view(nx,ny,nz) = 0.0; + } else { + Real int_value = 0.0; + for (int i=0; i<2*ng; i++) { + for (int j=0; j<2*ng; j++) { + for (int k=0; k<2*ng; k++) { + Real iwght = iwghts.d_view(nx,ny,nz,i,0)* + iwghts.d_view(nx,ny,nz,j,1)*iwghts.d_view(nx,ny,nz,k,2); + int_value += iwght*val(ii0,index,ii3-(ng-k-ks)+1, + ii2-(ng-j-js)+1,ii1-(ng-i-is)+1); + } + } + } + ivals.d_view(nx,ny,nz) = int_value; + } + }); + + // sync dual arrays + interp_vals.template modify(); + interp_vals.template sync(); + + return; +} + diff --git a/src/utils/cart_grid.hpp b/src/utils/cart_grid.hpp new file mode 100644 index 00000000..01c53d23 --- /dev/null +++ b/src/utils/cart_grid.hpp @@ -0,0 +1,51 @@ +#ifndef UTILS_CART_GRID_HPP_ +#define UTILS_CART_GRID_HPP_ + +//======================================================================================== +// AthenaXXX astrophysical plasma code +// Copyright(C) 2020 James M. Stone and the Athena code team +// Licensed under the 3-clause BSD License (the "LICENSE") +//======================================================================================== +//! \file cart_grid.hpp +// \brief definitions for SphericalGrid class + +#include "athena.hpp" + +// Forward declarations +class MeshBlockPack; + +//---------------------------------------------------------------------------------------- +//! \class CartesianGrid + +class CartesianGrid { + public: + // Creates a geodesic grid with refinement level nlev and radius rad + CartesianGrid(MeshBlockPack *pmy_pack, Real center[3], + Real extend[3], int numpoints[3], bool is_cheb = false); + ~CartesianGrid(); + + // parameters for the grid + Real center_x1, center_x2, center_x3; // grid centers + Real min_x1, min_x2, min_x3; // min for xyz + Real max_x1, max_x2, max_x3; // max value for xyz + Real d_x1, d_x2, d_x3; // resolution + int nx1, nx2, nx3; // number of points + Real extend_x1, extend_x2, extend_x3; + + // dump on chebyshev or uniform grid, default is uniform + bool is_cheby; + + // For simplicity, unravell all points into a 1d array + DualArray3D interp_vals; // container for data interpolated to sphere + void InterpolateToGrid(int nvars, DvceArray5D &val); // interpolate to sphere + void ResetCenter(Real center[3]); // set indexing for interpolation + + private: + MeshBlockPack* pmy_pack; // ptr to MeshBlockPack containing this Hydro + DualArray4D interp_indcs; // indices of MeshBlock and zones therein for interp + DualArray5D interp_wghts; // weights for interpolation + void SetInterpolationIndices(); // set indexing for interpolation + void SetInterpolationWeights(); // set weights for interpolation +}; + +#endif // UTILS_CART_GRID_HPP_ diff --git a/src/z4c/compact_object_tracker.hpp b/src/z4c/compact_object_tracker.hpp index 0dfa7db3..d681db42 100644 --- a/src/z4c/compact_object_tracker.hpp +++ b/src/z4c/compact_object_tracker.hpp @@ -55,11 +55,11 @@ class CompactObjectTracker { inline Real GetRadius() const { return radius; } + Real pos[NDIM]; private: bool owns_compact_object; CompactObjectType type; - Real pos[NDIM]; Real vel[NDIM]; int reflevel; // requested minimum refinement level (-1 for infinity) Real radius; // nominal radius of the object (for the AMR driver) diff --git a/src/z4c/horizon_dump.cpp b/src/z4c/horizon_dump.cpp new file mode 100644 index 00000000..9e61a1d2 --- /dev/null +++ b/src/z4c/horizon_dump.cpp @@ -0,0 +1,284 @@ +//======================================================================================== +// AthenaXXX astrophysical plasma code +// Copyright(C) 2020 James M. Stone and the Athena code team +// Licensed under the 3-clause BSD License (the "LICENSE") +//======================================================================================== + +#include +#include +#include // mkdir + +#include +#include +#include +#include +#include +#include +#include +#include + +#if MPI_PARALLEL_ENABLED +#include +#endif + +#include "horizon_dump.hpp" +#include "athena.hpp" +#include "globals.hpp" +#include "mesh/mesh.hpp" +#include "parameter_input.hpp" +#include "utils/cart_grid.hpp" +#include "coordinates/adm.hpp" +#include "mhd/mhd.hpp" +#include "z4c/z4c.hpp" + +//---------------------------------------------------------------------------------------- +HorizonDump::HorizonDump(MeshBlockPack *pmbp, ParameterInput *pin, int n, int is_common): + common_horizon{is_common}, pos{NAN, NAN, NAN}, + pmbp{pmbp}, horizon_ind{n} { + std::string nstr = std::to_string(n); + + pos[0] = pin->GetOrAddReal("z4c", "co_" + nstr + "_x", 0.0); + pos[1] = pin->GetOrAddReal("z4c", "co_" + nstr + "_y", 0.0); + pos[2] = pin->GetOrAddReal("z4c", "co_" + nstr + "_z", 0.0); + + horizon_extent = pin->GetOrAddReal("z4c", "co_" + nstr + "_dump_radius", 2.0); + horizon_nx = pin->GetOrAddInteger("z4c", "horizon_" + + std::to_string(n)+"_Nx",10); + horizon_dt = pin->GetOrAddReal("z4c", "horizon_dt", 1.0); + output_count = 0; + + Real extend[3] = {horizon_extent,horizon_extent,horizon_extent}; + int Nx[3] = {horizon_nx,horizon_nx,horizon_nx}; + pcat_grid = new CartesianGrid(pmbp, pos, extend, Nx); + + // Initializing variables that will be dumped + // The order is alpha, betax, betay, betaz, + // gxx, gxy, gxz, gyy, gyz, gzz + // Kxx, Kxy, Kxz, Kyy, Kyz, Kzz + variable_to_dump.push_back(std::make_pair(pmbp->pz4c->I_Z4C_ALPHA, true)); + variable_to_dump.push_back(std::make_pair(pmbp->pz4c->I_Z4C_BETAX, true)); + variable_to_dump.push_back(std::make_pair(pmbp->pz4c->I_Z4C_BETAY, true)); + variable_to_dump.push_back(std::make_pair(pmbp->pz4c->I_Z4C_BETAZ, true)); + variable_to_dump.push_back(std::make_pair(pmbp->padm->I_ADM_GXX, false)); + variable_to_dump.push_back(std::make_pair(pmbp->padm->I_ADM_GXY, false)); + variable_to_dump.push_back(std::make_pair(pmbp->padm->I_ADM_GXZ, false)); + variable_to_dump.push_back(std::make_pair(pmbp->padm->I_ADM_GYY, false)); + variable_to_dump.push_back(std::make_pair(pmbp->padm->I_ADM_GYZ, false)); + variable_to_dump.push_back(std::make_pair(pmbp->padm->I_ADM_GZZ, false)); + variable_to_dump.push_back(std::make_pair(pmbp->padm->I_ADM_KXX, false)); + variable_to_dump.push_back(std::make_pair(pmbp->padm->I_ADM_KXY, false)); + variable_to_dump.push_back(std::make_pair(pmbp->padm->I_ADM_KXZ, false)); + variable_to_dump.push_back(std::make_pair(pmbp->padm->I_ADM_KYY, false)); + variable_to_dump.push_back(std::make_pair(pmbp->padm->I_ADM_KYZ, false)); + variable_to_dump.push_back(std::make_pair(pmbp->padm->I_ADM_KZZ, false)); +} + +//---------------------------------------------------------------------------------------- +HorizonDump::~HorizonDump() { + delete pcat_grid; +} + +void HorizonDump::SetGridAndInterpolate(Real center[NDIM]) { + // update center location + pos[0] = center[0]; + pos[1] = center[1]; + pos[2] = center[2]; + + pcat_grid->ResetCenter(pos); + // Real* data_out = new Real []; + // swap out to 1d array + // Real data_out[horizon_nx][horizon_nx][horizon_nx][16]; + + // Define the size of each dimension + int count = horizon_nx * horizon_nx * horizon_nx * 16; + + // Dynamically allocate memory for the 4D array flattened into 1D + Real* data_out = new Real[count]; + + for(int nvar=0; nvar<16; nvar++) { + // Interpolate here + if (variable_to_dump[nvar].second) { + pcat_grid->InterpolateToGrid(variable_to_dump[nvar].first,pmbp->pz4c->u0); + } else { + pcat_grid->InterpolateToGrid(variable_to_dump[nvar].first,pmbp->padm->u_adm); + } + for (int nx = 0; nx < horizon_nx; nx ++) + for (int ny = 0; ny < horizon_nx; ny ++) + for (int nz = 0; nz < horizon_nx; nz ++) { + data_out[nvar * horizon_nx * horizon_nx * horizon_nx + // Section for nvar + nx * horizon_nx * horizon_nx + // Slice for nx + ny * horizon_nx + // Row for ny + nz] // Column for nz + = pcat_grid->interp_vals.h_view(nx, ny, nz); // Value being assigned + } + } + + // MPI reduce here + // Reduction to the master rank for data_out + #if MPI_PARALLEL_ENABLED + if (0 == global_variable::my_rank) { + MPI_Reduce(MPI_IN_PLACE, data_out, count, + MPI_ATHENA_REAL, MPI_SUM, 0, MPI_COMM_WORLD); + } else { + MPI_Reduce(data_out, data_out, count, MPI_ATHENA_REAL, MPI_SUM, 0, MPI_COMM_WORLD); + } + #endif + // Then write output file + // Open the file in binary write mode + std::string foldername = "horizon_"+std::to_string(horizon_ind) + +"/output_"+std::to_string(output_count); + mkdir(foldername.c_str(),0775); + + if (0 == global_variable::my_rank) { + std::string fname = foldername + "/etk_output_file.dat"; + FILE* etk_output_file = fopen(fname.c_str(), "wb"); + if (etk_output_file == nullptr) { + perror("Error opening file"); + return; + } + fwrite(&common_horizon, sizeof(int), 1, etk_output_file); + fwrite(&pmbp->pmesh->time, sizeof(Real), 1, etk_output_file); + // Write the 4D array to the binary file + size_t elementsWritten = fwrite(data_out, sizeof(Real), count, etk_output_file); + if (elementsWritten != count) { + perror("Error writing to file"); + } + // Close the file + fclose(etk_output_file); + + char parfilename[100]; + ETK_setup_parfile(0.5, parfilename); + output_count++; + // delete dataout + delete[] data_out; + } +} + +void HorizonDump::ETK_setup_parfile(const Real BH_radius_guess, char parfilename[100]) { + std::string foldername; + if (common_horizon == 0) { + foldername = "horizon_"+std::to_string(horizon_ind)+ + "/output_"+std::to_string(output_count); + } else { + foldername = "horizon_common/output_"+std::to_string(output_count); + } + std::string fname = foldername + "/ET_analyze_BHaH_data_horizon.par"; + + sprintf(parfilename, "%s", fname.c_str()); // NOLINT + + FILE *etk_parfile = fopen(parfilename, "w"); + fprintf(etk_parfile, + "ActiveThorns = \"PUGH SymBase CartGrid3D\"\n" + "cactus::cctk_itlast = 0\n" + "#cactus::cctk_show_schedule = \"yes\" # //" + "Disables initial scheduler printout.\n" + "cactus::cctk_show_schedule = \"no\" # //" + "Disables initial scheduler printout.\n" + "cactus::cctk_show_banners = \"no\" # // Disables banners.\n" + "Driver::ghost_size = 0\n" + "Driver::global_nsize = %d\n" + "Driver::info = load\n" + "grid::type = byrange\n" + "\n" + "grid::xmin = %e\n" + "grid::xmax = %e\n" + "grid::ymin = %e\n" + "grid::ymax = %e\n" + "grid::zmin = %e\n" + "grid::zmax = %e\n" + "ActiveThorns = ADMBase\n" + "#ActiveThorns = \"AHFinderDirect SphericalSurface SpaceMask StaticConformal" + " IOUtil AEILocalInterp PUGHInterp PUGHReduce QuasiLocalMeasures IOBasic" + " TmunuBase ADMCoupling ADMMacros LocalReduce\"\n" + "ActiveThorns = \"AHFinderDirect SphericalSurface SpaceMask StaticConformal" + " IOUtil AEILocalInterp PUGHInterp PUGHReduce QuasiLocalMeasures IOBasic" + " TmunuBase LocalReduce\"\n" + "ActiveThorns = \"readBHaHdata\"\n" + "ADMBase::metric_type = \"physical\"\n" + "AHFinderDirect::find_every = 1\n" + "AHFinderDirect::geometry_interpolator_name =" + " \"Lagrange polynomial interpolation\"\n" + "AHFinderDirect::geometry_interpolator_pars = \"order=4\"\n" + "AHFinderDirect::max_Newton_iterations__initial = 100\n" + "AHFinderDirect::max_Newton_iterations__subsequent = 10\n" + "AHFinderDirect::N_horizons = 1\n" + "AHFinderDirect::output_BH_diagnostics = \"yes\"\n" + "AHFinderDirect::reset_horizon_after_not_finding[1] = \"no\"\n" + "AHFinderDirect::set_mask_for_individual_horizon[1] = \"no\"\n" + "AHFinderDirect::surface_interpolator_name =" + " \"Lagrange polynomial interpolation\"\n" + "AHFinderDirect::surface_interpolator_pars = \"order=4\"\n" + "AHFinderDirect::verbose_level = \"physics details\"\n" + "#AHFinderDirect::verbose_level =" + " \"algorithm details\"\n" + "AHFinderDirect::which_surface_to_store_info[1] = 0\n" + "AHFinderDirect::run_at_CCTK_POSTSTEP = false\n" + "AHFinderDirect::run_at_CCTK_ANALYSIS = true\n" + "\n" + "# Parameters of thorn QuasiLocalMeasures (implementing QuasiLocalMeasures)\n" + "QuasiLocalMeasures::interpolator =" + " \"Lagrange polynomial interpolation\"\n" + "QuasiLocalMeasures::interpolator_options = \"order=4\"\n" + "# QuasiLocalMeasures::interpolator =" + " \"Hermite polynomial interpolation\"\n" + "# QuasiLocalMeasures::interpolator_options = \"order=3\"\n" + "QuasiLocalMeasures::num_surfaces = 1\n" + "QuasiLocalMeasures::spatial_order = 2\n" + "QuasiLocalMeasures::surface_index[0] = 0\n" + "QuasiLocalMeasures::verbose = yes\n" + "#QuasiLocalMeasures::veryverbose = yes\n" + "SphericalSurface::nsurfaces = 1\n" + "# You may find benefit using super high SphericalSurface " + "resolutions with very high spin BHs\n" + "# SphericalSurface::maxntheta = 301\n" + "# SphericalSurface::maxnphi = 504\n" + "# SphericalSurface::ntheta [0] = 301\n" + "# SphericalSurface::nphi [0] = 504\n" + "SphericalSurface::maxntheta = 161\n" + "SphericalSurface::maxnphi = 324\n" + "SphericalSurface::ntheta [0] = 161\n" + "SphericalSurface::nphi [0] = 324\n" + "SphericalSurface::nghoststheta[0] = 2\n" + "SphericalSurface::nghostsphi [0] = 2\n" + "# SphericalSurface::set_spherical[1]= yes\n" + "# SphericalSurface::radius [1]= 40\n" + "# SphericalSurface::radius [2]= 80\n" + "IOBasic::outInfo_every = 1\n" + "IOBasic::outInfo_vars = \"\n" + " QuasiLocalMeasures::qlm_scalars\n" + " QuasiLocalMeasures::qlm_spin[0]\n" + " QuasiLocalMeasures::qlm_radius[0]\n" + " QuasiLocalMeasures::qlm_mass[0]\n" + " QuasiLocalMeasures::qlm_3det[0] \"\n", horizon_nx, + -horizon_extent, horizon_extent, -horizon_extent, horizon_extent, + -horizon_extent, horizon_extent); + + if(horizon_ind==0) + fprintf(etk_parfile, + "IOUtil::out_dir = \"AHET_out_horizon_BH_0_ahf_ihf_diags\"\n" + "readBHaHdata::outfilename = \"horizon_BH_0_ahf_ihf_diags.txt\"\n" + "readBHaHdata::recent_ah_radius_max_filename = \"ah_radius_max_BH_0.txt\"\n"); + if(horizon_ind==1) + fprintf(etk_parfile, + "IOUtil::out_dir = \"AHET_out_horizon_BH_1_ahf_ihf_diags\"\n" + "readBHaHdata::outfilename = \"horizon_BH_1_ahf_ihf_diags.txt\"\n" + "readBHaHdata::recent_ah_radius_max_filename = \"ah_radius_max_BH_1.txt\"\n"); + + //if(commondata->time == 0) { + fprintf(etk_parfile, + "AHFinderDirect::initial_guess_method[1]" + " = \"coordinate sphere\"\n" + "AHFinderDirect::initial_guess__coord_sphere__radius[1] = %e\n", + BH_radius_guess); + /*} else { + if(BH==LESSMASSIVE_BH_PREMERGER) + fprintf(etk_parfile, + "AHFinderDirect::initial_guess_method[1] = \"read from named file\"\n" + "AHFinderDirect::initial_guess__read_from_named_file__file_name[1] = \"AHET_out_horizon_BH_m_ahf_ihf_diags/latest_ah_surface.gp\"\n"); + if(BH==MOREMASSIVE_BH_PREMERGER || BH==BH_POSTMERGER) + fprintf(etk_parfile, + "AHFinderDirect::initial_guess_method[1] = \"read from named file\"\n" + "AHFinderDirect::initial_guess__read_from_named_file__file_name[1] = \"AHET_out_horizon_BH_M_ahf_ihf_diags/latest_ah_surface.gp\"\n"); + }*/ + fclose(etk_parfile); +} diff --git a/src/z4c/horizon_dump.hpp b/src/z4c/horizon_dump.hpp new file mode 100644 index 00000000..bcde47b0 --- /dev/null +++ b/src/z4c/horizon_dump.hpp @@ -0,0 +1,62 @@ +//======================================================================================== +// AthenaXXX astrophysical plasma code +// Copyright(C) 2020 James M. Stone and the Athena code team +// Licensed under the 3-clause BSD License (the "LICENSE") +//======================================================================================== + +#ifndef Z4C_HORIZON_DUMP_HPP_ +#define Z4C_HORIZON_DUMP_HPP_ + +#include +#include +#include +#include +#include + +#include "athena.hpp" +#include "mesh/mesh.hpp" +#include "z4c_macros.hpp" + +// Forward declaration +class Mesh; +class MeshBlockPack; +class ParameterInput; +class CartesianGrid; + +//! \class CompactObjectTracker +//! \brief Tracks a single puncture +class HorizonDump { + public: + //! Initialize a tracker + HorizonDump(MeshBlockPack *pmbp, ParameterInput *pin, int n, int common_horizon); + //! Destructor (will close output file) + ~HorizonDump(); + + // Interpolate field to Cartesian Grid centered at the puncture locations + void SetGridAndInterpolate(Real center[NDIM]); + + //! Write parameter file for Einstein Toolkit + void ETK_setup_parfile(const Real BH_radius_guess, char parfilename[100]); + + int horizon_nx; // number of points in each direction + int common_horizon; // common horizon or not, triggering when to start dumping data + int horizon_ind; // indices for horizon + // TODO(hzhu) : check if this works with rst + int output_count; // counting the output number (for naming subfolders) + + Real horizon_dt; + Real horizon_last_output_time; + Real horizon_extent; // radius for dumping data in a cube + CartesianGrid *pcat_grid=nullptr; // pointer to cartesian grid + Real pos[NDIM]; // position of the puncture + // Real data[horizon_nx][horizon_nx][horizon_nx][16]; + + private: + MeshBlockPack const *pmbp; + Real radius; // nominal radius of the object (for the AMR driver) + // first element store variable index second + // store whether from z4c (true) or adm (false) array + std::vector> variable_to_dump; +}; + +#endif // Z4C_HORIZON_DUMP_HPP_ diff --git a/src/z4c/z4c.cpp b/src/z4c/z4c.cpp index ad5bb5e2..8edbdc3d 100644 --- a/src/z4c/z4c.cpp +++ b/src/z4c/z4c.cpp @@ -21,9 +21,11 @@ #include "mesh/mesh.hpp" #include "bvals/bvals.hpp" #include "z4c/compact_object_tracker.hpp" +#include "z4c/horizon_dump.hpp" #include "z4c/z4c.hpp" #include "z4c/z4c_amr.hpp" #include "coordinates/adm.hpp" +#include "utils/cart_grid.hpp" namespace z4c { @@ -120,6 +122,7 @@ Z4c::Z4c(MeshBlockPack *ppack, ParameterInput *pin) : opt.chi_psi_power = pin->GetOrAddReal("z4c", "chi_psi_power", -4.0); opt.chi_div_floor = pin->GetOrAddReal("z4c", "chi_div_floor", -1000.0); + opt.chi_min_floor = pin->GetOrAddReal("z4c", "chi_min_floor", 1e-12); opt.diss = pin->GetOrAddReal("z4c", "diss", 0.0); opt.eps_floor = pin->GetOrAddReal("z4c", "eps_floor", 1e-12); opt.damp_kappa1 = pin->GetOrAddReal("z4c", "damp_kappa1", 0.0); @@ -129,9 +132,13 @@ Z4c::Z4c(MeshBlockPack *ppack, ParameterInput *pin) : opt.lapse_harmonic = pin->GetOrAddReal("z4c", "lapse_harmonic", 0.0); opt.lapse_oplog = pin->GetOrAddReal("z4c", "lapse_oplog", 2.0); opt.lapse_advect = pin->GetOrAddReal("z4c", "lapse_advect", 1.0); + opt.slow_start_lapse = pin->GetOrAddBoolean("z4c", "slow_start_lapse", false); + opt.ssl_damping_amp = pin->GetOrAddReal("z4c", "ssl_damping_amp", 0.6); + opt.ssl_damping_time = pin->GetOrAddReal("z4c", "ssl_damping_time", 20.0); + opt.ssl_damping_index = pin->GetOrAddInteger("z4c", "ssl_damping_index", 1); + opt.shift_ggamma = pin->GetOrAddReal("z4c", "shift_Gamma", 1.0); opt.shift_advect = pin->GetOrAddReal("z4c", "shift_advect", 1.0); - opt.shift_alpha2ggamma = pin->GetOrAddReal("z4c", "shift_alpha2Gamma", 0.0); opt.shift_hh = pin->GetOrAddReal("z4c", "shift_H", 0.0); @@ -186,12 +193,42 @@ Z4c::Z4c(MeshBlockPack *ppack, ParameterInput *pin) : int n = 0; while (true) { if (pin->DoesParameterExist("z4c", "co_" + std::to_string(n) + "_type")) { - ptracker.emplace_back(pmy_pack->pmesh, pin, n); + ptracker.push_back(std::make_unique(pmy_pack->pmesh, pin, n)); n++; } else { break; } } + // Construct the Cartesian data grid for dumping horizon data + n = 0; + while (true) { + if (pin->GetOrAddBoolean("z4c", "dump_horizon_" + std::to_string(n),false)) { + // phorizon_dump.emplace_back(pmy_pack, pin, n,false); + phorizon_dump.push_back(std::make_unique(pmy_pack, pin, n, 0)); + std::string foldername = "horizon_"+std::to_string(n); + mkdir(foldername.c_str(),0775); + n++; + } else { + break; + } + } + /* + horizon_dt = pin->GetOrAddReal("z4c", "horizon_dt", 1); + horizon_last_output_time = 0; + n = 0; + for (auto & pt : ptracker) { + if (pin->GetOrAddBoolean("z4c", "dump_horizon_" + std::to_string(n),false)) { + horizon_extent.push_back(pin->GetOrAddReal("z4c", "horizon_" + + std::to_string(n)+"_radius",3)); + Real extend[3] = {horizon_extent[n],horizon_extent[n],horizon_extent[n]}; + horizon_nx.push_back(pin->GetOrAddInteger("z4c", "horizon_" + + std::to_string(n)+"_Nx",100)); + int Nx[3] = {horizon_nx[n],horizon_nx[n],horizon_nx[n]}; + horizon_dump.emplace_back(pmy_pack, pt.pos, extend, Nx); + n++; + } + } + */ } //---------------------------------------------------------------------------------------- diff --git a/src/z4c/z4c.hpp b/src/z4c/z4c.hpp index 928a8b61..a97d20d6 100644 --- a/src/z4c/z4c.hpp +++ b/src/z4c/z4c.hpp @@ -15,6 +15,7 @@ #include #include "athena.hpp" #include "utils/finite_diff.hpp" +#include "utils/cart_grid.hpp" #include "parameter_input.hpp" #include "tasklist/task_list.hpp" #include "bvals/bvals.hpp" @@ -26,39 +27,7 @@ class Coordinates; class Driver; class CompactObjectTracker; - -//---------------------------------------------------------------------------------------- -//! \struct Z4cTaskIDs -// \brief container to hold TaskIDs of all z4c tasks - -struct Z4cTaskIDs { - TaskID irecv; - TaskID irecvweyl; - TaskID copyu; - TaskID crhs; - TaskID sombc; - TaskID expl; - TaskID sendu; - TaskID recvu; - TaskID newdt; - TaskID bcs; - TaskID prol; - TaskID algc; - TaskID z4tad; - TaskID admc; - TaskID csend; - TaskID crecv; - TaskID restu; - TaskID ptrck; - TaskID weyl_scalar; - TaskID wave_extr; - TaskID weyl_rest; - TaskID weyl_send; - TaskID weyl_prol; - TaskID weyl_recv; - TaskID csendweyl; - TaskID crecvweyl; -}; +class HorizonDump; namespace z4c { class Z4c_AMR; @@ -173,6 +142,8 @@ class Z4c { // puncture's floor value for chi, use max(chi, chi_div_floor) // in non-differentiated chi Real chi_div_floor; + Real chi_min_floor; // minimum of chi, only used in slow-start-lapse + // where a square root is necessary. Real diss; // amount of numerical dissipation Real eps_floor; // a small number O(10^-12) // Constraint damping parameters @@ -183,12 +154,20 @@ class Z4c { Real lapse_harmonicf; Real lapse_harmonic; Real lapse_advect; + // slow start lapse condition + bool slow_start_lapse; + Real ssl_damping_amp; + Real ssl_damping_time; + Real ssl_damping_index; // Gauge condition for the shift Real shift_ggamma; Real shift_alpha2ggamma; Real shift_hh; Real shift_advect; Real shift_eta; + // turn on shift damping smoothly + bool slow_roll_eta; + Real turn_on_time; // Enable BSSN if false (disable theta) bool use_z4c; // Apply the Sommerfeld condition for user BCs. @@ -207,8 +186,6 @@ class Z4c { // following only used for time-evolving flow Real dtnew; - // container to hold names of TaskIDs - Z4cTaskIDs id; // geodesic grid for wave extr std::vector> spherical_grids; @@ -218,6 +195,8 @@ class Z4c { Real last_output_time; int nrad; // number of radii to perform wave extraction + // dump data cube at horizon + // functions void QueueZ4cTasks(); TaskStatus InitRecv(Driver *d, int stage); @@ -247,6 +226,7 @@ class Z4c { TaskStatus TrackCompactObjects(Driver *d, int stage); TaskStatus CalcWeylScalar(Driver *d, int stage); TaskStatus CalcWaveForm(Driver *d, int stage); + TaskStatus DumpHorizons(Driver *d, int stage); template TaskStatus CalcRHS(Driver *d, int stage); @@ -262,8 +242,19 @@ class Z4c { void AlgConstr(MeshBlockPack *pmbp); Z4c_AMR *pamr; - std::list ptracker; + std::vector> ptracker; + std::vector> phorizon_dump; + /* + std::list horizon_dump; + Real horizon_dt; + Real horizon_last_output_time; + std::vector horizon_extent; // radius for dumping data in a cube + std::vector horizon_nx; // number of points in each direction + */ + // TODO(@hzhu): think about how to automatically trigger common horizon + // maybe have a horizon dump object to save all the space here + // same for the waveform. private: MeshBlockPack* pmy_pack; // ptr to MeshBlockPack containing this Z4c }; diff --git a/src/z4c/z4c_amr.cpp b/src/z4c/z4c_amr.cpp index 7efac0a7..ac18c513 100644 --- a/src/z4c/z4c_amr.cpp +++ b/src/z4c/z4c_amr.cpp @@ -93,25 +93,25 @@ void Z4c_AMR::RefineTracker(MeshBlockPack *pmbp) { flag.clear(); for (auto & pt : pmbp->pz4c->ptracker) { Real d2[8] = { - SQ(x1min - pt.GetPos(0)) + SQ(x2min - pt.GetPos(1)) + SQ(x3min - pt.GetPos(2)), - SQ(x1max - pt.GetPos(0)) + SQ(x2min - pt.GetPos(1)) + SQ(x3min - pt.GetPos(2)), - SQ(x1min - pt.GetPos(0)) + SQ(x2max - pt.GetPos(1)) + SQ(x3min - pt.GetPos(2)), - SQ(x1max - pt.GetPos(0)) + SQ(x2max - pt.GetPos(1)) + SQ(x3min - pt.GetPos(2)), - SQ(x1min - pt.GetPos(0)) + SQ(x2min - pt.GetPos(1)) + SQ(x3max - pt.GetPos(2)), - SQ(x1max - pt.GetPos(0)) + SQ(x2min - pt.GetPos(1)) + SQ(x3max - pt.GetPos(2)), - SQ(x1min - pt.GetPos(0)) + SQ(x2max - pt.GetPos(1)) + SQ(x3max - pt.GetPos(2)), - SQ(x1max - pt.GetPos(0)) + SQ(x2max - pt.GetPos(1)) + SQ(x3max - pt.GetPos(2)), + SQ(x1min - pt->GetPos(0)) + SQ(x2min - pt->GetPos(1)) + SQ(x3min - pt->GetPos(2)), + SQ(x1max - pt->GetPos(0)) + SQ(x2min - pt->GetPos(1)) + SQ(x3min - pt->GetPos(2)), + SQ(x1min - pt->GetPos(0)) + SQ(x2max - pt->GetPos(1)) + SQ(x3min - pt->GetPos(2)), + SQ(x1max - pt->GetPos(0)) + SQ(x2max - pt->GetPos(1)) + SQ(x3min - pt->GetPos(2)), + SQ(x1min - pt->GetPos(0)) + SQ(x2min - pt->GetPos(1)) + SQ(x3max - pt->GetPos(2)), + SQ(x1max - pt->GetPos(0)) + SQ(x2min - pt->GetPos(1)) + SQ(x3max - pt->GetPos(2)), + SQ(x1min - pt->GetPos(0)) + SQ(x2max - pt->GetPos(1)) + SQ(x3max - pt->GetPos(2)), + SQ(x1max - pt->GetPos(0)) + SQ(x2max - pt->GetPos(1)) + SQ(x3max - pt->GetPos(2)), }; Real dmin2 = *std::min_element(&d2[0], &d2[8]); bool iscontained = - (pt.GetPos(0) >= x1min && pt.GetPos(0) <= x1max) && - (pt.GetPos(1) >= x2min && pt.GetPos(1) <= x2max) && - (pt.GetPos(2) >= x3min && pt.GetPos(2) <= x3max); + (pt->GetPos(0) >= x1min && pt->GetPos(0) <= x1max) && + (pt->GetPos(1) >= x2min && pt->GetPos(1) <= x2max) && + (pt->GetPos(2) >= x3min && pt->GetPos(2) <= x3max); - if (dmin2 < SQ(pt.GetRadius()) || iscontained) { - if (pt.GetReflevel() < 0 || level < pt.GetReflevel()) { + if (dmin2 < SQ(pt->GetRadius()) || iscontained) { + if (pt->GetReflevel() < 0 || level < pt->GetReflevel()) { flag.push_back(1); - } else if (level == pt.GetReflevel()) { + } else if (level == pt->GetReflevel()) { flag.push_back(0); } else { flag.push_back(-1); diff --git a/src/z4c/z4c_tasks.cpp b/src/z4c/z4c_tasks.cpp index 551cbcde..69384fbf 100644 --- a/src/z4c/z4c_tasks.cpp +++ b/src/z4c/z4c_tasks.cpp @@ -19,6 +19,7 @@ #include "mesh/mesh.hpp" #include "bvals/bvals.hpp" #include "z4c/compact_object_tracker.hpp" +#include "z4c/horizon_dump.hpp" #include "z4c/z4c.hpp" #include "tasklist/numerical_relativity.hpp" @@ -28,7 +29,7 @@ namespace z4c { //! \fn void Z4c::QueueZ4cTasks //! \brief queue Z4c tasks into NumericalRelativity void Z4c::QueueZ4cTasks() { - Kokkos::printf("AssembleZ4cTasks\n"); + printf("AssembleZ4cTasks\n"); using namespace mhd; // NOLINT(build/namespaces) using namespace numrel; // NOLINT(build/namespaces) NumericalRelativity *pnr = pmy_pack->pnr; @@ -95,6 +96,8 @@ void Z4c::QueueZ4cTasks() { pnr->QueueTask(&Z4c::CalcWaveForm, this, Z4c_Wave, "Z4c_Wave", Task_End, {Z4c_ClearRW}); pnr->QueueTask(&Z4c::TrackCompactObjects, this, Z4c_PT, "Z4c_PT", Task_End, {Z4c_Wave}); + pnr->QueueTask(&Z4c::DumpHorizons, this, Z4c_DumpHorizon, "Z4c_DumpHorizon", + Task_End, {Z4c_PT}); } //---------------------------------------------------------------------------------------- @@ -281,9 +284,9 @@ TaskStatus Z4c::ApplyPhysicalBCs(Driver *pdrive, int stage) { TaskStatus Z4c::TrackCompactObjects(Driver *pdrive, int stage) { if (stage == pdrive->nexp_stages) { for (auto & pt : ptracker) { - pt.InterpolateVelocity(pmy_pack); - pt.EvolveTracker(); - pt.WriteTracker(); + pt->InterpolateVelocity(pmy_pack); + pt->EvolveTracker(); + pt->WriteTracker(); } } return TaskStatus::complete; @@ -456,4 +459,26 @@ TaskStatus Z4c::ClearSendWeyl(Driver *pdrive, int stage) { } } +TaskStatus Z4c::DumpHorizons(Driver *pdrive, int stage) { + if (pmy_pack->pz4c->phorizon_dump.size() == 0 || stage != pdrive->nexp_stages) { + return TaskStatus::complete; + } else { + float time_32 = static_cast(pmy_pack->pmesh->time); + float next_32 = static_cast(pmy_pack->pz4c->phorizon_dump[0] + ->horizon_last_output_time + +pmy_pack->pz4c->phorizon_dump[0]->horizon_dt); + if (((time_32 >= next_32) || (time_32 == 0))) { + int i = 0; + for (auto & hd : phorizon_dump) { + hd->horizon_last_output_time = time_32; + hd->SetGridAndInterpolate(pmy_pack->pz4c->ptracker[i]->pos); + i++; + } + } + return TaskStatus::complete; + } + + return TaskStatus::complete; +} + } // namespace z4c From 58935fead49548c9f56bffbbc40a52dde19daa4e Mon Sep 17 00:00:00 2001 From: Hengrui Zhu Date: Mon, 25 Nov 2024 15:52:16 -0500 Subject: [PATCH 2/2] addressing comments and add kerr schild initial data --- src/pgen/gr_analytic/kerr_schild.cpp | 125 +++++++++++++++++++++++++++ src/utils/cart_grid.cpp | 6 +- src/z4c/compact_object_tracker.hpp | 2 +- src/z4c/horizon_dump.cpp | 15 ++-- src/z4c/horizon_dump.hpp | 5 +- src/z4c/z4c_tasks.cpp | 2 +- 6 files changed, 139 insertions(+), 16 deletions(-) create mode 100644 src/pgen/gr_analytic/kerr_schild.cpp diff --git a/src/pgen/gr_analytic/kerr_schild.cpp b/src/pgen/gr_analytic/kerr_schild.cpp new file mode 100644 index 00000000..da53557f --- /dev/null +++ b/src/pgen/gr_analytic/kerr_schild.cpp @@ -0,0 +1,125 @@ +//======================================================================================== +// AthenaXXX astrophysical plasma code +// Copyright(C) 2020 James M. Stone and the Athena code team +// Licensed under the 3-clause BSD License (the "LICENSE") +//======================================================================================== +//! \file z4c_one_puncture.cpp +// \brief Problem generator for a single puncture placed at the origin of the domain + +#include +#include +#include +#include +#include // endl +#include // numeric_limits::max() +#include +#include // c_str(), string +#include + +#include "athena.hpp" +#include "parameter_input.hpp" +#include "globals.hpp" +#include "mesh/mesh.hpp" +#include "z4c/z4c.hpp" +#include "z4c/z4c_amr.hpp" +#include "coordinates/adm.hpp" +#include "coordinates/cell_locations.hpp" +#include "coordinates/cartesian_ks.hpp" + + +void KerrSchild(MeshBlockPack *pmbp, ParameterInput *pin); +void RefinementCondition(MeshBlockPack* pmbp); +KOKKOS_INLINE_FUNCTION +bool invertSymmetricMatrix4x4(const double g_dd[4][4], double g_dd_inv[4][4]); + +//---------------------------------------------------------------------------------------- +//! \fn ProblemGenerator::UserProblem_() +//! \brief Problem Generator for single puncture +void ProblemGenerator::UserProblem(ParameterInput *pin, const bool restart) { + user_ref_func = RefinementCondition; + MeshBlockPack *pmbp = pmy_mesh_->pmb_pack; + auto &indcs = pmy_mesh_->mb_indcs; + + if (pmbp->pz4c == nullptr) { + std::cout << "### FATAL ERROR in " << __FILE__ << " at line " << __LINE__ << std::endl + << "One Puncture test can only be run in Z4c, but no block " + << "in input file" << std::endl; + exit(EXIT_FAILURE); + } + + KerrSchild(pmbp, pin); + switch (indcs.ng) { + case 2: pmbp->pz4c->ADMToZ4c<2>(pmbp, pin); + break; + case 3: pmbp->pz4c->ADMToZ4c<3>(pmbp, pin); + break; + case 4: pmbp->pz4c->ADMToZ4c<4>(pmbp, pin); + break; + } + pmbp->pz4c->Z4cToADM(pmbp); + pmbp->pz4c->GaugePreCollapsedLapse(pmbp, pin); + switch (indcs.ng) { + case 2: pmbp->pz4c->ADMConstraints<2>(pmbp); + break; + case 3: pmbp->pz4c->ADMConstraints<3>(pmbp); + break; + case 4: pmbp->pz4c->ADMConstraints<4>(pmbp); + break; + } + std::cout<<"Kerr Schild initialized."<GetOrAddReal("problem", "punc_ADM_mass", 1.); + Real center_x1 = pin->GetOrAddReal("problem", "punc_center_x1", 0.); + Real center_x2 = pin->GetOrAddReal("problem", "punc_center_x2", 0.); + Real center_x3 = pin->GetOrAddReal("problem", "punc_center_x3", 0.); + Real a = 0; + bool minkowski = false; + + // capture variables for the kernel + auto &indcs = pmbp->pmesh->mb_indcs; + auto &size = pmbp->pmb->mb_size; + int &is = indcs.is; int &ie = indcs.ie; + int &js = indcs.js; int &je = indcs.je; + int &ks = indcs.ks; int &ke = indcs.ke; + int nmb = pmbp->nmb_thispack; + auto &adm = pmbp->padm->adm; + int &ng = indcs.ng; + int n1 = indcs.nx1 + 2*ng; + int n2 = indcs.nx2 + 2*ng; + int n3 = indcs.nx3 + 2*ng; + par_for("pgen_adm_vars", DevExeSpace(), 0,nmb-1,0,(n3-1),0,(n2-1),0,(n1-1), + KOKKOS_LAMBDA(int m, int k, int j, int i) { + Real &x1min = size.d_view(m).x1min; + Real &x1max = size.d_view(m).x1max; + Real x1v = CellCenterX(i-is, indcs.nx1, x1min, x1max); + + Real &x2min = size.d_view(m).x2min; + Real &x2max = size.d_view(m).x2max; + Real x2v = CellCenterX(j-js, indcs.nx2, x2min, x2max); + + Real &x3min = size.d_view(m).x3min; + Real &x3max = size.d_view(m).x3max; + Real x3v = CellCenterX(k-ks, indcs.nx3, x3min, x3max); + + ComputeADMDecomposition(x1v, x2v, x3v, minkowski, a, + &adm.alpha(m,k,j,i), + &adm.beta_u(m,0,k,j,i), &adm.beta_u(m,1,k,j,i), &adm.beta_u(m,2,k,j,i), + &adm.psi4(m,k,j,i), + &adm.g_dd(m,0,0,k,j,i), &adm.g_dd(m,0,1,k,j,i), &adm.g_dd(m,0,2,k,j,i), + &adm.g_dd(m,1,1,k,j,i), &adm.g_dd(m,1,2,k,j,i), &adm.g_dd(m,2,2,k,j,i), + &adm.vK_dd(m,0,0,k,j,i), &adm.vK_dd(m,0,1,k,j,i), &adm.vK_dd(m,0,2,k,j,i), + &adm.vK_dd(m,1,1,k,j,i), &adm.vK_dd(m,1,2,k,j,i), &adm.vK_dd(m,2,2,k,j,i)); + }); +} + +// how decide the refinement +void RefinementCondition(MeshBlockPack* pmbp) { + pmbp->pz4c->pamr->Refine(pmbp); +} diff --git a/src/utils/cart_grid.cpp b/src/utils/cart_grid.cpp index c62b6bc3..cb89cdca 100644 --- a/src/utils/cart_grid.cpp +++ b/src/utils/cart_grid.cpp @@ -108,9 +108,9 @@ void CartesianGrid::SetInterpolationIndices() { Real x2 = min_x2 + ny * d_x2; Real x3 = min_x3 + nz * d_x3; if (is_cheby) { - x1 = center_x1 + extend_x1*std::cos((2*nx+1)*M_PI/(2*nx1)); - x2 = center_x2 + extend_x2*std::cos((2*ny+1)*M_PI/(2*nx2)); - x3 = center_x3 + extend_x3*std::cos((2*nz+1)*M_PI/(2*nx3)); + x1 = center_x1 + extend_x1*std::cos(nx*M_PI/(nx1-1)); + x2 = center_x2 + extend_x2*std::cos(ny*M_PI/(nx2-1)); + x3 = center_x3 + extend_x3*std::cos(nz*M_PI/(nx3-1)); } // indices default to -1 if point does not reside in this MeshBlockPack iindcs.h_view(nx,ny,nz,0) = -1; diff --git a/src/z4c/compact_object_tracker.hpp b/src/z4c/compact_object_tracker.hpp index d681db42..61f20bc9 100644 --- a/src/z4c/compact_object_tracker.hpp +++ b/src/z4c/compact_object_tracker.hpp @@ -55,7 +55,6 @@ class CompactObjectTracker { inline Real GetRadius() const { return radius; } - Real pos[NDIM]; private: bool owns_compact_object; @@ -66,6 +65,7 @@ class CompactObjectTracker { Mesh const *pmesh; int out_every; std::ofstream ofile; + Real pos[NDIM]; }; #endif // Z4C_COMPACT_OBJECT_TRACKER_HPP_ diff --git a/src/z4c/horizon_dump.cpp b/src/z4c/horizon_dump.cpp index 9e61a1d2..2286168c 100644 --- a/src/z4c/horizon_dump.cpp +++ b/src/z4c/horizon_dump.cpp @@ -43,8 +43,9 @@ HorizonDump::HorizonDump(MeshBlockPack *pmbp, ParameterInput *pin, int n, int is horizon_extent = pin->GetOrAddReal("z4c", "co_" + nstr + "_dump_radius", 2.0); horizon_nx = pin->GetOrAddInteger("z4c", "horizon_" - + std::to_string(n)+"_Nx",10); + + nstr+"_Nx",10); horizon_dt = pin->GetOrAddReal("z4c", "horizon_dt", 1.0); + r_guess = pin->GetOrAddReal("z4c", "horizon" + nstr + "r_guess", 0.5); output_count = 0; Real extend[3] = {horizon_extent,horizon_extent,horizon_extent}; @@ -146,15 +147,15 @@ void HorizonDump::SetGridAndInterpolate(Real center[NDIM]) { // Close the file fclose(etk_output_file); - char parfilename[100]; - ETK_setup_parfile(0.5, parfilename); + // Write input script for Einstein Toolkit + ETK_setup_parfile(); output_count++; // delete dataout delete[] data_out; } } -void HorizonDump::ETK_setup_parfile(const Real BH_radius_guess, char parfilename[100]) { +void HorizonDump::ETK_setup_parfile() { std::string foldername; if (common_horizon == 0) { foldername = "horizon_"+std::to_string(horizon_ind)+ @@ -164,9 +165,7 @@ void HorizonDump::ETK_setup_parfile(const Real BH_radius_guess, char parfilename } std::string fname = foldername + "/ET_analyze_BHaH_data_horizon.par"; - sprintf(parfilename, "%s", fname.c_str()); // NOLINT - - FILE *etk_parfile = fopen(parfilename, "w"); + FILE *etk_parfile = fopen(fname.c_str(), "w"); fprintf(etk_parfile, "ActiveThorns = \"PUGH SymBase CartGrid3D\"\n" "cactus::cctk_itlast = 0\n" @@ -269,7 +268,7 @@ void HorizonDump::ETK_setup_parfile(const Real BH_radius_guess, char parfilename "AHFinderDirect::initial_guess_method[1]" " = \"coordinate sphere\"\n" "AHFinderDirect::initial_guess__coord_sphere__radius[1] = %e\n", - BH_radius_guess); + r_guess); /*} else { if(BH==LESSMASSIVE_BH_PREMERGER) fprintf(etk_parfile, diff --git a/src/z4c/horizon_dump.hpp b/src/z4c/horizon_dump.hpp index bcde47b0..384b7f1e 100644 --- a/src/z4c/horizon_dump.hpp +++ b/src/z4c/horizon_dump.hpp @@ -36,7 +36,7 @@ class HorizonDump { void SetGridAndInterpolate(Real center[NDIM]); //! Write parameter file for Einstein Toolkit - void ETK_setup_parfile(const Real BH_radius_guess, char parfilename[100]); + void ETK_setup_parfile(); int horizon_nx; // number of points in each direction int common_horizon; // common horizon or not, triggering when to start dumping data @@ -49,11 +49,10 @@ class HorizonDump { Real horizon_extent; // radius for dumping data in a cube CartesianGrid *pcat_grid=nullptr; // pointer to cartesian grid Real pos[NDIM]; // position of the puncture - // Real data[horizon_nx][horizon_nx][horizon_nx][16]; + Real r_guess; // nominal radius of the object (for the AMR driver) private: MeshBlockPack const *pmbp; - Real radius; // nominal radius of the object (for the AMR driver) // first element store variable index second // store whether from z4c (true) or adm (false) array std::vector> variable_to_dump; diff --git a/src/z4c/z4c_tasks.cpp b/src/z4c/z4c_tasks.cpp index 69384fbf..db50221f 100644 --- a/src/z4c/z4c_tasks.cpp +++ b/src/z4c/z4c_tasks.cpp @@ -471,7 +471,7 @@ TaskStatus Z4c::DumpHorizons(Driver *pdrive, int stage) { int i = 0; for (auto & hd : phorizon_dump) { hd->horizon_last_output_time = time_32; - hd->SetGridAndInterpolate(pmy_pack->pz4c->ptracker[i]->pos); + hd->SetGridAndInterpolate(pmy_pack->pz4c->ptracker[i]->GetPos()); i++; } }