Dumping data on worldtube for Cauchy-Characteristic Extraction (#618)

* initiate cce.

* Add CCE dump task.

* Ctors for cce dump.

* Add cce compilation.

* Minor.

* Implement cce interpolation.

* Declare cce dump task.

* Fix declaration.

* Add hdf5 path.

* Fix compilation bugs.

* Fix bugs.

* Add dt param for dump freq.

* Comment out debugs.

* Add macros for enabling/disabling cce dump.

* Add compiling option for cce and hdf5 lib.

* Fix macro bug.

* Add header.

* Make SortByIndex and SortByRank a local struct to prevent linking error.

* Fix a leack.

* Add timer for profiling.

* Fix timer bug.

* init script for converting Athenak CCE output to Specter.

* add gauss-legendre infrastructure

* add gitignore for output

* fix compilation error

* now compiles

* template for new cce infra

* factor out hdf5 and clean up template

* design reader and data struct to save data.

* add sketch of Fourier time derivative.

* more time derive.

* near complete version

* minor.

* impl: fft time derivs.

* fix: fft.

* style.

* test: is time derivative working?

* chore: add todo and more gitignore.

* chore: rm test output.

* test: f derives.

* perf: optimize FFT

* horizon find and cartesian grid dump

* feat: start adding radial derivs.

* fix: complete cheb. derives and add more test.

* feat: add Cheb U interpolation at R=r.

* perf: provide support for mp.

* feat: more mp.

* style.

* fix: mp.

* feat: add h5 write.

* fix: write
chore: add todo.

* doc: try to figure out spec h5 format.

* add: data attrs.

* doc: add comments.

* fix: bug in radial derivs & max_l.

* feat: fin. diff time derives.

* fix: bugs.

* addressing comments and add kerr schild initial data

* running version on cpu

* working version

* delete unused variables

* fix linting and remove redundant files

* revert to an earlier working version

* fix bug

* fix bug

* leave out post analysis scripts for now

* fix lint

* correct version of kokkos

* restore unintentional files to main

---------

Co-authored-by: Alireza <[email protected]>
Co-authored-by: Hengrui Zhu <[email protected]>
Co-authored-by: Hengrui Zhu <[email protected]>
Co-authored-by: Hengrui Zhu <[email protected]>

.gitignore

@@ -1,5 +1,5 @@
 # build directories
-build/
+build*/
 
 # Athena output files
 *.tab

src/CMakeLists.txt

@@ -56,6 +56,7 @@ add_executable(
 
         geodesic-grid/geodesic_grid.cpp
         geodesic-grid/spherical_grid.cpp
+	geodesic-grid/gauss_legendre.cpp
 
         hydro/hydro.cpp
         hydro/hydro_fluxes.cpp
@@ -163,6 +164,7 @@ add_executable(
         z4c/z4c_calculate_weyl_scalars.cpp
         z4c/z4c_wave_extr.cpp
         z4c/z4c_amr.cpp
+        z4c/cce/cce.cpp
 )
 
 # custom problem generator to be included in compile

src/geodesic-grid/gauss_legendre.cpp

@@ -0,0 +1,263 @@
+//========================================================================================
+// AthenaXXX astrophysical plasma code
+// Copyright(C) 2020 James M. Stone <[email protected]> and the Athena code team
+// Licensed under the 3-clause BSD License (the "LICENSE")
+//========================================================================================
+//! \file gauss_legendre.cpp
+//  \brief Initializes a Gauss-Legendra grid to interpolate data onto
+
+#include <cmath>
+#include <iostream>
+#include <list>
+
+#include "athena.hpp"
+#include "coordinates/cell_locations.hpp"
+#include "mesh/mesh.hpp"
+#include "hydro/hydro.hpp"
+#include "mhd/mhd.hpp"
+#include "coordinates/coordinates.hpp"
+#include "gauss_legendre.hpp"
+#include "utils/spherical_harm.hpp"
+#include "utils/legendre_roots.hpp"
+
+//----------------------------------------------------------------------------------------
+// constructor, initializes data structures and parameters
+
+GaussLegendreGrid::GaussLegendreGrid(MeshBlockPack *pmy_pack, int ntheta, Real rad):
+  pmy_pack(pmy_pack),
+  radius(rad),
+  ntheta(ntheta),
+  int_weights("int_weights",1),
+  polar_pos("polar_pos",1,1),
+  cart_pos("cart_pos",1,1),
+  interp_indcs("interp_indcs",1,1),
+  interp_wghts("interp_wghts",1,1,1),
+  interp_vals("interp_vals",1) {
+  // reallocate and set interpolation coordinates, indices, and weights
+  int &ng = pmy_pack->pmesh->mb_indcs.ng;
+  nangles = 2*ntheta*ntheta;
+
+  Kokkos::realloc(int_weights,nangles);
+  Kokkos::realloc(polar_pos,nangles,2);
+  Kokkos::realloc(cart_pos,nangles,3);
+  Kokkos::realloc(interp_indcs,nangles,4);
+  Kokkos::realloc(interp_wghts,nangles,2*ng,3);
+
+  InitializeAngleAndWeights();
+  InitializeRadius();
+  SetInterpolationIndices();
+  SetInterpolationWeights();
+  return;
+}
+
+//----------------------------------------------------------------------------------------
+//! \brief GaussLegendreGrid destructor
+
+GaussLegendreGrid::~GaussLegendreGrid() {
+}
+
+void GaussLegendreGrid::InitializeAngleAndWeights() {
+  // calculate roots and weights for the legendre polynomial
+  auto roots_and_weights = RootsAndWeights(ntheta);
+
+  // order for storing angle, first theta, then phi
+  for (int n=0; n<nangles; ++n) {
+    // save the weights
+    int_weights.h_view(n) =  roots_and_weights[1][n%ntheta]*M_PI/ntheta;
+    // calculate theta
+    polar_pos.h_view(n,0) = acos(roots_and_weights[0][n%ntheta]);
+    // calculate phi
+    polar_pos.h_view(n,1) = 2*M_PI/(2*ntheta)*(n/ntheta);
+  }
+
+  // sync to device
+  polar_pos.template modify<HostMemSpace>();
+  polar_pos.template sync<DevExeSpace>();
+
+  int_weights.template modify<HostMemSpace>();
+  int_weights.template sync<DevExeSpace>();
+}
+
+void GaussLegendreGrid::InitializeRadius() {
+  for (int n=0; n<nangles; ++n) {
+    Real &theta = polar_pos.h_view(n,0);
+    Real &phi = polar_pos.h_view(n,1);
+    cart_pos.h_view(n,0) = radius*cos(phi)*sin(theta);
+    cart_pos.h_view(n,1) = radius*sin(phi)*sin(theta);
+    cart_pos.h_view(n,2) = radius*cos(theta);
+  }
+  cart_pos.template modify<HostMemSpace>();
+  cart_pos.template sync<DevExeSpace>();
+}
+
+//----------------------------------------------------------------------------------------
+//! \fn void GaussLegendreGrid::SetInterpolationIndices
+//! \brief determine which MeshBlocks and MeshBlock zones therein that will be used in
+//         interpolation onto the sphere
+
+void GaussLegendreGrid::SetInterpolationIndices() {
+  auto &size = pmy_pack->pmb->mb_size;
+
+  int nmb1 = pmy_pack->nmb_thispack - 1;
+  int nang1 = nangles - 1;
+  auto &rcoord = cart_pos;
+  auto &iindcs = interp_indcs;
+  for (int n=0; n<=nang1; ++n) {
+    // indices default to -1 if angle does not reside in this MeshBlockPack
+    iindcs.h_view(n,0) = -1;
+    iindcs.h_view(n,1) = -1;
+    iindcs.h_view(n,2) = -1;
+    iindcs.h_view(n,3) = -1;
+    for (int m=0; m<=nmb1; ++m) {
+      // extract MeshBlock bounds
+      Real &x1min = size.h_view(m).x1min;
+      Real &x1max = size.h_view(m).x1max;
+      Real &x2min = size.h_view(m).x2min;
+      Real &x2max = size.h_view(m).x2max;
+      Real &x3min = size.h_view(m).x3min;
+      Real &x3max = size.h_view(m).x3max;
+
+      // extract MeshBlock grid cell spacings
+      Real &dx1 = size.h_view(m).dx1;
+      Real &dx2 = size.h_view(m).dx2;
+      Real &dx3 = size.h_view(m).dx3;
+
+      // save MeshBlock and zone indicies for nearest position to spherical patch center
+      // if this angle position resides in this MeshBlock
+      if ((rcoord.h_view(n,0) >= x1min && rcoord.h_view(n,0) <= x1max) &&
+          (rcoord.h_view(n,1) >= x2min && rcoord.h_view(n,1) <= x2max) &&
+          (rcoord.h_view(n,2) >= x3min && rcoord.h_view(n,2) <= x3max)) {
+        iindcs.h_view(n,0) = m;
+        iindcs.h_view(n,1) = static_cast<int>(std::floor((rcoord.h_view(n,0)-
+                                                          (x1min+dx1/2.0))/dx1));
+        iindcs.h_view(n,2) = static_cast<int>(std::floor((rcoord.h_view(n,1)-
+                                                          (x2min+dx2/2.0))/dx2));
+        iindcs.h_view(n,3) = static_cast<int>(std::floor((rcoord.h_view(n,2)-
+                                                          (x3min+dx3/2.0))/dx3));
+      }
+    }
+  }
+
+  // sync dual arrays
+  interp_indcs.template modify<HostMemSpace>();
+  interp_indcs.template sync<DevExeSpace>();
+
+  return;
+}
+
+//----------------------------------------------------------------------------------------
+//! \fn void GaussLegendreGrid::SetInterpolationWeights
+//! \brief set weights used by Lagrangian interpolation
+
+void GaussLegendreGrid::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 n=0; n<nangles; ++n) {
+    // extract indices
+    int &ii0 = iindcs.h_view(n,0);
+    int &ii1 = iindcs.h_view(n,1);
+    int &ii2 = iindcs.h_view(n,2);
+    int &ii3 = iindcs.h_view(n,3);
+
+    if (ii0==-1) {  // angle not on this rank
+      for (int i=0; i<2*ng; ++i) {
+        iwghts.h_view(n,i,0) = 0.0;
+        iwghts.h_view(n,i,1) = 0.0;
+        iwghts.h_view(n,i,2) = 0.0;
+      }
+    } else {
+      // extract spherical grid positions
+      Real &x0 = cart_pos.h_view(n,0);
+      Real &y0 = cart_pos.h_view(n,1);
+      Real &z0 = cart_pos.h_view(n,2);
+
+      // extract MeshBlock bounds
+      Real &x1min = size.h_view(ii0).x1min;
+      Real &x1max = size.h_view(ii0).x1max;
+      Real &x2min = size.h_view(ii0).x2min;
+      Real &x2max = size.h_view(ii0).x2max;
+      Real &x3min = size.h_view(ii0).x3min;
+      Real &x3max = size.h_view(ii0).x3max;
+
+      // set interpolation weights
+      for (int i=0; i<2*ng; ++i) {
+        iwghts.h_view(n,i,0) = 1.;
+        iwghts.h_view(n,i,1) = 1.;
+        iwghts.h_view(n,i,2) = 1.;
+        for (int j=0; j<2*ng; ++j) {
+          if (j != i) {
+            Real x1vpi1 = CellCenterX(ii1-ng+i+1, indcs.nx1, x1min, x1max);
+            Real x1vpj1 = CellCenterX(ii1-ng+j+1, indcs.nx1, x1min, x1max);
+            iwghts.h_view(n,i,0) *= (x0-x1vpj1)/(x1vpi1-x1vpj1);
+            Real x2vpi1 = CellCenterX(ii2-ng+i+1, indcs.nx2, x2min, x2max);
+            Real x2vpj1 = CellCenterX(ii2-ng+j+1, indcs.nx2, x2min, x2max);
+            iwghts.h_view(n,i,1) *= (y0-x2vpj1)/(x2vpi1-x2vpj1);
+            Real x3vpi1 = CellCenterX(ii3-ng+i+1, indcs.nx3, x3min, x3max);
+            Real x3vpj1 = CellCenterX(ii3-ng+j+1, indcs.nx3, x3min, x3max);
+            iwghts.h_view(n,i,2) *= (z0-x3vpj1)/(x3vpi1-x3vpj1);
+          }
+        }
+      }
+    }
+  }
+
+  // sync dual arrays
+  interp_wghts.template modify<HostMemSpace>();
+  interp_wghts.template sync<DevExeSpace>();
+
+  return;
+}
+//----------------------------------------------------------------------------------------
+//! \fn void GaussLegendreGrid::InterpolateToSphere
+//! \brief interpolate Cartesian data to surface of sphere
+
+void GaussLegendreGrid::InterpolateToSphere(int var_ind, DvceArray5D<Real> &val) {
+  // reinitialize interpolation indices and weights if AMR
+  if (pmy_pack->pmesh->adaptive) {
+    SetInterpolationIndices();
+    SetInterpolationWeights();
+  }
+  auto &indcs = pmy_pack->pmesh->mb_indcs;
+  int &is = indcs.is; int &js = indcs.js; int &ks = indcs.ks;
+  int &ng = indcs.ng;
+  int nang1 = nangles - 1;
+  int v = var_ind;
+
+  // reallocate container
+  Kokkos::realloc(interp_vals,nangles);
+  auto &iindcs = interp_indcs;
+  auto &iwghts = interp_wghts;
+  auto &ivals = interp_vals;
+  par_for("int2sph",DevExeSpace(),0,nang1,
+  KOKKOS_LAMBDA(int n) {
+    int &ii0 = iindcs.d_view(n,0);
+    int &ii1 = iindcs.d_view(n,1);
+    int &ii2 = iindcs.d_view(n,2);
+    int &ii3 = iindcs.d_view(n,3);
+
+    if (ii0==-1) {  // angle not on this rank
+      ivals.d_view(n) = 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(n,i,0)*iwghts.d_view(n,j,1)*iwghts.d_view(n,k,2);
+            int_value += iwght*val(ii0,v,ii3-(ng-k-ks)+1,ii2-(ng-j-js)+1,ii1-(ng-i-is)+1);
+          }
+        }
+      }
+      ivals.d_view(n) = int_value;
+    }
+  });
+
+  // sync dual arrays
+  interp_vals.template modify<DevExeSpace>();
+  interp_vals.template sync<HostMemSpace>();
+
+  return;
+}

src/geodesic-grid/gauss_legendre.hpp

@@ -0,0 +1,48 @@
+#ifndef GEODESIC_GRID_GAUSS_LEGENDRE_HPP_
+#define GEODESIC_GRID_GAUSS_LEGENDRE_HPP_
+//========================================================================================
+// AthenaXXX astrophysical plasma code
+// Copyright(C) 2020 James M. Stone <[email protected]> and the Athena code team
+// Licensed under the 3-clause BSD License (the "LICENSE")
+//========================================================================================
+//! \file geodesic_grid.hpp
+//  \brief definitions for GaussLegendreGrid class
+
+#include "athena.hpp"
+#include "athena_tensor.hpp"
+
+// Forward declarations
+class MeshBlockPack;
+
+//----------------------------------------------------------------------------------------
+//! \class GaussLegendreGrid
+
+class GaussLegendreGrid {
+ public:
+  GaussLegendreGrid(MeshBlockPack *pmy_pack, int ntheta, Real rad);
+  ~GaussLegendreGrid();
+    int nangles;  // total number of gridpoints
+    int ntheta;  // number of gridpoints along theta direction, nphi = 2ntheta
+    Real radius; // radius to initialize the sphere
+    DualArray1D<Real> int_weights;         // weights for quadrature integration
+    DualArray2D<Real> cart_pos;             // coord position (cartesian) at gridpoints
+    DualArray2D<Real> polar_pos;
+    DualArray1D<Real> interp_vals;   // container for data interpolated to sphere
+
+    // functions
+    void InitializeAngleAndWeights();
+    void InitializeRadius();
+
+    // interpolate scalar field to sphere
+    void InterpolateToSphere(int nvars, DvceArray5D<Real> &val);
+    DualArray2D<int> interp_indcs;   // indices of MeshBlock and zones therein for interp
+    DualArray3D<Real> interp_wghts;  // weights for interpolation
+
+    void SetInterpolationCoordinates();  // set indexing for interpolation
+    void SetInterpolationIndices();      // set indexing for interpolation
+    void SetInterpolationWeights();      // set weights for interpolation
+
+ private:
+    MeshBlockPack* pmy_pack;  // ptr to MeshBlockPack containing this Hydro
+};
+#endif // GEODESIC_GRID_GAUSS_LEGENDRE_HPP_

src/mesh/meshblock_pack.cpp

@@ -23,6 +23,7 @@
 #include "tasklist/numerical_relativity.hpp"
 #include "z4c/z4c.hpp"
 #include "dyn_grmhd/dyn_grmhd.hpp"
+#include "z4c/cce/cce.hpp"
 #include "diffusion/viscosity.hpp"
 #include "diffusion/resistivity.hpp"
 #include "radiation/radiation.hpp"
@@ -61,7 +62,14 @@ MeshBlockPack::~MeshBlockPack() {
   if (pnr    != nullptr) {delete pnr;}
   if (pturb  != nullptr) {delete pturb;}
   if (punit  != nullptr) {delete punit;}
-  if (pz4c   != nullptr) {delete pz4c;}
+  if (pz4c   != nullptr) {
+    delete pz4c;
+    // cce dump
+    for (auto cce : pz4c_cce) {
+      delete cce;
+    }
+    pz4c_cce.resize(0);
+  }
   if (ppart  != nullptr) {delete ppart;}
   // must be last, since it calls ~BoundaryValues() which (MPI) uses pmy_pack->pmb->nnghbr
   delete pmb;
@@ -186,6 +194,14 @@ void MeshBlockPack::AddPhysics(ParameterInput *pin) {
     pz4c = new z4c::Z4c(this, pin);
     padm = new adm::ADM(this, pin);
     ptmunu = nullptr;
+    // init cce dump
+    pz4c_cce.reserve(0);
+    int ncce = pin->GetOrAddInteger("cce", "num_radii", 0);
+    pz4c_cce.reserve(ncce);// 10 different components for each radius
+    for(int n = 0; n < ncce; ++n) {
+      // NOTE: these names are used for pittnull code, so DON'T change the convention
+      pz4c_cce.push_back(new z4c::CCE(pmesh, pin,n));
+    }
     nphysics++;
   } else {
     pz4c = nullptr;

src/mesh/meshblock_pack.hpp

@@ -31,6 +31,7 @@ namespace numrel {class NumericalRelativity;}
 class TurbulenceDriver;
 namespace radiation {class Radiation;}
 namespace z4c {class Z4c;}
+namespace z4c {class CCE;}
 namespace adm {class ADM;}
 namespace particles {class Particles;}
 namespace units {class Units;}
@@ -71,6 +72,7 @@ class MeshBlockPack {
   ion_neutral::IonNeutral *pionn=nullptr;
   TurbulenceDriver *pturb=nullptr;
   radiation::Radiation *prad=nullptr;
+  std::vector<z4c::CCE *> pz4c_cce;
   particles::Particles *ppart=nullptr;
 
   // units (needed to convert code units to cgs for, e.g., cooling or radiation)