![]() |
deal.II version GIT relicensing-3356-g4636aadadd 2025-05-22 18:40:00+00:00
|
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):
![]() | ![]() | ![]() | ![]() |
![]() | ![]() | ![]() | ![]() |
![]() | ![]() | ![]() | ![]() |
![]() | ![]() | ![]() | ![]() |
![]() | ![]() | ![]() | ![]() |
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].
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.
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).
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).
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.
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.
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.
[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.
FLAGS
LOG FOR LEVEL SET FROM -1 to 1
////////////////////////////////////////////////////// ////////////////// 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.
//////////////////// INITIAL CONDITIONS ////////////////////
///////////////////// BOUNDARY CONDITIONS /////////////////////
////////////// SET VELOCITY //////////////
/////////////////// SET AND GET ALPHA ///////////////////
/////////////// NTH TIME STEP ///////////////
/////// SETUP ///////
//////////////////////////////////// ASSEMBLE MASS (and other) MATRICES ////////////////////////////////////
////////////////////////////////// LOW ORDER METHOD (DiJ Viscosity) //////////////////////////////////
/////////////////// ENTROPY VISCOSITY ///////////////////
/////////////////////// FOR MAXIMUM PRINCIPLE ///////////////////////
/////////////////// COMPUTE SOLUTIONS ///////////////////
/////////// UTILITIES ///////////
/////////////////// MY PETSC WRAPPERS ///////////////////
FINITE ELEMENT SPACE
OPERATORS times SOLUTION VECTOR
MASS MATRIX
BOUNDARIES
////////// 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
SOME PARAMETERS
UTILITIES
/////////////////////////////////////////////////////// /////////////////// 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 //////// /////////////////////////////////////
////////////////////////////// //////// SET VELOCITY //////// //////////////////////////////
SAVE OLD SOLUTION
update velocity
SAVE OLD SOLUTION
update velocity
/////////////////////////////// //////// SET AND GET U //////// ///////////////////////////////
COMPUTE SOLUTION
BOUNDARY CONDITIONS
CHECK MAXIMUM PRINCIPLE
pcout << "*********************************************************************... " << unp1.min() << ", " << unp1.max() << std::endl;
//////////////////////// SETUP FOR DOF HANDLERS //////////////////////// setup system LS
setup system U
////////////////// INIT CONSTRAINTS //////////////////
///////////////////// NON-GHOSTED VECTORS /////////////////////
vectors for lumped mass matrix
operators times solution
LIMITERS (FCT)
///////////////////////////////////////////////////// GHOSTED VECTORS (used within some assemble process) /////////////////////////////////////////////////////
init vectors for vx
init vectors for vy
init vectors for vz
LIMITERS (FCT)
//////////////// SETUP MATRICES //////////////// MATRICES
COMPUTE MASS MATRICES (AND OTHERS) FOR FIRST TIME STEP
get mat for DOFs between Q1 and Q2
distribute
compress
loop on locally owned i-DOFs (rows)
distribute
compress
Distribute
COMPRESS
loop on cells
compute cell_K_times_solution
distribute
K_times_solution=0;
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
loop on sparsity pattern of i-th DOF
ith_K_times_solution += soln[j]*(vj*C);
low order dissipative matrix
high order dissipative matrix (entropy viscosity)
high order compression matrix
save K times solution vector K_times_solution(gi)=ith_K_times_solution; save i-th row of matrices on global matrices
Restore matrices after reading rows
compress K_times_solution.compress(VectorOperation::insert);
get matrices times vector
get solutions at quadrature points
Distribute
ENTROPY NORM FACTOR
loop on locally owned i-DOFs (rows)
get solution at DOFs on the sparsity pattern of i-th DOF
compute bounds, ith row of flux matrix, P vectors
bounds
compute min and max vectors
NON-GHOSTED VECTORS: MPP_uL_solution, NMPP_uH_solution GHOSTED VECTORS: un_solution
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
Array for i-th row of matrices
compute bounds, ith row of flux matrix, P vectors
bounds
i-th row of flux matrix A
compute P vectors
save i-th row of flux matrix A
compute Q vectors
compute R vectors
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
Array for i-th row of A_times_L matrix
loop in sparsity pattern of i-th DOF
save i-th row of LxA
restore A matrix after reading it
loop in num of FCT iterations
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
compute bounds, ith row of flux matrix, P vectors
bounds
compute P vectors
compute Q vectors
compute R vectors
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
Array for i-th row of LxAkp1 matrix
save i-th row of LxA
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)
loop on DOFs
get i-th row of mass matrix (dummy, I just need the indices gj)
all vectors are NON-GHOSTED
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
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
/////////////////////// FOR TRANSPORT PROBLEM /////////////////////// TIME_INTEGRATION
PROBLEM
/////////////////////////////////////////////////// /////////////////// MAIN CLASS //////////////////// ///////////////////////////////////////////////////
SOLUTION VECTORS
BOUNDARY VECTORS
FOR NAVIER STOKES
FOR TRANSPORT
///////////////////////////////////// /////////////// SETUP /////////////// /////////////////////////////////////
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
Initial conditions init condition for phi
init condition for u=0
init condition for v
init condition for p
NO-SLIP CONDITION
LEFT
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
output_vectors();
//////////////////// GENERAL PARAMETERS ////////////////////
////////////////////////////////////////// PARAMETERS FOR THE NAVIER STOKES PROBLEM //////////////////////////////////////////
PROBLEM=FILLING_TANK; PROBLEM=SMALL_WAVE_PERTURBATION; PROBLEM=FALLING_DROP;
////////////////////////////////// PARAMETERS FOR TRANSPORT PROBLEM //////////////////////////////////
TRANSPORT_TIME_INTEGRATION=FORWARD_EULER;
ALGORITHM = "MPP_u1"; ALGORITHM = "NMPP_uH";
ADJUST PARAMETERS ACCORDING TO PROBLEM
////////// GEOMETRY //////////
SETUP
PARAMETERS FOR TIME STEPPING
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
TIME STEPPING
GET NAVIER STOKES VELOCITY
GET LEVEL SET SOLUTION
///////////////////////////////////////////////////////////// /////////////////// NAVIER STOKES SOLVER //////////////////// /////////////////////////////////////////////////////////////
constructor for using LEVEL SET
constructor for NOT LEVEL SET
rho and nu functions
initial conditions
boundary conditions
DO STEPS
SETUP
SETUP AND INITIAL CONDITION
ASSEMBLE SYSTEMS
SOLVERS
GET DIFFERENT FIELDS
OTHERS
CONSTRUCTOR FOR LEVEL SET
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)
get rho, nu
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 /////////// ///////////////////////////////////////////////////
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
sparsity pattern for S
sparsity pattern for M
grl constraints
constraints for dpsi
if (constraints_psi.can_store_line(0)) constraints_psi.add_line(0); //constraint u0 = 0
/////////////////////////////////////////////////// //////////////// ASSEMBLE SYSTEMS ///////////////// ///////////////////////////////////////////////////
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;
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
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];
FORCE TERMS
-2./3*time_step*rho*nonlinearity_u
-2./3*time_step*rho*nonlinearity_v
-2./3*time_step*rho*nonlinearity_w
distribute
BOUNDARY CONDITIONS
PRECONDITIONERS
get function values for LS
get function grads for u and v
Distribute
/////////////////////////////////////////////////// ///////////////////// SOLVERS ///////////////////// ///////////////////////////////////////////////////
PETScWrappers::SolverCG solver(solver_control, mpi_communicator); PETScWrappers::SolverGMRES solver(solver_control, mpi_communicator); PETScWrappers::SolverChebychev solver(solver_control, mpi_communicator);
PETScWrappers::SolverGMRES solver(solver_control, mpi_communicator);
/////////////////////////////////////////////////// ////////////// get different fields /////////////// ///////////////////////////////////////////////////
GET DPSI
SOLVE Q
UPDATE THE PRESSURE
/////////////////////////////////////////////////// ///////////////////// DO STEPS //////////////////// ///////////////////////////////////////////////////
/////////////////////////////////////////////////// ////////////////////// OTHERS ///////////////////// ///////////////////////////////////////////////////
/////////////////////// FOR TRANSPORT PROBLEM /////////////////////// TIME_INTEGRATION
PROBLEM
OTHER FLAGS
/////////////////////////////////////////////////// /////////////////// MAIN CLASS //////////////////// ///////////////////////////////////////////////////
BOUNDARY
VELOCITY
SETUP AND INIT CONDITIONS
POST PROCESSING
SOLUTION VECTORS
BOUNDARY VECTORS
GENERAL
FOR TRANSPORT
FOR RECONSTRUCTION OF MATERIAL FIELDS
MASS MATRIX
VELOCITY //////////
velocity in x
velocity in y
////////// BOUNDARY //////////
/////////////////////////////// SETUP AND INITIAL CONDITIONS //////////////////////////////
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
Initial conditions init condition for phi
Functions::ZeroFunction<dim>(),
init condition for u=0
init condition for v
///////////////// POST PROCESSING /////////////////
error for phi
//////////////////// GENERAL PARAMETERS ////////////////////
PROBLEM=DIAGONAL_ADVECTION;
////////////////////////////////// PARAMETERS FOR TRANSPORT PROBLEM //////////////////////////////////
TRANSPORT_TIME_INTEGRATION=FORWARD_EULER;
ALGORITHM = "MPP_u1";
ALGORITHM = "MPP_uH";
////////// GEOMETRY //////////
GridGenerator::hyper_rectangle(triangulation, Point<dim>(0.0,0.0), Point<dim>(1.0,1.0), true);
/////// SETUP ///////
for Reconstruction of MATERIAL FIELDS
GET TIME STEP
////////////////// TRANSPORT SOLVER //////////////////
/////////////////// INITIAL CONDITION ///////////////////
///////////////////////////// BOUNDARY CONDITIONS FOR PHI /////////////////////////////
OUTPUT DATA REGARDING TIME STEPPING AND MESH
TIME STEPPING
////////////// GET VELOCITY // (NS or interpolate from a function) at current time tn //////////////
SET VELOCITY TO LEVEL SET SOLVER
//////////////////////// GET LEVEL SET SOLUTION // (at tnp1) ////////////////////////
///////////// UPDATE TIME /////////////
//////// OUTPUT ////////
/////////////////////////////////////////////////// /////////////////// MAIN CLASS //////////////////// ///////////////////////////////////////////////////
TimerOutput timer;
timer(std::cout, TimerOutput::summary, TimerOutput::wall_times),
///////////////////////////////////// /////////////// SETUP /////////////// /////////////////////////////////////
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
error for u
error for v
error for w
error for p
PARAMETERS FOR THE NAVIER STOKES PROBLEM
if (cycle==0)
set INITIAL CONDITION within TRANSPORT PROBLEM
TIME STEPPING
/////////////// GET TIME_STEP ///////////////
///////////// FORCE TERMS /////////////
///////////////////////////// DENSITY AND VISCOSITY FIELD /////////////////////////////
///////////////////// BOUNDARY CONDITIONS /////////////////////
////////////// GET SOLUTION //////////////
////////////// FIX PRESSURE //////////////
///////////// UPDATE TIME /////////////
//////// OUTPUT ////////
///////////////////////////////////////////////// ////////////////// INITIAL PHI ////////////////// /////////////////////////////////////////////////
/////////////////////////////////////////////////// ////////////////// FORCE TERMS ///// ////////////// ///////////////////////////////////////////////////
///////////////////////////////////////////////// ////////////////// BOUNDARY PHI ///////////////// /////////////////////////////////////////////////
////////////////////////////////////////////////////// ////////////////// BOUNDARY VELOCITY ///////////////// //////////////////////////////////////////////////////
////////////////// FILLING THE TANK ////////////////// boundary for filling the tank (inlet)
boundary for filling the tank (outlet)
/////////////////////////////////////////////////// ///////////////// POST-PROCESSING ///////////////// ///////////////////////////////////////////////////
/////////////////////////////////////////////////// ////////////////// INITIAL PHI ////////////////// ///////////////////////////////////////////////////
///////////////////////////////////////////////// ////////////////// BOUNDARY PHI ///////////////// /////////////////////////////////////////////////
/////////////////////////////////////////////////// ////////////////// EXACT VELOCITY ///////////////// ///////////////////////////////////////////////////
PROBLEM = 3D_DIAGONAL_ADVECTION
/////////////////////////////////////////////////// ////////// EXACT SOLUTION RHO TO TEST NS ////////// ///////////////////////////////////////////////////
////////////////////////////////////////////////////////////// ///////////////// EXACT SOLUTION U to TEST NS //////////////// //////////////////////////////////////////////////////////////
THIS IS USED JUST FOR TESTING NS
/////////////////////////////////////////////////// ///////// EXACT SOLUTION FOR p TO TEST NS ///////// ///////////////////////////////////////////////////
////////////////////////////////////////////////////////////// ////////////////// FORCE TERMS to TEST NS //////////////////// //////////////////////////////////////////////////////////////
force in x
force in y
force in x