deal.II version GIT relicensing-3356-g4636aadadd 2025-05-22 18:40:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
The 'Two phase flow interaction ' code gallery program

This program was contributed by Manuel Quezada de Luna <manuel.quezada.dl@gmail.com>.
It comes without any warranty or support by its authors or the authors of deal.II.

This program is part of the deal.II code gallery and consists of the following files (click to inspect):

Pictures from this code gallery program

Annotated version of Readme.md

Two Phase Flow

General description of the problem

We consider the problem of two-phase incompressible flow. We start with an initial state of two phases (fluids) that define density and viscosity fields. Using these fields we solve the incompressible Navier-Stokes equations to obtain a velocity field.

We use the initial state to define a representation of the interface via a Level Set function \(\phi\in[-1, 1]\). The zero level set \(\{\phi=0\}\) defines the interface of the phases. Positive values of the level set function represent water while negative values represent air.

Using the velocity field from the Navier-Stokes equations we transport the level set function. To do this we assume the velocity is divergence free and write the transport equation in conservation form.

Using the advected level set function we reconstruct density and viscosity fields. We repeat the process until the final desired time.

The Navier-Stokes equations are solved using a projection scheme based on [1]. To solve the level set we use continuous Galerkin Finite Elements with high-order stabilization based on the entropy residual of the solution [2] and artificial compression inspired by [3] and [4].


General description of the code

Driver code: MultiPhase

The driver code of the simulation is the run function within MultiPhase.cc. The general idea is to define here everything that has to do with the problem, set all the (physical and numerical) parameters and perform the time loop. The run function does the following: Set some physical parameters like final time, density and viscosity coefficients, etc. and numerical parameters like cfl, numerical constants, algorithms to be used, etc. Creates the geometry for the specified problem. Currently we have the following problems: Breaking Dam problem in 2D. Filling a tank in 2D. Small wave perturbation in 2D. Falling drop in 2D. Creates an object of the class NavierStokesSolver and an object of the class LevelSetSolver.
Set the initial condition for each of the solvers. Performs the time loop. Within the time loop we do the following: Pass the current level set function to the Navier Stokes Solver. Ask the Navier Stokes Solver to perform one time step. Get the velocity field from the Navier Stokes Solver. Pass the velocity field to the Level Set Solver. Ask the Level Set Solver to perform one time step. Get the level set function from the Level Set Solver. Repeat until the final time. Output the solution at the requested times.

Navier Stokes Solver

The NavierStokesSolver class is responsible for solving the Navier Stokes equation for just one time step. It requires density and viscosity information. This information can be passed by either a function or by passing a vector containing the DOFs of the level set function. For this reason the class contains the following two constructors: First constructor. Here we have to pass density and viscosity constants for the two phases. In addition, we have to pass a vector of DOFs defining the level set function. This constructor is meant to be used during the two-phase flow simulations. Second constructor. Here we have to pass functions to define the viscosity and density fields. This is meant to test the convergence properties of the method (and to validate the implementation).

Level Set Solver

The LevelSetSolver.cc code is responsible for solving the Level Set for just one time step. It requires information about the velocity field and provides the transported level set function. The velocity field can be interpolated (outside of this class) from a given function to test the method (and to validate the implementation). Alternatively, the velocity can be provided from the solution of the Navier-Stokes equations (for the two phase flow simulations).

Testing the Navier Stokes Solver

The TestNavierStokes.cc code is used to test the convergence (in time) of the Navier-Stokes solver. To run it uncomment the line SET(TARGET "TestNavierStokes") within CMakeLists.txt (and make sure to comment SET(TARGET "TestLevelSet") and SET(TARGET "MultiPhase"). Then cmake and compile. The convergence can be done in 2 or 3 dimensions. Different exact solutions (and force terms) are used in each case. The dimension can be set in the line TestNavierStokes<2> test_navier_stokes(degree_LS, degree_U) within the main function.

Testing the Level Set Solver

The TestLevelSet.cc code is used to test the level set solver. To run it uncomment the corresponding line within CMakeLists.txt. Then cmake and compile. There are currently just two problems implemented: diagonal advection and circular rotation. If the velocity is independent of time set the flag VARIABLE_VELOCITY to zero to avoid interpolating the velocity field at every time step.

Utility files

The files utilities.cc, utilities_test_LS.cc and utilities_test_NS.cc contain functions required in MultiPhase.cc, TestLevelSet.cc and TestNavierStokes.cc respectively. The script clean.sh ereases all files created by cmake, compile and run any example.


References

[1] J.-L. Guermond and A. Salgado. A splitting method for incompressible flows with variable density based on a pressure Poisson equation. Journal of Computational Physics, 228(8):2834–2846, 2009.

[2] J.-L. Guermond, R. Pasquetti, and B. Popov. Entropy viscosity method for nonlinear conservation laws. Journal of Computational Physics, 230(11):4248– 4267, 2011.

[3] A. Harten. The artificial compression method for computation of shocks and contact discontinuities. I. Single conservation laws. Communications on Pure and Applied Mathematics, 30(5):611–638, 1977.

[4] A. Harten. The artificial compression method for computation of shocks and contact discontinuities. III. Self-adjusting hybrid schemes. Mathematics of Computation, 32:363–389, 1978.

Annotated version of LevelSetSolver.cc

  /* -----------------------------------------------------------------------------
  *
  * SPDX-License-Identifier: LGPL-2.1-or-later
  * Copyright (C) 2016 Manuel Quezada de Luna
  *
  * This file is part of the deal.II code gallery.
  *
  * -----------------------------------------------------------------------------
  */
  #include <deal.II/base/function.h>
  #include <deal.II/lac/affine_constraints.h>
  #include <deal.II/lac/vector.h>
  #include <deal.II/lac/petsc_vector.h>
  #include <deal.II/dofs/dof_handler.h>
  #include <deal.II/dofs/dof_tools.h>
  #include <deal.II/fe/fe_values.h>
  #include <deal.II/fe/fe_q.h>
  #include <deal.II/numerics/vector_tools.h>
  #include <deal.II/numerics/data_out.h>
  #include <deal.II/numerics/error_estimator.h>
  #include <deal.II/base/utilities.h>
  #include <deal.II/base/index_set.h>
  #include <deal.II/distributed/tria.h>
  #include <deal.II/distributed/grid_refinement.h>
  #include <deal.II/lac/vector.h>
  #include <deal.II/base/timer.h>
  #include <deal.II/grid/grid_tools.h>
  #include <deal.II/fe/mapping_q.h>
  #include <mpi.h>
  using namespace dealii;

FLAGS

LOG FOR LEVEL SET FROM -1 to 1

  #define ENTROPY_GRAD(phi,phix) 2*phi*phix*((1-phi*phi>=0) ? -1 : 1)/(std::abs(1-phi*phi)+1E-14)
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)

////////////////////////////////////////////////////// ////////////////// TRANSPORT SOLVER ////////////////// ////////////////////////////////////////////////////// This is a solver for the transpor solver. We assume the velocity is divergence free and solve the equation in conservation form. /////////////////////////////// -------— NOTATION -------— /////////////////////////////// We use notation popular in the literature of conservation laws. For this reason the solution is denoted as u, unm1, unp1, etc. and the velocity is treated as vx, vy and vz.

  template <int dim>
  {
  public:

//////////////////// INITIAL CONDITIONS ////////////////////

///////////////////// BOUNDARY CONDITIONS /////////////////////

  void set_boundary_conditions(std::vector<types::global_dof_index> &boundary_values_id_u,
  std::vector<double> boundary_values_u);

////////////// SET VELOCITY //////////////

/////////////////// SET AND GET ALPHA ///////////////////

/////////////// NTH TIME STEP ///////////////

/////// SETUP ///////

  void setup();
  LevelSetSolver (const unsigned int degree_LS,
  const unsigned int degree_U,
  const double time_step,
  const double cK,
  const double cE,
  const bool verbose,
  std::string ALGORITHM,
  const unsigned int TIME_INTEGRATION,
  MPI_Comm &mpi_communicator);
  private:
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

//////////////////////////////////// ASSEMBLE MASS (and other) MATRICES ////////////////////////////////////

  void assemble_ML();
  void invert_ML();
  void assemble_MC();

////////////////////////////////// LOW ORDER METHOD (DiJ Viscosity) //////////////////////////////////

/////////////////// ENTROPY VISCOSITY ///////////////////

/////////////////////// FOR MAXIMUM PRINCIPLE ///////////////////////

/////////////////// COMPUTE SOLUTIONS ///////////////////

/////////// UTILITIES ///////////

  void get_sparsity_pattern();
  void solve(const AffineConstraints<double> &constraints,
  std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,

/////////////////// MY PETSC WRAPPERS ///////////////////

  const std::vector<types::global_dof_index> &indices,
  std::vector<PetscScalar> &values);
  const std::vector<types::global_dof_index> &indices,
  std::map<types::global_dof_index, types::global_dof_index> &map_from_Q1_to_Q2,
  std::vector<PetscScalar> &values);
  MPI_Comm mpi_communicator;

FINITE ELEMENT SPACE

OPERATORS times SOLUTION VECTOR

MASS MATRIX

  std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> MC_preconditioner;

BOUNDARIES

  std::vector<types::global_dof_index> boundary_values_id_u;
  std::vector<double> boundary_values_u;

////////// MATRICES ////////// FOR FIRST ORDER VISCOSITY

FOR ENTROPY VISCOSITY

FOR FCT (flux and limited flux)

FOR ITERATIVE FCT

GHOSTED VECTORS

NON-GHOSTED VECTORS

LUMPED MASS MATRIX

CONSTRAINTS

TIME STEPPING

  double time_step;

SOME PARAMETERS

  double cE, cK;
  double solver_tolerance;

UTILITIES

  bool verbose;
  std::string ALGORITHM;
  unsigned int TIME_INTEGRATION;
  std::map<types::global_dof_index, types::global_dof_index> map_from_Q1_to_Q2;
  std::map<types::global_dof_index, std::vector<types::global_dof_index> > sparsity_pattern;
  };
  template <int dim>
  const unsigned int degree_U,
  const double time_step,
  const double cK,
  const double cE,
  const bool verbose,
  std::string ALGORITHM,
  const unsigned int TIME_INTEGRATION,
  MPI_Comm &mpi_communicator)
  :
  mpi_communicator (mpi_communicator),
  cE(cE),
  cK(cK),
  pcout (std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator)== 0))
  {
  pcout << "********** LEVEL SET SETUP **********" << std::endl;
  setup();
  }
  template <int dim>
  {
  }
constexpr void clear()
friend class Tensor
Definition tensor.h:866
STL namespace.

/////////////////////////////////////////////////////// /////////////////// PUBLIC FUNCTIONS ////////////////// /////////////////////////////////////////////////////// //////////////////////////////////// //////// INITIAL CONDITIONS //////// ////////////////////////////////////

initialize old vectors with current solution, this just happens the first time

initialize old vectors with current solution, this just happens the first time

///////////////////////////////////// //////// BOUNDARY CONDITIONS //////// /////////////////////////////////////

  template <int dim>
  void LevelSetSolver<dim>::set_boundary_conditions(std::vector<types::global_dof_index> &boundary_values_id_u,
  std::vector<double> boundary_values_u)
  {
  this->boundary_values_id_u = boundary_values_id_u;
  this->boundary_values_u = boundary_values_u;
  }

////////////////////////////// //////// SET VELOCITY //////// //////////////////////////////

SAVE OLD SOLUTION

update velocity

SAVE OLD SOLUTION

update velocity

  this->locally_relevant_solution_vx=locally_relevant_solution_vx;
  this->locally_relevant_solution_vy=locally_relevant_solution_vy;
  this->locally_relevant_solution_vz=locally_relevant_solution_vz;
  }

/////////////////////////////// //////// SET AND GET U //////// ///////////////////////////////


---------------------------— COMPUTE SOLUTIONS ---------------------------—

COMPUTE SOLUTION

BOUNDARY CONDITIONS

CHECK MAXIMUM PRINCIPLE

pcout << "*********************************************************************... " << unp1.min() << ", " << unp1.max() << std::endl;


---------------------------— SETUP ---------------------------—

  template <int dim>
  {
  solver_tolerance=1E-6;
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)

//////////////////////// SETUP FOR DOF HANDLERS //////////////////////// setup system LS

  dof_handler_LS.distribute_dofs (fe_LS);
  locally_owned_dofs_LS = dof_handler_LS.locally_owned_dofs ();
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)

setup system U

////////////////// INIT CONSTRAINTS //////////////////

  constraints.clear ();
  constraints.reinit (locally_relevant_dofs_LS);
  constraints.close ();
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)

///////////////////// NON-GHOSTED VECTORS /////////////////////

  MPP_uL_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
  NMPP_uH_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
  RHS.reinit(locally_owned_dofs_LS,mpi_communicator);
  uStage1_nonGhosted.reinit (locally_owned_dofs_LS,mpi_communicator);
  uStage2_nonGhosted.reinit (locally_owned_dofs_LS,mpi_communicator);
  unp1.reinit (locally_owned_dofs_LS,mpi_communicator);
  MPP_uH_solution.reinit (locally_owned_dofs_LS,mpi_communicator);

vectors for lumped mass matrix

  ML_vector.reinit(locally_owned_dofs_LS,mpi_communicator);
  inverse_ML_vector.reinit(locally_owned_dofs_LS,mpi_communicator);
  ones_vector.reinit(locally_owned_dofs_LS,mpi_communicator);

operators times solution

  K_times_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
  DL_times_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
  DH_times_solution.reinit(locally_owned_dofs_LS,mpi_communicator);

LIMITERS (FCT)

  umin_vector.reinit (locally_owned_dofs_LS,mpi_communicator);
  umax_vector.reinit (locally_owned_dofs_LS,mpi_communicator);

///////////////////////////////////////////////////// GHOSTED VECTORS (used within some assemble process) /////////////////////////////////////////////////////

init vectors for vx

init vectors for vy

init vectors for vz

LIMITERS (FCT)

//////////////// SETUP MATRICES //////////////// MATRICES

  dof_handler_LS.locally_owned_dofs(),
  mpi_communicator,
  MC_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
  dof_handler_LS.locally_owned_dofs(),
  dsp, mpi_communicator);
  Cx_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
  dof_handler_LS.locally_owned_dofs(),
  dsp, mpi_communicator);
  CTx_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
  dof_handler_LS.locally_owned_dofs(),
  dsp, mpi_communicator);
  Cy_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
  dof_handler_LS.locally_owned_dofs(),
  dsp, mpi_communicator);
  CTy_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
  dof_handler_LS.locally_owned_dofs(),
  dsp, mpi_communicator);
  if (dim==3)
  {
  Cz_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
  dof_handler_LS.locally_owned_dofs(),
  dsp, mpi_communicator);
  CTz_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
  dof_handler_LS.locally_owned_dofs(),
  dsp, mpi_communicator);
  }
  dLij_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
  dof_handler_LS.locally_owned_dofs(),
  dsp, mpi_communicator);
  EntRes_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
  dof_handler_LS.locally_owned_dofs(),
  dsp, mpi_communicator);
  SuppSize_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
  dof_handler_LS.locally_owned_dofs(),
  dsp, mpi_communicator);
  dCij_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
  dof_handler_LS.locally_owned_dofs(),
  dsp, mpi_communicator);
  A_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
  dof_handler_LS.locally_owned_dofs(),
  dsp, mpi_communicator);
  LxA_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
  dof_handler_LS.locally_owned_dofs(),
  dsp, mpi_communicator);
  Akp1_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
  dof_handler_LS.locally_owned_dofs(),
  dsp, mpi_communicator);
  LxAkp1_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
  dof_handler_LS.locally_owned_dofs(),
  dsp, mpi_communicator);
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm mpi_comm, const IndexSet &locally_relevant_rows)

COMPUTE MASS MATRICES (AND OTHERS) FOR FIRST TIME STEP

get mat for DOFs between Q1 and Q2

  get_sparsity_pattern();
  }

---------------------------— MASS MATRICES ---------------------------—

  template<int dim>
  {
  const unsigned int dofs_per_cell = fe_LS.dofs_per_cell;
  const unsigned int n_q_points = quadrature_formula.size();
  Vector<double> cell_ML (dofs_per_cell);
  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
  cell_LS = dof_handler_LS.begin_active(),
  if (cell_LS->is_locally_owned())
  {
  cell_ML = 0;
  for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
  {
  const double JxW = fe_values_LS.JxW(q_point);
  for (unsigned int i=0; i<dofs_per_cell; ++i)
  cell_ML (i) += fe_values_LS.shape_value(i,q_point)*JxW;
  }
typename ActiveSelector::active_cell_iterator active_cell_iterator
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.

distribute

  cell_LS->get_dof_indices (local_dof_indices);
  constraints.distribute_local_to_global (cell_ML,local_dof_indices,ML_vector);
  }

compress

loop on locally owned i-DOFs (rows)

  {
  int gi = *idofs_iter;
  }
  }
  template<int dim>
  {
  const unsigned int dofs_per_cell = fe_LS.dofs_per_cell;
  const unsigned int n_q_points = quadrature_formula.size();
  FullMatrix<double> cell_MC (dofs_per_cell, dofs_per_cell);
  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
  std::vector<double> shape_values(dofs_per_cell);
  cell_LS = dof_handler_LS.begin_active(),
  if (cell_LS->is_locally_owned())
  {
  cell_MC = 0;
  for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
  {
  const double JxW = fe_values_LS.JxW(q_point);
  for (unsigned int i=0; i<dofs_per_cell; ++i)
  shape_values[i] = fe_values_LS.shape_value(i,q_point);
  for (unsigned int i=0; i<dofs_per_cell; ++i)
  for (unsigned int j=0; j<dofs_per_cell; ++j)
  cell_MC(i,j) += shape_values[i]*shape_values[j]*JxW;
  }

distribute

  cell_LS->get_dof_indices (local_dof_indices);
  constraints.distribute_local_to_global (cell_MC,local_dof_indices,MC_matrix);
  }

compress


---------------------------— LO METHOD (Dij Viscosity) ---------------------------—

  template <int dim>
  {
  const unsigned int dofs_per_cell_LS = fe_LS.dofs_per_cell;
  const unsigned int n_q_points = quadrature_formula.size();
  std::vector<Tensor<1, dim> > shape_grads_LS(dofs_per_cell_LS);
  std::vector<double> shape_values_LS(dofs_per_cell_LS);
  std::vector<types::global_dof_index> local_dof_indices_LS (dofs_per_cell_LS);
  cell_LS = dof_handler_LS.begin_active();
  if (cell_LS->is_locally_owned())
  {
  if (dim==3)
  {
  }
  cell_LS->get_dof_indices (local_dof_indices_LS);
  for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
  {
  const double JxW = fe_values_LS.JxW(q_point);
  for (unsigned int i=0; i<dofs_per_cell_LS; ++i)
  {
  shape_values_LS[i] = fe_values_LS.shape_value(i,q_point);
  shape_grads_LS [i] = fe_values_LS.shape_grad (i,q_point);
  }
  for (unsigned int i=0; i<dofs_per_cell_LS; ++i)
  for (unsigned int j=0; j < dofs_per_cell_LS; ++j)
  {
  if (dim==3)
  {
  }
  }
  }

Distribute

  constraints.distribute_local_to_global(cell_Cij_x,local_dof_indices_LS,Cx_matrix);
  constraints.distribute_local_to_global(cell_Cji_x,local_dof_indices_LS,CTx_matrix);
  constraints.distribute_local_to_global(cell_Cij_y,local_dof_indices_LS,Cy_matrix);
  constraints.distribute_local_to_global(cell_Cji_y,local_dof_indices_LS,CTy_matrix);
  if (dim==3)
  {
  constraints.distribute_local_to_global(cell_Cij_z,local_dof_indices_LS,Cz_matrix);
  constraints.distribute_local_to_global(cell_Cji_z,local_dof_indices_LS,CTz_matrix);
  }
  }

COMPRESS

  if (dim==3)
  {
  }
  }
  template<int dim>
  {
  const unsigned int dofs_per_cell = fe_LS.dofs_per_cell;
  const unsigned int n_q_points = quadrature_formula.size();
  std::vector<Tensor<1,dim> > un_grads (n_q_points);
  std::vector<double> old_vx_values (n_q_points);
  std::vector<double> old_vy_values (n_q_points);
  std::vector<double> old_vz_values (n_q_points);
  std::vector<double> shape_values(dofs_per_cell);
  std::vector<Tensor<1,dim> > shape_grads(dofs_per_cell);
  Vector<double> un_dofs(dofs_per_cell);
  std::vector<types::global_dof_index> indices_LS (dofs_per_cell);

loop on cells

  cell_LS = dof_handler_LS.begin_active(),
  cell_U = dof_handler_U.begin_active();
  if (cell_LS->is_locally_owned())
  {
  cell_LS->get_dof_indices (indices_LS);
  fe_values_LS.get_function_gradients(solution,un_grads);
  if (dim==3) fe_values_U.get_function_values(locally_relevant_solution_vz,old_vz_values);

compute cell_K_times_solution

  for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
  {
  if (dim==3) v[2] = old_vz_values[q_point]; //dim=3
  for (unsigned int i=0; i<dofs_per_cell; ++i)
  }

distribute

K_times_solution=0;

  const PetscInt *gj;
  const PetscScalar *Cxi, *Cyi, *Czi, *CTxi, *CTyi, *CTzi;
  const PetscScalar *EntResi, *SuppSizei, *MCi;
  double solni;

loop on locally owned i-DOFs (rows)

double ith_K_times_solution = 0;

read velocity of i-th DOF

get i-th row of C matrices

get vector values for column indices

Array for i-th row of matrices

  std::vector<double> dLi(ncolumns), dCi(ncolumns);
  double dLii = 0, dCii = 0;

loop on sparsity pattern of i-th DOF

  for (int j =0; j < ncolumns; ++j)
  {
  C[0] = Cxi[j];
  C[1] = Cyi[j];
  CT[0]= CTxi[j];
  CT[1]= CTyi[j];
  vj[0] = vx[j];
  vj[1] = vy[j];
  if (dim==3)
  {
  C[2] = Czi[j];
  CT[2] = CTzi[j];
  vj[2] = vz[j];
  }

ith_K_times_solution += soln[j]*(vj*C);

  if (gi!=gj[j])
  {

low order dissipative matrix

high order dissipative matrix (entropy viscosity)

  double dEij = -std::min(-dLi[j],
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)

high order compression matrix

  double Compij = cK*std::max(1-std::pow(0.5*(solni+soln[j]),2),0.0)/(std::abs(solni-soln[j])+1E-14);
  dCi[j] = dEij*std::max(1-Compij,0.0);
  dCii -= dCi[j];
  }
  }
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)

save K times solution vector K_times_solution(gi)=ith_K_times_solution; save i-th row of matrices on global matrices

  MatSetValuesRow(dLij_matrix,gi,&dLi[0]); // BTW: there is a dealii wrapper for this
  MatSetValuesRow(dCij_matrix,gi,&dCi[0]); // BTW: there is a dealii wrapper for this

Restore matrices after reading rows

compress K_times_solution.compress(VectorOperation::insert);

get matrices times vector

  }

---------------------------— ENTROPY VISCOSITY ---------------------------—

  template <int dim>
  {
  const unsigned int dofs_per_cell_LS = fe_LS.dofs_per_cell;
  const unsigned int n_q_points = quadrature_formula.size();
  std::vector<double> uqn (n_q_points); // un at q point
  std::vector<double> uqnm1 (n_q_points);
  std::vector<Tensor<1,dim> > guqn (n_q_points); //grad of uqn
  std::vector<Tensor<1,dim> > guqnm1 (n_q_points);
  std::vector<double> vxqn (n_q_points);
  std::vector<double> vyqn (n_q_points);
  std::vector<double> vzqn (n_q_points);
  std::vector<double> vxqnm1 (n_q_points);
  std::vector<double> vyqnm1 (n_q_points);
  std::vector<double> vzqnm1 (n_q_points);
  std::vector<Tensor<1, dim> > shape_grads_LS(dofs_per_cell_LS);
  std::vector<double> shape_values_LS(dofs_per_cell_LS);
  std::vector<types::global_dof_index> local_dof_indices_LS (dofs_per_cell_LS);
  cell_LS = dof_handler_LS.begin_active();
  cell_U = dof_handler_U.begin_active();
  double Rk;
  double cell_volume_double, volume=0;
  if (cell_LS->is_locally_owned())
  {

get solutions at quadrature points

  cell_LS->get_dof_indices (local_dof_indices_LS);
  fe_values_LS.get_function_values(un,uqn);
  fe_values_LS.get_function_values(unm1,uqnm1);
  fe_values_LS.get_function_gradients(un,guqn);
  fe_values_LS.get_function_gradients(unm1,guqnm1);
  if (dim==3) fe_values_U.get_function_values(locally_relevant_solution_vz,vzqn);
  if (dim==3) fe_values_U.get_function_values(locally_relevant_solution_vz_old,vzqnm1);
  for (unsigned int q=0; q<n_q_points; ++q)
  {
  if (dim==3)
  const double JxW = fe_values_LS.JxW(q);
  for (unsigned int i=0; i<dofs_per_cell_LS; ++i)
  {
  shape_values_LS[i] = fe_values_LS.shape_value(i,q);
  shape_grads_LS [i] = fe_values_LS.shape_grad (i,q);
  }
  for (unsigned int i=0; i<dofs_per_cell_LS; ++i)
  for (unsigned int j=0; j < dofs_per_cell_LS; ++j)
  {
  cell_volume (i,j) += JxW;
  }
  }
double volume(const Triangulation< dim, spacedim > &tria)

Distribute

  constraints.distribute_local_to_global(cell_EntRes,local_dof_indices_LS,EntRes_matrix);
  constraints.distribute_local_to_global(cell_volume,local_dof_indices_LS,SuppSize_matrix);
  }

ENTROPY NORM FACTOR

  volume = Utilities::MPI::sum(volume,mpi_communicator);
  entropy_mass = Utilities::MPI::sum(entropy_mass,mpi_communicator)/volume;
  }
T sum(const T &t, const MPI_Comm mpi_communicator)
T max(const T &t, const MPI_Comm mpi_communicator)
T min(const T &t, const MPI_Comm mpi_communicator)

---------------------------— TO CHECK MAX PRINCIPLE ---------------------------—

loop on locally owned i-DOFs (rows)

get solution at DOFs on the sparsity pattern of i-th DOF

  std::vector<types::global_dof_index> gj_indices = sparsity_pattern[gi];
  std::vector<double> soln(gj_indices.size());

compute bounds, ith row of flux matrix, P vectors

  double mini=1E10, maxi=-1E10;
  for (unsigned int j =0; j < gj_indices.size(); ++j)
  {

bounds

compute min and max vectors

  const unsigned int dofs_per_cell = fe_LS.dofs_per_cell;
  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
  double tol=1e-10;
  cell_LS = dof_handler_LS.begin_active(),
  if (cell_LS->is_locally_owned() && !cell_LS->at_boundary())
  {
  cell_LS->get_dof_indices(local_dof_indices);
  for (unsigned int i=0; i<dofs_per_cell; ++i)
  if (locally_owned_dofs_LS.is_element(local_dof_indices[i]))
  {
  double solni = unp1_solution(local_dof_indices[i]);
  if (solni - umin_vector(local_dof_indices[i]) < -tol || umax_vector(local_dof_indices[i]) - solni < -tol)
  {
  pcout << "MAX Principle violated" << std::endl;
  abort();
  }
  }
  }
  }

---------------------------— COMPUTE SOLUTIONS ---------------------------—

NON-GHOSTED VECTORS: MPP_uL_solution, NMPP_uH_solution GHOSTED VECTORS: un_solution

  NMPP_uH_solution=un_solution; // to start iterative solver at un_solution (instead of zero)

assemble RHS VECTORS

///////////////////////// COMPUTE MPP u1 solution /////////////////////////

////////////////////////////// COMPUTE GALERKIN u2 solution //////////////////////////////

loop on locally owned i-DOFs (rows)

read vectors at i-th DOF

get i-th row of matrices

get vector values for support of i-th DOF

  const std::vector<types::global_dof_index> gj_indices (gj,gj+ncolumns);
  std::vector<double> soln(ncolumns);
  std::vector<double> solH(ncolumns);

Array for i-th row of matrices

  std::vector<double> Ai(ncolumns);

compute bounds, ith row of flux matrix, P vectors

  double mini=1E10, maxi=-1E10;
  double Pposi=0 ,Pnegi=0;
  for (int j =0; j < ncolumns; ++j)
  {

bounds

i-th row of flux matrix A

  Ai[j] = (((gi==gj[j]) ? 1 : 0)*mi - MCi[j])*(solH[j]-soln[j] - (solHi-solni))

compute P vectors

  Pposi += Ai[j]*((Ai[j] > 0) ? 1. : 0.);
  Pnegi += Ai[j]*((Ai[j] < 0) ? 1. : 0.);
  }

save i-th row of flux matrix A

compute Q vectors

  double Qposi = mi*(maxi-solLi);
  double Qnegi = mi*(mini-solLi);

compute R vectors

  R_pos_vector_nonGhosted(gi) = ((Pposi==0) ? 1. : std::min(1.0,Qposi/Pposi));
  R_neg_vector_nonGhosted(gi) = ((Pnegi==0) ? 1. : std::min(1.0,Qnegi/Pnegi));

Restore matrices after reading rows

compress A matrix

compress R vectors

update ghost values for R vectors

compute limiters. NOTE: this is a different loop due to need of i- and j-th entries of R vectors

get i-th row of A matrix

get vector values for column indices

  const std::vector<types::global_dof_index> gj_indices (gj,gj+ncolumns);
  std::vector<double> Rpos(ncolumns);
  std::vector<double> Rneg(ncolumns);

Array for i-th row of A_times_L matrix

  std::vector<double> LxAi(ncolumns);

loop in sparsity pattern of i-th DOF

  for (int j =0; j < ncolumns; ++j)
  LxAi[j] = Ai[j] * ((Ai[j]>0) ? std::min(Rposi,Rneg[j]) : std::min(Rnegi,Rpos[j]));

save i-th row of LxA

  MatSetValuesRow(LxA_matrix,gi,&LxAi[0]); // BTW: there is a dealii wrapper for this

restore A matrix after reading it

loop in num of FCT iterations

  const PetscInt *gj;
  const PetscScalar *Akp1i;
  double mi;
  for (int iter=0; iter<NUM_ITER; ++iter)
  {
  Akp1_matrix.add(-1.0, LxAkp1_matrix); //new matrix to limit: A-LxA

loop on locally owned i-DOFs (rows)

read vectors at i-th DOF

get i-th row of matrices

get vector values for support of i-th DOF

  const std::vector<types::global_dof_index> gj_indices (gj,gj+ncolumns);
  std::vector<double> soln(ncolumns);

compute bounds, ith row of flux matrix, P vectors

  double mini=1E10, maxi=-1E10;
  double Pposi=0 ,Pnegi=0;
  for (int j =0; j < ncolumns; ++j)
  {

bounds

compute P vectors

  Pposi += Akp1i[j]*((Akp1i[j] > 0) ? 1. : 0.);
  Pnegi += Akp1i[j]*((Akp1i[j] < 0) ? 1. : 0.);
  }

compute Q vectors

  double Qposi = mi*(maxi-solLi);
  double Qnegi = mi*(mini-solLi);

compute R vectors

  R_pos_vector_nonGhosted(gi) = ((Pposi==0) ? 1. : std::min(1.0,Qposi/Pposi));
  R_neg_vector_nonGhosted(gi) = ((Pnegi==0) ? 1. : std::min(1.0,Qnegi/Pnegi));

Restore matrices after reading rows

compress R vectors

update ghost values for R vectors

compute limiters. NOTE: this is a different loop due to need of i- and j-th entries of R vectors

get i-th row of Akp1 matrix

get vector values for column indices

  const std::vector<types::global_dof_index> gj_indices(gj,gj+ncolumns);
  std::vector<double> Rpos(ncolumns);
  std::vector<double> Rneg(ncolumns);

Array for i-th row of LxAkp1 matrix

  std::vector<double> LxAkp1i(ncolumns);
  for (int j =0; j < ncolumns; ++j)
  LxAkp1i[j] = Akp1i[j] * ((Akp1i[j]>0) ? std::min(Rposi,Rneg[j]) : std::min(Rnegi,Rpos[j]));

save i-th row of LxA

  MatSetValuesRow(LxAkp1_matrix,gi,&LxAkp1i[0]); // BTW: there is a dealii wrapper for this

restore A matrix after reading it

COMPUTE MPP LOW-ORDER SOLN and NMPP HIGH-ORDER SOLN

GHOSTED VECTORS: un NON-GHOSTED VECTORS: unp1

///////////// FIRST STAGE ///////////// u1=un-dt*RH*un

////////////// SECOND STAGE ////////////// u2=3/4*un+1/4*(u1-dt*RH*u1)

///////////// THIRD STAGE ///////////// unp1=1/3*un+2/3*(u2-dt*RH*u2)


---------------------------— UTILITIES ---------------------------—

loop on DOFs

get i-th row of mass matrix (dummy, I just need the indices gj)

  sparsity_pattern[gi] = std::vector<types::global_dof_index>(gj,gj+ncolumns);
  }
  }
  template<int dim>
  {
  const unsigned int dofs_per_cell_LS = fe_LS.dofs_per_cell;
  std::vector<types::global_dof_index> local_dof_indices_LS (dofs_per_cell_LS);
  const unsigned int dofs_per_cell_U = fe_U.dofs_per_cell;
  std::vector<types::global_dof_index> local_dof_indices_U (dofs_per_cell_U);
  cell_LS = dof_handler_LS.begin_active(),
  cell_U = dof_handler_U.begin_active();
  if (!cell_LS->is_artificial()) // loop on ghost cells as well
  {
  cell_LS->get_dof_indices(local_dof_indices_LS);
  cell_U->get_dof_indices(local_dof_indices_U);
  for (unsigned int i=0; i<dofs_per_cell_LS; ++i)
  }
  }
  template <int dim>
  std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
  {

all vectors are NON-GHOSTED

  SolverControl solver_control (dof_handler_LS.n_dofs(), solver_tolerance);
  PETScWrappers::SolverCG solver(solver_control, mpi_communicator);
  constraints.distribute (completely_distributed_solution);
  solver.solve (Matrix, completely_distributed_solution, rhs, *preconditioner);
  constraints.distribute (completely_distributed_solution);
  if (verbose==true) pcout << " Solved in " << solver_control.last_step() << " iterations." << std::endl;
  }
  template <int dim>
  {
  }
  template <int dim>
  {
  if (dim==3)
  }

---------------------------— MY PETSC WRAPPERS ---------------------------—

  template<int dim>
  const std::vector<types::global_dof_index> &indices,
  std::vector<PetscScalar> &values)
  {

PETSc wrapper to get sets of values from a petsc vector. we assume the vector is ghosted We need to figure out which elements we own locally. Then get a pointer to the elements that are stored here (both the ones we own as well as the ghost elements). In this array, the locally owned elements come first followed by the ghost elements whose position we can get from an index set

  ghost_indices.subtract_set(locally_owned_dofs_LS);
  PetscInt n_idx, begin, end, i;
  n_idx = indices.size();
  VecGetOwnershipRange (vector, &begin, &end);
  Vec solution_in_local_form = nullptr;
  PetscScalar *soln;
  for (i = 0; i < n_idx; ++i)
  {
  int index = indices[i];
  if (index >= begin && index < end)
  values[i] = *(soln+index-begin);
  else //ghost
  {
  const unsigned int ghostidx = ghost_indices.index_within_set(index);
  values[i] = *(soln+ghostidx+end-begin);
  }
  }
  }
  template<int dim>
  const std::vector<types::global_dof_index> &indices,
  std::map<types::global_dof_index, types::global_dof_index> &map_from_Q1_to_Q2,
  std::vector<PetscScalar> &values)
  {

THIS IS MEANT TO BE USED WITH VELOCITY VECTORS PETSc wrapper to get sets of values from a petsc vector. we assume the vector is ghosted We need to figure out which elements we own locally. Then get a pointer to the elements that are stored here (both the ones we own as well as the ghost elements). In this array, the locally owned elements come first followed by the ghost elements whose position we can get from an index set

  ghost_indices.subtract_set(locally_owned_dofs_U);
  PetscInt n_idx, begin, end, i;
  n_idx = indices.size();
  VecGetOwnershipRange (vector, &begin, &end);
  Vec solution_in_local_form = nullptr;
  PetscScalar *soln;
  for (i = 0; i < n_idx; ++i)
  {
  int index = map_from_Q1_to_Q2[indices[i]];
  if (index >= begin && index < end)
  values[i] = *(soln+index-begin);
  else //ghost
  {
  const unsigned int ghostidx = ghost_indices.index_within_set(index);
  values[i] = *(soln+ghostidx+end-begin);
  }
  }
  }

Annotated version of MultiPhase.cc

  /* -----------------------------------------------------------------------------
  *
  * SPDX-License-Identifier: LGPL-2.1-or-later
  * Copyright (C) 2016 Manuel Quezada de Luna
  *
  * This file is part of the deal.II code gallery.
  *
  * -----------------------------------------------------------------------------
  */
  #include <deal.II/base/function.h>
  #include <deal.II/lac/affine_constraints.h>
  #include <deal.II/lac/vector.h>
  #include <deal.II/lac/petsc_vector.h>
  #include <deal.II/dofs/dof_handler.h>
  #include <deal.II/dofs/dof_tools.h>
  #include <deal.II/fe/fe_values.h>
  #include <deal.II/fe/fe_q.h>
  #include <deal.II/numerics/vector_tools.h>
  #include <deal.II/numerics/data_out.h>
  #include <deal.II/numerics/error_estimator.h>
  #include <deal.II/base/utilities.h>
  #include <deal.II/base/index_set.h>
  #include <deal.II/distributed/tria.h>
  #include <deal.II/distributed/grid_refinement.h>
  #include <deal.II/base/timer.h>
  #include <deal.II/grid/grid_tools.h>
  #include <deal.II/fe/mapping_q.h>
  using namespace dealii;

/////////////////////// FOR TRANSPORT PROBLEM /////////////////////// TIME_INTEGRATION

  #define FORWARD_EULER 0

PROBLEM

  #include "NavierStokesSolver.cc"
  #include "LevelSetSolver.cc"
  #include "utilities.cc"

/////////////////////////////////////////////////// /////////////////// MAIN CLASS //////////////////// ///////////////////////////////////////////////////

SOLUTION VECTORS

BOUNDARY VECTORS

  std::vector<types::global_dof_index> boundary_values_id_u;
  std::vector<types::global_dof_index> boundary_values_id_v;
  std::vector<types::global_dof_index> boundary_values_id_phi;
  std::vector<double> boundary_values_u;
  std::vector<double> boundary_values_v;
  std::vector<double> boundary_values_phi;
  double time;
  double time_step;
  double final_time;
  unsigned int timestep_number;
  double cfl;
  double umax;
  double min_h;
  double sharpness;
  unsigned int n_refinement;
  unsigned int output_number;
  double output_time;
  bool verbose;

FOR NAVIER STOKES

  double rho_fluid;
  double nu_fluid;
  double rho_air;
  double nu_air;
  double nu;
  double eps;

FOR TRANSPORT

  double cK; //compression coeff
  double cE; //entropy-visc coeff
  std::string ALGORITHM;
  unsigned int PROBLEM;
  };
  template <int dim>
  const unsigned int degree_U)
  :
  mpi_communicator (MPI_COMM_WORLD),
  triangulation (mpi_communicator,
  typename Triangulation<dim>::MeshSmoothing
  (Triangulation<dim>::smoothing_on_refinement |
  Triangulation<dim>::smoothing_on_coarsening)),
  pcout (std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator)== 0))
  {}
  {
  dof_handler_LS.clear ();
  }

///////////////////////////////////// /////////////// SETUP /////////////// /////////////////////////////////////

  template <int dim>
  {

setup system LS

setup system U

setup system P

init vectors for phi

init vectors for u

init vectors for v

init vectors for p

INIT CONSTRAINTS

  }
  template <int dim>
  {
  time=0;

Initial conditions init condition for phi

  constraints.distribute (completely_distributed_solution_phi);
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask={}, const unsigned int level=numbers::invalid_unsigned_int)

init condition for u=0

init condition for v

init condition for p

  constraints.distribute (completely_distributed_solution_p);
  }
  template <int dim>
  {
  constraints.clear ();
  constraints.reinit (locally_relevant_dofs_LS);
  constraints.close ();
  }
  template <int dim>
  {
  std::map<types::global_dof_index, double> map_boundary_values_u;
  std::map<types::global_dof_index, double> map_boundary_values_v;
  std::map<types::global_dof_index, double> map_boundary_values_w;

NO-SLIP CONDITION

LEFT

void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})

RIGHT

BOTTOM

TOP

no slip in bottom and top and slip in left and right LEFT

RIGHT

BOTTOM

TOP

LEFT: entry in x, zero in y

RIGHT: no-slip condition

BOTTOM: non-slip

TOP: exit in y, zero in x

  }
  else
  {
  pcout << "Error in type of PROBLEM at Boundary Conditions" << std::endl;
  abort();
  }
  std::map<types::global_dof_index,double>::const_iterator boundary_value_u =map_boundary_values_u.begin();
  std::map<types::global_dof_index,double>::const_iterator boundary_value_v =map_boundary_values_v.begin();
  {
  }
  {
  }
  }
  template <int dim>
  {
  const QGauss<dim-1> face_quadrature_formula(1); // center of the face
  const unsigned int n_face_q_points = face_quadrature_formula.size();
  std::vector<double> u_value (n_face_q_points);
  std::vector<double> v_value (n_face_q_points);
  cell_U = dof_handler_U.begin_active(),
  for (; cell_U!=endc_U; ++cell_U)
  if (cell_U->is_locally_owned())
  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
  if (cell_U->face(face)->at_boundary())
  {
  fe_face_values.reinit(cell_U,face);
  fe_face_values.get_function_values(locally_relevant_solution_u,u_value);
  fe_face_values.get_function_values(locally_relevant_solution_v,v_value);
  u[0]=u_value[0];
  u[1]=v_value[0];
  if (fe_face_values.normal_vector(0)*u < -1e-14)
  cell_U->face(face)->set_boundary_id(10); // SET ID 10 to inlet BOUNDARY (10 is an arbitrary number)
  }
  }
  template <int dim>
  void MultiPhase<dim>::get_boundary_values_phi(std::vector<types::global_dof_index> &boundary_values_id_phi,
  std::vector<double> &boundary_values_phi)
  {
  std::map<types::global_dof_index, double> map_boundary_values_phi;
  unsigned int boundary_id=0;
  boundary_id=10; // inlet
  std::map<types::global_dof_index,double>::const_iterator boundary_value_phi = map_boundary_values_phi.begin();
  {
  }
  }
  template<int dim>
  {
@ update_normal_vectors
Normal vectors.
unsigned int boundary_id
Definition types.h:161

output_vectors();

  }
  template <int dim>
  {
  DataOut<dim> data_out;
  data_out.attach_dof_handler (dof_handler_LS);
  data_out.add_data_vector (locally_relevant_solution_phi, "phi");
  data_out.build_patches ();
  const std::string filename = ("sol_vectors-" +
  "." +
  (triangulation.locally_owned_subdomain(), 4));
  std::ofstream output ((filename + ".vtu").c_str());
  data_out.write_vtu (output);
  if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
  {
  std::vector<std::string> filenames;
  for (unsigned int i=0;
  i<Utilities::MPI::n_mpi_processes(mpi_communicator);
  ++i)
  filenames.push_back ("sol_vectors-" +
  "." +
  ".vtu");
  std::ofstream master_output ((filename + ".pvtu").c_str());
  data_out.write_pvtu_record (master_output, filenames);
  }
  }
  template <int dim>
  {
  DataOut<dim> data_out;
  data_out.attach_dof_handler (dof_handler_LS);
  data_out.add_data_vector (locally_relevant_solution_phi, postprocessor);
  data_out.build_patches ();
  const std::string filename = ("sol_rho-" +
  "." +
  (triangulation.locally_owned_subdomain(), 4));
  std::ofstream output ((filename + ".vtu").c_str());
  data_out.write_vtu (output);
  if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
  {
  std::vector<std::string> filenames;
  for (unsigned int i=0;
  i<Utilities::MPI::n_mpi_processes(mpi_communicator);
  ++i)
  filenames.push_back ("sol_rho-" +
  "." +
  ".vtu");
  std::ofstream master_output ((filename + ".pvtu").c_str());
  data_out.write_pvtu_record (master_output, filenames);
  }
  }
  template <int dim>
  {
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:99
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:114
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:466

//////////////////// GENERAL PARAMETERS ////////////////////

  umax=1;
  cfl=0.1;
  verbose = true;
  get_output = true;
  final_time = 10.0;

////////////////////////////////////////// PARAMETERS FOR THE NAVIER STOKES PROBLEM //////////////////////////////////////////

  rho_fluid = 1000.;
  nu_fluid = 1.0;
  rho_air = 1.0;
  nu_air = 1.8e-2;

PROBLEM=FILLING_TANK; PROBLEM=SMALL_WAVE_PERTURBATION; PROBLEM=FALLING_DROP;

  ForceTerms<dim> force_function(std::vector<double> {0.0,-1.0});

////////////////////////////////// PARAMETERS FOR TRANSPORT PROBLEM //////////////////////////////////

  cK = 1.0;
  cE = 1.0;
  sharpness_integer=10; //this will be multiplied by min_h

TRANSPORT_TIME_INTEGRATION=FORWARD_EULER;

ALGORITHM = "MPP_u1"; ALGORITHM = "NMPP_uH";

  ALGORITHM = "MPP_uH";

ADJUST PARAMETERS ACCORDING TO PROBLEM

////////// GEOMETRY //////////

  Point<dim>(0.0,0.0), Point<dim>(0.4,0.4), true);
  {
  std::vector< unsigned int > repetitions;
  repetitions.push_back(2);
  repetitions.push_back(1);
  (triangulation, repetitions, Point<dim>(0.0,0.0), Point<dim>(1.0,0.5), true);
  }
  {
  std::vector< unsigned int > repetitions;
  repetitions.push_back(1);
  repetitions.push_back(4);
  (triangulation, repetitions, Point<dim>(0.0,0.0), Point<dim>(0.3,0.9), true);
  }
  triangulation.refine_global (n_refinement);
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)

SETUP

  setup();

PARAMETERS FOR TIME STEPPING

  eps=1.*min_h; //For reconstruction of density in Navier Stokes
  sharpness=sharpness_integer*min_h; //adjust value of sharpness (for init cond of phi)
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)

INITIAL CONDITIONS

NAVIER STOKES SOLVER

BOUNDARY CONDITIONS FOR NAVIER STOKES

set INITIAL CONDITION within NAVIER STOKES

TRANSPORT SOLVER

BOUNDARY CONDITIONS FOR PHI

set INITIAL CONDITION within TRANSPORT PROBLEM

NO BOUNDARY CONDITIONS for LEVEL SET

  pcout << "Cfl: " << cfl << "; umax: " << umax << "; min h: " << min_h
  << "; time step: " << time_step << std::endl;
  pcout << " Number of active cells: "
  << triangulation.n_global_active_cells() << std::endl
  << " Number of degrees of freedom: " << std::endl
  << " U: " << dofs_U << std::endl
  << " P: " << dofs_P << std::endl
  << " LS: " << dofs_LS << std::endl
  << " TOTAL: " << dofs_TOTAL
  << std::endl;

TIME STEPPING

  for (timestep_number=1, time=time_step; time<=final_time;
  {
  pcout << "Time step " << timestep_number
  << " at t=" << time
  << std::endl;

GET NAVIER STOKES VELOCITY

GET LEVEL SET SOLUTION

  transport_solver.nth_time_step();
  }
  }
  int main(int argc, char *argv[])
  {
  try
  {
  using namespace dealii;
  PetscInitialize(&argc, &argv, nullptr, nullptr);
  {
  unsigned int degree_LS = 1;
  unsigned int degree_U = 2;
  }
  }
  catch (std::exception &exc)
  {
  std::cerr << std::endl << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Exception on processing: " << std::endl
  << exc.what() << std::endl
  << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
  return 1;
  }
  catch (...)
  {
  std::cerr << std::endl << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Unknown exception!" << std::endl
  << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
  return 1;
  }
  return 0;
  }
unsigned int depth_console(const unsigned int n)
Definition logstream.cc:351
LogStream deallog
Definition logstream.cc:38

Annotated version of NavierStokesSolver.cc

  /* -----------------------------------------------------------------------------
  *
  * SPDX-License-Identifier: LGPL-2.1-or-later
  * Copyright (C) 2016 Manuel Quezada de Luna
  *
  * This file is part of the deal.II code gallery.
  *
  * -----------------------------------------------------------------------------
  */
  #include <deal.II/base/function.h>
  #include <deal.II/lac/affine_constraints.h>
  #include <deal.II/lac/vector.h>
  #include <deal.II/lac/petsc_vector.h>
  #include <deal.II/dofs/dof_handler.h>
  #include <deal.II/dofs/dof_tools.h>
  #include <deal.II/fe/fe_values.h>
  #include <deal.II/fe/fe_q.h>
  #include <deal.II/numerics/vector_tools.h>
  #include <deal.II/numerics/data_out.h>
  #include <deal.II/numerics/error_estimator.h>
  #include <deal.II/base/utilities.h>
  #include <deal.II/base/index_set.h>
  #include <deal.II/distributed/tria.h>
  #include <deal.II/distributed/grid_refinement.h>
  #include <deal.II/lac/petsc_vector.h>
  #include <deal.II/base/timer.h>
  #include <deal.II/grid/grid_tools.h>
  #include <deal.II/fe/mapping_q.h>
  using namespace dealii;

///////////////////////////////////////////////////////////// /////////////////// NAVIER STOKES SOLVER //////////////////// /////////////////////////////////////////////////////////////

  template<int dim>
  {
  public:

constructor for using LEVEL SET

  NavierStokesSolver(const unsigned int degree_LS,
  const unsigned int degree_U,
  const double time_step,
  const double eps,
  const double rho_air,
  const double nu_air,
  const double rho_fluid,
  const double nu_fluid,
  const bool verbose,
  MPI_Comm &mpi_communicator);

constructor for NOT LEVEL SET

rho and nu functions

initial conditions

boundary conditions

DO STEPS

SETUP

  void setup();
  private:

SETUP AND INITIAL CONDITION

  void setup_DOF();

ASSEMBLE SYSTEMS

SOLVERS

  std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
  std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,

GET DIFFERENT FIELDS

  void get_rho_and_nu(double phi);

OTHERS

  MPI_Comm &mpi_communicator;
  double rho_air;
  double nu_air;
  double rho_fluid;
  double nu_fluid;
  double time_step;
  double eps;
  bool verbose;
  unsigned int LEVEL_SET;
  unsigned int RHO_TIMES_RHS;
  double rho_min;
  double rho_value;
  double nu_value;
  double h;
  double umax;
  std::vector<types::global_dof_index> boundary_values_id_u;
  std::vector<types::global_dof_index> boundary_values_id_v;
  std::vector<types::global_dof_index> boundary_values_id_w;
  std::vector<double> boundary_values_u;
  std::vector<double> boundary_values_v;
  std::vector<double> boundary_values_w;
  std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_Matrix_u;
  std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_Matrix_v;
  std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_Matrix_w;
  std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_S;
  std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_M;
  };

CONSTRUCTOR FOR LEVEL SET

  template<int dim>
  const unsigned int degree_U,
  const double time_step,
  const double eps,
  const double rho_air,
  const double nu_air,
  const double rho_fluid,
  const double nu_fluid,
  const bool verbose,
  MPI_Comm &mpi_communicator)
  :
  mpi_communicator(mpi_communicator),

This is dummy since rho and nu functions won't be used

CONSTRUCTOR NOT FOR LEVEL SET

///////////////////////////////////////////////////////// ////////////////// SETTERS AND GETTERS ////////////////// /////////////////////////////////////////////////////////

set old vectors to the initial condition (just for first time step)

set old vectors to the initial condition (just for first time step)

  }
  template<int dim>
  std::vector<types::global_dof_index> boundary_values_id_v,
  std::vector<double> boundary_values_u,
  std::vector<double> boundary_values_v)
  {
  this->boundary_values_id_u=boundary_values_id_u;
  this->boundary_values_id_v=boundary_values_id_v;
  this->boundary_values_u=boundary_values_u;
  this->boundary_values_v=boundary_values_v;
  }
  template<int dim>
  std::vector<types::global_dof_index> boundary_values_id_v,
  std::vector<types::global_dof_index> boundary_values_id_w,
  std::vector<double> boundary_values_u,
  std::vector<double> boundary_values_v,
  std::vector<double> boundary_values_w)
  {
  this->boundary_values_id_u=boundary_values_id_u;
  this->boundary_values_id_v=boundary_values_id_v;
  this->boundary_values_id_w=boundary_values_id_w;
  this->boundary_values_u=boundary_values_u;
  this->boundary_values_v=boundary_values_v;
  this->boundary_values_w=boundary_values_w;
  }
  template<int dim>
  {
  this->locally_relevant_solution_u=locally_relevant_solution_u;
  this->locally_relevant_solution_v=locally_relevant_solution_v;
  }
  template<int dim>
  {
  this->locally_relevant_solution_u=locally_relevant_solution_u;
  this->locally_relevant_solution_v=locally_relevant_solution_v;
  this->locally_relevant_solution_w=locally_relevant_solution_w;
  }
  template<int dim>
  {
  this->locally_relevant_solution_phi=locally_relevant_solution_phi;
  }
  template<int dim>
  {
  double H=0;

get rho, nu

  if (phi>eps)
  H=1;
  else if (phi<-eps)
  H=-1;
  else
  H=phi/eps;
  rho_value=rho_fluid*(1+H)/2.+rho_air*(1-H)/2.;
  nu_value=nu_fluid*(1+H)/2.+nu_air*(1-H)/2.;

rho_value=rho_fluid*(1+phi)/2.+rho_air*(1-phi)/2.; nu_value=nu_fluid*(1+phi)/2.+nu_air*(1-phi)/2.;

/////////////////////////////////////////////////// /////////// SETUP AND INITIAL CONDITION /////////// ///////////////////////////////////////////////////

  template<int dim>
  {
  pcout<<"***** SETUP IN NAVIER STOKES SOLVER *****"<<std::endl;
  }
  template<int dim>
  {
  rho_min = 1.;

setup system LS

setup system U

setup system P

init vectors for phi

init vectors for u

init vectors for u_old

init vectors for v

init vectors for v_old

init vectors for w

init vectors for w_old

init vectors for dpsi

init vectors for dpsi old

init vectors for q

init vectors for psi

init vectors for p

//////////////////////// Initialize constraints ////////////////////////

////////////////// Sparsity pattern ////////////////// sparsity pattern for A

  dof_handler_U.locally_owned_dofs(),
  mpi_communicator,
  system_Matrix_u.reinit(dof_handler_U.locally_owned_dofs(),
  dof_handler_U.locally_owned_dofs(),
  mpi_communicator);
  system_Matrix_v.reinit(dof_handler_U.locally_owned_dofs(),
  dof_handler_U.locally_owned_dofs(),
  mpi_communicator);
  system_Matrix_w.reinit(dof_handler_U.locally_owned_dofs(),
  dof_handler_U.locally_owned_dofs(),
  mpi_communicator);

sparsity pattern for S

sparsity pattern for M

grl constraints

  constraints.clear();
  constraints.reinit(locally_relevant_dofs_U);
  constraints.close();

constraints for dpsi

if (constraints_psi.can_store_line(0)) constraints_psi.add_line(0); //constraint u0 = 0

  }

/////////////////////////////////////////////////// //////////////// ASSEMBLE SYSTEMS ///////////////// ///////////////////////////////////////////////////

  template<int dim>
  {
  if (rebuild_Matrix_U==true)
  {
  }
  const unsigned int dofs_per_cell=fe_U.dofs_per_cell;
  const unsigned int n_q_points=quadrature_formula.size();
  FullMatrix<double> cell_A_u(dofs_per_cell,dofs_per_cell);
  Vector<double> cell_rhs_u(dofs_per_cell);
  Vector<double> cell_rhs_v(dofs_per_cell);
  Vector<double> cell_rhs_w(dofs_per_cell);
  std::vector<double> phiqnp1(n_q_points);
  std::vector<double> uqn(n_q_points);
  std::vector<double> uqnm1(n_q_points);
  std::vector<double> vqn(n_q_points);
  std::vector<double> vqnm1(n_q_points);
  std::vector<double> wqn(n_q_points);
  std::vector<double> wqnm1(n_q_points);

FOR Explicit nonlinearity std::vector<Tensor<1, dim> > grad_un(n_q_points); std::vector<Tensor<1, dim> > grad_vn(n_q_points); std::vector<Tensor<1, dim> > grad_wn(n_q_points); Tensor<1, dim> Un;

  std::vector<Tensor<1, dim> > grad_pqn(n_q_points);
  std::vector<Tensor<1, dim> > grad_psiqn(n_q_points);
  std::vector<Tensor<1, dim> > grad_psiqnm1(n_q_points);
  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
  std::vector<Tensor<1, dim> > shape_grad(dofs_per_cell);
  std::vector<double> shape_value(dofs_per_cell);
  double force_u;
  double force_v;
  double force_w;
  double u_star=0;
  double v_star=0;
  double w_star=0;
  double rho_star;
  double rho;
  cell_U=dof_handler_U.begin_active(), endc_U=dof_handler_U.end();
  if (cell_U->is_locally_owned())
  {

get function values for LS

get function values for U

For explicit nonlinearity get gradient values for U fe_values_U.get_function_gradients(locally_relevant_solution_u,grad_un); fe_values_U.get_function_gradients(locally_relevant_solution_v,grad_vn); if (dim==3) fe_values_U.get_function_gradients(locally_relevant_solution_w,grad_wn);

get values and gradients for p and dpsi

  for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
  {
  const double JxW=fe_values_U.JxW(q_point);
  for (unsigned int i=0; i<dofs_per_cell; ++i)
  {
  shape_grad[i]=fe_values_U.shape_grad(i,q_point);
  shape_value[i]=fe_values_U.shape_value(i,q_point);
  }
  if (dim==3)
  if (LEVEL_SET==1) // use level set to define rho and nu
  else // rho and nu are defined through functions
  {
  rho_value=rho_function.value(fe_values_U.quadrature_point(q_point));
  nu_value=nu_function.value(fe_values_U.quadrature_point(q_point));
  }

Non-linearity: for semi-implicit

for explicit nonlinearity Un[0] = uqn[q_point]; Un[1] = vqn[q_point]; if (dim==3) Un[2] = wqn[q_point];

double nonlinearity_u = Un*grad_un[q_point]; double nonlinearity_v = Un*grad_vn[q_point]; double nonlinearity_w = 0; if (dim==3) nonlinearity_w = Un*grad_wn[q_point];

  rho_star=rho_value; // This is because we consider rho*u_t instead of (rho*u)_t

FORCE TERMS

  force_function.vector_value(fe_values_U.quadrature_point(q_point),force_terms);
  if (dim==3)
  {
  if (dim==3)
  }
  for (unsigned int i=0; i<dofs_per_cell; ++i)
  {

-2./3*time_step*rho*nonlinearity_u

  )*shape_value[i])*JxW;

-2./3*time_step*rho*nonlinearity_v

  )*shape_value[i])*JxW;
  if (dim==3)

-2./3*time_step*rho*nonlinearity_w

  )*shape_value[i])*JxW;
  if (rebuild_Matrix_U==true)
  for (unsigned int j=0; j<dofs_per_cell; ++j)
  {
  if (dim==2)
  cell_A_u(i,j)+=(rho_star*shape_value[i]*shape_value[j]
  +2./3*time_step*nu_value*(shape_grad[i]*shape_grad[j])
  +2./3*time_step*rho*shape_value[i]
  *(u_star*shape_grad[j][0]+v_star*shape_grad[j][1]) // semi-implicit NL
  )*JxW;
  else //dim==3
  cell_A_u(i,j)+=(rho_star*shape_value[i]*shape_value[j]
  +2./3*time_step*nu_value*(shape_grad[i]*shape_grad[j])
  +2./3*time_step*rho*shape_value[i]
  *(u_star*shape_grad[j][0]+v_star*shape_grad[j][1]+w_star*shape_grad[j][2]) // semi-implicit NL
  )*JxW;
  }
  }
  }
  cell_U->get_dof_indices(local_dof_indices);

distribute

  if (rebuild_Matrix_U==true)
  constraints.distribute_local_to_global(cell_A_u,local_dof_indices,system_Matrix_u);
  constraints.distribute_local_to_global(cell_rhs_u,local_dof_indices,system_rhs_u);
  constraints.distribute_local_to_global(cell_rhs_v,local_dof_indices,system_rhs_v);
  if (dim==3)
  constraints.distribute_local_to_global(cell_rhs_w,local_dof_indices,system_rhs_w);
  }
  if (dim==3) system_rhs_w.compress(VectorOperation::add);
  if (rebuild_Matrix_U==true)
  {
  if (dim==3)
  }

BOUNDARY CONDITIONS

PRECONDITIONERS

  if (dim==3)
  }
  }
  }
  template<int dim>
  {
  if (rebuild_S_M==true)
  {
  }
  const unsigned int dofs_per_cell=fe_P.dofs_per_cell;
  const unsigned int n_q_points=quadrature_formula.size();
  FullMatrix<double> cell_S(dofs_per_cell,dofs_per_cell);
  FullMatrix<double> cell_M(dofs_per_cell,dofs_per_cell);
  Vector<double> cell_rhs_q(dofs_per_cell);
  std::vector<double> phiqnp1(n_q_points);
  std::vector<Tensor<1, dim> > gunp1(n_q_points);
  std::vector<Tensor<1, dim> > gvnp1(n_q_points);
  std::vector<Tensor<1, dim> > gwnp1(n_q_points);
  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
  std::vector<double> shape_value(dofs_per_cell);
  std::vector<Tensor<1, dim> > shape_grad(dofs_per_cell);
  cell_P=dof_handler_P.begin_active(), endc_P=dof_handler_P.end();
  if (cell_P->is_locally_owned())
  {

get function values for LS

get function grads for u and v

  if (dim==3)
  for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
  {
  const double JxW=fe_values_P.JxW(q_point);
  double divU = gunp1[q_point][0]+gvnp1[q_point][1];
  if (dim==3) divU += gwnp1[q_point][2];
  for (unsigned int i=0; i<dofs_per_cell; ++i)
  {
  shape_value[i]=fe_values_P.shape_value(i,q_point);
  shape_grad[i]=fe_values_P.shape_grad(i,q_point);
  }
  if (LEVEL_SET==1) // use level set to define rho and nu
  else // rho and nu are defined through functions
  nu_value=nu_function.value(fe_values_U.quadrature_point(q_point));
  for (unsigned int i=0; i<dofs_per_cell; ++i)
  {
  cell_rhs_psi(i)+=-3./2./time_step*rho_min*divU*shape_value[i]*JxW;
  cell_rhs_q(i)-=nu_value*divU*shape_value[i]*JxW;
  if (rebuild_S_M==true)
  {
  for (unsigned int j=0; j<dofs_per_cell; ++j)
  if (i==j)
  {
  cell_S(i,j)+=shape_grad[i]*shape_grad[j]*JxW+1E-10;
  cell_M(i,j)+=shape_value[i]*shape_value[j]*JxW;
  }
  else
  {
  cell_S(i,j)+=shape_grad[i]*shape_grad[j]*JxW;
  cell_M(i,j)+=shape_value[i]*shape_value[j]*JxW;
  }
  }
  }
  }
  cell_P->get_dof_indices(local_dof_indices);

Distribute

/////////////////////////////////////////////////// ///////////////////// SOLVERS ///////////////////// ///////////////////////////////////////////////////

  template<int dim>
  std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
  {
  SolverControl solver_control(dof_handler_U.n_dofs(),1e-6);

PETScWrappers::SolverCG solver(solver_control, mpi_communicator); PETScWrappers::SolverGMRES solver(solver_control, mpi_communicator); PETScWrappers::SolverChebychev solver(solver_control, mpi_communicator);

  PETScWrappers::SolverBicgstab solver(solver_control,mpi_communicator);
  constraints.distribute(completely_distributed_solution);
  solver.solve(Matrix,completely_distributed_solution,rhs,*preconditioner);
  constraints.distribute(completely_distributed_solution);
  if (solver_control.last_step() > MAX_NUM_ITER_TO_RECOMPUTE_PRECONDITIONER)
  if (verbose==true)
  pcout<<" Solved U in "<<solver_control.last_step()<<" iterations."<<std::endl;
  }
  template<int dim>
  std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
  {
  SolverControl solver_control(dof_handler_P.n_dofs(),1e-6);
  PETScWrappers::SolverCG solver(solver_control,mpi_communicator);

PETScWrappers::SolverGMRES solver(solver_control, mpi_communicator);

  constraints.distribute(completely_distributed_solution);
  solver.solve(Matrix,completely_distributed_solution,rhs,*preconditioner);
  constraints.distribute(completely_distributed_solution);
  if (solver_control.last_step() > MAX_NUM_ITER_TO_RECOMPUTE_PRECONDITIONER)
  if (verbose==true)
  pcout<<" Solved P in "<<solver_control.last_step()<<" iterations."<<std::endl;
  }

/////////////////////////////////////////////////// ////////////// get different fields /////////////// ///////////////////////////////////////////////////

GET DPSI

SOLVE Q

UPDATE THE PRESSURE

/////////////////////////////////////////////////// ///////////////////// DO STEPS //////////////////// ///////////////////////////////////////////////////

/////////////////////////////////////////////////// ////////////////////// OTHERS ///////////////////// ///////////////////////////////////////////////////

Annotated version of TestLevelSet.cc

  /* -----------------------------------------------------------------------------
  *
  * SPDX-License-Identifier: LGPL-2.1-or-later
  * Copyright (C) 2016 Manuel Quezada de Luna
  *
  * This file is part of the deal.II code gallery.
  *
  * -----------------------------------------------------------------------------
  */
  #include <deal.II/base/function.h>
  #include <deal.II/lac/affine_constraints.h>
  #include <deal.II/lac/vector.h>
  #include <deal.II/lac/petsc_vector.h>
  #include <deal.II/dofs/dof_handler.h>
  #include <deal.II/dofs/dof_tools.h>
  #include <deal.II/fe/fe_values.h>
  #include <deal.II/fe/fe_q.h>
  #include <deal.II/numerics/vector_tools.h>
  #include <deal.II/numerics/data_out.h>
  #include <deal.II/numerics/error_estimator.h>
  #include <deal.II/base/utilities.h>
  #include <deal.II/base/index_set.h>
  #include <deal.II/distributed/tria.h>
  #include <deal.II/distributed/grid_refinement.h>
  #include <deal.II/lac/petsc_vector.h>
  #include <deal.II/base/timer.h>
  #include <deal.II/grid/grid_tools.h>
  #include <deal.II/fe/mapping_q.h>
  #include <deal.II/fe/fe_system.h>
  using namespace dealii;

/////////////////////// FOR TRANSPORT PROBLEM /////////////////////// TIME_INTEGRATION

  #define FORWARD_EULER 0

PROBLEM

OTHER FLAGS

  #include "utilities_test_LS.cc"
  #include "LevelSetSolver.cc"

/////////////////////////////////////////////////// /////////////////// MAIN CLASS //////////////////// ///////////////////////////////////////////////////

  template <int dim>
  {
  public:
  TestLevelSet (const unsigned int degree_LS,
  const unsigned int degree_U);
  void run ();
  private:

BOUNDARY

  void get_boundary_values_phi(std::vector<unsigned int> &boundary_values_id_phi,
  std::vector<double> &boundary_values_phi);

VELOCITY

SETUP AND INIT CONDITIONS

  void setup();

POST PROCESSING

SOLUTION VECTORS

BOUNDARY VECTORS

  std::vector<unsigned int> boundary_values_id_phi;
  std::vector<double> boundary_values_phi;

GENERAL

FOR TRANSPORT

  double cK; //compression coeff
  double cE; //entropy-visc coeff
  std::string ALGORITHM;
  unsigned int PROBLEM;

FOR RECONSTRUCTION OF MATERIAL FIELDS

  double eps, rho_air, rho_fluid;

MASS MATRIX

  std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_MC;
  };
  template <int dim>
  const unsigned int degree_U)
  :
  mpi_communicator (MPI_COMM_WORLD),
  triangulation (mpi_communicator,
  typename Triangulation<dim>::MeshSmoothing
  (Triangulation<dim>::smoothing_on_refinement |
  Triangulation<dim>::smoothing_on_coarsening)),
  pcout (std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator)== 0))
  {}
  {
  }
Definition fe_q.h:554

VELOCITY //////////

velocity in x

velocity in y

////////// BOUNDARY //////////

  template <int dim>
  {
  const QGauss<dim-1> face_quadrature_formula(1); // center of the face
  const unsigned int n_face_q_points = face_quadrature_formula.size();
  std::vector<double> u_value (n_face_q_points);
  std::vector<double> v_value (n_face_q_points);
  std::vector<double> w_value (n_face_q_points);
  cell_U = dof_handler_U.begin_active(),
  for (; cell_U!=endc_U; ++cell_U)
  if (cell_U->is_locally_owned())
  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
  if (cell_U->face(face)->at_boundary())
  {
  fe_face_values.reinit(cell_U,face);
  fe_face_values.get_function_values(locally_relevant_solution_u,u_value);
  fe_face_values.get_function_values(locally_relevant_solution_v,v_value);
  if (dim==3)
  fe_face_values.get_function_values(locally_relevant_solution_w,w_value);
  u[0]=u_value[0];
  u[1]=v_value[0];
  if (dim==3)
  u[2]=w_value[0];
  if (fe_face_values.normal_vector(0)*u < -1e-14)
  cell_U->face(face)->set_boundary_id(10);
  }
  }
  template <int dim>
  std::vector<double> &boundary_values_phi)
  {
  std::map<unsigned int, double> map_boundary_values_phi;
  unsigned int boundary_id=0;
  boundary_id=10; // inlet
  boundary_id,BoundaryPhi<dim>(),
  std::map<unsigned int,double>::const_iterator boundary_value_phi = map_boundary_values_phi.begin();
  {
  }
  }

/////////////////////////////// SETUP AND INITIAL CONDITIONS //////////////////////////////

  template <int dim>
  {

setup system LS

setup system U

setup system U for disp field

init vectors for phi

init vectors for u

init vectors for v

init vectors for w

MASS MATRIX

  dof_handler_LS.n_locally_owned_dofs_per_processor(),
  mpi_communicator,
  matrix_MC.reinit (mpi_communicator,
  dof_handler_LS.n_locally_owned_dofs_per_processor(),
  dof_handler_LS.n_locally_owned_dofs_per_processor(),
  matrix_MC_tnm1.reinit (mpi_communicator,
  dof_handler_LS.n_locally_owned_dofs_per_processor(),
  dof_handler_LS.n_locally_owned_dofs_per_processor(),
  }
  template <int dim>
  {
  time=0;

Initial conditions init condition for phi

Functions::ZeroFunction<dim>(),

init condition for u=0

init condition for v

///////////////// POST PROCESSING /////////////////

error for phi

  solution,
  solution,
  double u_L2_error = difference_per_cell.l2_norm();
  pcout << "L1 error: " << u_L1_error << std::endl;
  pcout << "L2 error: " << u_L2_error << std::endl;
  }
  template<int dim>
  {
  }
  template <int dim>
  {
  DataOut<dim> data_out;
  data_out.attach_dof_handler(dof_handler_LS);
  data_out.add_data_vector (locally_relevant_solution_phi, "phi");
  data_out.build_patches();
  const std::string filename = ("solution-" +
  "." +
  (triangulation.locally_owned_subdomain(), 4));
  std::ofstream output ((filename + ".vtu").c_str());
  data_out.write_vtu (output);
  if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
  {
  std::vector<std::string> filenames;
  for (unsigned int i=0;
  i<Utilities::MPI::n_mpi_processes(mpi_communicator);
  ++i)
  filenames.push_back ("solution-" +
  "." +
  ".vtu");
  std::ofstream master_output ((filename + ".pvtu").c_str());
  data_out.write_pvtu_record (master_output, filenames);
  }
  }
  template <int dim>
  {
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition tensor.h:3021
void integrate_difference(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const ReadVector< Number > &fe_function, const Function< spacedim, Number > &exact_solution, OutVector &difference, const Quadrature< dim > &q, const NormType &norm, const Function< spacedim, double > *weight=nullptr, const double exponent=2.)

//////////////////// GENERAL PARAMETERS ////////////////////

  cfl=0.1;
  verbose = false;
  get_output = true;
  final_time = 1.0;
Definition timer.h:117

PROBLEM=DIAGONAL_ADVECTION;

  double umax = 0;
  else
constexpr double PI
Definition numbers.h:240

////////////////////////////////// PARAMETERS FOR TRANSPORT PROBLEM //////////////////////////////////

  cK = 1.0; // compression constant
  cE = 1.0; // entropy viscosity constant
  sharpness_integer=1; //this will be multiplied by min_h

TRANSPORT_TIME_INTEGRATION=FORWARD_EULER;

ALGORITHM = "MPP_u1";

  ALGORITHM = "NMPP_uH";

ALGORITHM = "MPP_uH";

////////// GEOMETRY //////////

void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)

GridGenerator::hyper_rectangle(triangulation, Point<dim>(0.0,0.0), Point<dim>(1.0,1.0), true);

  triangulation.refine_global (n_refinement);

/////// SETUP ///////

  setup();

for Reconstruction of MATERIAL FIELDS

  eps=1*min_h; //For reconstruction of density in Navier Stokes
  sharpness=sharpness_integer*min_h; //adjust value of sharpness (for init cond of phi)
  rho_fluid = 1000;
  rho_air = 1;

GET TIME STEP

////////////////// TRANSPORT SOLVER //////////////////

/////////////////// INITIAL CONDITION ///////////////////

///////////////////////////// BOUNDARY CONDITIONS FOR PHI /////////////////////////////

OUTPUT DATA REGARDING TIME STEPPING AND MESH

  int dofs_LS = dof_handler_LS.n_dofs();
  pcout << "Cfl: " << cfl << std::endl;
  pcout << " Number of active cells: "
  << triangulation.n_global_active_cells() << std::endl
  << " Number of degrees of freedom: " << std::endl
  << " LS: " << dofs_LS << std::endl;

TIME STEPPING

  time=0;
  while (time<final_time)
  {
  if (time+time_step > final_time)
  {
  pcout << "FINAL TIME STEP... " << std::endl;
  time_step = final_time-time;
  }
  pcout << "Time step " << timestep_number
  << "\twith dt=" << time_step
  << "\tat tn=" << time << std::endl;

////////////// GET VELOCITY // (NS or interpolate from a function) at current time tn //////////////

SET VELOCITY TO LEVEL SET SOLVER

//////////////////////// GET LEVEL SET SOLUTION // (at tnp1) ////////////////////////

  level_set.nth_time_step();

///////////// UPDATE TIME /////////////

  time+=time_step; // time tnp1

//////// OUTPUT ////////

  {
  level_set.get_unp1(locally_relevant_solution_phi);
  }
  }
  pcout << "FINAL TIME T=" << time << std::endl;
  }
  int main(int argc, char *argv[])
  {
  try
  {
  using namespace dealii;
  {
  unsigned int degree = 1;
  TestLevelSet<2> multiphase(degree, degree);
  multiphase.run();
  }
  }
  catch (std::exception &exc)
  {
  std::cerr << std::endl << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Exception on processing: " << std::endl
  << exc.what() << std::endl
  << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
  return 1;
  }
  catch (...)
  {
  std::cerr << std::endl << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Unknown exception!" << std::endl
  << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
  return 1;
  }
  return 0;
  }

Annotated version of TestNavierStokes.cc

  /* -----------------------------------------------------------------------------
  *
  * SPDX-License-Identifier: LGPL-2.1-or-later
  * Copyright (C) 2016 Manuel Quezada de Luna
  *
  * This file is part of the deal.II code gallery.
  *
  * -----------------------------------------------------------------------------
  */
  #include <deal.II/base/function.h>
  #include <deal.II/lac/affine_constraints.h>
  #include <deal.II/lac/vector.h>
  #include <deal.II/lac/petsc_vector.h>
  #include <deal.II/dofs/dof_handler.h>
  #include <deal.II/dofs/dof_tools.h>
  #include <deal.II/fe/fe_values.h>
  #include <deal.II/fe/fe_q.h>
  #include <deal.II/numerics/vector_tools.h>
  #include <deal.II/numerics/data_out.h>
  #include <deal.II/numerics/error_estimator.h>
  #include <deal.II/base/utilities.h>
  #include <deal.II/base/index_set.h>
  #include <deal.II/distributed/tria.h>
  #include <deal.II/distributed/grid_refinement.h>
  #include <deal.II/lac/petsc_vector.h>
  #include <deal.II/base/timer.h>
  #include <deal.II/grid/grid_tools.h>
  #include <deal.II/fe/mapping_q.h>
  #include <deal.II/base/function.h>
  using namespace dealii;
  #include "utilities_test_NS.cc"
  #include "NavierStokesSolver.cc"

/////////////////////////////////////////////////// /////////////////// MAIN CLASS //////////////////// ///////////////////////////////////////////////////

  template <int dim>
  {
  public:
  TestNavierStokes (const unsigned int degree_LS,
  const unsigned int degree_U);
  void run ();
  private:
  void get_boundary_values_U(double t);
  void process_solution(const unsigned int cycle);
  void setup();
  std::vector<unsigned int> boundary_values_id_u;
  std::vector<unsigned int> boundary_values_id_v;
  std::vector<unsigned int> boundary_values_id_w;
  std::vector<double> boundary_values_u;
  std::vector<double> boundary_values_v;
  std::vector<double> boundary_values_w;
  double rho_fluid;
  double nu_fluid;
  double rho_air;
  double nu_air;
  MPI_Comm mpi_communicator;

TimerOutput timer;

  double time;
  double time_step;
  double final_time;
  unsigned int timestep_number;
  double cfl;
  double min_h;
  unsigned int n_cycles;
  unsigned int n_refinement;
  unsigned int output_number;
  double output_time;
  double h;
  double umax;
  bool verbose;
  double nu;
  };
  template <int dim>
  const unsigned int degree_U)
  :
  mpi_communicator (MPI_COMM_WORLD),
  triangulation (mpi_communicator,
  typename Triangulation<dim>::MeshSmoothing
  (Triangulation<dim>::smoothing_on_refinement |
  Triangulation<dim>::smoothing_on_coarsening)),
  fe_P (degree_U-1), //TODO: change this to be degree_Q-1

timer(std::cout, TimerOutput::summary, TimerOutput::wall_times),

///////////////////////////////////// /////////////// SETUP /////////////// /////////////////////////////////////

  template <int dim>
  {

setup system LS

setup system U

setup system P

init vectors for rho

init vectors for u

init vectors for v

init vectors for w

init vectors for p

Initial conditions init condition for rho

init condition for u

init condition for v

init condition for w

init condition for p

fix the constant in the pressure

  0);
  if (dim==2)
  completely_distributed_solution_p.add(-mean_value+std::sin(1)*(std::cos(time)-cos(1+time)));
  else
  completely_distributed_solution_p.add(-mean_value+8*std::pow(std::sin(0.5),3)*std::sin(1.5+time));
  }
  template <int dim>
  {
  DataOut<dim> data_out;
  data_out.attach_dof_handler (dof_handler_U);
  data_out.add_data_vector (locally_relevant_solution_u, "u");
  data_out.add_data_vector (locally_relevant_solution_v, "v");
  if (dim==3) data_out.add_data_vector (locally_relevant_solution_w, "w");
  for (unsigned int i=0; i<subdomain.size(); ++i)
  subdomain(i) = triangulation.locally_owned_subdomain();
  data_out.add_data_vector (subdomain, "subdomain");
  data_out.build_patches ();
  const std::string filename = ("solution-" +
  "." +
  (triangulation.locally_owned_subdomain(), 4));
  std::ofstream output ((filename + ".vtu").c_str());
  data_out.write_vtu (output);
  if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
  {
  std::vector<std::string> filenames;
  for (unsigned int i=0;
  i<Utilities::MPI::n_mpi_processes(mpi_communicator);
  ++i)
  filenames.push_back ("solution-" +
  "." +
  ".vtu");
  std::ofstream master_output ((filename + ".pvtu").c_str());
  data_out.write_pvtu_record (master_output, filenames);
  }
  }
  template <int dim>
  void TestNavierStokes<dim>::process_solution(const unsigned int cycle)
  {
Number compute_mean_value(const hp::MappingCollection< dim, spacedim > &mapping_collection, const DoFHandler< dim, spacedim > &dof, const hp::QCollection< dim > &q_collection, const ReadVector< Number > &v, const unsigned int component)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)

error for u

error for v

error for w

error for p

  double p_L2_error = difference_per_cell.l2_norm();
  mpi_communicator));
  double p_H1_error = difference_per_cell.l2_norm();
  mpi_communicator));
  const unsigned int n_active_cells=triangulation.n_active_cells();
  const unsigned int n_dofs_U=dof_handler_U.n_dofs();
  const unsigned int n_dofs_P=dof_handler_P.n_dofs();
  convergence_table.add_value("cycle", cycle);
  convergence_table.add_value("cells", n_active_cells);
  convergence_table.add_value("dofs_U", n_dofs_U);
  convergence_table.add_value("dofs_P", n_dofs_P);
  convergence_table.add_value("dt", time_step);
  convergence_table.add_value("u L2", u_L2_error);
  convergence_table.add_value("u H1", u_H1_error);
  convergence_table.add_value("v L2", v_L2_error);
  convergence_table.add_value("v H1", v_H1_error);
  if (dim==3)
  {
  convergence_table.add_value("w L2", w_L2_error);
  convergence_table.add_value("w H1", w_H1_error);
  }
  convergence_table.add_value("p L2", p_L2_error);
  convergence_table.add_value("p H1", p_H1_error);
  }
  template <int dim>
  {
  std::map<unsigned int, double> map_boundary_values_u;
  std::map<unsigned int, double> map_boundary_values_v;
  std::map<unsigned int,double>::const_iterator boundary_value_u =map_boundary_values_u.begin();
  std::map<unsigned int,double>::const_iterator boundary_value_v =map_boundary_values_v.begin();
  if (dim==3)
  {
  std::map<unsigned int, double> map_boundary_values_w;
  std::map<unsigned int,double>::const_iterator boundary_value_w =map_boundary_values_w.begin();
  {
  }
  }
  {
  }
  {
  }
  }
  template <int dim>
  {
  if (Utilities::MPI::this_mpi_process(mpi_communicator)== 0)
  {
  std::cout << "***** CONVERGENCE TEST FOR NS *****" << std::endl;
  std::cout << "DEGREE LS: " << degree_LS << std::endl;
  std::cout << "DEGREE U: " << degree_U << std::endl;
  }

PARAMETERS FOR THE NAVIER STOKES PROBLEM

  final_time = 1.0;
  n_cycles=6;
  bool get_output = false;
  bool get_error = true;
  verbose = true;
  for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
  {
  if (cycle == 0)
  {
  triangulation.refine_global (n_refinement);
  setup();
  }
  else
  {
  triangulation.refine_global(1);
  setup();
  }

if (cycle==0)

set INITIAL CONDITION within TRANSPORT PROBLEM

  if (dim==2)
  else //dim=3
  pcout << "Cycle " << cycle << ':' << std::endl;
  pcout << " Cycle " << cycle
  << " Number of active cells: "
  << triangulation.n_global_active_cells() << std::endl
  << " Number of degrees of freedom (velocity): "
  << dof_handler_U.n_dofs() << std::endl
  << std::endl;

TIME STEPPING

  time=0;
  while (time<final_time)
  {

/////////////// GET TIME_STEP ///////////////

  if (time+time_step > final_time-1E-10)
  {
  pcout << "FINAL TIME STEP..." << std::endl;
  time_step=final_time-time;
  }
  pcout << "Time step " << timestep_number
  << "\twith dt=" << time_step
  << "\tat tn=" << time
  << std::endl;

///////////// FORCE TERMS /////////////

  force_function.set_time(time+time_step);

///////////////////////////// DENSITY AND VISCOSITY FIELD /////////////////////////////

  rho_function.set_time(time+time_step);
  nu_function.set_time(time+time_step);

///////////////////// BOUNDARY CONDITIONS /////////////////////

////////////// GET SOLUTION //////////////

////////////// FIX PRESSURE //////////////

///////////// UPDATE TIME /////////////

  time+=time_step;

//////// OUTPUT ////////

  if (get_output && time-(output_number)*output_time>=1E-10)
  }
  pcout << "FINAL TIME: " << time << std::endl;
  {
  convergence_table.set_precision("u L2", 2);
  convergence_table.set_precision("u H1", 2);
  convergence_table.set_scientific("u L2",true);
  convergence_table.set_scientific("u H1",true);
  convergence_table.set_precision("v L2", 2);
  convergence_table.set_precision("v H1", 2);
  convergence_table.set_scientific("v L2",true);
  convergence_table.set_scientific("v H1",true);
  if (dim==3)
  {
  convergence_table.set_precision("w L2", 2);
  convergence_table.set_precision("w H1", 2);
  convergence_table.set_scientific("w L2",true);
  convergence_table.set_scientific("w H1",true);
  }
  convergence_table.set_precision("p L2", 2);
  convergence_table.set_precision("p H1", 2);
  convergence_table.set_scientific("p L2",true);
  convergence_table.set_scientific("p H1",true);
  convergence_table.set_tex_format("cells","r");
  convergence_table.set_tex_format("dofs_U","r");
  convergence_table.set_tex_format("dofs_P","r");
  convergence_table.set_tex_format("dt","r");
  if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
  {
  std::cout << std::endl;
  convergence_table.write_text(std::cout);
  }
  }
  }
  }
  int main(int argc, char *argv[])
  {
  try
  {
  using namespace dealii;
  {
  unsigned int degree_LS = 1;
  unsigned int degree_U = 2;
  }
  }
  catch (std::exception &exc)
  {
  std::cerr << std::endl << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Exception on processing: " << std::endl
  << exc.what() << std::endl
  << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
  return 1;
  }
  catch (...)
  {
  std::cerr << std::endl << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Unknown exception!" << std::endl
  << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
  return 1;
  }
  return 0;
  }

Annotated version of clean.sh

rm -rf CMakeFiles CMakeCache.txt Makefile cmake_install.cmake *~
rm -f MultiPhase TestLevelSet TestNavierStokes
rm -f sol*
rm -f *#*
rm -f *.visit

Annotated version of utilities.cc

  /* -----------------------------------------------------------------------------
  *
  * SPDX-License-Identifier: LGPL-2.1-or-later
  * Copyright (C) 2016 Manuel Quezada de Luna
  *
  * This file is part of the deal.II code gallery.
  *
  * -----------------------------------------------------------------------------
  */

///////////////////////////////////////////////// ////////////////// INITIAL PHI ////////////////// /////////////////////////////////////////////////

  template <int dim>
  class InitialPhi : public Function <dim>
  {
  public:
  InitialPhi (unsigned int PROBLEM, double sharpness=0.005) : Function<dim>(),
  virtual double value (const Point<dim> &p, const unsigned int component=0) const override;
  double sharpness;
  unsigned int PROBLEM;
  };
  template <int dim>
  const unsigned int) const
  {
  double x = p[0];
  double y = p[1];
  double pi=numbers::PI;
  return 0.5*(-std::tanh((y-0.3)/sharpness)*std::tanh((y-0.35)/sharpness)+1)
  *(-std::tanh((x-0.02)/sharpness)+1)-1;
  return 0.5*(-std::tanh((x-0.35)/sharpness)*std::tanh((x-0.65)/sharpness)+1)
  *(1-std::tanh((y-0.35)/sharpness))-1;
  {
  double x0=0.15;
  double y0=0.75;
  double r0=0.1;
  double r = std::sqrt(std::pow(x-x0,2)+std::pow(y-y0,2));
  return 1-(std::tanh((r-r0)/sharpness)+std::tanh((y-0.3)/sharpness));
  }
  {
  double wave = 0.1*std::sin(pi*x)+0.25;
  }
  else
  {
  std::cout << "Error in type of PROBLEM" << std::endl;
  }
  }
void abort(const ExceptionBase &exc) noexcept
inline ::VectorizedArray< Number, width > tanh(const ::VectorizedArray< Number, width > &x)

/////////////////////////////////////////////////// ////////////////// FORCE TERMS ///// ////////////// ///////////////////////////////////////////////////

  template <int dim>
  {
  public:
  ForceTerms (const std::vector<double> values) : Functions::ConstantFunction<dim>(values) {}
  };
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
Definition tensor.h:852

///////////////////////////////////////////////// ////////////////// BOUNDARY PHI ///////////////// /////////////////////////////////////////////////

  template <int dim>
  {
  public:
  BoundaryPhi (const double value, const unsigned int n_components=1) : Functions::ConstantFunction<dim>(value,n_components) {}
  };

////////////////////////////////////////////////////// ////////////////// BOUNDARY VELOCITY ///////////////// //////////////////////////////////////////////////////

  template <int dim>
  class BoundaryU : public Function <dim>
  {
  public:
  BoundaryU (unsigned int PROBLEM, double t=0) : Function<dim>(), PROBLEM(PROBLEM) {this->set_time(t);}
  virtual double value (const Point<dim> &p, const unsigned int component=0) const override;
  unsigned PROBLEM;
  };
  template <int dim>
  double BoundaryU<dim>::value (const Point<dim> &p, const unsigned int) const
  {

////////////////// FILLING THE TANK ////////////////// boundary for filling the tank (inlet)

  double x = p[0];
  double y = p[1];
  {
  if (x==0 && y>=0.3 && y<=0.35)
  return 0.25;
  else
  return 0.0;
  }
  else
  {
  std::cout << "Error in PROBLEM definition" << std::endl;
  abort();
  }
  }
  template <int dim>
  class BoundaryV : public Function <dim>
  {
  public:
  BoundaryV (unsigned int PROBLEM, double t=0) : Function<dim>(), PROBLEM(PROBLEM) {this->set_time(t);}
  virtual double value (const Point<dim> &p, const unsigned int component=0) const override;
  unsigned int PROBLEM;
  };
  template <int dim>
  double BoundaryV<dim>::value (const Point<dim> &p, const unsigned int) const
  {

boundary for filling the tank (outlet)

  double x = p[0];
  double y = p[1];
  double return_value = 0;
  {
  if (y==0.4 && x>=0.3 && x<=0.35)
  return_value = 0.25;
  }
  return return_value;
  }

/////////////////////////////////////////////////// ///////////////// POST-PROCESSING ///////////////// ///////////////////////////////////////////////////

  template <int dim>
  {
  public:
  Postprocessor(double eps, double rho_air, double rho_fluid)
  :
  {
  this->eps=eps;
  this->rho_air=rho_air;
  this->rho_fluid=rho_fluid;
  }
  virtual
  void
  evaluate_scalar_field (const DataPostprocessorInputs::Scalar<dim> &input_data,
  std::vector<Vector<double> > &computed_quantities) const override;
  double eps;
  double rho_air;
  double rho_fluid;
  };
  template <int dim>
  void
  evaluate_scalar_field (const DataPostprocessorInputs::Scalar<dim> &input_data,
  std::vector<Vector<double> > &computed_quantities) const
  {
  const unsigned int n_quadrature_points = input_data.solution_values.size();
  for (unsigned int q=0; q<n_quadrature_points; ++q)
  {
  double H;
  double rho_value;
  double phi_value=input_data.solution_values[q];
  if (phi_value > eps)
  H=1;
  else if (phi_value < -eps)
  H=-1;
  else
  rho_value = rho_fluid*(1+H)/2. + rho_air*(1-H)/2.;
  }
  }

Annotated version of utilities_test_LS.cc

  /* -----------------------------------------------------------------------------
  *
  * SPDX-License-Identifier: LGPL-2.1-or-later
  * Copyright (C) 2016 Manuel Quezada de Luna
  *
  * This file is part of the deal.II code gallery.
  *
  * -----------------------------------------------------------------------------
  */

/////////////////////////////////////////////////// ////////////////// INITIAL PHI ////////////////// ///////////////////////////////////////////////////

  template <int dim>
  class InitialPhi : public Function <dim>
  {
  public:
  InitialPhi (unsigned int PROBLEM, double sharpness=0.005) : Function<dim>(),
  virtual double value (const Point<dim> &p, const unsigned int component=0) const;
  double sharpness;
  unsigned int PROBLEM;
  };
  template <int dim>
  const unsigned int) const
  {
  double x = p[0];
  double y = p[1];
  double return_value = -1.;
  {
  double x0=0.5;
  double y0=0.75;
  double r0=0.15;
  double r = std::sqrt(std::pow(x-x0,2)+std::pow(y-y0,2));
  return_value = -std::tanh((r-r0)/sharpness);
  }
  else // (PROBLEM==DIAGONAL_ADVECTION)
  {
  double x0=0.25;
  double y0=0.25;
  double r0=0.15;
  double r=0;
  if (dim==2)
  else
  {
  double z0=0.25;
  double z=p[2];
  }
  return_value = -std::tanh((r-r0)/sharpness);
  }
  return return_value;
  }

///////////////////////////////////////////////// ////////////////// BOUNDARY PHI ///////////////// /////////////////////////////////////////////////

  template <int dim>
  class BoundaryPhi : public Function <dim>
  {
  public:
  BoundaryPhi (double t=0)
  :
  Function<dim>()
  {this->set_time(t);}
  virtual double value (const Point<dim> &p, const unsigned int component=0) const;
  };
  template <int dim>
  double BoundaryPhi<dim>::value (const Point<dim> &, const unsigned int) const
  {
  return -1.0;
  }

/////////////////////////////////////////////////// ////////////////// EXACT VELOCITY ///////////////// ///////////////////////////////////////////////////

  template <int dim>
  class ExactU : public Function <dim>
  {
  public:
  ExactU (unsigned int PROBLEM, double time=0) : Function<dim>(), PROBLEM(PROBLEM), time(time) {}
  virtual double value (const Point<dim> &p, const unsigned int component=0) const;
  void set_time(double time) {this->time=time;};
  unsigned PROBLEM;
  double time;
  };
  template <int dim>
  double ExactU<dim>::value (const Point<dim> &p, const unsigned int) const
  {
  return -2*numbers::PI*(p[1]-0.5);
  else // (PROBLEM==DIAGONAL_ADVECTION)
  return 1.0;
  }
  template <int dim>
  class ExactV : public Function <dim>
  {
  public:
  ExactV (unsigned int PROBLEM, double time=0) : Function<dim>(), PROBLEM(PROBLEM), time(time) {}
  virtual double value (const Point<dim> &p, const unsigned int component=0) const;
  void set_time(double time) {this->time=time;};
  unsigned int PROBLEM;
  double time;
  };
  template <int dim>
  double ExactV<dim>::value (const Point<dim> &p, const unsigned int) const
  {
  return 2*numbers::PI*(p[0]-0.5);
  else // (PROBLEM==DIAGONAL_ADVECTION)
  return 1.0;
  }
  template <int dim>
  class ExactW : public Function <dim>
  {
  public:
  ExactW (unsigned int PROBLEM, double time=0) : Function<dim>(), PROBLEM(PROBLEM), time(time) {}
  virtual double value (const Point<dim> &p, const unsigned int component=0) const;
  void set_time(double time) {this->time=time;};
  unsigned int PROBLEM;
  double time;
  };
  template <int dim>
  double ExactW<dim>::value (const Point<dim> &, const unsigned int) const
  {

PROBLEM = 3D_DIAGONAL_ADVECTION

  return 1.0;
  }

Annotated version of utilities_test_NS.cc

  /* -----------------------------------------------------------------------------
  *
  * SPDX-License-Identifier: LGPL-2.1-or-later
  * Copyright (C) 2016 Manuel Quezada de Luna
  *
  * This file is part of the deal.II code gallery.
  *
  * -----------------------------------------------------------------------------
  */

/////////////////////////////////////////////////// ////////// EXACT SOLUTION RHO TO TEST NS ////////// ///////////////////////////////////////////////////

  template <int dim>
  class RhoFunction : public Function <dim>
  {
  public:
  RhoFunction (double t=0) : Function<dim>() {this->set_time(t);}
  virtual double value (const Point<dim> &p, const unsigned int component=0) const;
  };
  template <int dim>
  const unsigned int) const
  {
  double t = this->get_time();
  double return_value = 0;
  if (dim==2)
  return_value = std::pow(std::sin(p[0]+p[1]+t),2)+1;
  else //dim=3
  return_value = std::pow(std::sin(p[0]+p[1]+p[2]+t),2)+1;
  return return_value;
  }
  template <int dim>
  class NuFunction : public Function <dim>
  {
  public:
  NuFunction (double t=0) : Function<dim>() {this->set_time(t);}
  virtual double value (const Point<dim> &p, const unsigned int component=0) const;
  };
  template <int dim>
  double NuFunction<dim>::value (const Point<dim> &, const unsigned int) const
  {
  return 1.;
  }
std::string get_time()

////////////////////////////////////////////////////////////// ///////////////// EXACT SOLUTION U to TEST NS //////////////// //////////////////////////////////////////////////////////////

  template <int dim>
  class ExactSolution_and_BC_U : public Function <dim>
  {
  public:
  ExactSolution_and_BC_U (double t=0, int field=0)
  :
  Function<dim>(),
  {
  this->set_time(t);
  }
  virtual double value (const Point<dim> &p, const unsigned int component=1) const;
  virtual Tensor<1,dim> gradient (const Point<dim> &p, const unsigned int component=1) const;
  virtual void set_field(int field) {this->field=field;}
  int field;
  unsigned int type_simulation;
  };
  template <int dim>
  const unsigned int) const
  {
  double t = this->get_time();
  double return_value = 0;
  double Pi = numbers::PI;
  double x = p[0];
  double y = p[1];
  double z = 0;
  if (dim == 2)
  if (field == 0)
  return_value = std::sin(x)*std::sin(y+t);
  else
  return_value = std::cos(x)*std::cos(y+t);
  else //dim=3
  {
  z = p[2];
  if (field == 0)
  return_value = std::cos(t)*std::cos(Pi*y)*std::cos(Pi*z)*std::sin(Pi*x);
  else if (field == 1)
  return_value = std::cos(t)*std::cos(Pi*x)*std::cos(Pi*z)*std::sin(Pi*y);
  else
  return_value = -2*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::sin(Pi*z);
  }
  return return_value;
  }
  template <int dim>
  const unsigned int) const
  {

THIS IS USED JUST FOR TESTING NS

  Tensor<1,dim> return_value;
  double t = this->get_time();
  double Pi = numbers::PI;
  double x = p[0];
  double y = p[1];
  double z = 0;
  if (dim == 2)
  if (field == 0)
  {
  return_value[0] = std::cos(x)*std::sin(y+t);
  return_value[1] = std::sin(x)*std::cos(y+t);
  }
  else
  {
  return_value[0] = -std::sin(x)*std::cos(y+t);
  return_value[1] = -std::cos(x)*std::sin(y+t);
  }
  else //dim=3
  {
  z=p[2];
  if (field == 0)
  {
  return_value[0] = Pi*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::cos(Pi*z);
  return_value[1] = -(Pi*std::cos(t)*std::cos(Pi*z)*std::sin(Pi*x)*std::sin(Pi*y));
  return_value[2] = -(Pi*std::cos(t)*std::cos(Pi*y)*std::sin(Pi*x)*std::sin(Pi*z));
  }
  else if (field == 1)
  {
  return_value[0] = -(Pi*std::cos(t)*std::cos(Pi*z)*std::sin(Pi*x)*std::sin(Pi*y));
  return_value[1] = Pi*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::cos(Pi*z);
  return_value[2] = -(Pi*std::cos(t)*std::cos(Pi*x)*std::sin(Pi*y)*std::sin(Pi*z));
  }
  else
  {
  return_value[0] = 2*Pi*std::cos(t)*std::cos(Pi*y)*std::sin(Pi*x)*std::sin(Pi*z);
  return_value[1] = 2*Pi*std::cos(t)*std::cos(Pi*x)*std::sin(Pi*y)*std::sin(Pi*z);
  return_value[2] = -2*Pi*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::cos(Pi*z);
  }
  }
  return return_value;
  }

/////////////////////////////////////////////////// ///////// EXACT SOLUTION FOR p TO TEST NS ///////// ///////////////////////////////////////////////////

  template <int dim>
  class ExactSolution_p : public Function <dim>
  {
  public:
  ExactSolution_p (double t=0) : Function<dim>() {this->set_time(t);}
  virtual double value (const Point<dim> &p, const unsigned int component=0) const;
  virtual Tensor<1,dim> gradient (const Point<dim> &p, const unsigned int component = 0) const;
  };
  template <int dim>
  double ExactSolution_p<dim>::value (const Point<dim> &p, const unsigned int) const
  {
  double t = this->get_time();
  double return_value = 0;
  if (dim == 2)
  return_value = std::cos(p[0])*std::sin(p[1]+t);
  else //dim=3
  return_value = std::sin(p[0]+p[1]+p[2]+t);
  return return_value;
  }
  template <int dim>
  Tensor<1,dim> ExactSolution_p<dim>::gradient (const Point<dim> &p, const unsigned int) const
  {
  Tensor<1,dim> return_value;
  double t = this->get_time();
  if (dim == 2)
  {
  return_value[0] = -std::sin(p[0])*std::sin(p[1]+t);
  return_value[1] = std::cos(p[0])*std::cos(p[1]+t);
  }
  else //dim=3
  {
  return_value[0] = std::cos(t+p[0]+p[1]+p[2]);
  return_value[1] = std::cos(t+p[0]+p[1]+p[2]);
  return_value[2] = std::cos(t+p[0]+p[1]+p[2]);
  }
  return return_value;
  }

////////////////////////////////////////////////////////////// ////////////////// FORCE TERMS to TEST NS //////////////////// //////////////////////////////////////////////////////////////

  template <int dim>
  class ForceTerms : public Function <dim>
  {
  public:
  ForceTerms (double t=0)
  :
  Function<dim>()
  {
  this->set_time(t);
  nu = 1.;
  }
  virtual void vector_value (const Point<dim> &p, Vector<double> &values) const;
  double nu;
  };
  template <int dim>
  {
  double x = p[0];
  double y = p[1];
  double z = 0;
  double t = this->get_time();
  double Pi = numbers::PI;
  if (dim == 2)
  {

force in x

  values[0] = std::cos(t+y)*std::sin(x)*(1+std::pow(std::sin(t+x+y),2)) // time derivative
  +2*nu*std::sin(x)*std::sin(t+y) // viscosity
  +std::cos(x)*std::sin(x)*(1+std::pow(std::sin(t+x+y),2)) // non-linearity
  -std::sin(x)*std::sin(y+t); // pressure

force in y

  values[1] = -(std::cos(x)*std::sin(t+y)*(1+std::pow(std::sin(t+x+y),2))) // time derivative
  +2*nu*std::cos(x)*std::cos(t+y) // viscosity
  -(std::sin(2*(t+y))*(1+std::pow(std::sin(t+x+y),2)))/2. // non-linearity
  +std::cos(x)*std::cos(y+t); // pressure
  }
  else //3D
  {
  z = p[2];

force in x

  values[0]=
  -(std::cos(Pi*y)*std::cos(Pi*z)*std::sin(t)*std::sin(Pi*x)*(1+std::pow(std::sin(t+x+y+z),2))) //time der.
  +3*std::pow(Pi,2)*std::cos(t)*std::cos(Pi*y)*std::cos(Pi*z)*std::sin(Pi*x) //viscosity
  -(Pi*std::pow(std::cos(t),2)*(-3+std::cos(2*(t+x+y+z)))*std::sin(2*Pi*x)*(std::cos(2*Pi*y)+std::pow(std::sin(Pi*z),2)))/4. //NL
  +std::cos(t+x+y+z); // pressure
  values[1]=
  -(std::cos(Pi*x)*std::cos(Pi*z)*std::sin(t)*std::sin(Pi*y)*(1+std::pow(std::sin(t+x+y+z),2))) //time der
  +3*std::pow(Pi,2)*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*z)*std::sin(Pi*y) //viscosity
  -(Pi*std::pow(std::cos(t),2)*(-3+std::cos(2*(t+x+y+z)))*std::sin(2*Pi*y)*(std::cos(2*Pi*x)+std::pow(std::sin(Pi*z),2)))/4. //NL
  +std::cos(t+x+y+z); // pressure
  values[2]=
  2*std::cos(Pi*x)*std::cos(Pi*y)*std::sin(t)*std::sin(Pi*z)*(1+std::pow(std::sin(t+x+y+z),2)) //time der
  -6*std::pow(Pi,2)*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::sin(Pi*z) //viscosity
  -(Pi*std::pow(std::cos(t),2)*(2+std::cos(2*Pi*x)+std::cos(2*Pi*y))*(-3+std::cos(2*(t+x+y+z)))*std::sin(2*Pi*z))/4. //NL
  +std::cos(t+x+y+z); // pressure
  }
  }