NavierStokesExample Struct Reference

[Example Navier Stokes] More...

Collaboration diagram for NavierStokesExample:

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::Interface &	mField
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::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< PostProcVol >	fePostProcPtr
boost::shared_ptr< PostProcFace >	fePostProcDragPtr

Detailed Description

[Example Navier Stokes]

Examples

navier_stokes.cpp.

Definition at line 22 of file navier_stokes.cpp.

Constructor & Destructor Documentation

NavierStokesExample()

NavierStokesExample::NavierStokesExample ( MoFEM::Interface &  m_field )

inline

Definition at line 24 of file navier_stokes.cpp.

                                               : mField(m_field) {
    commonData = boost::make_shared<NavierStokesElement::CommonData>(m_field);
  }

~NavierStokesExample()

NavierStokesExample::~NavierStokesExample ( )

inline

Definition at line 27 of file navier_stokes.cpp.

                         {
    CHKERR SNESDestroy(&snes);
    CHKERR DMDestroy(&dM);
  }

Member Function Documentation

defineFiniteElements()

MoFEMErrorCode NavierStokesExample::defineFiniteElements ( )

private

[Setup fields]

[Define finite elements]

Examples

navier_stokes.cpp.

Definition at line 334 of file navier_stokes.cpp.

                                                         {
  MoFEMFunctionBegin;
 
  // 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
  CHKERR MetaNeumannForces::addNeumannBCElements(mField, "VELOCITY");
 
  // build finite elements
  CHKERR mField.build_finite_elements();
  // build adjacencies between elements and degrees of freedom
  CHKERR mField.build_adjacencies(bitLevel);
 
  MoFEMFunctionReturn(0);
}

findBlocks()

MoFEMErrorCode NavierStokesExample::findBlocks ( )

private

[Read input]

[Find blocks]

Examples

navier_stokes.cpp.

Definition at line 182 of file navier_stokes.cpp.

                                               {
  MoFEMFunctionBegin;
 
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, bit)) {
    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;
    }
  }
 
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, bit)) {
    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");
      }
    }
  }
 
  MoFEMFunctionReturn(0);
}

postProcess()

MoFEMErrorCode NavierStokesExample::postProcess ( )

private

[Solve problem]

[Post process]

Examples

navier_stokes.cpp.

Definition at line 649 of file navier_stokes.cpp.

                                                {
  MoFEMFunctionBegin;
 
  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->writeFile(out_file_name.c_str());
  }
 
  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) {
    MoFEMFunctionBegin;
    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);
    MoFEMFunctionReturn(0);
  };
 
  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->writeFile(out_file_name);
  }
  MoFEMFunctionReturn(0);
}

readInput()

MoFEMErrorCode NavierStokesExample::readInput ( )

private

[Run problem]

[Read input]

Examples

navier_stokes.cpp.

Definition at line 112 of file navier_stokes.cpp.

                                              {
  MoFEMFunctionBegin;
  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);
 
  MoFEMFunctionReturn(0);
}

runProblem()

MoFEMErrorCode NavierStokesExample::runProblem

(

)

[Example Navier Stokes]

[Run problem]

Examples

navier_stokes.cpp.

Definition at line 96 of file navier_stokes.cpp.

                                               {
  MoFEMFunctionBegin;
  CHKERR readInput();
  CHKERR findBlocks();
  CHKERR setupFields();
  CHKERR defineFiniteElements();
  CHKERR setupDiscreteManager();
  CHKERR setupAlgebraicStructures();
  CHKERR setupElementInstances();
  CHKERR setupSNES();
  CHKERR solveProblem();
  MoFEMFunctionReturn(0);
}

setupAlgebraicStructures()

MoFEMErrorCode NavierStokesExample::setupAlgebraicStructures ( )

private

[Setup discrete manager]

[Setup algebraic structures]

Examples

navier_stokes.cpp.

Definition at line 381 of file navier_stokes.cpp.

                                                             {
  MoFEMFunctionBegin;
 
  D = createDMVector(dM);
  F = vectorDuplicate(D);
  Aij = createDMMatrix(dM);
 
  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);
 
  MoFEMFunctionReturn(0);
}

setupDiscreteManager()

MoFEMErrorCode NavierStokesExample::setupDiscreteManager ( )

private

[Define finite elements]

[Setup discrete manager]

Examples

navier_stokes.cpp.

Definition at line 357 of file navier_stokes.cpp.

                                                         {
  MoFEMFunctionBegin;
  DMType dm_name = "DM_NAVIER_STOKES";
  // Register DM problem
  CHKERR DMRegister_MoFEM(dm_name);
  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");
  CHKERR DMMoFEMAddElement(dM, "DRAG");
  if (mField.check_finite_element("PRESSURE_FE"))
    CHKERR DMMoFEMAddElement(dM, "PRESSURE_FE");
  CHKERR DMMoFEMSetIsPartitioned(dM, isPartitioned);
  // setup the DM
  CHKERR DMSetUp(dM);
  MoFEMFunctionReturn(0);
}

setupElementInstances()

MoFEMErrorCode NavierStokesExample::setupElementInstances ( )

private

[Setup algebraic structures]

[Setup element instances]

Examples

navier_stokes.cpp.

Definition at line 404 of file navier_stokes.cpp.

                                                          {
  MoFEMFunctionBegin;
 
  feLhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
      new VolumeElementForcesAndSourcesCore(mField));
  feRhsPtr = boost::shared_ptr<VolumeElementForcesAndSourcesCore>(
      new VolumeElementForcesAndSourcesCore(mField));
 
  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>(
      new FaceElementForcesAndSourcesCore(mField));
  feDragSidePtr = boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide>(
      new VolumeElementForcesAndSourcesCoreOnSide(mField));
 
  feDragPtr->getRuleHook = NavierStokesElement::FaceRule();
  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *feDragSidePtr, true, false, false,
                     true);
 
  if (isStokesFlow) {
    CHKERR NavierStokesElement::setStokesOperators(
        feRhsPtr, feLhsPtr, "VELOCITY", "PRESSURE", commonData);
  } else {
    CHKERR NavierStokesElement::setNavierStokesOperators(
        feRhsPtr, feLhsPtr, "VELOCITY", "PRESSURE", commonData);
  }
 
  NavierStokesElement::setCalcDragOperators(feDragPtr, feDragSidePtr,
                                            "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
  CHKERR MetaNeumannForces::setMomentumFluxOperators(mField, neumannForces,
                                                     NULL, "VELOCITY");
 
  // for postprocessing:
  fePostProcPtr = boost::make_shared<PostProcVol>(mField);
  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", *fePostProcPtr, true, false, false,
                     true);
 
  auto v_ptr = boost::make_shared<MatrixDouble>();
  auto grad_ptr = boost::make_shared<MatrixDouble>();
  auto pos_ptr = boost::make_shared<MatrixDouble>();
  auto p_ptr = boost::make_shared<VectorDouble>();
 
  fePostProcPtr->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<3>("VELOCITY", v_ptr));
  fePostProcPtr->getOpPtrVector().push_back(
      new OpCalculateVectorFieldGradient<3, 3>("VELOCITY", grad_ptr));
  fePostProcPtr->getOpPtrVector().push_back(
      new OpCalculateScalarFieldValues("PRESSURE", p_ptr));
  fePostProcPtr->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS", pos_ptr));
 
  using OpPPMap = OpPostProcMapInMoab<3, 3>;
 
  fePostProcPtr->getOpPtrVector().push_back(
 
      new OpPPMap(
 
          fePostProcPtr->getPostProcMesh(), fePostProcPtr->getMapGaussPts(),
 
          {{"PRESSURE", p_ptr}},
 
          {{"VELOCITY", v_ptr}, {"MESH_NODE_POSITIONS", pos_ptr}},
 
          {{"VELOCITY_GRAD", grad_ptr}},
 
          {}
 
          )
 
  );
 
  fePostProcDragPtr = boost::make_shared<PostProcFace>(mField);
  fePostProcDragPtr->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<3>("MESH_NODE_POSITIONS", pos_ptr));
  fePostProcDragPtr->getOpPtrVector().push_back(
 
      new OpPPMap(
 
          fePostProcDragPtr->getPostProcMesh(),
          fePostProcDragPtr->getMapGaussPts(),
 
          {},
 
          {{"MESH_NODE_POSITIONS", pos_ptr}},
 
          {},
 
          {}
 
          )
 
  );
 
  CHKERR NavierStokesElement::setPostProcDragOperators(
      fePostProcDragPtr, feDragSidePtr, "NAVIER_STOKES", "VELOCITY", "PRESSURE",
      commonData);
 
  MoFEMFunctionReturn(0);
}

setupFields()

MoFEMErrorCode NavierStokesExample::setupFields ( )

private

[Find blocks]

[Setup fields]

Examples

navier_stokes.cpp.

Definition at line 239 of file navier_stokes.cpp.

                                                {
  MoFEMFunctionBegin;
 
  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);
 
  MoFEMFunctionReturn(0);
}

setupSNES()

MoFEMErrorCode NavierStokesExample::setupSNES ( )

private

[Setup element instances]

[Setup SNES]

Examples

navier_stokes.cpp.

Definition at line 520 of file navier_stokes.cpp.

                                              {
  MoFEMFunctionBegin;
 
  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);
  CHKERR DMMoFEMSNESSetFunction(dM, DM_NO_ELEMENT, 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);
  CHKERR DMMoFEMSNESSetJacobian(dM, DM_NO_ELEMENT, 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);
 
  MoFEMFunctionReturn(0);
}

solveProblem()

MoFEMErrorCode NavierStokesExample::solveProblem ( )

private

[Setup SNES]

[Solve problem]

Examples

navier_stokes.cpp.

Definition at line 556 of file navier_stokes.cpp.

                                                 {
  MoFEMFunctionBegin;
 
  auto scale_problem = [&](double U, double L, double P) {
    MoFEMFunctionBegin;
    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");
    MoFEMFunctionReturn(0);
  };
 
  pressureScale = viscosity * velocityScale / lengthScale;
  NavierStokesElement::LoadScale::lambda = 1.0 / velocityScale;
  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++;
  }
 
  MoFEMFunctionReturn(0);
}

Member Data Documentation

Aij

SmartPetscObj<Mat> NavierStokesExample::Aij

private

Definition at line 68 of file navier_stokes.cpp.

bitLevel

BitRefLevel NavierStokesExample::bitLevel

private

Definition at line 59 of file navier_stokes.cpp.

commonData

boost::shared_ptr<NavierStokesElement::CommonData> NavierStokesExample::commonData

private

Definition at line 64 of file navier_stokes.cpp.

D

SmartPetscObj<Vec> NavierStokesExample::D

private

Definition at line 66 of file navier_stokes.cpp.

density

double NavierStokesExample::density

private

Definition at line 49 of file navier_stokes.cpp.

dirichletBcPtr

boost::shared_ptr<DirichletDisplacementBc> NavierStokesExample::dirichletBcPtr

private

Definition at line 76 of file navier_stokes.cpp.

dM

DM NavierStokesExample::dM

private

Definition at line 61 of file navier_stokes.cpp.

F

SmartPetscObj<Vec> NavierStokesExample::F

private

Definition at line 67 of file navier_stokes.cpp.

feDragPtr

boost::shared_ptr<FaceElementForcesAndSourcesCore> NavierStokesExample::feDragPtr

private

Definition at line 73 of file navier_stokes.cpp.

feDragSidePtr

boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> NavierStokesExample::feDragSidePtr

private

Definition at line 74 of file navier_stokes.cpp.

feLhsPtr

boost::shared_ptr<VolumeElementForcesAndSourcesCore> NavierStokesExample::feLhsPtr

private

Definition at line 70 of file navier_stokes.cpp.

fePostProcDragPtr

boost::shared_ptr<PostProcFace> NavierStokesExample::fePostProcDragPtr

private

Definition at line 80 of file navier_stokes.cpp.

fePostProcPtr

boost::shared_ptr<PostProcVol> NavierStokesExample::fePostProcPtr

private

Definition at line 79 of file navier_stokes.cpp.

feRhsPtr

boost::shared_ptr<VolumeElementForcesAndSourcesCore> NavierStokesExample::feRhsPtr

private

Definition at line 71 of file navier_stokes.cpp.

isDiscontPressure

PetscBool NavierStokesExample::isDiscontPressure = PETSC_FALSE

private

Definition at line 42 of file navier_stokes.cpp.

isHoGeometry

PetscBool NavierStokesExample::isHoGeometry = PETSC_TRUE

private

Definition at line 43 of file navier_stokes.cpp.

isPartitioned

PetscBool NavierStokesExample::isPartitioned = PETSC_FALSE

private

Definition at line 37 of file navier_stokes.cpp.

isStokesFlow

PetscBool NavierStokesExample::isStokesFlow = PETSC_FALSE

private

Definition at line 51 of file navier_stokes.cpp.

lambda

double NavierStokesExample::lambda = 0.0

private

Definition at line 55 of file navier_stokes.cpp.

lambdaStep

double NavierStokesExample::lambdaStep = 1.0

private

Definition at line 54 of file navier_stokes.cpp.

lengthScale

double NavierStokesExample::lengthScale = 1.0

private

Definition at line 46 of file navier_stokes.cpp.

mField

MoFEM::Interface& NavierStokesExample::mField

private

Definition at line 35 of file navier_stokes.cpp.

neumannForces

boost::ptr_map<std::string, NeumannForcesSurface> NavierStokesExample::neumannForces

private

Definition at line 77 of file navier_stokes.cpp.

numHoLevels

int NavierStokesExample::numHoLevels = 0

private

Definition at line 41 of file navier_stokes.cpp.

numSteps

int NavierStokesExample::numSteps = 1

private

Definition at line 53 of file navier_stokes.cpp.

orderPressure

int NavierStokesExample::orderPressure = 1

private

Definition at line 40 of file navier_stokes.cpp.

orderVelocity

int NavierStokesExample::orderVelocity = 2

private

Definition at line 39 of file navier_stokes.cpp.

pressureScale

double NavierStokesExample::pressureScale = 1.0

private

Definition at line 45 of file navier_stokes.cpp.

reNumber

double NavierStokesExample::reNumber = 1.0

private

Definition at line 48 of file navier_stokes.cpp.

snes

SNES NavierStokesExample::snes

private

Definition at line 62 of file navier_stokes.cpp.

solidFaces

Range NavierStokesExample::solidFaces

private

Definition at line 58 of file navier_stokes.cpp.

step

int NavierStokesExample::step

private

Definition at line 56 of file navier_stokes.cpp.

velocityScale

double NavierStokesExample::velocityScale = 1.0

private

Definition at line 47 of file navier_stokes.cpp.

viscosity

double NavierStokesExample::viscosity

private

Definition at line 50 of file navier_stokes.cpp.

---

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

• users_modules/tutorials/cor-10/navier_stokes.cpp