v0.14.0
unsaturated_transport.cpp
Go to the documentation of this file.
1 /** \file unsaturated_transport.cpp
2 \brief Implementation of unsaturated flow problem
3 \example unsaturated_transport.cpp
4 
5 \ingroup mofem_mix_transport_elem
6 */
7 
8 
9 
10 
11 #include <BasicFiniteElements.hpp>
12 #include <MixTransportElement.hpp>
13 #include <UnsaturatedFlow.hpp>
14 #include <MaterialUnsaturatedFlow.hpp>
15 
16 #ifndef WITH_ADOL_C
17 #error "MoFEM need to be compiled with ADOL-C"
18 #endif
19 
20 using namespace MoFEM;
21 using namespace MixTransport;
22 static char help[] = "...\n\n";
23 
24 double GenericMaterial::ePsilon0 = 0;
25 double GenericMaterial::ePsilon1 = 0;
26 double GenericMaterial::scaleZ = 1;
27 
28 map<std::string, CommonMaterialData::RegisterHook>
29  RegisterMaterials::mapOfRegistredMaterials;
30 
31 int main(int argc, char *argv[]) {
32 
33  const string default_options = "-ksp_type fgmres\n"
34  "-pc_type lu \n"
35  "-pc_factor_mat_solver_type mumps \n"
36  "-mat_mumps_icntl_20 0 \n";
37 
38  string param_file = "param_file.petsc";
39  if (!static_cast<bool>(ifstream(param_file))) {
40  std::ofstream file(param_file.c_str(), std::ios::ate);
41  if (file.is_open()) {
42  file << default_options;
43  file.close();
44  }
45  }
46 
47  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
48 
49  // Register DM Manager
50  CHKERR DMRegister_MoFEM("DMMOFEM"); // register MoFEM DM in PETSc
51  // Register materials
52  CHKERR RegisterMaterials()();
53 
54  PetscBool test_mat = PETSC_FALSE;
55  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_mat", &test_mat,
56  PETSC_NULL);
57 
58  if (test_mat == PETSC_TRUE) {
59  CommonMaterialData data;
60  // Testing van_genuchten
61  MaterialVanGenuchten van_genuchten(data);
62  van_genuchten.printMatParameters(-1, "testing");
63  van_genuchten.printTheta(-1e-16, -1e4, -1e-12, "hTheta");
64  van_genuchten.printKappa(-1e-16, -1, -1e-12, "hK");
65  van_genuchten.printC(-1e-16, -1, -1e-12, "hC");
66  CHKERR PetscFinalize();
67  return 0;
68  }
69 
70  try {
71 
72  moab::Core mb_instance;
73  moab::Interface &moab = mb_instance;
74 
75  // get file name form command line
76  PetscBool flg_mesh_file_name = PETSC_TRUE;
77  char mesh_file_name[255];
78  PetscBool flg_conf_file_name = PETSC_TRUE;
79  char conf_file_name[255];
80 
81  int order = 1;
82 
83  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Unsaturated flow options",
84  "none");
85  CHKERRG(ierr);
86  ierr = PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
87  mesh_file_name, 255, &flg_mesh_file_name);
88  CHKERRG(ierr);
89  ierr = PetscOptionsString(
90  "-configure", "material and bc configuration file name", "",
91  "unsaturated.cfg", conf_file_name, 255, &flg_conf_file_name);
92  CHKERRG(ierr);
93  ierr = PetscOptionsInt("-my_order", "default approximation order", "",
94  order, &order, PETSC_NULL);
95  CHKERRG(ierr);
96  ierr = PetscOptionsEnd();
97  CHKERRG(ierr);
98 
99  if (flg_mesh_file_name != PETSC_TRUE) {
100  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
101  "File name not set (e.g. -my_name mesh.h5m)");
102  }
103  if (flg_conf_file_name != PETSC_TRUE) {
104  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
105  "File name not set (e.g. -config unsaturated.cfg)");
106  }
107 
108  // Create MOAB communicator
109  auto pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
110  auto moab_comm_wrap =
111  boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
112  if (pcomm == NULL)
113  pcomm =
114  new ParallelComm(&moab, moab_comm_wrap->get_comm());
115 
116  const char *option;
117  option = "PARALLEL=READ_PART;"
118  "PARALLEL_RESOLVE_SHARED_ENTS;"
119  "PARTITION=PARALLEL_PARTITION;";
120  CHKERR moab.load_file(mesh_file_name, 0, option);
121 
122  // Create mofem interface
123  MoFEM::Core core(moab);
124  MoFEM::Interface &m_field = core;
125 
126  // Add meshsets with material and boundary conditions
127  MeshsetsManager *meshsets_manager_ptr;
128  CHKERR m_field.getInterface(meshsets_manager_ptr);
129  CHKERR meshsets_manager_ptr->setMeshsetFromFile();
130 
131  PetscPrintf(PETSC_COMM_WORLD,
132  "Read meshsets and added meshsets from bc.cfg\n");
133  for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, it)) {
134  PetscPrintf(PETSC_COMM_WORLD, "%s",
135  static_cast<std::ostringstream &>(
136  std::ostringstream().seekp(0) << *it << endl)
137  .str()
138  .c_str());
139  }
140 
141  // Set entities bit level
142  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
143  0, 3, BitRefLevel().set(0));
144 
145  UnsaturatedFlowElement uf(m_field);
146 
147  ifstream ini_file(conf_file_name);
148  po::variables_map vm;
149  po::options_description config_file_options;
150  Range domain_ents, bc_boundary_ents;
151 
152  map<int, CommonMaterialData> material_blocks;
153  map<int, double> head_blocks;
154  map<int, double> flux_blocks;
155  // Set material blocks
156  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
157  if (it->getName().compare(0, 4, "SOIL") != 0)
158  continue;
159  // get block id
160  const int block_id = it->getMeshsetId();
161  std::string block_name =
162  "mat_block_" + boost::lexical_cast<std::string>(block_id);
163  material_blocks[block_id].blockId = block_id;
164  material_blocks[block_id].addOptions(config_file_options, block_name);
165  }
166  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
167  if (it->getName().compare(0, 4, "HEAD") != 0)
168  continue;
169  // get block id
170  const int block_id = it->getMeshsetId();
171  std::string block_name =
172  "head_block_" + boost::lexical_cast<std::string>(block_id);
173  config_file_options.add_options()(
174  (block_name + ".head").c_str(),
175  po::value<double>(&head_blocks[block_id])->default_value(0));
176  }
177  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
178  if (it->getName().compare(0, 4, "FLUX") != 0)
179  continue;
180  // get block id
181  const int block_id = it->getMeshsetId();
182  std::string block_name =
183  "flux_block_" + boost::lexical_cast<std::string>(block_id);
184  config_file_options.add_options()(
185  (block_name + ".flux").c_str(),
186  po::value<double>(&head_blocks[block_id])->default_value(0));
187  }
188  po::parsed_options parsed =
189  parse_config_file(ini_file, config_file_options, true);
190  store(parsed, vm);
191  po::notify(vm);
192  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
193  if (it->getName().compare(0, 4, "SOIL") != 0)
194  continue;
195  const int block_id = it->getMeshsetId();
196  uf.dMatMap[block_id] = RegisterMaterials::mapOfRegistredMaterials.at(
197  material_blocks[block_id].matName)(material_blocks.at(block_id));
198  if (!uf.dMatMap.at(block_id)) {
199  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
200  "Material block not set");
201  }
202  // get block test
203  CHKERR m_field.get_moab().get_entities_by_type(
204  it->meshset, MBTET, uf.dMatMap.at(block_id)->tEts, true);
205  domain_ents.merge(uf.dMatMap.at(block_id)->tEts);
206  uf.dMatMap.at(block_id)->printMatParameters(block_id, "Read material");
207  }
208  std::vector<std::string> additional_parameters;
209  additional_parameters =
210  collect_unrecognized(parsed.options, po::include_positional);
211  for (std::vector<std::string>::iterator vit = additional_parameters.begin();
212  vit != additional_parameters.end(); vit++) {
213  CHKERR PetscPrintf(m_field.get_comm(),
214  "** WARNING Unrecognized option %s\n", vit->c_str());
215  }
216  // Set capillary pressure bc data
217  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
218  if (it->getName().compare(0, 4, "HEAD") != 0)
219  continue;
220  // get block id
221  const int block_id = it->getMeshsetId();
222  // create block data instance
223  uf.bcValueMap[block_id] =
224  boost::shared_ptr<UnsaturatedFlowElement::BcData>(
225  new UnsaturatedFlowElement::BcData());
226  // get bc value
227  std::vector<double> attributes;
228  CHKERR it->getAttributes(attributes);
229  if (attributes.size() != 1)
230  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
231  "Wrong number of head parameters %d", attributes.size());
232  uf.bcValueMap[block_id]->fixValue = attributes[0];
233  std::string block_name =
234  "head_block_" + boost::lexical_cast<std::string>(block_id);
235  if (vm.count((block_name) + ".head")) {
236  uf.bcValueMap[block_id]->fixValue = head_blocks[block_id];
237  }
238  // cerr << uf.bcValueMap[block_id]->fixValue << endl;
239  // CHKERR it->printAttributes(std::cout);
240  // get faces in the block
241  CHKERR m_field.get_moab().get_entities_by_type(
242  it->meshset, MBTRI, uf.bcValueMap[block_id]->eNts, true);
243  bc_boundary_ents.merge(uf.bcValueMap[block_id]->eNts);
244  }
245 
246  int max_flux_id = 0;
247  // Set water flux bc data
248  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
249  if (it->getName().compare(0, 4, "FLUX") != 0)
250  continue;
251  // get block id
252  const int block_id = it->getMeshsetId();
253  // create block data instance
254  uf.bcFluxMap[block_id] =
255  boost::shared_ptr<UnsaturatedFlowElement::BcData>(
256  new UnsaturatedFlowElement::BcData());
257  // get bc value
258  std::vector<double> attributes;
259  CHKERR it->getAttributes(attributes);
260  // CHKERR it->printAttributes(std::cout);
261  uf.bcFluxMap[block_id]->fixValue = attributes[0];
262  std::string block_name =
263  "flux_block_" + boost::lexical_cast<std::string>(block_id);
264  if (vm.count((block_name) + ".flux")) {
265  uf.bcValueMap[block_id]->fixValue = head_blocks[block_id];
266  }
267  // get faces in the block
268  CHKERR m_field.get_moab().get_entities_by_type(
269  it->meshset, MBTRI, uf.bcFluxMap[block_id]->eNts, true);
270  bc_boundary_ents.merge(uf.bcFluxMap[block_id]->eNts);
271  max_flux_id = max_flux_id > block_id ? max_flux_id : block_id + 1;
272  }
273 
274  // Add zero flux on the rest of the boundary
275  auto get_zero_flux_entities = [&m_field,
276  &pcomm](const auto &domain_ents,
277  const auto &bc_boundary_ents) {
278  Range zero_flux_ents;
279  Skinner skin(&m_field.get_moab());
280  Range domain_skin;
281  CHKERR skin.find_skin(0, domain_ents, false, domain_skin);
282  // filter not owned entities, those are not on boundary
283  CHKERR pcomm->filter_pstatus(domain_skin,
284  PSTATUS_SHARED | PSTATUS_MULTISHARED,
285  PSTATUS_NOT, -1, &zero_flux_ents);
286  zero_flux_ents = subtract(zero_flux_ents, bc_boundary_ents);
287  return zero_flux_ents;
288  };
289 
290  CHKERR uf.addFields("VALUES", "FLUXES", order);
291  CHKERR uf.addFiniteElements("FLUXES", "VALUES");
292  CHKERR uf.buildProblem(
293  get_zero_flux_entities(domain_ents, bc_boundary_ents));
294  CHKERR uf.createMatrices();
295  CHKERR uf.setFiniteElements();
296  CHKERR uf.calculateEssentialBc();
297  CHKERR uf.calculateInitialPc();
298  CHKERR uf.solveProblem();
299  CHKERR uf.destroyMatrices();
300 
301  }
302  CATCH_ERRORS;
303 
304  CHKERR MoFEM::Core::Finalize();
305 
306  return 0;
307 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
main
int main(int argc, char *argv[])
Definition: unsaturated_transport.cpp:31
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
MixTransport::RegisterMaterials
Definition: MaterialUnsaturatedFlow.hpp:350
MixTransport::UnsaturatedFlowElement::BcData
Class storing information about boundary condition.
Definition: UnsaturatedFlow.hpp:145
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MixTransport::CommonMaterialData::printMatParameters
void printMatParameters(const int id, const std::string &prefix) const
Definition: MaterialUnsaturatedFlow.hpp:77
MixTransport::CommonMaterialData
Definition: MaterialUnsaturatedFlow.hpp:13
help
static char help[]
Definition: unsaturated_transport.cpp:22
MixTransport::MaterialVanGenuchten::printKappa
void printKappa(const double b, const double e, double s, const std::string &prefix)
Definition: MaterialUnsaturatedFlow.hpp:328
MixTransport::MaterialVanGenuchten::printTheta
void printTheta(const double b, const double e, double s, const std::string &prefix)
Definition: MaterialUnsaturatedFlow.hpp:315
BasicFiniteElements.hpp
MixTransport::MaterialVanGenuchten
Definition: MaterialUnsaturatedFlow.hpp:254
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MixTransport::MaterialVanGenuchten::printC
void printC(const double b, const double e, double s, const std::string &prefix)
Definition: MaterialUnsaturatedFlow.hpp:339
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
UnsaturatedFlow.hpp
Mix implementation of transport element.
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::MeshsetsManager::setMeshsetFromFile
MoFEMErrorCode setMeshsetFromFile(const string file_name, const bool clean_file_options=true)
add blocksets reading config file
Definition: MeshsetsManager.cpp:788
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
MixTransport::UnsaturatedFlowElement
Implementation of operators, problem and finite elements for unsaturated flow.
Definition: UnsaturatedFlow.hpp:102
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
Range
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:71
MixTransport
Definition: MaterialUnsaturatedFlow.hpp:11
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
MixTransportElement.hpp
Mix implementation of transport element.
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MaterialUnsaturatedFlow.hpp
Mix implementation of transport element.
_IT_CUBITMESHSETS_FOR_LOOP_
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
Definition: MeshsetsManager.hpp:34
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
Generated by Doxygen 1.8.17 and hosted at University of Glasgow