v0.13.1
navier_stokes.cpp
Go to the documentation of this file.
1/** \file navier_stokes.cpp
2 * \example navier_stokes.cpp
3 *
4 * Example of viscous fluid flow problem
5 *
6 */
7
8/* MoFEM is free software: you can redistribute it and/or modify it under
9 * the terms of the GNU Lesser General Public License as published by the
10 * Free Software Foundation, either version 3 of the License, or (at your
11 * option) any later version.
12 *
13 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 * License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>.
20 * */
21
22#include <BasicFiniteElements.hpp>
23
24using namespace boost::numeric;
25using namespace MoFEM;
26using namespace std;
27
28static char help[] = "Navier-Stokes Example\n";
29
31
32//! [Example Navier Stokes]
34
36 commonData = boost::make_shared<NavierStokesElement::CommonData>(m_field);
37 }
39 CHKERR SNESDestroy(&snes);
40 CHKERR DMDestroy(&dM);
41 }
42
44
45private:
47
48 PetscBool isPartitioned = PETSC_FALSE;
49
50 int orderVelocity = 2; // default approximation orderVelocity
51 int orderPressure = 1; // default approximation orderPressure
52 int numHoLevels = 0;
53 PetscBool isDiscontPressure = PETSC_FALSE;
54 PetscBool isHoGeometry = PETSC_TRUE;
55
56 double pressureScale = 1.0;
57 double lengthScale = 1.0;
58 double velocityScale = 1.0;
59 double reNumber = 1.0;
60 double density;
61 double viscosity;
62 PetscBool isStokesFlow = PETSC_FALSE;
63
64 int numSteps = 1;
65 double lambdaStep = 1.0;
66 double lambda = 0.0;
67 int step;
68
71
72 DM dM;
73 SNES snes;
74
75 boost::shared_ptr<NavierStokesElement::CommonData> commonData;
76
80
81 boost::shared_ptr<VolumeElementForcesAndSourcesCore> feLhsPtr;
82 boost::shared_ptr<VolumeElementForcesAndSourcesCore> feRhsPtr;
83
84 boost::shared_ptr<FaceElementForcesAndSourcesCore> feDragPtr;
85 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feDragSidePtr;
86
87 boost::shared_ptr<DirichletDisplacementBc> dirichletBcPtr;
88 boost::ptr_map<std::string, NeumannForcesSurface> neumannForces;
89
90 boost::shared_ptr<PostProcVolumeOnRefinedMesh> fePostProcPtr;
91 boost::shared_ptr<PostProcFaceOnRefinedMesh> fePostProcDragPtr;
92
103};
104//! [Example Navier Stokes]
105
106//! [Run problem]
119}
120//! [Run problem]
121
122//! [Read input]
125 char mesh_file_name[255];
126 PetscBool flg_mesh_file;
127
128 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "NAVIER_STOKES problem",
129 "none");
130
131 CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
132 mesh_file_name, 255, &flg_mesh_file);
133 CHKERR PetscOptionsBool("-my_is_partitioned", "is partitioned?", "",
134 isPartitioned, &isPartitioned, PETSC_NULL);
135
136 CHKERR PetscOptionsInt("-my_order_u", "approximation orderVelocity", "",
137 orderVelocity, &orderVelocity, PETSC_NULL);
138 CHKERR PetscOptionsInt("-my_order_p", "approximation orderPressure", "",
139 orderPressure, &orderPressure, PETSC_NULL);
140 CHKERR PetscOptionsInt("-my_num_ho_levels",
141 "number of higher order boundary levels", "",
142 numHoLevels, &numHoLevels, PETSC_NULL);
143 CHKERR PetscOptionsBool("-my_discont_pressure", "discontinuous pressure", "",
145 CHKERR PetscOptionsBool("-my_ho_geometry", "use second order geometry", "",
146 isHoGeometry, &isHoGeometry, PETSC_NULL);
147
148 CHKERR PetscOptionsScalar("-my_length_scale", "length scale", "", lengthScale,
149 &lengthScale, PETSC_NULL);
150 CHKERR PetscOptionsScalar("-my_velocity_scale", "velocity scale", "",
151 velocityScale, &velocityScale, PETSC_NULL);
152 CHKERR PetscOptionsBool("-my_stokes_flow", "stokes flow", "", isStokesFlow,
153 &isStokesFlow, PETSC_NULL);
154
155 CHKERR PetscOptionsInt("-my_step_num", "number of steps", "", numSteps,
156 &numSteps, PETSC_NULL);
157
158 ierr = PetscOptionsEnd();
159
160 if (flg_mesh_file != PETSC_TRUE) {
161 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
162 }
163
164 // Read whole mesh or part of it if partitioned
165 if (isPartitioned == PETSC_TRUE) {
166 // This is a case of distributed mesh and algebra. In that case each
167 // processor keeps only part of the problem.
168 const char *option;
169 option = "PARALLEL=READ_PART;"
170 "PARALLEL_RESOLVE_SHARED_ENTS;"
171 "PARTITION=PARALLEL_PARTITION;";
172 CHKERR mField.get_moab().load_file(mesh_file_name, 0, option);
173 } else {
174 // In this case, we have distributed algebra, i.e. assembly of vectors and
175 // matrices is in parallel, but whole mesh is stored on all processors.
176 // snes and matrix scales well, however problem set-up of problem is
177 // not fully parallel.
178 const char *option;
179 option = "";
180 CHKERR mField.get_moab().load_file(mesh_file_name, 0, option);
181 }
183
184 bitLevel.set(0);
185 CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(0, 3,
186 bitLevel);
187
189}
190//! [Read input]
191
192//! [Find blocks]
195
197 if (bit->getName().compare(0, 9, "MAT_FLUID") == 0) {
198 const int id = bit->getMeshsetId();
199 CHKERR mField.get_moab().get_entities_by_dimension(
200 bit->getMeshset(), 3, commonData->setOfBlocksData[id].eNts, true);
201
202 std::vector<double> attributes;
203 bit->getAttributes(attributes);
204 if (attributes.size() < 2) {
205 SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
206 "should be at least 2 attributes but is %d",
207 attributes.size());
208 }
209 commonData->setOfBlocksData[id].iD = id;
210 commonData->setOfBlocksData[id].fluidViscosity = attributes[0];
211 commonData->setOfBlocksData[id].fluidDensity = attributes[1];
212
213 viscosity = commonData->setOfBlocksData[id].fluidViscosity;
214 density = commonData->setOfBlocksData[id].fluidDensity;
215 }
216 }
217
219 if (bit->getName().compare(0, 9, "INT_SOLID") == 0) {
220 Range tets, tet;
221 const int id = bit->getMeshsetId();
222 CHKERR mField.get_moab().get_entities_by_type(
223 bit->getMeshset(), MBTRI, commonData->setOfFacesData[id].eNts, true);
224 solidFaces.merge(commonData->setOfFacesData[id].eNts);
225 CHKERR mField.get_moab().get_adjacencies(
226 commonData->setOfFacesData[id].eNts, 3, true, tets,
227 moab::Interface::UNION);
228 tet = Range(tets.front(), tets.front());
229 for (auto &bit : commonData->setOfBlocksData) {
230 if (bit.second.eNts.contains(tet)) {
231 commonData->setOfFacesData[id].fluidViscosity =
232 bit.second.fluidViscosity;
233 commonData->setOfFacesData[id].fluidDensity = bit.second.fluidDensity;
234 commonData->setOfFacesData[id].iD = id;
235 break;
236 }
237 }
238 if (commonData->setOfFacesData[id].fluidViscosity < 0) {
239 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
240 "Cannot find a fluid block adjacent to a given solid face");
241 }
242 }
243 }
244
246}
247//! [Find blocks]
248
249//! [Setup fields]
252
254 if (isDiscontPressure) {
256 } else {
258 }
259 CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
260 3);
261
262 CHKERR mField.add_ents_to_field_by_dim(0, 3, "VELOCITY");
263 CHKERR mField.add_ents_to_field_by_dim(0, 3, "PRESSURE");
264 CHKERR mField.add_ents_to_field_by_dim(0, 3, "MESH_NODE_POSITIONS");
265
266 CHKERR mField.set_field_order(0, MBVERTEX, "VELOCITY", 1);
267 CHKERR mField.set_field_order(0, MBEDGE, "VELOCITY", orderVelocity);
268 CHKERR mField.set_field_order(0, MBTRI, "VELOCITY", orderVelocity);
269 CHKERR mField.set_field_order(0, MBTET, "VELOCITY", orderVelocity);
270
271 if (!isDiscontPressure) {
272 CHKERR mField.set_field_order(0, MBVERTEX, "PRESSURE", 1);
273 CHKERR mField.set_field_order(0, MBEDGE, "PRESSURE", orderPressure);
274 CHKERR mField.set_field_order(0, MBTRI, "PRESSURE", orderPressure);
275 }
276 CHKERR mField.set_field_order(0, MBTET, "PRESSURE", orderPressure);
277
278 if (numHoLevels > 0) {
279 std::vector<Range> levels(numHoLevels);
280 Range ents;
281 ents.merge(solidFaces);
282 for (int ll = 0; ll != numHoLevels; ll++) {
283 Range verts;
284 CHKERR mField.get_moab().get_connectivity(ents, verts, true);
285 for (auto d : {1, 2, 3}) {
286 CHKERR mField.get_moab().get_adjacencies(verts, d, false, ents,
287 moab::Interface::UNION);
288 }
289 levels[ll] = subtract(ents, ents.subset_by_type(MBVERTEX));
290 }
291 for (int ll = numHoLevels - 1; ll >= 1; ll--) {
292 levels[ll] = subtract(levels[ll], levels[ll - 1]);
293 }
294
295 int add_order = 1;
296 for (int ll = numHoLevels - 1; ll >= 0; ll--) {
297 if (isPartitioned)
298 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
299 levels[ll]);
300
301 CHKERR mField.set_field_order(levels[ll], "VELOCITY",
302 orderVelocity + add_order);
304 CHKERR mField.set_field_order(levels[ll], "PRESSURE",
305 orderPressure + add_order);
306 else
307 CHKERR mField.set_field_order(levels[ll].subset_by_type(MBTET),
308 "PRESSURE", orderPressure + add_order);
309 ++add_order;
310 }
311 }
312
313 CHKERR mField.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
314 // Set 2nd order of approximation of geometry
315 if (isHoGeometry) {
316 Range ents, edges;
317 CHKERR mField.get_moab().get_entities_by_dimension(0, 3, ents);
318 CHKERR mField.get_moab().get_adjacencies(ents, 1, false, edges,
319 moab::Interface::UNION);
320 if (isPartitioned)
321 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(edges);
322 CHKERR mField.set_field_order(edges, "MESH_NODE_POSITIONS", 2);
323 }
324
325 if (isPartitioned) {
326 CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
327 "VELOCITY");
328 CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
329 "PRESSURE");
330 CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
331 "MESH_NODE_POSITIONS");
332 }
333
335
336 Projection10NodeCoordsOnField ent_method_material(mField,
337 "MESH_NODE_POSITIONS");
338 CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
339
341}
342//! [Setup fields]
343
344//! [Define finite elements]
347
348 // Add finite element (this defines element, declaration comes later)
349 CHKERR NavierStokesElement::addElement(mField, "NAVIER_STOKES", "VELOCITY",
350 "PRESSURE", "MESH_NODE_POSITIONS");
351
352 CHKERR NavierStokesElement::addElement(mField, "DRAG", "VELOCITY", "PRESSURE",
353 "MESH_NODE_POSITIONS", 2, &solidFaces);
354
355 // setup elements for loading
357
358 // build finite elements
360 // build adjacencies between elements and degrees of freedom
362
364}
365//! [Define finite elements]
366
367//! [Setup discrete manager]
370 DMType dm_name = "DM_NAVIER_STOKES";
371 // Register DM problem
372 CHKERR DMRegister_MoFEM(dm_name);
373 CHKERR DMCreate(PETSC_COMM_WORLD, &dM);
374 CHKERR DMSetType(dM, dm_name);
375 // Create DM instance
377 // Configure DM form line command options s
378 CHKERR DMSetFromOptions(dM);
379 // Add elements to dM (only one here)
380 CHKERR DMMoFEMAddElement(dM, "NAVIER_STOKES");
381 CHKERR DMMoFEMAddElement(dM, "DRAG");
382 if (mField.check_finite_element("PRESSURE_FE"))
383 CHKERR DMMoFEMAddElement(dM, "PRESSURE_FE");
385 // setup the DM
386 CHKERR DMSetUp(dM);
388}
389//! [Setup discrete manager]
390
391//! [Setup algebraic structures]
394
398
399 CHKERR VecZeroEntries(F);
400 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
401 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
402
403 CHKERR VecZeroEntries(D);
404 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
405 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
406
407 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
408 CHKERR MatZeroEntries(Aij);
409
411}
412//! [Setup algebraic structures]
413
414//! [Setup element instances]
417
418 feLhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
420 feRhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
422
423 feLhsPtr->getRuleHook = NavierStokesElement::VolRule();
424 feRhsPtr->getRuleHook = NavierStokesElement::VolRule();
425 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *feLhsPtr, true, false, false, true);
426 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *feRhsPtr, true, false, false, true);
427
428 feDragPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
430 feDragSidePtr = boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide>(
432
434 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *feDragSidePtr, true, false, false,
435 true);
436
437 if (isStokesFlow) {
439 feRhsPtr, feLhsPtr, "VELOCITY", "PRESSURE", commonData);
440 } else {
442 feRhsPtr, feLhsPtr, "VELOCITY", "PRESSURE", commonData);
443 }
444
446 "NAVIER_STOKES", "VELOCITY",
447 "PRESSURE", commonData);
448
449 dirichletBcPtr = boost::shared_ptr<DirichletDisplacementBc>(
450 new DirichletDisplacementBc(mField, "VELOCITY", Aij, D, F));
451 dirichletBcPtr->methodsOp.push_back(new NavierStokesElement::LoadScale());
452 // dirichletBcPtr->snes_ctx = FEMethod::CTX_SNESNONE;
453 dirichletBcPtr->snes_x = D;
454
455 // Assemble pressure and traction forces
457 NULL, "VELOCITY");
458
459 // for postprocessing:
460 fePostProcPtr = boost::make_shared<PostProcVolumeOnRefinedMesh>(mField);
461 CHKERR fePostProcPtr->generateReferenceElementMesh();
462 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *fePostProcPtr, true, false, false,
463 true);
464 CHKERR fePostProcPtr->addFieldValuesPostProc("VELOCITY");
465 CHKERR fePostProcPtr->addFieldValuesPostProc("PRESSURE");
466 CHKERR fePostProcPtr->addFieldValuesPostProc("MESH_NODE_POSITIONS");
467 CHKERR fePostProcPtr->addFieldValuesGradientPostProc("VELOCITY");
468
469 fePostProcDragPtr = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
470 CHKERR fePostProcDragPtr->generateReferenceElementMesh();
471 CHKERR fePostProcDragPtr->addFieldValuesPostProc("MESH_NODE_POSITIONS");
473 fePostProcDragPtr, feDragSidePtr, "NAVIER_STOKES", "VELOCITY", "PRESSURE",
474 commonData);
475
477}
478//! [Setup element instances]
479
480//! [Setup SNES]
483
484 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
485 neumannForces.begin();
486 for (; mit != neumannForces.end(); mit++) {
487 CHKERR DMMoFEMSNESSetFunction(dM, mit->first.c_str(),
488 &mit->second->getLoopFe(), NULL, NULL);
489 }
490
491 boost::shared_ptr<FEMethod> null_fe;
492 CHKERR DMMoFEMSNESSetFunction(dM, "NAVIER_STOKES", feRhsPtr, null_fe,
493 null_fe);
496
498 null_fe);
499 CHKERR DMMoFEMSNESSetJacobian(dM, "NAVIER_STOKES", feLhsPtr, null_fe,
500 null_fe);
503
504 SnesCtx *snes_ctx;
505 // create snes nonlinear solver
506 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
507 CHKERR DMMoFEMGetSnesCtx(dM, &snes_ctx);
508 CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
509 CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, snes_ctx);
510 CHKERR SNESSetFromOptions(snes);
511
513}
514//! [Setup SNES]
515
516//! [Solve problem]
519
520 auto scale_problem = [&](double U, double L, double P) {
522 CHKERR mField.getInterface<FieldBlas>()->fieldScale(L,
523 "MESH_NODE_POSITIONS");
524
525 ProjectionFieldOn10NodeTet ent_method_on_10nodeTet(
526 mField, "MESH_NODE_POSITIONS", true, true);
527 CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_on_10nodeTet);
528 ent_method_on_10nodeTet.setNodes = false;
529 CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_on_10nodeTet);
530
531 CHKERR mField.getInterface<FieldBlas>()->fieldScale(U, "VELOCITY");
532 CHKERR mField.getInterface<FieldBlas>()->fieldScale(P, "PRESSURE");
534 };
535
538 CHKERR scale_problem(1.0 / velocityScale, 1.0 / lengthScale,
539 1.0 / pressureScale);
540
541 step = 0;
542
543 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Viscosity: %6.4e\n", viscosity);
544 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Density: %6.4e\n", density);
545 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Length scale: %6.4e\n", lengthScale);
546 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Velocity scale: %6.4e\n",
548 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Pressure scale: %6.4e\n",
550 if (isStokesFlow) {
551 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Re number: 0 (Stokes flow)\n");
552 } else {
553 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Re number: %6.4e\n",
555 }
556
557 lambdaStep = 1.0 / numSteps;
558
559 while (lambda < 1.0 - 1e-12) {
560
562
563 if (isStokesFlow) {
564 reNumber = 0.0;
565 for (auto &bit : commonData->setOfBlocksData) {
566 bit.second.inertiaCoef = 0.0;
567 bit.second.viscousCoef = 1.0;
568 }
569 } else {
571 for (auto &bit : commonData->setOfBlocksData) {
572 bit.second.inertiaCoef = reNumber;
573 bit.second.viscousCoef = 1.0;
574 }
575 }
576
577 CHKERR PetscPrintf(
578 PETSC_COMM_WORLD,
579 "Step: %d | Lambda: %6.4e | Inc: %6.4e | Re number: %6.4e \n", step,
581
583
584 CHKERR VecAssemblyBegin(D);
585 CHKERR VecAssemblyEnd(D);
586 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
587 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
588
589 CHKERR SNESSolve(snes, PETSC_NULL, D);
590
591 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
592 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
593 CHKERR DMoFEMMeshToGlobalVector(dM, D, INSERT_VALUES, SCATTER_REVERSE);
594
596
598
599 CHKERR scale_problem(1.0 / velocityScale, 1.0 / lengthScale,
600 1.0 / pressureScale);
601
602 step++;
603 }
604
606}
607//! [Solve problem]
608
609//! [Post process]
612
613 string out_file_name;
614
616 {
617 std::ostringstream stm;
618 stm << "out_" << step << ".h5m";
619 out_file_name = stm.str();
620 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Write file %s\n",
621 out_file_name.c_str());
622 CHKERR fePostProcPtr->postProcMesh.write_file(out_file_name.c_str(), "MOAB",
623 "PARALLEL=WRITE_PART");
624 }
625
626 CHKERR VecZeroEntries(commonData->pressureDragForceVec);
627 CHKERR VecZeroEntries(commonData->shearDragForceVec);
628 CHKERR VecZeroEntries(commonData->totalDragForceVec);
630
631 auto get_vec_data = [&](auto vec, std::array<double, 3> &data) {
633 CHKERR VecAssemblyBegin(vec);
634 CHKERR VecAssemblyEnd(vec);
635 const double *array;
636 CHKERR VecGetArrayRead(vec, &array);
637 if (mField.get_comm_rank() == 0) {
638 for (int i : {0, 1, 2})
639 data[i] = array[i];
640 }
641 CHKERR VecRestoreArrayRead(vec, &array);
643 };
644
645 std::array<double, 3> pressure_drag;
646 std::array<double, 3> shear_drag;
647 std::array<double, 3> total_drag;
648
649 CHKERR get_vec_data(commonData->pressureDragForceVec, pressure_drag);
650 CHKERR get_vec_data(commonData->shearDragForceVec, shear_drag);
651 CHKERR get_vec_data(commonData->totalDragForceVec, total_drag);
652
653 if (mField.get_comm_rank() == 0) {
654 MOFEM_LOG_C("SELF", Sev::inform,
655 "Pressure drag force: [ %6.4e, %6.4e, %6.4e ]",
656 pressure_drag[0], pressure_drag[1], pressure_drag[2]);
657 MOFEM_LOG_C("SELF", Sev::inform,
658 "Shear drag force: [ %6.4e, %6.4e, %6.4e ]", shear_drag[0],
659 shear_drag[1], shear_drag[2]);
660 MOFEM_LOG_C("SELF", Sev::inform,
661 "Total drag force: [ %6.4e, %6.4e, %6.4e ]", total_drag[0],
662 total_drag[1], total_drag[2]);
663 }
664
666 {
667 std::ostringstream stm;
668 stm << "out_drag_" << step << ".h5m";
669 out_file_name = stm.str();
670 CHKERR PetscPrintf(PETSC_COMM_WORLD, "out file %s\n",
671 out_file_name.c_str());
672 CHKERR fePostProcDragPtr->postProcMesh.write_file(
673 out_file_name.c_str(), "MOAB", "PARALLEL=WRITE_PART");
674 }
676}
677//! [Post process]
678
679//! [Main function]
680int main(int argc, char *argv[]) {
681
682 const char param_file[] = "param_file.petsc";
683 // Initialise MoFEM
685
686
687
688 try {
689
690 // Create mesh database
691 moab::Core mb_instance; // create database
692 moab::Interface &moab = mb_instance; // create interface to database
693
694 // Create MoFEM database and link it to MoAB
695 MoFEM::Core core(moab);
696 MoFEM::Interface &m_field = core;
697
698 NavierStokesExample ex(m_field);
699 CHKERR ex.runProblem();
700 }
702
703 // finish work cleaning memory, getting statistics, etc
705 return 0;
706}
707//! [Main function]
static Index< 'L', 3 > L
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:22
std::string param_file
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:314
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ BLOCKSET
Definition: definitions.h:161
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
auto smartCreateDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:965
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1082
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition: DMMMoFEM.cpp:677
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMMoFEM.cpp:130
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:462
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition: DMMMoFEM.cpp:718
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:545
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMMoFEM.cpp:1053
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:494
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:977
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:505
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
#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.
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
char mesh_file_name[255]
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
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 SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition: SnesCtx.cpp:148
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition: SnesCtx.cpp:39
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
STL namespace.
int main(int argc, char *argv[])
[Post process]
static char help[]
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:69
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
Managing BitRefLevels.
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual int get_comm_rank() const =0
virtual MoFEMErrorCode rebuild_database(int verb=DEFAULT_VERBOSITY)=0
Clear database and initialize it once again.
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.
Basic algebra on fields.
Definition: FieldBlas.hpp:31
Projection of edge entities with one mid-node on hierarchical basis.
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:27
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Base volume element used to integrate on skeleton.
Set integration rule to volume elements.
static MoFEMErrorCode setCalcDragOperators(boost::shared_ptr< FaceElementForcesAndSourcesCore > dragFe, boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideDragFe, std::string side_fe_name, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data)
Setting up operators for calculating drag force on the solid surface.
static MoFEMErrorCode setPostProcDragOperators(boost::shared_ptr< PostProcFaceOnRefinedMesh > postProcDragPtr, boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideDragFe, std::string side_fe_name, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data)
Setting up operators for post processing output of drag traction.
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const string element_name, const string velocity_field_name, const string pressure_field_name, const string mesh_field_name, const int dim=3, Range *ents=nullptr)
Setting up elements.
static MoFEMErrorCode setNavierStokesOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhs, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for solving Navier-Stokes equations.
static MoFEMErrorCode setStokesOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhs, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for solving Stokes equations.
[Example Navier Stokes]
MoFEMErrorCode setupElementInstances()
[Setup algebraic structures]
MoFEMErrorCode setupDiscreteManager()
[Define finite elements]
boost::shared_ptr< FaceElementForcesAndSourcesCore > feDragPtr
MoFEMErrorCode defineFiniteElements()
[Setup fields]
boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhsPtr
MoFEMErrorCode setupFields()
[Find blocks]
MoFEMErrorCode setupSNES()
[Setup element instances]
MoFEMErrorCode runProblem()
[Example Navier Stokes]
SmartPetscObj< Vec > F
boost::shared_ptr< NavierStokesElement::CommonData > commonData
boost::shared_ptr< DirichletDisplacementBc > dirichletBcPtr
PetscBool isDiscontPressure
boost::shared_ptr< PostProcVolumeOnRefinedMesh > fePostProcPtr
MoFEMErrorCode setupAlgebraicStructures()
[Setup discrete manager]
NavierStokesExample(MoFEM::Interface &m_field)
boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > feDragSidePtr
boost::shared_ptr< PostProcFaceOnRefinedMesh > fePostProcDragPtr
MoFEMErrorCode findBlocks()
[Read input]
boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhsPtr
boost::ptr_map< std::string, NeumannForcesSurface > neumannForces
MoFEM::Interface & mField
MoFEMErrorCode solveProblem()
[Setup SNES]
MoFEMErrorCode postProcess()
[Solve problem]
MoFEMErrorCode readInput()
[Run problem]
SmartPetscObj< Vec > D
SmartPetscObj< Mat > Aij