v0.10.0
Functions | Variables
unsaturated_transport.cpp File Reference

Implementation of unsaturated flow problem. More...

#include <boost/program_options.hpp>
#include <BasicFiniteElements.hpp>
#include <MixTransportElement.hpp>
#include <UnsaturatedFlow.hpp>
#include <MaterialUnsaturatedFlow.hpp>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 

Detailed Description

Implementation of unsaturated flow problem.

Definition in file unsaturated_transport.cpp.

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)
Examples
unsaturated_transport.cpp.

Definition at line 46 of file unsaturated_transport.cpp.

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

Variable Documentation

◆ help

char help[] = "...\n\n"
static
Examples
unsaturated_transport.cpp.

Definition at line 37 of file unsaturated_transport.cpp.

MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:26
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:292
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:127
help
static char help[]
Definition: unsaturated_transport.cpp:37
MixTransport::MaterialVanGenuchten
Definition: MaterialUnsaturatedFlow.hpp:254
MixTransport::UnsaturatedFlowElement::BcData
Class storing information about boundary condition.
Definition: UnsaturatedFlow.hpp:158
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MixTransport::UnsaturatedFlowElement
Implementation of operators, problem and finite elements for unsaturated flow.
Definition: UnsaturatedFlow.hpp:115
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:604
param_file
std::string param_file
Definition: DefaultParams.hpp:29
MixTransport::RegisterMaterials
Definition: MaterialUnsaturatedFlow.hpp:350
MixTransport::CommonMaterialData
Definition: MaterialUnsaturatedFlow.hpp:13
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:60
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
StdRDOperators::order
const int order
Definition: rd_stdOperators.hpp:28
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:77
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:35
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
_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:91
default_options
const std::string default_options
Definition: DefaultParams.hpp:1
BLOCKSET
@ BLOCKSET
Definition: definitions.h:217
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:34
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:123
_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:54
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:193
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1129
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:552
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
Definition: UnknownInterface.hpp:130
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:128
MoFEM::MeshsetsManager::setMeshsetFromFile
MoFEMErrorCode setMeshsetFromFile(const string file_name, const bool clean_file_options=true)
add blocksets reading config file