diff --git a/CHANGELOG.md b/CHANGELOG.md index 3abf18c52538..00e37aa88f07 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,7 @@ ## Current develop ### Added (new features/APIs/variables/...) +- [[PR 1210]](https://github.com/parthenon-hpc-lab/parthenon/pull/1210) Add cycle based output - [[PR 1103]](https://github.com/parthenon-hpc-lab/parthenon/pull/1103) Add sparsity to vector wave equation test - [[PR 1185]](https://github.com/parthenon-hpc-lab/parthenon/pull/1185) Bugfix to particle defragmentation - [[PR 1184]](https://github.com/parthenon-hpc-lab/parthenon/pull/1184) Fix swarm block neighbor indexing in 1D, 2D diff --git a/doc/sphinx/src/outputs.rst b/doc/sphinx/src/outputs.rst index fa44e27da298..d8d861b540e4 100644 --- a/doc/sphinx/src/outputs.rst +++ b/doc/sphinx/src/outputs.rst @@ -6,10 +6,20 @@ Outputs Outputs from Parthenon are controlled via ```` blocks, where ``*`` should be replaced by a unique integer for each block. +The frequency of outputs can be controlled for each block separately +and can be triggered by either (simulation) time or cycle, i.e., + +- ``dt = 0.1`` means that the output for the block is written every 0.1 + in simulation time. +- ``dn = 100`` means that the output for the block is written every 100 + cycles. + +Note that only one option can be chosen for a given block. To disable an output block without removing it from the input file set -the block's ``dt < 0.0``. +the block's ``dt < 0.0`` and ``dn < 0`` (which is also happening by default +if the paramter is not provided in the input file). -In addition to time base outputs, two additional options to trigger +In addition to time or cycle based outputs, two additional options to trigger outputs (applies to HDF5, restart and histogram outputs) exist. - Signaling: If ``Parthenon`` catches a signal, e.g., ``SIGALRM`` which diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 37883ddcb7af..825014341aef 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -207,7 +207,6 @@ add_library(parthenon outputs/restart.hpp outputs/restart_hdf5.cpp outputs/restart_hdf5.hpp - outputs/vtk.cpp parthenon/driver.hpp parthenon/package.hpp diff --git a/src/outputs/ascent.cpp b/src/outputs/ascent.cpp index bab4632485f1..5392ef815b69 100644 --- a/src/outputs/ascent.cpp +++ b/src/outputs/ascent.cpp @@ -210,9 +210,15 @@ void AscentOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm, // advance output parameters output_params.file_number++; - output_params.next_time += output_params.dt; pin->SetInteger(output_params.block_name, "file_number", output_params.file_number); - pin->SetReal(output_params.block_name, "next_time", output_params.next_time); + if (output_params.dt > 0.0) { + output_params.next_time += output_params.dt; + pin->SetReal(output_params.block_name, "next_time", output_params.next_time); + } + if (output_params.dn > 0) { + output_params.next_n += output_params.dn; + pin->SetInteger(output_params.block_name, "next_n", output_params.next_n); + } } } // namespace parthenon diff --git a/src/outputs/histogram.cpp b/src/outputs/histogram.cpp index f983df0018e9..9d2da54831ca 100644 --- a/src/outputs/histogram.cpp +++ b/src/outputs/histogram.cpp @@ -559,9 +559,15 @@ void HistogramOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm // Only applies to default time-based data dumps, so that writing "now" and "final" // outputs does not change the desired output numbering. output_params.file_number++; - output_params.next_time += output_params.dt; pin->SetInteger(output_params.block_name, "file_number", output_params.file_number); - pin->SetReal(output_params.block_name, "next_time", output_params.next_time); + if (output_params.dt > 0.0) { + output_params.next_time += output_params.dt; + pin->SetReal(output_params.block_name, "next_time", output_params.next_time); + } + if (output_params.dn > 0) { + output_params.next_n += output_params.dn; + pin->SetInteger(output_params.block_name, "next_n", output_params.next_n); + } } } diff --git a/src/outputs/history.cpp b/src/outputs/history.cpp index 074bea4feb0c..89663d424074 100644 --- a/src/outputs/history.cpp +++ b/src/outputs/history.cpp @@ -211,9 +211,15 @@ void HistoryOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm, // advance output parameters output_params.file_number++; - output_params.next_time += output_params.dt; pin->SetInteger(output_params.block_name, "file_number", output_params.file_number); - pin->SetReal(output_params.block_name, "next_time", output_params.next_time); + if (output_params.dt > 0.0) { + output_params.next_time += output_params.dt; + pin->SetReal(output_params.block_name, "next_time", output_params.next_time); + } + if (output_params.dn > 0) { + output_params.next_n += output_params.dn; + pin->SetInteger(output_params.block_name, "next_n", output_params.next_n); + } } } // namespace parthenon diff --git a/src/outputs/output_parameters.hpp b/src/outputs/output_parameters.hpp index 8f97b0350e1f..dfc3abf10c40 100644 --- a/src/outputs/output_parameters.hpp +++ b/src/outputs/output_parameters.hpp @@ -44,6 +44,7 @@ struct OutputParameters { std::string data_format; std::vector packages; Real next_time, dt; + int next_n, dn; int file_number; bool include_ghost_zones, cartesian_vector; bool single_precision_output; @@ -53,7 +54,7 @@ struct OutputParameters { bool write_swarm_xdmf; // TODO(felker): some of the parameters in this class are not initialized in constructor OutputParameters() - : block_number(0), next_time(0.0), dt(-1.0), file_number(0), + : block_number(0), next_time(0.0), dt(-1.0), next_n(0), dn(-1), file_number(0), include_ghost_zones(false), cartesian_vector(false), single_precision_output(false), sparse_seed_nans(false), hdf5_compression_level(5), write_xdmf(false), write_swarm_xdmf(false) {} diff --git a/src/outputs/outputs.cpp b/src/outputs/outputs.cpp index 430ccdf3026c..feee23da478a 100644 --- a/src/outputs/outputs.cpp +++ b/src/outputs/outputs.cpp @@ -117,19 +117,31 @@ Outputs::Outputs(Mesh *pm, ParameterInput *pin, SimTime *tm) { op.block_name.assign(pib->block_name); Real dt = 0.0; // default value == 0 means that initial data is written by default - // for temporal drivers, setting dt to tlim ensures a final output is also written + int dn = -1; if (tm != nullptr) { - dt = pin->GetOrAddReal(op.block_name, "dt", tm->tlim); + dn = pin->GetOrAddInteger(op.block_name, "dn", -1.0); + + // If this is a dn controlled output (dn >= 0), soft disable dt based triggering + // (-> dt = -1), otherwise setting dt to tlim ensures a final output is also + // written for temporal drivers. + const auto tlim = dn >= 0 ? -1 : tm->tlim; + dt = pin->GetOrAddReal(op.block_name, "dt", tlim); } // if this output is "soft-disabled" (negative value) skip processing - if (dt < 0.0) { + if (dt < 0.0 && dn < 0) { pib = pib->pnext; // move to next input block name continue; } + + PARTHENON_REQUIRE_THROWS(!(dt >= 0.0 && dn >= 0), + "dt and dn are enabled for the same output block, which " + "is not supported. Please set at most one value >= 0."); // set time of last output, time between outputs if (tm != nullptr) { op.next_time = pin->GetOrAddReal(op.block_name, "next_time", tm->time); op.dt = dt; + op.next_n = pin->GetOrAddInteger(op.block_name, "next_n", tm->ncycle); + op.dn = dn; } // set file number, basename, id, and format @@ -259,8 +271,6 @@ Outputs::Outputs(Mesh *pm, ParameterInput *pin, SimTime *tm) { if (op.file_type == "hst") { pnew_type = new HistoryOutput(op); num_hst_outputs++; - } else if (op.file_type == "vtk") { - pnew_type = new VTKOutput(op); } else if (op.file_type == "ascent") { pnew_type = new AscentOutput(op); } else if (op.file_type == "histogram") { @@ -438,9 +448,19 @@ void Outputs::MakeOutputs(Mesh *pm, ParameterInput *pin, SimTime *tm, OutputType *ptype = pfirst_type_; while (ptype != nullptr) { if ((tm == nullptr) || - ((ptype->output_params.dt >= 0.0) && - ((tm->ncycle == 0) || (tm->time >= ptype->output_params.next_time) || - (tm->time >= tm->tlim) || (signal == SignalHandler::OutputSignal::now) || + // output is not soft disabled and + (((ptype->output_params.dt >= 0.0) || (ptype->output_params.dn >= 0)) && + // either dump initial data + ((tm->ncycle == 0) || + // or by triggering time or cycle based conditions + ((ptype->output_params.dt >= 0.0) && + ((tm->time >= ptype->output_params.next_time) || + (tm->tlim > 0.0 && tm->time >= tm->tlim))) || + ((ptype->output_params.dn >= 0) && + ((tm->ncycle >= ptype->output_params.next_n) || + (tm->nlim > 0 && tm->ncycle >= tm->nlim))) || + // or by manual triggers + (signal == SignalHandler::OutputSignal::now) || (signal == SignalHandler::OutputSignal::final) || (signal == SignalHandler::OutputSignal::analysis && ptype->output_params.analysis_flag)))) { diff --git a/src/outputs/outputs.hpp b/src/outputs/outputs.hpp index 9d37245eabc8..c2a351654063 100644 --- a/src/outputs/outputs.hpp +++ b/src/outputs/outputs.hpp @@ -143,18 +143,6 @@ class HistoryOutput : public OutputType { const SignalHandler::OutputSignal signal) override; }; -//---------------------------------------------------------------------------------------- -//! \class VTKOutput -// \brief derived OutputType class for vtk dumps - -class VTKOutput : public OutputType { - public: - explicit VTKOutput(const OutputParameters &oparams) : OutputType(oparams) {} - void WriteContainer(SimTime &tm, Mesh *pm, ParameterInput *pin, bool flag) override; - void WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm, - const SignalHandler::OutputSignal signal) override; -}; - //---------------------------------------------------------------------------------------- //! \class AscentOutput // \brief derived OutputType class for Ascent in situ situ visualization and analysis diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index cd5009d490d9..3068c07a1048 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -597,9 +597,15 @@ std::string PHDF5Output::GenerateFilename_(ParameterInput *pin, SimTime *tm, // Only applies to default time-based data dumps, so that writing "now" and "final" // outputs does not change the desired output numbering. output_params.file_number++; - output_params.next_time += output_params.dt; pin->SetInteger(output_params.block_name, "file_number", output_params.file_number); - pin->SetReal(output_params.block_name, "next_time", output_params.next_time); + if (output_params.dt > 0.0) { + output_params.next_time += output_params.dt; + pin->SetReal(output_params.block_name, "next_time", output_params.next_time); + } + if (output_params.dn > 0) { + output_params.next_n += output_params.dn; + pin->SetInteger(output_params.block_name, "next_n", output_params.next_n); + } } return filename; } diff --git a/src/outputs/vtk.cpp b/src/outputs/vtk.cpp deleted file mode 100644 index 668b41c3b105..000000000000 --- a/src/outputs/vtk.cpp +++ /dev/null @@ -1,240 +0,0 @@ -//======================================================================================== -// Athena++ astrophysical MHD code -// Copyright(C) 2014 James M. Stone and other code contributors -// Licensed under the 3-clause BSD License, see LICENSE file for details -//======================================================================================== -// (C) (or copyright) 2020-2021. Triad National Security, LLC. All rights reserved. -// -// This program was produced under U.S. Government contract 89233218CNA000001 for Los -// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC -// for the U.S. Department of Energy/National Nuclear Security Administration. All rights -// in the program are reserved by Triad National Security, LLC, and the U.S. Department -// of Energy/National Nuclear Security Administration. The Government is granted for -// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide -// license in this material to reproduce, prepare derivative works, distribute copies to -// the public, perform publicly and display publicly, and to permit others to do so. -//======================================================================================== -//! \file vtk.cpp -// \brief writes output data in (legacy) vtk format. -// Data is written in RECTILINEAR_GRID geometry, in BINARY format, and in FLOAT type -// Writes one file per MeshBlock. - -#include -#include -#include -#include -#include -#include -#include -#include - -#include "coordinates/coordinates.hpp" -#include "defs.hpp" -#include "mesh/mesh.hpp" -#include "mesh/meshblock.hpp" -#include "outputs/outputs.hpp" -#include "parthenon_arrays.hpp" -#include "utils/error_checking.hpp" - -namespace parthenon { - -//---------------------------------------------------------------------------------------- -// Functions to detect big endian machine, and to byte-swap 32-bit words. The vtk -// legacy format requires data to be stored as big-endian. - -int IsBigEndian() { - std::int32_t n = 1; - // careful! although int -> char * -> int round-trip conversion is safe, - // an arbitrary char* may not be converted to int* - char *ep = reinterpret_cast(&n); - return (*ep == 0); // Returns 1 (true) on a big endian machine -} - -// namespace { -// inline void Swap4Bytes(void *vdat) { -// char tmp, *dat = static_cast(vdat); -// tmp = dat[0]; -// dat[0] = dat[3]; -// dat[3] = tmp; -// tmp = dat[1]; -// dat[1] = dat[2]; -// dat[2] = tmp; -// } -// } // namespace - -//---------------------------------------------------------------------------------------- -//! \fn void VTKOutput:::WriteOutputFile(Mesh *pm) -// \brief Cycles over all MeshBlocks and writes OutputData in (legacy) vtk format, one -// MeshBlock per file - -void VTKOutput::WriteContainer(SimTime &tm, Mesh *pm, ParameterInput *pin, bool flag) { - throw std::runtime_error(std::string(__func__) + " is not implemented"); - /* - int big_end = IsBigEndian(); // =1 on big endian machine - - // Loop over MeshBlocks - for (auto &pmb : pm->block_list) { - // set start/end array indices depending on whether ghost zones are included - IndexDomain domain = IndexDomain::interior; - if (output_params.include_ghost_zones) { - domain = IndexDomain::entire; - } - // build doubly linked list of OutputData nodes (setting data ptrs to appropriate - // quantity on MeshBlock for each node), then slice/sum as needed - // create filename: "file_basename"+ "."+"blockid"+"."+"file_id"+"."+XXXXX+".vtk", - // where XXXXX = 5-digit file_number - std::string fname; - char number[6]; - std::snprintf(number, sizeof(number), "%05d", output_params.file_number); - char blockid[12]; - std::snprintf(blockid, sizeof(blockid), "block%d", pmb->gid); - - fname.assign(output_params.file_basename); - fname.append(".N."); - fname.append(blockid); - fname.append("."); - fname.append(output_params.file_id); - fname.append("."); - fname.append(number); - fname.append(".vtk"); - - // open file for output - FILE *pfile; - std::stringstream msg; - if ((pfile = std::fopen(fname.c_str(), "w")) == nullptr) { - msg << "### FATAL ERROR in function [VTKOutput::WriteOutputFile]" << std::endl - << "Output file '" << fname << "' could not be opened" << std::endl; - PARTHENON_FAIL(msg); - } - - // There are five basic parts to the VTK "legacy" file format. - // 1. Write file version and identifier - std::fprintf(pfile, "# vtk DataFile Version 2.0\n"); - - // 2. Header - std::fprintf(pfile, "# Athena++ data at time=%e", tm.time); - std::fprintf(pfile, " cycle=%d", tm.ncycle); - std::fprintf(pfile, " variables=%s \n", output_params.variable.c_str()); - - // 3. File format - std::fprintf(pfile, "BINARY\n"); - - // 4. Dataset structure - int ncells1 = pmb->cellbounds.ncellsi(domain); - int ncells2 = pmb->cellbounds.ncellsj(domain); - int ncells3 = pmb->cellbounds.ncellsk(domain); - int ncoord1 = ncells1; - if (ncells1 > 1) ncoord1++; - int ncoord2 = ncells2; - if (ncells2 > 1) ncoord2++; - int ncoord3 = ncells3; - if (ncells3 > 1) ncoord3++; - - float *data; - int ndata = std::max(ncoord1, ncoord2); - ndata = std::max(ndata, ncoord3); - data = new float[3 * ndata + 1]; - - // Specify the type of data, dimensions, and coordinates. If N>1, then write N+1 - // cell faces as binary floats. If N=1, then write 1 cell center position. - std::fprintf(pfile, "DATASET RECTILINEAR_GRID\n"); - std::fprintf(pfile, "DIMENSIONS %d %d %d\n", ncoord1, ncoord2, ncoord3); - - IndexRange out_ib = pmb->cellbounds.GetBoundsI(domain); - IndexRange out_jb = pmb->cellbounds.GetBoundsJ(domain); - IndexRange out_kb = pmb->cellbounds.GetBoundsK(domain); - // write x1-coordinates as binary float in big endian order - std::fprintf(pfile, "X_COORDINATES %d float\n", ncoord1); - if (ncells1 == 1) { - data[0] = static_cast(pmb->coords.Xc<1>(out_ib.s)); - } else { - for (int i = out_ib.s; i <= out_ib.e + 1; ++i) { - data[i - out_ib.s] = static_cast(pmb->coords.Xf<1>(i)); - } - } - if (!big_end) { - for (int i = 0; i < ncoord1; ++i) - Swap4Bytes(&data[i]); - } - std::fwrite(data, sizeof(float), static_cast(ncoord1), pfile); - - // write x2-coordinates as binary float in big endian order - std::fprintf(pfile, "\nY_COORDINATES %d float\n", ncoord2); - if (ncells2 == 1) { - data[0] = static_cast(pmb->coords.Xc<2>(out_jb.s)); - } else { - for (int j = out_jb.s; j <= out_jb.e + 1; ++j) { - data[j - out_jb.s] = static_cast(pmb->coords.Xf<2>(j)); - } - } - if (!big_end) { - for (int i = 0; i < ncoord2; ++i) - Swap4Bytes(&data[i]); - } - std::fwrite(data, sizeof(float), static_cast(ncoord2), pfile); - - // write x3-coordinates as binary float in big endian order - std::fprintf(pfile, "\nZ_COORDINATES %d float\n", ncoord3); - if (ncells3 == 1) { - data[0] = static_cast(pmb->coords.Xc<3>(out_kb.s)); - } else { - for (int k = out_kb.s; k <= out_kb.e + 1; ++k) { - data[k - out_kb.s] = static_cast(pmb->coords.Xf<3>(k)); - } - } - if (!big_end) { - for (int i = 0; i < ncoord3; ++i) - Swap4Bytes(&data[i]); - } - std::fwrite(data, sizeof(float), static_cast(ncoord3), pfile); - - // 5. Data. An arbitrary number of scalars and vectors can be written (every node - // in the OutputData doubly linked lists), all in binary floats format - - std::fprintf(pfile, "\nCELL_DATA %d", pmb->cellbounds.GetTotal(domain)); - // reset container iterator to point to current block data - auto ci = - MeshBlockDataIterator(pmb->meshblock_data.Get(), {Metadata::Graphics}); - for (auto &v : ci.vars) { - if (!data) { - std::cout << "____________________SKIPPPING:" << v->label() << std::endl; - continue; - } - std::fprintf(pfile, "\nLOOKUP_TABLE default\n"); - for (int k = out_kb.s; k <= out_kb.e; k++) { - for (int j = out_jb.s; j <= out_jb.e; j++) { - int index = 0; - for (int i = out_ib.s; i <= out_ib.e; i++, index++) { - data[(i - out_ib.s) + index] = (*v)(k, j, i); - } - - // write data in big endian order - if (!big_end) { - for (int i = 0; i < (ncells1); ++i) - Swap4Bytes(&data[i]); - } - std::fwrite(data, sizeof(float), static_cast(ncells1), pfile); - } - } - } - - // don't forget to close the output file and clean up ptrs to data in OutputData - std::fclose(pfile); - delete[] data; - } // end loop over MeshBlocks - - // increment counters - output_params.file_number++; - output_params.next_time += output_params.dt; - pin->SetInteger(output_params.block_name, "file_number", output_params.file_number); - pin->SetReal(output_params.block_name, "next_time", output_params.next_time); - - return; - */ -} -void VTKOutput::WriteOutputFile(Mesh *pm, ParameterInput *pin, SimTime *tm, - const SignalHandler::OutputSignal signal) { - throw std::runtime_error(std::string(__func__) + " is not implemented"); -} - -} // namespace parthenon