solution_vector_comparison_atom.cpp

Go to the documentation of this file.

/** @file solution_vector_comparison_atom.cpp
 * @brief Compare complete PETSc solution vectors from independent solves.
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, nullptr,
                          "Compare saved solution vectors in order.\n");
  try {
    char left_file[PETSC_MAX_PATH_LEN] = "";
    char right_file[PETSC_MAX_PATH_LEN] = "";
    PetscInt count = 1;
    PetscReal tolerance = 1.e-9;
    CHKERR PetscOptionsGetString(nullptr, nullptr, "-compare_left", left_file,
                                 sizeof(left_file), nullptr);
    CHKERR PetscOptionsGetString(nullptr, nullptr, "-compare_right", right_file,
                                 sizeof(right_file), nullptr);
    CHKERR PetscOptionsGetInt(nullptr, nullptr, "-compare_count", &count,
                              nullptr);
    CHKERR PetscOptionsGetReal(nullptr, nullptr, "-compare_tolerance",
                               &tolerance, nullptr);
    if (!left_file[0] || !right_file[0] || count < 1 ||
        !std::isfinite(tolerance) || tolerance <= 0.)
      SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
              "Specify both vector files, a positive count and tolerance");
    PetscViewer left_viewer, right_viewer;
    CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD, left_file, FILE_MODE_READ,
                                 &left_viewer);
    CHKERR PetscViewerBinaryOpen(PETSC_COMM_WORLD, right_file, FILE_MODE_READ,
                                 &right_viewer);
    for (PetscInt step = 0; step != count; ++step) {
      Vec left_raw, right_raw;
      CHKERR VecCreate(PETSC_COMM_WORLD, &left_raw);
      CHKERR VecCreate(PETSC_COMM_WORLD, &right_raw);
      SmartPetscObj<Vec> left(left_raw), right(right_raw);
      CHKERR VecLoad(left, left_viewer);
      CHKERR VecLoad(right, right_viewer);
      PetscInt left_size, right_size;
      CHKERR VecGetSize(left, &left_size);
      CHKERR VecGetSize(right, &right_size);
      if (left_size != right_size || !left_size)
        SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                "Solution vectors have different or empty layouts");
      PetscReal scale, difference;
      CHKERR VecNorm(left, NORM_INFINITY, &scale);
      CHKERR VecAXPY(right, -1., left);
      CHKERR VecNorm(right, NORM_INFINITY, &difference);
      if (!std::isfinite(scale) || !std::isfinite(difference) ||
          difference > tolerance * std::max(1., scale))
        SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                "Solution vector %" PetscInt_FMT
                " differs by %.16e (scale %.16e)",
                step, difference, scale);
      CHKERR PetscPrintf(PETSC_COMM_WORLD,
                         "Solution vector %" PetscInt_FMT ": %" PetscInt_FMT
                         " coefficients, maximum difference %.16e\n",
                         step, left_size, difference);
    }
    for (auto viewer : {left_viewer, right_viewer}) {
      PetscInt extra, count_read;
      CHKERR PetscViewerBinaryRead(viewer, &extra, 1, &count_read, PETSC_INT);
      if (count_read)
        SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                "Vector file contains more than the expected number of states");
    }
    CHKERR PetscViewerDestroy(&left_viewer);
    CHKERR PetscViewerDestroy(&right_viewer);
  }
  CATCH_ERRORS;
  MoFEM::Core::Finalize();
  return 0;
}