diff --git a/include/bout/petsc_preconditioner.hxx b/include/bout/petsc_preconditioner.hxx index be283413fc..9fe33f7330 100644 --- a/include/bout/petsc_preconditioner.hxx +++ b/include/bout/petsc_preconditioner.hxx @@ -10,8 +10,11 @@ #ifndef BOUT_PETSC_PRECONDITIONER_H #define BOUT_PETSC_PRECONDITIONER_H +#include "bout/bout_enum_class.hxx" #include "bout/build_defines.hxx" +BOUT_ENUM_CLASS(PetscMatrixExportFormat, binary, ascii); + #if BOUT_HAS_PETSC #include "bout/petsc_interface.hxx" @@ -22,6 +25,8 @@ #include #include +#include + class Options; class Field3D; @@ -73,6 +78,13 @@ public: Mat jacobian() const { return Jfd; } MatFDColoring coloring() const { return fdcoloring; } + static PetscErrorCode + saveMatrix(Mat matrix, const std::string& filename, + PetscMatrixExportFormat format = PetscMatrixExportFormat::binary); + PetscErrorCode + saveMatrix(const std::string& filename, + PetscMatrixExportFormat format = PetscMatrixExportFormat::binary) const; + void reset(); private: @@ -86,6 +98,9 @@ private: // unconditionally in PETSc-enabled compilation units. class PetscPreconditioner { public: + void saveMatrix( + const std::string& UNUSED(filename), + PetscMatrixExportFormat UNUSED(format) = PetscMatrixExportFormat::binary) const {} void reset() {} }; diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index 0d4b330ec6..10a3d294f1 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -38,6 +38,7 @@ #include "bout/build_defines.hxx" +#include "bout/bout_enum_class.hxx" #include "bout/bout_types.hxx" #include "bout/boutexception.hxx" #include "bout/globals.hxx" @@ -101,6 +102,8 @@ constexpr auto SOLVERRKGENERIC = "rkgeneric"; enum class FieldCategories : std::uint8_t { VARS, DERIVS, MMS }; enum class SOLVER_VAR_OP : std::uint8_t { LOAD, SET_ID, SAVE }; +BOUT_ENUM_CLASS(JacobianExportKind, system, scaled, rhs); + /// A type to set where in the list monitors are added enum class MonitorPosition { BACK, FRONT }; @@ -368,6 +371,25 @@ public: protected: friend class SundialsNVectorInterface; + struct JacobianVariableMetadata { + int offset{0}; + std::string name; + std::string location; + bool evolve_bndry{false}; + bool constraint{false}; + std::string description; + }; + + struct JacobianMetadata { + int format_version{1}; + std::string solver_name; + int n2d{0}; + int n3d{0}; + std::vector variables_2d; + std::vector variables_3d; + std::string ordering; + }; + /// Number of command-line arguments static int* pargc; /// Command-line arguments @@ -606,6 +628,12 @@ protected: /// Returns a Field3D containing the global indices Field3D globalIndex(int localStart); + Field3D jacobianIndexBase(int localStart = 0); + std::vector getJacobianMetadata2D() const; + std::vector getJacobianMetadata3D() const; + JacobianMetadata getJacobianMetadata(const std::string& solver_name) const; + void writeJacobianMetadataJson(const std::string& filename, + const std::string& solver_name) const; /// Maximum internal timestep BoutReal max_dt{-1.0}; @@ -670,6 +698,7 @@ private: std::string run_restart_from = "yyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyy"; /// Save `run_id` and `run_restart_from` every output bool save_repeat_run_id{false}; + bool save_jacobian_index_base{false}; /// Current iteration (output time-step) number int iteration{0}; diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index e7d367331a..8c40cf5722 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -21,6 +21,7 @@ #include #include #include +#include #include #include "petscerror.h" @@ -61,6 +62,15 @@ PetscErrorCode FormFunctionForColoring(void* UNUSED(snes), Vec x, Vec f, void* c return static_cast(ctx)->snes_function(x, f, true); } +PetscErrorCode FormRawFunctionForColoring(void* UNUSED(snes), Vec x, Vec f, void* ctx) { + return static_cast(ctx)->raw_rhs_function(x, f, true); +} + +PetscErrorCode FormScaledFunctionForColoring(void* UNUSED(snes), Vec x, Vec f, + void* ctx) { + return static_cast(ctx)->scaled_rhs_function(x, f, true); +} + PetscErrorCode snesPCapply(PC pc, Vec x, Vec y) { // Get the context SNESSolver* s; @@ -71,6 +81,8 @@ PetscErrorCode snesPCapply(PC pc, Vec x, Vec y) { PetscErrorCode ComputeJacobianScaledColor(SNES snes, Vec x1, Mat Jac, Mat Jac_new, void* ctx); +PetscErrorCode ComputeJacobianDefaultMaybeExport(SNES snes, Vec x1, Mat Jac, Mat Jac_new, + void* ctx); } // namespace PetscErrorCode SNESSolver::FDJinitialise() { @@ -111,9 +123,9 @@ PetscErrorCode SNESSolver::FDJinitialise() { nullptr, &Jfd); if (matrix_free_operator) { - SNESSetJacobian(snes, Jmf, Jfd, SNESComputeJacobianDefault, this); + SNESSetJacobian(snes, Jmf, Jfd, ComputeJacobianDefaultMaybeExport, this); } else { - SNESSetJacobian(snes, Jfd, Jfd, SNESComputeJacobianDefault, this); + SNESSetJacobian(snes, Jfd, Jfd, ComputeJacobianDefaultMaybeExport, this); } MatSetOption(Jfd, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); @@ -364,7 +376,99 @@ SNESSolver::SNESSolver(Options* opts) .withDefault(100.)), asinh_vars((*options)["asinh_vars"] .doc("Apply asinh() to all variables?") - .withDefault(false)) {} + .withDefault(false)), + save_jacobian((*options)["save_jacobian"] + .doc("Save Jacobian matrices for diagnostics?") + .withDefault(false)), + jacobian_export_kind((*options)["jacobian_export_kind"] + .doc("Which Jacobian to save: system, scaled, rhs") + .withDefault(JacobianExportKind::system)), + jacobian_export_prefix( + (*options)["jacobian_export_prefix"] + .doc("Prefix for saved Jacobian matrix and metadata files") + .withDefault("jacobian")), + jacobian_export_format((*options)["jacobian_export_format"] + .doc("Format for saved Jacobian matrices: binary, ascii") + .withDefault(PetscMatrixExportFormat::binary)) {} + +std::string SNESSolver::getJacobianExportStem(JacobianExportKind kind) { + return fmt::format("{}_{}_{:06d}", jacobian_export_prefix, toString(kind), + jacobian_export_counter++); +} + +std::string SNESSolver::getJacobianMatrixFilename(const std::string& stem) const { + return stem + + (jacobian_export_format == PetscMatrixExportFormat::binary ? ".dat" : ".txt"); +} + +PetscErrorCode +SNESSolver::exportMatrixAndMetadata(const PetscPreconditioner& preconditioner, + const std::string& stem) { + if (!jacobian_metadata_written) { + writeJacobianMetadataJson(jacobian_export_prefix + "_metadata.json", "snes"); + jacobian_metadata_written = true; + } + + PetscCall( + preconditioner.saveMatrix(getJacobianMatrixFilename(stem), jacobian_export_format)); + PetscFunctionReturn(PETSC_SUCCESS); +} + +PetscErrorCode SNESSolver::saveDiagnosticJacobian(JacobianExportKind kind, Vec x_solver) { + PetscPreconditioner diagnostic_preconditioner; + Field3D index = globalIndex(0); + PetscCall(diagnostic_preconditioner.createJacobianPattern( + index, *options, nlocal, n2Dvars(), n3Dvars(), BoutComm::get())); + + if (kind == JacobianExportKind::rhs) { + PetscCall(diagnostic_preconditioner.updateColoring(FormRawFunctionForColoring, this)); + } else { + PetscCall( + diagnostic_preconditioner.updateColoring(FormScaledFunctionForColoring, this)); + } + + Vec x_evaluate = x_solver; + Vec physical_x{nullptr}; + if (kind == JacobianExportKind::rhs) { + PetscCall(VecDuplicate(x_solver, &physical_x)); + PetscCall(toPhysicalState(x_solver, physical_x)); + x_evaluate = physical_x; + } + + Mat diagnostic_jacobian = diagnostic_preconditioner.jacobian(); + PetscCall(MatZeroEntries(diagnostic_jacobian)); + PetscCall(SNESComputeJacobianDefaultColor(snes, x_evaluate, diagnostic_jacobian, + diagnostic_jacobian, + diagnostic_preconditioner.coloring())); + PetscCall( + exportMatrixAndMetadata(diagnostic_preconditioner, getJacobianExportStem(kind))); + + if (physical_x != nullptr) { + PetscCall(VecDestroy(&physical_x)); + } + + PetscFunctionReturn(PETSC_SUCCESS); +} + +PetscErrorCode SNESSolver::maybeExportJacobian(Mat system_jacobian, Vec x_solver) { + if (!save_jacobian) { + PetscFunctionReturn(PETSC_SUCCESS); + } + + if (jacobian_export_kind == JacobianExportKind::system) { + if (!jacobian_metadata_written) { + writeJacobianMetadataJson(jacobian_export_prefix + "_metadata.json", "snes"); + jacobian_metadata_written = true; + } + PetscCall(PetscPreconditioner::saveMatrix( + system_jacobian, + getJacobianMatrixFilename(getJacobianExportStem(jacobian_export_kind)), + jacobian_export_format)); + PetscFunctionReturn(PETSC_SUCCESS); + } + + PetscFunctionReturn(saveDiagnosticJacobian(jacobian_export_kind, x_solver)); +} int SNESSolver::init() { Solver::init(); @@ -1128,13 +1232,13 @@ PetscErrorCode SNESSolver::updateResiduals(Vec x) { const BoutReal* current_residual = nullptr; if (diagnose) { // Call RHS function to get time derivatives - PetscCall(rhs_function(x, deriv, false)); + PetscCall(scaled_rhs_function(x, deriv, false)); // Reading the residual vectors PetscCall(VecGetArrayRead(deriv, ¤t_residual)); } else { // Call RHS function to get time derivatives - PetscCall(rhs_function(x, snes_f, false)); + PetscCall(scaled_rhs_function(x, snes_f, false)); // Reading the residual vectors PetscCall(VecGetArrayRead(snes_f, ¤t_residual)); @@ -1415,34 +1519,34 @@ BoutReal SNESSolver::updatePseudoTimestep(BoutReal previous_timestep, throw BoutException("SNESSolver::updatePseudoTimestep invalid BoutPTCStrategy"); } -PetscErrorCode SNESSolver::rhs_function(Vec x, Vec f, bool linear) { - // Get data from PETSc into BOUT++ fields +PetscErrorCode SNESSolver::toPhysicalState(Vec x, Vec physical_x) { if (scale_vars) { - // scaled_x <- x * var_scaling_factors - PetscCall(VecPointwiseMult(scaled_x, x, var_scaling_factors)); - } else if (asinh_vars) { - PetscCall(VecCopy(x, scaled_x)); + PetscCall(VecPointwiseMult(physical_x, x, var_scaling_factors)); } else { - scaled_x = x; + PetscCall(VecCopy(x, physical_x)); } if (asinh_vars) { PetscInt size; - PetscCall(VecGetLocalSize(scaled_x, &size)); + PetscCall(VecGetLocalSize(physical_x, &size)); - BoutReal* scaled_data = nullptr; - PetscCall(VecGetArray(scaled_x, &scaled_data)); + BoutReal* physical_data = nullptr; + PetscCall(VecGetArray(physical_x, &physical_data)); for (PetscInt i = 0; i != size; ++i) { - scaled_data[i] = asinh_scale * std::sinh(scaled_data[i]); + physical_data[i] = asinh_scale * std::sinh(physical_data[i]); } - PetscCall(VecRestoreArray(scaled_x, &scaled_data)); + PetscCall(VecRestoreArray(physical_x, &physical_data)); } + return PETSC_SUCCESS; +} + +PetscErrorCode SNESSolver::raw_rhs_function(Vec x, Vec f, bool linear) { const BoutReal* xdata = nullptr; - PetscCall(VecGetArrayRead(scaled_x, &xdata)); + PetscCall(VecGetArrayRead(x, &xdata)); // const_cast needed due to load_vars API. Not writing to xdata. load_vars(const_cast(xdata)); - PetscCall(VecRestoreArrayRead(scaled_x, &xdata)); + PetscCall(VecRestoreArrayRead(x, &xdata)); try { // Call RHS function @@ -1460,6 +1564,18 @@ PetscErrorCode SNESSolver::rhs_function(Vec x, Vec f, bool linear) { BoutReal* fdata = nullptr; PetscCall(VecGetArray(f, &fdata)); save_derivs(fdata); + PetscCall(VecRestoreArray(f, &fdata)); + + return PETSC_SUCCESS; +} + +PetscErrorCode SNESSolver::scaled_rhs_function(Vec x, Vec f, bool linear) { + if (!scale_vars && !asinh_vars) { + return raw_rhs_function(x, f, linear); + } + + PetscCall(toPhysicalState(x, scaled_x)); + PetscCall(raw_rhs_function(scaled_x, f, linear)); if (asinh_vars) { // Modify time-derivatives for asinh(var) using chain rule @@ -1472,14 +1588,15 @@ PetscErrorCode SNESSolver::rhs_function(Vec x, Vec f, bool linear) { PetscCall(VecGetLocalSize(f, &size)); const BoutReal* scaled_data = nullptr; PetscCall(VecGetArrayRead(scaled_x, &scaled_data)); + BoutReal* fdata = nullptr; + PetscCall(VecGetArray(f, &fdata)); for (PetscInt i = 0; i != size; ++i) { fdata[i] /= std::sqrt(SQ(scaled_data[i]) + SQ(asinh_scale)); } + PetscCall(VecRestoreArray(f, &fdata)); PetscCall(VecRestoreArrayRead(scaled_x, &scaled_data)); } - PetscCall(VecRestoreArray(f, &fdata)); - if (scale_vars) { PetscCall(VecPointwiseDivide(f, f, var_scaling_factors)); } @@ -1490,7 +1607,7 @@ PetscErrorCode SNESSolver::rhs_function(Vec x, Vec f, bool linear) { PetscErrorCode SNESSolver::snes_function(Vec x, Vec f, bool linear) { // Call the RHS function - if (rhs_function(x, f, linear) != PETSC_SUCCESS) { + if (scaled_rhs_function(x, f, linear) != PETSC_SUCCESS) { // Tell SNES that the input was out of domain SNESSetFunctionDomainError(snes); // Note: Returning non-zero error here leaves vectors in locked state @@ -1664,7 +1781,20 @@ PetscErrorCode ComputeJacobianScaledColor(SNES snes, Vec x1, Mat Jac, Mat Jac_ne CHKERRQ(err); // Call the SNESSolver function - return fctx->scaleJacobian(Jac_new); + PetscCall(fctx->scaleJacobian(Jac_new)); + PetscFunctionReturn(fctx->maybeExportJacobian(Jac_new, x1)); +} + +PetscErrorCode ComputeJacobianDefaultMaybeExport(SNES snes, Vec x1, Mat Jac, Mat Jac_new, + void* ctx) { + PetscErrorCode err = SNESComputeJacobianDefault(snes, x1, Jac, Jac_new, ctx); + CHKERRQ(err); + + if ((err != 0) or (ctx == nullptr)) { + return err; + } + + PetscFunctionReturn(static_cast(ctx)->maybeExportJacobian(Jac_new, x1)); } } // namespace diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index c8042cf2dc..d45ab1cdaf 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -110,6 +110,18 @@ public: /// finite difference approximated Jacobian. PetscErrorCode scaleJacobian(Mat Jac_new); + /// Convert solver coordinates into the physical variables used by the model. + PetscErrorCode toPhysicalState(Vec x, Vec physical_x); + + /// Call the physics model RHS function on a vector of physical variables. + PetscErrorCode raw_rhs_function(Vec x, Vec f, bool linear); + + /// Apply solver-coordinate transforms, call the raw RHS, and transform the + /// resulting derivatives back into solver coordinates. + PetscErrorCode scaled_rhs_function(Vec x, Vec f, bool linear); + + PetscErrorCode maybeExportJacobian(Mat system_jacobian, Vec x_solver); + /// Save diagnostics to output void outputVars(Options& output_options, bool save_repeat = true) override; @@ -121,13 +133,11 @@ private: /// Rescale state (snes_x) so that all quantities are around 1. If /// quantities are near zero then RTOL is used. PetscErrorCode rescale(); - - /// Call the physics model RHS function - /// - /// @param[in] x The state vector. Will be scaled if scale_vars=true - /// @param[out] f The vector for the result f(x) - /// @param[in] linear Specifies that the SNES solver is in a linear (KSP) inner loop - PetscErrorCode rhs_function(Vec x, Vec f, bool linear); + PetscErrorCode saveDiagnosticJacobian(JacobianExportKind kind, Vec x_solver); + std::string getJacobianExportStem(JacobianExportKind kind); + std::string getJacobianMatrixFilename(const std::string& stem) const; + PetscErrorCode exportMatrixAndMetadata(const PetscPreconditioner& preconditioner, + const std::string& stem); BoutSnesOutput output_trigger; ///< Sets when outputs are written @@ -278,6 +288,13 @@ private: bool asinh_vars; ///< Evolve asinh(vars) to compress magnitudes while preserving signs const BoutReal asinh_scale = 1e-5; // Scale below which asinh response becomes ~linear + bool save_jacobian; ///< Save Jacobian matrices for diagnostics + JacobianExportKind jacobian_export_kind; ///< Which Jacobian to save + std::string jacobian_export_prefix; ///< Prefix for Jacobian matrix/metadata outputs + PetscMatrixExportFormat jacobian_export_format; ///< Output format for matrix save + int jacobian_export_counter{0}; ///< Running counter for saved Jacobians + bool jacobian_metadata_written{false}; ///< Has the shared JSON metadata been written? + std::vector resid_2d; ///< Storage for residuals of SNES solve, unpacked from snes_f std::vector diff --git a/src/solver/petsc_preconditioner.cxx b/src/solver/petsc_preconditioner.cxx index 508d001fd3..063d85812e 100644 --- a/src/solver/petsc_preconditioner.cxx +++ b/src/solver/petsc_preconditioner.cxx @@ -6,6 +6,7 @@ #include "bout/assert.hxx" #include "bout/boutcomm.hxx" +#include "bout/boutexception.hxx" #include "bout/field3d.hxx" #include "bout/globals.hxx" #include "bout/mesh.hxx" @@ -16,6 +17,7 @@ #include #include #include +#include #include #include @@ -71,6 +73,31 @@ void PetscPreconditioner::reset() { } } +PetscErrorCode PetscPreconditioner::saveMatrix(Mat matrix, const std::string& filename, + PetscMatrixExportFormat format) { + if (matrix == nullptr) { + throw BoutException("Cannot save Jacobian matrix: matrix has not been created yet"); + } + + PetscViewer viewer{nullptr}; + if (format == PetscMatrixExportFormat::binary) { + PetscCall(PetscViewerBinaryOpen(BoutComm::get(), filename.c_str(), FILE_MODE_WRITE, + &viewer)); + } else { + PetscCall(PetscViewerASCIIOpen(BoutComm::get(), filename.c_str(), &viewer)); + } + + PetscCall(MatView(matrix, viewer)); + PetscCall(PetscViewerDestroy(&viewer)); + + PetscFunctionReturn(PETSC_SUCCESS); +} + +PetscErrorCode PetscPreconditioner::saveMatrix(const std::string& filename, + PetscMatrixExportFormat format) const { + return saveMatrix(Jfd, filename, format); +} + PetscErrorCode PetscPreconditioner::createJacobianPattern(Field3D& index, Options& options, PetscInt nlocal, int n2d, diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx index 1f7b654427..133451c9e7 100644 --- a/src/solver/solver.cxx +++ b/src/solver/solver.cxx @@ -51,6 +51,7 @@ #include #include +#include #include #include #include @@ -79,6 +80,44 @@ int* Solver::pargc = nullptr; char*** Solver::pargv = nullptr; +namespace { +std::string jsonEscape(const std::string& input) { + std::string escaped; + escaped.reserve(input.size()); + + for (const char ch : input) { + switch (ch) { + case '\\': + escaped += "\\\\"; + break; + case '"': + escaped += "\\\""; + break; + case '\b': + escaped += "\\b"; + break; + case '\f': + escaped += "\\f"; + break; + case '\n': + escaped += "\\n"; + break; + case '\r': + escaped += "\\r"; + break; + case '\t': + escaped += "\\t"; + break; + default: + escaped += ch; + break; + } + } + + return escaped; +} +} // namespace + /************************************************************************** * Constructor **************************************************************************/ @@ -94,6 +133,10 @@ Solver::Solver(Options* opts) "timestep, to make it easier to concatenate output " "data sets in time") .withDefault(false)), + save_jacobian_index_base( + (*options)["save_jacobian_index_base"] + .doc("Write the base global index field used for Jacobian diagnostics") + .withDefault(false)), is_nonsplit_model_diffusive( (*options)["is_nonsplit_model_diffusive"] .doc("If not a split operator, treat RHS as diffusive?") @@ -705,6 +748,11 @@ void Solver::outputVars(Options& output_options, bool save_repeat) { "or the previous run did not have a run_id.") .assignRepeat(run_restart_from, "t", save_repeat and save_repeat_run_id, "Solver"); + if (save_jacobian_index_base) { + output_options["jacobian_index_base"].assignRepeat(jacobianIndexBase(), "t", + save_repeat, "Solver"); + } + // Add 2D and 3D evolving fields to output file for (const auto& f : f2d) { // Add to dump file (appending) @@ -1207,6 +1255,91 @@ Field3D Solver::globalIndex(int localStart) { return index; } +Field3D Solver::jacobianIndexBase(int localStart) { return globalIndex(localStart); } + +std::vector Solver::getJacobianMetadata2D() const { + std::vector metadata; + metadata.reserve(f2d.size()); + + for (int i = 0; i < static_cast(f2d.size()); ++i) { + metadata.push_back(JacobianVariableMetadata{i, f2d[i].name, toString(f2d[i].location), + f2d[i].evolve_bndry, f2d[i].constraint, + f2d[i].description}); + } + + return metadata; +} + +std::vector Solver::getJacobianMetadata3D() const { + std::vector metadata; + metadata.reserve(f3d.size()); + + for (int i = 0; i < static_cast(f3d.size()); ++i) { + metadata.push_back(JacobianVariableMetadata{i, f3d[i].name, toString(f3d[i].location), + f3d[i].evolve_bndry, f3d[i].constraint, + f3d[i].description}); + } + + return metadata; +} + +Solver::JacobianMetadata +Solver::getJacobianMetadata(const std::string& solver_name) const { + return JacobianMetadata{1, + solver_name, + n2Dvars(), + n3Dvars(), + getJacobianMetadata2D(), + getJacobianMetadata3D(), + "For each (x,y): 2D variables at z=0, then 3D variables for " + "z=0..Nz-1; evolved boundary points precede RGN_NOBNDRY"}; +} + +void Solver::writeJacobianMetadataJson(const std::string& filename, + const std::string& solver_name) const { + if (MYPE != 0) { + return; + } + + const auto metadata = getJacobianMetadata(solver_name); + std::ofstream json_file(filename); + if (!json_file.is_open()) { + throw BoutException("Failed to open Jacobian metadata file '{}'", filename); + } + + auto write_variables = [&](const std::vector& variables, + const char* name) { + json_file << " \"" << name << "\": [\n"; + for (std::size_t i = 0; i < variables.size(); ++i) { + const auto& variable = variables[i]; + json_file << " {\"offset\": " << variable.offset << ", " + << "\"name\": \"" << jsonEscape(variable.name) << "\", " + << "\"location\": \"" << jsonEscape(variable.location) << "\", " + << "\"evolve_bndry\": " << (variable.evolve_bndry ? "true" : "false") + << ", " + << "\"constraint\": " << (variable.constraint ? "true" : "false") << ", " + << "\"description\": \"" << jsonEscape(variable.description) << "\"}"; + if (i + 1 != variables.size()) { + json_file << ","; + } + json_file << "\n"; + } + json_file << " ]"; + }; + + json_file << "{\n" + << " \"format_version\": " << metadata.format_version << ",\n" + << " \"solver\": \"" << jsonEscape(metadata.solver_name) << "\",\n" + << " \"n2d\": " << metadata.n2d << ",\n" + << " \"n3d\": " << metadata.n3d << ",\n"; + write_variables(metadata.variables_2d, "variables_2d"); + json_file << ",\n"; + write_variables(metadata.variables_3d, "variables_3d"); + json_file << ",\n" + << " \"ordering\": \"" << jsonEscape(metadata.ordering) << "\"\n" + << "}\n"; +} + /************************************************************************** * Running user-supplied functions **************************************************************************/