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