v0.13.0
navier_stokes.cpp

Example of viscous fluid flow problem

/** \file navier_stokes.cpp
* \example navier_stokes.cpp
*
* Example of viscous fluid flow problem
*
*/
/* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>.
* */
#include <BasicFiniteElements.hpp>
using namespace boost::numeric;
using namespace MoFEM;
using namespace std;
static char help[] = "Navier-Stokes Example\n";
//! [Example Navier Stokes]
NavierStokesExample(MoFEM::Interface &m_field) : mField(m_field) {
commonData = boost::make_shared<NavierStokesElement::CommonData>(m_field);
}
CHKERR SNESDestroy(&snes);
CHKERR DMDestroy(&dM);
}
MoFEMErrorCode runProblem();
private:
PetscBool isPartitioned = PETSC_FALSE;
int orderVelocity = 2; // default approximation orderVelocity
int orderPressure = 1; // default approximation orderPressure
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::CommonData> commonData;
SmartPetscObj<Vec> D;
SmartPetscObj<Vec> F;
SmartPetscObj<Mat> Aij;
boost::shared_ptr<VolumeElementForcesAndSourcesCore> feLhsPtr;
boost::shared_ptr<VolumeElementForcesAndSourcesCore> feRhsPtr;
boost::shared_ptr<FaceElementForcesAndSourcesCore> feDragPtr;
boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feDragSidePtr;
boost::shared_ptr<DirichletDisplacementBc> dirichletBcPtr;
boost::ptr_map<std::string, NeumannForcesSurface> neumannForces;
boost::shared_ptr<PostProcVolumeOnRefinedMesh> fePostProcPtr;
boost::shared_ptr<PostProcFaceOnRefinedMesh> fePostProcDragPtr;
MoFEMErrorCode readInput();
MoFEMErrorCode findBlocks();
MoFEMErrorCode setupFields();
MoFEMErrorCode defineFiniteElements();
MoFEMErrorCode setupDiscreteManager();
MoFEMErrorCode setupAlgebraicStructures();
MoFEMErrorCode setupElementInstances();
MoFEMErrorCode setupSNES();
MoFEMErrorCode solveProblem();
MoFEMErrorCode postProcess();
};
//! [Example Navier Stokes]
//! [Run problem]
CHKERR readInput();
CHKERR findBlocks();
CHKERR setupFields();
CHKERR defineFiniteElements();
CHKERR setupDiscreteManager();
CHKERR setupAlgebraicStructures();
CHKERR setupElementInstances();
CHKERR setupSNES();
CHKERR solveProblem();
}
//! [Run problem]
//! [Read input]
char mesh_file_name[255];
PetscBool flg_mesh_file;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "NAVIER_STOKES problem",
"none");
CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
mesh_file_name, 255, &flg_mesh_file);
CHKERR PetscOptionsBool("-my_is_partitioned", "is partitioned?", "",
isPartitioned, &isPartitioned, PETSC_NULL);
CHKERR PetscOptionsInt("-my_order_u", "approximation orderVelocity", "",
orderVelocity, &orderVelocity, PETSC_NULL);
CHKERR PetscOptionsInt("-my_order_p", "approximation orderPressure", "",
orderPressure, &orderPressure, PETSC_NULL);
CHKERR PetscOptionsInt("-my_num_ho_levels",
"number of higher order boundary levels", "",
numHoLevels, &numHoLevels, PETSC_NULL);
CHKERR PetscOptionsBool("-my_discont_pressure", "discontinuous pressure", "",
isDiscontPressure, &isDiscontPressure, PETSC_NULL);
CHKERR PetscOptionsBool("-my_ho_geometry", "use second order geometry", "",
isHoGeometry, &isHoGeometry, PETSC_NULL);
CHKERR PetscOptionsScalar("-my_length_scale", "length scale", "", lengthScale,
&lengthScale, PETSC_NULL);
CHKERR PetscOptionsScalar("-my_velocity_scale", "velocity scale", "",
velocityScale, &velocityScale, PETSC_NULL);
CHKERR PetscOptionsBool("-my_stokes_flow", "stokes flow", "", isStokesFlow,
&isStokesFlow, PETSC_NULL);
CHKERR PetscOptionsInt("-my_step_num", "number of steps", "", numSteps,
&numSteps, PETSC_NULL);
ierr = PetscOptionsEnd();
if (flg_mesh_file != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
}
// Read whole mesh or part of it if partitioned
if (isPartitioned == PETSC_TRUE) {
// This is a case of distributed mesh and algebra. In that case each
// processor keeps only part of the problem.
const char *option;
option = "PARALLEL=READ_PART;"
"PARALLEL_RESOLVE_SHARED_ENTS;"
"PARTITION=PARALLEL_PARTITION;";
CHKERR mField.get_moab().load_file(mesh_file_name, 0, option);
} else {
// In this case, we have distributed algebra, i.e. assembly of vectors and
// matrices is in parallel, but whole mesh is stored on all processors.
// snes and matrix scales well, however problem set-up of problem is
// not fully parallel.
const char *option;
option = "";
CHKERR mField.get_moab().load_file(mesh_file_name, 0, option);
}
CHKERR mField.rebuild_database();
bitLevel.set(0);
CHKERR mField.getInterface<BitRefManager>()->setBitRefLevelByDim(0, 3,
bitLevel);
}
//! [Read input]
//! [Find blocks]
if (bit->getName().compare(0, 9, "MAT_FLUID") == 0) {
const int id = bit->getMeshsetId();
CHKERR mField.get_moab().get_entities_by_dimension(
bit->getMeshset(), 3, commonData->setOfBlocksData[id].eNts, true);
std::vector<double> attributes;
bit->getAttributes(attributes);
if (attributes.size() < 2) {
SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"should be at least 2 attributes but is %d",
attributes.size());
}
commonData->setOfBlocksData[id].iD = id;
commonData->setOfBlocksData[id].fluidViscosity = attributes[0];
commonData->setOfBlocksData[id].fluidDensity = attributes[1];
viscosity = commonData->setOfBlocksData[id].fluidViscosity;
density = commonData->setOfBlocksData[id].fluidDensity;
}
}
if (bit->getName().compare(0, 9, "INT_SOLID") == 0) {
Range tets, tet;
const int id = bit->getMeshsetId();
CHKERR mField.get_moab().get_entities_by_type(
bit->getMeshset(), MBTRI, commonData->setOfFacesData[id].eNts, true);
solidFaces.merge(commonData->setOfFacesData[id].eNts);
CHKERR mField.get_moab().get_adjacencies(
commonData->setOfFacesData[id].eNts, 3, true, tets,
moab::Interface::UNION);
tet = Range(tets.front(), tets.front());
for (auto &bit : commonData->setOfBlocksData) {
if (bit.second.eNts.contains(tet)) {
commonData->setOfFacesData[id].fluidViscosity =
bit.second.fluidViscosity;
commonData->setOfFacesData[id].fluidDensity = bit.second.fluidDensity;
commonData->setOfFacesData[id].iD = id;
break;
}
}
if (commonData->setOfFacesData[id].fluidViscosity < 0) {
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Cannot find a fluid block adjacent to a given solid face");
}
}
}
}
//! [Find blocks]
//! [Setup fields]
CHKERR mField.add_field("VELOCITY", H1, AINSWORTH_LEGENDRE_BASE, 3);
if (isDiscontPressure) {
CHKERR mField.add_field("PRESSURE", L2, AINSWORTH_LEGENDRE_BASE, 1);
} else {
CHKERR mField.add_field("PRESSURE", H1, AINSWORTH_LEGENDRE_BASE, 1);
}
CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
3);
CHKERR mField.add_ents_to_field_by_dim(0, 3, "VELOCITY");
CHKERR mField.add_ents_to_field_by_dim(0, 3, "PRESSURE");
CHKERR mField.add_ents_to_field_by_dim(0, 3, "MESH_NODE_POSITIONS");
CHKERR mField.set_field_order(0, MBVERTEX, "VELOCITY", 1);
CHKERR mField.set_field_order(0, MBEDGE, "VELOCITY", orderVelocity);
CHKERR mField.set_field_order(0, MBTRI, "VELOCITY", orderVelocity);
CHKERR mField.set_field_order(0, MBTET, "VELOCITY", orderVelocity);
if (!isDiscontPressure) {
CHKERR mField.set_field_order(0, MBVERTEX, "PRESSURE", 1);
CHKERR mField.set_field_order(0, MBEDGE, "PRESSURE", orderPressure);
CHKERR mField.set_field_order(0, MBTRI, "PRESSURE", orderPressure);
}
CHKERR mField.set_field_order(0, MBTET, "PRESSURE", orderPressure);
if (numHoLevels > 0) {
std::vector<Range> levels(numHoLevels);
Range ents;
ents.merge(solidFaces);
for (int ll = 0; ll != numHoLevels; ll++) {
Range verts;
CHKERR mField.get_moab().get_connectivity(ents, verts, true);
for (auto d : {1, 2, 3}) {
CHKERR mField.get_moab().get_adjacencies(verts, d, false, ents,
moab::Interface::UNION);
}
levels[ll] = subtract(ents, ents.subset_by_type(MBVERTEX));
}
for (int ll = numHoLevels - 1; ll >= 1; ll--) {
levels[ll] = subtract(levels[ll], levels[ll - 1]);
}
int add_order = 1;
for (int ll = numHoLevels - 1; ll >= 0; ll--) {
if (isPartitioned)
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
levels[ll]);
CHKERR mField.set_field_order(levels[ll], "VELOCITY",
orderVelocity + add_order);
if (!isDiscontPressure)
CHKERR mField.set_field_order(levels[ll], "PRESSURE",
orderPressure + add_order);
else
CHKERR mField.set_field_order(levels[ll].subset_by_type(MBTET),
"PRESSURE", orderPressure + add_order);
++add_order;
}
}
CHKERR mField.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
// Set 2nd order of approximation of geometry
if (isHoGeometry) {
Range ents, edges;
CHKERR mField.get_moab().get_entities_by_dimension(0, 3, ents);
CHKERR mField.get_moab().get_adjacencies(ents, 1, false, edges,
moab::Interface::UNION);
if (isPartitioned)
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(edges);
CHKERR mField.set_field_order(edges, "MESH_NODE_POSITIONS", 2);
}
if (isPartitioned) {
CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
"VELOCITY");
CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
"PRESSURE");
CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
"MESH_NODE_POSITIONS");
}
CHKERR mField.build_fields();
Projection10NodeCoordsOnField ent_method_material(mField,
"MESH_NODE_POSITIONS");
CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
}
//! [Setup fields]
//! [Define finite elements]
// Add finite element (this defines element, declaration comes later)
CHKERR NavierStokesElement::addElement(mField, "NAVIER_STOKES", "VELOCITY",
"PRESSURE", "MESH_NODE_POSITIONS");
CHKERR NavierStokesElement::addElement(mField, "DRAG", "VELOCITY", "PRESSURE",
"MESH_NODE_POSITIONS", 2, &solidFaces);
// setup elements for loading
// build finite elements
CHKERR mField.build_finite_elements();
// build adjacencies between elements and degrees of freedom
CHKERR mField.build_adjacencies(bitLevel);
}
//! [Define finite elements]
//! [Setup discrete manager]
DMType dm_name = "DM_NAVIER_STOKES";
// Register DM problem
CHKERR DMCreate(PETSC_COMM_WORLD, &dM);
CHKERR DMSetType(dM, dm_name);
// Create DM instance
CHKERR DMMoFEMCreateMoFEM(dM, &mField, dm_name, bitLevel);
// Configure DM form line command options s
CHKERR DMSetFromOptions(dM);
// Add elements to dM (only one here)
CHKERR DMMoFEMAddElement(dM, "NAVIER_STOKES");
if (mField.check_finite_element("PRESSURE_FE"))
CHKERR DMMoFEMAddElement(dM, "PRESSURE_FE");
CHKERR DMMoFEMSetIsPartitioned(dM, isPartitioned);
// setup the DM
CHKERR DMSetUp(dM);
}
//! [Setup discrete manager]
//! [Setup algebraic structures]
CHKERR VecZeroEntries(F);
CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecZeroEntries(D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR MatSetOption(Aij, MAT_SPD, PETSC_TRUE);
CHKERR MatZeroEntries(Aij);
}
//! [Setup algebraic structures]
//! [Setup element instances]
feLhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
feRhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
feLhsPtr->getRuleHook = NavierStokesElement::VolRule();
feRhsPtr->getRuleHook = NavierStokesElement::VolRule();
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *feLhsPtr, true, false, false, true);
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *feRhsPtr, true, false, false, true);
feDragPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
feDragSidePtr = boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide>(
feDragPtr->getRuleHook = NavierStokesElement::FaceRule();
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *feDragSidePtr, true, false, false,
true);
if (isStokesFlow) {
feRhsPtr, feLhsPtr, "VELOCITY", "PRESSURE", commonData);
} else {
feRhsPtr, feLhsPtr, "VELOCITY", "PRESSURE", commonData);
}
"NAVIER_STOKES", "VELOCITY",
"PRESSURE", commonData);
dirichletBcPtr = boost::shared_ptr<DirichletDisplacementBc>(
new DirichletDisplacementBc(mField, "VELOCITY", Aij, D, F));
dirichletBcPtr->methodsOp.push_back(new NavierStokesElement::LoadScale());
// dirichletBcPtr->snes_ctx = FEMethod::CTX_SNESNONE;
dirichletBcPtr->snes_x = D;
// Assemble pressure and traction forces
NULL, "VELOCITY");
// for postprocessing:
fePostProcPtr = boost::make_shared<PostProcVolumeOnRefinedMesh>(mField);
CHKERR fePostProcPtr->generateReferenceElementMesh();
CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *fePostProcPtr, true, false, false,
true);
CHKERR fePostProcPtr->addFieldValuesPostProc("VELOCITY");
CHKERR fePostProcPtr->addFieldValuesPostProc("PRESSURE");
CHKERR fePostProcPtr->addFieldValuesPostProc("MESH_NODE_POSITIONS");
CHKERR fePostProcPtr->addFieldValuesGradientPostProc("VELOCITY");
fePostProcDragPtr = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
CHKERR fePostProcDragPtr->generateReferenceElementMesh();
CHKERR fePostProcDragPtr->addFieldValuesPostProc("MESH_NODE_POSITIONS");
fePostProcDragPtr, feDragSidePtr, "NAVIER_STOKES", "VELOCITY", "PRESSURE",
commonData);
}
//! [Setup element instances]
//! [Setup SNES]
boost::ptr_map<std::string, NeumannForcesSurface>::iterator mit =
neumannForces.begin();
for (; mit != neumannForces.end(); mit++) {
CHKERR DMMoFEMSNESSetFunction(dM, mit->first.c_str(),
&mit->second->getLoopFe(), NULL, NULL);
}
boost::shared_ptr<FEMethod> null_fe;
CHKERR DMMoFEMSNESSetFunction(dM, "NAVIER_STOKES", feRhsPtr, null_fe,
null_fe);
dirichletBcPtr);
CHKERR DMMoFEMSNESSetJacobian(dM, DM_NO_ELEMENT, null_fe, dirichletBcPtr,
null_fe);
CHKERR DMMoFEMSNESSetJacobian(dM, "NAVIER_STOKES", feLhsPtr, null_fe,
null_fe);
dirichletBcPtr);
SnesCtx *snes_ctx;
// create snes nonlinear solver
CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
CHKERR DMMoFEMGetSnesCtx(dM, &snes_ctx);
CHKERR SNESSetFunction(snes, F, SnesRhs, snes_ctx);
CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, snes_ctx);
CHKERR SNESSetFromOptions(snes);
}
//! [Setup SNES]
//! [Solve problem]
auto scale_problem = [&](double U, double L, double P) {
CHKERR mField.getInterface<FieldBlas>()->fieldScale(L,
"MESH_NODE_POSITIONS");
ProjectionFieldOn10NodeTet ent_method_on_10nodeTet(
mField, "MESH_NODE_POSITIONS", true, true);
CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_on_10nodeTet);
ent_method_on_10nodeTet.setNodes = false;
CHKERR mField.loop_dofs("MESH_NODE_POSITIONS", ent_method_on_10nodeTet);
CHKERR mField.getInterface<FieldBlas>()->fieldScale(U, "VELOCITY");
CHKERR mField.getInterface<FieldBlas>()->fieldScale(P, "PRESSURE");
};
pressureScale = viscosity * velocityScale / lengthScale;
CHKERR scale_problem(1.0 / velocityScale, 1.0 / lengthScale,
1.0 / pressureScale);
step = 0;
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Viscosity: %6.4e\n", viscosity);
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Density: %6.4e\n", density);
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Length scale: %6.4e\n", lengthScale);
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Velocity scale: %6.4e\n",
velocityScale);
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Pressure scale: %6.4e\n",
pressureScale);
if (isStokesFlow) {
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Re number: 0 (Stokes flow)\n");
} else {
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Re number: %6.4e\n",
density / viscosity * velocityScale * lengthScale);
}
lambdaStep = 1.0 / numSteps;
while (lambda < 1.0 - 1e-12) {
lambda += lambdaStep;
if (isStokesFlow) {
reNumber = 0.0;
for (auto &bit : commonData->setOfBlocksData) {
bit.second.inertiaCoef = 0.0;
bit.second.viscousCoef = 1.0;
}
} else {
reNumber = density / viscosity * velocityScale * lengthScale * lambda;
for (auto &bit : commonData->setOfBlocksData) {
bit.second.inertiaCoef = reNumber;
bit.second.viscousCoef = 1.0;
}
}
CHKERR PetscPrintf(
PETSC_COMM_WORLD,
"Step: %d | Lambda: %6.4e | Inc: %6.4e | Re number: %6.4e \n", step,
lambda, lambdaStep, reNumber);
CHKERR DMoFEMPreProcessFiniteElements(dM, dirichletBcPtr.get());
CHKERR VecAssemblyBegin(D);
CHKERR VecAssemblyEnd(D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR SNESSolve(snes, PETSC_NULL, D);
CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToGlobalVector(dM, D, INSERT_VALUES, SCATTER_REVERSE);
CHKERR scale_problem(velocityScale, lengthScale, pressureScale);
CHKERR postProcess();
CHKERR scale_problem(1.0 / velocityScale, 1.0 / lengthScale,
1.0 / pressureScale);
step++;
}
}
//! [Solve problem]
//! [Post process]
string out_file_name;
CHKERR DMoFEMLoopFiniteElements(dM, "NAVIER_STOKES", fePostProcPtr);
{
std::ostringstream stm;
stm << "out_" << step << ".h5m";
out_file_name = stm.str();
CHKERR PetscPrintf(PETSC_COMM_WORLD, "Write file %s\n",
out_file_name.c_str());
CHKERR fePostProcPtr->postProcMesh.write_file(out_file_name.c_str(), "MOAB",
"PARALLEL=WRITE_PART");
}
CHKERR VecZeroEntries(commonData->pressureDragForceVec);
CHKERR VecZeroEntries(commonData->shearDragForceVec);
CHKERR VecZeroEntries(commonData->totalDragForceVec);
CHKERR DMoFEMLoopFiniteElements(dM, "DRAG", feDragPtr);
auto get_vec_data = [&](auto vec, std::array<double, 3> &data) {
CHKERR VecAssemblyBegin(vec);
CHKERR VecAssemblyEnd(vec);
const double *array;
CHKERR VecGetArrayRead(vec, &array);
if (mField.get_comm_rank() == 0) {
for (int i : {0, 1, 2})
data[i] = array[i];
}
CHKERR VecRestoreArrayRead(vec, &array);
};
std::array<double, 3> pressure_drag;
std::array<double, 3> shear_drag;
std::array<double, 3> total_drag;
CHKERR get_vec_data(commonData->pressureDragForceVec, pressure_drag);
CHKERR get_vec_data(commonData->shearDragForceVec, shear_drag);
CHKERR get_vec_data(commonData->totalDragForceVec, total_drag);
if (mField.get_comm_rank() == 0) {
MOFEM_LOG_C("SELF", Sev::inform,
"Pressure drag force: [ %6.4e, %6.4e, %6.4e ]",
pressure_drag[0], pressure_drag[1], pressure_drag[2]);
MOFEM_LOG_C("SELF", Sev::inform,
"Shear drag force: [ %6.4e, %6.4e, %6.4e ]", shear_drag[0],
shear_drag[1], shear_drag[2]);
MOFEM_LOG_C("SELF", Sev::inform,
"Total drag force: [ %6.4e, %6.4e, %6.4e ]", total_drag[0],
total_drag[1], total_drag[2]);
}
CHKERR DMoFEMLoopFiniteElements(dM, "DRAG", fePostProcDragPtr);
{
std::ostringstream stm;
stm << "out_drag_" << step << ".h5m";
out_file_name = stm.str();
CHKERR PetscPrintf(PETSC_COMM_WORLD, "out file %s\n",
out_file_name.c_str());
CHKERR fePostProcDragPtr->postProcMesh.write_file(
out_file_name.c_str(), "MOAB", "PARALLEL=WRITE_PART");
}
}
//! [Post process]
//! [Main function]
int main(int argc, char *argv[]) {
const char param_file[] = "param_file.petsc";
// Initialise MoFEM
try {
// Create mesh database
moab::Core mb_instance; // create database
moab::Interface &moab = mb_instance; // create interface to database
// Create MoFEM database and link it to MoAB
MoFEM::Core core(moab);
MoFEM::Interface &m_field = core;
NavierStokesExample ex(m_field);
CHKERR ex.runProblem();
}
// finish work cleaning memory, getting statistics, etc
return 0;
}
//! [Main function]
static Index< 'L', 3 > L
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:22
std::string param_file
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:314
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ 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
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ BLOCKSET
Definition: definitions.h:161
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
constexpr double lambda
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1061
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:676
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:461
auto smartCreateDMMatrix
Get smart matrix from DM.
Definition: DMMoFEM.hpp:944
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:717
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
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:544
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMMoFEM.cpp:1032
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:493
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:504
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreOnSideSwitch< 0 > VolumeElementForcesAndSourcesCoreOnSide
Volume element used to integrate on skeleton.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
char mesh_file_name[255]
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
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
const double D
diffusivity
int main(int argc, char *argv[])
[Post process]
static char help[]
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:69
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeumannForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
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.
[Example Navier Stokes]
MoFEMErrorCode setupElementInstances()
[Setup algebraic structures]
MoFEMErrorCode setupDiscreteManager()
[Define finite elements]
MoFEMErrorCode defineFiniteElements()
[Setup fields]
MoFEMErrorCode setupFields()
[Find blocks]
MoFEMErrorCode setupSNES()
[Setup element instances]
MoFEMErrorCode runProblem()
[Example Navier Stokes]
MoFEMErrorCode setupAlgebraicStructures()
[Setup discrete manager]
MoFEMErrorCode findBlocks()
[Read input]
MoFEMErrorCode solveProblem()
[Setup SNES]
MoFEMErrorCode postProcess()
[Solve problem]
MoFEMErrorCode readInput()
[Run problem]