466 *
void get_sparsity_pattern();
470 *
std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
483 *
const std::vector<types::global_dof_index> &indices,
484 *
std::vector<PetscScalar> &values);
486 *
const std::vector<types::global_dof_index> &indices,
488 *
std::vector<PetscScalar> &values);
620 *
double solver_tolerance;
635 *
std::map<types::global_dof_index, std::vector<types::global_dof_index> > sparsity_pattern;
650 *
mpi_communicator (mpi_communicator),
665 *
pcout <<
"********** LEVEL SET SETUP **********" << std::endl;
799 * -------------------------------------------------------------------------------
800 * ------------------------------
COMPUTE SOLUTIONS ------------------------------
801 * -------------------------------------------------------------------------------
836 *
pcout <<
"*********************************************************************... "
837 * <<
unp1.min() <<
", " <<
unp1.max() << std::endl;
845 * --------------------------------------------------------------------
846 * ------------------------------
SETUP ------------------------------
847 * --------------------------------------------------------------------
853 *
solver_tolerance=1
E-6;
881 *
constraints.
clear ();
884 *
constraints.close ();
941 * init vectors
for vx
948 * init vectors
for vy
955 * init vectors
for vz
983 *
dsp, mpi_communicator);
986 *
dsp, mpi_communicator);
989 *
dsp, mpi_communicator);
992 *
dsp, mpi_communicator);
995 *
dsp, mpi_communicator);
1000 *
dsp, mpi_communicator);
1003 *
dsp, mpi_communicator);
1007 *
dsp, mpi_communicator);
1010 *
dsp, mpi_communicator);
1013 *
dsp, mpi_communicator);
1016 *
dsp, mpi_communicator);
1019 *
dsp, mpi_communicator);
1022 *
dsp, mpi_communicator);
1025 *
dsp, mpi_communicator);
1028 *
dsp, mpi_communicator);
1044 *
get_sparsity_pattern();
1049 * ----------------------------------------------------------------------------
1050 * ------------------------------
MASS MATRICES ------------------------------
1051 * ----------------------------------------------------------------------------
1065 *
const unsigned int dofs_per_cell =
fe_LS.dofs_per_cell;
1069 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1076 *
if (
cell_LS->is_locally_owned())
1083 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1091 *
cell_LS->get_dof_indices (local_dof_indices);
1092 *
constraints.distribute_local_to_global (
cell_ML,local_dof_indices,
ML_vector);
1130 *
const unsigned int dofs_per_cell =
fe_LS.dofs_per_cell;
1134 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1135 *
std::vector<double> shape_values(dofs_per_cell);
1142 *
if (
cell_LS->is_locally_owned())
1149 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1152 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1153 *
for (
unsigned int j=0;
j<dofs_per_cell; ++
j)
1154 *
cell_MC(i,
j) += shape_values[i]*shape_values[
j]*JxW;
1161 *
cell_LS->get_dof_indices (local_dof_indices);
1162 *
constraints.distribute_local_to_global (
cell_MC,local_dof_indices,
MC_matrix);
1175 * ---------------------------------------------------------------------------------------
1177 * ---------------------------------------------------------------------------------------
1180 *
template <
int dim>
1216 *
if (
cell_LS->is_locally_owned())
1299 *
const unsigned int dofs_per_cell =
fe_LS.dofs_per_cell;
1304 *
std::vector<Tensor<1,dim> >
un_grads (n_q_points);
1309 *
std::vector<double> shape_values(dofs_per_cell);
1310 *
std::vector<Tensor<1,dim> > shape_grads(dofs_per_cell);
1314 *
std::vector<types::global_dof_index>
indices_LS (dofs_per_cell);
1329 *
if (
cell_LS->is_locally_owned())
1353 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1367 *
template <
int dim>
1415 * get i-
th row
of C matrices
1433 * get vector
values for column indices
1507 * save
K times solution vector
1509 * save i-
th row
of matrices
on global matrices
1545 * get matrices
times vector
1554 * --------------------------------------------------------------------------------------
1555 * ------------------------------
ENTROPY VISCOSITY ------------------------------
1556 * --------------------------------------------------------------------------------------
1559 *
template <
int dim>
1579 *
std::vector<double>
uqn (n_q_points);
1580 *
std::vector<double>
uqnm1 (n_q_points);
1581 *
std::vector<Tensor<1,dim> >
guqn (n_q_points);
1582 *
std::vector<Tensor<1,dim> >
guqnm1 (n_q_points);
1584 *
std::vector<double>
vxqn (n_q_points);
1585 *
std::vector<double>
vyqn (n_q_points);
1586 *
std::vector<double>
vzqn (n_q_points);
1587 *
std::vector<double>
vxqnm1 (n_q_points);
1588 *
std::vector<double>
vyqnm1 (n_q_points);
1589 *
std::vector<double>
vzqnm1 (n_q_points);
1612 *
if (
cell_LS->is_locally_owned())
1641 *
for (
unsigned int q=0;
q<n_q_points; ++
q)
1697 * ------------------------------------------------------------------------------------
1699 * ------------------------------------------------------------------------------------
1722 *
std::vector<types::global_dof_index>
gj_indices = sparsity_pattern[
gi];
1756 *
const unsigned int dofs_per_cell =
fe_LS.dofs_per_cell;
1757 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1767 *
cell_LS->get_dof_indices(local_dof_indices);
1768 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1774 *
pcout <<
"MAX Principle violated" << std::endl;
1783 * -------------------------------------------------------------------------------
1784 * ------------------------------
COMPUTE SOLUTIONS ------------------------------
1785 * -------------------------------------------------------------------------------
1832 *
template <
int dim>
1857 * read vectors at i-
th DOF
1867 * get i-
th row
of matrices
1890 *
std::vector<double> Ai(
ncolumns);
1921 *
Pposi += Ai[
j]*((Ai[
j] > 0) ? 1. : 0.);
1922 *
Pnegi += Ai[
j]*((Ai[
j] < 0) ? 1. : 0.);
2000 * get vector
values for column indices
2064 *
const PetscScalar *
Akp1i;
2066 *
for (
int iter=0; iter<
NUM_ITER; ++iter)
2083 * read vectors at i-
th DOF
2091 * get i-
th row
of matrices
2189 * get vector
values for column indices
2243 *
else if (
algorithm.compare(
"NMPP_uH")==0)
2245 *
else if (
algorithm.compare(
"MPP_uH")==0)
2254 *
pcout <<
"Error in algorithm" << std::endl;
2310 * -----------------------------------------------------------------------
2311 * ------------------------------
UTILITIES ------------------------------
2312 * -----------------------------------------------------------------------
2326 *
const PetscScalar *
MCi;
2337 *
sparsity_pattern[
gi] = std::vector<types::global_dof_index>(
gj,
gj+
ncolumns);
2358 *
if (!
cell_LS->is_artificial())
2367 *
template <
int dim>
2370 *
std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
2384 *
if (
verbose==
true)
pcout <<
" Solved in " << solver_control.last_step() <<
" iterations." << std::endl;
2387 *
template <
int dim>
2394 *
template <
int dim>
2405 * -------------------------------------------------------------------------------
2406 * ------------------------------
MY PETSC WRAPPERS ------------------------------
2407 * -------------------------------------------------------------------------------
2412 *
const std::vector<types::global_dof_index> &indices,
2413 *
std::vector<PetscScalar> &values)
2435 *
n_idx = indices.size();
2442 *
PetscScalar *
soln;
2445 *
for (i = 0; i <
n_idx; ++i)
2447 *
int index = indices[i];
2448 *
if (index >= begin && index < end)
2452 *
const unsigned int ghostidx = ghost_indices.index_within_set(index);
2462 *
const std::vector<types::global_dof_index> &indices,
2463 *
std::map<types::global_dof_index, types::global_dof_index> &
map_from_Q1_to_Q2,
2464 *
std::vector<PetscScalar> &values)
2487 *
n_idx = indices.size();
2494 *
PetscScalar *
soln;
2497 *
for (i = 0; i <
n_idx; ++i)
2500 *
if (index >= begin && index < end)
2504 *
const unsigned int ghostidx = ghost_indices.index_within_set(index);
2515<a name=
"ann-MultiPhase.cc"></a>
2558 *
#
include <
deal.II/distributed/grid_refinement.h>
2567 *
using namespace dealii;
2589 *
#
include "NavierStokesSolver.cc"
2590 *
#
include "LevelSetSolver.cc"
2600 *
template <
int dim>
2673 *
double final_time;
2713 *
template <
int dim>
2733 *
template <
int dim>
2748 *
template <
int dim>
2777 * init vectors
for phi
2785 * init vectors
for u
2793 * init vectors
for v
2801 * init vectors
for p
2815 *
template <
int dim>
2866 *
template <
int dim>
2869 *
constraints.
clear ();
2872 *
constraints.close ();
2875 *
template <
int dim>
2981 *
pcout <<
"Error in type of PROBLEM at Boundary Conditions" << std::endl;
3003 *
template <
int dim>
3020 *
if (
cell_U->is_locally_owned())
3022 *
if (
cell_U->face(face)->at_boundary())
3024 *
fe_face_values.reinit(
cell_U,face);
3029 *
if (fe_face_values.normal_vector(0)*
u < -1e-14)
3030 *
cell_U->face(face)->set_boundary_id(10);
3034 *
template <
int dim>
3066 *
template <
int dim>
3072 *
data_out.build_patches ();
3074 *
const std::string
filename = (
"sol_vectors-" +
3080 *
data_out.write_vtu (output);
3085 *
for (
unsigned int i=0;
3099 *
template <
int dim>
3107 *
data_out.build_patches ();
3109 *
const std::string
filename = (
"sol_rho-" +
3115 *
data_out.write_vtu (output);
3120 *
for (
unsigned int i=0;
3134 *
template <
int dim>
3151 *
final_time = 10.0;
3300 *
mpi_communicator);
3328 *
<<
"; time step: " <<
time_step << std::endl;
3329 *
pcout <<
" Number of active cells: "
3331 *
<<
" Number of degrees of freedom: " << std::endl
3332 *
<<
" U: " <<
dofs_U << std::endl
3333 *
<<
" P: " <<
dofs_P << std::endl
3334 *
<<
" LS: " <<
dofs_LS << std::endl
3347 *
<<
" at t=" << time
3378 *
using namespace dealii;
3391 *
catch (std::exception &exc)
3393 *
std::cerr << std::endl << std::endl
3394 *
<<
"----------------------------------------------------"
3396 *
std::cerr <<
"Exception on processing: " << std::endl
3397 *
<< exc.what() << std::endl
3398 *
<<
"Aborting!" << std::endl
3399 *
<<
"----------------------------------------------------"
3405 *
std::cerr << std::endl << std::endl
3406 *
<<
"----------------------------------------------------"
3408 *
std::cerr <<
"Unknown exception!" << std::endl
3409 *
<<
"Aborting!" << std::endl
3410 *
<<
"----------------------------------------------------"
3419<a name=
"ann-NavierStokesSolver.cc"></a>
3462 *
#
include <
deal.II/distributed/grid_refinement.h>
3474 *
using namespace dealii;
3603 *
std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
3607 *
std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
3740 *
mpi_communicator(mpi_communicator),
3756 * rho_function(force_function),
3757 * nu_function(force_function),
3760 * rho_fluid(rho_fluid),
3761 * nu_fluid(nu_fluid),
3762 * time_step(time_step),
3767 * pcout(std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator)==0)),
3768 * rebuild_Matrix_U(true),
3769 * rebuild_S_M(true),
3770 * rebuild_Matrix_U_preconditioners(true),
3771 * rebuild_S_M_preconditioners(true)
3776 * CONSTRUCTOR NOT FOR LEVEL SET
3780 * NavierStokesSolver<dim>::NavierStokesSolver(const unsigned int degree_LS,
3781 * const unsigned int degree_U,
3782 * const double time_step,
3783 * Function<dim> &force_function,
3784 * Function<dim> &rho_function,
3785 * Function<dim> &nu_function,
3786 * const bool verbose,
3787 * parallel::distributed::Triangulation<dim> &triangulation,
3788 * MPI_Comm &mpi_communicator) :
3789 * mpi_communicator(mpi_communicator),
3790 * triangulation(triangulation),
3791 * degree_LS(degree_LS),
3792 * dof_handler_LS(triangulation),
3794 * degree_U(degree_U),
3795 * dof_handler_U(triangulation),
3797 * dof_handler_P(triangulation),
3799 * force_function(force_function),
3800 * rho_function(rho_function),
3801 * nu_function(nu_function),
3802 * time_step(time_step),
3806 * pcout(std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator)==0)),
3807 * rebuild_Matrix_U(true),
3808 * rebuild_S_M(true),
3809 * rebuild_Matrix_U_preconditioners(true),
3810 * rebuild_S_M_preconditioners(true)
3814 * NavierStokesSolver<dim>::~NavierStokesSolver()
3816 * dof_handler_LS.clear();
3817 * dof_handler_U.clear();
3818 * dof_handler_P.clear();
3823 * /////////////////////////////////////////////////////////
3824 * ////////////////// SETTERS AND GETTERS //////////////////
3825 * /////////////////////////////////////////////////////////
3829 * void NavierStokesSolver<dim>::set_rho_and_nu_functions(const Function<dim> &rho_function,
3830 * const Function<dim> &nu_function)
3832 * this->rho_function=rho_function;
3833 * this->nu_function=nu_function;
3837 * void NavierStokesSolver<dim>::initial_condition(PETScWrappers::MPI::Vector locally_relevant_solution_phi,
3838 * PETScWrappers::MPI::Vector locally_relevant_solution_u,
3839 * PETScWrappers::MPI::Vector locally_relevant_solution_v,
3840 * PETScWrappers::MPI::Vector locally_relevant_solution_p)
3842 * this->locally_relevant_solution_phi=locally_relevant_solution_phi;
3843 * this->locally_relevant_solution_u=locally_relevant_solution_u;
3844 * this->locally_relevant_solution_v=locally_relevant_solution_v;
3845 * this->locally_relevant_solution_p=locally_relevant_solution_p;
3848 * set old vectors to the initial condition (just for first time step)
3851 * save_old_solution();
3855 * void NavierStokesSolver<dim>::initial_condition(PETScWrappers::MPI::Vector locally_relevant_solution_phi,
3856 * PETScWrappers::MPI::Vector locally_relevant_solution_u,
3857 * PETScWrappers::MPI::Vector locally_relevant_solution_v,
3858 * PETScWrappers::MPI::Vector locally_relevant_solution_w,
3859 * PETScWrappers::MPI::Vector locally_relevant_solution_p)
3861 * this->locally_relevant_solution_phi=locally_relevant_solution_phi;
3862 * this->locally_relevant_solution_u=locally_relevant_solution_u;
3863 * this->locally_relevant_solution_v=locally_relevant_solution_v;
3864 * this->locally_relevant_solution_w=locally_relevant_solution_w;
3865 * this->locally_relevant_solution_p=locally_relevant_solution_p;
3868 * set old vectors to the initial condition (just for first time step)
3871 * save_old_solution();
3875 * void NavierStokesSolver<dim>::set_boundary_conditions(std::vector<types::global_dof_index> boundary_values_id_u,
3876 * std::vector<types::global_dof_index> boundary_values_id_v,
3877 * std::vector<double> boundary_values_u,
3878 * std::vector<double> boundary_values_v)
3880 * this->boundary_values_id_u=boundary_values_id_u;
3881 * this->boundary_values_id_v=boundary_values_id_v;
3882 * this->boundary_values_u=boundary_values_u;
3883 * this->boundary_values_v=boundary_values_v;
3887 * void NavierStokesSolver<dim>::set_boundary_conditions(std::vector<types::global_dof_index> boundary_values_id_u,
3888 * std::vector<types::global_dof_index> boundary_values_id_v,
3889 * std::vector<types::global_dof_index> boundary_values_id_w,
3890 * std::vector<double> boundary_values_u,
3891 * std::vector<double> boundary_values_v,
3892 * std::vector<double> boundary_values_w)
3894 * this->boundary_values_id_u=boundary_values_id_u;
3895 * this->boundary_values_id_v=boundary_values_id_v;
3896 * this->boundary_values_id_w=boundary_values_id_w;
3897 * this->boundary_values_u=boundary_values_u;
3898 * this->boundary_values_v=boundary_values_v;
3899 * this->boundary_values_w=boundary_values_w;
3903 * void NavierStokesSolver<dim>::set_velocity(PETScWrappers::MPI::Vector locally_relevant_solution_u,
3904 * PETScWrappers::MPI::Vector locally_relevant_solution_v)
3906 * this->locally_relevant_solution_u=locally_relevant_solution_u;
3907 * this->locally_relevant_solution_v=locally_relevant_solution_v;
3911 * void NavierStokesSolver<dim>::set_velocity(PETScWrappers::MPI::Vector locally_relevant_solution_u,
3912 * PETScWrappers::MPI::Vector locally_relevant_solution_v,
3913 * PETScWrappers::MPI::Vector locally_relevant_solution_w)
3915 * this->locally_relevant_solution_u=locally_relevant_solution_u;
3916 * this->locally_relevant_solution_v=locally_relevant_solution_v;
3917 * this->locally_relevant_solution_w=locally_relevant_solution_w;
3921 * void NavierStokesSolver<dim>::set_phi(PETScWrappers::MPI::Vector locally_relevant_solution_phi)
3923 * this->locally_relevant_solution_phi=locally_relevant_solution_phi;
3927 * void NavierStokesSolver<dim>::get_rho_and_nu(double phi)
3937 * else if (phi<-eps)
3941 * rho_value=rho_fluid*(1+H)/2.+rho_air*(1-H)/2.;
3942 * nu_value=nu_fluid*(1+H)/2.+nu_air*(1-H)/2.;
3945 * rho_value=rho_fluid*(1+phi)/2.+rho_air*(1-phi)/2.;
3946 * nu_value=nu_fluid*(1+phi)/2.+nu_air*(1-phi)/2.;
3952 * void NavierStokesSolver<dim>::get_pressure(PETScWrappers::MPI::Vector &locally_relevant_solution_p)
3954 * locally_relevant_solution_p=this->locally_relevant_solution_p;
3958 * void NavierStokesSolver<dim>::get_velocity(PETScWrappers::MPI::Vector &locally_relevant_solution_u,
3959 * PETScWrappers::MPI::Vector &locally_relevant_solution_v)
3961 * locally_relevant_solution_u=this->locally_relevant_solution_u;
3962 * locally_relevant_solution_v=this->locally_relevant_solution_v;
3966 * void NavierStokesSolver<dim>::get_velocity(PETScWrappers::MPI::Vector &locally_relevant_solution_u,
3967 * PETScWrappers::MPI::Vector &locally_relevant_solution_v,
3968 * PETScWrappers::MPI::Vector &locally_relevant_solution_w)
3970 * locally_relevant_solution_u=this->locally_relevant_solution_u;
3971 * locally_relevant_solution_v=this->locally_relevant_solution_v;
3972 * locally_relevant_solution_w=this->locally_relevant_solution_w;
3977 * ///////////////////////////////////////////////////
3978 * /////////// SETUP AND INITIAL CONDITION ///////////
3979 * ///////////////////////////////////////////////////
3983 * void NavierStokesSolver<dim>::setup()
3985 * pcout<<"***** SETUP IN NAVIER STOKES SOLVER *****"<<std::endl;
3987 * init_constraints();
3992 * void NavierStokesSolver<dim>::setup_DOF()
3995 * degree_MAX=std::max(degree_LS,degree_U);
4001 * dof_handler_LS.distribute_dofs(fe_LS);
4002 * locally_owned_dofs_LS = dof_handler_LS.locally_owned_dofs();
4003 * locally_relevant_dofs_LS = DoFTools::extract_locally_relevant_dofs(dof_handler_LS);
4009 * dof_handler_U.distribute_dofs(fe_U);
4010 * locally_owned_dofs_U = dof_handler_U.locally_owned_dofs();
4011 * locally_relevant_dofs_U = DoFTools::extract_locally_relevant_dofs(dof_handler_U);
4017 * dof_handler_P.distribute_dofs(fe_P);
4018 * locally_owned_dofs_P = dof_handler_P.locally_owned_dofs();
4019 * locally_relevant_dofs_P = DoFTools::extract_locally_relevant_dofs(dof_handler_P);
4023 * void NavierStokesSolver<dim>::setup_VECTORS()
4027 * init vectors for phi
4030 * locally_relevant_solution_phi.reinit(locally_owned_dofs_LS,locally_relevant_dofs_LS,
4031 * mpi_communicator);
4032 * locally_relevant_solution_phi=0;
4035 * init vectors for u
4038 * locally_relevant_solution_u.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4039 * mpi_communicator);
4040 * locally_relevant_solution_u=0;
4041 * completely_distributed_solution_u.reinit(locally_owned_dofs_U,mpi_communicator);
4042 * system_rhs_u.reinit(locally_owned_dofs_U,mpi_communicator);
4045 * init vectors for u_old
4048 * locally_relevant_solution_u_old.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4049 * mpi_communicator);
4050 * locally_relevant_solution_u_old=0;
4053 * init vectors for v
4056 * locally_relevant_solution_v.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4057 * mpi_communicator);
4058 * locally_relevant_solution_v=0;
4059 * completely_distributed_solution_v.reinit(locally_owned_dofs_U,mpi_communicator);
4060 * system_rhs_v.reinit(locally_owned_dofs_U,mpi_communicator);
4063 * init vectors for v_old
4066 * locally_relevant_solution_v_old.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4067 * mpi_communicator);
4068 * locally_relevant_solution_v_old=0;
4071 * init vectors for w
4074 * locally_relevant_solution_w.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4075 * mpi_communicator);
4076 * locally_relevant_solution_w=0;
4077 * completely_distributed_solution_w.reinit(locally_owned_dofs_U,mpi_communicator);
4078 * system_rhs_w.reinit(locally_owned_dofs_U,mpi_communicator);
4081 * init vectors for w_old
4084 * locally_relevant_solution_w_old.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4085 * mpi_communicator);
4086 * locally_relevant_solution_w_old=0;
4089 * init vectors for dpsi
4092 * locally_relevant_solution_psi.reinit(locally_owned_dofs_P,locally_relevant_dofs_P,
4093 * mpi_communicator);
4094 * locally_relevant_solution_psi=0;
4095 * system_rhs_psi.reinit(locally_owned_dofs_P,mpi_communicator);
4098 * init vectors for dpsi old
4101 * locally_relevant_solution_psi_old.reinit(locally_owned_dofs_P,locally_relevant_dofs_P,
4102 * mpi_communicator);
4103 * locally_relevant_solution_psi_old=0;
4106 * init vectors for q
4109 * completely_distributed_solution_q.reinit(locally_owned_dofs_P,mpi_communicator);
4110 * system_rhs_q.reinit(locally_owned_dofs_P,mpi_communicator);
4113 * init vectors for psi
4116 * completely_distributed_solution_psi.reinit(locally_owned_dofs_P,mpi_communicator);
4119 * init vectors for p
4122 * locally_relevant_solution_p.reinit(locally_owned_dofs_P,locally_relevant_dofs_P,
4123 * mpi_communicator);
4124 * locally_relevant_solution_p=0;
4125 * completely_distributed_solution_p.reinit(locally_owned_dofs_P,mpi_communicator);
4128 * ////////////////////////
4129 * Initialize constraints
4130 * ////////////////////////
4133 * init_constraints();
4136 * //////////////////
4138 * //////////////////
4139 * sparsity pattern for A
4142 * DynamicSparsityPattern dsp_Matrix(locally_relevant_dofs_U);
4143 * DoFTools::make_sparsity_pattern(dof_handler_U,dsp_Matrix,constraints,false);
4144 * SparsityTools::distribute_sparsity_pattern(dsp_Matrix,
4145 * dof_handler_U.locally_owned_dofs(),
4147 * locally_relevant_dofs_U);
4148 * system_Matrix_u.reinit(dof_handler_U.locally_owned_dofs(),
4149 * dof_handler_U.locally_owned_dofs(),
4151 * mpi_communicator);
4152 * system_Matrix_v.reinit(dof_handler_U.locally_owned_dofs(),
4153 * dof_handler_U.locally_owned_dofs(),
4155 * mpi_communicator);
4156 * system_Matrix_w.reinit(dof_handler_U.locally_owned_dofs(),
4157 * dof_handler_U.locally_owned_dofs(),
4159 * mpi_communicator);
4160 * rebuild_Matrix_U=true;
4163 * sparsity pattern for S
4166 * DynamicSparsityPattern dsp_S(locally_relevant_dofs_P);
4167 * DoFTools::make_sparsity_pattern(dof_handler_P,dsp_S,constraints_psi,false);
4168 * SparsityTools::distribute_sparsity_pattern(dsp_S,
4169 * dof_handler_P.locally_owned_dofs(),
4171 * locally_relevant_dofs_P);
4172 * system_S.reinit(dof_handler_P.locally_owned_dofs(),
4173 * dof_handler_P.locally_owned_dofs(),
4175 * mpi_communicator);
4178 * sparsity pattern for M
4181 * DynamicSparsityPattern dsp_M(locally_relevant_dofs_P);
4182 * DoFTools::make_sparsity_pattern(dof_handler_P,dsp_M,constraints_psi,false);
4183 * SparsityTools::distribute_sparsity_pattern(dsp_M,
4184 * dof_handler_P.locally_owned_dofs(),
4186 * locally_relevant_dofs_P);
4187 * system_M.reinit(dof_handler_P.locally_owned_dofs(),
4188 * dof_handler_P.locally_owned_dofs(),
4190 * mpi_communicator);
4195 * void NavierStokesSolver<dim>::init_constraints()
4202 * constraints.clear();
4203 * constraints.reinit(locally_relevant_dofs_U);
4204 * DoFTools::make_hanging_node_constraints(dof_handler_U,constraints);
4205 * constraints.close();
4208 * constraints for dpsi
4211 * constraints_psi.clear();
4212 * constraints_psi.reinit(locally_relevant_dofs_P);
4213 * DoFTools::make_hanging_node_constraints(dof_handler_P,constraints_psi);
4216 * if (constraints_psi.can_store_line(0))
4217 * constraints_psi.add_line(0); //constraint u0 = 0
4220 * constraints_psi.close();
4225 * ///////////////////////////////////////////////////
4226 * //////////////// ASSEMBLE SYSTEMS /////////////////
4227 * ///////////////////////////////////////////////////
4231 * void NavierStokesSolver<dim>::assemble_system_U()
4233 * if (rebuild_Matrix_U==true)
4235 * system_Matrix_u=0;
4236 * system_Matrix_v=0;
4237 * system_Matrix_w=0;
4243 * const QGauss<dim> quadrature_formula(degree_MAX+1);
4244 * FEValues<dim> fe_values_LS(fe_LS,quadrature_formula,
4245 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4246 * FEValues<dim> fe_values_U(fe_U,quadrature_formula,
4247 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4248 * FEValues<dim> fe_values_P(fe_P,quadrature_formula,
4249 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4251 * const unsigned int dofs_per_cell=fe_U.dofs_per_cell;
4252 * const unsigned int n_q_points=quadrature_formula.size();
4254 * FullMatrix<double> cell_A_u(dofs_per_cell,dofs_per_cell);
4255 * Vector<double> cell_rhs_u(dofs_per_cell);
4256 * Vector<double> cell_rhs_v(dofs_per_cell);
4257 * Vector<double> cell_rhs_w(dofs_per_cell);
4259 * std::vector<double> phiqnp1(n_q_points);
4261 * std::vector<double> uqn(n_q_points);
4262 * std::vector<double> uqnm1(n_q_points);
4263 * std::vector<double> vqn(n_q_points);
4264 * std::vector<double> vqnm1(n_q_points);
4265 * std::vector<double> wqn(n_q_points);
4266 * std::vector<double> wqnm1(n_q_points);
4270 * FOR Explicit nonlinearity
4271 * std::vector<Tensor<1, dim> > grad_un(n_q_points);
4272 * std::vector<Tensor<1, dim> > grad_vn(n_q_points);
4273 * std::vector<Tensor<1, dim> > grad_wn(n_q_points);
4274 * Tensor<1, dim> Un;
4280 * std::vector<Tensor<1, dim> > grad_pqn(n_q_points);
4281 * std::vector<Tensor<1, dim> > grad_psiqn(n_q_points);
4282 * std::vector<Tensor<1, dim> > grad_psiqnm1(n_q_points);
4284 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
4285 * std::vector<Tensor<1, dim> > shape_grad(dofs_per_cell);
4286 * std::vector<double> shape_value(dofs_per_cell);
4291 * double pressure_grad_u;
4292 * double pressure_grad_v;
4293 * double pressure_grad_w;
4299 * Vector<double> force_terms(dim);
4301 * typename DoFHandler<dim>::active_cell_iterator
4302 * cell_U=dof_handler_U.begin_active(), endc_U=dof_handler_U.end();
4303 * typename DoFHandler<dim>::active_cell_iterator cell_P=dof_handler_P.begin_active();
4304 * typename DoFHandler<dim>::active_cell_iterator cell_LS=dof_handler_LS.begin_active();
4306 * for (; cell_U!=endc_U; ++cell_U,++cell_P,++cell_LS)
4307 * if (cell_U->is_locally_owned())
4314 * fe_values_LS.reinit(cell_LS);
4315 * fe_values_U.reinit(cell_U);
4316 * fe_values_P.reinit(cell_P);
4320 * get function values for LS
4323 * fe_values_LS.get_function_values(locally_relevant_solution_phi,phiqnp1);
4326 * get function values for U
4329 * fe_values_U.get_function_values(locally_relevant_solution_u,uqn);
4330 * fe_values_U.get_function_values(locally_relevant_solution_u_old,uqnm1);
4331 * fe_values_U.get_function_values(locally_relevant_solution_v,vqn);
4332 * fe_values_U.get_function_values(locally_relevant_solution_v_old,vqnm1);
4335 * fe_values_U.get_function_values(locally_relevant_solution_w,wqn);
4336 * fe_values_U.get_function_values(locally_relevant_solution_w_old,wqnm1);
4340 * For explicit nonlinearity
4341 * get gradient values for U
4342 * fe_values_U.get_function_gradients(locally_relevant_solution_u,grad_un);
4343 * fe_values_U.get_function_gradients(locally_relevant_solution_v,grad_vn);
4345 * fe_values_U.get_function_gradients(locally_relevant_solution_w,grad_wn);
4349 * get values and gradients for p and dpsi
4352 * fe_values_P.get_function_gradients(locally_relevant_solution_p,grad_pqn);
4353 * fe_values_P.get_function_gradients(locally_relevant_solution_psi,grad_psiqn);
4354 * fe_values_P.get_function_gradients(locally_relevant_solution_psi_old,grad_psiqnm1);
4356 * for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
4358 * const double JxW=fe_values_U.JxW(q_point);
4359 * for (unsigned int i=0; i<dofs_per_cell; ++i)
4361 * shape_grad[i]=fe_values_U.shape_grad(i,q_point);
4362 * shape_value[i]=fe_values_U.shape_value(i,q_point);
4365 * pressure_grad_u=(grad_pqn[q_point][0]+4./3*grad_psiqn[q_point][0]-1./3*grad_psiqnm1[q_point][0]);
4366 * pressure_grad_v=(grad_pqn[q_point][1]+4./3*grad_psiqn[q_point][1]-1./3*grad_psiqnm1[q_point][1]);
4368 * pressure_grad_w=(grad_pqn[q_point][2]+4./3*grad_psiqn[q_point][2]-1./3*grad_psiqnm1[q_point][2]);
4370 * if (LEVEL_SET==1) // use level set to define rho and nu
4371 * get_rho_and_nu(phiqnp1[q_point]);
4372 * else // rho and nu are defined through functions
4374 * rho_value=rho_function.value(fe_values_U.quadrature_point(q_point));
4375 * nu_value=nu_function.value(fe_values_U.quadrature_point(q_point));
4380 * Non-linearity: for semi-implicit
4383 * u_star=2*uqn[q_point]-uqnm1[q_point];
4384 * v_star=2*vqn[q_point]-vqnm1[q_point];
4386 * w_star=2*wqn[q_point]-wqnm1[q_point];
4390 * for explicit nonlinearity
4391 * Un[0] = uqn[q_point];
4392 * Un[1] = vqn[q_point];
4394 * Un[2] = wqn[q_point];
4398 * double nonlinearity_u = Un*grad_un[q_point];
4399 * double nonlinearity_v = Un*grad_vn[q_point];
4400 * double nonlinearity_w = 0;
4402 * nonlinearity_w = Un*grad_wn[q_point];
4408 * rho_star=rho_value; // This is because we consider rho*u_t instead of (rho*u)_t
4416 * force_function.vector_value(fe_values_U.quadrature_point(q_point),force_terms);
4417 * force_u=force_terms[0];
4418 * force_v=force_terms[1];
4420 * force_w=force_terms[2];
4421 * if (RHO_TIMES_RHS==1)
4429 * for (unsigned int i=0; i<dofs_per_cell; ++i)
4431 * cell_rhs_u(i)+=((4./3*rho*uqn[q_point]-1./3*rho*uqnm1[q_point]
4432 * +2./3*time_step*(force_u-pressure_grad_u)
4435 * -2./3*time_step*rho*nonlinearity_u
4438 * )*shape_value[i])*JxW;
4439 * cell_rhs_v(i)+=((4./3*rho*vqn[q_point]-1./3*rho*vqnm1[q_point]
4440 * +2./3*time_step*(force_v-pressure_grad_v)
4443 * -2./3*time_step*rho*nonlinearity_v
4446 * )*shape_value[i])*JxW;
4448 * cell_rhs_w(i)+=((4./3*rho*wqn[q_point]-1./3*rho*wqnm1[q_point]
4449 * +2./3*time_step*(force_w-pressure_grad_w)
4452 * -2./3*time_step*rho*nonlinearity_w
4455 * )*shape_value[i])*JxW;
4456 * if (rebuild_Matrix_U==true)
4457 * for (unsigned int j=0; j<dofs_per_cell; ++j)
4460 * cell_A_u(i,j)+=(rho_star*shape_value[i]*shape_value[j]
4461 * +2./3*time_step*nu_value*(shape_grad[i]*shape_grad[j])
4462 * +2./3*time_step*rho*shape_value[i]
4463 * *(u_star*shape_grad[j][0]+v_star*shape_grad[j][1]) // semi-implicit NL
4466 * cell_A_u(i,j)+=(rho_star*shape_value[i]*shape_value[j]
4467 * +2./3*time_step*nu_value*(shape_grad[i]*shape_grad[j])
4468 * +2./3*time_step*rho*shape_value[i]
4469 * *(u_star*shape_grad[j][0]+v_star*shape_grad[j][1]+w_star*shape_grad[j][2]) // semi-implicit NL
4474 * cell_U->get_dof_indices(local_dof_indices);
4480 * if (rebuild_Matrix_U==true)
4481 * constraints.distribute_local_to_global(cell_A_u,local_dof_indices,system_Matrix_u);
4482 * constraints.distribute_local_to_global(cell_rhs_u,local_dof_indices,system_rhs_u);
4483 * constraints.distribute_local_to_global(cell_rhs_v,local_dof_indices,system_rhs_v);
4485 * constraints.distribute_local_to_global(cell_rhs_w,local_dof_indices,system_rhs_w);
4487 * system_rhs_u.compress(VectorOperation::add);
4488 * system_rhs_v.compress(VectorOperation::add);
4489 * if (dim==3) system_rhs_w.compress(VectorOperation::add);
4490 * if (rebuild_Matrix_U==true)
4492 * system_Matrix_u.compress(VectorOperation::add);
4493 * system_Matrix_v.copy_from(system_Matrix_u);
4495 * system_Matrix_w.copy_from(system_Matrix_u);
4499 * BOUNDARY CONDITIONS
4502 * system_rhs_u.set(boundary_values_id_u,boundary_values_u);
4503 * system_rhs_u.compress(VectorOperation::insert);
4504 * system_rhs_v.set(boundary_values_id_v,boundary_values_v);
4505 * system_rhs_v.compress(VectorOperation::insert);
4508 * system_rhs_w.set(boundary_values_id_w,boundary_values_w);
4509 * system_rhs_w.compress(VectorOperation::insert);
4511 * if (rebuild_Matrix_U)
4513 * system_Matrix_u.clear_rows(boundary_values_id_u,1);
4514 * system_Matrix_v.clear_rows(boundary_values_id_v,1);
4516 * system_Matrix_w.clear_rows(boundary_values_id_w,1);
4517 * if (rebuild_Matrix_U_preconditioners)
4524 * rebuild_Matrix_U_preconditioners=false;
4525 * preconditioner_Matrix_u.reset(new PETScWrappers::PreconditionBoomerAMG
4526 * (system_Matrix_u,PETScWrappers::PreconditionBoomerAMG::AdditionalData(false)));
4527 * preconditioner_Matrix_v.reset( new PETScWrappers::PreconditionBoomerAMG
4528 * (system_Matrix_v,PETScWrappers::PreconditionBoomerAMG::AdditionalData(false)));
4530 * preconditioner_Matrix_w.reset(new PETScWrappers::PreconditionBoomerAMG
4531 * (system_Matrix_w,PETScWrappers::PreconditionBoomerAMG::AdditionalData(false)));
4534 * rebuild_Matrix_U=true;
4538 * void NavierStokesSolver<dim>::assemble_system_dpsi_q()
4540 * if (rebuild_S_M==true)
4548 * const QGauss<dim> quadrature_formula(degree_MAX+1);
4550 * FEValues<dim> fe_values_U(fe_U,quadrature_formula,
4551 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4552 * FEValues<dim> fe_values_P(fe_P,quadrature_formula,
4553 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4554 * FEValues<dim> fe_values_LS(fe_LS,quadrature_formula,
4555 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4557 * const unsigned int dofs_per_cell=fe_P.dofs_per_cell;
4558 * const unsigned int n_q_points=quadrature_formula.size();
4560 * FullMatrix<double> cell_S(dofs_per_cell,dofs_per_cell);
4561 * FullMatrix<double> cell_M(dofs_per_cell,dofs_per_cell);
4562 * Vector<double> cell_rhs_psi(dofs_per_cell);
4563 * Vector<double> cell_rhs_q(dofs_per_cell);
4565 * std::vector<double> phiqnp1(n_q_points);
4566 * std::vector<Tensor<1, dim> > gunp1(n_q_points);
4567 * std::vector<Tensor<1, dim> > gvnp1(n_q_points);
4568 * std::vector<Tensor<1, dim> > gwnp1(n_q_points);
4570 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
4571 * std::vector<double> shape_value(dofs_per_cell);
4572 * std::vector<Tensor<1, dim> > shape_grad(dofs_per_cell);
4574 * typename DoFHandler<dim>::active_cell_iterator
4575 * cell_P=dof_handler_P.begin_active(), endc_P=dof_handler_P.end();
4576 * typename DoFHandler<dim>::active_cell_iterator cell_U=dof_handler_U.begin_active();
4577 * typename DoFHandler<dim>::active_cell_iterator cell_LS=dof_handler_LS.begin_active();
4579 * for (; cell_P!=endc_P; ++cell_P,++cell_U,++cell_LS)
4580 * if (cell_P->is_locally_owned())
4587 * fe_values_P.reinit(cell_P);
4588 * fe_values_U.reinit(cell_U);
4589 * fe_values_LS.reinit(cell_LS);
4593 * get function values for LS
4596 * fe_values_LS.get_function_values(locally_relevant_solution_phi,phiqnp1);
4600 * get function grads for u and v
4603 * fe_values_U.get_function_gradients(locally_relevant_solution_u,gunp1);
4604 * fe_values_U.get_function_gradients(locally_relevant_solution_v,gvnp1);
4606 * fe_values_U.get_function_gradients(locally_relevant_solution_w,gwnp1);
4608 * for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
4610 * const double JxW=fe_values_P.JxW(q_point);
4611 * double divU = gunp1[q_point][0]+gvnp1[q_point][1];
4612 * if (dim==3) divU += gwnp1[q_point][2];
4613 * for (unsigned int i=0; i<dofs_per_cell; ++i)
4615 * shape_value[i]=fe_values_P.shape_value(i,q_point);
4616 * shape_grad[i]=fe_values_P.shape_grad(i,q_point);
4618 * if (LEVEL_SET==1) // use level set to define rho and nu
4619 * get_rho_and_nu (phiqnp1[q_point]);
4620 * else // rho and nu are defined through functions
4621 * nu_value=nu_function.value(fe_values_U.quadrature_point(q_point));
4623 * for (unsigned int i=0; i<dofs_per_cell; ++i)
4625 * cell_rhs_psi(i)+=-3./2./time_step*rho_min*divU*shape_value[i]*JxW;
4626 * cell_rhs_q(i)-=nu_value*divU*shape_value[i]*JxW;
4627 * if (rebuild_S_M==true)
4629 * for (unsigned int j=0; j<dofs_per_cell; ++j)
4632 * cell_S(i,j)+=shape_grad[i]*shape_grad[j]*JxW+1E-10;
4633 * cell_M(i,j)+=shape_value[i]*shape_value[j]*JxW;
4637 * cell_S(i,j)+=shape_grad[i]*shape_grad[j]*JxW;
4638 * cell_M(i,j)+=shape_value[i]*shape_value[j]*JxW;
4643 * cell_P->get_dof_indices(local_dof_indices);
4649 * if (rebuild_S_M==true)
4651 * constraints_psi.distribute_local_to_global(cell_S,local_dof_indices,system_S);
4652 * constraints_psi.distribute_local_to_global(cell_M,local_dof_indices,system_M);
4654 * constraints_psi.distribute_local_to_global(cell_rhs_q,local_dof_indices,system_rhs_q);
4655 * constraints_psi.distribute_local_to_global(cell_rhs_psi,local_dof_indices,system_rhs_psi);
4657 * if (rebuild_S_M==true)
4659 * system_M.compress(VectorOperation::add);
4660 * system_S.compress(VectorOperation::add);
4661 * if (rebuild_S_M_preconditioners)
4663 * rebuild_S_M_preconditioners=false;
4664 * preconditioner_S.reset(new PETScWrappers::PreconditionBoomerAMG
4665 * (system_S,PETScWrappers::PreconditionBoomerAMG::AdditionalData(true)));
4666 * preconditioner_M.reset(new PETScWrappers::PreconditionBoomerAMG
4667 * (system_M,PETScWrappers::PreconditionBoomerAMG::AdditionalData(true)));
4670 * system_rhs_psi.compress(VectorOperation::add);
4671 * system_rhs_q.compress(VectorOperation::add);
4672 * rebuild_S_M=false;
4677 * ///////////////////////////////////////////////////
4678 * ///////////////////// SOLVERS /////////////////////
4679 * ///////////////////////////////////////////////////
4683 * void NavierStokesSolver<dim>::solve_U(const AffineConstraints<double> &constraints,
4684 * PETScWrappers::MPI::SparseMatrix &Matrix,
4685 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
4686 * PETScWrappers::MPI::Vector &completely_distributed_solution,
4687 * const PETScWrappers::MPI::Vector &rhs)
4689 * SolverControl solver_control(dof_handler_U.n_dofs(),1e-6);
4692 * PETScWrappers::SolverCG solver(solver_control, mpi_communicator);
4693 * PETScWrappers::SolverGMRES solver(solver_control, mpi_communicator);
4694 * PETScWrappers::SolverChebychev solver(solver_control, mpi_communicator);
4697 * PETScWrappers::SolverBicgstab solver(solver_control,mpi_communicator);
4698 * constraints.distribute(completely_distributed_solution);
4699 * solver.solve(Matrix,completely_distributed_solution,rhs,*preconditioner);
4700 * constraints.distribute(completely_distributed_solution);
4701 * if (solver_control.last_step() > MAX_NUM_ITER_TO_RECOMPUTE_PRECONDITIONER)
4702 * rebuild_Matrix_U_preconditioners=true;
4703 * if (verbose==true)
4704 * pcout<<" Solved U in "<<solver_control.last_step()<<" iterations."<<std::endl;
4708 * void NavierStokesSolver<dim>::solve_P(const AffineConstraints<double> &constraints,
4709 * PETScWrappers::MPI::SparseMatrix &Matrix,
4710 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
4711 * PETScWrappers::MPI::Vector &completely_distributed_solution,
4712 * const PETScWrappers::MPI::Vector &rhs)
4714 * SolverControl solver_control(dof_handler_P.n_dofs(),1e-6);
4715 * PETScWrappers::SolverCG solver(solver_control,mpi_communicator);
4718 * PETScWrappers::SolverGMRES solver(solver_control, mpi_communicator);
4721 * constraints.distribute(completely_distributed_solution);
4722 * solver.solve(Matrix,completely_distributed_solution,rhs,*preconditioner);
4723 * constraints.distribute(completely_distributed_solution);
4724 * if (solver_control.last_step() > MAX_NUM_ITER_TO_RECOMPUTE_PRECONDITIONER)
4725 * rebuild_S_M_preconditioners=true;
4726 * if (verbose==true)
4727 * pcout<<" Solved P in "<<solver_control.last_step()<<" iterations."<<std::endl;
4732 * ///////////////////////////////////////////////////
4733 * ////////////// get different fields ///////////////
4734 * ///////////////////////////////////////////////////
4738 * void NavierStokesSolver<dim>::get_velocity()
4740 * assemble_system_U();
4741 * save_old_solution();
4742 * solve_U(constraints,system_Matrix_u,preconditioner_Matrix_u,completely_distributed_solution_u,system_rhs_u);
4743 * locally_relevant_solution_u=completely_distributed_solution_u;
4744 * solve_U(constraints,system_Matrix_v,preconditioner_Matrix_v,completely_distributed_solution_v,system_rhs_v);
4745 * locally_relevant_solution_v=completely_distributed_solution_v;
4748 * solve_U(constraints,system_Matrix_w,preconditioner_Matrix_w,completely_distributed_solution_w,system_rhs_w);
4749 * locally_relevant_solution_w=completely_distributed_solution_w;
4754 * void NavierStokesSolver<dim>::get_pressure()
4761 * assemble_system_dpsi_q();
4762 * solve_P(constraints_psi,system_S,preconditioner_S,completely_distributed_solution_psi,system_rhs_psi);
4763 * locally_relevant_solution_psi=completely_distributed_solution_psi;
4769 * solve_P(constraints,system_M,preconditioner_M,completely_distributed_solution_q,system_rhs_q);
4772 * UPDATE THE PRESSURE
4775 * completely_distributed_solution_p.add(1,completely_distributed_solution_psi);
4776 * completely_distributed_solution_p.add(1,completely_distributed_solution_q);
4777 * locally_relevant_solution_p = completely_distributed_solution_p;
4782 * ///////////////////////////////////////////////////
4783 * ///////////////////// DO STEPS ////////////////////
4784 * ///////////////////////////////////////////////////
4788 * void NavierStokesSolver<dim>::nth_time_step()
4796 * ///////////////////////////////////////////////////
4797 * ////////////////////// OTHERS /////////////////////
4798 * ///////////////////////////////////////////////////
4802 * void NavierStokesSolver<dim>::save_old_solution()
4804 * locally_relevant_solution_u_old=locally_relevant_solution_u;
4805 * locally_relevant_solution_v_old=locally_relevant_solution_v;
4806 * locally_relevant_solution_w_old=locally_relevant_solution_w;
4807 * locally_relevant_solution_psi_old=locally_relevant_solution_psi;
4813<a name="ann-TestLevelSet.cc"></a>
4814<h1>Annotated version of TestLevelSet.cc</h1>
4820 * /* -----------------------------------------------------------------------------
4822 * * SPDX-License-Identifier: LGPL-2.1-or-later
4823 * * Copyright (C) 2016 Manuel Quezada de Luna
4825 * * This file is part of the deal.II code gallery.
4827 * * -----------------------------------------------------------------------------
4830 * #include <deal.II/base/quadrature_lib.h>
4831 * #include <deal.II/base/function.h>
4832 * #include <deal.II/lac/affine_constraints.h>
4833 * #include <deal.II/lac/vector.h>
4834 * #include <deal.II/lac/full_matrix.h>
4835 * #include <deal.II/lac/solver_cg.h>
4836 * #include <deal.II/lac/petsc_sparse_matrix.h>
4837 * #include <deal.II/lac/petsc_vector.h>
4838 * #include <deal.II/lac/petsc_solver.h>
4839 * #include <deal.II/lac/petsc_precondition.h>
4840 * #include <deal.II/grid/grid_generator.h>
4841 * #include <deal.II/grid/tria_accessor.h>
4842 * #include <deal.II/grid/tria_iterator.h>
4843 * #include <deal.II/dofs/dof_handler.h>
4844 * #include <deal.II/dofs/dof_accessor.h>
4845 * #include <deal.II/dofs/dof_tools.h>
4846 * #include <deal.II/fe/fe_values.h>
4847 * #include <deal.II/fe/fe_q.h>
4848 * #include <deal.II/numerics/vector_tools.h>
4849 * #include <deal.II/numerics/data_out.h>
4850 * #include <deal.II/numerics/error_estimator.h>
4851 * #include <deal.II/base/utilities.h>
4852 * #include <deal.II/base/conditional_ostream.h>
4853 * #include <deal.II/base/index_set.h>
4854 * #include <deal.II/lac/sparsity_tools.h>
4855 * #include <deal.II/distributed/tria.h>
4856 * #include <deal.II/distributed/grid_refinement.h>
4857 * #include <deal.II/lac/petsc_vector.h>
4858 * #include <deal.II/base/convergence_table.h>
4859 * #include <deal.II/base/timer.h>
4860 * #include <deal.II/base/parameter_handler.h>
4861 * #include <deal.II/grid/grid_tools.h>
4862 * #include <deal.II/fe/mapping_q.h>
4863 * #include <deal.II/fe/fe_system.h>
4865 * #include <fstream>
4866 * #include <iostream>
4869 * using namespace dealii;
4873 * ///////////////////////
4874 * FOR TRANSPORT PROBLEM
4875 * ///////////////////////
4879 * #define FORWARD_EULER 0
4886 * #define CIRCULAR_ROTATION 0
4887 * #define DIAGONAL_ADVECTION 1
4893 * #define VARIABLE_VELOCITY 0
4895 * #include "utilities_test_LS.cc"
4896 * #include "LevelSetSolver.cc"
4900 * ///////////////////////////////////////////////////
4901 * /////////////////// MAIN CLASS ////////////////////
4902 * ///////////////////////////////////////////////////
4905 * template <int dim>
4906 * class TestLevelSet
4909 * TestLevelSet (const unsigned int degree_LS,
4910 * const unsigned int degree_U);
4920 * void set_boundary_inlet();
4921 * void get_boundary_values_phi(std::vector<unsigned int> &boundary_values_id_phi,
4922 * std::vector<double> &boundary_values_phi);
4928 * void get_interpolated_velocity();
4931 * SETUP AND INIT CONDITIONS
4935 * void initial_condition();
4936 * void init_constraints();
4942 * void process_solution(parallel::distributed::Triangulation<dim> &triangulation,
4943 * DoFHandler<dim> &dof_handler_LS,
4944 * PETScWrappers::MPI::Vector &solution);
4945 * void output_results();
4946 * void output_solution();
4953 * PETScWrappers::MPI::Vector locally_relevant_solution_phi;
4954 * PETScWrappers::MPI::Vector locally_relevant_solution_u;
4955 * PETScWrappers::MPI::Vector locally_relevant_solution_v;
4956 * PETScWrappers::MPI::Vector locally_relevant_solution_w;
4957 * PETScWrappers::MPI::Vector completely_distributed_solution_phi;
4958 * PETScWrappers::MPI::Vector completely_distributed_solution_u;
4959 * PETScWrappers::MPI::Vector completely_distributed_solution_v;
4960 * PETScWrappers::MPI::Vector completely_distributed_solution_w;
4966 * std::vector<unsigned int> boundary_values_id_phi;
4967 * std::vector<double> boundary_values_phi;
4974 * MPI_Comm mpi_communicator;
4975 * parallel::distributed::Triangulation<dim> triangulation;
4979 * DoFHandler<dim> dof_handler_LS;
4981 * IndexSet locally_owned_dofs_LS;
4982 * IndexSet locally_relevant_dofs_LS;
4985 * DoFHandler<dim> dof_handler_U;
4987 * IndexSet locally_owned_dofs_U;
4988 * IndexSet locally_relevant_dofs_U;
4990 * DoFHandler<dim> dof_handler_U_disp_field;
4991 * FESystem<dim> fe_U_disp_field;
4992 * IndexSet locally_owned_dofs_U_disp_field;
4993 * IndexSet locally_relevant_dofs_U_disp_field;
4995 * AffineConstraints<double> constraints;
4996 * AffineConstraints<double> constraints_disp_field;
5000 * double final_time;
5001 * unsigned int timestep_number;
5006 * int sharpness_integer;
5008 * unsigned int n_refinement;
5009 * unsigned int output_number;
5010 * double output_time;
5014 * ConditionalOStream pcout;
5021 * double cK; //compression coeff
5022 * double cE; //entropy-visc coeff
5023 * unsigned int TRANSPORT_TIME_INTEGRATION;
5024 * std::string ALGORITHM;
5025 * unsigned int PROBLEM;
5029 * FOR RECONSTRUCTION OF MATERIAL FIELDS
5032 * double eps, rho_air, rho_fluid;
5039 * PETScWrappers::MPI::SparseMatrix matrix_MC, matrix_MC_tnm1;
5040 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_MC;
5044 * template <int dim>
5045 * TestLevelSet<dim>::TestLevelSet (const unsigned int degree_LS,
5046 * const unsigned int degree_U)
5048 * mpi_communicator (MPI_COMM_WORLD),
5049 * triangulation (mpi_communicator,
5050 * typename Triangulation<dim>::MeshSmoothing
5051 * (Triangulation<dim>::smoothing_on_refinement |
5052 * Triangulation<dim>::smoothing_on_coarsening)),
5053 * degree_LS(degree_LS),
5054 * dof_handler_LS (triangulation),
5055 * fe_LS (degree_LS),
5056 * degree_U(degree_U),
5057 * dof_handler_U (triangulation),
5059 * dof_handler_U_disp_field(triangulation),
5060 * fe_U_disp_field(FE_Q<dim>(degree_U),dim),
5061 * pcout (std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator)== 0))
5064 * template <int dim>
5065 * TestLevelSet<dim>::~TestLevelSet ()
5067 * dof_handler_U_disp_field.clear();
5068 * dof_handler_LS.clear ();
5069 * dof_handler_U.clear ();
5078 * template <int dim>
5079 * void TestLevelSet<dim>::get_interpolated_velocity()
5086 * completely_distributed_solution_u = 0;
5087 * VectorTools::interpolate(dof_handler_U,
5088 * ExactU<dim>(PROBLEM,time),
5089 * completely_distributed_solution_u);
5090 * constraints.distribute (completely_distributed_solution_u);
5091 * locally_relevant_solution_u = completely_distributed_solution_u;
5097 * completely_distributed_solution_v = 0;
5098 * VectorTools::interpolate(dof_handler_U,
5099 * ExactV<dim>(PROBLEM,time),
5100 * completely_distributed_solution_v);
5101 * constraints.distribute (completely_distributed_solution_v);
5102 * locally_relevant_solution_v = completely_distributed_solution_v;
5105 * completely_distributed_solution_w = 0;
5106 * VectorTools::interpolate(dof_handler_U,
5107 * ExactW<dim>(PROBLEM,time),
5108 * completely_distributed_solution_w);
5109 * constraints.distribute (completely_distributed_solution_w);
5110 * locally_relevant_solution_w = completely_distributed_solution_w;
5121 * template <int dim>
5122 * void TestLevelSet<dim>::set_boundary_inlet()
5124 * const QGauss<dim-1> face_quadrature_formula(1); // center of the face
5125 * FEFaceValues<dim> fe_face_values (fe_U,face_quadrature_formula,
5126 * update_values | update_quadrature_points |
5127 * update_normal_vectors);
5128 * const unsigned int n_face_q_points = face_quadrature_formula.size();
5129 * std::vector<double> u_value (n_face_q_points);
5130 * std::vector<double> v_value (n_face_q_points);
5131 * std::vector<double> w_value (n_face_q_points);
5133 * typename DoFHandler<dim>::active_cell_iterator
5134 * cell_U = dof_handler_U.begin_active(),
5135 * endc_U = dof_handler_U.end();
5138 * for (; cell_U!=endc_U; ++cell_U)
5139 * if (cell_U->is_locally_owned())
5140 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
5141 * if (cell_U->face(face)->at_boundary())
5143 * fe_face_values.reinit(cell_U,face);
5144 * fe_face_values.get_function_values(locally_relevant_solution_u,u_value);
5145 * fe_face_values.get_function_values(locally_relevant_solution_v,v_value);
5147 * fe_face_values.get_function_values(locally_relevant_solution_w,w_value);
5152 * if (fe_face_values.normal_vector(0)*u < -1e-14)
5153 * cell_U->face(face)->set_boundary_id(10);
5157 * template <int dim>
5158 * void TestLevelSet<dim>::get_boundary_values_phi(std::vector<unsigned int> &boundary_values_id_phi,
5159 * std::vector<double> &boundary_values_phi)
5161 * std::map<unsigned int, double> map_boundary_values_phi;
5162 * unsigned int boundary_id=0;
5164 * set_boundary_inlet();
5165 * boundary_id=10; // inlet
5166 * VectorTools::interpolate_boundary_values (dof_handler_LS,
5167 * boundary_id,BoundaryPhi<dim>(),
5168 * map_boundary_values_phi);
5170 * boundary_values_id_phi.resize(map_boundary_values_phi.size());
5171 * boundary_values_phi.resize(map_boundary_values_phi.size());
5172 * std::map<unsigned int,double>::const_iterator boundary_value_phi = map_boundary_values_phi.begin();
5173 * for (int i=0; boundary_value_phi !=map_boundary_values_phi.end(); ++boundary_value_phi, ++i)
5175 * boundary_values_id_phi[i]=boundary_value_phi->first;
5176 * boundary_values_phi[i]=boundary_value_phi->second;
5182 * ///////////////////////////////
5183 * SETUP AND INITIAL CONDITIONS
5184 * //////////////////////////////
5187 * template <int dim>
5188 * void TestLevelSet<dim>::setup()
5190 * degree = std::max(degree_LS,degree_U);
5196 * dof_handler_LS.distribute_dofs (fe_LS);
5197 * locally_owned_dofs_LS = dof_handler_LS.locally_owned_dofs ();
5198 * locally_relevant_dofs_LS = DoFTools::extract_locally_relevant_dofs (dof_handler_LS);
5204 * dof_handler_U.distribute_dofs (fe_U);
5205 * locally_owned_dofs_U = dof_handler_U.locally_owned_dofs ();
5206 * locally_relevant_dofs_U = DoFTools::extract_locally_relevant_dofs (dof_handler_U);
5209 * setup system U for disp field
5212 * dof_handler_U_disp_field.distribute_dofs (fe_U_disp_field);
5213 * locally_owned_dofs_U_disp_field = dof_handler_U_disp_field.locally_owned_dofs ();
5214 * locally_relevant_dofs_U_disp_field = DoFTools::extract_locally_relevant_dofs (dof_handler_U_disp_field);
5217 * init vectors for phi
5220 * locally_relevant_solution_phi.reinit(locally_owned_dofs_LS,
5221 * locally_relevant_dofs_LS,
5222 * mpi_communicator);
5223 * locally_relevant_solution_phi = 0;
5224 * completely_distributed_solution_phi.reinit(mpi_communicator,
5225 * dof_handler_LS.n_dofs(),
5226 * dof_handler_LS.n_locally_owned_dofs());
5229 * init vectors for u
5232 * locally_relevant_solution_u.reinit(locally_owned_dofs_U,
5233 * locally_relevant_dofs_U,
5234 * mpi_communicator);
5235 * locally_relevant_solution_u = 0;
5236 * completely_distributed_solution_u.reinit(mpi_communicator,
5237 * dof_handler_U.n_dofs(),
5238 * dof_handler_U.n_locally_owned_dofs());
5241 * init vectors for v
5244 * locally_relevant_solution_v.reinit(locally_owned_dofs_U,
5245 * locally_relevant_dofs_U,
5246 * mpi_communicator);
5247 * locally_relevant_solution_v = 0;
5248 * completely_distributed_solution_v.reinit(mpi_communicator,
5249 * dof_handler_U.n_dofs(),
5250 * dof_handler_U.n_locally_owned_dofs());
5253 * init vectors for w
5256 * locally_relevant_solution_w.reinit(locally_owned_dofs_U,
5257 * locally_relevant_dofs_U,
5258 * mpi_communicator);
5259 * locally_relevant_solution_w = 0;
5260 * completely_distributed_solution_w.reinit(mpi_communicator,
5261 * dof_handler_U.n_dofs(),
5262 * dof_handler_U.n_locally_owned_dofs());
5263 * init_constraints();
5269 * DynamicSparsityPattern dsp (locally_relevant_dofs_LS);
5270 * DoFTools::make_sparsity_pattern (dof_handler_LS,dsp,constraints,false);
5271 * SparsityTools::distribute_sparsity_pattern (dsp,
5272 * dof_handler_LS.n_locally_owned_dofs_per_processor(),
5274 * locally_relevant_dofs_LS);
5275 * matrix_MC.reinit (mpi_communicator,
5277 * dof_handler_LS.n_locally_owned_dofs_per_processor(),
5278 * dof_handler_LS.n_locally_owned_dofs_per_processor(),
5279 * Utilities::MPI::this_mpi_process(mpi_communicator));
5280 * matrix_MC_tnm1.reinit (mpi_communicator,
5282 * dof_handler_LS.n_locally_owned_dofs_per_processor(),
5283 * dof_handler_LS.n_locally_owned_dofs_per_processor(),
5284 * Utilities::MPI::this_mpi_process(mpi_communicator));
5287 * template <int dim>
5288 * void TestLevelSet<dim>::initial_condition()
5293 * Initial conditions
5294 * init condition for phi
5297 * completely_distributed_solution_phi = 0;
5298 * VectorTools::interpolate(dof_handler_LS,
5299 * InitialPhi<dim>(PROBLEM, sharpness),
5302 * Functions::ZeroFunction<dim>(),
5305 * completely_distributed_solution_phi);
5306 * constraints.distribute (completely_distributed_solution_phi);
5307 * locally_relevant_solution_phi = completely_distributed_solution_phi;
5310 * init condition for u=0
5313 * completely_distributed_solution_u = 0;
5314 * VectorTools::interpolate(dof_handler_U,
5315 * ExactU<dim>(PROBLEM,time),
5316 * completely_distributed_solution_u);
5317 * constraints.distribute (completely_distributed_solution_u);
5318 * locally_relevant_solution_u = completely_distributed_solution_u;
5321 * init condition for v
5324 * completely_distributed_solution_v = 0;
5325 * VectorTools::interpolate(dof_handler_U,
5326 * ExactV<dim>(PROBLEM,time),
5327 * completely_distributed_solution_v);
5328 * constraints.distribute (completely_distributed_solution_v);
5329 * locally_relevant_solution_v = completely_distributed_solution_v;
5332 * template <int dim>
5333 * void TestLevelSet<dim>::init_constraints()
5335 * constraints.clear ();
5336 * constraints.reinit (locally_relevant_dofs_LS);
5337 * DoFTools::make_hanging_node_constraints (dof_handler_LS, constraints);
5338 * constraints.close ();
5339 * constraints_disp_field.clear ();
5340 * constraints_disp_field.reinit (locally_relevant_dofs_LS);
5341 * DoFTools::make_hanging_node_constraints (dof_handler_LS, constraints_disp_field);
5342 * constraints_disp_field.close ();
5352 * template <int dim>
5353 * void TestLevelSet<dim>::process_solution(parallel::distributed::Triangulation<dim> &triangulation,
5354 * DoFHandler<dim> &dof_handler_LS,
5355 * PETScWrappers::MPI::Vector &solution)
5357 * Vector<double> difference_per_cell (triangulation.n_active_cells());
5363 * VectorTools::integrate_difference (dof_handler_LS,
5365 * InitialPhi<dim>(PROBLEM,sharpness),
5366 * difference_per_cell,
5367 * QGauss<dim>(degree_LS+3),
5368 * VectorTools::L1_norm);
5370 * double u_L1_error = difference_per_cell.l1_norm();
5371 * u_L1_error = std::sqrt(Utilities::MPI::sum(u_L1_error * u_L1_error, mpi_communicator));
5373 * VectorTools::integrate_difference (dof_handler_LS,
5375 * InitialPhi<dim>(PROBLEM,sharpness),
5376 * difference_per_cell,
5377 * QGauss<dim>(degree_LS+3),
5378 * VectorTools::L2_norm);
5379 * double u_L2_error = difference_per_cell.l2_norm();
5380 * u_L2_error = std::sqrt(Utilities::MPI::sum(u_L2_error * u_L2_error, mpi_communicator));
5382 * pcout << "L1 error: " << u_L1_error << std::endl;
5383 * pcout << "L2 error: " << u_L2_error << std::endl;
5387 * void TestLevelSet<dim>::output_results()
5389 * output_solution();
5393 * template <int dim>
5394 * void TestLevelSet<dim>::output_solution()
5396 * DataOut<dim> data_out;
5397 * data_out.attach_dof_handler(dof_handler_LS);
5398 * data_out.add_data_vector (locally_relevant_solution_phi, "phi");
5399 * data_out.build_patches();
5401 * const std::string filename = ("solution-" +
5402 * Utilities::int_to_string (output_number, 3) +
5404 * Utilities::int_to_string
5405 * (triangulation.locally_owned_subdomain(), 4));
5406 * std::ofstream output ((filename + ".vtu").c_str());
5407 * data_out.write_vtu (output);
5409 * if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
5411 * std::vector<std::string> filenames;
5412 * for (unsigned int i=0;
5413 * i<Utilities::MPI::n_mpi_processes(mpi_communicator);
5415 * filenames.push_back ("solution-" +
5416 * Utilities::int_to_string (output_number, 3) +
5418 * Utilities::int_to_string (i, 4) +
5421 * std::ofstream master_output ((filename + ".pvtu").c_str());
5422 * data_out.write_pvtu_record (master_output, filenames);
5426 * template <int dim>
5427 * void TestLevelSet<dim>::run()
5431 * ////////////////////
5432 * GENERAL PARAMETERS
5433 * ////////////////////
5438 * get_output = true;
5439 * output_number = 0;
5442 * output_time = 0.1;
5444 * PROBLEM=CIRCULAR_ROTATION;
5447 * PROBLEM=DIAGONAL_ADVECTION;
5451 * if (PROBLEM==CIRCULAR_ROTATION)
5452 * umax = std::sqrt(2)*numbers::PI;
5454 * umax = std::sqrt(2);
5458 * //////////////////////////////////
5459 * PARAMETERS FOR TRANSPORT PROBLEM
5460 * //////////////////////////////////
5463 * cK = 1.0; // compression constant
5464 * cE = 1.0; // entropy viscosity constant
5465 * sharpness_integer=1; //this will be multiplied by min_h
5468 * TRANSPORT_TIME_INTEGRATION=FORWARD_EULER;
5471 * TRANSPORT_TIME_INTEGRATION=SSP33;
5474 * ALGORITHM = "MPP_u1";
5477 * ALGORITHM = "NMPP_uH";
5480 * ALGORITHM = "MPP_uH";
5489 * if (PROBLEM==CIRCULAR_ROTATION || PROBLEM==DIAGONAL_ADVECTION)
5490 * GridGenerator::hyper_cube(triangulation);
5493 * GridGenerator::hyper_rectangle(triangulation, Point<dim>(0.0,0.0), Point<dim>(1.0,1.0), true);
5496 * triangulation.refine_global (n_refinement);
5509 * for Reconstruction of MATERIAL FIELDS
5512 * min_h = GridTools::minimal_cell_diameter(triangulation)/std::sqrt(dim)/degree;
5513 * eps=1*min_h; //For reconstruction of density in Navier Stokes
5514 * sharpness=sharpness_integer*min_h; //adjust value of sharpness (for init cond of phi)
5523 * time_step = cfl*min_h/umax;
5527 * //////////////////
5529 * //////////////////
5532 * LevelSetSolver<dim> level_set (degree_LS,degree_U,
5536 * TRANSPORT_TIME_INTEGRATION,
5538 * mpi_communicator);
5542 * ///////////////////
5544 * ///////////////////
5547 * initial_condition();
5550 * level_set.initial_condition(locally_relevant_solution_phi,
5551 * locally_relevant_solution_u,locally_relevant_solution_v);
5553 * level_set.initial_condition(locally_relevant_solution_phi,
5554 * locally_relevant_solution_u,locally_relevant_solution_v,locally_relevant_solution_w);
5558 * /////////////////////////////
5559 * BOUNDARY CONDITIONS FOR PHI
5560 * /////////////////////////////
5563 * get_boundary_values_phi(boundary_values_id_phi,boundary_values_phi);
5564 * level_set.set_boundary_conditions(boundary_values_id_phi,boundary_values_phi);
5568 * OUTPUT DATA REGARDING TIME STEPPING AND MESH
5571 * int dofs_LS = dof_handler_LS.n_dofs();
5572 * pcout << "Cfl: " << cfl << std::endl;
5573 * pcout << " Number of active cells: "
5574 * << triangulation.n_global_active_cells() << std::endl
5575 * << " Number of degrees of freedom: " << std::endl
5576 * << " LS: " << dofs_LS << std::endl;
5583 * timestep_number=0;
5585 * while (time<final_time)
5587 * timestep_number++;
5588 * if (time+time_step > final_time)
5590 * pcout << "FINAL TIME STEP... " << std::endl;
5591 * time_step = final_time-time;
5593 * pcout << "Time step " << timestep_number
5594 * << "\twith dt=" << time_step
5595 * << "\tat tn=" << time << std::endl;
5600 * GET VELOCITY // (NS or interpolate from a function) at current time tn
5604 * if (VARIABLE_VELOCITY)
5606 * get_interpolated_velocity();
5609 * SET VELOCITY TO LEVEL SET SOLVER
5612 * level_set.set_velocity(locally_relevant_solution_u,locally_relevant_solution_v);
5616 * ////////////////////////
5617 * GET LEVEL SET SOLUTION // (at tnp1)
5618 * ////////////////////////
5621 * level_set.nth_time_step();
5630 * time+=time_step; // time tnp1
5639 * if (get_output && time-(output_number)*output_time>=0)
5641 * level_set.get_unp1(locally_relevant_solution_phi);
5645 * pcout << "FINAL TIME T=" << time << std::endl;
5648 * int main(int argc, char *argv[])
5652 * using namespace dealii;
5653 * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
5654 * PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
5655 * deallog.depth_console (0);
5657 * unsigned int degree = 1;
5658 * TestLevelSet<2> multiphase(degree, degree);
5663 * catch (std::exception &exc)
5665 * std::cerr << std::endl << std::endl
5666 * << "----------------------------------------------------"
5668 * std::cerr << "Exception on processing: " << std::endl
5669 * << exc.what() << std::endl
5670 * << "Aborting!" << std::endl
5671 * << "----------------------------------------------------"
5677 * std::cerr << std::endl << std::endl
5678 * << "----------------------------------------------------"
5680 * std::cerr << "Unknown exception!" << std::endl
5681 * << "Aborting!" << std::endl
5682 * << "----------------------------------------------------"
5693<a name="ann-TestNavierStokes.cc"></a>
5694<h1>Annotated version of TestNavierStokes.cc</h1>
5700 * /* -----------------------------------------------------------------------------
5702 * * SPDX-License-Identifier: LGPL-2.1-or-later
5703 * * Copyright (C) 2016 Manuel Quezada de Luna
5705 * * This file is part of the deal.II code gallery.
5707 * * -----------------------------------------------------------------------------
5710 * #include <deal.II/base/quadrature_lib.h>
5711 * #include <deal.II/base/function.h>
5712 * #include <deal.II/lac/affine_constraints.h>
5713 * #include <deal.II/lac/vector.h>
5714 * #include <deal.II/lac/full_matrix.h>
5715 * #include <deal.II/lac/solver_cg.h>
5716 * #include <deal.II/lac/petsc_sparse_matrix.h>
5717 * #include <deal.II/lac/petsc_vector.h>
5718 * #include <deal.II/lac/petsc_solver.h>
5719 * #include <deal.II/lac/petsc_precondition.h>
5720 * #include <deal.II/grid/grid_generator.h>
5721 * #include <deal.II/grid/tria_accessor.h>
5722 * #include <deal.II/grid/tria_iterator.h>
5723 * #include <deal.II/dofs/dof_handler.h>
5724 * #include <deal.II/dofs/dof_accessor.h>
5725 * #include <deal.II/dofs/dof_tools.h>
5726 * #include <deal.II/fe/fe_values.h>
5727 * #include <deal.II/fe/fe_q.h>
5728 * #include <deal.II/numerics/vector_tools.h>
5729 * #include <deal.II/numerics/data_out.h>
5730 * #include <deal.II/numerics/error_estimator.h>
5731 * #include <deal.II/base/utilities.h>
5732 * #include <deal.II/base/conditional_ostream.h>
5733 * #include <deal.II/base/index_set.h>
5734 * #include <deal.II/lac/sparsity_tools.h>
5735 * #include <deal.II/distributed/tria.h>
5736 * #include <deal.II/distributed/grid_refinement.h>
5737 * #include <deal.II/lac/petsc_vector.h>
5738 * #include <deal.II/base/convergence_table.h>
5739 * #include <deal.II/base/timer.h>
5740 * #include <deal.II/base/parameter_handler.h>
5741 * #include <fstream>
5742 * #include <iostream>
5743 * #include <deal.II/grid/grid_tools.h>
5744 * #include <deal.II/fe/mapping_q.h>
5745 * #include <deal.II/base/function.h>
5747 * using namespace dealii;
5749 * #include "utilities_test_NS.cc"
5750 * #include "NavierStokesSolver.cc"
5754 * ///////////////////////////////////////////////////
5755 * /////////////////// MAIN CLASS ////////////////////
5756 * ///////////////////////////////////////////////////
5759 * template <int dim>
5760 * class TestNavierStokes
5763 * TestNavierStokes (const unsigned int degree_LS,
5764 * const unsigned int degree_U);
5765 * ~TestNavierStokes ();
5769 * void get_boundary_values_U(double t);
5770 * void fix_pressure();
5771 * void output_results();
5772 * void process_solution(const unsigned int cycle);
5774 * void initial_condition();
5775 * void init_constraints();
5777 * PETScWrappers::MPI::Vector locally_relevant_solution_rho;
5778 * PETScWrappers::MPI::Vector locally_relevant_solution_u;
5779 * PETScWrappers::MPI::Vector locally_relevant_solution_v;
5780 * PETScWrappers::MPI::Vector locally_relevant_solution_w;
5781 * PETScWrappers::MPI::Vector locally_relevant_solution_p;
5782 * PETScWrappers::MPI::Vector completely_distributed_solution_rho;
5783 * PETScWrappers::MPI::Vector completely_distributed_solution_u;
5784 * PETScWrappers::MPI::Vector completely_distributed_solution_v;
5785 * PETScWrappers::MPI::Vector completely_distributed_solution_w;
5786 * PETScWrappers::MPI::Vector completely_distributed_solution_p;
5788 * std::vector<unsigned int> boundary_values_id_u;
5789 * std::vector<unsigned int> boundary_values_id_v;
5790 * std::vector<unsigned int> boundary_values_id_w;
5791 * std::vector<double> boundary_values_u;
5792 * std::vector<double> boundary_values_v;
5793 * std::vector<double> boundary_values_w;
5800 * MPI_Comm mpi_communicator;
5801 * parallel::distributed::Triangulation<dim> triangulation;
5804 * DoFHandler<dim> dof_handler_LS;
5806 * IndexSet locally_owned_dofs_LS;
5807 * IndexSet locally_relevant_dofs_LS;
5810 * DoFHandler<dim> dof_handler_U;
5812 * IndexSet locally_owned_dofs_U;
5813 * IndexSet locally_relevant_dofs_U;
5815 * DoFHandler<dim> dof_handler_P;
5817 * IndexSet locally_owned_dofs_P;
5818 * IndexSet locally_relevant_dofs_P;
5820 * AffineConstraints<double> constraints;
5824 * TimerOutput timer;
5832 * double final_time;
5833 * unsigned int timestep_number;
5838 * unsigned int n_cycles;
5839 * unsigned int n_refinement;
5840 * unsigned int output_number;
5841 * double output_time;
5849 * ConditionalOStream pcout;
5850 * ConvergenceTable convergence_table;
5855 * template <int dim>
5856 * TestNavierStokes<dim>::TestNavierStokes (const unsigned int degree_LS,
5857 * const unsigned int degree_U)
5859 * mpi_communicator (MPI_COMM_WORLD),
5860 * triangulation (mpi_communicator,
5861 * typename Triangulation<dim>::MeshSmoothing
5862 * (Triangulation<dim>::smoothing_on_refinement |
5863 * Triangulation<dim>::smoothing_on_coarsening)),
5864 * degree_LS(degree_LS),
5865 * dof_handler_LS (triangulation),
5866 * fe_LS (degree_LS),
5867 * degree_U(degree_U),
5868 * dof_handler_U (triangulation),
5870 * dof_handler_P (triangulation),
5871 * fe_P (degree_U-1), //TODO: change this to be degree_Q-1
5874 * timer(std::cout, TimerOutput::summary, TimerOutput::wall_times),
5877 * pcout (std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator)== 0))
5880 * template <int dim>
5881 * TestNavierStokes<dim>::~TestNavierStokes ()
5883 * dof_handler_LS.clear ();
5884 * dof_handler_U.clear ();
5885 * dof_handler_P.clear ();
5890 * /////////////////////////////////////
5891 * /////////////// SETUP ///////////////
5892 * /////////////////////////////////////
5895 * template <int dim>
5896 * void TestNavierStokes<dim>::setup()
5903 * dof_handler_LS.distribute_dofs (fe_LS);
5904 * locally_owned_dofs_LS = dof_handler_LS.locally_owned_dofs ();
5905 * locally_relevant_dofs_LS = DoFTools::extract_locally_relevant_dofs (dof_handler_LS);
5911 * dof_handler_U.distribute_dofs (fe_U);
5912 * locally_owned_dofs_U = dof_handler_U.locally_owned_dofs ();
5913 * locally_relevant_dofs_U = DoFTools::extract_locally_relevant_dofs (dof_handler_U);
5919 * dof_handler_P.distribute_dofs (fe_P);
5920 * locally_owned_dofs_P = dof_handler_P.locally_owned_dofs ();
5921 * locally_relevant_dofs_P = DoFTools::extract_locally_relevant_dofs (dof_handler_P);
5922 * init_constraints();
5925 * init vectors for rho
5928 * locally_relevant_solution_rho.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
5929 * locally_relevant_solution_rho = 0;
5930 * completely_distributed_solution_rho.reinit(locally_owned_dofs_LS,mpi_communicator);
5933 * init vectors for u
5936 * locally_relevant_solution_u.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
5937 * locally_relevant_solution_u = 0;
5938 * completely_distributed_solution_u.reinit(locally_owned_dofs_U,mpi_communicator);
5941 * init vectors for v
5944 * locally_relevant_solution_v.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
5945 * locally_relevant_solution_v = 0;
5946 * completely_distributed_solution_v.reinit(locally_owned_dofs_U,mpi_communicator);
5949 * init vectors for w
5952 * locally_relevant_solution_w.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
5953 * locally_relevant_solution_w = 0;
5954 * completely_distributed_solution_w.reinit(locally_owned_dofs_U,mpi_communicator);
5957 * init vectors for p
5960 * locally_relevant_solution_p.reinit(locally_owned_dofs_P,locally_relevant_dofs_P,mpi_communicator);
5961 * locally_relevant_solution_p = 0;
5962 * completely_distributed_solution_p.reinit(locally_owned_dofs_P,mpi_communicator);
5965 * template <int dim>
5966 * void TestNavierStokes<dim>::initial_condition()
5971 * Initial conditions
5972 * init condition for rho
5975 * completely_distributed_solution_rho = 0;
5976 * VectorTools::interpolate(dof_handler_LS,
5977 * RhoFunction<dim>(0),
5978 * completely_distributed_solution_rho);
5979 * constraints.distribute (completely_distributed_solution_rho);
5980 * locally_relevant_solution_rho = completely_distributed_solution_rho;
5983 * init condition for u
5986 * completely_distributed_solution_u = 0;
5987 * VectorTools::interpolate(dof_handler_U,
5988 * ExactSolution_and_BC_U<dim>(0,0),
5989 * completely_distributed_solution_u);
5990 * constraints.distribute (completely_distributed_solution_u);
5991 * locally_relevant_solution_u = completely_distributed_solution_u;
5994 * init condition for v
5997 * completely_distributed_solution_v = 0;
5998 * VectorTools::interpolate(dof_handler_U,
5999 * ExactSolution_and_BC_U<dim>(0,1),
6000 * completely_distributed_solution_v);
6001 * constraints.distribute (completely_distributed_solution_v);
6002 * locally_relevant_solution_v = completely_distributed_solution_v;
6005 * init condition for w
6010 * completely_distributed_solution_w = 0;
6011 * VectorTools::interpolate(dof_handler_U,
6012 * ExactSolution_and_BC_U<dim>(0,2),
6013 * completely_distributed_solution_w);
6014 * constraints.distribute (completely_distributed_solution_w);
6015 * locally_relevant_solution_w = completely_distributed_solution_w;
6019 * init condition for p
6022 * completely_distributed_solution_p = 0;
6023 * VectorTools::interpolate(dof_handler_P,
6024 * ExactSolution_p<dim>(0),
6025 * completely_distributed_solution_p);
6026 * constraints.distribute (completely_distributed_solution_p);
6027 * locally_relevant_solution_p = completely_distributed_solution_p;
6030 * template <int dim>
6031 * void TestNavierStokes<dim>::init_constraints()
6033 * constraints.clear ();
6034 * constraints.reinit (locally_relevant_dofs_LS);
6035 * DoFTools::make_hanging_node_constraints (dof_handler_LS, constraints);
6036 * constraints.close ();
6040 * void TestNavierStokes<dim>::fix_pressure()
6044 * fix the constant in the pressure
6047 * completely_distributed_solution_p = locally_relevant_solution_p;
6048 * double mean_value = VectorTools::compute_mean_value(dof_handler_P,
6050 * locally_relevant_solution_p,
6053 * completely_distributed_solution_p.add(-mean_value+std::sin(1)*(std::cos(time)-cos(1+time)));
6055 * completely_distributed_solution_p.add(-mean_value+8*std::pow(std::sin(0.5),3)*std::sin(1.5+time));
6056 * locally_relevant_solution_p = completely_distributed_solution_p;
6059 * template <int dim>
6060 * void TestNavierStokes<dim>::output_results ()
6062 * DataOut<dim> data_out;
6063 * data_out.attach_dof_handler (dof_handler_U);
6064 * data_out.add_data_vector (locally_relevant_solution_u, "u");
6065 * data_out.add_data_vector (locally_relevant_solution_v, "v");
6066 * if (dim==3) data_out.add_data_vector (locally_relevant_solution_w, "w");
6068 * Vector<float> subdomain (triangulation.n_active_cells());
6069 * for (unsigned int i=0; i<subdomain.size(); ++i)
6070 * subdomain(i) = triangulation.locally_owned_subdomain();
6071 * data_out.add_data_vector (subdomain, "subdomain");
6073 * data_out.build_patches ();
6075 * const std::string filename = ("solution-" +
6076 * Utilities::int_to_string (output_number, 3) +
6078 * Utilities::int_to_string
6079 * (triangulation.locally_owned_subdomain(), 4));
6080 * std::ofstream output ((filename + ".vtu").c_str());
6081 * data_out.write_vtu (output);
6083 * if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
6085 * std::vector<std::string> filenames;
6086 * for (unsigned int i=0;
6087 * i<Utilities::MPI::n_mpi_processes(mpi_communicator);
6089 * filenames.push_back ("solution-" +
6090 * Utilities::int_to_string (output_number, 3) +
6092 * Utilities::int_to_string (i, 4) +
6095 * std::ofstream master_output ((filename + ".pvtu").c_str());
6096 * data_out.write_pvtu_record (master_output, filenames);
6101 * template <int dim>
6102 * void TestNavierStokes<dim>::process_solution(const unsigned int cycle)
6104 * Vector<double> difference_per_cell (triangulation.n_active_cells());
6110 * VectorTools::integrate_difference (dof_handler_U,
6111 * locally_relevant_solution_u,
6112 * ExactSolution_and_BC_U<dim>(time,0),
6113 * difference_per_cell,
6114 * QGauss<dim>(degree_U+1),
6115 * VectorTools::L2_norm);
6116 * double u_L2_error = difference_per_cell.l2_norm();
6118 * std::sqrt(Utilities::MPI::sum(u_L2_error * u_L2_error, mpi_communicator));
6119 * VectorTools::integrate_difference (dof_handler_U,
6120 * locally_relevant_solution_u,
6121 * ExactSolution_and_BC_U<dim>(time,0),
6122 * difference_per_cell,
6123 * QGauss<dim>(degree_U+1),
6124 * VectorTools::H1_norm);
6125 * double u_H1_error = difference_per_cell.l2_norm();
6127 * std::sqrt(Utilities::MPI::sum(u_H1_error * u_H1_error, mpi_communicator));
6133 * VectorTools::integrate_difference (dof_handler_U,
6134 * locally_relevant_solution_v,
6135 * ExactSolution_and_BC_U<dim>(time,1),
6136 * difference_per_cell,
6137 * QGauss<dim>(degree_U+1),
6138 * VectorTools::L2_norm);
6139 * double v_L2_error = difference_per_cell.l2_norm();
6141 * std::sqrt(Utilities::MPI::sum(v_L2_error * v_L2_error,
6142 * mpi_communicator));
6143 * VectorTools::integrate_difference (dof_handler_U,
6144 * locally_relevant_solution_v,
6145 * ExactSolution_and_BC_U<dim>(time,1),
6146 * difference_per_cell,
6147 * QGauss<dim>(degree_U+1),
6148 * VectorTools::H1_norm);
6149 * double v_H1_error = difference_per_cell.l2_norm();
6151 * std::sqrt(Utilities::MPI::sum(v_H1_error *
6152 * v_H1_error, mpi_communicator));
6158 * double w_L2_error = 0;
6159 * double w_H1_error = 0;
6162 * VectorTools::integrate_difference (dof_handler_U,
6163 * locally_relevant_solution_w,
6164 * ExactSolution_and_BC_U<dim>(time,2),
6165 * difference_per_cell,
6166 * QGauss<dim>(degree_U+1),
6167 * VectorTools::L2_norm);
6168 * w_L2_error = difference_per_cell.l2_norm();
6170 * std::sqrt(Utilities::MPI::sum(w_L2_error * w_L2_error,
6171 * mpi_communicator));
6172 * VectorTools::integrate_difference (dof_handler_U,
6173 * locally_relevant_solution_w,
6174 * ExactSolution_and_BC_U<dim>(time,2),
6175 * difference_per_cell,
6176 * QGauss<dim>(degree_U+1),
6177 * VectorTools::H1_norm);
6178 * w_H1_error = difference_per_cell.l2_norm();
6180 * std::sqrt(Utilities::MPI::sum(w_H1_error *
6181 * w_H1_error, mpi_communicator));
6188 * VectorTools::integrate_difference (dof_handler_P,
6189 * locally_relevant_solution_p,
6190 * ExactSolution_p<dim>(time),
6191 * difference_per_cell,
6192 * QGauss<dim>(degree_U+1),
6193 * VectorTools::L2_norm);
6194 * double p_L2_error = difference_per_cell.l2_norm();
6196 * std::sqrt(Utilities::MPI::sum(p_L2_error * p_L2_error,
6197 * mpi_communicator));
6198 * VectorTools::integrate_difference (dof_handler_P,
6199 * locally_relevant_solution_p,
6200 * ExactSolution_p<dim>(time),
6201 * difference_per_cell,
6202 * QGauss<dim>(degree_U+1),
6203 * VectorTools::H1_norm);
6204 * double p_H1_error = difference_per_cell.l2_norm();
6206 * std::sqrt(Utilities::MPI::sum(p_H1_error * p_H1_error,
6207 * mpi_communicator));
6209 * const unsigned int n_active_cells=triangulation.n_active_cells();
6210 * const unsigned int n_dofs_U=dof_handler_U.n_dofs();
6211 * const unsigned int n_dofs_P=dof_handler_P.n_dofs();
6213 * convergence_table.add_value("cycle", cycle);
6214 * convergence_table.add_value("cells", n_active_cells);
6215 * convergence_table.add_value("dofs_U", n_dofs_U);
6216 * convergence_table.add_value("dofs_P", n_dofs_P);
6217 * convergence_table.add_value("dt", time_step);
6218 * convergence_table.add_value("u L2", u_L2_error);
6219 * convergence_table.add_value("u H1", u_H1_error);
6220 * convergence_table.add_value("v L2", v_L2_error);
6221 * convergence_table.add_value("v H1", v_H1_error);
6224 * convergence_table.add_value("w L2", w_L2_error);
6225 * convergence_table.add_value("w H1", w_H1_error);
6227 * convergence_table.add_value("p L2", p_L2_error);
6228 * convergence_table.add_value("p H1", p_H1_error);
6231 * template <int dim>
6232 * void TestNavierStokes<dim>::get_boundary_values_U(double t)
6234 * std::map<unsigned int, double> map_boundary_values_u;
6235 * std::map<unsigned int, double> map_boundary_values_v;
6237 * VectorTools::interpolate_boundary_values (dof_handler_U,0,ExactSolution_and_BC_U<dim>(t,0),map_boundary_values_u);
6238 * VectorTools::interpolate_boundary_values (dof_handler_U,0,ExactSolution_and_BC_U<dim>(t,1),map_boundary_values_v);
6240 * boundary_values_id_u.resize(map_boundary_values_u.size());
6241 * boundary_values_id_v.resize(map_boundary_values_v.size());
6242 * boundary_values_u.resize(map_boundary_values_u.size());
6243 * boundary_values_v.resize(map_boundary_values_v.size());
6244 * std::map<unsigned int,double>::const_iterator boundary_value_u =map_boundary_values_u.begin();
6245 * std::map<unsigned int,double>::const_iterator boundary_value_v =map_boundary_values_v.begin();
6248 * std::map<unsigned int, double> map_boundary_values_w;
6249 * VectorTools::interpolate_boundary_values (dof_handler_U,0,ExactSolution_and_BC_U<dim>(t,2),map_boundary_values_w);
6250 * boundary_values_id_w.resize(map_boundary_values_w.size());
6251 * boundary_values_w.resize(map_boundary_values_w.size());
6252 * std::map<unsigned int,double>::const_iterator boundary_value_w =map_boundary_values_w.begin();
6253 * for (int i=0; boundary_value_w !=map_boundary_values_w.end(); ++boundary_value_w, ++i)
6255 * boundary_values_id_w[i]=boundary_value_w->first;
6256 * boundary_values_w[i]=boundary_value_w->second;
6259 * for (int i=0; boundary_value_u !=map_boundary_values_u.end(); ++boundary_value_u, ++i)
6261 * boundary_values_id_u[i]=boundary_value_u->first;
6262 * boundary_values_u[i]=boundary_value_u->second;
6264 * for (int i=0; boundary_value_v !=map_boundary_values_v.end(); ++boundary_value_v, ++i)
6266 * boundary_values_id_v[i]=boundary_value_v->first;
6267 * boundary_values_v[i]=boundary_value_v->second;
6271 * template <int dim>
6272 * void TestNavierStokes<dim>::run()
6274 * if (Utilities::MPI::this_mpi_process(mpi_communicator)== 0)
6276 * std::cout << "***** CONVERGENCE TEST FOR NS *****" << std::endl;
6277 * std::cout << "DEGREE LS: " << degree_LS << std::endl;
6278 * std::cout << "DEGREE U: " << degree_U << std::endl;
6282 * PARAMETERS FOR THE NAVIER STOKES PROBLEM
6289 * ForceTerms<dim> force_function;
6290 * RhoFunction<dim> rho_function;
6291 * NuFunction<dim> nu_function;
6295 * bool get_output = false;
6296 * bool get_error = true;
6299 * for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
6303 * GridGenerator::hyper_cube (triangulation);
6304 * triangulation.refine_global (n_refinement);
6306 * initial_condition();
6310 * triangulation.refine_global(1);
6312 * initial_condition();
6322 * NavierStokesSolver<dim> navier_stokes (degree_LS,
6330 * mpi_communicator);
6333 * set INITIAL CONDITION within TRANSPORT PROBLEM
6337 * navier_stokes.initial_condition(locally_relevant_solution_rho,
6338 * locally_relevant_solution_u,
6339 * locally_relevant_solution_v,
6340 * locally_relevant_solution_p);
6342 * navier_stokes.initial_condition(locally_relevant_solution_rho,
6343 * locally_relevant_solution_u,
6344 * locally_relevant_solution_v,
6345 * locally_relevant_solution_w,
6346 * locally_relevant_solution_p);
6348 * pcout << "Cycle " << cycle << ':
' << std::endl;
6349 * pcout << " Cycle " << cycle
6350 * << " Number of active cells: "
6351 * << triangulation.n_global_active_cells() << std::endl
6352 * << " Number of degrees of freedom (velocity): "
6353 * << dof_handler_U.n_dofs() << std::endl
6354 * << " min h=" << GridTools::minimal_cell_diameter(triangulation)/std::sqrt(2)/degree_U
6362 * timestep_number=0;
6364 * double time_step_backup=time_step;
6365 * while (time<final_time)
6367 * timestep_number++;
6375 * if (time+time_step > final_time-1E-10)
6377 * pcout << "FINAL TIME STEP..." << std::endl;
6378 * time_step_backup=time_step;
6379 * time_step=final_time-time;
6381 * pcout << "Time step " << timestep_number
6382 * << "\twith dt=" << time_step
6383 * << "\tat tn=" << time
6392 * force_function.set_time(time+time_step);
6395 * /////////////////////////////
6396 * DENSITY AND VISCOSITY FIELD
6397 * /////////////////////////////
6400 * rho_function.set_time(time+time_step);
6401 * nu_function.set_time(time+time_step);
6404 * /////////////////////
6405 * BOUNDARY CONDITIONS
6406 * /////////////////////
6409 * get_boundary_values_U(time+time_step);
6410 * if (dim==2) navier_stokes.set_boundary_conditions(boundary_values_id_u, boundary_values_id_v,
6411 * boundary_values_u, boundary_values_v);
6412 * else navier_stokes.set_boundary_conditions(boundary_values_id_u,
6413 * boundary_values_id_v,
6414 * boundary_values_id_w,
6415 * boundary_values_u, boundary_values_v, boundary_values_w);
6423 * navier_stokes.nth_time_step();
6425 * navier_stokes.get_velocity(locally_relevant_solution_u,locally_relevant_solution_v);
6427 * navier_stokes.get_velocity(locally_relevant_solution_u,
6428 * locally_relevant_solution_v,
6429 * locally_relevant_solution_w);
6430 * navier_stokes.get_pressure(locally_relevant_solution_p);
6457 * if (get_output && time-(output_number)*output_time>=1E-10)
6460 * pcout << "FINAL TIME: " << time << std::endl;
6461 * time_step=time_step_backup;
6463 * process_solution(cycle);
6467 * convergence_table.set_precision("u L2", 2);
6468 * convergence_table.set_precision("u H1", 2);
6469 * convergence_table.set_scientific("u L2",true);
6470 * convergence_table.set_scientific("u H1",true);
6472 * convergence_table.set_precision("v L2", 2);
6473 * convergence_table.set_precision("v H1", 2);
6474 * convergence_table.set_scientific("v L2",true);
6475 * convergence_table.set_scientific("v H1",true);
6479 * convergence_table.set_precision("w L2", 2);
6480 * convergence_table.set_precision("w H1", 2);
6481 * convergence_table.set_scientific("w L2",true);
6482 * convergence_table.set_scientific("w H1",true);
6485 * convergence_table.set_precision("p L2", 2);
6486 * convergence_table.set_precision("p H1", 2);
6487 * convergence_table.set_scientific("p L2",true);
6488 * convergence_table.set_scientific("p H1",true);
6490 * convergence_table.set_tex_format("cells","r");
6491 * convergence_table.set_tex_format("dofs_U","r");
6492 * convergence_table.set_tex_format("dofs_P","r");
6493 * convergence_table.set_tex_format("dt","r");
6495 * if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
6497 * std::cout << std::endl;
6498 * convergence_table.write_text(std::cout);
6504 * int main(int argc, char *argv[])
6508 * using namespace dealii;
6509 * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
6510 * PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
6511 * deallog.depth_console (0);
6514 * unsigned int degree_LS = 1;
6515 * unsigned int degree_U = 2;
6516 * TestNavierStokes<2> test_navier_stokes(degree_LS, degree_U);
6517 * test_navier_stokes.run();
6524 * catch (std::exception &exc)
6526 * std::cerr << std::endl << std::endl
6527 * << "----------------------------------------------------"
6529 * std::cerr << "Exception on processing: " << std::endl
6530 * << exc.what() << std::endl
6531 * << "Aborting!" << std::endl
6532 * << "----------------------------------------------------"
6539 * std::cerr << std::endl << std::endl
6540 * << "----------------------------------------------------"
6542 * std::cerr << "Unknown exception!" << std::endl
6543 * << "Aborting!" << std::endl
6544 * << "----------------------------------------------------"
6554<a name="ann-clean.sh"></a>
6555<h1>Annotated version of clean.sh</h1>
6557rm -rf CMakeFiles CMakeCache.txt Makefile cmake_install.cmake *~
6558rm -f MultiPhase TestLevelSet TestNavierStokes
6565<a name="ann-utilities.cc"></a>
6566<h1>Annotated version of utilities.cc</h1>
6572 * /* -----------------------------------------------------------------------------
6574 * * SPDX-License-Identifier: LGPL-2.1-or-later
6575 * * Copyright (C) 2016 Manuel Quezada de Luna
6577 * * This file is part of the deal.II code gallery.
6579 * * -----------------------------------------------------------------------------
6584 * /////////////////////////////////////////////////
6585 * ////////////////// INITIAL PHI //////////////////
6586 * /////////////////////////////////////////////////
6589 * template <int dim>
6590 * class InitialPhi : public Function <dim>
6593 * InitialPhi (unsigned int PROBLEM, double sharpness=0.005) : Function<dim>(),
6594 * sharpness(sharpness),
6595 * PROBLEM(PROBLEM) {}
6596 * virtual double value (const Point<dim> &p, const unsigned int component=0) const override;
6598 * unsigned int PROBLEM;
6600 * template <int dim>
6601 * double InitialPhi<dim>::value (const Point<dim> &p,
6602 * const unsigned int) const
6606 * double pi=numbers::PI;
6608 * if (PROBLEM==FILLING_TANK)
6609 * return 0.5*(-std::tanh((y-0.3)/sharpness)*std::tanh((y-0.35)/sharpness)+1)
6610 * *(-std::tanh((x-0.02)/sharpness)+1)-1;
6611 * else if (PROBLEM==BREAKING_DAM)
6612 * return 0.5*(-std::tanh((x-0.35)/sharpness)*std::tanh((x-0.65)/sharpness)+1)
6613 * *(1-std::tanh((y-0.35)/sharpness))-1;
6614 * else if (PROBLEM==FALLING_DROP)
6619 * double r = std::sqrt(std::pow(x-x0,2)+std::pow(y-y0,2));
6620 * return 1-(std::tanh((r-r0)/sharpness)+std::tanh((y-0.3)/sharpness));
6622 * else if (PROBLEM==SMALL_WAVE_PERTURBATION)
6624 * double wave = 0.1*std::sin(pi*x)+0.25;
6625 * return -std::tanh((y-wave)/sharpness);
6629 * std::cout << "Error in type of PROBLEM" << std::endl;
6636 * ///////////////////////////////////////////////////
6637 * ////////////////// FORCE TERMS ///// //////////////
6638 * ///////////////////////////////////////////////////
6641 * template <int dim>
6642 * class ForceTerms : public Functions::ConstantFunction <dim>
6645 * ForceTerms (const std::vector<double> values) : Functions::ConstantFunction<dim>(values) {}
6650 * /////////////////////////////////////////////////
6651 * ////////////////// BOUNDARY PHI /////////////////
6652 * /////////////////////////////////////////////////
6655 * template <int dim>
6656 * class BoundaryPhi : public Functions::ConstantFunction <dim>
6659 * BoundaryPhi (const double value, const unsigned int n_components=1) : Functions::ConstantFunction<dim>(value,n_components) {}
6664 * //////////////////////////////////////////////////////
6665 * ////////////////// BOUNDARY VELOCITY /////////////////
6666 * //////////////////////////////////////////////////////
6669 * template <int dim>
6670 * class BoundaryU : public Function <dim>
6673 * BoundaryU (unsigned int PROBLEM, double t=0) : Function<dim>(), PROBLEM(PROBLEM) {this->set_time(t);}
6674 * virtual double value (const Point<dim> &p, const unsigned int component=0) const override;
6677 * template <int dim>
6678 * double BoundaryU<dim>::value (const Point<dim> &p, const unsigned int) const
6682 * //////////////////
6684 * //////////////////
6685 * boundary for filling the tank (inlet)
6691 * if (PROBLEM==FILLING_TANK)
6693 * if (x==0 && y>=0.3 && y<=0.35)
6700 * std::cout << "Error in PROBLEM definition" << std::endl;
6705 * template <int dim>
6706 * class BoundaryV : public Function <dim>
6709 * BoundaryV (unsigned int PROBLEM, double t=0) : Function<dim>(), PROBLEM(PROBLEM) {this->set_time(t);}
6710 * virtual double value (const Point<dim> &p, const unsigned int component=0) const override;
6711 * unsigned int PROBLEM;
6713 * template <int dim>
6714 * double BoundaryV<dim>::value (const Point<dim> &p, const unsigned int) const
6718 * boundary for filling the tank (outlet)
6723 * double return_value = 0;
6725 * if (PROBLEM==FILLING_TANK)
6727 * if (y==0.4 && x>=0.3 && x<=0.35)
6728 * return_value = 0.25;
6730 * return return_value;
6735 * ///////////////////////////////////////////////////
6736 * ///////////////// POST-PROCESSING /////////////////
6737 * ///////////////////////////////////////////////////
6740 * template <int dim>
6741 * class Postprocessor : public DataPostprocessorScalar <dim>
6744 * Postprocessor(double eps, double rho_air, double rho_fluid)
6746 * DataPostprocessorScalar<dim>("Density",update_values)
6749 * this->rho_air=rho_air;
6750 * this->rho_fluid=rho_fluid;
6755 * evaluate_scalar_field (const DataPostprocessorInputs::Scalar<dim> &input_data,
6756 * std::vector<Vector<double> > &computed_quantities) const override;
6764 * template <int dim>
6766 * Postprocessor<dim>::
6767 * evaluate_scalar_field (const DataPostprocessorInputs::Scalar<dim> &input_data,
6768 * std::vector<Vector<double> > &computed_quantities) const
6770 * const unsigned int n_quadrature_points = input_data.solution_values.size();
6771 * for (unsigned int q=0; q<n_quadrature_points; ++q)
6775 * double phi_value=input_data.solution_values[q];
6776 * if (phi_value > eps)
6778 * else if (phi_value < -eps)
6782 * rho_value = rho_fluid*(1+H)/2. + rho_air*(1-H)/2.;
6783 * computed_quantities[q] = rho_value;
6790<a name="ann-utilities_test_LS.cc"></a>
6791<h1>Annotated version of utilities_test_LS.cc</h1>
6797 * /* -----------------------------------------------------------------------------
6799 * * SPDX-License-Identifier: LGPL-2.1-or-later
6800 * * Copyright (C) 2016 Manuel Quezada de Luna
6802 * * This file is part of the deal.II code gallery.
6804 * * -----------------------------------------------------------------------------
6809 * ///////////////////////////////////////////////////
6810 * ////////////////// INITIAL PHI //////////////////
6811 * ///////////////////////////////////////////////////
6814 * template <int dim>
6815 * class InitialPhi : public Function <dim>
6818 * InitialPhi (unsigned int PROBLEM, double sharpness=0.005) : Function<dim>(),
6819 * sharpness(sharpness),
6820 * PROBLEM(PROBLEM) {}
6821 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6823 * unsigned int PROBLEM;
6825 * template <int dim>
6826 * double InitialPhi<dim>::value (const Point<dim> &p,
6827 * const unsigned int) const
6831 * double return_value = -1.;
6833 * if (PROBLEM==CIRCULAR_ROTATION)
6838 * double r = std::sqrt(std::pow(x-x0,2)+std::pow(y-y0,2));
6839 * return_value = -std::tanh((r-r0)/sharpness);
6841 * else // (PROBLEM==DIAGONAL_ADVECTION)
6848 * r = std::sqrt(std::pow(x-x0,2)+std::pow(y-y0,2));
6853 * r = std::sqrt(std::pow(x-x0,2)+std::pow(y-y0,2)+std::pow(z-z0,2));
6855 * return_value = -std::tanh((r-r0)/sharpness);
6857 * return return_value;
6862 * /////////////////////////////////////////////////
6863 * ////////////////// BOUNDARY PHI /////////////////
6864 * /////////////////////////////////////////////////
6867 * template <int dim>
6868 * class BoundaryPhi : public Function <dim>
6871 * BoundaryPhi (double t=0)
6874 * {this->set_time(t);}
6875 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6878 * template <int dim>
6879 * double BoundaryPhi<dim>::value (const Point<dim> &, const unsigned int) const
6886 * ///////////////////////////////////////////////////
6887 * ////////////////// EXACT VELOCITY /////////////////
6888 * ///////////////////////////////////////////////////
6891 * template <int dim>
6892 * class ExactU : public Function <dim>
6895 * ExactU (unsigned int PROBLEM, double time=0) : Function<dim>(), PROBLEM(PROBLEM), time(time) {}
6896 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6897 * void set_time(double time) {this->time=time;};
6902 * template <int dim>
6903 * double ExactU<dim>::value (const Point<dim> &p, const unsigned int) const
6905 * if (PROBLEM==CIRCULAR_ROTATION)
6906 * return -2*numbers::PI*(p[1]-0.5);
6907 * else // (PROBLEM==DIAGONAL_ADVECTION)
6911 * template <int dim>
6912 * class ExactV : public Function <dim>
6915 * ExactV (unsigned int PROBLEM, double time=0) : Function<dim>(), PROBLEM(PROBLEM), time(time) {}
6916 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6917 * void set_time(double time) {this->time=time;};
6918 * unsigned int PROBLEM;
6922 * template <int dim>
6923 * double ExactV<dim>::value (const Point<dim> &p, const unsigned int) const
6925 * if (PROBLEM==CIRCULAR_ROTATION)
6926 * return 2*numbers::PI*(p[0]-0.5);
6927 * else // (PROBLEM==DIAGONAL_ADVECTION)
6931 * template <int dim>
6932 * class ExactW : public Function <dim>
6935 * ExactW (unsigned int PROBLEM, double time=0) : Function<dim>(), PROBLEM(PROBLEM), time(time) {}
6936 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6937 * void set_time(double time) {this->time=time;};
6938 * unsigned int PROBLEM;
6942 * template <int dim>
6943 * double ExactW<dim>::value (const Point<dim> &, const unsigned int) const
6947 * PROBLEM = 3D_DIAGONAL_ADVECTION
6956<a name="ann-utilities_test_NS.cc"></a>
6957<h1>Annotated version of utilities_test_NS.cc</h1>
6963 * /* -----------------------------------------------------------------------------
6965 * * SPDX-License-Identifier: LGPL-2.1-or-later
6966 * * Copyright (C) 2016 Manuel Quezada de Luna
6968 * * This file is part of the deal.II code gallery.
6970 * * -----------------------------------------------------------------------------
6975 * ///////////////////////////////////////////////////
6976 * ////////// EXACT SOLUTION RHO TO TEST NS //////////
6977 * ///////////////////////////////////////////////////
6980 * template <int dim>
6981 * class RhoFunction : public Function <dim>
6984 * RhoFunction (double t=0) : Function<dim>() {this->set_time(t);}
6985 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6987 * template <int dim>
6988 * double RhoFunction<dim>::value (const Point<dim> &p,
6989 * const unsigned int) const
6991 * double t = this->get_time();
6992 * double return_value = 0;
6994 * return_value = std::pow(std::sin(p[0]+p[1]+t),2)+1;
6996 * return_value = std::pow(std::sin(p[0]+p[1]+p[2]+t),2)+1;
6997 * return return_value;
7000 * template <int dim>
7001 * class NuFunction : public Function <dim>
7004 * NuFunction (double t=0) : Function<dim>() {this->set_time(t);}
7005 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
7007 * template <int dim>
7008 * double NuFunction<dim>::value (const Point<dim> &, const unsigned int) const
7015 * //////////////////////////////////////////////////////////////
7016 * ///////////////// EXACT SOLUTION U to TEST NS ////////////////
7017 * //////////////////////////////////////////////////////////////
7020 * template <int dim>
7021 * class ExactSolution_and_BC_U : public Function <dim>
7024 * ExactSolution_and_BC_U (double t=0, int field=0)
7029 * this->set_time(t);
7031 * virtual double value (const Point<dim> &p, const unsigned int component=1) const;
7032 * virtual Tensor<1,dim> gradient (const Point<dim> &p, const unsigned int component=1) const;
7033 * virtual void set_field(int field) {this->field=field;}
7035 * unsigned int type_simulation;
7037 * template <int dim>
7038 * double ExactSolution_and_BC_U<dim>::value (const Point<dim> &p,
7039 * const unsigned int) const
7041 * double t = this->get_time();
7042 * double return_value = 0;
7043 * double Pi = numbers::PI;
7050 * return_value = std::sin(x)*std::sin(y+t);
7052 * return_value = std::cos(x)*std::cos(y+t);
7057 * return_value = std::cos(t)*std::cos(Pi*y)*std::cos(Pi*z)*std::sin(Pi*x);
7058 * else if (field == 1)
7059 * return_value = std::cos(t)*std::cos(Pi*x)*std::cos(Pi*z)*std::sin(Pi*y);
7061 * return_value = -2*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::sin(Pi*z);
7063 * return return_value;
7065 * template <int dim>
7066 * Tensor<1,dim> ExactSolution_and_BC_U<dim>::gradient (const Point<dim> &p,
7067 * const unsigned int) const
7071 * THIS IS USED JUST FOR TESTING NS
7074 * Tensor<1,dim> return_value;
7075 * double t = this->get_time();
7076 * double Pi = numbers::PI;
7083 * return_value[0] = std::cos(x)*std::sin(y+t);
7084 * return_value[1] = std::sin(x)*std::cos(y+t);
7088 * return_value[0] = -std::sin(x)*std::cos(y+t);
7089 * return_value[1] = -std::cos(x)*std::sin(y+t);
7096 * return_value[0] = Pi*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::cos(Pi*z);
7097 * return_value[1] = -(Pi*std::cos(t)*std::cos(Pi*z)*std::sin(Pi*x)*std::sin(Pi*y));
7098 * return_value[2] = -(Pi*std::cos(t)*std::cos(Pi*y)*std::sin(Pi*x)*std::sin(Pi*z));
7100 * else if (field == 1)
7102 * return_value[0] = -(Pi*std::cos(t)*std::cos(Pi*z)*std::sin(Pi*x)*std::sin(Pi*y));
7103 * return_value[1] = Pi*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::cos(Pi*z);
7104 * return_value[2] = -(Pi*std::cos(t)*std::cos(Pi*x)*std::sin(Pi*y)*std::sin(Pi*z));
7108 * return_value[0] = 2*Pi*std::cos(t)*std::cos(Pi*y)*std::sin(Pi*x)*std::sin(Pi*z);
7109 * return_value[1] = 2*Pi*std::cos(t)*std::cos(Pi*x)*std::sin(Pi*y)*std::sin(Pi*z);
7110 * return_value[2] = -2*Pi*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::cos(Pi*z);
7113 * return return_value;
7118 * ///////////////////////////////////////////////////
7119 * ///////// EXACT SOLUTION FOR p TO TEST NS /////////
7120 * ///////////////////////////////////////////////////
7123 * template <int dim>
7124 * class ExactSolution_p : public Function <dim>
7127 * ExactSolution_p (double t=0) : Function<dim>() {this->set_time(t);}
7128 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
7129 * virtual Tensor<1,dim> gradient (const Point<dim> &p, const unsigned int component = 0) const;
7132 * template <int dim>
7133 * double ExactSolution_p<dim>::value (const Point<dim> &p, const unsigned int) const
7135 * double t = this->get_time();
7136 * double return_value = 0;
7138 * return_value = std::cos(p[0])*std::sin(p[1]+t);
7140 * return_value = std::sin(p[0]+p[1]+p[2]+t);
7141 * return return_value;
7144 * template <int dim>
7145 * Tensor<1,dim> ExactSolution_p<dim>::gradient (const Point<dim> &p, const unsigned int) const
7147 * Tensor<1,dim> return_value;
7148 * double t = this->get_time();
7151 * return_value[0] = -std::sin(p[0])*std::sin(p[1]+t);
7152 * return_value[1] = std::cos(p[0])*std::cos(p[1]+t);
7156 * return_value[0] = std::cos(t+p[0]+p[1]+p[2]);
7157 * return_value[1] = std::cos(t+p[0]+p[1]+p[2]);
7158 * return_value[2] = std::cos(t+p[0]+p[1]+p[2]);
7160 * return return_value;
7165 * //////////////////////////////////////////////////////////////
7166 * ////////////////// FORCE TERMS to TEST NS ////////////////////
7167 * //////////////////////////////////////////////////////////////
7170 * template <int dim>
7171 * class ForceTerms : public Function <dim>
7174 * ForceTerms (double t=0)
7178 * this->set_time(t);
7181 * virtual void vector_value (const Point<dim> &p, Vector<double> &values) const;
7185 * template <int dim>
7186 * void ForceTerms<dim>::vector_value (const Point<dim> &p, Vector<double> &values) const
7191 * double t = this->get_time();
7192 * double Pi = numbers::PI;
7201 * values[0] = std::cos(t+y)*std::sin(x)*(1+std::pow(std::sin(t+x+y),2)) // time derivative
7202 * +2*nu*std::sin(x)*std::sin(t+y) // viscosity
7203 * +std::cos(x)*std::sin(x)*(1+std::pow(std::sin(t+x+y),2)) // non-linearity
7204 * -std::sin(x)*std::sin(y+t); // pressure
7210 * values[1] = -(std::cos(x)*std::sin(t+y)*(1+std::pow(std::sin(t+x+y),2))) // time derivative
7211 * +2*nu*std::cos(x)*std::cos(t+y) // viscosity
7212 * -(std::sin(2*(t+y))*(1+std::pow(std::sin(t+x+y),2)))/2. // non-linearity
7213 * +std::cos(x)*std::cos(y+t); // pressure
7224 * -(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.
7225 * +3*std::pow(Pi,2)*std::cos(t)*std::cos(Pi*y)*std::cos(Pi*z)*std::sin(Pi*x) //viscosity
7226 * -(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
7227 * +std::cos(t+x+y+z); // pressure
7229 * -(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
7230 * +3*std::pow(Pi,2)*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*z)*std::sin(Pi*y) //viscosity
7231 * -(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
7232 * +std::cos(t+x+y+z); // pressure
7234 * 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
7235 * -6*std::pow(Pi,2)*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::sin(Pi*z) //viscosity
7236 * -(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
7237 * +std::cos(t+x+y+z); // pressure
unsigned int depth_console(const unsigned int n)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
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)
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
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)
@ matrix
Contents is actually a matrix.
constexpr types::blas_int zero
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(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)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
std::string compress(const std::string &input)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void abort(const ExceptionBase &exc) noexcept
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation