v0.13.1
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< PostProcVolumeOnRefinedMeshfePostProcPtr
 
boost::shared_ptr< PostProcFaceOnRefinedMeshfePostProcDragPtr
 

Detailed Description

[Example Navier Stokes]

Examples
navier_stokes.cpp.

Definition at line 33 of file navier_stokes.cpp.

Constructor & Destructor Documentation

◆ NavierStokesExample()

NavierStokesExample::NavierStokesExample ( MoFEM::Interface m_field)
Examples
navier_stokes.cpp.

Definition at line 35 of file navier_stokes.cpp.

35 : mField(m_field) {
36 commonData = boost::make_shared<NavierStokesElement::CommonData>(m_field);
37 }
boost::shared_ptr< NavierStokesElement::CommonData > commonData
MoFEM::Interface & mField

◆ ~NavierStokesExample()

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

Definition at line 38 of file navier_stokes.cpp.

38 {
39 CHKERR SNESDestroy(&snes);
40 CHKERR DMDestroy(&dM);
41 }
#define CHKERR
Inline error check.
Definition: definitions.h:548

Member Function Documentation

◆ defineFiniteElements()

MoFEMErrorCode NavierStokesExample::defineFiniteElements ( )
private

[Setup fields]

[Define finite elements]

Examples
navier_stokes.cpp.

Definition at line 345 of file navier_stokes.cpp.

345 {
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}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
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 193 of file navier_stokes.cpp.

193 {
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}
@ BLOCKSET
Definition: definitions.h:161
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#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 610 of file navier_stokes.cpp.

610 {
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}
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:314
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
FTensor::Index< 'i', SPACE_DIM > i
virtual int get_comm_rank() const =0
boost::shared_ptr< FaceElementForcesAndSourcesCore > feDragPtr
boost::shared_ptr< PostProcVolumeOnRefinedMesh > fePostProcPtr
boost::shared_ptr< PostProcFaceOnRefinedMesh > fePostProcDragPtr

◆ readInput()

MoFEMErrorCode NavierStokesExample::readInput ( )
private

[Run problem]

[Read input]

Examples
navier_stokes.cpp.

Definition at line 123 of file navier_stokes.cpp.

123 {
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}
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
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 107 of file navier_stokes.cpp.

107 {
119}
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 392 of file navier_stokes.cpp.

392 {
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}
auto smartCreateDMMatrix(DM dm)
Get smart matrix from DM.
Definition: DMMoFEM.hpp:965
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:977
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< 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 368 of file navier_stokes.cpp.

368 {
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}
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1082
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 DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
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 415 of file navier_stokes.cpp.

415 {
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}
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:69
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.
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 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 250 of file navier_stokes.cpp.

250 {
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}
@ 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
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.
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
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 481 of file navier_stokes.cpp.

481 {
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}
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:22
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 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 DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMMoFEM.cpp:1053
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
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:27

◆ solveProblem()

MoFEMErrorCode NavierStokesExample::solveProblem ( )
private

[Setup SNES]

[Solve problem]

Examples
navier_stokes.cpp.

Definition at line 517 of file navier_stokes.cpp.

517 {
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}
static Index< 'L', 3 > L
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:494
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:505
Basic algebra on fields.
Definition: FieldBlas.hpp:31
MoFEMErrorCode postProcess()
[Solve problem]

Member Data Documentation

◆ Aij

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

Definition at line 79 of file navier_stokes.cpp.

◆ bitLevel

BitRefLevel NavierStokesExample::bitLevel
private
Examples
navier_stokes.cpp.

Definition at line 70 of file navier_stokes.cpp.

◆ commonData

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

Definition at line 75 of file navier_stokes.cpp.

◆ D

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

Definition at line 77 of file navier_stokes.cpp.

◆ density

double NavierStokesExample::density
private
Examples
navier_stokes.cpp.

Definition at line 60 of file navier_stokes.cpp.

◆ dirichletBcPtr

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

Definition at line 87 of file navier_stokes.cpp.

◆ dM

DM NavierStokesExample::dM
private
Examples
navier_stokes.cpp.

Definition at line 72 of file navier_stokes.cpp.

◆ F

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

Definition at line 78 of file navier_stokes.cpp.

◆ feDragPtr

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

Definition at line 84 of file navier_stokes.cpp.

◆ feDragSidePtr

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

Definition at line 85 of file navier_stokes.cpp.

◆ feLhsPtr

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

Definition at line 81 of file navier_stokes.cpp.

◆ fePostProcDragPtr

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

Definition at line 91 of file navier_stokes.cpp.

◆ fePostProcPtr

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

Definition at line 90 of file navier_stokes.cpp.

◆ feRhsPtr

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

Definition at line 82 of file navier_stokes.cpp.

◆ isDiscontPressure

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

Definition at line 53 of file navier_stokes.cpp.

◆ isHoGeometry

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

Definition at line 54 of file navier_stokes.cpp.

◆ isPartitioned

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

Definition at line 48 of file navier_stokes.cpp.

◆ isStokesFlow

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

Definition at line 62 of file navier_stokes.cpp.

◆ lambda

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

Definition at line 66 of file navier_stokes.cpp.

◆ lambdaStep

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

Definition at line 65 of file navier_stokes.cpp.

◆ lengthScale

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

Definition at line 57 of file navier_stokes.cpp.

◆ mField

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

Definition at line 46 of file navier_stokes.cpp.

◆ neumannForces

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

Definition at line 88 of file navier_stokes.cpp.

◆ numHoLevels

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

Definition at line 52 of file navier_stokes.cpp.

◆ numSteps

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

Definition at line 64 of file navier_stokes.cpp.

◆ orderPressure

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

Definition at line 51 of file navier_stokes.cpp.

◆ orderVelocity

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

Definition at line 50 of file navier_stokes.cpp.

◆ pressureScale

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

Definition at line 56 of file navier_stokes.cpp.

◆ reNumber

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

Definition at line 59 of file navier_stokes.cpp.

◆ snes

SNES NavierStokesExample::snes
private
Examples
navier_stokes.cpp.

Definition at line 73 of file navier_stokes.cpp.

◆ solidFaces

Range NavierStokesExample::solidFaces
private
Examples
navier_stokes.cpp.

Definition at line 69 of file navier_stokes.cpp.

◆ step

int NavierStokesExample::step
private
Examples
navier_stokes.cpp.

Definition at line 67 of file navier_stokes.cpp.

◆ velocityScale

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

Definition at line 58 of file navier_stokes.cpp.

◆ viscosity

double NavierStokesExample::viscosity
private
Examples
navier_stokes.cpp.

Definition at line 61 of file navier_stokes.cpp.


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