246 * #include <deal.II/base/quadrature_lib.h>
247 * #include <deal.II/base/function.h>
248 * #include <deal.II/lac/affine_constraints.h>
249 * #include <deal.II/lac/vector.h>
250 * #include <deal.II/lac/full_matrix.h>
251 * #include <deal.II/lac/solver_cg.h>
252 * #include <deal.II/lac/petsc_sparse_matrix.h>
253 * #include <deal.II/lac/petsc_sparse_matrix.h>
254 * #include <deal.II/lac/petsc_vector.h>
255 * #include <deal.II/lac/petsc_solver.h>
256 * #include <deal.II/lac/petsc_precondition.h>
257 * #include <deal.II/grid/grid_generator.h>
258 * #include <deal.II/grid/tria_accessor.h>
259 * #include <deal.II/grid/tria_iterator.h>
260 * #include <deal.II/dofs/dof_handler.h>
261 * #include <deal.II/dofs/dof_accessor.h>
262 * #include <deal.II/dofs/dof_tools.h>
263 * #include <deal.II/fe/fe_values.h>
264 * #include <deal.II/fe/fe_q.h>
265 * #include <deal.II/numerics/vector_tools.h>
266 * #include <deal.II/numerics/data_out.h>
267 * #include <deal.II/numerics/error_estimator.h>
268 * #include <deal.II/base/utilities.h>
269 * #include <deal.II/base/conditional_ostream.h>
270 * #include <deal.II/base/index_set.h>
271 * #include <deal.II/lac/sparsity_tools.h>
272 * #include <deal.II/distributed/tria.h>
273 * #include <deal.II/distributed/grid_refinement.h>
274 * #include <deal.II/lac/vector.h>
275 * #include <deal.II/base/convergence_table.h>
276 * #include <deal.II/base/timer.h>
277 * #include <deal.II/base/parameter_handler.h>
278 * #include <deal.II/grid/grid_tools.h>
279 * #include <deal.II/fe/mapping_q.h>
284 * #include <iostream>
295 * #define CHECK_MAX_PRINCIPLE 0
299 * LOG FOR LEVEL SET FROM -1 to 1
303 * #define ENTROPY_GRAD(phi,phix) 2*phi*phix*((1-phi*phi>=0) ? -1 : 1)/(
std::
abs(1-phi*phi)+1
E-14)
310 * This is a solver for the transpor solver.
311 * We assume the velocity is divergence free
312 * and solve the equation in conservation form.
314 * ---------- NOTATION ----------
316 * We use notation popular in the literature of conservation laws.
317 * For this reason the solution is denoted as u, unm1, unp1, etc.
318 * and the velocity is treated as vx, vy and vz.
322 * class LevelSetSolver
342 * BOUNDARY CONDITIONS
346 *
void set_boundary_conditions(std::vector<types::global_dof_index> &boundary_values_id_u,
347 * std::vector<double> boundary_values_u);
375 *
void nth_time_step();
385 * LevelSetSolver (
const unsigned int degree_LS,
386 *
const unsigned int degree_U,
387 *
const double time_step,
390 *
const bool verbose,
391 * std::string ALGORITHM,
392 *
const unsigned int TIME_INTEGRATION,
401 * ASSEMBLE MASS (and other) MATRICES
405 *
void assemble_ML();
407 *
void assemble_MC();
411 * LOW ORDER METHOD (DiJ Viscosity)
415 *
void assemble_C_Matrix();
425 *
void assemble_EntRes_Matrix();
429 * FOR MAXIMUM PRINCIPLE
455 * std::string algorithm);
458 * std::string algorithm);
466 *
void get_sparsity_pattern();
467 *
void get_map_from_Q1_to_Q2();
470 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
473 *
void save_old_solution();
474 *
void save_old_vel_solution();
483 *
const std::vector<types::global_dof_index> &indices,
484 * std::vector<PetscScalar> &values);
486 *
const std::vector<types::global_dof_index> &indices,
487 * std::map<types::global_dof_index, types::global_dof_index> &map_from_Q1_to_Q2,
488 * std::vector<PetscScalar> &values);
494 * FINITE ELEMENT SPACE
502 *
IndexSet locally_relevant_dofs_LS;
512 * OPERATORS times SOLUTION VECTOR
525 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> MC_preconditioner;
532 * std::vector<types::global_dof_index> boundary_values_id_u;
533 * std::vector<double> boundary_values_u;
540 * FOR FIRST ORDER VISCOSITY
547 * FOR ENTROPY VISCOSITY
553 * FOR FCT (flux and limited flux)
582 * NON-GHOSTED VECTORS
620 *
double solver_tolerance;
621 *
double entropy_normalization_factor;
629 * std::string ALGORITHM;
630 *
unsigned int TIME_INTEGRATION;
634 * std::map<types::global_dof_index, types::global_dof_index> map_from_Q1_to_Q2;
635 * std::map<types::global_dof_index, std::vector<types::global_dof_index> > sparsity_pattern;
639 * LevelSetSolver<dim>::LevelSetSolver (
const unsigned int degree_LS,
640 *
const unsigned int degree_U,
641 *
const double time_step,
644 *
const bool verbose,
645 * std::string ALGORITHM,
646 *
const unsigned int TIME_INTEGRATION,
650 * mpi_communicator (mpi_communicator),
651 * degree_LS(degree_LS),
654 * degree_U(degree_U),
657 * time_step(time_step),
661 * ALGORITHM(ALGORITHM),
662 * TIME_INTEGRATION(TIME_INTEGRATION),
665 * pcout <<
"********** LEVEL SET SETUP **********" << std::endl;
670 * LevelSetSolver<dim>::~LevelSetSolver ()
672 * dof_handler_LS.clear ();
673 * dof_handler_U.clear ();
692 * this->locally_relevant_solution_vx = locally_relevant_solution_vx;
693 * this->locally_relevant_solution_vy = locally_relevant_solution_vy;
696 * initialize old vectors with current solution,
this just happens the
first time
700 * locally_relevant_solution_vx_old = locally_relevant_solution_vx;
701 * locally_relevant_solution_vy_old = locally_relevant_solution_vy;
711 * this->locally_relevant_solution_vx = locally_relevant_solution_vx;
712 * this->locally_relevant_solution_vy = locally_relevant_solution_vy;
713 * this->locally_relevant_solution_vz = locally_relevant_solution_vz;
716 * initialize old vectors with current solution,
this just happens the
first time
720 * locally_relevant_solution_vx_old = locally_relevant_solution_vx;
721 * locally_relevant_solution_vy_old = locally_relevant_solution_vy;
722 * locally_relevant_solution_vz_old = locally_relevant_solution_vz;
733 *
void LevelSetSolver<dim>::set_boundary_conditions(std::vector<types::global_dof_index> &boundary_values_id_u,
734 * std::vector<double> boundary_values_u)
736 * this->boundary_values_id_u = boundary_values_id_u;
737 * this->boundary_values_u = boundary_values_u;
756 * save_old_vel_solution();
762 * this->locally_relevant_solution_vx=locally_relevant_solution_vx;
763 * this->locally_relevant_solution_vy=locally_relevant_solution_vy;
776 * save_old_vel_solution();
782 * this->locally_relevant_solution_vx=locally_relevant_solution_vx;
783 * this->locally_relevant_solution_vy=locally_relevant_solution_vy;
784 * this->locally_relevant_solution_vz=locally_relevant_solution_vz;
799 * -------------------------------------------------------------------------------
800 * ------------------------------ COMPUTE SOLUTIONS ------------------------------
801 * -------------------------------------------------------------------------------
805 *
void LevelSetSolver<dim>::nth_time_step()
807 * assemble_EntRes_Matrix();
813 *
if (TIME_INTEGRATION==FORWARD_EULER)
814 * compute_solution(unp1,un,ALGORITHM);
816 * compute_solution_SSP33(unp1,un,ALGORITHM);
819 * BOUNDARY CONDITIONS
822 * unp1.set(boundary_values_id_u,boundary_values_u);
826 * CHECK MAXIMUM PRINCIPLE
829 *
if (CHECK_MAX_PRINCIPLE)
831 * compute_bounds(un);
832 * check_max_principle(unp1);
836 * pcout <<
"*********************************************************************... "
837 * << unp1.min() <<
", " << unp1.max() << std::endl;
840 * save_old_solution();
845 * --------------------------------------------------------------------
846 * ------------------------------ SETUP ------------------------------
847 * --------------------------------------------------------------------
851 *
void LevelSetSolver<dim>::setup()
853 * solver_tolerance=1
E-6;
854 * degree_MAX =
std::max(degree_LS,degree_U);
858 * SETUP FOR DOF HANDLERS
863 * dof_handler_LS.distribute_dofs (fe_LS);
864 * locally_owned_dofs_LS = dof_handler_LS.locally_owned_dofs ();
871 * dof_handler_U.distribute_dofs (fe_U);
872 * locally_owned_dofs_U = dof_handler_U.locally_owned_dofs ();
881 * constraints.clear ();
882 * constraints.reinit (locally_relevant_dofs_LS);
884 * constraints.close ();
888 * NON-GHOSTED VECTORS
892 * MPP_uL_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
893 * NMPP_uH_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
894 * RHS.reinit(locally_owned_dofs_LS,mpi_communicator);
895 * uStage1_nonGhosted.reinit (locally_owned_dofs_LS,mpi_communicator);
896 * uStage2_nonGhosted.reinit (locally_owned_dofs_LS,mpi_communicator);
897 * unp1.reinit (locally_owned_dofs_LS,mpi_communicator);
898 * MPP_uH_solution.reinit (locally_owned_dofs_LS,mpi_communicator);
901 * vectors
for lumped mass
matrix
904 * ML_vector.reinit(locally_owned_dofs_LS,mpi_communicator);
905 * inverse_ML_vector.reinit(locally_owned_dofs_LS,mpi_communicator);
906 * ones_vector.reinit(locally_owned_dofs_LS,mpi_communicator);
910 * operators times solution
913 * K_times_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
914 * DL_times_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
915 * DH_times_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
921 * R_pos_vector_nonGhosted.reinit (locally_owned_dofs_LS,mpi_communicator);
922 * R_neg_vector_nonGhosted.reinit (locally_owned_dofs_LS,mpi_communicator);
923 * umin_vector.reinit (locally_owned_dofs_LS,mpi_communicator);
924 * umax_vector.reinit (locally_owned_dofs_LS,mpi_communicator);
928 * GHOSTED VECTORS (used within some assemble process)
932 * uStage1.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
933 * uStage2.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
934 * unm1.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
935 * un.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
936 * MPP_uL_solution_ghosted.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
937 * MPP_uLkp1_solution_ghosted.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
938 * NMPP_uH_solution_ghosted.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
941 * init vectors
for vx
944 * locally_relevant_solution_vx.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
945 * locally_relevant_solution_vx_old.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
948 * init vectors
for vy
951 * locally_relevant_solution_vy.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
952 * locally_relevant_solution_vy_old.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
955 * init vectors
for vz
958 * locally_relevant_solution_vz.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
959 * locally_relevant_solution_vz_old.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
965 * R_pos_vector.reinit(locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
966 * R_neg_vector.reinit(locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
978 * dof_handler_LS.locally_owned_dofs(),
980 * locally_relevant_dofs_LS);
981 * MC_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
982 * dof_handler_LS.locally_owned_dofs(),
983 * dsp, mpi_communicator);
984 * Cx_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
985 * dof_handler_LS.locally_owned_dofs(),
986 * dsp, mpi_communicator);
987 * CTx_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
988 * dof_handler_LS.locally_owned_dofs(),
989 * dsp, mpi_communicator);
990 * Cy_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
991 * dof_handler_LS.locally_owned_dofs(),
992 * dsp, mpi_communicator);
993 * CTy_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
994 * dof_handler_LS.locally_owned_dofs(),
995 * dsp, mpi_communicator);
998 * Cz_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
999 * dof_handler_LS.locally_owned_dofs(),
1000 * dsp, mpi_communicator);
1001 * CTz_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1002 * dof_handler_LS.locally_owned_dofs(),
1003 * dsp, mpi_communicator);
1005 * dLij_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1006 * dof_handler_LS.locally_owned_dofs(),
1007 * dsp, mpi_communicator);
1008 * EntRes_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1009 * dof_handler_LS.locally_owned_dofs(),
1010 * dsp, mpi_communicator);
1011 * SuppSize_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1012 * dof_handler_LS.locally_owned_dofs(),
1013 * dsp, mpi_communicator);
1014 * dCij_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1015 * dof_handler_LS.locally_owned_dofs(),
1016 * dsp, mpi_communicator);
1017 * A_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1018 * dof_handler_LS.locally_owned_dofs(),
1019 * dsp, mpi_communicator);
1020 * LxA_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1021 * dof_handler_LS.locally_owned_dofs(),
1022 * dsp, mpi_communicator);
1023 * Akp1_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1024 * dof_handler_LS.locally_owned_dofs(),
1025 * dsp, mpi_communicator);
1026 * LxAkp1_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1027 * dof_handler_LS.locally_owned_dofs(),
1028 * dsp, mpi_communicator);
1031 * COMPUTE MASS MATRICES (AND OTHERS) FOR FIRST TIME STEP
1037 * assemble_C_Matrix();
1040 * get mat
for DOFs between Q1 and Q2
1043 * get_map_from_Q1_to_Q2();
1044 * get_sparsity_pattern();
1049 * ----------------------------------------------------------------------------
1050 * ------------------------------ MASS MATRICES ------------------------------
1051 * ----------------------------------------------------------------------------
1055 *
void LevelSetSolver<dim>::assemble_ML()
1059 *
const QGauss<dim> quadrature_formula(degree_MAX+1);
1065 *
const unsigned int dofs_per_cell = fe_LS.dofs_per_cell;
1066 *
const unsigned int n_q_points = quadrature_formula.size();
1069 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1072 * cell_LS = dof_handler_LS.begin_active(),
1073 * endc_LS = dof_handler_LS.end();
1075 *
for (; cell_LS!=endc_LS; ++cell_LS)
1076 *
if (cell_LS->is_locally_owned())
1079 * fe_values_LS.reinit (cell_LS);
1080 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
1082 *
const double JxW = fe_values_LS.JxW(q_point);
1083 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1084 * cell_ML (i) += fe_values_LS.shape_value(i,q_point)*JxW;
1091 * cell_LS->get_dof_indices (local_dof_indices);
1092 * constraints.distribute_local_to_global (cell_ML,local_dof_indices,ML_vector);
1103 *
void LevelSetSolver<dim>::invert_ML()
1107 *
loop on locally owned i-DOFs (rows)
1111 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
1113 *
int gi = *idofs_iter;
1114 * inverse_ML_vector(gi) = 1./ML_vector(gi);
1120 *
void LevelSetSolver<dim>::assemble_MC()
1124 *
const QGauss<dim> quadrature_formula(degree_MAX+1);
1130 *
const unsigned int dofs_per_cell = fe_LS.dofs_per_cell;
1131 *
const unsigned int n_q_points = quadrature_formula.size();
1134 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1135 * std::vector<double> shape_values(dofs_per_cell);
1138 * cell_LS = dof_handler_LS.begin_active(),
1139 * endc_LS = dof_handler_LS.end();
1141 *
for (; cell_LS!=endc_LS; ++cell_LS)
1142 *
if (cell_LS->is_locally_owned())
1145 * fe_values_LS.reinit (cell_LS);
1146 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
1148 *
const double JxW = fe_values_LS.JxW(q_point);
1149 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1150 * shape_values[i] = fe_values_LS.shape_value(i,q_point);
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 * ---------------------------------------------------------------------------------------
1176 * ------------------------------
LO METHOD (Dij Viscosity) ------------------------------
1177 * ---------------------------------------------------------------------------------------
1180 *
template <
int dim>
1181 *
void LevelSetSolver<dim>::assemble_C_Matrix ()
1190 *
const QGauss<dim> quadrature_formula(degree_MAX+1);
1196 *
const unsigned int dofs_per_cell_LS = fe_LS.dofs_per_cell;
1197 *
const unsigned int n_q_points = quadrature_formula.size();
1206 * std::vector<Tensor<1, dim> > shape_grads_LS(dofs_per_cell_LS);
1207 * std::vector<double> shape_values_LS(dofs_per_cell_LS);
1209 * std::vector<types::global_dof_index> local_dof_indices_LS (dofs_per_cell_LS);
1212 * cell_LS = dof_handler_LS.begin_active();
1213 * endc_LS = dof_handler_LS.end();
1215 *
for (; cell_LS!=endc_LS; ++cell_LS)
1216 *
if (cell_LS->is_locally_owned())
1228 * fe_values_LS.reinit (cell_LS);
1229 * cell_LS->get_dof_indices (local_dof_indices_LS);
1231 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
1233 *
const double JxW = fe_values_LS.JxW(q_point);
1234 *
for (
unsigned int i=0; i<dofs_per_cell_LS; ++i)
1236 * shape_values_LS[i] = fe_values_LS.shape_value(i,q_point);
1237 * shape_grads_LS [i] = fe_values_LS.shape_grad (i,q_point);
1240 *
for (
unsigned int i=0; i<dofs_per_cell_LS; ++i)
1241 *
for (
unsigned int j=0; j < dofs_per_cell_LS; ++j)
1243 * cell_Cij_x(i,j) += (shape_grads_LS[j][0])*shape_values_LS[i]*JxW;
1244 * cell_Cij_y(i,j) += (shape_grads_LS[j][1])*shape_values_LS[i]*JxW;
1245 * cell_Cji_x(i,j) += (shape_grads_LS[i][0])*shape_values_LS[j]*JxW;
1246 * cell_Cji_y(i,j) += (shape_grads_LS[i][1])*shape_values_LS[j]*JxW;
1249 * cell_Cij_z(i,j) += (shape_grads_LS[j][2])*shape_values_LS[i]*JxW;
1250 * cell_Cji_z(i,j) += (shape_grads_LS[i][2])*shape_values_LS[j]*JxW;
1259 * constraints.distribute_local_to_global(cell_Cij_x,local_dof_indices_LS,Cx_matrix);
1260 * constraints.distribute_local_to_global(cell_Cji_x,local_dof_indices_LS,CTx_matrix);
1261 * constraints.distribute_local_to_global(cell_Cij_y,local_dof_indices_LS,Cy_matrix);
1262 * constraints.distribute_local_to_global(cell_Cji_y,local_dof_indices_LS,CTy_matrix);
1265 * constraints.distribute_local_to_global(cell_Cij_z,local_dof_indices_LS,Cz_matrix);
1266 * constraints.distribute_local_to_global(cell_Cji_z,local_dof_indices_LS,CTz_matrix);
1288 * K_times_solution = 0;
1290 *
const QGauss<dim> quadrature_formula(degree_MAX+1);
1299 *
const unsigned int dofs_per_cell = fe_LS.dofs_per_cell;
1300 *
const unsigned int n_q_points = quadrature_formula.size();
1304 * std::vector<Tensor<1,dim> > un_grads (n_q_points);
1305 * std::vector<double> old_vx_values (n_q_points);
1306 * std::vector<double> old_vy_values (n_q_points);
1307 * std::vector<double> old_vz_values (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);
1322 * cell_LS = dof_handler_LS.begin_active(),
1323 * endc_LS = dof_handler_LS.end();
1325 * cell_U = dof_handler_U.begin_active();
1328 *
for (; cell_LS!=endc_LS; ++cell_U, ++cell_LS)
1329 *
if (cell_LS->is_locally_owned())
1331 * cell_K_times_solution=0;
1333 * fe_values_LS.reinit (cell_LS);
1334 * cell_LS->get_dof_indices (indices_LS);
1335 * fe_values_LS.get_function_gradients(solution,un_grads);
1337 * fe_values_U.reinit (cell_U);
1338 * fe_values_U.get_function_values(locally_relevant_solution_vx,old_vx_values);
1339 * fe_values_U.get_function_values(locally_relevant_solution_vy,old_vy_values);
1340 *
if (dim==3) fe_values_U.get_function_values(locally_relevant_solution_vz,old_vz_values);
1344 * compute cell_K_times_solution
1347 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
1349 * v[0] = old_vx_values[q_point];
1350 * v[1] = old_vy_values[q_point];
1351 *
if (dim==3) v[2] = old_vz_values[q_point];
1353 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1354 * cell_K_times_solution(i) += (v*un_grads[q_point])
1355 * *fe_values_LS.shape_value(i,q_point)*fe_values_LS.JxW(q_point);
1362 * constraints.distribute_local_to_global (cell_K_times_solution, indices_LS, K_times_solution);
1367 *
template <
int dim>
1368 *
void LevelSetSolver<dim>::assemble_K_DL_DH_times_vector
1373 * K_times_solution=0;
1376 * DL_times_solution=0;
1377 * DH_times_solution=0;
1381 * PetscInt ncolumns;
1382 *
const PetscInt *gj;
1383 *
const PetscScalar *Cxi, *Cyi, *Czi, *CTxi, *CTyi, *CTzi;
1384 *
const PetscScalar *EntResi, *SuppSizei, *MCi;
1391 *
loop on locally owned i-DOFs (rows)
1396 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
1398 * PetscInt gi = *idofs_iter;
1401 *
double ith_K_times_solution = 0;
1405 * read velocity of i-th DOF
1408 * vi[0] = locally_relevant_solution_vx(map_from_Q1_to_Q2[gi]);
1409 * vi[1] = locally_relevant_solution_vy(map_from_Q1_to_Q2[gi]);
1410 *
if (dim==3) vi[2] = locally_relevant_solution_vz(map_from_Q1_to_Q2[gi]);
1411 * solni = solution(gi);
1415 * get i-th row of
C matrices
1418 * MatGetRow(Cx_matrix,gi,&ncolumns,&gj,&Cxi);
1419 * MatGetRow(Cy_matrix,gi,&ncolumns,&gj,&Cyi);
1420 * MatGetRow(CTx_matrix,gi,&ncolumns,&gj,&CTxi);
1421 * MatGetRow(CTy_matrix,gi,&ncolumns,&gj,&CTyi);
1424 * MatGetRow(Cz_matrix,gi,&ncolumns,&gj,&Czi);
1425 * MatGetRow(CTz_matrix,gi,&ncolumns,&gj,&CTzi);
1427 * MatGetRow(EntRes_matrix,gi,&ncolumns,&gj,&EntResi);
1428 * MatGetRow(SuppSize_matrix,gi,&ncolumns,&gj,&SuppSizei);
1429 * MatGetRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
1433 * get vector
values for column indices
1436 *
const std::vector<types::global_dof_index> gj_indices (gj,gj+ncolumns);
1437 * std::vector<double> soln(ncolumns);
1438 * std::vector<double> vx(ncolumns);
1439 * std::vector<double> vy(ncolumns);
1440 * std::vector<double> vz(ncolumns);
1441 * get_vector_values(solution,gj_indices,soln);
1442 * get_vector_values(locally_relevant_solution_vx,gj_indices,map_from_Q1_to_Q2,vx);
1443 * get_vector_values(locally_relevant_solution_vy,gj_indices,map_from_Q1_to_Q2,vy);
1445 * get_vector_values(locally_relevant_solution_vz,gj_indices,map_from_Q1_to_Q2,vz);
1449 * Array
for i-th row of matrices
1452 * std::vector<double> dLi(ncolumns), dCi(ncolumns);
1453 *
double dLii = 0, dCii = 0;
1456 *
loop on sparsity pattern of i-th DOF
1459 *
for (
int j =0; j < ncolumns; ++j)
1476 * ith_K_times_solution += soln[j]*(vj*
C);
1483 * low order dissipative
matrix
1490 * high order dissipative
matrix (entropy viscosity)
1494 * cE*
std::abs(EntResi[j])/(entropy_normalization_factor*MCi[j]/SuppSizei[j]));
1497 * high order compression
matrix
1501 * dCi[j] = dEij*
std::max(1-Compij,0.0);
1507 * save
K times solution vector
1508 * K_times_solution(gi)=ith_K_times_solution;
1509 * save i-th row of matrices on global matrices
1512 * MatSetValuesRow(dLij_matrix,gi,&dLi[0]);
1513 * dLij_matrix.set(gi,gi,dLii);
1514 * MatSetValuesRow(dCij_matrix,gi,&dCi[0]);
1515 * dCij_matrix.set(gi,gi,dCii);
1519 * Restore matrices after reading rows
1522 * MatRestoreRow(Cx_matrix,gi,&ncolumns,&gj,&Cxi);
1523 * MatRestoreRow(Cy_matrix,gi,&ncolumns,&gj,&Cyi);
1524 * MatRestoreRow(CTx_matrix,gi,&ncolumns,&gj,&CTxi);
1525 * MatRestoreRow(CTy_matrix,gi,&ncolumns,&gj,&CTyi);
1528 * MatRestoreRow(Cz_matrix,gi,&ncolumns,&gj,&Czi);
1529 * MatRestoreRow(CTz_matrix,gi,&ncolumns,&gj,&CTzi);
1531 * MatRestoreRow(EntRes_matrix,gi,&ncolumns,&gj,&EntResi);
1532 * MatRestoreRow(SuppSize_matrix,gi,&ncolumns,&gj,&SuppSizei);
1533 * MatRestoreRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
1545 * get matrices times vector
1548 * dLij_matrix.vmult(DL_times_solution,solution);
1549 * dCij_matrix.vmult(DH_times_solution,solution);
1554 * --------------------------------------------------------------------------------------
1555 * ------------------------------ ENTROPY VISCOSITY ------------------------------
1556 * --------------------------------------------------------------------------------------
1559 *
template <
int dim>
1560 *
void LevelSetSolver<dim>::assemble_EntRes_Matrix ()
1563 * entropy_normalization_factor=0;
1564 * SuppSize_matrix=0;
1566 *
const QGauss<dim> quadrature_formula(degree_MAX+1);
1576 *
const unsigned int dofs_per_cell_LS = fe_LS.dofs_per_cell;
1577 *
const unsigned int n_q_points = quadrature_formula.size();
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);
1594 * std::vector<Tensor<1, dim> > shape_grads_LS(dofs_per_cell_LS);
1595 * std::vector<double> shape_values_LS(dofs_per_cell_LS);
1597 * std::vector<types::global_dof_index> local_dof_indices_LS (dofs_per_cell_LS);
1600 * cell_LS = dof_handler_LS.begin_active();
1601 * endc_LS = dof_handler_LS.end();
1603 * cell_U = dof_handler_U.begin_active();
1606 *
double max_entropy=-1E10, min_entropy=1E10;
1607 *
double cell_max_entropy, cell_min_entropy;
1608 *
double cell_entropy_mass, entropy_mass=0;
1609 *
double cell_volume_double,
volume=0;
1611 *
for (; cell_LS!=endc_LS; ++cell_LS, ++cell_U)
1612 *
if (cell_LS->is_locally_owned())
1614 * cell_entropy_mass = 0;
1615 * cell_volume_double = 0;
1616 * cell_max_entropy = -1E10;
1617 * cell_min_entropy = 1E10;
1623 * get solutions at quadrature points
1626 * fe_values_LS.reinit(cell_LS);
1627 * cell_LS->get_dof_indices (local_dof_indices_LS);
1628 * fe_values_LS.get_function_values(un,uqn);
1629 * fe_values_LS.get_function_values(unm1,uqnm1);
1630 * fe_values_LS.get_function_gradients(un,guqn);
1631 * fe_values_LS.get_function_gradients(unm1,guqnm1);
1633 * fe_values_U.reinit(cell_U);
1634 * fe_values_U.get_function_values(locally_relevant_solution_vx,vxqn);
1635 * fe_values_U.get_function_values(locally_relevant_solution_vy,vyqn);
1636 *
if (dim==3) fe_values_U.get_function_values(locally_relevant_solution_vz,vzqn);
1637 * fe_values_U.get_function_values(locally_relevant_solution_vx_old,vxqnm1);
1638 * fe_values_U.get_function_values(locally_relevant_solution_vy_old,vyqnm1);
1639 *
if (dim==3) fe_values_U.get_function_values(locally_relevant_solution_vz_old,vzqnm1);
1641 *
for (
unsigned int q=0; q<n_q_points; ++q)
1643 * Rk = 1./time_step*(ENTROPY(uqn[q])-ENTROPY(uqnm1[q]))
1644 * +(vxqn[q]*ENTROPY_GRAD(uqn[q],guqn[q][0])+vyqn[q]*ENTROPY_GRAD(uqn[q],guqn[q][1]))/2.
1645 * +(vxqnm1[q]*ENTROPY_GRAD(uqnm1[q],guqnm1[q][0])+vyqnm1[q]*ENTROPY_GRAD(uqnm1[q],guqnm1[q][1]))/2.;
1647 * Rk += 0.5*(vzqn[q]*ENTROPY_GRAD(uqn[q],guqn[q][2])+vzqnm1[q]*ENTROPY_GRAD(uqnm1[q],guqnm1[q][2]));
1649 *
const double JxW = fe_values_LS.JxW(q);
1650 *
for (
unsigned int i=0; i<dofs_per_cell_LS; ++i)
1652 * shape_values_LS[i] = fe_values_LS.shape_value(i,q);
1653 * shape_grads_LS [i] = fe_values_LS.shape_grad (i,q);
1656 *
for (
unsigned int i=0; i<dofs_per_cell_LS; ++i)
1657 *
for (
unsigned int j=0; j < dofs_per_cell_LS; ++j)
1659 * cell_EntRes (i,j) += Rk*shape_values_LS[i]*shape_values_LS[j]*JxW;
1660 * cell_volume (i,j) += JxW;
1662 * cell_entropy_mass += ENTROPY(uqn[q])*JxW;
1663 * cell_volume_double += JxW;
1665 * cell_min_entropy =
std::min(cell_min_entropy,ENTROPY(uqn[q]));
1666 * cell_max_entropy =
std::max(cell_max_entropy,ENTROPY(uqn[q]));
1668 * entropy_mass += cell_entropy_mass;
1669 *
volume += cell_volume_double;
1671 * min_entropy =
std::min(min_entropy,cell_min_entropy);
1672 * max_entropy =
std::max(max_entropy,cell_max_entropy);
1678 * constraints.distribute_local_to_global(cell_EntRes,local_dof_indices_LS,EntRes_matrix);
1679 * constraints.distribute_local_to_global(cell_volume,local_dof_indices_LS,SuppSize_matrix);
1685 * ENTROPY NORM FACTOR
1697 * ------------------------------------------------------------------------------------
1698 * ------------------------------ TO CHECK MAX PRINCIPLE ------------------------------
1699 * ------------------------------------------------------------------------------------
1709 *
loop on locally owned i-DOFs (rows)
1713 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
1715 *
int gi = *idofs_iter;
1719 * get solution at DOFs on the sparsity pattern of i-th DOF
1722 * std::vector<types::global_dof_index> gj_indices = sparsity_pattern[gi];
1723 * std::vector<double> soln(gj_indices.size());
1724 * get_vector_values(un_solution,gj_indices,soln);
1727 * compute bounds, ith row of flux
matrix, P vectors
1730 *
double mini=1E10, maxi=-1E10;
1731 *
for (
unsigned int j =0; j < gj_indices.size(); ++j)
1741 * umin_vector(gi) = mini;
1742 * umax_vector(gi) = maxi;
1753 * compute
min and
max vectors
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);
1761 * cell_LS = dof_handler_LS.begin_active(),
1762 * endc_LS = dof_handler_LS.end();
1764 *
for (; cell_LS!=endc_LS; ++cell_LS)
1765 *
if (cell_LS->is_locally_owned() && !cell_LS->at_boundary())
1767 * cell_LS->get_dof_indices(local_dof_indices);
1768 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1769 *
if (locally_owned_dofs_LS.is_element(local_dof_indices[i]))
1771 *
double solni = unp1_solution(local_dof_indices[i]);
1772 *
if (solni - umin_vector(local_dof_indices[i]) < -tol || umax_vector(local_dof_indices[i]) - solni < -tol)
1774 * pcout <<
"MAX Principle violated" << std::endl;
1783 * -------------------------------------------------------------------------------
1784 * ------------------------------ COMPUTE SOLUTIONS ------------------------------
1785 * -------------------------------------------------------------------------------
1789 *
void LevelSetSolver<dim>::compute_MPP_uL_and_NMPP_uH
1796 * NON-GHOSTED VECTORS: MPP_uL_solution, NMPP_uH_solution
1797 * GHOSTED VECTORS: un_solution
1800 * MPP_uL_solution=un_solution;
1801 * NMPP_uH_solution=un_solution;
1807 * assemble_K_times_vector(un_solution);
1808 * assemble_K_DL_DH_times_vector(un_solution);
1812 * COMPUTE MPP u1 solution
1816 * MPP_uL_solution.scale(ML_vector);
1817 * MPP_uL_solution.add(-time_step,K_times_solution);
1818 * MPP_uL_solution.add(-time_step,DL_times_solution);
1819 * MPP_uL_solution.scale(inverse_ML_vector);
1823 * COMPUTE GALERKIN u2 solution
1827 * MC_matrix.vmult(RHS,un_solution);
1828 * RHS.add(-time_step,K_times_solution,-time_step,DH_times_solution);
1829 * solve(constraints,MC_matrix,MC_preconditioner,NMPP_uH_solution,RHS);
1832 *
template <
int dim>
1833 *
void LevelSetSolver<dim>::compute_MPP_uH
1839 * MPP_uH_solution=0;
1842 *
loop on locally owned i-DOFs (rows)
1847 * PetscInt ncolumns;
1848 *
const PetscInt *gj;
1849 *
const PetscScalar *MCi, *dLi, *dCi;
1850 *
double solni, mi, solLi, solHi;
1852 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
1854 *
int gi = *idofs_iter;
1857 * read vectors at i-th DOF
1860 * solni=solution(gi);
1861 * solHi=NMPP_uH_solution_ghosted(gi);
1862 * solLi=MPP_uL_solution_ghosted(gi);
1867 * get i-th row of matrices
1870 * MatGetRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
1871 * MatGetRow(dLij_matrix,gi,&ncolumns,&gj,&dLi);
1872 * MatGetRow(dCij_matrix,gi,&ncolumns,&gj,&dCi);
1876 * get vector
values for support of i-th DOF
1879 *
const std::vector<types::global_dof_index> gj_indices (gj,gj+ncolumns);
1880 * std::vector<double> soln(ncolumns);
1881 * std::vector<double> solH(ncolumns);
1882 * get_vector_values(solution,gj_indices,soln);
1883 * get_vector_values(NMPP_uH_solution_ghosted,gj_indices,solH);
1887 * Array
for i-th row of matrices
1890 * std::vector<double> Ai(ncolumns);
1893 * compute bounds, ith row of flux
matrix, P vectors
1896 *
double mini=1E10, maxi=-1E10;
1897 *
double Pposi=0 ,Pnegi=0;
1898 *
for (
int j =0; j < ncolumns; ++j)
1910 * i-th row of flux
matrix A
1913 * Ai[j] = (((gi==gj[j]) ? 1 : 0)*mi - MCi[j])*(solH[j]-soln[j] - (solHi-solni))
1914 * +time_step*(dLi[j]-dCi[j])*(soln[j]-solni);
1921 * Pposi += Ai[j]*((Ai[j] > 0) ? 1. : 0.);
1922 * Pnegi += Ai[j]*((Ai[j] < 0) ? 1. : 0.);
1926 * save i-th row of flux
matrix A
1929 * MatSetValuesRow(A_matrix,gi,&Ai[0]);
1936 *
double Qposi = mi*(maxi-solLi);
1937 *
double Qnegi = mi*(mini-solLi);
1944 * R_pos_vector_nonGhosted(gi) = ((Pposi==0) ? 1. :
std::
min(1.0,Qposi/Pposi));
1945 * R_neg_vector_nonGhosted(gi) = ((Pnegi==0) ? 1. :
std::
min(1.0,Qnegi/Pnegi));
1949 * Restore matrices after reading rows
1952 * MatRestoreRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
1953 * MatRestoreRow(dLij_matrix,gi,&ncolumns,&gj,&dLi);
1954 * MatRestoreRow(dCij_matrix,gi,&ncolumns,&gj,&dCi);
1971 * update ghost
values for R vectors
1974 * R_pos_vector = R_pos_vector_nonGhosted;
1975 * R_neg_vector = R_neg_vector_nonGhosted;
1979 * compute limiters. NOTE:
this is a different
loop due to need of i- and j-th entries of R vectors
1983 *
double Rposi, Rnegi;
1984 * idofs_iter=locally_owned_dofs_LS.begin();
1985 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
1987 *
int gi = *idofs_iter;
1988 * Rposi = R_pos_vector(gi);
1989 * Rnegi = R_neg_vector(gi);
1993 * get i-th row of A
matrix
1996 * MatGetRow(A_matrix,gi,&ncolumns,&gj,&Ai);
2000 * get vector
values for column indices
2003 *
const std::vector<types::global_dof_index> gj_indices (gj,gj+ncolumns);
2004 * std::vector<double> Rpos(ncolumns);
2005 * std::vector<double> Rneg(ncolumns);
2006 * get_vector_values(R_pos_vector,gj_indices,Rpos);
2007 * get_vector_values(R_neg_vector,gj_indices,Rneg);
2011 * Array
for i-th row of A_times_L
matrix
2014 * std::vector<double> LxAi(ncolumns);
2017 *
loop in sparsity pattern of i-th DOF
2020 *
for (
int j =0; j < ncolumns; ++j)
2021 * LxAi[j] = Ai[j] * ((Ai[j]>0) ?
std::min(Rposi,Rneg[j]) :
std::
min(Rnegi,Rpos[j]));
2025 * save i-th row of LxA
2028 * MatSetValuesRow(LxA_matrix,gi,&LxAi[0]);
2031 * restore A
matrix after reading it
2034 * MatRestoreRow(A_matrix,gi,&ncolumns,&gj,&Ai);
2037 * LxA_matrix.vmult(MPP_uH_solution,ones_vector);
2038 * MPP_uH_solution.scale(inverse_ML_vector);
2039 * MPP_uH_solution.add(1.0,MPP_uL_solution_ghosted);
2043 *
void LevelSetSolver<dim>::compute_MPP_uH_with_iterated_FCT
2049 * MPP_uH_solution=0;
2050 * compute_MPP_uH(MPP_uH_solution,MPP_uL_solution_ghosted,NMPP_uH_solution_ghosted,un_solution);
2054 * Akp1_matrix.copy_from(A_matrix);
2055 * LxAkp1_matrix.copy_from(LxA_matrix);
2059 *
loop in num of FCT iterations
2062 * PetscInt ncolumns;
2063 *
const PetscInt *gj;
2064 *
const PetscScalar *Akp1i;
2066 *
for (
int iter=0; iter<NUM_ITER; ++iter)
2068 * MPP_uLkp1_solution_ghosted = MPP_uH_solution;
2069 * Akp1_matrix.add(-1.0, LxAkp1_matrix);
2073 *
loop on locally owned i-DOFs (rows)
2077 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
2079 *
int gi = *idofs_iter;
2083 * read vectors at i-th DOF
2087 *
double solLi = MPP_uLkp1_solution_ghosted(gi);
2091 * get i-th row of matrices
2094 * MatGetRow(Akp1_matrix,gi,&ncolumns,&gj,&Akp1i);
2097 * get vector
values for support of i-th DOF
2100 *
const std::vector<types::global_dof_index> gj_indices (gj,gj+ncolumns);
2101 * std::vector<double> soln(ncolumns);
2102 * get_vector_values(un_solution,gj_indices,soln);
2106 * compute bounds, ith row of flux
matrix, P vectors
2109 *
double mini=1E10, maxi=-1E10;
2110 *
double Pposi=0 ,Pnegi=0;
2111 *
for (
int j =0; j < ncolumns; ++j)
2126 * Pposi += Akp1i[j]*((Akp1i[j] > 0) ? 1. : 0.);
2127 * Pnegi += Akp1i[j]*((Akp1i[j] < 0) ? 1. : 0.);
2134 *
double Qposi = mi*(maxi-solLi);
2135 *
double Qnegi = mi*(mini-solLi);
2142 * R_pos_vector_nonGhosted(gi) = ((Pposi==0) ? 1. :
std::
min(1.0,Qposi/Pposi));
2143 * R_neg_vector_nonGhosted(gi) = ((Pnegi==0) ? 1. :
std::
min(1.0,Qnegi/Pnegi));
2147 * Restore matrices after reading rows
2150 * MatRestoreRow(Akp1_matrix,gi,&ncolumns,&gj,&Akp1i);
2161 * update ghost
values for R vectors
2164 * R_pos_vector = R_pos_vector_nonGhosted;
2165 * R_neg_vector = R_neg_vector_nonGhosted;
2169 * compute limiters. NOTE:
this is a different
loop due to need of i- and j-th entries of R vectors
2172 *
double Rposi, Rnegi;
2173 * idofs_iter=locally_owned_dofs_LS.begin();
2174 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
2176 *
int gi = *idofs_iter;
2177 * Rposi = R_pos_vector(gi);
2178 * Rnegi = R_neg_vector(gi);
2182 * get i-th row of Akp1
matrix
2185 * MatGetRow(Akp1_matrix,gi,&ncolumns,&gj,&Akp1i);
2189 * get vector
values for column indices
2192 *
const std::vector<types::global_dof_index> gj_indices(gj,gj+ncolumns);
2193 * std::vector<double> Rpos(ncolumns);
2194 * std::vector<double> Rneg(ncolumns);
2195 * get_vector_values(R_pos_vector,gj_indices,Rpos);
2196 * get_vector_values(R_neg_vector,gj_indices,Rneg);
2200 * Array
for i-th row of LxAkp1
matrix
2203 * std::vector<double> LxAkp1i(ncolumns);
2204 *
for (
int j =0; j < ncolumns; ++j)
2205 * LxAkp1i[j] = Akp1i[j] * ((Akp1i[j]>0) ?
std::min(Rposi,Rneg[j]) :
std::
min(Rnegi,Rpos[j]));
2209 * save i-th row of LxA
2212 * MatSetValuesRow(LxAkp1_matrix,gi,&LxAkp1i[0]);
2215 * restore A
matrix after reading it
2218 * MatRestoreRow(Akp1_matrix,gi,&ncolumns,&gj,&Akp1i);
2221 * LxAkp1_matrix.vmult(MPP_uH_solution,ones_vector);
2222 * MPP_uH_solution.scale(inverse_ML_vector);
2223 * MPP_uH_solution.add(1.0,MPP_uLkp1_solution_ghosted);
2231 * std::string algorithm)
2236 * COMPUTE MPP LOW-ORDER SOLN and NMPP HIGH-ORDER SOLN
2239 * compute_MPP_uL_and_NMPP_uH(MPP_uL_solution,NMPP_uH_solution,un);
2241 *
if (algorithm.compare(
"MPP_u1")==0)
2242 * unp1=MPP_uL_solution;
2243 *
else if (algorithm.compare(
"NMPP_uH")==0)
2244 * unp1=NMPP_uH_solution;
2245 *
else if (algorithm.compare(
"MPP_uH")==0)
2247 * MPP_uL_solution_ghosted = MPP_uL_solution;
2248 * NMPP_uH_solution_ghosted=NMPP_uH_solution;
2249 * compute_MPP_uH_with_iterated_FCT(MPP_uH_solution,MPP_uL_solution_ghosted,NMPP_uH_solution_ghosted,un);
2250 * unp1=MPP_uH_solution;
2254 * pcout <<
"Error in algorithm" << std::endl;
2262 * std::string algorithm)
2266 * GHOSTED VECTORS: un
2267 * NON-GHOSTED VECTORS: unp1
2271 * uStage1=0., uStage2=0.;
2272 * uStage1_nonGhosted=0., uStage2_nonGhosted=0.;
2281 * compute_solution(uStage1_nonGhosted,un,algorithm);
2282 * uStage1=uStage1_nonGhosted;
2288 * u2=3/4*un+1/4*(u1-dt*RH*u1)
2291 * compute_solution(uStage2_nonGhosted,uStage1,algorithm);
2292 * uStage2_nonGhosted*=1./4;
2293 * uStage2_nonGhosted.add(3./4,un);
2294 * uStage2=uStage2_nonGhosted;
2300 * unp1=1/3*un+2/3*(u2-dt*RH*u2)
2303 * compute_solution(unp1,uStage2,algorithm);
2305 * unp1.add(1./3,un);
2310 * -----------------------------------------------------------------------
2311 * ------------------------------ UTILITIES ------------------------------
2312 * -----------------------------------------------------------------------
2316 *
void LevelSetSolver<dim>::get_sparsity_pattern()
2324 * PetscInt ncolumns;
2325 *
const PetscInt *gj;
2326 *
const PetscScalar *MCi;
2328 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
2330 * PetscInt gi = *idofs_iter;
2333 * get i-th row of mass
matrix (dummy, I just need the indices gj)
2336 * MatGetRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
2337 * sparsity_pattern[gi] = std::vector<types::global_dof_index>(gj,gj+ncolumns);
2338 * MatRestoreRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
2343 *
void LevelSetSolver<dim>::get_map_from_Q1_to_Q2()
2345 * map_from_Q1_to_Q2.clear();
2346 *
const unsigned int dofs_per_cell_LS = fe_LS.dofs_per_cell;
2347 * std::vector<types::global_dof_index> local_dof_indices_LS (dofs_per_cell_LS);
2348 *
const unsigned int dofs_per_cell_U = fe_U.dofs_per_cell;
2349 * std::vector<types::global_dof_index> local_dof_indices_U (dofs_per_cell_U);
2352 * cell_LS = dof_handler_LS.begin_active(),
2353 * endc_LS = dof_handler_LS.end();
2355 * cell_U = dof_handler_U.begin_active();
2357 *
for (; cell_LS!=endc_LS; ++cell_LS, ++cell_U)
2358 *
if (!cell_LS->is_artificial())
2360 * cell_LS->get_dof_indices(local_dof_indices_LS);
2361 * cell_U->get_dof_indices(local_dof_indices_U);
2362 *
for (
unsigned int i=0; i<dofs_per_cell_LS; ++i)
2363 * map_from_Q1_to_Q2[local_dof_indices_LS[i]] = local_dof_indices_U[i];
2367 *
template <
int dim>
2370 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
2376 * all vectors are NON-GHOSTED
2379 *
SolverControl solver_control (dof_handler_LS.n_dofs(), solver_tolerance);
2381 * constraints.distribute (completely_distributed_solution);
2382 * solver.solve (Matrix, completely_distributed_solution, rhs, *preconditioner);
2383 * constraints.distribute (completely_distributed_solution);
2384 *
if (verbose==
true) pcout <<
" Solved in " << solver_control.last_step() <<
" iterations." << std::endl;
2387 *
template <
int dim>
2388 *
void LevelSetSolver<dim>::save_old_solution()
2394 *
template <
int dim>
2395 *
void LevelSetSolver<dim>::save_old_vel_solution()
2397 * locally_relevant_solution_vx_old = locally_relevant_solution_vx;
2398 * locally_relevant_solution_vy_old = locally_relevant_solution_vy;
2400 * locally_relevant_solution_vz_old = locally_relevant_solution_vz;
2405 * -------------------------------------------------------------------------------
2406 * ------------------------------ MY PETSC WRAPPERS ------------------------------
2407 * -------------------------------------------------------------------------------
2412 *
const std::vector<types::global_dof_index> &indices,
2413 * std::vector<PetscScalar> &values)
2417 *
PETSc wrapper to get sets of
values from a petsc vector.
2418 * we assume the vector is ghosted
2419 * We need to figure out which elements we
2420 * own locally. Then get a pointer to the
2421 * elements that are stored here (both the
2422 * ones we own as well as the ghost elements).
2423 * In
this array, the locally owned elements
2424 * come
first followed by the ghost elements whose
2425 * position we can get from an
index set
2431 *
IndexSet ghost_indices = locally_relevant_dofs_LS;
2435 * n_idx = indices.size();
2437 * VecGetOwnershipRange (vector, &begin, &end);
2439 * Vec solution_in_local_form =
nullptr;
2440 * VecGhostGetLocalForm(vector, &solution_in_local_form);
2442 * PetscScalar *soln;
2443 * VecGetArray(solution_in_local_form, &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);
2456 * VecRestoreArray(solution_in_local_form, &soln);
2457 * VecGhostRestoreLocalForm(vector, &solution_in_local_form);
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)
2468 * THIS IS MEANT TO BE USED WITH VELOCITY VECTORS
2469 *
PETSc wrapper to get sets of
values from a petsc vector.
2470 * we assume the vector is ghosted
2471 * We need to figure out which elements we
2472 * own locally. Then get a pointer to the
2473 * elements that are stored here (both the
2474 * ones we own as well as the ghost elements).
2475 * In
this array, the locally owned elements
2476 * come
first followed by the ghost elements whose
2477 * position we can get from an
index set
2483 *
IndexSet ghost_indices = locally_relevant_dofs_U;
2487 * n_idx = indices.size();
2489 * VecGetOwnershipRange (vector, &begin, &end);
2491 * Vec solution_in_local_form =
nullptr;
2492 * VecGhostGetLocalForm(vector, &solution_in_local_form);
2494 * PetscScalar *soln;
2495 * VecGetArray(solution_in_local_form, &soln);
2497 *
for (i = 0; i < n_idx; ++i)
2499 *
int index = map_from_Q1_to_Q2[indices[i]];
2500 *
if (index >= begin && index < end)
2504 *
const unsigned int ghostidx = ghost_indices.index_within_set(index);
2508 * VecRestoreArray(solution_in_local_form, &soln);
2509 * VecGhostRestoreLocalForm(vector, &solution_in_local_form);
2515<a name=
"ann-MultiPhase.cc"></a>
2516<h1>Annotated version of MultiPhase.cc</h1>
2532 * #include <deal.II/base/quadrature_lib.h>
2533 * #include <deal.II/base/function.h>
2534 * #include <deal.II/lac/affine_constraints.h>
2535 * #include <deal.II/lac/vector.h>
2536 * #include <deal.II/lac/full_matrix.h>
2537 * #include <deal.II/lac/solver_cg.h>
2538 * #include <deal.II/lac/petsc_sparse_matrix.h>
2539 * #include <deal.II/lac/petsc_vector.h>
2540 * #include <deal.II/lac/petsc_solver.h>
2541 * #include <deal.II/lac/petsc_precondition.h>
2542 * #include <deal.II/grid/grid_generator.h>
2543 * #include <deal.II/grid/tria_accessor.h>
2544 * #include <deal.II/grid/tria_iterator.h>
2545 * #include <deal.II/dofs/dof_handler.h>
2546 * #include <deal.II/dofs/dof_accessor.h>
2547 * #include <deal.II/dofs/dof_tools.h>
2548 * #include <deal.II/fe/fe_values.h>
2549 * #include <deal.II/fe/fe_q.h>
2550 * #include <deal.II/numerics/vector_tools.h>
2551 * #include <deal.II/numerics/data_out.h>
2552 * #include <deal.II/numerics/error_estimator.h>
2553 * #include <deal.II/base/utilities.h>
2554 * #include <deal.II/base/conditional_ostream.h>
2555 * #include <deal.II/base/index_set.h>
2556 * #include <deal.II/lac/sparsity_tools.h>
2557 * #include <deal.II/distributed/tria.h>
2558 * #include <deal.II/distributed/grid_refinement.h>
2559 * #include <deal.II/base/convergence_table.h>
2560 * #include <deal.II/base/timer.h>
2561 * #include <deal.II/base/parameter_handler.h>
2562 * #include <fstream>
2563 * #include <iostream>
2564 * #include <deal.II/grid/grid_tools.h>
2565 * #include <deal.II/fe/mapping_q.h>
2567 *
using namespace dealii;
2572 * FOR TRANSPORT PROBLEM
2584 * #define FILLING_TANK 0
2585 * #define BREAKING_DAM 1
2586 * #define FALLING_DROP 2
2587 * #define SMALL_WAVE_PERTURBATION 3
2589 * #include
"NavierStokesSolver.cc"
2590 * #include
"LevelSetSolver.cc"
2591 * #include
"utilities.cc"
2600 *
template <
int dim>
2604 * MultiPhase (
const unsigned int degree_LS,
2605 *
const unsigned int degree_U);
2610 *
void set_boundary_inlet();
2611 *
void get_boundary_values_U();
2612 *
void get_boundary_values_phi(std::vector<types::global_dof_index> &boundary_values_id_phi,
2613 * std::vector<double> &boundary_values_phi);
2614 *
void output_results();
2615 *
void output_vectors();
2616 *
void output_rho();
2618 *
void initial_condition();
2619 *
void init_constraints();
2628 *
IndexSet locally_relevant_dofs_LS;
2635 *
IndexSet locally_relevant_dofs_U;
2640 *
IndexSet locally_relevant_dofs_P;
2662 * std::vector<types::global_dof_index> boundary_values_id_u;
2663 * std::vector<types::global_dof_index> boundary_values_id_v;
2664 * std::vector<types::global_dof_index> boundary_values_id_phi;
2665 * std::vector<double> boundary_values_u;
2666 * std::vector<double> boundary_values_v;
2667 * std::vector<double> boundary_values_phi;
2673 *
double final_time;
2674 *
unsigned int timestep_number;
2680 *
int sharpness_integer;
2682 *
unsigned int n_refinement;
2683 *
unsigned int output_number;
2684 *
double output_time;
2708 *
unsigned int TRANSPORT_TIME_INTEGRATION;
2709 * std::string ALGORITHM;
2710 *
unsigned int PROBLEM;
2713 *
template <
int dim>
2714 * MultiPhase<dim>::MultiPhase (
const unsigned int degree_LS,
2715 *
const unsigned int degree_U)
2717 * mpi_communicator (MPI_COMM_WORLD),
2722 * degree_LS(degree_LS),
2724 * fe_LS (degree_LS),
2725 * degree_U(degree_U),
2729 * fe_P (degree_U-1),
2733 *
template <
int dim>
2734 * MultiPhase<dim>::~MultiPhase ()
2736 * dof_handler_LS.clear ();
2737 * dof_handler_U.clear ();
2738 * dof_handler_P.clear ();
2748 *
template <
int dim>
2749 *
void MultiPhase<dim>::setup()
2756 * dof_handler_LS.distribute_dofs (fe_LS);
2757 * locally_owned_dofs_LS = dof_handler_LS.locally_owned_dofs ();
2764 * dof_handler_U.distribute_dofs (fe_U);
2765 * locally_owned_dofs_U = dof_handler_U.locally_owned_dofs ();
2772 * dof_handler_P.distribute_dofs (fe_P);
2773 * locally_owned_dofs_P = dof_handler_P.locally_owned_dofs ();
2777 * init vectors
for phi
2780 * locally_relevant_solution_phi.reinit(locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
2781 * locally_relevant_solution_phi = 0;
2782 * completely_distributed_solution_phi.reinit (locally_owned_dofs_P,mpi_communicator);
2785 * init vectors
for u
2788 * locally_relevant_solution_u.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
2789 * locally_relevant_solution_u = 0;
2790 * completely_distributed_solution_u.reinit (locally_owned_dofs_U,mpi_communicator);
2793 * init vectors
for v
2796 * locally_relevant_solution_v.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
2797 * locally_relevant_solution_v = 0;
2798 * completely_distributed_solution_v.reinit (locally_owned_dofs_U,mpi_communicator);
2801 * init vectors
for p
2804 * locally_relevant_solution_p.reinit (locally_owned_dofs_P,locally_relevant_dofs_P,mpi_communicator);
2805 * locally_relevant_solution_p = 0;
2806 * completely_distributed_solution_p.reinit (locally_owned_dofs_P,mpi_communicator);
2812 * init_constraints();
2815 *
template <
int dim>
2816 *
void MultiPhase<dim>::initial_condition()
2821 * Initial conditions
2822 * init condition
for phi
2825 * completely_distributed_solution_phi = 0;
2827 * InitialPhi<dim>(PROBLEM, sharpness),
2828 * completely_distributed_solution_phi);
2829 * constraints.distribute (completely_distributed_solution_phi);
2830 * locally_relevant_solution_phi = completely_distributed_solution_phi;
2833 * init condition
for u=0
2836 * completely_distributed_solution_u = 0;
2839 * completely_distributed_solution_u);
2840 * constraints.distribute (completely_distributed_solution_u);
2841 * locally_relevant_solution_u = completely_distributed_solution_u;
2844 * init condition
for v
2847 * completely_distributed_solution_v = 0;
2850 * completely_distributed_solution_v);
2851 * constraints.distribute (completely_distributed_solution_v);
2852 * locally_relevant_solution_v = completely_distributed_solution_v;
2855 * init condition
for p
2858 * completely_distributed_solution_p = 0;
2861 * completely_distributed_solution_p);
2862 * constraints.distribute (completely_distributed_solution_p);
2863 * locally_relevant_solution_p = completely_distributed_solution_p;
2866 *
template <
int dim>
2867 *
void MultiPhase<dim>::init_constraints()
2869 * constraints.clear ();
2870 * constraints.reinit (locally_relevant_dofs_LS);
2872 * constraints.close ();
2875 *
template <
int dim>
2876 *
void MultiPhase<dim>::get_boundary_values_U()
2878 * std::map<types::global_dof_index, double> map_boundary_values_u;
2879 * std::map<types::global_dof_index, double> map_boundary_values_v;
2880 * std::map<types::global_dof_index, double> map_boundary_values_w;
2887 *
if (PROBLEM==BREAKING_DAM || PROBLEM==FALLING_DROP)
2918 *
else if (PROBLEM==SMALL_WAVE_PERTURBATION)
2922 * no slip in bottom and top and slip in left and right
2948 *
else if (PROBLEM==FILLING_TANK)
2952 * LEFT: entry in x, zero in y
2959 * RIGHT: no-slip condition
2973 * TOP: exit in y, zero in x
2981 * pcout <<
"Error in type of PROBLEM at Boundary Conditions" << std::endl;
2984 * boundary_values_id_u.resize(map_boundary_values_u.size());
2985 * boundary_values_id_v.resize(map_boundary_values_v.size());
2986 * boundary_values_u.resize(map_boundary_values_u.size());
2987 * boundary_values_v.resize(map_boundary_values_v.size());
2988 * std::map<types::global_dof_index,double>::const_iterator boundary_value_u =map_boundary_values_u.begin();
2989 * std::map<types::global_dof_index,double>::const_iterator boundary_value_v =map_boundary_values_v.begin();
2991 *
for (
int i=0; boundary_value_u !=map_boundary_values_u.end(); ++boundary_value_u, ++i)
2993 * boundary_values_id_u[i]=boundary_value_u->first;
2994 * boundary_values_u[i]=boundary_value_u->second;
2996 *
for (
int i=0; boundary_value_v !=map_boundary_values_v.end(); ++boundary_value_v, ++i)
2998 * boundary_values_id_v[i]=boundary_value_v->first;
2999 * boundary_values_v[i]=boundary_value_v->second;
3003 *
template <
int dim>
3004 *
void MultiPhase<dim>::set_boundary_inlet()
3006 *
const QGauss<dim-1> face_quadrature_formula(1);
3010 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
3011 * std::vector<double> u_value (n_face_q_points);
3012 * std::vector<double> v_value (n_face_q_points);
3015 * cell_U = dof_handler_U.begin_active(),
3016 * endc_U = dof_handler_U.end();
3019 *
for (; cell_U!=endc_U; ++cell_U)
3020 *
if (cell_U->is_locally_owned())
3021 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
3022 *
if (cell_U->face(face)->at_boundary())
3024 * fe_face_values.reinit(cell_U,face);
3025 * fe_face_values.get_function_values(locally_relevant_solution_u,u_value);
3026 * fe_face_values.get_function_values(locally_relevant_solution_v,v_value);
3029 *
if (fe_face_values.normal_vector(0)*u < -1e-14)
3030 * cell_U->face(face)->set_boundary_id(10);
3034 *
template <
int dim>
3035 *
void MultiPhase<dim>::get_boundary_values_phi(std::vector<types::global_dof_index> &boundary_values_id_phi,
3036 * std::vector<double> &boundary_values_phi)
3038 * std::map<types::global_dof_index, double> map_boundary_values_phi;
3041 * set_boundary_inlet();
3044 * boundary_values_id_phi.resize(map_boundary_values_phi.size());
3045 * boundary_values_phi.resize(map_boundary_values_phi.size());
3046 * std::map<types::global_dof_index,double>::const_iterator boundary_value_phi = map_boundary_values_phi.begin();
3047 *
for (
int i=0; boundary_value_phi !=map_boundary_values_phi.end(); ++boundary_value_phi, ++i)
3049 * boundary_values_id_phi[i]=boundary_value_phi->first;
3050 * boundary_values_phi[i]=boundary_value_phi->second;
3055 *
void MultiPhase<dim>::output_results()
3066 *
template <
int dim>
3067 *
void MultiPhase<dim>::output_vectors()
3071 * data_out.add_data_vector (locally_relevant_solution_phi,
"phi");
3072 * data_out.build_patches ();
3074 *
const std::string filename = (
"sol_vectors-" +
3079 * std::ofstream output ((filename +
".vtu").c_str());
3080 * data_out.write_vtu (output);
3084 * std::vector<std::string> filenames;
3085 *
for (
unsigned int i=0;
3088 * filenames.push_back (
"sol_vectors-" +
3094 * std::ofstream master_output ((filename +
".pvtu").c_str());
3095 * data_out.write_pvtu_record (master_output, filenames);
3099 *
template <
int dim>
3100 *
void MultiPhase<dim>::output_rho()
3102 * Postprocessor<dim> postprocessor(eps,rho_air,rho_fluid);
3105 * data_out.add_data_vector (locally_relevant_solution_phi, postprocessor);
3107 * data_out.build_patches ();
3109 *
const std::string filename = (
"sol_rho-" +
3114 * std::ofstream output ((filename +
".vtu").c_str());
3115 * data_out.write_vtu (output);
3119 * std::vector<std::string> filenames;
3120 *
for (
unsigned int i=0;
3123 * filenames.push_back (
"sol_rho-" +
3129 * std::ofstream master_output ((filename +
".pvtu").c_str());
3130 * data_out.write_pvtu_record (master_output, filenames);
3134 *
template <
int dim>
3135 *
void MultiPhase<dim>::run()
3140 * GENERAL PARAMETERS
3147 * get_output =
true;
3148 * output_number = 0;
3150 * output_time = 0.1;
3151 * final_time = 10.0;
3155 * PARAMETERS FOR THE NAVIER STOKES PROBLEM
3159 * rho_fluid = 1000.;
3163 * PROBLEM=BREAKING_DAM;
3166 * PROBLEM=FILLING_TANK;
3167 * PROBLEM=SMALL_WAVE_PERTURBATION;
3168 * PROBLEM=FALLING_DROP;
3174 * ForceTerms<dim> force_function(std::vector<double> {0.0,-1.0});
3178 * PARAMETERS FOR TRANSPORT PROBLEM
3184 * sharpness_integer=10;
3190 * TRANSPORT_TIME_INTEGRATION=SSP33;
3193 * ALGORITHM =
"MPP_u1";
3194 * ALGORITHM =
"NMPP_uH";
3197 * ALGORITHM =
"MPP_uH";
3201 * ADJUST PARAMETERS ACCORDING TO PROBLEM
3204 *
if (PROBLEM==FALLING_DROP)
3214 *
if (PROBLEM==FILLING_TANK)
3217 *
else if (PROBLEM==BREAKING_DAM || PROBLEM==SMALL_WAVE_PERTURBATION)
3219 * std::vector< unsigned int > repetitions;
3220 * repetitions.push_back(2);
3221 * repetitions.push_back(1);
3225 *
else if (PROBLEM==FALLING_DROP)
3227 * std::vector< unsigned int > repetitions;
3228 * repetitions.push_back(1);
3229 * repetitions.push_back(4);
3243 * PARAMETERS FOR TIME STEPPING
3247 * time_step = cfl*min_h/umax;
3249 * sharpness=sharpness_integer*min_h;
3253 * INITIAL CONDITIONS
3256 * initial_condition();
3261 * NAVIER STOKES SOLVER
3264 * NavierStokesSolver<dim> navier_stokes (degree_LS,degree_U,
3267 * rho_fluid,nu_fluid,
3273 * BOUNDARY CONDITIONS FOR NAVIER STOKES
3276 * get_boundary_values_U();
3277 * navier_stokes.set_boundary_conditions(boundary_values_id_u, boundary_values_id_v,
3278 * boundary_values_u, boundary_values_v);
3282 * set INITIAL CONDITION within NAVIER STOKES
3285 * navier_stokes.initial_condition(locally_relevant_solution_phi,
3286 * locally_relevant_solution_u,
3287 * locally_relevant_solution_v,
3288 * locally_relevant_solution_p);
3294 * LevelSetSolver<dim> transport_solver (degree_LS,degree_U,
3298 * TRANSPORT_TIME_INTEGRATION,
3300 * mpi_communicator);
3303 * BOUNDARY CONDITIONS FOR PHI
3306 * get_boundary_values_phi(boundary_values_id_phi,boundary_values_phi);
3307 * transport_solver.set_boundary_conditions(boundary_values_id_phi,boundary_values_phi);
3311 * set INITIAL CONDITION within TRANSPORT PROBLEM
3314 * transport_solver.initial_condition(locally_relevant_solution_phi,
3315 * locally_relevant_solution_u,
3316 * locally_relevant_solution_v);
3317 *
int dofs_U = 2*dof_handler_U.n_dofs();
3318 *
int dofs_P = 2*dof_handler_P.n_dofs();
3319 *
int dofs_LS = dof_handler_LS.n_dofs();
3320 *
int dofs_TOTAL = dofs_U+dofs_P+dofs_LS;
3324 * NO BOUNDARY CONDITIONS
for LEVEL SET
3327 * pcout <<
"Cfl: " << cfl <<
"; umax: " << umax <<
"; min h: " << min_h
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
3335 * <<
" TOTAL: " << dofs_TOTAL
3343 *
for (timestep_number=1, time=time_step; time<=final_time;
3344 * time+=time_step,++timestep_number)
3346 * pcout <<
"Time step " << timestep_number
3347 * <<
" at t=" << time
3351 * GET NAVIER STOKES VELOCITY
3354 * navier_stokes.set_phi(locally_relevant_solution_phi);
3355 * navier_stokes.nth_time_step();
3356 * navier_stokes.get_velocity(locally_relevant_solution_u,locally_relevant_solution_v);
3357 * transport_solver.set_velocity(locally_relevant_solution_u,locally_relevant_solution_v);
3360 * GET LEVEL SET SOLUTION
3363 * transport_solver.nth_time_step();
3364 * transport_solver.get_unp1(locally_relevant_solution_phi);
3365 *
if (get_output && time-(output_number)*output_time>0)
3368 * navier_stokes.get_velocity(locally_relevant_solution_u, locally_relevant_solution_v);
3369 * transport_solver.get_unp1(locally_relevant_solution_phi);
3374 *
int main(
int argc,
char *argv[])
3378 *
using namespace dealii;
3380 * PetscInitialize(&argc, &argv,
nullptr,
nullptr);
3383 *
unsigned int degree_LS = 1;
3384 *
unsigned int degree_U = 2;
3385 * MultiPhase<2> multi_phase(degree_LS, degree_U);
3386 * multi_phase.run();
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>
3420<h1>Annotated version of NavierStokesSolver.cc</h1>
3436 * #include <deal.II/base/quadrature_lib.h>
3437 * #include <deal.II/base/function.h>
3438 * #include <deal.II/lac/affine_constraints.h>
3439 * #include <deal.II/lac/vector.h>
3440 * #include <deal.II/lac/full_matrix.h>
3441 * #include <deal.II/lac/solver_cg.h>
3442 * #include <deal.II/lac/petsc_sparse_matrix.h>
3443 * #include <deal.II/lac/petsc_vector.h>
3444 * #include <deal.II/lac/petsc_solver.h>
3445 * #include <deal.II/lac/petsc_precondition.h>
3446 * #include <deal.II/grid/grid_generator.h>
3447 * #include <deal.II/grid/tria_accessor.h>
3448 * #include <deal.II/grid/tria_iterator.h>
3449 * #include <deal.II/dofs/dof_handler.h>
3450 * #include <deal.II/dofs/dof_accessor.h>
3451 * #include <deal.II/dofs/dof_tools.h>
3452 * #include <deal.II/fe/fe_values.h>
3453 * #include <deal.II/fe/fe_q.h>
3454 * #include <deal.II/numerics/vector_tools.h>
3455 * #include <deal.II/numerics/data_out.h>
3456 * #include <deal.II/numerics/error_estimator.h>
3457 * #include <deal.II/base/utilities.h>
3458 * #include <deal.II/base/conditional_ostream.h>
3459 * #include <deal.II/base/index_set.h>
3460 * #include <deal.II/lac/sparsity_tools.h>
3461 * #include <deal.II/distributed/tria.h>
3462 * #include <deal.II/distributed/grid_refinement.h>
3463 * #include <deal.II/lac/petsc_vector.h>
3464 * #include <deal.II/base/convergence_table.h>
3465 * #include <deal.II/base/timer.h>
3466 * #include <deal.II/base/parameter_handler.h>
3467 * #include <deal.II/grid/grid_tools.h>
3468 * #include <deal.II/fe/mapping_q.h>
3470 * #include <fstream>
3471 * #include <iostream>
3474 *
using namespace dealii;
3476 * #define MAX_NUM_ITER_TO_RECOMPUTE_PRECONDITIONER 10
3486 *
class NavierStokesSolver
3491 * constructor
for using LEVEL SET
3494 * NavierStokesSolver(
const unsigned int degree_LS,
3495 *
const unsigned int degree_U,
3496 *
const double time_step,
3498 *
const double rho_air,
3499 *
const double nu_air,
3500 *
const double rho_fluid,
3501 *
const double nu_fluid,
3503 *
const bool verbose,
3508 * constructor
for NOT LEVEL SET
3511 * NavierStokesSolver(
const unsigned int degree_LS,
3512 *
const unsigned int degree_U,
3513 *
const double time_step,
3517 *
const bool verbose,
3526 *
void set_rho_and_nu_functions(
const Function<dim> &rho_function,
3544 * boundary conditions
3547 *
void set_boundary_conditions(std::vector<types::global_dof_index> boundary_values_id_u,
3548 * std::vector<types::global_dof_index> boundary_values_id_v, std::vector<double> boundary_values_u,
3549 * std::vector<double> boundary_values_v);
3550 *
void set_boundary_conditions(std::vector<types::global_dof_index> boundary_values_id_u,
3551 * std::vector<types::global_dof_index> boundary_values_id_v,
3552 * std::vector<types::global_dof_index> boundary_values_id_w, std::vector<double> boundary_values_u,
3553 * std::vector<double> boundary_values_v, std::vector<double> boundary_values_w);
3571 *
void nth_time_step();
3579 * ~NavierStokesSolver();
3584 * SETUP AND INITIAL CONDITION
3588 *
void setup_VECTORS();
3589 *
void init_constraints();
3595 *
void assemble_system_U();
3596 *
void assemble_system_dpsi_q();
3603 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
3607 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
3612 * GET DIFFERENT FIELDS
3615 *
void get_rho_and_nu(
double phi);
3616 *
void get_velocity();
3617 *
void get_pressure();
3623 *
void save_old_solution();
3632 *
IndexSet locally_relevant_dofs_LS;
3638 *
IndexSet locally_relevant_dofs_U;
3643 *
IndexSet locally_relevant_dofs_P;
3658 *
unsigned int LEVEL_SET;
3659 *
unsigned int RHO_TIMES_RHS;
3675 * std::vector<types::global_dof_index> boundary_values_id_u;
3676 * std::vector<types::global_dof_index> boundary_values_id_v;
3677 * std::vector<types::global_dof_index> boundary_values_id_w;
3678 * std::vector<double> boundary_values_u;
3679 * std::vector<double> boundary_values_v;
3680 * std::vector<double> boundary_values_w;
3685 *
bool rebuild_Matrix_U;
3686 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_Matrix_u;
3687 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_Matrix_v;
3688 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_Matrix_w;
3690 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_S;
3692 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_M;
3694 *
bool rebuild_Matrix_U_preconditioners;
3695 *
bool rebuild_S_M_preconditioners;
3723 * CONSTRUCTOR FOR LEVEL SET
3727 * NavierStokesSolver<dim>::NavierStokesSolver(
const unsigned int degree_LS,
3728 *
const unsigned int degree_U,
3729 *
const double time_step,
3731 *
const double rho_air,
3732 *
const double nu_air,
3733 *
const double rho_fluid,
3734 *
const double nu_fluid,
3736 *
const bool verbose,
3740 * mpi_communicator(mpi_communicator),
3742 * degree_LS(degree_LS),
3745 * degree_U(degree_U),
3750 * force_function(force_function),
3753 * This is dummy since rho and nu
functions won
't be used
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
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
void subtract_set(const IndexSet &other)
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.
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > E(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