sim_utils.cc

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
// SPDX-License-Identifier: LicenseRef-APL-1.1
 
#include "backend.h"
#include "basics.h"
#include "sim_utils.h"
#include "fem.h"
#include "physics.h"
#include "async_io.h"
#include "SF_init.h"
#include "opencarp_types.h"
#ifdef WITH_LEADFIELD
#include "leadfield.h"
#endif
#ifdef WITH_POWERCAPPING
#include <vector>
#include "powercapping.h"
#endif
 
#include <cctype>
#include <cstdlib>
#include <libgen.h>
#include <iomanip>
#include <map>
#include <set>
#include <sstream>
#include <string>
#include <sys/wait.h>
#include <unistd.h>
#include <vector>
 
#include "openCARP_schema.hpp"
#include "runtime.hpp"
#include "snapshot_file_io.hpp"
 
namespace opencarp {
 
namespace {
 
enum class ParserCompareMode {
  Off,
  Warn,
  Strict,
};
 
enum class ParserFallbackMode {
  Off,
  Legacy,
};
 
struct RuntimeCompatOptions {
  ParserFallbackMode fallback_mode = ParserFallbackMode::Off;
};
 
struct LegacyCompareInput {
  bool                     available = false;
  bool                     warn_when_unavailable = true;
  std::vector<std::string> runtime_args;
  std::string              unavailable_reason;
};
 
struct LegacySnapshotHelperRunResult {
  bool exited = false;
  int  exit_status = -1;
  bool signaled = false;
  int  signal = 0;
};
 
std::vector<paramschema::CitationSuggestion> active_citation_suggestions;
 
std::string trim_copy(const std::string& value);
std::string to_lower_ascii(std::string value);
bool build_legacy_compare_input(int argc, char** argv, LegacyCompareInput* input, std::string* error);
 
std::string join_paths_for_message(const std::vector<std::string>& paths)
{
  if (paths.empty()) {
    return std::string();
  }
 
  std::ostringstream joined;
  for (std::size_t i = 0; i < paths.size(); ++i) {
    if (i != 0) {
      joined << ", ";
    }
    joined << paths[i];
  }
  return joined.str();
}
 
bool match_long_option(const std::string& token, const char* option, std::string* attached_value)
{
  *attached_value = std::string();
  if (token == option) {
    return true;
  }
 
  const std::string prefix = std::string(option) + "=";
  if (token.size() > prefix.size() && token.compare(0, prefix.size(), prefix) == 0) {
    *attached_value = token.substr(prefix.size());
    return true;
  }
 
  return false;
}
 
bool is_help_topic_candidate(const char* token)
{
  return token != NULL && token[0] != '\0' && token[0] != '-' && token[0] != '+';
}
 
bool is_long_option_argument_error(const std::string& token,
                                   const char*        option,
                                   const std::string& attached_value,
                                   std::string*       error)
{
  if (attached_value.empty()) {
    return false;
  }
 
  *error = "Unexpected argument for " + std::string(option) + " in '" + token + "'";
  return true;
}
 
bool normalize_runtime_args(int                       argc,
                            char**                    argv,
                            std::vector<std::string>* normalized,
                            RuntimeCompatOptions*      compat,
                            std::string*              error)
{
  normalized->clear();
  compat->fallback_mode = ParserFallbackMode::Off;
 
  if (argc <= 0 || argv == NULL || argv[0] == NULL) {
    *error = "Missing program name";
    return false;
  }
 
  normalized->push_back(argv[0]);
 
  for (int i = 1; i < argc; ++i) {
    const std::string token = argv[i];
    if (token == "+") {
      break;
    }
 
    std::string attached_value;
 
    if (token == "+Help" || match_long_option(token, "--help", &attached_value)) {
      std::string topic = "PrM";
      if (!attached_value.empty()) {
        topic = attached_value;
      } else if (i + 1 < argc && is_help_topic_candidate(argv[i + 1])) {
        topic = argv[++i];
      }
      normalized->push_back("+Help");
      normalized->push_back(topic);
      break;
    }
 
    if (token == "+Doc" || match_long_option(token, "--doc", &attached_value)) {
      std::string topic = "ALL";
      if (!attached_value.empty()) {
        topic = attached_value;
      } else if (i + 1 < argc && is_help_topic_candidate(argv[i + 1])) {
        topic = argv[++i];
      }
      normalized->push_back("+Doc");
      normalized->push_back(topic);
      break;
    }
 
    if (token == "+Default" || match_long_option(token, "--default", &attached_value)) {
      if (token != "+Default" && is_long_option_argument_error(token, "--default", attached_value, error)) {
        return false;
      }
      normalized->push_back("+Default");
      break;
    }
 
    if (token == "+Run" || match_long_option(token, "--run", &attached_value)) {
      if (token != "+Run" && is_long_option_argument_error(token, "--run", attached_value, error)) {
        return false;
      }
      normalized->push_back("+Run");
      continue;
    }
 
    if (token == "+I" || match_long_option(token, "--interactive", &attached_value)) {
      if (!attached_value.empty()) {
        *error = "Unexpected argument for --interactive in '" + token + "'";
      } else {
        *error = "Unsupported option " + token + " (interactive mode is not available)";
      }
      return false;
    }
 
    if (match_long_option(token, "--param-fallback", &attached_value)) {
      std::string mode = attached_value;
      if (mode.empty()) {
        if (i + 1 >= argc) {
          *error = "Missing argument after --param-fallback";
          return false;
        }
        mode = argv[++i];
      }
 
      if (to_lower_ascii(trim_copy(mode)) != "legacy") {
        *error = "Unsupported value '" + mode + "' for --param-fallback (expected legacy)";
        return false;
      }
 
      compat->fallback_mode = ParserFallbackMode::Legacy;
      continue;
    }
 
    if (token == "+F" || match_long_option(token, "--file", &attached_value)) {
      std::string filename = attached_value;
      if (filename.empty()) {
        if (i + 1 >= argc) {
          *error = "Missing filename after " + token;
          return false;
        }
        filename = argv[++i];
      }
      normalized->push_back("+F");
      normalized->push_back(filename);
      continue;
    }
 
    if (token == "+Save" || match_long_option(token, "--save", &attached_value)) {
      std::string filename = attached_value;
      if (filename.empty()) {
        if (i + 1 >= argc) {
          *error = "Missing argument after " + token;
          return false;
        }
        filename = argv[++i];
      }
      normalized->push_back("+Save");
      normalized->push_back(filename);
      continue;
    }
 
    normalized->push_back(token);
  }
 
  return true;
}
 
bool parse_option_argument(int*                    index,
                           int                     argc,
                           char**                  argv,
                           const std::string&      token,
                           const std::string&      attached_value,
                           const char*             option_name,
                           std::string*            value,
                           std::string*            error)
{
  if (!attached_value.empty()) {
    *value = attached_value;
    return true;
  }
  if (*index + 1 >= argc) {
    *error = "Missing argument after " + std::string(option_name);
    return false;
  }
  *value = argv[++(*index)];
  return true;
}
 
bool filename_has_suffix(const std::string& filename, const char* suffix)
{
  const std::string normalized_filename = to_lower_ascii(trim_copy(filename));
  const std::string normalized_suffix = to_lower_ascii(std::string(suffix == NULL ? "" : suffix));
  return normalized_filename.size() >= normalized_suffix.size() &&
         normalized_filename.compare(normalized_filename.size() - normalized_suffix.size(),
                                     normalized_suffix.size(),
                                     normalized_suffix) == 0;
}
 
bool build_legacy_compare_input(int argc, char** argv, LegacyCompareInput* input, std::string* error)
{
  input->available = false;
  input->warn_when_unavailable = true;
  input->runtime_args.clear();
  input->unavailable_reason.clear();
 
  if (argc <= 0 || argv == NULL || argv[0] == NULL) {
    *error = "Missing program name";
    return false;
  }
 
  input->runtime_args.push_back(argv[0]);
  bool saw_passthrough_token_before_file = false;
 
  for (int i = 1; i < argc; ++i) {
    const std::string token = argv[i];
    if (token == "+") {
      break;
    }
 
    std::string attached_value;
 
    if (match_long_option(token, "--param-fallback", &attached_value)) {
      std::string ignored;
      if (!parse_option_argument(&i, argc, argv, token, attached_value, "--param-fallback", &ignored, error)) {
        return false;
      }
      continue;
    }
 
    if (token == "+Save" || match_long_option(token, "--save", &attached_value)) {
      std::string ignored;
      if (!parse_option_argument(&i, argc, argv, token, attached_value, token.c_str(), &ignored, error)) {
        return false;
      }
      continue;
    }
 
    if (token == "+Help" || match_long_option(token, "--help", &attached_value)) {
      std::string topic = "PrM";
      if (!attached_value.empty()) {
        topic = attached_value;
      } else if (i + 1 < argc && is_help_topic_candidate(argv[i + 1])) {
        topic = argv[++i];
      }
      input->runtime_args.push_back("+Help");
      input->runtime_args.push_back(topic);
      input->available = true;
      return true;
    }
 
    if (token == "+Doc" || match_long_option(token, "--doc", &attached_value)) {
      std::string topic = "ALL";
      if (!attached_value.empty()) {
        topic = attached_value;
      } else if (i + 1 < argc && is_help_topic_candidate(argv[i + 1])) {
        topic = argv[++i];
      }
      input->runtime_args.push_back("+Doc");
      input->runtime_args.push_back(topic);
      input->available = true;
      return true;
    }
 
    if (token == "+Default" || match_long_option(token, "--default", &attached_value)) {
      if (token != "+Default" && is_long_option_argument_error(token, "--default", attached_value, error)) {
        return false;
      }
      input->runtime_args.push_back("+Default");
      continue;
    }
 
    if (token == "+Run" || match_long_option(token, "--run", &attached_value)) {
      if (token != "+Run" && is_long_option_argument_error(token, "--run", attached_value, error)) {
        return false;
      }
      input->runtime_args.push_back("+Run");
      continue;
    }
 
    if (token == "+F" || match_long_option(token, "--file", &attached_value)) {
      std::string filename;
      if (!parse_option_argument(&i, argc, argv, token, attached_value, token.c_str(), &filename, error)) {
        return false;
      }
 
      if (saw_passthrough_token_before_file) {
        input->unavailable_reason =
            "legacy compare is unavailable when direct parameter arguments precede a .par input file";
        input->runtime_args.clear();
        input->runtime_args.push_back(argv[0]);
        return true;
      }
 
      if (!filename_has_suffix(filename, ".par")) {
        input->unavailable_reason =
            "legacy compare is unavailable for non-.par input file '" + filename + "'";
        input->runtime_args.clear();
        input->runtime_args.push_back(argv[0]);
        return true;
      }
 
      input->runtime_args.push_back("+F");
      input->runtime_args.push_back(filename);
      continue;
    }
 
    input->runtime_args.push_back(token);
    saw_passthrough_token_before_file = true;
  }
 
  input->available = input->runtime_args.size() > 1;
  if (!input->available && input->unavailable_reason.empty()) {
    input->unavailable_reason = "unable to reconstruct a legacy-compatible parameter input";
  }
  return true;
}
 
std::string trim_copy(const std::string& value)
{
  std::string::size_type first = 0;
  while (first < value.size() && std::isspace(static_cast<unsigned char>(value[first]))) {
    ++first;
  }
 
  std::string::size_type last = value.size();
  while (last > first && std::isspace(static_cast<unsigned char>(value[last - 1]))) {
    --last;
  }
 
  return value.substr(first, last - first);
}
 
std::string to_lower_ascii(std::string value)
{
  for (std::string::size_type i = 0; i < value.size(); ++i) {
    value[i] = static_cast<char>(std::tolower(static_cast<unsigned char>(value[i])));
  }
  return value;
}
 
ParserCompareMode parser_compare_mode()
{
  const char* raw = std::getenv("OPENCARP_PARAM_COMPARE");
  if (raw == NULL) {
    return ParserCompareMode::Strict;
  }
 
  const std::string normalized = to_lower_ascii(trim_copy(raw));
  if (normalized.empty() || normalized == "1" || normalized == "on" || normalized == "true" ||
      normalized == "yes" || normalized == "strict" || normalized == "fail" || normalized == "error") {
    return ParserCompareMode::Strict;
  }
  if (normalized == "warn") {
    return ParserCompareMode::Warn;
  }
  if (normalized == "0" || normalized == "off" || normalized == "false" || normalized == "no") {
    return ParserCompareMode::Off;
  }
 
  return ParserCompareMode::Strict;
}
 
ParserFallbackMode parser_fallback_mode()
{
  const char* raw = std::getenv("OPENCARP_PARAM_FALLBACK");
  if (raw == NULL) {
    return ParserFallbackMode::Off;
  }
 
  const std::string normalized = to_lower_ascii(trim_copy(raw));
  if (normalized.empty() || normalized == "0" || normalized == "off" || normalized == "false" ||
      normalized == "no") {
    return ParserFallbackMode::Off;
  }
  if (normalized == "legacy") {
    return ParserFallbackMode::Legacy;
  }
 
  return ParserFallbackMode::Off;
}
 
void populate_arg_pointers(const std::vector<std::string>& values, std::vector<char*>* argv)
{
  argv->assign(values.size(), NULL);
  for (std::size_t i = 0; i < values.size(); ++i) {
    (*argv)[i] = const_cast<char*>(values[i].c_str());
  }
}
 
void print_lines(FILE* stream, const char* label, const std::vector<std::string>& lines)
{
  for (std::size_t i = 0; i < lines.size(); ++i) {
    fprintf(stream, "%s%s\n", label, lines[i].c_str());
  }
}
 
void clear_active_citation_suggestions()
{
  active_citation_suggestions.clear();
}
 
void store_active_citation_suggestions(const std::vector<paramschema::CitationSuggestion>& suggestions)
{
  active_citation_suggestions = suggestions;
}
 
void print_active_citation_suggestions()
{
  if (active_citation_suggestions.empty() || get_rank() != 0) {
    return;
  }
 
  const std::vector<std::string> lines =
      paramschema::format_citation_suggestions(active_citation_suggestions);
  if (lines.empty()) {
    return;
  }
 
  std::fprintf(stdout, "\nIf you publish studies based on this simulation, the following references are likely relevant:\n");
  for (std::size_t i = 0; i < lines.size(); ++i) {
    std::fprintf(stdout, "  %s\n", lines[i].c_str());
  }
  std::fprintf(stdout, "  You can use https://citation.doi.org to format references in your preferred format or https://www.doi2bib.org to turn the DOIs into BibTeX format.\n");
  std::fprintf(stdout, "\nWe'd be happy to also see your experiment published to boost its impact and foster reproducibility.\nSee https://opencarp.org/community/upload-experiment for details.\n\n");
 
  std::fflush(stdout);
}
 
void print_parser_runtime_state(const std::vector<std::string>& runtime_args)
{
  const auto&       schema = paramschema::openCARP_schema();
  const std::string program_name =
      runtime_args.empty() || runtime_args[0].empty() ? "openCARP" : runtime_args[0];
 
  std::string              rendered;
  std::vector<std::string> errors;
  if (!paramschema::render_save_text(schema, program_name, &rendered, &errors)) {
    print_lines(stderr, "parameter parser warning: ", errors);
    return;
  }
 
  fputs(rendered.c_str(), stderr);
}
 
paramschema::ExecutionResult execute_parser_runtime_args(const std::vector<std::string>& runtime_args,
                                                         const bool                      allow_save)
{
  std::vector<char*> argv;
  populate_arg_pointers(runtime_args, &argv);
 
  const auto& schema = paramschema::openCARP_schema();
 
  paramschema::ExecutionOptions options;
  options.allow_save = allow_save;
  return paramschema::execute_legacy_cli(schema, static_cast<int>(argv.size()), argv.data(), options);
}
 
bool apply_parser_runtime_args(const std::vector<std::string>& runtime_args,
                               const bool                      allow_save,
                               paramschema::ExecutionResult*   executed_out)
{
  clear_active_citation_suggestions();
  const paramschema::ExecutionResult executed = execute_parser_runtime_args(runtime_args, allow_save);
  if (executed_out != NULL) {
    *executed_out = executed;
  }
  print_lines(stderr, "parameter parser warning: ", executed.warnings);
  if (!executed.rendered_output.empty()) {
    fputs(executed.rendered_output.c_str(), stdout);
  }
  if (!executed.errors.empty() || executed.status == paramschema::ExecutionStatus::Fatal) {
    print_lines(stderr, "parameter parser error: ", executed.errors);
    return false;
  }
 
  if (param_globals::output_setup) {
    print_parser_runtime_state(runtime_args);
  }
 
  if (executed.status == paramschema::ExecutionStatus::Help) {
    exit(EXIT_SUCCESS);
  }
 
  store_active_citation_suggestions(executed.citations);
  return true;
}
 
std::string parent_directory(const std::string& path)
{
  const std::string::size_type slash = path.rfind('/');
  if (slash == std::string::npos) {
    return ".";
  }
  if (slash == 0) {
    return "/";
  }
  return path.substr(0, slash);
}
 
std::string join_path(const std::string& left, const std::string& right)
{
  if (left.empty() || left == ".") {
    return right;
  }
  if (!left.empty() && left[left.size() - 1] == '/') {
    return left + right;
  }
  return left + "/" + right;
}
 
bool find_executable_on_path(const std::string& name, std::string* resolved_path)
{
  if (name.empty() || name.find('/') != std::string::npos) {
    return false;
  }
 
  const char* path_env = std::getenv("PATH");
  if (path_env == NULL || path_env[0] == '\0') {
    return false;
  }
 
  const std::string path_list = path_env;
  std::string::size_type start = 0;
  while (start <= path_list.size()) {
    std::string::size_type end = path_list.find(':', start);
    if (end == std::string::npos) {
      end = path_list.size();
    }
 
    const std::string directory = path_list.substr(start, end - start);
    const std::string candidate = join_path(directory.empty() ? "." : directory, name);
    if (access(candidate.c_str(), X_OK) == 0) {
      *resolved_path = candidate;
      return true;
    }
 
    if (end == path_list.size()) {
      break;
    }
    start = end + 1;
  }
 
  return false;
}
 
bool resolve_legacy_snapshot_helper(const std::string& program_path, std::string* helper_path)
{
  std::vector<std::string> resolved_program_paths;
  if (!program_path.empty()) {
    resolved_program_paths.push_back(program_path);
  }
 
  if (program_path.find('/') == std::string::npos) {
    std::string resolved_program_path;
    if (find_executable_on_path(program_path, &resolved_program_path)) {
      resolved_program_paths.push_back(resolved_program_path);
    }
  }
 
  for (std::size_t i = 0; i < resolved_program_paths.size(); ++i) {
    const std::string program_dir = parent_directory(resolved_program_paths[i]);
    const std::string parent_dir = parent_directory(program_dir);
 
    const std::vector<std::string> candidates = {
        join_path(program_dir, "param-parser-legacy-snapshot"),
        join_path(join_path(parent_dir, "simulator"), "param-parser-legacy-snapshot"),
    };
 
    for (std::size_t j = 0; j < candidates.size(); ++j) {
      if (access(candidates[j].c_str(), X_OK) == 0) {
        *helper_path = candidates[j];
        return true;
      }
    }
  }
 
  return find_executable_on_path("param-parser-legacy-snapshot", helper_path);
}
 
bool create_temp_output_path(const char* suffix, std::string* path)
{
  char temp_path[128];
  if (suffix != NULL && suffix[0] != '\0') {
    std::snprintf(temp_path, sizeof temp_path, "/tmp/opencarp-parser-compare-XXXXXX%s", suffix);
  } else {
    std::snprintf(temp_path, sizeof temp_path, "/tmp/opencarp-parser-compare-XXXXXX");
  }
 
  const int fd = suffix != NULL && suffix[0] != '\0' ? mkstemps(temp_path, static_cast<int>(std::strlen(suffix))) :
                 mkstemp(temp_path);
  if (fd < 0) {
    return false;
  }
  close(fd);
  *path = temp_path;
  return true;
}
 
bool capture_current_parser_snapshot(paramschema::SnapshotResult* snapshot)
{
  *snapshot = paramschema::snapshot_schema_state(paramschema::openCARP_schema());
  print_lines(stderr, "parameter compare warning: ", snapshot->warnings);
  if (!snapshot->errors.empty()) {
    print_lines(stderr, "parameter compare error: ", snapshot->errors);
    return false;
  }
  return true;
}
 
void maybe_force_test_mismatch(paramschema::SnapshotResult* snapshot)
{
  const char* raw = std::getenv("OPENCARP_PARAM_TEST_FORCE_MISMATCH");
  if (raw == NULL) {
    return;
  }
 
  const std::string normalized = to_lower_ascii(trim_copy(raw));
  if (normalized.empty() || normalized == "0" || normalized == "off" || normalized == "false" ||
      normalized == "no") {
    return;
  }
 
  if (!snapshot->entries.empty()) {
    snapshot->entries[0].value += "__forced_parser_compare_mismatch__";
    return;
  }
 
  paramschema::SnapshotEntry entry;
  entry.path = "buildinfo";
  entry.value = "__forced_parser_compare_mismatch__";
  snapshot->entries.push_back(entry);
}
 
LegacySnapshotHelperRunResult run_legacy_snapshot_helper(const std::string&             helper_path,
                                                         const std::vector<std::string>& runtime_args,
                                                         const std::string&             output_path)
{
  LegacySnapshotHelperRunResult result;
  std::vector<std::string> child_values;
  child_values.reserve(runtime_args.size() + 4);
  child_values.push_back(helper_path);
  child_values.push_back("--snapshot-out");
  child_values.push_back(output_path);
  child_values.push_back("--");
  child_values.insert(child_values.end(), runtime_args.begin(), runtime_args.end());
 
  std::vector<char*> child_argv;
  child_argv.reserve(child_values.size() + 1);
  for (std::size_t i = 0; i < child_values.size(); ++i) {
    child_argv.push_back(const_cast<char*>(child_values[i].c_str()));
  }
  child_argv.push_back(NULL);
 
  const pid_t pid = fork();
  if (pid < 0) {
    std::perror("fork");
    return result;
  }
 
  if (pid == 0) {
    execv(helper_path.c_str(), child_argv.data());
    std::perror("execv");
    _exit(127);
  }
 
  int status = 0;
  if (waitpid(pid, &status, 0) < 0) {
    std::perror("waitpid");
    return result;
  }
 
  if (WIFSIGNALED(status)) {
    result.signaled = true;
    result.signal = WTERMSIG(status);
    return result;
  }
 
  if (WIFEXITED(status)) {
    result.exited = true;
    result.exit_status = WEXITSTATUS(status);
  }
 
  return result;
}
 
bool legacy_snapshot_helper_rejected_input(const LegacySnapshotHelperRunResult& result)
{
  return result.exited && (result.exit_status == 3 || result.exit_status == 4);
}
 
bool legacy_fallback_requested(const RuntimeCompatOptions& compat)
{
  return compat.fallback_mode == ParserFallbackMode::Legacy ||
         parser_fallback_mode() == ParserFallbackMode::Legacy;
}
 
bool should_print_global_message()
{
  if (mpi_runtime_ready()) {
    return get_rank(MPI_COMM_WORLD) == 0;
  }
 
  const char* const rank_variables[] = {
      "OMPI_COMM_WORLD_RANK",
      "PMI_RANK",
      "PMIX_RANK",
      "MV2_COMM_WORLD_RANK",
      "SLURM_PROCID",
      "ALPS_APP_PE",
  };
  for (std::size_t i = 0; i < sizeof rank_variables / sizeof rank_variables[0]; ++i) {
    const char* raw_rank = std::getenv(rank_variables[i]);
    if (raw_rank == NULL || raw_rank[0] == '\0') {
      continue;
    }
 
    char* end = NULL;
    const long rank = std::strtol(raw_rank, &end, 10);
    if (end != raw_rank && end != NULL && end[0] == '\0') {
      return rank == 0;
    }
  }
 
  return true;
}
 
void print_legacy_fallback_workaround()
{
  std::fprintf(stderr,
               "parameter compare error: rerun with OPENCARP_PARAM_FALLBACK=legacy to continue with the legacy parameter state\n");
  std::fprintf(stderr,
               "parameter compare error: or add --param-fallback=legacy to the openCARP command line\n");
}
 
bool restore_legacy_snapshot_state(const paramschema::SnapshotResult& legacy_snapshot)
{
  const paramschema::SnapshotRestoreResult restored =
      paramschema::restore_snapshot_state(paramschema::openCARP_schema(), legacy_snapshot);
  print_lines(stderr, "parameter compare warning: ", restored.warnings);
  if (!restored.errors.empty()) {
    print_lines(stderr, "parameter compare error: ", restored.errors);
    return false;
  }
  return true;
}
 
bool run_parser_legacy_compare(const std::vector<std::string>& runtime_args,
                               const LegacyCompareInput&       legacy_input,
                               const RuntimeCompatOptions&     compat)
{
  const ParserCompareMode mode = parser_compare_mode();
  if (mode == ParserCompareMode::Off) {
    return true;
  }
  const bool fallback_to_legacy = legacy_fallback_requested(compat);
 
  if (runtime_args.empty()) {
    std::fprintf(stderr, "parameter compare error: unable to reconstruct simulator argv\n");
    if (fallback_to_legacy) {
      print_legacy_fallback_workaround();
    }
    return mode != ParserCompareMode::Strict;
  }
 
  paramschema::SnapshotResult parser_snapshot;
  if (!capture_current_parser_snapshot(&parser_snapshot)) {
    std::fprintf(stderr, "parameter compare error: unable to snapshot the parser runtime state\n");
    if (fallback_to_legacy) {
      print_legacy_fallback_workaround();
    }
    return mode != ParserCompareMode::Strict;
  }
  maybe_force_test_mismatch(&parser_snapshot);
 
  std::string helper_path;
  if (!resolve_legacy_snapshot_helper(runtime_args[0], &helper_path)) {
    std::fprintf(stderr, "parameter compare error: unable to locate param-parser-legacy-snapshot\n");
    if (fallback_to_legacy) {
      print_legacy_fallback_workaround();
    } else {
      std::fprintf(stderr,
                   "parameter compare error: rerun with OPENCARP_PARAM_COMPARE=0 to bypass this temporary compatibility gate\n");
    }
    return mode != ParserCompareMode::Strict;
  }
 
  if (!legacy_input.available) {
    const bool print_message = should_print_global_message();
    if (print_message && legacy_input.warn_when_unavailable && !legacy_input.unavailable_reason.empty()) {
      std::fprintf(stderr, "parameter compare warning: %s\n", legacy_input.unavailable_reason.c_str());
    }
    if (print_message && fallback_to_legacy) {
      std::fprintf(stderr,
                   "parameter compare warning: legacy fallback is unavailable because no legacy baseline exists for this input\n");
    }
    return true;
  }
 
  std::string snapshot_path;
  if (!create_temp_output_path("", &snapshot_path)) {
    std::fprintf(stderr, "parameter compare error: unable to create temporary snapshot file\n");
    if (fallback_to_legacy) {
      print_legacy_fallback_workaround();
    } else {
      std::fprintf(stderr,
                   "parameter compare error: rerun with OPENCARP_PARAM_COMPARE=0 to bypass this temporary compatibility gate\n");
    }
    return mode != ParserCompareMode::Strict;
  }
 
  const LegacySnapshotHelperRunResult helper_result =
      run_legacy_snapshot_helper(helper_path, legacy_input.runtime_args, snapshot_path);
  const bool helper_ok = helper_result.exited && helper_result.exit_status == 0;
  paramschema::SnapshotResult legacy_snapshot;
  std::string                 read_error;
  const bool loaded = helper_ok &&
                      paramschema::snapshotio::read_snapshot_file(snapshot_path, &legacy_snapshot, &read_error);
  unlink(snapshot_path.c_str());
 
  if (!loaded) {
    if (legacy_snapshot_helper_rejected_input(helper_result)) {
      std::fprintf(stderr,
                   "parameter compare error: legacy param() rejected the normalized input while the parser runtime accepted it\n");
      std::fprintf(stderr,
                   "parameter compare error: this indicates a parser/legacy validation mismatch, not snapshot helper infrastructure\n");
      std::fprintf(stderr,
                   "parameter compare error: please open an issue and include the triggering command line and parameter files\n");
      std::fprintf(stderr,
                   "https://git.opencarp.org/openCARP/openCARP/-/issues/new?type=ISSUE&initialCreationContext=list-route\n");
      if (fallback_to_legacy) {
        std::fprintf(stderr,
                     "parameter compare error: legacy fallback is unavailable because legacy produced no valid parameter state\n");
      }
    } else if (!read_error.empty()) {
      std::fprintf(stderr, "parameter compare error: %s\n", read_error.c_str());
    } else if (helper_result.signaled) {
      std::fprintf(stderr, "parameter compare error: legacy snapshot helper terminated with signal %d\n",
                   helper_result.signal);
    } else if (helper_result.exited) {
      std::fprintf(stderr, "parameter compare error: legacy snapshot helper exited with status %d\n",
                   helper_result.exit_status);
    } else {
      std::fprintf(stderr, "parameter compare error: unable to capture the legacy parameter state\n");
    }
    if (fallback_to_legacy && !legacy_snapshot_helper_rejected_input(helper_result)) {
      print_legacy_fallback_workaround();
    } else if (!legacy_snapshot_helper_rejected_input(helper_result)) {
      std::fprintf(stderr,
                   "parameter compare error: rerun with OPENCARP_PARAM_COMPARE=0 to bypass this temporary compatibility gate\n");
    }
    return mode != ParserCompareMode::Strict;
  }
 
  const paramschema::SnapshotComparisonResult comparison =
      paramschema::compare_snapshot_results(paramschema::openCARP_schema(), parser_snapshot, legacy_snapshot);
  if (!comparison.errors.empty() || !comparison.mismatches.empty()) {
    print_lines(stderr, "", paramschema::format_snapshot_comparison_report(comparison, "parser compare"));
    std::fprintf(stderr,
                 "parameter compare error: the parser runtime and legacy param() produced different parameter states\n");
    std::fprintf(stderr,
                 "parameter compare error: please open an issue and include the triggering command line and parameter files\n");
    std::fprintf(stderr,
                 "https://git.opencarp.org/openCARP/openCARP/-/issues/new?type=ISSUE&initialCreationContext=list-route\n");
 
    if (fallback_to_legacy) {
      if (!restore_legacy_snapshot_state(legacy_snapshot)) {
        print_legacy_fallback_workaround();
        return false;
      }
      clear_active_citation_suggestions();
      std::fprintf(stderr, "parameter compare warning: continuing with the legacy parameter state\n");
      return true;
    }
 
    print_legacy_fallback_workaround();
    return mode != ParserCompareMode::Strict;
  }
 
  return true;
}
 
}  // namespace
 
static char  input_dir[1024],        // directory from which to read input
             output_dir[1024],       // directory to which to write results
             postproc_dir[1024],     // postprocessing directory
             current_dir[1024];      // current directory
 
namespace {
 
void localize_gridout_nodes(sf_mesh& mesh,
                            SF::vector<mesh_int_t>& selected_nodes,
                            SF::vector<mesh_int_t>& output_idx,
                            bool async)
{
  const int mpi_rank = get_rank();
  const int mpi_size = get_size();
  const SF::vector<mesh_int_t>& nbr = mesh.get_numbering(SF::NBR_SUBMESH);
 
  for(mesh_int_t& node : selected_nodes)
    node = nbr[node];
 
  binary_sort(selected_nodes);
  unique_resize(selected_nodes);
 
  output_idx.resize(0);
 
  if(async) {
    const SF::vector<mesh_int_t>& alg_nod = mesh.pl.algebraic_nodes();
    hashmap::unordered_map<mesh_int_t, mesh_int_t> alg_map;
    alg_map.reserve(alg_nod.size());
 
    for(mesh_int_t local_node : alg_nod)
      alg_map[nbr[local_node]] = local_node;
 
    SF::vector<mesh_int_t> recv_nodes;
    size_t buffsize = 0;
    for(int pid = 0; pid < mpi_size; pid++) {
      if(mpi_rank == pid) {
        recv_nodes = selected_nodes;
        buffsize = recv_nodes.size();
      }
 
      MPI_Bcast(&buffsize, sizeof(size_t), MPI_BYTE, pid, mesh.comm);
      recv_nodes.resize(buffsize);
      MPI_Bcast(recv_nodes.data(), buffsize * sizeof(mesh_int_t), MPI_BYTE, pid, mesh.comm);
 
      for(mesh_int_t node : recv_nodes) {
        auto it = alg_map.find(node);
        if(it != alg_map.end())
          output_idx.push_back(it->second);
      }
    }
  } else {
    const SF::vector<mesh_int_t>& layout = mesh.pl.algebraic_layout();
    const mesh_int_t start = layout[mpi_rank];
    const mesh_int_t stop  = layout[mpi_rank + 1];
 
    SF::vector<mesh_int_t> recv_nodes;
    size_t buffsize = 0;
    for(int pid = 0; pid < mpi_size; pid++) {
      if(mpi_rank == pid) {
        recv_nodes = selected_nodes;
        buffsize = recv_nodes.size();
      }
 
      MPI_Bcast(&buffsize, sizeof(size_t), MPI_BYTE, pid, mesh.comm);
      recv_nodes.resize(buffsize);
      MPI_Bcast(recv_nodes.data(), buffsize * sizeof(mesh_int_t), MPI_BYTE, pid, mesh.comm);
 
      for(mesh_int_t node : recv_nodes) {
        if(node >= start && node < stop)
          output_idx.push_back(node - start);
      }
    }
  }
 
  binary_sort(output_idx);
  unique_resize(output_idx);
}
 
bool mesh_has_gridout_tags(const sf_mesh& mesh,
                           const hashmap::unordered_set<int>& output_tags)
{
  bool local_has_tag = false;
  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
    if(output_tags.count(mesh.tag[eidx]) != 0) {
      local_has_tag = true;
      break;
    }
  }
 
  int local = local_has_tag ? 1 : 0;
  int global = 0;
  MPI_Allreduce(&local, &global, 1, MPI_INT, MPI_MAX, mesh.comm);
  return global != 0;
}
 
void extract_gridout_tag_mesh(const sf_mesh& mesh,
                              const hashmap::unordered_set<int>& output_tags,
                              sf_mesh& out_mesh)
{
  SF::vector<bool> keep_elem(mesh.l_numelem, false);
  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++)
    keep_elem[eidx] = output_tags.count(mesh.tag[eidx]) != 0;
 
  out_mesh.name = mesh.name;
  extract_mesh(keep_elem, mesh, out_mesh);
  if(out_mesh.g_numelem == 0) {
    log_msg(0, 5, ECHO, "gridout_tags selected no elements in \"%s\".", mesh.name.c_str());
    EXIT(EXIT_FAILURE);
  }
 
  rebalance_mesh(out_mesh);
 
  SF::vector<mesh_real_t> pts(mesh.xyz);
  SF::vector<mesh_int_t> ptsidx(mesh.get_numbering(SF::NBR_REF));
  std::list<sf_mesh*> meshlist;
  meshlist.push_back(&out_mesh);
  insert_points(pts, ptsidx, meshlist);
 
  SF::submesh_numbering<mesh_int_t, mesh_real_t> numbering;
  numbering(out_mesh);
  out_mesh.generate_par_layout();
}
 
void extract_gridout_surface_mesh(const sf_mesh& mesh, sf_mesh& out_mesh)
{
  compute_surface_mesh(mesh, SF::NBR_REF, out_mesh);
  out_mesh.name = mesh.name;
 
  // the surface connectivity refers to the nodes of the parent mesh. Turn it into a
  // node numbering of its own, so that the surface has points independent of the parent.
  out_mesh.localize(SF::NBR_REF);
 
  // surface elements are created here and carry no reference numbering of their own
  SF::vector<mesh_int_t>& ref_eidx = out_mesh.register_numbering(SF::NBR_ELEM_REF);
  ref_eidx.resize(out_mesh.l_numelem);
 
  long int num_local = out_mesh.l_numelem, offset = 0;
  MPI_Exscan(&num_local, &offset, 1, MPI_LONG, MPI_SUM, out_mesh.comm);
  if(get_rank() == 0) offset = 0;
  for(size_t eidx = 0; eidx < out_mesh.l_numelem; eidx++)
    ref_eidx[eidx] = offset + eidx;
 
  SF::vector<mesh_real_t> pts(mesh.xyz);
  SF::vector<mesh_int_t> ptsidx(mesh.get_numbering(SF::NBR_REF));
  std::list<sf_mesh*> meshlist;
  meshlist.push_back(&out_mesh);
  insert_points(pts, ptsidx, meshlist);
 
  SF::submesh_numbering<mesh_int_t, mesh_real_t> numbering;
  numbering(out_mesh);
  out_mesh.generate_par_layout();
}
 
void write_gridout_surface(const sf_mesh& mesh, const std::string& basename)
{
  sf_mesh surf;
  extract_gridout_surface_mesh(mesh, surf);
 
  write_points_parallel(surf, false, basename);
 
  sf_mesh surf_global = surf;
  const SF::vector<mesh_int_t>& nbr = surf.get_numbering(SF::NBR_SUBMESH);
  for(size_t i = 0; i < surf_global.con.size(); i++)
    surf_global.con[i] = nbr[surf_global.con[i]];
 
  write_surface(surf_global, basename + ".surf");
}
 
} // namespace
 
void parse_params_cpy(int argc, char** argv)
{
  LegacyCompareInput        legacy_input;
  std::vector<std::string> normalized_args;
  RuntimeCompatOptions     compat;
  std::string              normalize_error;
  if (!build_legacy_compare_input(argc, argv, &legacy_input, &normalize_error) ||
      !normalize_runtime_args(argc, argv, &normalized_args, &compat, &normalize_error)) {
    fprintf(stderr, "\n*** %s\n\n", normalize_error.c_str());
    exit(EXIT_FAILURE);
  }
 
  paramschema::ExecutionResult executed;
  if (!apply_parser_runtime_args(normalized_args, true, &executed)) {
    exit(EXIT_FAILURE);
  }
 
  if (legacy_input.available) {
    for (std::size_t i = 0; i < executed.validation.assignments.size(); ++i) {
      if (!executed.validation.assignments[i].synthesized) {
        continue;
      }
      legacy_input.available = false;
      legacy_input.runtime_args.clear();
      legacy_input.runtime_args.push_back(normalized_args[0]);
      legacy_input.unavailable_reason =
          "legacy compare is unavailable because the parser inferred optional controller counts from the original input";
      break;
    }
  }
 
  if (legacy_input.available && !executed.validation.legacy_compare_incompatible_paths.empty()) {
    legacy_input.available = false;
    legacy_input.warn_when_unavailable = false;
    legacy_input.runtime_args.clear();
    legacy_input.runtime_args.push_back(normalized_args[0]);
    legacy_input.unavailable_reason =
        "legacy compare is unavailable because the original input uses aggregate ID syntax without a legacy baseline for " +
        join_paths_for_message(executed.validation.legacy_compare_incompatible_paths);
  }
 
  if (!run_parser_legacy_compare(normalized_args, legacy_input, compat)) {
    exit(EXIT_FAILURE);
  }
}
 
 
void register_physics()
{
  // here all the physics can be registered to the physics registry
  // then they should be processed automatically
  if(phys_defined(PHYSREG_LAPLACE))
    user_globals::physics_reg[laplace_phys] = new Laplace();
  if(phys_defined(PHYSREG_EMI)) {
#if WITH_EMI_MODEL
    user_globals::physics_reg[emi_phys] = new EMI();
#else
    log_msg(NULL, 5, ECHO, "The EMI model was not compiled for this binary.\n");
    exit(EXIT_FAILURE);
#endif
  }
 
  if (phys_defined(PHYSREG_EIKONAL))
    user_globals::physics_reg[eikonal_phys] = new Eikonal();
 
  if(phys_defined(PHYSREG_INTRA_ELEC))
    user_globals::physics_reg[elec_phys] = new Electrics();
}
 
void initialize_physics()
{
  log_msg(0,0,0, "\n    *** Initializing physics ***\n");
 
  //load in the external imp modules
#ifdef HAVE_DLOPEN
  for (int ii = 0; ii < param_globals::num_external_imp; ii++) {
    int loading_succeeded = limpet::load_ionic_module(param_globals::external_imp[ii]);
    assert(loading_succeeded);
  }
#else
  if(param_globals::num_external_imp)
    log_msg(NULL, 4, ECHO,"Loading of external LIMPET modules not enabled.\n"
                          "Recompile with DLOPEN set.\n" );
#endif
 
  // init physics via Basic_physic interface
  for(auto it : user_globals::physics_reg) {
    Basic_physic* p = it.second;
    log_msg(NULL, 0, 0, "Initializing %s ..", p->name);
    p->initialize();
  }
}
 
void destroy_physics()
{
  log_msg(0,0,0, "\n    *** Destroying physics ***\n");
 
  for(auto it : user_globals::physics_reg) {
    Basic_physic* p = it.second;
    log_msg(NULL, 0, 0, "Destroying %s ..", p->name);
    p->destroy();
  }
}
 
// ignore_extracellular stim must be moved to stimulate.cc to be able
// to use all defined set operations instead of defines
 
void ignore_extracellular_stim(Stimulus *st, int ns, int ignore)
{
  // needs to be switch to stim enum types defined in stimulate.h
  for ( int i=0; i<ns; i++ ) {
    int turn_off = 0;
    turn_off += (st[i].stimtype == Extracellular_Ground) && (ignore & NO_EXTRA_GND);
    turn_off += (IsExtraV(st[i])) && (ignore & NO_EXTRA_V);
    turn_off += (st[i].stimtype==Extracellular_I) && (ignore & NO_EXTRA_I);
 
    if (turn_off)  {
      st[i].stimtype = Ignore_Stim;
      log_msg( NULL, 1, 0, "Extracellular stimulus %d ignored for monodomain", i );
    } else if ( st[i].stimtype==Intracellular_I ) {
      st[i].stimtype = Transmembrane_I;
      log_msg( NULL, 1, 0, "Intracellular stimulus %d converted to transmembrane", i );
    }
  }
}
 
int set_ignore_flags( int mode )
{
  if(mode==MONODOMAIN)
    return STM_IGNORE_MONODOMAIN;
  if(mode==BIDOMAIN)
    return STM_IGNORE_BIDOMAIN;
  if(mode==PSEUDO_BIDM)
    return STM_IGNORE_PSEUDO_BIDM;
 
  return IGNORE_NONE;
}
 
 
void check_nullspace_ok()
{
  Stimulus* s = param_globals::stimulus;
 
  for(int i=0; i < param_globals::num_stim; i++) {
    if(s[i].stimtype == Extracellular_Ground ||
       s[i].stimtype == Extracellular_V      ||
       s[i].stimtype == Extracellular_V_OL)
      return;
  }
 
  // for now we only warn, although we should actually stop the run
  log_msg( NULL, 4, 0,"Elliptic system is singular!\n"
                      "Use an explicit ground:voltage (stimulus[X].stimtype=3)\n"
                      "Do not trust the elliptic solution of this simulation run!\n");
}
 
void show_build_info()
{
  const char* real_precision = OPENCARP_REAL_BITS == 32 ? "single" : "double";
  printf("\n""*** GIT tag:            %s\n",   GIT_COMMIT_TAG);
  printf(    "*** GIT hash:           %s\n",   GIT_COMMIT_HASH);
  printf(    "*** GIT repo:           %s\n",   GIT_PATH);
  printf(    "*** local index bits:   %d\n",   OPENCARP_LOCAL_INDEX_BITS);
  printf(    "*** global index bits:  %d\n",   OPENCARP_GLOBAL_INDEX_BITS);
  printf(    "*** real precision:     %s (%d-bit)\n", real_precision, OPENCARP_REAL_BITS);
  printf(    "*** dependency commits: %s\n\n", SUBREPO_COMMITS);
}
 
void check_and_convert_params()
{
  // convert time steps to milliseconds
  param_globals::dt      /= 1000.;
 
  // check parab-solve solution method
  // explicit integration inverts the mass matrix directly, which is only feasible if it is lumped
  if(param_globals::mass_lumping == 0 && param_globals::parab_solve==0)  {
    log_msg(0,5,0, "parab_solve = 0 (explicit) requires mass_lumping = 1. "
                   "Either enable mass lumping or choose an implicit parab_solve method.");
    EXIT(1);
  }
 
  // check if we have to modify stimuli based on used bidomain setting
  if(!param_globals::extracell_monodomain_stim)
    ignore_extracellular_stim(param_globals::stimulus, param_globals::num_stim,
            set_ignore_flags(param_globals::bidomain));
 
  // check nullspace if necessary
  // if((param_globals::bidomain==BIDOMAIN) ||
  //    (param_globals::bidomain==PSEUDO_BIDM))
  //   check_nullspace_ok();
 
  if(param_globals::t_sentinel > 0 && param_globals::sentinel_ID < 0 ) {
    log_msg(0,4,0, "Warning: t_sentinel is set but no sentinel_ID has been specified; check_quiescence() behavior may not be as expected");
  }
 
  if(param_globals::num_external_imp > 0 ) {
    for(int ext_imp_i = 0; ext_imp_i < param_globals::num_external_imp; ext_imp_i++) {
      if(param_globals::external_imp[ext_imp_i][0] != '/') {
        log_msg(0,5,0, "external_imp[%d] error: absolute paths must be used for .so file loading (\'%s\')",
                ext_imp_i, param_globals::external_imp[ext_imp_i]);
        EXIT(1);
      }
    }
  }
 
  if(param_globals::experiment == EXP_LAPLACE && param_globals::bidomain != 1) {
    log_msg(0,4,0, "Warning: Laplace experiment mode requires bidomain = 1. Setting bidomain = 1.");
    param_globals::bidomain = 1;
  }
 
  if(param_globals::num_phys_regions == 0) {
    log_msg(0,4,0, "Warning: No physics region defined! Please set phys_region parameters to correctly define physics.");
 
    if(param_globals::experiment != EXP_LAPLACE) {
      log_msg(0,4,0, "Intra-elec and Extra-elec domains will be derived from fibers.\n");
      param_globals::num_phys_regions = param_globals::bidomain ? 2 : 1;
      param_globals::phys_region = (p_region*) calloc(param_globals::num_phys_regions, sizeof(p_region));
      param_globals::phys_region[0].ptype  = PHYSREG_INTRA_ELEC;
      param_globals::phys_region[0].name   = strdup("Autogenerated intracellular Electrics");
      param_globals::phys_region[0].num_IDs = 0;
 
      if(param_globals::bidomain) {
        param_globals::phys_region[1].ptype   = PHYSREG_EXTRA_ELEC;
        param_globals::phys_region[1].name    = strdup("Autogenerated extracellular Electrics");
        param_globals::phys_region[1].num_IDs = 0;
      }
    } else {
      log_msg(0,4,0, "Laplace domain will be derived from fibers.\n");
      param_globals::num_phys_regions = 1;
      param_globals::phys_region = (p_region*) calloc(param_globals::num_phys_regions, sizeof(p_region));
      param_globals::phys_region[0].ptype  = PHYSREG_LAPLACE;
      param_globals::phys_region[0].name   = strdup("Autogenerated Laplace");
      param_globals::phys_region[0].num_IDs = 0;
    }
  }
 
  if(param_globals::experiment == EXP_LAPLACE && !phys_defined(PHYSREG_LAPLACE)) {
    log_msg(0,4,0, "Warning: Laplace experiment mode requires a laplace physics regions defined.");
 
    int idx = -1;
    if((idx = get_phys_index(PHYSREG_EXTRA_ELEC)) > -1) {
      log_msg(0,4,0, "Converting the defined extracellular-electrics-region to laplace-region.");
      param_globals::phys_region[idx].ptype = PHYSREG_LAPLACE;
    } else if ((idx = get_phys_index(PHYSREG_INTRA_ELEC)) > -1) {
      log_msg(0,4,0, "Converting the defined intracellular-electrics-region to laplace-region.");
      param_globals::phys_region[idx].ptype = PHYSREG_LAPLACE;
    } else {
      param_globals::num_phys_regions += 1;
      param_globals::phys_region = (p_region*) realloc(param_globals::phys_region, param_globals::num_phys_regions * sizeof(p_region));
 
      param_globals::phys_region[param_globals::num_phys_regions - 1].ptype  = PHYSREG_LAPLACE;
      param_globals::phys_region[param_globals::num_phys_regions - 1].name   = strdup("Autogenerated Laplace");
      param_globals::phys_region[param_globals::num_phys_regions - 1].num_IDs = 0;
    }
  }
 
#ifndef WITH_PARMETIS
  if(param_globals::pstrat == 1) {
    log_msg(0,3,0, "openCARP was built without Parmetis support. Swithing to KDtree.");
    param_globals::pstrat = 2;
  }
#endif
 
  // check if we have the legacy stimuli or the new stimuli defined by the user
  bool legacy_stim_set = false, new_stim_set = false;
 
  for(int i=0; i<param_globals::num_stim; i++) {
    Stimulus & legacy_stim = param_globals::stimulus[i];
    Stim & new_stim = param_globals::stim[i];
 
    if(legacy_stim.stimtype || legacy_stim.strength)
      legacy_stim_set = true;
 
    if(new_stim.crct.type || new_stim.pulse.strength)
      new_stim_set = true;
  }
 
  if(legacy_stim_set || new_stim_set) {
    if(legacy_stim_set && new_stim_set) {
      log_msg(0,4,0, "Warning: Legacy stimuli and default stimuli are defined. Only default stimuli will be used!");
    }
    else if (legacy_stim_set) {
      log_msg(0,1,0, "Warning: Legacy stimuli defined. Please consider switching to stimulus definition \"stim[]\"!");
      user_globals::using_legacy_stimuli = true;
    }
  }
  else {
    log_msg(0,4,0, "Warning: No potential or current stimuli found!");
  }
}
 
void log_runtime_apl_notice()
{
  if(mpi_runtime_ready() && get_rank() != 0)
    return;
 
  log_msg(0, 0, 0, "");
  log_msg(0, 0, 0, "LICENSE NOTICE:");
  log_msg(0, 0, 0,
          "This run uses openCARP components licensed under the Academic Public License v1.1.");
  log_msg(0, 0, 0,
          "That includes code from simulator/**, param/**, and APL-licensed parts of physics/**.");
  log_msg(0, 0, FLUSH,
          "Commercial use is not allowed without a separate commercial license.");
}
 
void set_io_dirs(char *sim_ID, char *pp_ID, IO_t init)
{
  int flg = 0, err = 0, rank = get_rank();
 
  char *ptr = getcwd(current_dir, 1024);
  if (ptr == NULL) err++;
  ptr = getcwd(input_dir, 1024);
  if (ptr == NULL) err++;
  //if (param_globals::experiment == 4 && post_processing_opts == MECHANIC_POSTPROCESS)
  //  {  sim_ID = param_globals::ppID; param_globals::ppID = "POSTPROC_DIR"; }
 
  // output directory
  if (rank == 0)  {
    if (strcmp(sim_ID, "OUTPUT_DIR")) {
      if (mkdir(sim_ID, 0775)) {             // rwxrwxr-x
        if (errno == EEXIST )  {
          log_msg(NULL, 2, 0, "Output directory exists: %s\n", sim_ID);
        } else  {
          log_msg(NULL, 5, 0, "Unable to make output directory\n");
          flg = 1;
        }
      }
    } else if (mkdir(sim_ID, 0775) &&  errno != EEXIST) {
      log_msg(NULL, 5, 0, "Unable to make output directory\n");
      flg = 1;
    }
  }
 
  // terminate?
  if(get_global(flg, MPI_SUM)) { EXIT(-1); }
 
  err += chdir(sim_ID);
  ptr = getcwd(output_dir, 1024);
  if (ptr == NULL) err++;
 
  // terminate?
  if(get_global(err, MPI_SUM)) { EXIT(-1); }
 
  err += chdir(output_dir);
 
  // postprocessing directory
  if (rank == 0 && ((param_globals::experiment==EXP_POSTPROCESS) || (param_globals::post_processing_opts & LEADFIELD)))  {
 
    if (strcmp(param_globals::ppID, "POSTPROC_DIR")) {
      if (mkdir(param_globals::ppID, 0775)) {             // rwxrwxr-x
        if (errno == EEXIST )  {
          log_msg(NULL, 2, ECHO, "Postprocessing directory exists: %s\n\n", param_globals::ppID);
        } else  {
          log_msg(NULL, 5, ECHO, "Unable to make postprocessing directory\n\n");
          flg = 1;
        }
      }
    } else if (mkdir(param_globals::ppID, 0775) &&  errno != EEXIST) {
      log_msg(NULL, 5, ECHO, "Unable to make postprocessing directory\n\n");
      flg = 1;
    }
 
  }
 
  if(get_global(flg, MPI_SUM)) { EXIT(-1); }
 
  err += chdir(param_globals::ppID);
  ptr = getcwd(postproc_dir, 1024);
  if (ptr == NULL) err++;
  err = chdir(output_dir);
  if(get_global(err, MPI_SUM)) { EXIT(-1); }
 
  err = set_dir(init);
  if(get_global(err, MPI_SUM)) { EXIT(-1); }
}
 
bool setup_IO(int argc, char **argv)
{
  bool io_node = false;
  int psize = get_size(), prank = get_rank();
 
  if (param_globals::num_io_nodes > 0) {
    // Can't do async IO with only one core
    if (get_size() == 1) {
      log_msg(NULL, 5, 0, "You cannot run with async IO on only one core.\n");
      EXIT(EXIT_FAILURE);
    }
    // Can't do async IO with more IO cores than compute cores
    if (2 * param_globals::num_io_nodes >= psize) {
      log_msg(NULL, 5, 0, "The number of IO cores be less " "than the number of compute cores.");
      EXIT(EXIT_FAILURE);
    }
#if 0
    if (param_globals::num_PS_nodes && param_globals::num_io_nodes > param_globals::num_PS_nodes) {
      LOG_MSG(NULL, 5, 0,
              "The number of IO cores (%d) should not "
              "exceed the number of PS compute cores (%d).\n",
               param_globals::num_io_nodes, param_globals::num_PS_nodes);
      EXIT(-1);
    }
#endif
    // root IO node is global node 0
    io_node = prank < param_globals::num_io_nodes;
 
    MPI_Comm comm;
    MPI_Comm_split(SF_COMM, io_node, get_rank(), &comm);
    MPI_Comm_set_name(comm, io_node ? "IO" : "compute");
 
    SF_COMM = comm;    // either the compute world or IO world
#ifdef WITH_PETSC
    if(petsc_active()) petsc_set_comm_world(comm);
#endif
 
    prank = get_rank();
 
    MPI_Intercomm_create(comm, 0, MPI_COMM_WORLD, io_node ? param_globals::num_io_nodes : 0,
                         6667, &user_globals::IO_Intercomm);
 
    if(prank != get_rank(user_globals::IO_Intercomm))
      log_msg(NULL, 4, 0, "Global node %d, Comm rank %d != Intercomm rank %d\n",
              get_rank(MPI_COMM_WORLD), get_rank(SF_COMM),
              get_rank(user_globals::IO_Intercomm));
  } else
    MPI_Comm_set_name(SF_COMM, "compute");
 
  set_io_dirs(param_globals::simID, param_globals::ppID, OUTPUT);
 
  int flg = 0;
  if((io_node || !param_globals::num_io_nodes) && !prank)
    flg = !output_parameter_file("parameters.par", argc, argv);
  if(get_global(flg, MPI_SUM)) EXIT(EXIT_FAILURE);
 
  return io_node;
}
void update_cwd()
{
  getcwd(current_dir, 1024);
}
 
int set_dir(IO_t dest)
{
  int err;
 
  if      (dest==OUTPUT)   err = chdir(output_dir);
  else if (dest==POSTPROC) err = chdir(postproc_dir);
  else if (dest==CURDIR)   err = chdir(current_dir);
  else                     err = chdir(input_dir);
 
  return err;
}
 
void basic_timer_setup()
{
  // if we restart from a checkpoint, the timer_manager will be notified at a later stage
  double start_time = 0.0;
  user_globals::tm_manager = new timer_manager(param_globals::dt, start_time, param_globals::tend);
 
  double end_time    = param_globals::tend;
  timer_manager & tm = *user_globals::tm_manager;
 
  if(param_globals::experiment == EXP_LAPLACE) {
    tm.initialize_singlestep_timer(tm.time, 0, iotm_console,   "IO (console)",  nullptr);
    tm.initialize_singlestep_timer(tm.time, 0, iotm_state_var, "IO (state vars)", nullptr);
    tm.initialize_singlestep_timer(tm.time, 0, iotm_spacedt,   "IO (spacedt)",  nullptr);
  }
  else {
    tm.initialize_eq_timer(tm.time, end_time, 0, param_globals::timedt,  0, iotm_console,   "IO (console)");
    tm.initialize_eq_timer(tm.time, end_time, 0, param_globals::spacedt, 0, iotm_state_var, "IO (state vars)");
    tm.initialize_eq_timer(tm.time, end_time, 0, param_globals::spacedt, 0, iotm_spacedt,   "IO (spacedt)");
  }
 
  if(param_globals::num_tsav) {
    std::vector<double> trig(param_globals::num_tsav);
    for(size_t i=0; i<trig.size(); i++) trig[i] = param_globals::tsav[i];
 
    tm.initialize_neq_timer(trig, 0, iotm_chkpt_list, "instance checkpointing");
  }
 
  if(param_globals::chkpt_intv)
    tm.initialize_eq_timer(param_globals::chkpt_start, param_globals::chkpt_stop, 0,
                           param_globals::chkpt_intv, 0, iotm_chkpt_intv, "interval checkpointing");
 
  if(param_globals::num_trace)
    tm.initialize_eq_timer(tm.time, end_time, 0, param_globals::tracedt, 0, iotm_trace, "IO (node trace)");
}
 
#ifdef WITH_POWERCAPPING
void basic_powercapping_setup()
{
  user_globals::pc_manager = new powercapping_manager();
}
 
void basic_powercapping_cleanup()
{
  /* Must be called explicitly to ensure the orderly termination of potential
   * background daemons */
  delete user_globals::pc_manager; user_globals::pc_manager = nullptr;
}
#endif
 
void get_protocol_column_widths(std::vector<int> & col_width, std::vector<int> & used_timer_ids)
{
  char buff[256];
  const short padding = 4;
  Electrics* elec = (Electrics*) get_physics(elec_phys, false);
 
  do {
    snprintf(buff, sizeof buff, "%.3lf", user_globals::tm_manager->time);
    if(col_width[0] < int(strlen(buff)+padding))
      col_width[0] = strlen(buff)+padding;
 
    snprintf(buff, sizeof buff, "%.3ld", user_globals::tm_manager->d_time);
    if(col_width[1] < int(strlen(buff)+padding))
      col_width[1] = strlen(buff)+padding;
 
    int col = 2;
    for (size_t tid = 0; tid < used_timer_ids.size(); tid++)
    {
      int timer_id  = used_timer_ids[tid];
      base_timer* t = user_globals::tm_manager->timers[timer_id];
 
      if(t->d_trigger_dur && elec) {
        // figure out value of signal linked to this timer
        double val = 0.;
 
        // determine timer linked to which physics, for now we deal with electrics only
        val = elec->timer_val(timer_id);
 
        snprintf(buff, sizeof buff, "%.3lf", val);
        if(col_width[col] < int(strlen(buff)+padding))
          col_width[col] = strlen(buff)+padding;
      }
      col++;
    }
 
    // advance time
    user_globals::tm_manager->update_timers();
  } while (!user_globals::tm_manager->elapsed());
 
  user_globals::tm_manager->reset_timers();
}
int plot_protocols(const char *fname)
{
  int   err = {0};
  std::ofstream fh;
  const char* smpl_endl = "\n";
 
  if(!get_rank()) {
    fh.open(fname);
 
    // If we couldn't open the output file stream for writing
    if (!fh) {
        // Print an error and exit
        log_msg(0,5,0,"Protocol file %s could not be opened for writing!\n", fname);
        err = -1;
    }
  }
 
  // broadcast and return if err
  if(get_global(err, MPI_SUM))
    return err;
 
  // only rank 0 writes
  if(!get_rank()) {
 
    // collect timer information, label, short label, unit
    std::vector<std::string> col_labels       = {"time", "tick"};
    std::vector<std::string> col_short_labels = {"A",  "B"};
    std::vector<std::string> col_unit_labels  = {"ms", "--" };
    std::vector<int> col_width = {4, 4};
 
    char c_label = {'C'};
    std::string label = {""};
    std::string unit  = {""};
 
    // here we store the IDs of the timers that we care about. currently this are the IO and TS timers
    // and the electricts timers
    std::vector<int> used_timer_ids;
    std::vector<int> used_stim_ids;
 
    Electrics* elec = (Electrics*) get_physics(elec_phys, false);
    if(elec) {
      int sidx = 0;
      for(const stimulus & s : elec->stimuli) {
        if(s.ptcl.timer_id > -1) {
          base_timer* t = user_globals::tm_manager->timers[s.ptcl.timer_id];
          if(t) {
            used_timer_ids.push_back(s.ptcl.timer_id);
            used_stim_ids.push_back(sidx);
          }
        }
 
        sidx++;
      }
    }
 
    // determine longest timer label
    int mx_llen = 0;
    for (size_t tid = 0; tid < used_timer_ids.size(); tid++)
    {
      int timer_id  = used_timer_ids[tid];
      base_timer* t = user_globals::tm_manager->timers[timer_id];
 
      int llen = strlen(t->name);
      mx_llen = llen > mx_llen ? llen : mx_llen;
    }
 
    for (size_t tid = 0; tid < used_timer_ids.size(); tid++)
    {
      int timer_id  = used_timer_ids[tid];
      base_timer* t = user_globals::tm_manager->timers[timer_id];
 
      col_labels.push_back(t->name);
      label = c_label;
      col_short_labels.push_back(label);
 
      if(elec) {
        // search physics for signals linked to timer
        unit = elec->timer_unit(timer_id);
        if(unit.empty()) unit = "--";
        col_unit_labels.push_back(unit);
        col_width.push_back(4);
      }
      c_label++;
    }
 
    get_protocol_column_widths(col_width, used_timer_ids);
 
    // print header + legend first
    fh << "# Protocol header\n#\n" << "# Legend:\n";
    for(size_t i = 0; i<col_short_labels.size(); i++)
    {
      fh << "#" << std::setw(2) << col_short_labels[i] << " = " << std::setw(mx_llen) << col_labels[i];
      fh << " [" << std::setw(10) << col_unit_labels[i] << "]";
 
      if(i >= 2 && used_stim_ids[i-2] > -1) {
        stimulus & s = elec->stimuli[used_stim_ids[i-2]];
 
        if (is_potential(s.phys.type)) {
          if(s.phys.type == GND_ex)
            fh << " ground stim" << smpl_endl;
          else
            fh << " applied: " << std::to_string(s.pulse.strength) << smpl_endl;
        } else {
          fh << smpl_endl;
        }
      } else {
        fh << smpl_endl;
      }
    }
 
    // plot column short labels
    fh << "#";
    for(size_t i = 0; i<col_short_labels.size(); i++)
      fh << std::setw(col_width[i] - 3) <<  col_short_labels[i].c_str() << std::setw(3) << " ";
 
    // plot column units
    fh << smpl_endl << "#";
    for(size_t i = 0; i<col_unit_labels.size(); i++)
      fh << "[" << std::setw(col_width[i]-2) <<  col_unit_labels[i].c_str() << "]";
 
    // step through simulated time period
    fh << smpl_endl << std::fixed;
    do {
      // time and discrete time
      fh << std::setw(col_width[0]) << std::setprecision(3) << user_globals::tm_manager->time;
      fh << std::setw(col_width[1]) << user_globals::tm_manager->d_time;
 
      // iterate over all timers
      int col = 2;
      for (size_t tid = 0; tid < used_timer_ids.size(); tid++)
      {
        int timer_id  = used_timer_ids[tid];
        base_timer* t = user_globals::tm_manager->timers[timer_id];
 
        // type of timer: plain trigger or trigger linked to signal
        if(!t->d_trigger_dur) {
          int On = t->triggered ? 1 : 0;
          fh << std::setw(col_width[col]) << On;
        } else if(elec) {
          // figure out value of signal linked to this timer
          double val = 0.;
 
          // determine timer linked to which physics, for now we deal with electrics only
          val = elec->timer_val(timer_id);
 
          fh << std::setw(col_width[col]) << std::setprecision(3) << val;
        }
        col++;
      }
 
      fh << smpl_endl;
 
      // advance time
      user_globals::tm_manager->update_timers();
    } while (!user_globals::tm_manager->elapsed());
 
    fh.close();
 
    // reset timer to start before actual simulation
    user_globals::tm_manager->reset_timers();
  }
 
  return err;
}
 
void init_console_output(const timer_manager & tm, prog_stats & p)
{
  const char*  h1_prog = "PROG\t----- \t----\t-------\t-------|";
  const char*  h2_prog = "time\t%%comp\ttime\t ctime \t  ETA  |";
  const char*  h1_wc   = "\tELAPS |";
  const char*  h2_wc   = "\twc    |";
 
  p.start = get_time();
  p.last  = p.start;
 
  log_msg(NULL, 0, 0, "%s", h1_prog );
  log_msg(NULL, 0, NONL, "%s", h2_prog );
  log_msg(NULL, 0, 0, "" );
}
 
 
void time_to_string(float time, char* str, short str_size)
{
  int req_hours =  ((int)(time)) / 3600;
  int req_min   = (((int)(time)) % 3600) / 60;
  int req_sec   = (((int)(time)) % 3600) % 60;
 
  snprintf(str, str_size, "%d:%02d:%02d", req_hours, req_min, req_sec);
}
 
void update_console_output(const timer_manager & tm, prog_stats & p)
{
 
  float progress = 100.*(tm.time - tm.start) / (tm.end - tm.start);
  float elapsed_time = timing(p.curr, p.start);
  float req_time = (elapsed_time / progress) * (100.0f - progress);
 
  if(progress == 0.0f)
    req_time = 0.0f;
 
  char elapsed_time_str[256];
  char req_time_str[256];
  time_to_string(elapsed_time, elapsed_time_str, 255);
  time_to_string(req_time,     req_time_str,     255);
 
  log_msg( NULL, 0, NONL, "%.2f\t%.1f\t%.1f\t%s\t%s",
      tm.time,
      progress,
      (float)(p.curr - p.last),
      elapsed_time_str,
      req_time_str);
 
  p.last = p.curr;
 
  // we add an empty string for newline and flush
  log_msg( NULL, 0, ECHO | FLUSH, "");
}
 
void simulate()
{
  // here we want to include all time dependent physics, to check if we have any of those
  bool have_timedependent_phys = (phys_defined(PHYSREG_INTRA_ELEC) || phys_defined(PHYSREG_EIKONAL) || phys_defined(PHYSREG_EMI));
 
  if(!have_timedependent_phys) {
    log_msg(0,0,0, "\n no time-dependent physics region registered, skipping simulate loop..\n");
    return;
  }
 
  log_msg(0,0,0, "\n    *** Launching simulation ***\n");
 
  set_dir(OUTPUT);
 
  if(param_globals::dump_protocol)
    plot_protocols("protocol.trc");
 
  prog_stats prog;
  timer_manager & tm = *user_globals::tm_manager;
  init_console_output(tm, prog);
 
#ifdef WITH_POWERCAPPING
  basic_powercapping_setup();
  powercapping_manager *pc = user_globals::pc_manager;
#endif
 
  // main loop
  do {
    // console output
    if(tm.trigger(iotm_console)) {
      // print console
      update_console_output(tm, prog);
    }
 
#ifdef WITH_POWERCAPPING
    std::vector<int> flops_per_rank;
 
    // TODO: fill vector with expected flops per MPI rank for upcoming iteration
#endif
 
#ifdef WITH_POWERCAPPING
    pc->iteration_begin(flops_per_rank);
#endif
 
    // in order to be closer to carpentry we first do output and then compute the solution
    // for the next time slice ..
    if (tm.trigger(iotm_spacedt)) {
      for(const auto & it : user_globals::physics_reg) {
        it.second->output_step();
      }
    }
#ifdef WITH_POWERCAPPING
    pc->sample("output step");
#endif
 
    // compute step
    for(const auto & it : user_globals::physics_reg) {
      Basic_physic* p = it.second;
      if (tm.trigger(p->timer_idx))
        p->compute_step();
#ifdef WITH_POWERCAPPING
      pc->sample(p->name);
#endif
    }
 
#ifdef WITH_POWERCAPPING
    pc->iteration_end(flops_per_rank);
#endif
 
    // advance time
    tm.update_timers();
  } while (!tm.elapsed());
 
  log_runtime_apl_notice();
  log_msg(0,0,0, "\n\nTimings of individual physics:");
  log_msg(0,0,0, "------------------------------\n");
 
  for(const auto & it : user_globals::physics_reg) {
    Basic_physic* p = it.second;
    p->output_timings();
  }
 
#ifdef WITH_POWERCAPPING
  pc = nullptr;
  basic_powercapping_cleanup();
#endif
}
 
void post_process()
{
  if(param_globals::post_processing_opts & RECOVER_PHIE)  {
    log_msg(NULL,0,ECHO,"\nPOSTPROCESSOR: Recovering Phie ...");
    log_msg(NULL,0,ECHO,  "----------------------------------\n");
 
    // do postprocessing
    int err = postproc_recover_phie();
 
    if(!err)  {
      log_msg(NULL,0,ECHO,"\n-----------------------------------------");
      log_msg(NULL,0,ECHO,  "POSTPROCESSOR: Successfully recoverd Phie.\n");
    }
  }
 
  if(param_globals::post_processing_opts & LEADFIELD)  {
#ifdef WITH_LEADFIELD
    log_msg(NULL,0,ECHO,"\nPOSTPROCESSOR: Computing leadfields ...");
    log_msg(NULL,0,ECHO,  "-------------------------------------\n");
 
    Electrics* elec = static_cast<Electrics*>(get_physics(elec_phys));
    if(!elec) {
      log_msg(NULL, 5, 0, "Error: Leadfield requires active EP physics. Aborting.");
      EXIT(1);
    }
    Leadfield leadfield;
    int err = leadfield.run(*elec);
 
    if(!err)  {
      log_msg(NULL,0,ECHO,"\n------------------------------------------");
      log_msg(NULL,0,ECHO,  "POSTPROCESSOR: Successfully computed leadfields.\n");
    }
#else
    log_msg(NULL, 5, 0, "Error: leadfield support not compiled in. Rebuild with -DENABLE_LEADFIELD=ON.");
    EXIT(1);
#endif
  }
}
 
Basic_physic* get_physics(physic_t p, bool error_if_missing)
{
  auto it = user_globals::physics_reg.find(p);
 
  if(it != user_globals::physics_reg.end()) {
    return it->second;
  } else {
    if(error_if_missing) {
      log_msg(0,5,0, "%s error: required physic is not active! Usually this is due to an inconsistent experiment configuration. Aborting!", __func__);
      EXIT(EXIT_FAILURE);
    }
 
    return NULL;
  }
}
 
sf_vec* get_data(datavec_t d)
{
  sf_vec* ret = NULL;
 
  if(user_globals::datavec_reg.count(d))
    ret = user_globals::datavec_reg[d];
 
  return ret;
}
 
void register_data(sf_vec* dat, datavec_t d)
{
  if(user_globals::datavec_reg.count(d) == 0) {
    user_globals::datavec_reg[d] = dat;
  }
  else {
    log_msg(0,5,0, "%s warning: trying to register already registered data vector.", __func__);
  }
}
 
void parse_mesh_types()
{
  std::map<mesh_t, sf_mesh> & mesh_registry = user_globals::mesh_reg;
 
  // This is the initial grid we read the hard-disk data into
  mesh_registry[reference_msh] = sf_mesh();
  // we specify the MPI communicator for the reference mesh,
  // all derived meshes will get this comminicator automatically
  mesh_registry[reference_msh].comm = SF_COMM;
 
  auto register_new_mesh = [&] (mesh_t mt, int pidx) {
    if(!mesh_registry.count(mt)) {
      mesh_registry[mt] = sf_mesh();
      mesh_registry[mt].name = param_globals::phys_region[pidx].name;
    }
    return &mesh_registry[mt];
  };
 
  // based on cli parameters we determine which grids need to be defined
  for(int i=0; i<param_globals::num_phys_regions; i++)
  {
    sf_mesh* curmesh = NULL;
    // register mesh type
    switch(param_globals::phys_region[i].ptype) {
      case PHYSREG_EIKONAL:
      case PHYSREG_INTRA_ELEC:
        curmesh = register_new_mesh(intra_elec_msh, i);
        break;
 
      case PHYSREG_LAPLACE:
      case PHYSREG_EXTRA_ELEC:
        curmesh = register_new_mesh(extra_elec_msh, i);
        break;
#if WITH_EMI_MODEL
      case PHYSREG_EMI:
        curmesh = register_new_mesh(emi_msh, i);
        break;
#endif
 
      default:
        log_msg(0,5,0, "Unsupported mesh type %d! Aborting!", param_globals::phys_region[i].ptype);
        EXIT(EXIT_FAILURE);
    }
 
    if(curmesh) {
      // set mesh unique tags
      for(int j=0; j<param_globals::phys_region[i].num_IDs; j++)
        curmesh->extr_tag.insert(param_globals::phys_region[i].ID[j]);
    }
  }
}
 
void retag_elements(sf_mesh & mesh, TagRegion *tagRegs, int ntr)
{
  if(ntr == 0) return;
  // checkTagRegDefs(ntr, tagRegs);
 
  SF::vector<mesh_int_t> & ref_eidx = mesh.get_numbering(SF::NBR_ELEM_REF);
 
  for (int i=0; i<ntr; i++) {
    tagreg_t type = tagreg_t(tagRegs[i].type);
    SF::vector<mesh_int_t> elem_indices;
 
    if (type == tagreg_list)
      read_indices(elem_indices, tagRegs[i].elemfile, ref_eidx, mesh.comm);
    else {
      geom_shape shape;
      shape.type   = geom_shape::shape_t(tagRegs[i].type);
      shape.p0.x   = tagRegs[i].p0[0];
      shape.p0.y   = tagRegs[i].p0[1];
      shape.p0.z   = tagRegs[i].p0[2];
      shape.p1.x   = tagRegs[i].p1[0];
      shape.p1.y   = tagRegs[i].p1[1];
      shape.p1.z   = tagRegs[i].p1[2];
      shape.radius = tagRegs[i].radius;
 
      bool nodal = false;
      indices_from_geom_shape(elem_indices, mesh, shape, nodal);
    }
 
    if(get_global((long int)elem_indices.size(), MPI_SUM, mesh.comm) == 0)
      log_msg(0,3,0,"Tag region %d is empty", i);
 
    for(size_t j=0; j<elem_indices.size(); j++)
      mesh.tag[elem_indices[j]] = tagRegs[i].tag;
  }
 
  // output the vector?
  if(strlen(param_globals::retagfile))
  {
    update_cwd();
    set_dir(OUTPUT);
 
    int dpn = 1;
    SF::write_data_ascii(mesh.comm, ref_eidx, mesh.tag, param_globals::retagfile, dpn);
 
    // Set dir back to what is was prior to retagfile output
    set_dir(CURDIR);
  }
}
 
size_t renormalise_fibres(SF::vector<mesh_real_t> &fib, size_t l_numelem)
{
  size_t renormalised_count = 0;
 
  // using pragma omp without global OMP control can lead to massive compute stalls,
  // as all cores may be already occupied by MPI, thus they become oversubscribed. Once
  // there is a global OMP control in place, we can activate this parallel for again. -Aurel, 20.01.2022
  // #pragma omp parallel for schedule(static) reduction(+ : renormalised_count)
  for (size_t i = 0; i < l_numelem; i++)
  {
    const mesh_real_t f0 = fib[3*i+0], f1 = fib[3*i+1], f2 = fib[3*i+2];
    mesh_real_t fibre_len = sqrt(f0*f0 + f1*f1 + f2*f2);
 
    if (fibre_len && fabs(fibre_len - 1) > 1e-3) {
      fib[3 * i + 0] /= fibre_len;
      fib[3 * i + 1] /= fibre_len;
      fib[3 * i + 2] /= fibre_len;
      renormalised_count++;
    }
  }
 
  return renormalised_count;
}
 
void setup_meshes(bool require_fibers=true)
{
  log_msg(0,0,0, "\n    *** Processing meshes ***\n");
 
  const std::string basename = param_globals::meshname;
  const int verb = param_globals::output_level;
  std::map<mesh_t, sf_mesh> & mesh_registry = user_globals::mesh_reg;
  assert(mesh_registry.count(reference_msh) == 1); // There must be a reference mesh
 
  set_dir(INPUT);
 
  // we always read into the reference mesh
  sf_mesh & ref_mesh = mesh_registry[reference_msh];
  MPI_Comm comm = ref_mesh.comm;
 
  int size, rank;
  double t1, t2, s1, s2;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
 
  SF::vector<mesh_real_t>    pts;
  SF::vector<mesh_int_t> ptsidx;
 
  // we add pointers to the meshes that need vertex cooridnates to this list
  std::list< sf_mesh* > ptsread_list;
 
  // read element mesh data
  t1 = MPI_Wtime(); s1 = t1;
  if(verb) log_msg(NULL, 0, 0,"Reading reference mesh: %s.*", basename.c_str());
 
  SF::read_elements(ref_mesh, basename, require_fibers);
  SF::read_points(basename, comm, pts, ptsidx);
 
  if (strlen(param_globals::tagfile)) {
    if(verb) log_msg(NULL, 0, 0, "Overriding element tags from: %s", param_globals::tagfile);
    SF::read_element_tags(ref_mesh, param_globals::tagfile);
  }
 
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  bool check_fibre_normality = true;
  if (check_fibre_normality and ref_mesh.fib.size()>0) {
    t1 = MPI_Wtime();
 
    // make sure that all fibre vectors have unit length
    size_t l_num_fixed_fib = renormalise_fibres(ref_mesh.fib, ref_mesh.l_numelem);
 
    size_t l_num_fixed_she = 0;
    if (ref_mesh.she.size() > 0)
      l_num_fixed_she = renormalise_fibres(ref_mesh.she, ref_mesh.l_numelem);
 
    unsigned long fixed[2] = {(unsigned long) l_num_fixed_fib, (unsigned long) l_num_fixed_she};
    MPI_Allreduce(MPI_IN_PLACE, fixed, 2, MPI_UNSIGNED_LONG, MPI_SUM, comm);
 
    if (fixed[0] + fixed[1] > 0)
      log_msg(NULL, 0, 0, "Renormalised %ld longitudinal and %ld sheet-transverse fibre vectors.", fixed[0], fixed[1]);
 
    t2 = MPI_Wtime();
    if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
  }
 
  if(param_globals::numtagreg > 0) {
    log_msg(0, 0, 0, "Re-tagging reference mesh");
 
    // the retagging requires vertex coordinates, as such we need to read them into
    // the reference mesh
    ptsread_list.push_back(&ref_mesh);
    SF::insert_points(pts, ptsidx, ptsread_list);
 
    retag_elements(ref_mesh, param_globals::tagreg, param_globals::numtagreg);
 
    // we clear the list of meshet to receive vertices
    ptsread_list.clear();
  }
 
  if(verb) log_msg(NULL, 0, 0, "Processing submeshes");
 
  bool have_emi_mesh = false;
  bool have_non_emi_mesh = false;
 
  for(auto it = mesh_registry.begin(); it != mesh_registry.end(); ++it) {
    mesh_t grid_type = it->first;
    sf_mesh & submesh = it->second;
 
    if(grid_type != reference_msh) {
      if(grid_type == emi_msh)
        have_emi_mesh = true;
      else
        have_non_emi_mesh = true;
 
      if(verb > 1) log_msg(NULL, 0, 0, "\nSubmesh name: %s", submesh.name.c_str());
      t1 = MPI_Wtime();
 
      if(submesh.extr_tag.size() && grid_type != emi_msh)
        extract_tagbased(ref_mesh, submesh);
      else {
        // all submeshes should be defined on sets of tags, for backwards compatibility
        // we do a fiber based intra_elec_msh extraction if no tags are provided. Also, we
        // could do special treatments of any other physics type here. It would defeat
        // the purpose of the tag-based design, though. -Aurel
        switch(grid_type) {
          case emi_msh:
          case intra_elec_msh: extract_myocardium(ref_mesh, submesh, require_fibers); break;
          default:        extract_tagbased(ref_mesh, submesh);   break;
        }
      }
      t2 = MPI_Wtime();
      if(verb > 1) log_msg(NULL, 0, 0, "Extraction done in %f sec.", float(t2 - t1));
 
      ptsread_list.push_back(&submesh);
    }
  }
 
  if(have_emi_mesh && have_non_emi_mesh) {
    log_msg(NULL, 5, ECHO, "EMI and non-EMI submeshes cannot be mixed during mesh setup.");
    EXIT(EXIT_FAILURE);
  }
 
  // KDtree partitioning requires the coordinates to be present in the mesh data
  if(param_globals::pstrat == 2 && have_non_emi_mesh)
    SF::insert_points(pts, ptsidx, ptsread_list);
 
  for(auto it = mesh_registry.begin(); it != mesh_registry.end(); ++it)
  {
    mesh_t grid_type = it->first;
    sf_mesh & submesh = it->second;
    if(grid_type != reference_msh && grid_type!=emi_msh) {
      if(verb > 2) log_msg(NULL, 0, 0, "\nSubmesh name: %s", submesh.name.c_str());
      SF::vector<mesh_int_t> part;
 
      // generate partitioning vector
      t1 = MPI_Wtime();
      switch(param_globals::pstrat) {
        case 0:
          if(verb > 2) log_msg(NULL, 0, 0, "Using linear partitioning ..");
          break;
 
#ifdef WITH_PARMETIS
        case 1:
        {
          if(verb > 2) log_msg(NULL, 0, 0, "Using Parmetis partitioner ..");
          SF::parmetis_partitioner<mesh_int_t, mesh_real_t> partitioner(param_globals::pstrat_imbalance, 2);
          partitioner(submesh, part);
          break;
        }
#endif
        default:
        case 2: {
          if(verb > 2) log_msg(NULL, 0, 0, "Using KDtree partitioner ..");
          SF::kdtree_partitioner<mesh_int_t, mesh_real_t> partitioner;
          partitioner(submesh, part);
          break;
        }
      }
      t2 = MPI_Wtime();
      if(verb > 2) log_msg(NULL, 0, 0, "Partitioning done in %f sec.", float(t2 - t1));
 
      if(param_globals::pstrat > 0) {
        if(param_globals::gridout_p) {
          std::string out_name = get_basename(param_globals::meshname);
          if(grid_type == intra_elec_msh)      out_name += "_i.part.dat";
          else if(grid_type == extra_elec_msh) out_name += "_e.part.dat";
 
          set_dir(OUTPUT);
          log_msg(0,0,0, "Writing \"%s\" partitioning to: %s", submesh.name.c_str(), out_name.c_str());
          write_data_ascii(submesh.comm, submesh.get_numbering(SF::NBR_ELEM_REF), part, out_name);
        }
 
        t1 = MPI_Wtime();
        SF::redistribute_elements(submesh, part);
        t2 = MPI_Wtime();
        if(verb > 2) log_msg(NULL, 0, 0, "Redistributing done in %f sec.", float(t2 - t1));
      }
 
      t1 = MPI_Wtime();
      SF::submesh_numbering<mesh_int_t, mesh_real_t> sm_numbering;
      sm_numbering(submesh);
      t2 = MPI_Wtime();
      if(verb > 2) log_msg(NULL, 0, 0, "Canonical numbering done in %f sec.", float(t2 - t1));
 
      t1 = MPI_Wtime();
      submesh.generate_par_layout();
      SF::petsc_numbering<mesh_int_t, mesh_real_t> p_numbering(submesh.pl, param_globals::renumbering);
      p_numbering(submesh);
      t2 = MPI_Wtime();
      if(verb > 2) log_msg(NULL, 0, 0, "PETSc numbering done in %f sec.", float(t2 - t1));
      if(verb > 2) print_DD_info(submesh);
    }
  }
 
  if(have_non_emi_mesh)
    SF::insert_points(pts, ptsidx, ptsread_list);
 
  ref_mesh.clear_data();
 
  s2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "All done in %f sec.", float(s2 - s1));
}
 
bool parse_gridout_tags(hashmap::unordered_set<int>& output_tags)
{
  static const char* parameter_name = "gridout_tags";
  const std::string spec = param_globals::gridout_tags ? param_globals::gridout_tags : "";
 
  std::vector<int> tags;
  std::string error;
  if(!paramschema::parse_idset_spec(spec, &tags, &error)) {
    log_msg(0, 5, ECHO, "Could not parse %s: %s.", parameter_name, error.c_str());
    EXIT(EXIT_FAILURE);
  }
 
  output_tags.clear();
  output_tags.insert(tags.begin(), tags.end());
 
  if(output_tags.size() == 0) return false;
 
  log_msg(0, 0, 0, "Restricting grid output to %zu tag(s) from %s.",
          output_tags.size(), parameter_name);
  return true;
}
 
void build_tagged_nodal_output_restriction(sf_mesh& mesh,
                                           const hashmap::unordered_set<int>& output_tags,
                                           SF::vector<mesh_int_t>& output_idx,
                                           bool async)
{
  SF::vector<mesh_int_t> selected_nodes;
 
  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
    if(output_tags.count(mesh.tag[eidx]) == 0) continue;
 
    for(mesh_int_t j = mesh.dsp[eidx]; j < mesh.dsp[eidx + 1]; j++)
      selected_nodes.push_back(mesh.con[j]);
  }
 
  localize_gridout_nodes(mesh, selected_nodes, output_idx, async);
}
 
void validate_gridout_tags(hashmap::unordered_set<int>& output_tags,
                           const SF::vector<mesh_t>& mesh_ids)
{
  hashmap::unordered_set<int> local_seen;
 
  for(mesh_t mesh_id : mesh_ids) {
    if(!mesh_is_registered(mesh_id)) continue;
 
    const sf_mesh& mesh = get_mesh(mesh_id);
    for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
      if(output_tags.count(mesh.tag[eidx]) != 0)
        local_seen.insert(mesh.tag[eidx]);
    }
  }
 
  SF::vector<int> sorted_tags;
  sorted_tags.assign(output_tags.begin(), output_tags.end());
  binary_sort(sorted_tags);
 
  SF::vector<int> missing_tags;
  for(int tag : sorted_tags) {
    int local_found = local_seen.count(tag) ? 1 : 0;
    int global_found = 0;
    MPI_Allreduce(&local_found, &global_found, 1, MPI_INT, MPI_MAX, SF_COMM);
    if(global_found == 0)
      missing_tags.push_back(tag);
  }
 
  if(missing_tags.size()) {
    binary_sort(missing_tags);
 
    std::stringstream msg;
    for(size_t i = 0; i < missing_tags.size(); i++) {
      if(i) msg << ", ";
      msg << missing_tags[i];
    }
 
    log_msg(0, 3, ECHO, "Warning: ignoring gridout_tags not present in the selected output mesh(es): %s.",
            msg.str().c_str());
 
    for(int tag : missing_tags)
      output_tags.erase(tag);
 
    if(output_tags.size() == 0) {
      log_msg(0, 5, ECHO, "gridout_tags did not match any tag in the selected output mesh(es).");
      EXIT(EXIT_FAILURE);
    }
  }
}
 
void intersect_output_restriction(SF::vector<mesh_int_t>*& restr,
                                  const SF::vector<mesh_int_t>& selected_idx)
{
  if(restr == NULL) {
    restr = new SF::vector<mesh_int_t>(selected_idx);
    return;
  }
 
  binary_sort(*restr);
  SF::vector<mesh_int_t> selected(selected_idx);
  binary_sort(selected);
 
  SF::vector<mesh_int_t> intersection;
  size_t lhs = 0, rhs = 0;
  while(lhs < restr->size() && rhs < selected.size()) {
    if((*restr)[lhs] == selected[rhs]) {
      intersection.push_back((*restr)[lhs]);
      lhs++;
      rhs++;
    } else if((*restr)[lhs] < selected[rhs]) {
      lhs++;
    } else {
      rhs++;
    }
  }
 
  *restr = intersection;
}
 
 
void output_meshes()
{
  bool write_intra_elec = mesh_is_registered(intra_elec_msh) && param_globals::gridout_i;
  bool write_extra_elec = mesh_is_registered(extra_elec_msh) && param_globals::gridout_e;
 
  set_dir(OUTPUT);
  std::string output_base = get_basename(param_globals::meshname);
  hashmap::unordered_set<int> output_tags;
  const bool restrict_gridout =
      (write_intra_elec || write_extra_elec) && parse_gridout_tags(output_tags);
  if(restrict_gridout) {
    SF::vector<mesh_t> mesh_ids;
    if(write_intra_elec) mesh_ids.push_back(intra_elec_msh);
    if(write_extra_elec) mesh_ids.push_back(extra_elec_msh);
    validate_gridout_tags(output_tags, mesh_ids);
  }
 
  if(write_intra_elec) {
    sf_mesh & mesh = get_mesh(intra_elec_msh);
    sf_mesh restricted_mesh;
    sf_mesh* output_mesh = &mesh;
    if(restrict_gridout) {
      if(!mesh_has_gridout_tags(mesh, output_tags)) {
        log_msg(0, 5, ECHO, "gridout_tags selected no intracellular grid elements.");
        EXIT(EXIT_FAILURE);
      }
      extract_gridout_tag_mesh(mesh, output_tags, restricted_mesh);
      output_mesh = &restricted_mesh;
    }
 
    if(param_globals::gridout_i & 1) {
      if(param_globals::output_level > 1)
        log_msg(0,0,0, "Computing \"%s\" surface ..", output_mesh->name.c_str());
 
      std::string output_file = output_base + "_i.surf";
      log_msg(0,0,0, "Writing \"%s\" surface: %s", output_mesh->name.c_str(), output_file.c_str());
 
      if(param_globals::gridout_i & 2) {
        // the surface indexes into the points of the volume mesh written below
        sf_mesh surfmesh;
        compute_surface_mesh(*output_mesh, SF::NBR_SUBMESH, surfmesh);
        write_surface(surfmesh, output_file);
      }
      else {
        // no volume mesh is written, so the surface needs points of its own
        write_gridout_surface(*output_mesh, output_base + "_i");
      }
    }
    if(param_globals::gridout_i & 2) {
      bool write_binary = false;
 
      std::string output_file = output_base + "_i";
      log_msg(0,0,0, "Writing \"%s\" mesh: %s", output_mesh->name.c_str(), output_file.c_str());
      write_mesh_parallel(*output_mesh, write_binary, output_file.c_str());
    }
  }
 
  if(write_extra_elec) {
    sf_mesh & mesh = get_mesh(extra_elec_msh);
    sf_mesh restricted_mesh;
    sf_mesh* output_mesh = &mesh;
    if(restrict_gridout) {
      if(!mesh_has_gridout_tags(mesh, output_tags)) {
        log_msg(0, 5, ECHO, "gridout_tags selected no extracellular grid elements.");
        EXIT(EXIT_FAILURE);
      }
      extract_gridout_tag_mesh(mesh, output_tags, restricted_mesh);
      output_mesh = &restricted_mesh;
    }
 
    if(param_globals::gridout_e & 1) {
      if(param_globals::output_level > 1)
        log_msg(0,0,0, "Computing \"%s\" surface ..", output_mesh->name.c_str());
 
      std::string output_file = output_base + "_e.surf";
      log_msg(0,0,0, "Writing \"%s\" surface: %s", output_mesh->name.c_str(), output_file.c_str());
 
      if(param_globals::gridout_e & 2) {
        // the surface indexes into the points of the volume mesh written below
        sf_mesh surfmesh;
        compute_surface_mesh(*output_mesh, SF::NBR_SUBMESH, surfmesh);
        write_surface(surfmesh, output_file);
      }
      else {
        // no volume mesh is written, so the surface needs points of its own
        write_gridout_surface(*output_mesh, output_base + "_e");
      }
    }
    if(param_globals::gridout_e & 2) {
      bool write_binary = false;
      std::string output_file = output_base + "_e";
      log_msg(0,0,0, "Writing \"%s\" mesh: %s", output_mesh->name.c_str(), output_file.c_str());
      write_mesh_parallel(*output_mesh, write_binary, output_file.c_str());
    }
  }
}
 
[[noreturn]] void cleanup_and_exit()
{
  destroy_physics();
  user_globals::scatter_reg.free_scatterings();
  print_active_citation_suggestions();
 
  const paramschema::ResetResult reset = paramschema::reset_schema_state(paramschema::openCARP_schema());
  print_lines(stderr, "parameter cleanup warning: ", reset.warnings);
  print_lines(stderr, "parameter cleanup error: ", reset.errors);
  finalize_backend();
 
  // close petsc error FD
  if(user_globals::petsc_error_fd)
    fclose(user_globals::petsc_error_fd);
 
  exit(EXIT_SUCCESS);
}
 
char* get_file_dir(const char* file)
{
  char* filecopy = dupstr(file);
  char* dir      = dupstr(dirname(filecopy));
 
  free(filecopy);
  return dir;
}
 
void setup_petsc_err_log()
{
  set_dir(OUTPUT);
 
#ifdef WITH_PETSC
  if(!petsc_active()) return;
 
  if(get_rank() == 0) {
    // we open a error log file handle and set it as petsc stderr
    user_globals::petsc_error_fd = fopen("petsc_err_log.txt", "w");
    petsc_redirect_errors(user_globals::petsc_error_fd);
  }
  else {
    petsc_redirect_errors(NULL);
  }
#endif
}
 
short get_mesh_dim(mesh_t id)
{
  const sf_mesh & mesh = get_mesh(id);
  SF::element_view<mesh_int_t, mesh_real_t> view(mesh, SF::NBR_PETSC);
  short mindim = 3;
 
  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
    view.set_elem(eidx);
    short cdim = view.dimension();
    if(mindim < cdim) mindim = cdim;
  }
 
  mindim = get_global(mindim, MPI_MIN, mesh.comm);
 
  return mindim;
}
 
void igb_output_manager::register_output_sync(sf_vec* inp_data,
                                              const mesh_t inp_meshid,
                                              const int dpn,
                                              const char* name,
                                              const char* units,
                                              const SF::vector<mesh_int_t>* idx,
                                              bool elem_data)
{
  sync_io_item IO;
 
  IO.data      = inp_data;
  IO.elem_flag = elem_data;
  IO.restr_idx = idx;
  if(idx != NULL) {
    const sf_mesh& mesh = get_mesh(inp_meshid);
    const SF::vector<mesh_int_t>& layout = elem_data ? mesh.epl.algebraic_layout() :
                                                       mesh.pl.algebraic_layout();
    const mesh_int_t start = layout[get_rank()];
 
    IO.restr_global_idx.resize(idx->size());
    for(size_t i = 0; i < idx->size(); i++)
      IO.restr_global_idx[i] = start + (*idx)[i];
  }
 
  IGBheader regigb;
  const timer_manager & tm = *user_globals::tm_manager;
  const int num_io = tm.timers[iotm_spacedt]->numIOs;
  int       err    = 0;
 
  int gsize = inp_data->gsize();
 
  // if we are restricting, we have to compute the restricted global size
  if(idx != NULL)
    gsize = get_global(int(idx->size()), MPI_SUM) * dpn;
 
  regigb.x(gsize / dpn);
  regigb.dim_x(regigb.x()-1);
  regigb.inc_x(1);
 
  regigb.y(1); regigb.z(1);
  regigb.t(num_io);
  regigb.dim_t(tm.end-tm.start);
 
  switch(dpn) {
    default:
    case 1: regigb.type(IGB_FLOAT); break;
    case 3: regigb.type(IGB_VEC3_f); break;
    case 4: regigb.type(IGB_VEC4_f); break;
    case 9: regigb.type(IGB_VEC9_f); break;
  }
 
  regigb.unites_x("um"); regigb.unites_y("um"); regigb.unites_z("um");
  regigb.unites_t("ms");
  regigb.unites(units);
  regigb.org_t(tm.time);
  regigb.inc_t(param_globals::spacedt);
 
  if(get_rank() == 0) {
    FILE_SPEC file = f_open(name, "w");
    if(file != NULL) {
      regigb.fileptr(file->fd);
      regigb.write();
      delete file;
    }
    else err++;
  }
 
  err = get_global(err, MPI_SUM);
  if(err) {
    log_msg(0,5,0, "%s error: Could not set up data output! Aborting!", __func__);
    EXIT(1);
  }
 
  IO.igb = regigb;
 
  SF::mixed_tuple<mesh_t, int> mesh_spec = {inp_meshid, dpn};
  IO.spec = mesh_spec;
 
  if(elem_data) {
    if(buffmap_elem.find(mesh_spec) == buffmap_elem.end()) {
      sf_vec *inp_copy; SF::init_vector(&inp_copy, inp_data);
      buffmap_elem[mesh_spec] = inp_copy;
    }
  } else {
    if(buffmap.find(mesh_spec) == buffmap.end()) {
      sf_vec *inp_copy; SF::init_vector(&inp_copy, inp_data);
      buffmap[mesh_spec] = inp_copy;
    }
  }
 
  this->sync_IOs.push_back(IO);
}
 
void igb_output_manager::register_output_async(sf_vec* inp_data,
                                               const mesh_t inp_meshid,
                                               const int dpn,
                                               const char* name,
                                               const char* units,
                                               const SF::vector<mesh_int_t>* idx,
                                               bool elem_data)
{
  sf_mesh & mesh = get_mesh(inp_meshid);
  SF::vector<mesh_int_t> ioidx;
  int rank = get_rank();
 
  async_io_item IO;
  IO.data = inp_data;
  IO.restr_idx = idx;
 
  if(elem_data) {
    const SF::vector<mesh_int_t> & nbr = mesh.get_numbering(SF::NBR_ELEM_SUBMESH);
    ioidx.resize(mesh.l_numelem);
    for(size_t i=0; i<mesh.l_numelem; i++)
      ioidx[i] = nbr[i];
  } else {
    const SF::vector<mesh_int_t> & alg_nod = mesh.pl.algebraic_nodes();
    const SF::vector<mesh_int_t> & nbr     = mesh.get_numbering(SF::NBR_SUBMESH);
 
    if(idx == NULL) {
      ioidx.resize(alg_nod.size());
 
      for(size_t i=0; i<alg_nod.size(); i++)
        ioidx[i] = nbr[alg_nod[i]];
    } else {
      ioidx.resize(idx->size());
      IO.restr_petsc_idx.resize(idx->size());
 
      for(size_t i=0; i<idx->size(); i++) {
        mesh_int_t loc_nodal = (*idx)[i];
        ioidx[i] = nbr[loc_nodal];
        IO.restr_petsc_idx[i] = local_nodal_to_local_petsc(mesh, rank, loc_nodal);
      }
    }
  }
 
  int id = async::COMPUTE_register_output(ioidx, dpn, name, units);
  IO.IO_id = id;
 
  this->async_IOs.push_back(IO);
}
 
void igb_output_manager::register_output(sf_vec* inp_data,
                                         const mesh_t inp_meshid,
                                         const int dpn,
                                         const char* name,
                                         const char* units,
                                         const SF::vector<mesh_int_t>* idx,
                                         bool elem_data)
{
  if(param_globals::num_io_nodes == 0)
    register_output_sync(inp_data, inp_meshid, dpn, name, units, idx, elem_data);
  else
    register_output_async(inp_data, inp_meshid, dpn, name, units, idx, elem_data);
}
 
sf_vec* igb_output_manager::fill_output_buffer(const sync_io_item & it)
{
  const SF::mixed_tuple<mesh_t, int> & cspec = it.spec;
  sf_vec* data_vec = it.data;
 
  bool have_perm = it.elem_flag ? have_permutation(cspec.v1, ELEM_PETSC_TO_CANONICAL, cspec.v2):
                                  have_permutation(cspec.v1, PETSC_TO_CANONICAL, cspec.v2);
 
  if(have_perm) {
    sf_vec*   perm_vec = it.elem_flag ? this->buffmap_elem[cspec] : this->buffmap[cspec];
    SF::scattering* sc       = it.elem_flag ? get_permutation(cspec.v1, ELEM_PETSC_TO_CANONICAL, cspec.v2):
                                              get_permutation(cspec.v1, PETSC_TO_CANONICAL, cspec.v2);
    sc->forward(*data_vec, *perm_vec);
    return perm_vec;
  } else {
    return data_vec;
  }
}
 
void igb_output_manager::write_data()
{
  SF::vector<float>      restr_buff;
  int rank = get_rank();
  // loop over registered datasets and root-write one by one
  //
  for (sync_io_item & it : sync_IOs) {
    // write to associated file descriptor
    FILE* fd = static_cast<FILE*>(it.igb.fileptr());
 
    // fill the output buffer
    sf_vec* buff = fill_output_buffer(it);
 
    if(it.restr_idx == NULL) {
      buff->write_binary<float>(fd);
    } else {
      const SF::vector<mesh_int_t> & idx = *it.restr_idx;
      const int dpn = it.spec.v2;
      SF_real* p = buff->ptr();
 
      restr_buff.resize(idx.size()); restr_buff.resize(0);
 
      for(mesh_int_t ii : idx) {
        const mesh_int_t offset = ii * dpn;
        for(int j = 0; j < dpn; j++)
          restr_buff.push_back(p[offset + j]);
      }
 
      if(dpn == 1) {
        root_write_ordered(fd, it.restr_global_idx, restr_buff, SF_COMM);
      } else {
        SF::vector<mesh_int_t> cnt(it.restr_global_idx.size(), dpn);
        root_write_ordered(fd, it.restr_global_idx, cnt, restr_buff, SF_COMM);
      }
      buff->release_ptr(p);
    }
  }
 
  // do all A-synchronous output:
  // loop over IDs received from the IO nodes and trigger async output
  //
  for (async_io_item & it : async_IOs) {
    SF_real* p    = it.data->ptr();
    int        ls = it.data->lsize();
    int        id = it.IO_id;
 
    if(it.restr_idx == NULL)
      async::COMPUTE_do_output(p, ls, id);
    else {
      async::COMPUTE_do_output(p, it.restr_petsc_idx, id);
    }
 
    it.data->release_ptr(p);
  }
}
 
void igb_output_manager::close_files_and_cleanup()
{
  if(get_rank() == 0) {
    // loop over registered datasets and close fd
    for(sync_io_item & it : sync_IOs) {
      FILE* fd = static_cast<FILE*>(it.igb.fileptr());
      fclose(fd);
    }
  }
 
  for(auto it = buffmap.begin(); it != buffmap.end(); ++it)
    delete it->second;
 
  for(auto it = buffmap_elem.begin(); it != buffmap_elem.end(); ++it)
    delete it->second;
 
  // we resize the arrays and clear the maps so that we are safe when calling
  // close_files_and_cleanup multiple times.
  sync_IOs.resize(0);
  async_IOs.resize(0);
  buffmap.clear();
  buffmap_elem.clear();
}
 
IGBheader* igb_output_manager::get_igb_header(const sf_vec* vec)
{
  for(sync_io_item & it : sync_IOs) {
    if(it.data == vec)
      return &it.igb;
  }
 
  return NULL;
}
 
namespace {
 
bool render_output_parameter_arguments(int                       argc,
                                       char**                    argv,
                                       std::vector<std::string>* rendered_values)
{
  rendered_values->assign(static_cast<std::size_t>(argc), std::string());
 
  for (int i = 1; i < argc; ++i) {
    const std::string argument(argv[i]);
    if (argument == "+F" || argument.find("_options_file") != std::string::npos) {
      if (i + 1 >= argc) {
        std::fprintf(stderr, "parameter file error: missing value for '%s'\n", argument.c_str());
        return false;
      }
      if (argument.find("_options_file") != std::string::npos) {
        std::string render_error;
        if (!render_parameter_file_value(argv[i + 1], false, &(*rendered_values)[i], &render_error)) {
          std::fprintf(stderr, "parameter file error: cannot render '%s': %s\n",
                       argument.substr(1).c_str(), render_error.c_str());
          return false;
        }
      }
      ++i;
      continue;
    }
 
    if (argument.empty() || argument[0] != '-') {
      continue;
    }
 
    const bool prelast = (i == argc - 1);
    const bool parameter_follows = !prelast &&
                                   ((argv[i + 1][0] != '-') ||
                                    ((argv[i + 1][1] >= '0') && (argv[i + 1][1] <= '9')));
    if (!parameter_follows) {
      continue;
    }
 
    const std::string raw_value(argv[i + 1]);
    std::string::size_type front = 0;
    while (front < raw_value.size() && raw_value[front] == ' ') {
      ++front;
    }
 
    bool        aggregate_value = false;
    std::string value = raw_value.substr(front);
    if (front < raw_value.size() && raw_value[front] == '{') {
      aggregate_value = true;
      ++front;
      while (front < raw_value.size() && raw_value[front] == ' ') {
        ++front;
      }
 
      std::string::size_type back = raw_value.size();
      while (back > front && raw_value[back - 1] == ' ') {
        --back;
      }
      if (back > front && raw_value[back - 1] == '}') {
        --back;
        value = raw_value.substr(front, back - front);
      } else {
        value = raw_value.substr(front);
      }
    }
 
    std::string render_error;
    if (!render_parameter_file_value(value, aggregate_value, &(*rendered_values)[i], &render_error)) {
      std::fprintf(stderr, "parameter file error: cannot render '%s': %s\n",
                   argv[i] + 1, render_error.c_str());
      return false;
    }
    ++i;
  }
 
  return true;
}
 
}  // namespace
 
bool output_parameter_file(const char *fname, int argc, char **argv)
{
  const int   max_line_len = 128;
  const char* file_sep = "#=======================================================";
 
  // make sure only root executes this function
  if(mpi_runtime_ready() && get_rank() != 0)
    return true;
 
  std::vector<std::string> rendered_values;
  if (!render_output_parameter_arguments(argc, argv, &rendered_values)) {
    return false;
  }
 
  FILE_SPEC out = f_open(fname, "w");
  fprintf(out->fd, "# CARP GIT commit hash: %s\n", GIT_COMMIT_HASH);
  fprintf(out->fd, "# dependency hashes:    %s\n", SUBREPO_COMMITS);
  fprintf(out->fd, "\n");
 
  // output the command line
  char line[8196] = "# ";
 
  for (int j=0; j<argc; j++) {
    strcat(line, argv[j]);
    if(strlen(line) > max_line_len) {
      fprintf(out->fd, "%s\n", line);
      strcpy(line, "# ");
    } else
      strcat(line, " ");
  }
 
  fprintf(out->fd, "%s\n\n", line);
  set_dir(INPUT);
 
  // convert command line to a par file
  for (int i=1; i<argc; i++) {
    std::string argument(argv[i]);
    if (argument == "+F" || argument.find("_options_file")!= std::string::npos) {
 
      std::string init = "";
      if (argument.find("_options_file")!= std::string::npos) {
        fprintf(out->fd, "%s = %s\n", argument.substr(1).c_str(), rendered_values[i].c_str());
        init = "#";
      }
      fprintf(out->fd, "%s>>\n", file_sep);
      // import par files
      i++;
      fprintf(out->fd, "## %s ##\n", argv[i]);
      FILE *in = fopen(argv[i], "r");
      while (fgets(line, 8196, in))
        fprintf(out->fd, "%s%s", init.c_str(), line);
      fclose(in);
      fprintf(out->fd, "\n##END of %s\n", argv[i]);
      fprintf(out->fd, "%s<<\n\n", file_sep);
    }
    else if(argv[i][0] == '-')
    {
      bool prelast       = (i==argc-1);
      bool paramFollows  = !prelast && ((argv[i+1][0] != '-') ||
                           ((argv[i+1][1] >= '0') && (argv[i+1][1] <= '9')));
 
      // strip leading hyphens from command line opts
      // assume options do not start with numbers
      if(paramFollows) {
        fprintf(out->fd, "%-40s= %s\n", argv[i]+1, rendered_values[i].c_str());
        i++;
      }
      else // a flag was specified
        fprintf(out->fd, "%-40s= 1\n", argv[i]);
    }
  }
  f_close(out);
  return true;
}
 
bool render_parameter_file_value(const std::string& value,
                                 bool               aggregate_value,
                                 std::string*       rendered,
                                 std::string*       error)
{
  if (aggregate_value) {
    *rendered = value;
    if (error != NULL) {
      error->clear();
    }
    return true;
  }
  return paramschema::quote_legacy_par_value(value, rendered, error);
}
 
void savequit()
{
  if(!get_physics(elec_phys)) return;
 
  set_dir(OUTPUT);
 
  double time = user_globals::tm_manager->time;
  char save_fn[512];
 
  snprintf(save_fn, sizeof save_fn, "exit.save.%.3f.roe", time);
  log_msg(NULL, 0, 0, "savequit called at time %g\n", time);
 
  Electrics* elec = (Electrics*) get_physics(elec_phys);
  elec->ion.miif->dump_state(save_fn, time, intra_elec_msh, false, GIT_COMMIT_COUNT);
 
  cleanup_and_exit();
}
 
}  // namespace opencarp