MoFEMManager.hpp

Go to the documentation of this file.

/* This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
namespace MoFEMManager {
 
struct MFManager {
 
  int order;
  int saveEveryNthStep;
  int atomTest;
  PetscBool isQuasiStatic;
  PetscBool isExplicit;
  PetscBool isPartitioned;
  PetscBool doPrintMaxMin;
 
  int lagJacFromNthSTep;
 
  MoFEM::Interface &mField;
  SmartPetscObj<TS> tSolver;
  SmartPetscObj<DM> dM;
  
  boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
 
  // FIXME: should we have a global (static) post-processing on volume and skin,
  // probably yes
  //  boost::shared_ptr<PostProcEle> postProc;
  //  boost::shared_ptr<PostProcSkinEle> postProcSkin;
 
  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> vScatter[3];
  boost::shared_ptr<FEMethod> monitorPtr;
 
  MFManager(MoFEM::Interface &m_field) : mField(m_field) {
    monitorPtr = boost::make_shared<FEMethod>();
    saveEveryNthStep = 1;
    order = 1;
    isPartitioned = PETSC_FALSE;
    atomTest = -1;
    doPrintMaxMin = PETSC_TRUE;
  }
 
  MoFEMErrorCode getParams();
  MoFEMErrorCode setMainField();
  MoFEMErrorCode defineDM();
  MoFEMErrorCode tsSetup();
  MoFEMErrorCode tsSolve();
 
  MoFEMErrorCode
  printMaxMin(std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> &tuple,
              const std::string msg, Vec &ts_u, double time);
 
  auto getTSType();
};
 
auto MFManager::getTSType() {
  auto t_type = GenericElementInterface::IM2;
  if (isQuasiStatic)
    t_type = GenericElementInterface::IM;
  if (isExplicit)
    t_type = GenericElementInterface::EX;
  return t_type;
}
 
MoFEMErrorCode MFManager::getParams() {
  MoFEMFunctionBegin;
  isQuasiStatic = PETSC_TRUE;
  isExplicit = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULL, "-is_quasi_static", &isQuasiStatic,
                             PETSC_NULL);
  CHKERR PetscOptionsGetBool(PETSC_NULL, "-is_explicit", &isExplicit,
                             PETSC_NULL);
  CHKERR PetscOptionsGetInt(PETSC_NULL, "-order", &order, PETSC_NULL);
  CHKERR PetscOptionsGetInt(PETSC_NULL, "-output_every", &saveEveryNthStep,
                            PETSC_NULL);                       
  CHKERR PetscOptionsGetInt(PETSC_NULL, "-atom_test", &atomTest, PETSC_NULL);
 
  CHKERR PetscOptionsGetBool(PETSC_NULL, "-print_max_min", &doPrintMaxMin,
                             PETSC_NULL);
  CHKERR PetscOptionsGetBool(PETSC_NULL, "-is_partitioned", &isPartitioned,
                             PETSC_NULL);
  MOFEM_LOG("WORLD", Sev::inform)
      << "Mesh Partition Flag Status: " << isPartitioned;
 
  MoFEMFunctionReturn(0);
}
 
 
MoFEMErrorCode MFManager::setMainField() {
  MoFEMFunctionBegin;
  // MoFEMFunctionReturnHot(0);
 
  // Select base
  enum bases { AINSWORTH, DEMKOWICZ, BERNSTEIN, LASBASETOPT };
  const char *list_bases[] = {"ainsworth", "demkowicz", "bernstein"};
  PetscInt choice_base_value = AINSWORTH;
  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
                              LASBASETOPT, &choice_base_value, PETSC_NULL);
 
  FieldApproximationBase base;
  switch (choice_base_value) {
  case AINSWORTH:
    base = AINSWORTH_LEGENDRE_BASE;
    MOFEM_LOG("WORLD", Sev::inform)
        << "Set AINSWORTH_LEGENDRE_BASE for displacements";
    break;
  case DEMKOWICZ:
    base = DEMKOWICZ_JACOBI_BASE;
    MOFEM_LOG("WORLD", Sev::inform)
        << "Set DEMKOWICZ_JACOBI_BASE for displacements";
    break;
  case BERNSTEIN:
    base = AINSWORTH_BERNSTEIN_BEZIER_BASE;
    MOFEM_LOG("WORLD", Sev::inform)
        << "Set AINSWORTH_BERNSTEIN_BEZIER_BASE for displacements";
    break;
  default:
    base = LASTBASE;
    break;
  }
 
  // Add displacement field
  CHKERR mField.add_field("U", H1, base, 3);
 
  // Add field representing ho-geometry
  CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
                          3);
 
  // Add entities to field
  CHKERR mField.add_ents_to_field_by_type(0, MBTET, "U");
  CHKERR mField.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
 
  // Set approximation order to entities
  if (base == AINSWORTH_BERNSTEIN_BEZIER_BASE)
    CHKERR mField.set_field_order(0, MBVERTEX, "U", order);
  else
    CHKERR mField.set_field_order(0, MBVERTEX, "U", 1);
  CHKERR mField.set_field_order(0, MBEDGE, "U", order);
  CHKERR mField.set_field_order(0, MBTRI, "U", order);
  CHKERR mField.set_field_order(0, MBTET, "U", order);
 
  CHKERR mField.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
  // CHKERR mField.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
  CHKERR mField.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 1);
  CHKERR mField.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 1);
  CHKERR mField.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 1);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MFManager::defineDM() {
  MoFEMFunctionBegin;
 
  auto simple = mField.getInterface<Simple>();
  DMType dm_name = "DMMOFEM";
  // Register DM problem
  CHKERR DMRegister_MoFEM(dm_name);
  dM = createSmartDM(mField.get_comm(), dm_name);
  CHKERR DMSetType(dM, dm_name);
  CHKERR DMMoFEMSetIsPartitioned(dM, isPartitioned);
  // Create DM instance
  CHKERR DMMoFEMCreateMoFEM(dM, &mField, simple->getProblemName().c_str(),
                            simple->getBitRefLevel());
  CHKERR DMSetFromOptions(dM);
 
  // mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
 
  MoFEMFunctionReturn(0);
}
 
 
MoFEMErrorCode MFManager::tsSetup() {
  MoFEMFunctionBegin;
 
  Simple *simple = mField.getInterface<Simple>();
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
  ISManager *is_manager = mField.getInterface<ISManager>();
 
  auto set_time_monitor = [&](auto solver) {
    MoFEMFunctionBegin;
    boost::shared_ptr<ForcesAndSourcesCore> null;
    CHKERR DMMoFEMTSSetMonitor(dM, solver, simple->getDomainFEName(),
                               monitorPtr, null, null);
    MoFEMFunctionReturn(0);
  };
 
  auto set_section_monitor = [&](auto solver) {
    MoFEMFunctionBegin;
    SNES snes;
    CHKERR TSGetSNES(solver, &snes);
    PetscViewerAndFormat *vf;
    CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
                                      PETSC_VIEWER_DEFAULT, &vf);
    CHKERR SNESMonitorSet(
        snes,
        (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
        vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
    MoFEMFunctionReturn(0);
  };
 
 
  auto set_dm_section = [&](auto dm) {
    MoFEMFunctionBeginHot;
    auto section = mField.getInterface<ISManager>()->sectionCreate(
        simple->getProblemName());
    CHKERR DMSetSection(dm, section);
    MoFEMFunctionReturnHot(0);
  };
 
  auto scatter_create = [&](auto coeff, auto D) {
    SmartPetscObj<IS> is;
    CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
                                                   ROW, "U", coeff, coeff, is);
    int loc_size;
    CHKERR ISGetLocalSize(is, &loc_size);
    Vec v;
    CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
    VecScatter scatter;
    CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
    return std::make_tuple(SmartPetscObj<Vec>(v),
                           SmartPetscObj<VecScatter>(scatter));
  };
 
  CHKERR set_dm_section(dM);
 
  auto D = smartCreateDMVector(dM);
  tSolver = MoFEM::createTS(mField.get_comm());
  
  //FIXME: for explicit you need to implement proper mass matrix for rhs operators
  if (isQuasiStatic || isExplicit) {
    CHKERR TSSetSolution(tSolver, D);
  } else {
    CHKERR TSSetType(tSolver, TSALPHA2);
    auto DD = smartVectorDuplicate(D);
    CHKERR TS2SetSolution(tSolver, D, DD);
  }
  CHKERR set_section_monitor(tSolver);
  CHKERR TSSetDM(tSolver, dM);
 
  //default max time is 1
  CHKERR TSSetMaxTime(tSolver, 1.0);
  CHKERR TSSetExactFinalTime(tSolver, TS_EXACTFINALTIME_MATCHSTEP);
  CHKERR TSSetFromOptions(tSolver);
 
  for (int i = 0; i != 3; ++i)
    vScatter[i] = scatter_create(i, D);
  CHKERR set_time_monitor(tSolver);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MFManager::tsSolve() {
  MoFEMFunctionBegin;
 
  CHKERR TSSetUp(tSolver);
  CHKERR TSSolve(tSolver, NULL);
 
  // CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
  // CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
  // CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode MFManager::printMaxMin(
    std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> &tuple,
    const std::string msg, Vec &ts_u, double time) {
  MoFEMFunctionBegin;
  double min_disp, max_disp;
 
  CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
                         INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
                       INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max_disp);
  CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min_disp);
  MOFEM_LOG_C("WORLD", Sev::inform, "%s time %3.4e min %3.4e max %3.4e\n",
              msg.c_str(), time, min_disp, max_disp);
 
  if (atomTest == 3 && time == 1.00)
    if (msg == "Ux") {
      MOFEM_LOG("WORLD", Sev::inform)
          << "Atom test min displacement X:" << min_disp;
      if (abs(min_disp + 0.2220) > 1e-4)
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "atom test diverged!");
    }
 
  MoFEMFunctionReturn(0);
 
}
} // namespace MoFEMManager