v0.14.0
Loading...
Searching...
No Matches
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
13#include <UnsaturatedFlow.hpp>
15
16#ifndef WITH_ADOL_C
17#error "MoFEM need to be compiled with ADOL-C"
18#endif
19
20using namespace MoFEM;
21using namespace MixTransport;
22static char help[] = "...\n\n";
23
27
28map<std::string, CommonMaterialData::RegisterHook>
30
31int 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
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) {
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");
86 ierr = PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
87 mesh_file_name, 255, &flg_mesh_file_name);
89 ierr = PetscOptionsString(
90 "-configure", "material and bc configuration file name", "",
91 "unsaturated.cfg", conf_file_name, 255, &flg_conf_file_name);
93 ierr = PetscOptionsInt("-my_order", "default approximation order", "",
94 order, &order, PETSC_NULL);
96 ierr = PetscOptionsEnd();
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
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 }
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 }
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);
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
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>(
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
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>(
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 }
303
305
306 return 0;
307}
const std::string default_options
std::string param_file
Mix implementation of transport element.
Mix implementation of transport element.
Mix implementation of transport element.
int main()
Definition: adol-c_atom.cpp:46
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
#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:76
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
void printMatParameters(const int id, const std::string &prefix) const
static double ePsilon0
Regularization parameter.
static double scaleZ
Scale z direction.
static double ePsilon1
Regularization parameter.
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)
static map< std::string, CommonMaterialData::RegisterHook > mapOfRegistredMaterials
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:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
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.
static char help[]