v0.13.2
Loading...
Searching...
No Matches
Public Member Functions | Private Member Functions | Private Attributes | List of all members
NavierStokesExample Struct Reference

[Example Navier Stokes] More...

Collaboration diagram for NavierStokesExample:
[legend]

Public Member Functions

 NavierStokesExample (MoFEM::Interface &m_field)
 
 ~NavierStokesExample ()
 
MoFEMErrorCode runProblem ()
 [Example Navier Stokes] More...
 

Private Member Functions

MoFEMErrorCode readInput ()
 [Run problem] More...
 
MoFEMErrorCode findBlocks ()
 [Read input] More...
 
MoFEMErrorCode setupFields ()
 [Find blocks] More...
 
MoFEMErrorCode defineFiniteElements ()
 [Setup fields] More...
 
MoFEMErrorCode setupDiscreteManager ()
 [Define finite elements] More...
 
MoFEMErrorCode setupAlgebraicStructures ()
 [Setup discrete manager] More...
 
MoFEMErrorCode setupElementInstances ()
 [Setup algebraic structures] More...
 
MoFEMErrorCode setupSNES ()
 [Setup element instances] More...
 
MoFEMErrorCode solveProblem ()
 [Setup SNES] More...
 
MoFEMErrorCode postProcess ()
 [Solve problem] More...
 

Private Attributes

MoFEM::InterfacemField
 
PetscBool isPartitioned = PETSC_FALSE
 
int orderVelocity = 2
 
int orderPressure = 1
 
int numHoLevels = 0
 
PetscBool isDiscontPressure = PETSC_FALSE
 
PetscBool isHoGeometry = PETSC_TRUE
 
double pressureScale = 1.0
 
double lengthScale = 1.0
 
double velocityScale = 1.0
 
double reNumber = 1.0
 
double density
 
double viscosity
 
PetscBool isStokesFlow = PETSC_FALSE
 
int numSteps = 1
 
double lambdaStep = 1.0
 
double lambda = 0.0
 
int step
 
Range solidFaces
 
BitRefLevel bitLevel
 
DM dM
 
SNES snes
 
boost::shared_ptr< NavierStokesElement::CommonDatacommonData
 
SmartPetscObj< Vec > D
 
SmartPetscObj< Vec > F
 
SmartPetscObj< Mat > Aij
 
boost::shared_ptr< VolumeElementForcesAndSourcesCorefeLhsPtr
 
boost::shared_ptr< VolumeElementForcesAndSourcesCorefeRhsPtr
 
boost::shared_ptr< FaceElementForcesAndSourcesCorefeDragPtr
 
boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSidefeDragSidePtr
 
boost::shared_ptr< DirichletDisplacementBcdirichletBcPtr
 
boost::ptr_map< std::string, NeumannForcesSurfaceneumannForces
 
boost::shared_ptr< PostProcVolfePostProcPtr
 
boost::shared_ptr< PostProcFacefePostProcDragPtr
 

Detailed Description

[Example Navier Stokes]

Examples
navier_stokes.cpp.

Definition at line 22 of file navier_stokes.cpp.

Constructor & Destructor Documentation

◆ NavierStokesExample()

NavierStokesExample::NavierStokesExample ( MoFEM::Interface m_field)
inline

Definition at line 24 of file navier_stokes.cpp.

24 : mField(m_field) {
25 commonData = boost::make_shared<NavierStokesElement::CommonData>(m_field);
26 }
boost::shared_ptr< NavierStokesElement::CommonData > commonData
MoFEM::Interface & mField

◆ ~NavierStokesExample()

NavierStokesExample::~NavierStokesExample ( )
inline
Examples
navier_stokes.cpp.

Definition at line 27 of file navier_stokes.cpp.

27 {
28 CHKERR SNESDestroy(&snes);
29 CHKERR DMDestroy(&dM);
30 }
#define CHKERR
Inline error check.
Definition: definitions.h:535

Member Function Documentation

◆ defineFiniteElements()

MoFEMErrorCode NavierStokesExample::defineFiniteElements ( )
private

[Setup fields]

[Define finite elements]

Examples
navier_stokes.cpp.

Definition at line 334 of file navier_stokes.cpp.

334 {
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}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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.

◆ findBlocks()

MoFEMErrorCode NavierStokesExample::findBlocks ( )
private

[Read input]

[Find blocks]

Examples
navier_stokes.cpp.

Definition at line 182 of file navier_stokes.cpp.

182 {
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}
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#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
virtual moab::Interface & get_moab()=0

◆ postProcess()

MoFEMErrorCode NavierStokesExample::postProcess ( )
private

[Solve problem]

[Post process]

Examples
navier_stokes.cpp.

Definition at line 649 of file navier_stokes.cpp.

649 {
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}
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
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
FTensor::Index< 'i', SPACE_DIM > i
char out_file_name[255]
virtual int get_comm_rank() const =0
boost::shared_ptr< FaceElementForcesAndSourcesCore > feDragPtr
boost::shared_ptr< PostProcFace > fePostProcDragPtr
boost::shared_ptr< PostProcVol > fePostProcPtr

◆ readInput()

MoFEMErrorCode NavierStokesExample::readInput ( )
private

[Run problem]

[Read input]

Examples
navier_stokes.cpp.

Definition at line 112 of file navier_stokes.cpp.

112 {
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}
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
Managing BitRefLevels.
virtual MoFEMErrorCode rebuild_database(int verb=DEFAULT_VERBOSITY)=0
Clear database and initialize it once again.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
PetscBool isDiscontPressure

◆ runProblem()

MoFEMErrorCode NavierStokesExample::runProblem ( )

[Example Navier Stokes]

[Run problem]

Examples
navier_stokes.cpp.

Definition at line 96 of file navier_stokes.cpp.

96 {
108}
MoFEMErrorCode setupElementInstances()
[Setup algebraic structures]
MoFEMErrorCode setupDiscreteManager()
[Define finite elements]
MoFEMErrorCode defineFiniteElements()
[Setup fields]
MoFEMErrorCode setupFields()
[Find blocks]
MoFEMErrorCode setupSNES()
[Setup element instances]
MoFEMErrorCode setupAlgebraicStructures()
[Setup discrete manager]
MoFEMErrorCode findBlocks()
[Read input]
MoFEMErrorCode solveProblem()
[Setup SNES]
MoFEMErrorCode readInput()
[Run problem]

◆ setupAlgebraicStructures()

MoFEMErrorCode NavierStokesExample::setupAlgebraicStructures ( )
private

[Setup discrete manager]

[Setup algebraic structures]

Examples
navier_stokes.cpp.

Definition at line 381 of file navier_stokes.cpp.

381 {
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}
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
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
SmartPetscObj< Vec > F
SmartPetscObj< Vec > D
SmartPetscObj< Mat > Aij

◆ setupDiscreteManager()

MoFEMErrorCode NavierStokesExample::setupDiscreteManager ( )
private

[Define finite elements]

[Setup discrete manager]

Examples
navier_stokes.cpp.

Definition at line 357 of file navier_stokes.cpp.

357 {
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}
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 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 DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.

◆ setupElementInstances()

MoFEMErrorCode NavierStokesExample::setupElementInstances ( )
private

[Setup algebraic structures]

[Setup element instances]

Examples
navier_stokes.cpp.

Definition at line 404 of file navier_stokes.cpp.

404 {
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}
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:57
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.
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.
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< 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.
boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhsPtr
boost::shared_ptr< DirichletDisplacementBc > dirichletBcPtr
boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > feDragSidePtr
boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhsPtr
boost::ptr_map< std::string, NeumannForcesSurface > neumannForces

◆ setupFields()

MoFEMErrorCode NavierStokesExample::setupFields ( )
private

[Find blocks]

[Setup fields]

Examples
navier_stokes.cpp.

Definition at line 239 of file navier_stokes.cpp.

239 {
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}
@ 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
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.
Managing BitRefLevels.
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.
Projection of edge entities with one mid-node on hierarchical basis.

◆ setupSNES()

MoFEMErrorCode NavierStokesExample::setupSNES ( )
private

[Setup element instances]

[Setup SNES]

Examples
navier_stokes.cpp.

Definition at line 520 of file navier_stokes.cpp.

520 {
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}
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
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 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 DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMoFEM.cpp:1080
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:136
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
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:13

◆ solveProblem()

MoFEMErrorCode NavierStokesExample::solveProblem ( )
private

[Setup SNES]

[Solve problem]

Examples
navier_stokes.cpp.

Definition at line 556 of file navier_stokes.cpp.

556 {
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}
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
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEMErrorCode postProcess()
[Solve problem]

Member Data Documentation

◆ Aij

SmartPetscObj<Mat> NavierStokesExample::Aij
private
Examples
navier_stokes.cpp.

Definition at line 68 of file navier_stokes.cpp.

◆ bitLevel

BitRefLevel NavierStokesExample::bitLevel
private
Examples
navier_stokes.cpp.

Definition at line 59 of file navier_stokes.cpp.

◆ commonData

boost::shared_ptr<NavierStokesElement::CommonData> NavierStokesExample::commonData
private
Examples
navier_stokes.cpp.

Definition at line 64 of file navier_stokes.cpp.

◆ D

SmartPetscObj<Vec> NavierStokesExample::D
private
Examples
navier_stokes.cpp.

Definition at line 66 of file navier_stokes.cpp.

◆ density

double NavierStokesExample::density
private
Examples
navier_stokes.cpp.

Definition at line 49 of file navier_stokes.cpp.

◆ dirichletBcPtr

boost::shared_ptr<DirichletDisplacementBc> NavierStokesExample::dirichletBcPtr
private
Examples
navier_stokes.cpp.

Definition at line 76 of file navier_stokes.cpp.

◆ dM

DM NavierStokesExample::dM
private
Examples
navier_stokes.cpp.

Definition at line 61 of file navier_stokes.cpp.

◆ F

SmartPetscObj<Vec> NavierStokesExample::F
private
Examples
navier_stokes.cpp.

Definition at line 67 of file navier_stokes.cpp.

◆ feDragPtr

boost::shared_ptr<FaceElementForcesAndSourcesCore> NavierStokesExample::feDragPtr
private
Examples
navier_stokes.cpp.

Definition at line 73 of file navier_stokes.cpp.

◆ feDragSidePtr

boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> NavierStokesExample::feDragSidePtr
private
Examples
navier_stokes.cpp.

Definition at line 74 of file navier_stokes.cpp.

◆ feLhsPtr

boost::shared_ptr<VolumeElementForcesAndSourcesCore> NavierStokesExample::feLhsPtr
private
Examples
navier_stokes.cpp.

Definition at line 70 of file navier_stokes.cpp.

◆ fePostProcDragPtr

boost::shared_ptr<PostProcFace> NavierStokesExample::fePostProcDragPtr
private
Examples
navier_stokes.cpp.

Definition at line 80 of file navier_stokes.cpp.

◆ fePostProcPtr

boost::shared_ptr<PostProcVol> NavierStokesExample::fePostProcPtr
private
Examples
navier_stokes.cpp.

Definition at line 79 of file navier_stokes.cpp.

◆ feRhsPtr

boost::shared_ptr<VolumeElementForcesAndSourcesCore> NavierStokesExample::feRhsPtr
private
Examples
navier_stokes.cpp.

Definition at line 71 of file navier_stokes.cpp.

◆ isDiscontPressure

PetscBool NavierStokesExample::isDiscontPressure = PETSC_FALSE
private
Examples
navier_stokes.cpp.

Definition at line 42 of file navier_stokes.cpp.

◆ isHoGeometry

PetscBool NavierStokesExample::isHoGeometry = PETSC_TRUE
private
Examples
navier_stokes.cpp.

Definition at line 43 of file navier_stokes.cpp.

◆ isPartitioned

PetscBool NavierStokesExample::isPartitioned = PETSC_FALSE
private
Examples
navier_stokes.cpp.

Definition at line 37 of file navier_stokes.cpp.

◆ isStokesFlow

PetscBool NavierStokesExample::isStokesFlow = PETSC_FALSE
private
Examples
navier_stokes.cpp.

Definition at line 51 of file navier_stokes.cpp.

◆ lambda

double NavierStokesExample::lambda = 0.0
private
Examples
navier_stokes.cpp.

Definition at line 55 of file navier_stokes.cpp.

◆ lambdaStep

double NavierStokesExample::lambdaStep = 1.0
private
Examples
navier_stokes.cpp.

Definition at line 54 of file navier_stokes.cpp.

◆ lengthScale

double NavierStokesExample::lengthScale = 1.0
private
Examples
navier_stokes.cpp.

Definition at line 46 of file navier_stokes.cpp.

◆ mField

MoFEM::Interface& NavierStokesExample::mField
private
Examples
navier_stokes.cpp.

Definition at line 35 of file navier_stokes.cpp.

◆ neumannForces

boost::ptr_map<std::string, NeumannForcesSurface> NavierStokesExample::neumannForces
private
Examples
navier_stokes.cpp.

Definition at line 77 of file navier_stokes.cpp.

◆ numHoLevels

int NavierStokesExample::numHoLevels = 0
private
Examples
navier_stokes.cpp.

Definition at line 41 of file navier_stokes.cpp.

◆ numSteps

int NavierStokesExample::numSteps = 1
private
Examples
navier_stokes.cpp.

Definition at line 53 of file navier_stokes.cpp.

◆ orderPressure

int NavierStokesExample::orderPressure = 1
private
Examples
navier_stokes.cpp.

Definition at line 40 of file navier_stokes.cpp.

◆ orderVelocity

int NavierStokesExample::orderVelocity = 2
private
Examples
navier_stokes.cpp.

Definition at line 39 of file navier_stokes.cpp.

◆ pressureScale

double NavierStokesExample::pressureScale = 1.0
private
Examples
navier_stokes.cpp.

Definition at line 45 of file navier_stokes.cpp.

◆ reNumber

double NavierStokesExample::reNumber = 1.0
private
Examples
navier_stokes.cpp.

Definition at line 48 of file navier_stokes.cpp.

◆ snes

SNES NavierStokesExample::snes
private
Examples
navier_stokes.cpp.

Definition at line 62 of file navier_stokes.cpp.

◆ solidFaces

Range NavierStokesExample::solidFaces
private
Examples
navier_stokes.cpp.

Definition at line 58 of file navier_stokes.cpp.

◆ step

int NavierStokesExample::step
private
Examples
navier_stokes.cpp.

Definition at line 56 of file navier_stokes.cpp.

◆ velocityScale

double NavierStokesExample::velocityScale = 1.0
private
Examples
navier_stokes.cpp.

Definition at line 47 of file navier_stokes.cpp.

◆ viscosity

double NavierStokesExample::viscosity
private
Examples
navier_stokes.cpp.

Definition at line 50 of file navier_stokes.cpp.


The documentation for this struct was generated from the following file: