v0.14.0
Loading...
Searching...
No Matches
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
9
12
13using namespace boost::numeric;
14using namespace MoFEM;
15using namespace std;
16
17static char help[] = "Navier-Stokes Example\n";
18
20
21//! [Example Navier Stokes]
23
25 commonData = boost::make_shared<NavierStokesElement::CommonData>(m_field);
26 }
28 CHKERR SNESDestroy(&snes);
29 CHKERR DMDestroy(&dM);
30 }
31
33
34private:
36
37 PetscBool isPartitioned = PETSC_FALSE;
38
39 int orderVelocity = 2; // default approximation orderVelocity
40 int orderPressure = 1; // default approximation orderPressure
41 int numHoLevels = 0;
42 PetscBool isDiscontPressure = PETSC_FALSE;
43 PetscBool isHoGeometry = PETSC_TRUE;
44
45 double pressureScale = 1.0;
46 double lengthScale = 1.0;
47 double velocityScale = 1.0;
48 double reNumber = 1.0;
49 double density;
50 double viscosity;
51 PetscBool isStokesFlow = PETSC_FALSE;
52
53 int numSteps = 1;
54 double lambdaStep = 1.0;
55 double lambda = 0.0;
56 int step;
57
60
61 DM dM;
62 SNES snes;
63
64 boost::shared_ptr<NavierStokesElement::CommonData> commonData;
65
69
70 boost::shared_ptr<VolumeElementForcesAndSourcesCore> feLhsPtr;
71 boost::shared_ptr<VolumeElementForcesAndSourcesCore> feRhsPtr;
72
73 boost::shared_ptr<FaceElementForcesAndSourcesCore> feDragPtr;
74 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feDragSidePtr;
75
76 boost::shared_ptr<DirichletDisplacementBc> dirichletBcPtr;
77 boost::ptr_map<std::string, NeumannForcesSurface> neumannForces;
78
79 boost::shared_ptr<PostProcVol> fePostProcPtr;
80 boost::shared_ptr<PostProcFace> fePostProcDragPtr;
81
92};
93//! [Example Navier Stokes]
94
95//! [Run problem]
109//! [Run problem]
110
111//! [Read input]
114 char mesh_file_name[255];
115 PetscBool flg_mesh_file;
116
117 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "NAVIER_STOKES problem",
118 "none");
119
120 CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
121 mesh_file_name, 255, &flg_mesh_file);
122 CHKERR PetscOptionsBool("-my_is_partitioned", "is partitioned?", "",
123 isPartitioned, &isPartitioned, PETSC_NULL);
124
125 CHKERR PetscOptionsInt("-my_order_u", "approximation orderVelocity", "",
126 orderVelocity, &orderVelocity, PETSC_NULL);
127 CHKERR PetscOptionsInt("-my_order_p", "approximation orderPressure", "",
128 orderPressure, &orderPressure, PETSC_NULL);
129 CHKERR PetscOptionsInt("-my_num_ho_levels",
130 "number of higher order boundary levels", "",
131 numHoLevels, &numHoLevels, PETSC_NULL);
132 CHKERR PetscOptionsBool("-my_discont_pressure", "discontinuous pressure", "",
134 CHKERR PetscOptionsBool("-my_ho_geometry", "use second order geometry", "",
135 isHoGeometry, &isHoGeometry, PETSC_NULL);
136
137 CHKERR PetscOptionsScalar("-my_length_scale", "length scale", "", lengthScale,
138 &lengthScale, PETSC_NULL);
139 CHKERR PetscOptionsScalar("-my_velocity_scale", "velocity scale", "",
140 velocityScale, &velocityScale, PETSC_NULL);
141 CHKERR PetscOptionsBool("-my_stokes_flow", "stokes flow", "", isStokesFlow,
142 &isStokesFlow, PETSC_NULL);
143
144 CHKERR PetscOptionsInt("-my_step_num", "number of steps", "", numSteps,
145 &numSteps, PETSC_NULL);
146
147 ierr = PetscOptionsEnd();
148
149 if (flg_mesh_file != PETSC_TRUE) {
150 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
151 }
152
153 // Read whole mesh or part of it if partitioned
154 if (isPartitioned == PETSC_TRUE) {
155 // This is a case of distributed mesh and algebra. In that case each
156 // processor keeps only part of the problem.
157 const char *option;
158 option = "PARALLEL=READ_PART;"
159 "PARALLEL_RESOLVE_SHARED_ENTS;"
160 "PARTITION=PARALLEL_PARTITION;";
161 CHKERR mField.get_moab().load_file(mesh_file_name, 0, option);
162 } else {
163 // In this case, we have distributed algebra, i.e. assembly of vectors and
164 // matrices is in parallel, but whole mesh is stored on all processors.
165 // snes and matrix scales well, however problem set-up of problem is
166 // not fully parallel.
167 const char *option;
168 option = "";
169 CHKERR mField.get_moab().load_file(mesh_file_name, 0, option);
170 }
172
173 bitLevel.set(0);
174 CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(0, 3,
175 bitLevel);
176
178}
179//! [Read input]
180
181//! [Find blocks]
184
186 if (bit->getName().compare(0, 9, "MAT_FLUID") == 0) {
187 const int id = bit->getMeshsetId();
188 CHKERR mField.get_moab().get_entities_by_dimension(
189 bit->getMeshset(), 3, commonData->setOfBlocksData[id].eNts, true);
190
191 std::vector<double> attributes;
192 bit->getAttributes(attributes);
193 if (attributes.size() < 2) {
194 SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
195 "should be at least 2 attributes but is %d",
196 attributes.size());
197 }
198 commonData->setOfBlocksData[id].iD = id;
199 commonData->setOfBlocksData[id].fluidViscosity = attributes[0];
200 commonData->setOfBlocksData[id].fluidDensity = attributes[1];
201
202 viscosity = commonData->setOfBlocksData[id].fluidViscosity;
203 density = commonData->setOfBlocksData[id].fluidDensity;
204 }
205 }
206
208 if (bit->getName().compare(0, 9, "INT_SOLID") == 0) {
209 Range tets, tet;
210 const int id = bit->getMeshsetId();
211 CHKERR mField.get_moab().get_entities_by_type(
212 bit->getMeshset(), MBTRI, commonData->setOfFacesData[id].eNts, true);
213 solidFaces.merge(commonData->setOfFacesData[id].eNts);
214 CHKERR mField.get_moab().get_adjacencies(
215 commonData->setOfFacesData[id].eNts, 3, true, tets,
216 moab::Interface::UNION);
217 tet = Range(tets.front(), tets.front());
218 for (auto &bit : commonData->setOfBlocksData) {
219 if (bit.second.eNts.contains(tet)) {
220 commonData->setOfFacesData[id].fluidViscosity =
221 bit.second.fluidViscosity;
222 commonData->setOfFacesData[id].fluidDensity = bit.second.fluidDensity;
223 commonData->setOfFacesData[id].iD = id;
224 break;
225 }
226 }
227 if (commonData->setOfFacesData[id].fluidViscosity < 0) {
228 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
229 "Cannot find a fluid block adjacent to a given solid face");
230 }
231 }
232 }
233
235}
236//! [Find blocks]
237
238//! [Setup fields]
241
243 if (isDiscontPressure) {
245 } else {
247 }
248 CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
249 3);
250
251 CHKERR mField.add_ents_to_field_by_dim(0, 3, "VELOCITY");
252 CHKERR mField.add_ents_to_field_by_dim(0, 3, "PRESSURE");
253 CHKERR mField.add_ents_to_field_by_dim(0, 3, "MESH_NODE_POSITIONS");
254
255 CHKERR mField.set_field_order(0, MBVERTEX, "VELOCITY", 1);
256 CHKERR mField.set_field_order(0, MBEDGE, "VELOCITY", orderVelocity);
257 CHKERR mField.set_field_order(0, MBTRI, "VELOCITY", orderVelocity);
258 CHKERR mField.set_field_order(0, MBTET, "VELOCITY", orderVelocity);
259
260 if (!isDiscontPressure) {
261 CHKERR mField.set_field_order(0, MBVERTEX, "PRESSURE", 1);
262 CHKERR mField.set_field_order(0, MBEDGE, "PRESSURE", orderPressure);
263 CHKERR mField.set_field_order(0, MBTRI, "PRESSURE", orderPressure);
264 }
265 CHKERR mField.set_field_order(0, MBTET, "PRESSURE", orderPressure);
266
267 if (numHoLevels > 0) {
268 std::vector<Range> levels(numHoLevels);
269 Range ents;
270 ents.merge(solidFaces);
271 for (int ll = 0; ll != numHoLevels; ll++) {
272 Range verts;
273 CHKERR mField.get_moab().get_connectivity(ents, verts, true);
274 for (auto d : {1, 2, 3}) {
275 CHKERR mField.get_moab().get_adjacencies(verts, d, false, ents,
276 moab::Interface::UNION);
277 }
278 levels[ll] = subtract(ents, ents.subset_by_type(MBVERTEX));
279 }
280 for (int ll = numHoLevels - 1; ll >= 1; ll--) {
281 levels[ll] = subtract(levels[ll], levels[ll - 1]);
282 }
283
284 int add_order = 1;
285 for (int ll = numHoLevels - 1; ll >= 0; ll--) {
286 if (isPartitioned)
287 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
288 levels[ll]);
289
290 CHKERR mField.set_field_order(levels[ll], "VELOCITY",
291 orderVelocity + add_order);
293 CHKERR mField.set_field_order(levels[ll], "PRESSURE",
294 orderPressure + add_order);
295 else
296 CHKERR mField.set_field_order(levels[ll].subset_by_type(MBTET),
297 "PRESSURE", orderPressure + add_order);
298 ++add_order;
299 }
300 }
301
302 CHKERR mField.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
303 // Set 2nd order of approximation of geometry
304 if (isHoGeometry) {
305 Range ents, edges;
306 CHKERR mField.get_moab().get_entities_by_dimension(0, 3, ents);
307 CHKERR mField.get_moab().get_adjacencies(ents, 1, false, edges,
308 moab::Interface::UNION);
309 if (isPartitioned)
310 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(edges);
311 CHKERR mField.set_field_order(edges, "MESH_NODE_POSITIONS", 2);
312 }
313
314 if (isPartitioned) {
315 CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
316 "VELOCITY");
317 CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
318 "PRESSURE");
319 CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
320 "MESH_NODE_POSITIONS");
321 }
322
324
325 Projection10NodeCoordsOnField ent_method_material(mField,
326 "MESH_NODE_POSITIONS");
327 CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
328
330}
331//! [Setup fields]
332
333//! [Define finite elements]
336
337 // Add finite element (this defines element, declaration comes later)
338 CHKERR NavierStokesElement::addElement(mField, "NAVIER_STOKES", "VELOCITY",
339 "PRESSURE", "MESH_NODE_POSITIONS");
340
341 CHKERR NavierStokesElement::addElement(mField, "DRAG", "VELOCITY", "PRESSURE",
342 "MESH_NODE_POSITIONS", 2, &solidFaces);
343
344 // setup elements for loading
346
347 // build finite elements
349 // build adjacencies between elements and degrees of freedom
351
353}
354//! [Define finite elements]
355
356//! [Setup discrete manager]
359 DMType dm_name = "DM_NAVIER_STOKES";
360 // Register DM problem
361 CHKERR DMRegister_MoFEM(dm_name);
362 CHKERR DMCreate(PETSC_COMM_WORLD, &dM);
363 CHKERR DMSetType(dM, dm_name);
364 // Create DM instance
366 // Configure DM form line command options s
367 CHKERR DMSetFromOptions(dM);
368 // Add elements to dM (only one here)
369 CHKERR DMMoFEMAddElement(dM, "NAVIER_STOKES");
370 CHKERR DMMoFEMAddElement(dM, "DRAG");
371 if (mField.check_finite_element("PRESSURE_FE"))
372 CHKERR DMMoFEMAddElement(dM, "PRESSURE_FE");
374 // setup the DM
375 CHKERR DMSetUp(dM);
377}
378//! [Setup discrete manager]
379
380//! [Setup algebraic structures]
383
387
388 CHKERR VecZeroEntries(F);
389 CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
390 CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
391
392 CHKERR VecZeroEntries(D);
393 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
394 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
395
396 CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
397 CHKERR MatZeroEntries(Aij);
398
400}
401//! [Setup algebraic structures]
402
403//! [Setup element instances]
406
407 feLhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
409 feRhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
411
412 feLhsPtr->getRuleHook = NavierStokesElement::VolRule();
413 feRhsPtr->getRuleHook = NavierStokesElement::VolRule();
414 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *feLhsPtr, true, false, false,
415 true);
416 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *feRhsPtr, true, false, false,
417 true);
418
419 feDragPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
421 feDragSidePtr = boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide>(
423
425 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *feDragSidePtr, true, false, false,
426 true);
427
428 if (isStokesFlow) {
430 feRhsPtr, feLhsPtr, "VELOCITY", "PRESSURE", commonData);
431 } else {
433 feRhsPtr, feLhsPtr, "VELOCITY", "PRESSURE", commonData);
434 }
435
437 "NAVIER_STOKES", "VELOCITY",
438 "PRESSURE", commonData);
439
440 dirichletBcPtr = boost::shared_ptr<DirichletDisplacementBc>(
441 new DirichletDisplacementBc(mField, "VELOCITY", Aij, D, F));
442 dirichletBcPtr->methodsOp.push_back(new NavierStokesElement::LoadScale());
443 // dirichletBcPtr->snes_ctx = FEMethod::CTX_SNESNONE;
444 dirichletBcPtr->snes_x = D;
445
446 // Assemble pressure and traction forces
448 NULL, "VELOCITY");
449
450 // for postprocessing:
451 fePostProcPtr = boost::make_shared<PostProcVol>(mField);
452 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *fePostProcPtr, true, false, false,
453 true);
454
455 auto v_ptr = boost::make_shared<MatrixDouble>();
456 auto grad_ptr = boost::make_shared<MatrixDouble>();
457 auto pos_ptr = boost::make_shared<MatrixDouble>();
458 auto p_ptr = boost::make_shared<VectorDouble>();
459
460 fePostProcPtr->getOpPtrVector().push_back(
461 new OpCalculateVectorFieldValues<3>("VELOCITY", v_ptr));
462 fePostProcPtr->getOpPtrVector().push_back(
463 new OpCalculateVectorFieldGradient<3, 3>("VELOCITY", grad_ptr));
464 fePostProcPtr->getOpPtrVector().push_back(
465 new OpCalculateScalarFieldValues("PRESSURE", p_ptr));
466 fePostProcPtr->getOpPtrVector().push_back(
467 new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS", pos_ptr));
468
470
471 fePostProcPtr->getOpPtrVector().push_back(
472
473 new OpPPMap(
474
475 fePostProcPtr->getPostProcMesh(), fePostProcPtr->getMapGaussPts(),
476
477 {{"PRESSURE", p_ptr}},
478
479 {{"VELOCITY", v_ptr}, {"MESH_NODE_POSITIONS", pos_ptr}},
480
481 {{"VELOCITY_GRAD", grad_ptr}},
482
483 {}
484
485 )
486
487 );
488
489 fePostProcDragPtr = boost::make_shared<PostProcFace>(mField);
490 fePostProcDragPtr->getOpPtrVector().push_back(
491 new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS", pos_ptr));
492 fePostProcDragPtr->getOpPtrVector().push_back(
493
494 new OpPPMap(
495
496 fePostProcDragPtr->getPostProcMesh(),
497 fePostProcDragPtr->getMapGaussPts(),
498
499 {},
500
501 {{"MESH_NODE_POSITIONS", pos_ptr}},
502
503 {},
504
505 {}
506
507 )
508
509 );
510
512 fePostProcDragPtr, feDragSidePtr, "NAVIER_STOKES", "VELOCITY", "PRESSURE",
513 commonData);
514
516}
517//! [Setup element instances]
518
519//! [Setup SNES]
522
523 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
524 neumannForces.begin();
525 for (; mit != neumannForces.end(); mit++) {
526 CHKERR DMMoFEMSNESSetFunction(dM, mit->first.c_str(),
527 &mit->second->getLoopFe(), NULL, NULL);
528 }
529
530 boost::shared_ptr<FEMethod> null_fe;
531 CHKERR DMMoFEMSNESSetFunction(dM, "NAVIER_STOKES", feRhsPtr, null_fe,
532 null_fe);
535
537 null_fe);
538 CHKERR DMMoFEMSNESSetJacobian(dM, "NAVIER_STOKES", feLhsPtr, null_fe,
539 null_fe);
542
543 SnesCtx *snes_ctx;
544 // create snes nonlinear solver
545 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
546 CHKERR DMMoFEMGetSnesCtx(dM, &snes_ctx);
547 CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
548 CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, snes_ctx);
549 CHKERR SNESSetFromOptions(snes);
550
552}
553//! [Setup SNES]
554
555//! [Solve problem]
558
559 auto scale_problem = [&](double U, double L, double P) {
561 CHKERR mField.getInterface<FieldBlas>()->fieldScale(L,
562 "MESH_NODE_POSITIONS");
563
564 ProjectionFieldOn10NodeTet ent_method_on_10nodeTet(
565 mField, "MESH_NODE_POSITIONS", true, true);
566 CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_on_10nodeTet);
567 ent_method_on_10nodeTet.setNodes = false;
568 CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_on_10nodeTet);
569
570 CHKERR mField.getInterface<FieldBlas>()->fieldScale(U, "VELOCITY");
571 CHKERR mField.getInterface<FieldBlas>()->fieldScale(P, "PRESSURE");
573 };
574
577 CHKERR scale_problem(1.0 / velocityScale, 1.0 / lengthScale,
578 1.0 / pressureScale);
579
580 step = 0;
581
582 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Viscosity: %6.4e\n", viscosity);
583 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Density: %6.4e\n", density);
584 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Length scale: %6.4e\n", lengthScale);
585 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Velocity scale: %6.4e\n",
587 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Pressure scale: %6.4e\n",
589 if (isStokesFlow) {
590 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Re number: 0 (Stokes flow)\n");
591 } else {
592 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Re number: %6.4e\n",
594 }
595
596 lambdaStep = 1.0 / numSteps;
597
598 while (lambda < 1.0 - 1e-12) {
599
601
602 if (isStokesFlow) {
603 reNumber = 0.0;
604 for (auto &bit : commonData->setOfBlocksData) {
605 bit.second.inertiaCoef = 0.0;
606 bit.second.viscousCoef = 1.0;
607 }
608 } else {
610 for (auto &bit : commonData->setOfBlocksData) {
611 bit.second.inertiaCoef = reNumber;
612 bit.second.viscousCoef = 1.0;
613 }
614 }
615
616 CHKERR PetscPrintf(
617 PETSC_COMM_WORLD,
618 "Step: %d | Lambda: %6.4e | Inc: %6.4e | Re number: %6.4e \n", step,
620
622
623 CHKERR VecAssemblyBegin(D);
624 CHKERR VecAssemblyEnd(D);
625 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
626 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
627
628 CHKERR SNESSolve(snes, PETSC_NULL, D);
629
630 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
631 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
632 CHKERR DMoFEMMeshToGlobalVector(dM, D, INSERT_VALUES, SCATTER_REVERSE);
633
635
637
638 CHKERR scale_problem(1.0 / velocityScale, 1.0 / lengthScale,
639 1.0 / pressureScale);
640
641 step++;
642 }
643
645}
646//! [Solve problem]
647
648//! [Post process]
651
652 string out_file_name;
653
655 {
656 std::ostringstream stm;
657 stm << "out_" << step << ".h5m";
658 out_file_name = stm.str();
659 CHKERR PetscPrintf(PETSC_COMM_WORLD, "Write file %s\n",
660 out_file_name.c_str());
661 CHKERR fePostProcPtr->writeFile(out_file_name.c_str());
662 }
663
664 CHKERR VecZeroEntries(commonData->pressureDragForceVec);
665 CHKERR VecZeroEntries(commonData->shearDragForceVec);
666 CHKERR VecZeroEntries(commonData->totalDragForceVec);
668
669 auto get_vec_data = [&](auto vec, std::array<double, 3> &data) {
671 CHKERR VecAssemblyBegin(vec);
672 CHKERR VecAssemblyEnd(vec);
673 const double *array;
674 CHKERR VecGetArrayRead(vec, &array);
675 if (mField.get_comm_rank() == 0) {
676 for (int i : {0, 1, 2})
677 data[i] = array[i];
678 }
679 CHKERR VecRestoreArrayRead(vec, &array);
681 };
682
683 std::array<double, 3> pressure_drag;
684 std::array<double, 3> shear_drag;
685 std::array<double, 3> total_drag;
686
687 CHKERR get_vec_data(commonData->pressureDragForceVec, pressure_drag);
688 CHKERR get_vec_data(commonData->shearDragForceVec, shear_drag);
689 CHKERR get_vec_data(commonData->totalDragForceVec, total_drag);
690
691 if (mField.get_comm_rank() == 0) {
692 MOFEM_LOG_C("SELF", Sev::inform,
693 "Pressure drag force: [ %6.4e, %6.4e, %6.4e ]",
694 pressure_drag[0], pressure_drag[1], pressure_drag[2]);
695 MOFEM_LOG_C("SELF", Sev::inform,
696 "Shear drag force: [ %6.4e, %6.4e, %6.4e ]", shear_drag[0],
697 shear_drag[1], shear_drag[2]);
698 MOFEM_LOG_C("SELF", Sev::inform,
699 "Total drag force: [ %6.4e, %6.4e, %6.4e ]", total_drag[0],
700 total_drag[1], total_drag[2]);
701 }
702
704 {
705 std::ostringstream stm;
706 stm << "out_drag_" << step << ".h5m";
707 out_file_name = stm.str();
708 CHKERR PetscPrintf(PETSC_COMM_WORLD, "out file %s\n",
709 out_file_name.c_str());
711 }
713}
714//! [Post process]
715
716//! [Main function]
717int main(int argc, char *argv[]) {
718
719 const char param_file[] = "param_file.petsc";
720 // Initialise MoFEM
722
723 try {
724
725 // Create mesh database
726 moab::Core mb_instance; // create database
727 moab::Interface &moab = mb_instance; // create interface to database
728
729 // Create MoFEM database and link it to MoAB
730 MoFEM::Core core(moab);
731 MoFEM::Interface &m_field = core;
732
733 NavierStokesExample ex(m_field);
734 CHKERR ex.runProblem();
735 }
737
738 // finish work cleaning memory, getting statistics, etc
740 return 0;
741}
742//! [Main function]
#define DM_NO_ELEMENT
Definition DMMoFEM.hpp:10
std::string param_file
#define MOFEM_LOG_C(channel, severity, format,...)
int main()
PostProcEleByDim< SPACE_DIM >::PostProcFace PostProcFace
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ L2
field with C-1 continuity
Definition definitions.h:88
@ H1
continuous field
Definition definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ BLOCKSET
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition DMMoFEM.cpp:1109
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:483
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 DMMoFEM.cpp:704
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 DMMoFEM.cpp:118
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 DMMoFEM.cpp:745
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:572
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1003
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:988
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition DMMoFEM.cpp:1080
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition DMMoFEM.cpp:521
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:532
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 out_file_name[255]
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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 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:139
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:27
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
static char help[]
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Set Dirichlet boundary conditions on displacements.
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: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.
Basic algebra on fields.
Definition FieldBlas.hpp:21
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Projection of edge entities with one mid-node on hierarchical basis.
intrusive_ptr for managing petsc objects
Interface for nonlinear (SNES) solver.
Definition SnesCtx.hpp:13
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
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 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 setPostProcDragOperators(boost::shared_ptr< PostProcFace > 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 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< PostProcFace > fePostProcDragPtr
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
MoFEMErrorCode setupAlgebraicStructures()
[Setup discrete manager]
NavierStokesExample(MoFEM::Interface &m_field)
boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > feDragSidePtr
MoFEMErrorCode findBlocks()
[Read input]
boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhsPtr
boost::ptr_map< std::string, NeumannForcesSurface > neumannForces
MoFEM::Interface & mField
boost::shared_ptr< PostProcVol > fePostProcPtr
MoFEMErrorCode solveProblem()
[Setup SNES]
MoFEMErrorCode postProcess()
[Solve problem]
MoFEMErrorCode readInput()
[Run problem]
SmartPetscObj< Vec > D
SmartPetscObj< Mat > Aij