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