13using namespace boost::numeric;
17static char help[] =
"Navier-Stokes Example\n";
25 commonData = boost::make_shared<NavierStokesElement::CommonData>(m_field);
64 boost::shared_ptr<NavierStokesElement::CommonData>
commonData;
70 boost::shared_ptr<VolumeElementForcesAndSourcesCore>
feLhsPtr;
71 boost::shared_ptr<VolumeElementForcesAndSourcesCore>
feRhsPtr;
73 boost::shared_ptr<FaceElementForcesAndSourcesCore>
feDragPtr;
74 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide>
feDragSidePtr;
115 PetscBool flg_mesh_file;
117 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"NAVIER_STOKES problem",
120 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
122 CHKERR PetscOptionsBool(
"-my_is_partitioned",
"is partitioned?",
"",
125 CHKERR PetscOptionsInt(
"-my_order_u",
"approximation orderVelocity",
"",
127 CHKERR PetscOptionsInt(
"-my_order_p",
"approximation orderPressure",
"",
129 CHKERR PetscOptionsInt(
"-my_num_ho_levels",
130 "number of higher order boundary levels",
"",
132 CHKERR PetscOptionsBool(
"-my_discont_pressure",
"discontinuous pressure",
"",
134 CHKERR PetscOptionsBool(
"-my_ho_geometry",
"use second order geometry",
"",
139 CHKERR PetscOptionsScalar(
"-my_velocity_scale",
"velocity scale",
"",
144 CHKERR PetscOptionsInt(
"-my_step_num",
"number of steps",
"",
numSteps,
147 ierr = PetscOptionsEnd();
149 if (flg_mesh_file != PETSC_TRUE) {
150 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR -my_file (MESH FILE NEEDED)");
158 option =
"PARALLEL=READ_PART;"
159 "PARALLEL_RESOLVE_SHARED_ENTS;"
160 "PARTITION=PARALLEL_PARTITION;";
186 if (
bit->getName().compare(0, 9,
"MAT_FLUID") == 0) {
187 const int id =
bit->getMeshsetId();
189 bit->getMeshset(), 3,
commonData->setOfBlocksData[
id].eNts,
true);
191 std::vector<double> attributes;
192 bit->getAttributes(attributes);
193 if (attributes.size() < 2) {
195 "should be at least 2 attributes but is %d",
199 commonData->setOfBlocksData[id].fluidViscosity = attributes[0];
200 commonData->setOfBlocksData[id].fluidDensity = attributes[1];
208 if (
bit->getName().compare(0, 9,
"INT_SOLID") == 0) {
210 const int id =
bit->getMeshsetId();
212 bit->getMeshset(), MBTRI,
commonData->setOfFacesData[
id].eNts,
true);
215 commonData->setOfFacesData[
id].eNts, 3,
true, tets,
216 moab::Interface::UNION);
217 tet =
Range(tets.front(), tets.front());
219 if (
bit.second.eNts.contains(tet)) {
220 commonData->setOfFacesData[id].fluidViscosity =
221 bit.second.fluidViscosity;
222 commonData->setOfFacesData[id].fluidDensity =
bit.second.fluidDensity;
227 if (
commonData->setOfFacesData[
id].fluidViscosity < 0) {
229 "Cannot find a fluid block adjacent to a given solid face");
274 for (
auto d : {1, 2, 3}) {
276 moab::Interface::UNION);
278 levels[ll] = subtract(ents, ents.subset_by_type(MBVERTEX));
281 levels[ll] = subtract(levels[ll], levels[ll - 1]);
308 moab::Interface::UNION);
320 "MESH_NODE_POSITIONS");
326 "MESH_NODE_POSITIONS");
339 "PRESSURE",
"MESH_NODE_POSITIONS");
359 DMType dm_name =
"DM_NAVIER_STOKES";
362 CHKERR DMCreate(PETSC_COMM_WORLD, &
dM);
389 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
390 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
393 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
394 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
396 CHKERR MatSetOption(
Aij, MAT_SPD, PETSC_TRUE);
407 feLhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
409 feRhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
419 feDragPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
421 feDragSidePtr = boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide>(
437 "NAVIER_STOKES",
"VELOCITY",
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>();
477 {{
"PRESSURE", p_ptr}},
479 {{
"VELOCITY", v_ptr}, {
"MESH_NODE_POSITIONS", pos_ptr}},
481 {{
"VELOCITY_GRAD", grad_ptr}},
489 fePostProcDragPtr = boost::make_shared<PostProcFace>(mField);
490 fePostProcDragPtr->getOpPtrVector().push_back(
492 fePostProcDragPtr->getOpPtrVector().push_back(
496 fePostProcDragPtr->getPostProcMesh(),
497 fePostProcDragPtr->getMapGaussPts(),
501 {{
"MESH_NODE_POSITIONS", pos_ptr}},
512 fePostProcDragPtr, feDragSidePtr,
"NAVIER_STOKES",
"VELOCITY",
"PRESSURE",
523 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
527 &mit->second->getLoopFe(), NULL, NULL);
530 boost::shared_ptr<FEMethod> null_fe;
559 auto scale_problem = [&](
double U,
double L,
double P) {
562 "MESH_NODE_POSITIONS");
565 mField,
"MESH_NODE_POSITIONS",
true,
true);
567 ent_method_on_10nodeTet.
setNodes =
false;
583 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Density: %6.4e\n",
density);
585 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Velocity scale: %6.4e\n",
587 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Pressure scale: %6.4e\n",
590 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Re number: 0 (Stokes flow)\n");
592 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Re number: %6.4e\n",
598 while (
lambda < 1.0 - 1e-12) {
605 bit.second.inertiaCoef = 0.0;
606 bit.second.viscousCoef = 1.0;
612 bit.second.viscousCoef = 1.0;
618 "Step: %d | Lambda: %6.4e | Inc: %6.4e | Re number: %6.4e \n",
step,
625 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
626 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
630 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
631 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
656 std::ostringstream stm;
657 stm <<
"out_" <<
step <<
".h5m";
659 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Write file %s\n",
669 auto get_vec_data = [&](
auto vec, std::array<double, 3> &data) {
671 CHKERR VecAssemblyBegin(vec);
672 CHKERR VecAssemblyEnd(vec);
674 CHKERR VecGetArrayRead(vec, &array);
676 for (
int i : {0, 1, 2})
679 CHKERR VecRestoreArrayRead(vec, &array);
683 std::array<double, 3> pressure_drag;
684 std::array<double, 3> shear_drag;
685 std::array<double, 3> total_drag;
693 "Pressure drag force: [ %6.4e, %6.4e, %6.4e ]",
694 pressure_drag[0], pressure_drag[1], pressure_drag[2]);
696 "Shear drag force: [ %6.4e, %6.4e, %6.4e ]", shear_drag[0],
697 shear_drag[1], shear_drag[2]);
699 "Total drag force: [ %6.4e, %6.4e, %6.4e ]", total_drag[0],
700 total_drag[1], total_drag[2]);
705 std::ostringstream stm;
706 stm <<
"out_drag_" <<
step <<
".h5m";
708 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"out file %s\n",
717int main(
int argc,
char *argv[]) {
726 moab::Core mb_instance;
727 moab::Interface &moab = mb_instance;
#define MOFEM_LOG_C(channel, severity, format,...)
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ L2
field with C-1 continuity
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#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)
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
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.
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm)
Get smart vector from DM.
auto createDMMatrix(DM dm)
Get smart matrix from DM.
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
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.
FTensor::Index< 'i', SPACE_DIM > i
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.
implementation of Data Operators for Forces and Sources
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.
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.
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)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Set Dirichlet boundary conditions on displacements.
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.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
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.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
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 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.
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]
boost::shared_ptr< NavierStokesElement::CommonData > commonData
boost::shared_ptr< DirichletDisplacementBc > dirichletBcPtr
PetscBool isDiscontPressure
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]