unsaturate_flow.cpp

Go to the documentation of this file.

#include <stdlib.h>
#include <BasicFiniteElements.hpp>
#include <UnsaturatedFlowOperators.hpp>
 
using namespace MoFEM;
using namespace UFOperators;
 
static char help[] = "...\n\n";
 
struct UFProblem {
public:
  UFProblem(moab::Core &mb_instance, MoFEM::Core &core, const int order,
            const int n_species)
      : moab(mb_instance), m_field(core), order(order), nb_species(n_species),
        cOmm(m_field.get_comm()), rAnk(m_field.get_comm_rank()) {
    vol_ele_stiff_rhs = boost::shared_ptr<VolEle>(new VolEle(m_field));
    vol_ele_stiff_lhs = boost::shared_ptr<VolEle>(new VolEle(m_field));
 
    boundary_ele_rhs = boost::shared_ptr<FaceEle>(new FaceEle(m_field));
 
    post_proc = boost::shared_ptr<PostProcVolumeOnRefinedMesh>(
        new PostProcVolumeOnRefinedMesh(m_field));
 
    data.resize(nb_species);
    values_ptr.resize(nb_species);
    grads_ptr.resize(nb_species);
    dots_ptr.resize(nb_species);
    inner_surface.resize(nb_species);
 
    for (int i = 0; i < nb_species; ++i) {
      data[i] = boost::shared_ptr<PreviousData>(new PreviousData());
      grads_ptr[i] = boost::shared_ptr<MatrixDouble>(data[i], &data[i]->grads);
      values_ptr[i] =
          boost::shared_ptr<VectorDouble>(data[i], &data[i]->values);
      dots_ptr[i] =
          boost::shared_ptr<VectorDouble>(data[i], &data[i]->dot_values);
    }
  }
 
  // UFProblem(const int order) : order(order){}
  MoFEMErrorCode run_analysis();
 
private:
  MoFEMErrorCode setup_system();
  MoFEMErrorCode add_fe(std::string field_name);
  MoFEMErrorCode set_blockData(std::map<int, BlockData> &block_data_map);
 
  MoFEMErrorCode set_initial_values(std::string field_name, int block_id,
                                    Range &surface, double &init_val);
 
  MoFEMErrorCode push_mass_ele(std::string field_name);
 
  MoFEMErrorCode update_vol_fe(boost::shared_ptr<VolEle> &vol_ele,
                               boost::shared_ptr<PreviousData> &data);
  MoFEMErrorCode update_stiff_rhs(std::string field_name,
                                  boost::shared_ptr<VectorDouble> &values_ptr,
                                  boost::shared_ptr<MatrixDouble> &grads_ptr,
                                  boost::shared_ptr<VectorDouble> &dots_ptr);
 
  MoFEMErrorCode push_stiff_rhs(std::string field_name,
                                boost::shared_ptr<PreviousData> &data,
                                std::map<int, BlockData> &block_map);
 
  MoFEMErrorCode update_stiff_lhs(std::string field_name,
                                  boost::shared_ptr<VectorDouble> &values_ptr,
                                  boost::shared_ptr<MatrixDouble> &grads_ptr);
  MoFEMErrorCode push_stiff_lhs(std::string field_name,
                                boost::shared_ptr<PreviousData> &data,
                                std::map<int, BlockData> &block_map);
 
  MoFEMErrorCode set_integration_rule();
 
  MoFEMErrorCode set_fe_in_loop();
  MoFEMErrorCode post_proc_fields(std::string field_name);
  MoFEMErrorCode output_result();
  MoFEMErrorCode solve();
 
  MoFEM::Interface &m_field;
  Simple *simple_interface;
 
  moab::Interface &moab;
 
  SmartPetscObj<DM> dm;
  SmartPetscObj<TS> ts;
 
  Range natural_bdry_ents;
 
  std::vector<Range> inner_surface;
 
  double global_error;
 
  MPI_Comm cOmm;
  const int rAnk;
 
  int order;
  int nb_species;
 
  std::map<int, BlockData> material_blocks;
 
  boost::shared_ptr<VolEle> vol_ele_stiff_rhs;
  boost::shared_ptr<VolEle> vol_ele_stiff_lhs;
  boost::shared_ptr<FaceEle> boundary_ele_rhs;
 
  boost::shared_ptr<VolEle> vol_mass_ele;
 
  boost::shared_ptr<PostProcVolumeOnRefinedMesh> post_proc;
  boost::shared_ptr<Monitor> monitor_ptr;
 
  std::vector<boost::shared_ptr<PreviousData>> data;
 
  std::vector<boost::shared_ptr<MatrixDouble>> grads_ptr;
  std::vector<boost::shared_ptr<VectorDouble>> values_ptr;
  std::vector<boost::shared_ptr<VectorDouble>> dots_ptr;
 
  boost::shared_ptr<ForcesAndSourcesCore> null;
};
 
MoFEMErrorCode UFProblem::setup_system() {
  MoFEMFunctionBegin;
  CHKERR m_field.getInterface(simple_interface);
  CHKERR simple_interface->getOptions();
  CHKERR simple_interface->loadFile();
  MoFEMFunctionReturn(0);
}
MoFEMErrorCode UFProblem::add_fe(std::string field_name) {
  MoFEMFunctionBegin;
  CHKERR simple_interface->addDomainField(field_name, H1,
                                          AINSWORTH_LEGENDRE_BASE, 1);
 
  CHKERR simple_interface->addBoundaryField(field_name, H1,
                                            AINSWORTH_LEGENDRE_BASE, 1);
  CHKERR simple_interface->setFieldOrder(field_name, order);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode UFProblem::set_blockData(std::map<int, BlockData> &block_map) {
  MoFEMFunctionBegin;
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
    string name = it->getName();
    const int id = it->getMeshsetId();
    if (name.compare(0, 14, "REGION1") == 0) {
      CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
          id, BLOCKSET, 3, block_map[id].block_ents, true);
 
      block_map[id].block_id = id;
      block_map[id].K_s = 1.000;
      block_map[id].h_s = 0.000;
      block_map[id].theta_s = 0.43000;
      block_map[id].theta_m = 0.43000;
      block_map[id].theta_r = 0.04500;
      block_map[id].alpha = 1.812500;
      block_map[id].nn = 5.3800;
 
    } else if (name.compare(0, 14, "REGION2") == 0) {
      CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
          id, BLOCKSET, 3, block_map[id].block_ents, true);
 
      block_map[id].block_id = id;
    } else if (name.compare(0, 14, "REGION3") == 0) {
      CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
          id, BLOCKSET, 3, block_map[id].block_ents, true);
 
      block_map[id].block_id = id;
    } else if (name.compare(0, 14, "REGION4") == 0) {
      CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
          id, BLOCKSET, 3, block_map[id].block_ents, true);
 
      block_map[id].block_id = id;
    } else if (name.compare(0, 14, "REGION5") == 0) {
      CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
          id, BLOCKSET, 3, block_map[id].block_ents, true);
 
      block_map[id].block_id = id;
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode UFProblem::set_initial_values(std::string field_name,
                                             int block_id, Range &surface,
                                             double &init_val) {
  MoFEMFunctionBegin;
  if (m_field.getInterface<MeshsetsManager>()->checkMeshset(block_id,
                                                            BLOCKSET)) {
    CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
        block_id, BLOCKSET, 3, surface, true);
  }
  if (!surface.empty()) {
    Range surface_verts;
 
    CHKERR moab.get_connectivity(surface, surface_verts, false);
    CHKERR m_field.getInterface<FieldBlas>()->setField(
        init_val, MBVERTEX, surface_verts, field_name);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode UFProblem::update_vol_fe(boost::shared_ptr<VolEle> &vol_ele,
                                        boost::shared_ptr<PreviousData> &data) {
  MoFEMFunctionBegin;
  auto det_ptr = boost::make_shared<VectorDouble>();
  auto jac_ptr = boost::make_shared<MatrixDouble>();
  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
  vol_ele->getOpPtrVector().push_back(new OpCalculateHOJacForFace(jac_ptr));
  vol_ele->getOpPtrVector().push_back(
      new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
  vol_ele->getOpPtrVector().push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
UFProblem::update_stiff_rhs(std::string field_name,
                            boost::shared_ptr<VectorDouble> &values_ptr,
                            boost::shared_ptr<MatrixDouble> &grads_ptr,
                            boost::shared_ptr<VectorDouble> &dots_ptr) {
 
  MoFEMFunctionBegin;
  vol_ele_stiff_rhs->getOpPtrVector().push_back(
      new OpCalculateScalarFieldValues(field_name, values_ptr));
  vol_ele_stiff_rhs->getOpPtrVector().push_back(
      new OpCalculateScalarValuesDot(field_name, dots_ptr));
  vol_ele_stiff_rhs->getOpPtrVector().push_back(
      new OpCalculateScalarFieldGradient<2>(field_name, grads_ptr));
  MoFEMFunctionReturn(0);
}
MoFEMErrorCode UFProblem::push_stiff_rhs(std::string field_name,
                                         boost::shared_ptr<PreviousData> &data,
                                         std::map<int, BlockData> &block_map) {
  MoFEMFunctionBegin;
 
  vol_ele_stiff_rhs->getOpPtrVector().push_back(
      new OpAssembleStiffRhs<3>(field_name, data, block_map));
  // boundary_ele_rhs->getOpPtrVector().push_back(
  //     new OpAssembleNaturalBCRhs(field_name,natural_bdry_ents));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
UFProblem::update_stiff_lhs(std::string field_name,
                            boost::shared_ptr<VectorDouble> &values_ptr,
                            boost::shared_ptr<MatrixDouble> &grads_ptr) {
  MoFEMFunctionBegin;
  vol_ele_stiff_lhs->getOpPtrVector().push_back(
      new OpCalculateScalarFieldValues(field_name, values_ptr));
 
  vol_ele_stiff_lhs->getOpPtrVector().push_back(
      new OpCalculateScalarFieldGradient<3>(field_name, grads_ptr));
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode UFProblem::push_stiff_lhs(std::string field_name,
                                         boost::shared_ptr<PreviousData> &data,
                                         std::map<int, BlockData> &block_map) {
  MoFEMFunctionBegin;
 
  vol_ele_stiff_lhs->getOpPtrVector().push_back(
      new OpAssembleStiffLhs<3>(field_name, data, block_map));
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode UFProblem::set_integration_rule() {
  MoFEMFunctionBegin;
  auto vol_rule = [](int, int, int p) -> int { return 2 * p; };
 
  vol_ele_stiff_rhs->getRuleHook = vol_rule;
 
  vol_ele_stiff_lhs->getRuleHook = vol_rule;
  boundary_ele_rhs->getRuleHook = vol_rule;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode UFProblem::set_fe_in_loop() {
  MoFEMFunctionBegin;
  CHKERR TSSetType(ts, TSARKIMEX);
  CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
 
  CHKERR DMMoFEMTSSetIJacobian(dm, simple_interface->getDomainFEName(),
                               vol_ele_stiff_lhs, null, null);
 
  CHKERR DMMoFEMTSSetIFunction(dm, simple_interface->getDomainFEName(),
                               vol_ele_stiff_rhs, null, null);
 
  CHKERR DMMoFEMTSSetIFunction(dm, simple_interface->getBoundaryFEName(),
                               boundary_ele_rhs, null, null);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode UFProblem::post_proc_fields(std::string field_name) {
  MoFEMFunctionBegin;
  post_proc->addFieldValuesPostProc(field_name);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode UFProblem::output_result() {
  MoFEMFunctionBegin;
  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple_interface->getDomainFEName(),
                             monitor_ptr, null, null);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode UFProblem::solve() {
  MoFEMFunctionBegin;
  // Create solution vector
  SmartPetscObj<Vec> X;
  CHKERR DMCreateGlobalVector_MoFEM(dm, X);
  CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_FORWARD);
  // Solve problem
  double ftime = 1;
  CHKERR TSSetDM(ts, dm);
  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
  CHKERR TSSetSolution(ts, X);
  CHKERR TSSetFromOptions(ts);
  CHKERR TSSolve(ts, X);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode UFProblem::run_analysis() {
  MoFEMFunctionBegin;
  global_error = 0;
  std::vector<std::string> mass_names(nb_species);
 
  for (int i = 0; i < nb_species; ++i) {
    mass_names[i] = "h" + boost::lexical_cast<std::string>(i + 1);
  }
  CHKERR setup_system();
  for (int i = 0; i < nb_species; ++i) {
    add_fe(mass_names[i]);
  }
 
  CHKERR simple_interface->setUp();
 
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
    string name = it->getName();
    if (name.compare(0, 14, "ESSENTIAL") == 0) {
      CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 2,
                                                 natural_bdry_ents, true);
    }
  }
  Range face_edges;
 
  CHKERR moab.get_adjacencies(natural_bdry_ents, 2, false, face_edges,
                              moab::Interface::UNION);
 
  Range face_edges_verts;
  CHKERR moab.get_connectivity(face_edges, face_edges_verts, false);
 
  Range bdry_ents;
  bdry_ents = unite(natural_bdry_ents, face_edges_verts);
  bdry_ents = unite(bdry_ents, face_edges_verts);
 
  bdry_ents.print();
  cout << "size: " << endl;
  cout << bdry_ents.size() << endl;
  cout << bdry_ents << endl;
 
  CHKERR m_field.getInterface<ProblemsManager>()->removeDofsOnEntities(
      "SimpleProblem", mass_names[0], bdry_ents);
 
  CHKERR set_blockData(material_blocks);
 
  VectorDouble initVals;
  initVals.resize(3, false);
  initVals.clear();
 
  initVals[0] = -0.8;
  initVals[1] = 0.0;
  initVals[2] = 0.0;
 
  for (int i = 0; i < nb_species; ++i) {
    CHKERR set_initial_values(mass_names[i], i + 2, inner_surface[i],
                              initVals[i]);
  }
 
  if (!natural_bdry_ents.empty()) {
    Range surface_verts;
 
    CHKERR moab.get_connectivity(natural_bdry_ents, surface_verts, false);
    CHKERR m_field.getInterface<FieldBlas>()->setField(
        0.0, MBVERTEX, surface_verts, mass_names[0]);
  }
 
  dm = simple_interface->getDM();
 
  // CHKERR update_vol_fe(vol_ele_stiff_rhs, data[0]);
 
  for (int i = 0; i < nb_species; ++i) {
    CHKERR update_stiff_rhs(mass_names[i], values_ptr[i], grads_ptr[i],
                            dots_ptr[i]);
    CHKERR push_stiff_rhs(mass_names[i], data[i], material_blocks);
  }
 
  // CHKERR update_vol_fe(vol_ele_stiff_lhs, data[0]);
 
  for (int i = 0; i < nb_species; ++i) {
    CHKERR update_stiff_lhs(mass_names[i], values_ptr[i], grads_ptr[i]);
    CHKERR push_stiff_lhs(mass_names[i], data[i],
                          material_blocks); // nb_species times
  }
  // vol_ele_stiff_lhs->getOpPtrVector().push_back(
  //   new OpError(ExactFunction(), ExactFunctionLap(), ExactFunctionGrad(),
  //   data[0], material_blocks, global_error));
 
  CHKERR set_integration_rule();
 
  ts = createTS(m_field.get_comm());
 
  CHKERR set_fe_in_loop();
 
  post_proc->generateReferenceElementMesh(); // only once
 
  for (int i = 0; i < nb_species; ++i) {
    CHKERR post_proc_fields(mass_names[i]);
  }
 
  monitor_ptr = boost::shared_ptr<Monitor>(
      new Monitor(cOmm, rAnk, dm, post_proc, global_error)); // nb_species times
  CHKERR output_result();                                    // only once
 
  CHKERR solve(); // only once
 
  MoFEMFunctionReturn(0);
}
 
int main(int argc, char *argv[]) {
  const char param_file[] = "param_file.petsc";
  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
  try {
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
    MoFEM::Core core(moab);
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    int order = 1;
    CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
    int nb_species = 1;
    UFProblem unsatu_flow_problem(mb_instance, core, order, nb_species);
    CHKERR unsatu_flow_problem.run_analysis();
  }
  CATCH_ERRORS;
  MoFEM::Core::Finalize();
  return 0;
}