22#include <BasicFiniteElements.hpp>
24using namespace boost::numeric;
28static char help[] =
"Navier-Stokes Example\n";
36 commonData = boost::make_shared<NavierStokesElement::CommonData>(m_field);
75 boost::shared_ptr<NavierStokesElement::CommonData>
commonData;
81 boost::shared_ptr<VolumeElementForcesAndSourcesCore>
feLhsPtr;
82 boost::shared_ptr<VolumeElementForcesAndSourcesCore>
feRhsPtr;
84 boost::shared_ptr<FaceElementForcesAndSourcesCore>
feDragPtr;
85 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide>
feDragSidePtr;
126 PetscBool flg_mesh_file;
128 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"NAVIER_STOKES problem",
131 CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"mesh.h5m",
133 CHKERR PetscOptionsBool(
"-my_is_partitioned",
"is partitioned?",
"",
136 CHKERR PetscOptionsInt(
"-my_order_u",
"approximation orderVelocity",
"",
138 CHKERR PetscOptionsInt(
"-my_order_p",
"approximation orderPressure",
"",
140 CHKERR PetscOptionsInt(
"-my_num_ho_levels",
141 "number of higher order boundary levels",
"",
143 CHKERR PetscOptionsBool(
"-my_discont_pressure",
"discontinuous pressure",
"",
145 CHKERR PetscOptionsBool(
"-my_ho_geometry",
"use second order geometry",
"",
150 CHKERR PetscOptionsScalar(
"-my_velocity_scale",
"velocity scale",
"",
155 CHKERR PetscOptionsInt(
"-my_step_num",
"number of steps",
"",
numSteps,
158 ierr = PetscOptionsEnd();
160 if (flg_mesh_file != PETSC_TRUE) {
161 SETERRQ(PETSC_COMM_SELF, 1,
"*** ERROR -my_file (MESH FILE NEEDED)");
169 option =
"PARALLEL=READ_PART;"
170 "PARALLEL_RESOLVE_SHARED_ENTS;"
171 "PARTITION=PARALLEL_PARTITION;";
197 if (
bit->getName().compare(0, 9,
"MAT_FLUID") == 0) {
198 const int id =
bit->getMeshsetId();
200 bit->getMeshset(), 3,
commonData->setOfBlocksData[
id].eNts,
true);
202 std::vector<double> attributes;
203 bit->getAttributes(attributes);
204 if (attributes.size() < 2) {
206 "should be at least 2 attributes but is %d",
210 commonData->setOfBlocksData[id].fluidViscosity = attributes[0];
211 commonData->setOfBlocksData[id].fluidDensity = attributes[1];
219 if (
bit->getName().compare(0, 9,
"INT_SOLID") == 0) {
221 const int id =
bit->getMeshsetId();
223 bit->getMeshset(), MBTRI,
commonData->setOfFacesData[
id].eNts,
true);
226 commonData->setOfFacesData[
id].eNts, 3,
true, tets,
227 moab::Interface::UNION);
228 tet = Range(tets.front(), tets.front());
230 if (
bit.second.eNts.contains(tet)) {
231 commonData->setOfFacesData[id].fluidViscosity =
232 bit.second.fluidViscosity;
233 commonData->setOfFacesData[id].fluidDensity =
bit.second.fluidDensity;
238 if (
commonData->setOfFacesData[
id].fluidViscosity < 0) {
240 "Cannot find a fluid block adjacent to a given solid face");
285 for (
auto d : {1, 2, 3}) {
287 moab::Interface::UNION);
289 levels[ll] = subtract(ents, ents.subset_by_type(MBVERTEX));
292 levels[ll] = subtract(levels[ll], levels[ll - 1]);
319 moab::Interface::UNION);
331 "MESH_NODE_POSITIONS");
337 "MESH_NODE_POSITIONS");
350 "PRESSURE",
"MESH_NODE_POSITIONS");
370 DMType dm_name =
"DM_NAVIER_STOKES";
373 CHKERR DMCreate(PETSC_COMM_WORLD, &
dM);
400 CHKERR VecGhostUpdateBegin(
F, INSERT_VALUES, SCATTER_FORWARD);
401 CHKERR VecGhostUpdateEnd(
F, INSERT_VALUES, SCATTER_FORWARD);
404 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
405 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
407 CHKERR MatSetOption(
Aij, MAT_SPD, PETSC_TRUE);
418 feLhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
420 feRhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
428 feDragPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
430 feDragSidePtr = boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide>(
446 "NAVIER_STOKES",
"VELOCITY",
484 boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
488 &mit->second->getLoopFe(), NULL, NULL);
491 boost::shared_ptr<FEMethod> null_fe;
520 auto scale_problem = [&](
double U,
double L,
double P) {
523 "MESH_NODE_POSITIONS");
526 mField,
"MESH_NODE_POSITIONS",
true,
true);
528 ent_method_on_10nodeTet.
setNodes =
false;
544 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Density: %6.4e\n",
density);
546 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Velocity scale: %6.4e\n",
548 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Pressure scale: %6.4e\n",
551 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Re number: 0 (Stokes flow)\n");
553 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"Re number: %6.4e\n",
559 while (
lambda < 1.0 - 1e-12) {
566 bit.second.inertiaCoef = 0.0;
567 bit.second.viscousCoef = 1.0;
573 bit.second.viscousCoef = 1.0;
579 "Step: %d | Lambda: %6.4e | Inc: %6.4e | Re number: %6.4e \n",
step,
586 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
587 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
591 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
592 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
613 string out_file_name;
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());
623 "PARALLEL=WRITE_PART");
631 auto get_vec_data = [&](
auto vec, std::array<double, 3> &data) {
633 CHKERR VecAssemblyBegin(vec);
634 CHKERR VecAssemblyEnd(vec);
636 CHKERR VecGetArrayRead(vec, &array);
638 for (
int i : {0, 1, 2})
641 CHKERR VecRestoreArrayRead(vec, &array);
645 std::array<double, 3> pressure_drag;
646 std::array<double, 3> shear_drag;
647 std::array<double, 3> total_drag;
655 "Pressure drag force: [ %6.4e, %6.4e, %6.4e ]",
656 pressure_drag[0], pressure_drag[1], pressure_drag[2]);
658 "Shear drag force: [ %6.4e, %6.4e, %6.4e ]", shear_drag[0],
659 shear_drag[1], shear_drag[2]);
661 "Total drag force: [ %6.4e, %6.4e, %6.4e ]", total_drag[0],
662 total_drag[1], total_drag[2]);
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());
673 out_file_name.c_str(),
"MOAB",
"PARALLEL=WRITE_PART");
680int main(
int argc,
char *argv[]) {
#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.
auto smartCreateDMMatrix(DM dm)
Get smart matrix from DM.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
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 DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
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.
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
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
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
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)
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 > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
DeprecatedCoreInterface Interface
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
int main(int argc, char *argv[])
[Post process]
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.
Projection of edge entities with one mid-node on hierarchical basis.
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 setPostProcDragOperators(boost::shared_ptr< PostProcFaceOnRefinedMesh > postProcDragPtr, boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideDragFe, std::string side_fe_name, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data)
Setting up operators for post processing output of drag traction.
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const string element_name, const string velocity_field_name, const string pressure_field_name, const string mesh_field_name, const int dim=3, Range *ents=nullptr)
Setting up elements.
static MoFEMErrorCode setNavierStokesOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhs, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for solving Navier-Stokes equations.
static MoFEMErrorCode setStokesOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhs, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for solving Stokes equations.
MoFEMErrorCode setupElementInstances()
[Setup algebraic structures]
MoFEMErrorCode setupDiscreteManager()
[Define finite elements]
boost::shared_ptr< FaceElementForcesAndSourcesCore > feDragPtr
MoFEMErrorCode defineFiniteElements()
[Setup fields]
boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhsPtr
MoFEMErrorCode setupFields()
[Find blocks]
MoFEMErrorCode setupSNES()
[Setup element instances]
MoFEMErrorCode runProblem()
[Example Navier Stokes]
boost::shared_ptr< NavierStokesElement::CommonData > commonData
boost::shared_ptr< DirichletDisplacementBc > dirichletBcPtr
PetscBool isDiscontPressure
boost::shared_ptr< PostProcVolumeOnRefinedMesh > fePostProcPtr
MoFEMErrorCode setupAlgebraicStructures()
[Setup discrete manager]
NavierStokesExample(MoFEM::Interface &m_field)
boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > feDragSidePtr
boost::shared_ptr< PostProcFaceOnRefinedMesh > fePostProcDragPtr
MoFEMErrorCode findBlocks()
[Read input]
boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhsPtr
boost::ptr_map< std::string, NeumannForcesSurface > neumannForces
MoFEM::Interface & mField
MoFEMErrorCode solveProblem()
[Setup SNES]
MoFEMErrorCode postProcess()
[Solve problem]
MoFEMErrorCode readInput()
[Run problem]