From c3f8919a764fec71cdb8f1e612d07a217e0010a2 Mon Sep 17 00:00:00 2001 From: David Radice Date: Sat, 25 Jan 2025 17:58:55 +0000 Subject: [PATCH 1/6] Fix interpolation on Chebyshev grid Fix bug when interpolating on Chebyshev grid in which the indices were computed with the Chebyshev grids, but the interpolation weights still used the Cartesian grid --- src/utils/cart_grid.cpp | 7 +++++-- src/utils/cart_grid.hpp | 3 +-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/utils/cart_grid.cpp b/src/utils/cart_grid.cpp index cb89cdca..ddb79bf6 100644 --- a/src/utils/cart_grid.cpp +++ b/src/utils/cart_grid.cpp @@ -74,8 +74,6 @@ CartesianGrid::CartesianGrid(MeshBlockPack *pmy_pack, Real center[3], return; } -CartesianGrid::~CartesianGrid() {} - void CartesianGrid::ResetCenter(Real center[3]) { // grid center center_x1 = center[0]; @@ -184,6 +182,11 @@ void CartesianGrid::SetInterpolationWeights() { Real x0 = min_x1 + nx * d_x1; Real y0 = min_x2 + ny * d_x2; Real z0 = min_x3 + nz * d_x3; + if (is_cheby) { + x0 = center_x1 + extend_x1*std::cos(nx*M_PI/(nx1-1)); + y0 = center_x2 + extend_x2*std::cos(ny*M_PI/(nx2-1)); + z0 = center_x3 + extend_x3*std::cos(nz*M_PI/(nx3-1)); + } // extract MeshBlock bounds Real &x1min = size.h_view(ii0).x1min; Real &x1max = size.h_view(ii0).x1max; diff --git a/src/utils/cart_grid.hpp b/src/utils/cart_grid.hpp index 01c53d23..ba70a24e 100644 --- a/src/utils/cart_grid.hpp +++ b/src/utils/cart_grid.hpp @@ -21,8 +21,7 @@ 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(); + Real extend[3], int numpoints[3], bool is_cheb = false); // parameters for the grid Real center_x1, center_x2, center_x3; // grid centers From c08e3ea182a049fb8728ed2171d0a61278a646f7 Mon Sep 17 00:00:00 2001 From: David Radice Date: Sat, 25 Jan 2025 18:02:10 +0000 Subject: [PATCH 2/6] Added new CartesianGridOuput output class `CartesianGridOutput` can output any variable in a Cartesian box (within the simulation domain) at a prescribed position and extent. To use it, add a block to the output file like this: ``` file_type = cart variable = hydro_w dt = 0.1 center_x = 0.0 center_y = 0.0 center_z = 0.0 extent_x = 0.1 extent_y = 0.1 extent_z = 0.2 chebyshev = 1 # 0:uniform grid; 1:Chebyshev-Lobatto grid ``` Note that the extent of the grid is actually half of the extent variable. This is because we are using the same convention as in `CartesianGrid`. This commit also adds `vis/python/cartgrid.py` to read the data produced by `CartesianGridOutput` --- src/CMakeLists.txt | 1 + src/outputs/cartgrid.cpp | 140 +++++++++++++++++++++++++++++++++++++++ src/outputs/outputs.cpp | 3 + src/outputs/outputs.hpp | 28 ++++++++ vis/python/cartgrid.py | 70 ++++++++++++++++++++ 5 files changed, 242 insertions(+) create mode 100644 src/outputs/cartgrid.cpp create mode 100644 vis/python/cartgrid.py diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c442f864..ead8caf4 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -88,6 +88,7 @@ add_executable( outputs/io_wrapper.cpp outputs/outputs.cpp outputs/basetype_output.cpp + outputs/cartgrid.cpp outputs/derived_variables.cpp outputs/binary.cpp outputs/eventlog.cpp diff --git a/src/outputs/cartgrid.cpp b/src/outputs/cartgrid.cpp new file mode 100644 index 00000000..7c6e4934 --- /dev/null +++ b/src/outputs/cartgrid.cpp @@ -0,0 +1,140 @@ +//======================================================================================== +// 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 cartgrid.cpp +//! \brief writes data on a Cartesian sub-grid in binary format + +#include // mkdir + +#include // snprintf +#include +#include +#include + +#include "athena.hpp" +#include "globals.hpp" +#include "mesh/mesh.hpp" +#include "outputs.hpp" +#include "parameter_input.hpp" +#include "utils/cart_grid.hpp" + +CartesianGridOutput::CartesianGridOutput(ParameterInput *pin, Mesh *pm, + OutputParameters op) + : BaseTypeOutput(pin, pm, op) { + Real center[3], extent[3]; + int numpoints[3]; + bool is_cheb; + + mkdir("cart", 0755); + + center[0] = pin->GetOrAddReal(op.block_name, "center_x", 0.0); + center[1] = pin->GetOrAddReal(op.block_name, "center_y", 0.0); + center[2] = pin->GetOrAddReal(op.block_name, "center_z", 0.0); + extent[0] = pin->GetOrAddReal(op.block_name, "extent_x", 2.0); + extent[1] = pin->GetOrAddReal(op.block_name, "extent_y", 2.0); + extent[2] = pin->GetOrAddReal(op.block_name, "extent_z", 2.0); + numpoints[0] = pin->GetOrAddInteger(op.block_name, "numpoints_x", 32); + numpoints[1] = pin->GetOrAddInteger(op.block_name, "numpoints_y", 32); + numpoints[2] = pin->GetOrAddInteger(op.block_name, "numpoints_z", 32); + is_cheb = pin->GetOrAddBoolean(op.block_name, "chebyshev", false); + + pcart = new CartesianGrid(pm->pmb_pack, center, extent, numpoints, is_cheb); + + for (int d = 0; d < 3; ++d) { + md.center[d] = center[d]; + md.extent[d] = extent[d]; + md.numpoints[d] = numpoints[d]; + } + md.is_cheb = is_cheb; +} + +CartesianGridOutput::~CartesianGridOutput() { delete pcart; } + +void CartesianGridOutput::LoadOutputData(Mesh *pm) { + int nout_vars = outvars.size(); + Kokkos::realloc(outarray, nout_vars, 1, md.numpoints[0], md.numpoints[1], + md.numpoints[2]); + for (int n = 0; n < nout_vars; ++n) { + pcart->InterpolateToGrid(outvars[n].data_index, *(outvars[n].data_ptr)); + auto v_slice = + Kokkos::subview(outarray, n, 0, Kokkos::ALL, Kokkos::ALL, Kokkos::ALL); + Kokkos::deep_copy(v_slice, pcart->interp_vals.h_view); + } +#if MPI_PARALLEL_ENABLED + // Note that InterpolateToGrid will set zero values for points not owned by + // current rank + int count = nout_vars * md.numpoints[0] * md.numpoints[1] * md.numpoints[2]; + if (0 == global_variable::my_rank) { + MPI_Reduce(MPI_IN_PLACE, outarray.data(), count, MPI_ATHENA_REAL, MPI_SUM, + 0, MPI_COMM_WORLD); + } else { + MPI_Reduce(outarray.data(), outarray.data(), count, MPI_ATHENA_REAL, + MPI_SUM, 0, MPI_COMM_WORLD); + } +#endif +} + +void CartesianGridOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin) { +#if MPI_PARALLEL_ENABLED + if (0 == global_variable::my_rank) { +#endif + + // Assemble filename + char fname[BUFSIZ]; + std::snprintf(fname, BUFSIZ, "cart/%s.%s.%05d.bin", + out_params.file_basename.c_str(), out_params.file_id.c_str(), + out_params.file_number); + + // Open file + std::ofstream ofile(fname, std::ios::binary); + + // Write metadata + md.cycle = pm->ncycle; + md.time = pm->time; + md.noutvars = outvars.size(); + ofile.write(reinterpret_cast(&md), sizeof(MetaData)); + + // Write list of variables + { + std::stringstream msg; + for (int n = 0; n < md.noutvars - 1; ++n) { + msg << outvars[n].label << " "; + } + msg << outvars[md.noutvars - 1].label; + std::string smsg = msg.str(); + int len = smsg.size(); + ofile.write(reinterpret_cast(&len), sizeof(int)); + ofile.write(smsg.c_str(), len); + } + + // Write actual data + for (int n = 0; n < md.noutvars; ++n) { + for (int k = 0; k < md.numpoints[2]; ++k) { + for (int j = 0; j < md.numpoints[1]; ++j) { + for (int i = 0; i < md.numpoints[0]; ++i) { + // Note that we are accessing the array with the convention of + // CartesianGrid which is opposite from the one used in the rest of + // the code, but we write the output as k, j, i + float var = outarray(n, 0, i, j, k); + ofile.write(reinterpret_cast(&var), sizeof(float)); + } + } + } + } + +#if MPI_PARALLEL_ENABLED + } +#endif + + // increment counters + out_params.file_number++; + if (out_params.last_time < 0.0) { + out_params.last_time = pm->time; + } else { + out_params.last_time += out_params.dt; + } + pin->SetInteger(out_params.block_name, "file_number", out_params.file_number); + pin->SetReal(out_params.block_name, "last_time", out_params.last_time); +} diff --git a/src/outputs/outputs.cpp b/src/outputs/outputs.cpp index e047f13b..36dd7b99 100644 --- a/src/outputs/outputs.cpp +++ b/src/outputs/outputs.cpp @@ -272,6 +272,9 @@ Outputs::Outputs(ParameterInput *pin, Mesh *pm) { } else if (opar.file_type.compare("bin") == 0) { pnode = new MeshBinaryOutput(pin,pm,opar); pout_list.insert(pout_list.begin(),pnode); + } else if (opar.file_type.compare("cart") == 0) { + pnode = new CartesianGridOutput(pin,pm,opar); + pout_list.insert(pout_list.begin(),pnode); } else if (opar.file_type.compare("rst") == 0) { // Add restarts to the tail end of BaseTypeOutput list, so file counters for other // output types are up-to-date in restart file diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index a8ebf4ac..c9784492 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -384,6 +384,34 @@ class RestartOutput : public BaseTypeOutput { void WriteOutputFile(Mesh *pm, ParameterInput *pin) override; }; +// Forward declaration +class CartesianGrid; + +//---------------------------------------------------------------------------------------- +//! \class CartesianGridOutput +// \brief derived BaseTypeOutput class for output on a Cartesian grid +class CartesianGridOutput : public BaseTypeOutput { + struct MetaData { + int cycle; + float time; + float center[3]; + float extent[3]; + int numpoints[3]; + bool is_cheb; + int noutvars; + }; + public: + CartesianGridOutput(ParameterInput *pin, Mesh *pm, OutputParameters oparams); + ~CartesianGridOutput(); + //! Interpolate the data on the Cartesian grid and handle MPI communication + void LoadOutputData(Mesh *pm) override; + //! Write the data to file + void WriteOutputFile(Mesh *pm, ParameterInput *pin) override; + private: + CartesianGrid *pcart; + MetaData md; +}; + //---------------------------------------------------------------------------------------- //! \class EventLogOutput // \brief derived BaseTypeOutput class for event counter data diff --git a/vis/python/cartgrid.py b/vis/python/cartgrid.py new file mode 100644 index 00000000..56be6ad7 --- /dev/null +++ b/vis/python/cartgrid.py @@ -0,0 +1,70 @@ +""" +Reads CartesianGrid binary format +""" + +import numpy as np +from struct import Struct + + +class CartesianGridData: + """Class representing a single CartesianGrid dump + + Members are: + * cycle:int simulation cycle + * time:float simulation time + * center[3]:float box center + * extent[3]:float box extent + * numpoints[3]:int number of grid points in each direction + * is_cheb:bool Chebyshev grid + * variables dictionary with all grid functions + """ + + def __init__(self, fname): + mdata = Struct("i7f3i?2i") + with open(fname, "rb") as f: + # Parse metadata + blob = f.read(mdata.size) + ncycle, time, cx, cy, cz, ex, ey, ez, nx, ny, nz, cheb, nout, nstr = ( + mdata.unpack(blob) + ) + self.cycle = ncycle + self.time = time + self.center = (cx, cy, cz) + self.extent = (ex, ey, ez) + self.numpoints = (nx, ny, nz) + self.is_cheb = cheb + self.variables = {} + + # Read variable names + blob = f.read(nstr).decode("ascii") + names = blob.split() + if len(names) != nout: + raise IOError(f"Could not read {fname}") + + # Read variables + for n in names: + self.variables[n] = ( + np.fromfile(f, dtype=np.float32, count=np.prod(self.numpoints)) + .reshape(self.numpoints) + .transpose() + ) + + def coords(self, d=None): + """Returns the coordinates""" + if d is None: + return self.coords(0), self.coords(1), self.coords(2) + if self.is_cheb: + from math import pi + + return self.center[d] + self.extent[d] * np.cos( + np.linspace(0, pi, self.numpoints[d]) + ) + else: + return self.center[d] + self.extent[d] * np.linspace( + -1, 1, self.numpoints[d] + ) + + def meshgrid(self): + """Returns a mesh grid with all the coordinates""" + x, y, z = self.coords() + return np.meshgrid(x, y, z, indexing="ij") From 8d49d194d20b7795d997147454897012933f8ce1 Mon Sep 17 00:00:00 2001 From: David Radice Date: Fri, 31 Jan 2025 15:16:33 +0000 Subject: [PATCH 3/6] Add support for derived variables in CartesianGridOutput --- src/outputs/cartgrid.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/outputs/cartgrid.cpp b/src/outputs/cartgrid.cpp index 7c6e4934..583ffa2d 100644 --- a/src/outputs/cartgrid.cpp +++ b/src/outputs/cartgrid.cpp @@ -56,6 +56,12 @@ void CartesianGridOutput::LoadOutputData(Mesh *pm) { int nout_vars = outvars.size(); Kokkos::realloc(outarray, nout_vars, 1, md.numpoints[0], md.numpoints[1], md.numpoints[2]); + + // Calculate derived variables, if required + if (out_params.contains_derived) { + ComputeDerivedVariable(out_params.variable, pm); + } + for (int n = 0; n < nout_vars; ++n) { pcart->InterpolateToGrid(outvars[n].data_index, *(outvars[n].data_ptr)); auto v_slice = From eea9a06185976b889173a0b2b135ff00fff517cb Mon Sep 17 00:00:00 2001 From: David Radice Date: Fri, 31 Jan 2025 16:11:29 +0000 Subject: [PATCH 4/6] Fix CartesianGridOutput with AMR Reset CartesianGrid interpolation indices and weights before interpolation if AMR is enabled --- src/outputs/cartgrid.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/outputs/cartgrid.cpp b/src/outputs/cartgrid.cpp index 583ffa2d..eb66cb5d 100644 --- a/src/outputs/cartgrid.cpp +++ b/src/outputs/cartgrid.cpp @@ -53,6 +53,12 @@ CartesianGridOutput::CartesianGridOutput(ParameterInput *pin, Mesh *pm, CartesianGridOutput::~CartesianGridOutput() { delete pcart; } void CartesianGridOutput::LoadOutputData(Mesh *pm) { + // If AMR is enabled we need to reset the CartesianGrid + if (pm->adaptive) { + pcart->SetInterpolationIndices(); + pcart->SetInterpolationWeights(); + } + int nout_vars = outvars.size(); Kokkos::realloc(outarray, nout_vars, 1, md.numpoints[0], md.numpoints[1], md.numpoints[2]); From 2a0853b1cf82002ed3e8b443847d82ec88d58cf2 Mon Sep 17 00:00:00 2001 From: David Radice Date: Fri, 31 Jan 2025 18:41:39 +0000 Subject: [PATCH 5/6] Fix bug introduced in previous commit --- src/utils/cart_grid.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/utils/cart_grid.hpp b/src/utils/cart_grid.hpp index ba70a24e..12fa4470 100644 --- a/src/utils/cart_grid.hpp +++ b/src/utils/cart_grid.hpp @@ -38,13 +38,13 @@ class CartesianGrid { 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 + void SetInterpolationIndices(); // set indexing for interpolation + void SetInterpolationWeights(); // set weights 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_ From 9a3a5149fae8a7e8fc072f444137c581e64970f0 Mon Sep 17 00:00:00 2001 From: David Radice Date: Sun, 2 Feb 2025 14:05:02 +0000 Subject: [PATCH 6/6] Fix typo introduced in merge --- src/utils/cart_grid.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/utils/cart_grid.cpp b/src/utils/cart_grid.cpp index ab7c5b75..6f87dd39 100644 --- a/src/utils/cart_grid.cpp +++ b/src/utils/cart_grid.cpp @@ -209,9 +209,9 @@ void CartesianGrid::SetInterpolationWeights() { Real y0 = min_x2 + ny * d_x2; Real z0 = min_x3 + nz * d_x3; if (is_cheby) { - x0 = center_x1 + extend_x1*std::cos(nx*M_PI/(nx1-1)); - y0 = center_x2 + extend_x2*std::cos(ny*M_PI/(nx2-1)); - z0 = center_x3 + extend_x3*std::cos(nz*M_PI/(nx3-1)); + x0 = center_x1 + extent_x1*std::cos(nx*M_PI/(nx1-1)); + y0 = center_x2 + extent_x2*std::cos(ny*M_PI/(nx2-1)); + z0 = center_x3 + extent_x3*std::cos(nz*M_PI/(nx3-1)); } // extract MeshBlock bounds Real &x1min = size.h_view(ii0).x1min;