247 * #include <deal.II/base/quadrature_lib.h>
248 * #include <deal.II/base/function.h>
249 * #include <deal.II/lac/affine_constraints.h>
250 * #include <deal.II/lac/vector.h>
251 * #include <deal.II/lac/full_matrix.h>
252 * #include <deal.II/lac/solver_cg.h>
253 * #include <deal.II/lac/petsc_sparse_matrix.h>
254 * #include <deal.II/lac/petsc_sparse_matrix.h>
255 * #include <deal.II/lac/petsc_vector.h>
256 * #include <deal.II/lac/petsc_solver.h>
257 * #include <deal.II/lac/petsc_precondition.h>
258 * #include <deal.II/grid/grid_generator.h>
259 * #include <deal.II/grid/tria_accessor.h>
260 * #include <deal.II/grid/tria_iterator.h>
261 * #include <deal.II/dofs/dof_handler.h>
262 * #include <deal.II/dofs/dof_accessor.h>
263 * #include <deal.II/dofs/dof_tools.h>
264 * #include <deal.II/fe/fe_values.h>
265 * #include <deal.II/fe/fe_q.h>
266 * #include <deal.II/numerics/vector_tools.h>
267 * #include <deal.II/numerics/data_out.h>
268 * #include <deal.II/numerics/error_estimator.h>
269 * #include <deal.II/base/utilities.h>
270 * #include <deal.II/base/conditional_ostream.h>
271 * #include <deal.II/base/index_set.h>
272 * #include <deal.II/lac/sparsity_tools.h>
273 * #include <deal.II/distributed/tria.h>
274 * #include <deal.II/distributed/grid_refinement.h>
275 * #include <deal.II/lac/vector.h>
276 * #include <deal.II/base/convergence_table.h>
277 * #include <deal.II/base/timer.h>
278 * #include <deal.II/base/parameter_handler.h>
279 * #include <deal.II/grid/grid_tools.h>
280 * #include <deal.II/fe/mapping_q.h>
285 * #include <iostream>
296 * #define CHECK_MAX_PRINCIPLE 0
300 * LOG FOR LEVEL SET FROM -1 to 1
304 * #define ENTROPY_GRAD(phi,phix) 2*phi*phix*((1-phi*phi>=0) ? -1 : 1)/(
std::
abs(1-phi*phi)+1
E-14)
311 * This is a solver for the transpor solver.
312 * We assume the velocity is divergence
free
313 * and solve the equation in conservation form.
315 * ---------- NOTATION ----------
317 * We use notation popular in the literature of conservation laws.
318 * For this reason the solution is denoted as u, unm1, unp1, etc.
319 * and the velocity is treated as vx, vy and vz.
323 * class LevelSetSolver
343 * BOUNDARY CONDITIONS
347 *
void set_boundary_conditions(std::vector<types::global_dof_index> &boundary_values_id_u,
348 * std::vector<double> boundary_values_u);
376 *
void nth_time_step();
386 * LevelSetSolver (
const unsigned int degree_LS,
387 *
const unsigned int degree_U,
388 *
const double time_step,
391 *
const bool verbose,
392 * std::string ALGORITHM,
393 *
const unsigned int TIME_INTEGRATION,
402 * ASSEMBLE MASS (and other) MATRICES
406 *
void assemble_ML();
408 *
void assemble_MC();
412 * LOW ORDER METHOD (DiJ Viscosity)
416 *
void assemble_C_Matrix();
426 *
void assemble_EntRes_Matrix();
430 * FOR MAXIMUM PRINCIPLE
456 * std::string algorithm);
459 * std::string algorithm);
467 *
void get_sparsity_pattern();
468 *
void get_map_from_Q1_to_Q2();
471 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
474 *
void save_old_solution();
475 *
void save_old_vel_solution();
484 *
const std::vector<types::global_dof_index> &indices,
485 * std::vector<PetscScalar> &values);
487 *
const std::vector<types::global_dof_index> &indices,
488 * std::map<types::global_dof_index, types::global_dof_index> &map_from_Q1_to_Q2,
489 * std::vector<PetscScalar> &values);
495 * FINITE ELEMENT SPACE
503 *
IndexSet locally_relevant_dofs_LS;
513 * OPERATORS times SOLUTION VECTOR
526 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> MC_preconditioner;
533 * std::vector<types::global_dof_index> boundary_values_id_u;
534 * std::vector<double> boundary_values_u;
541 * FOR FIRST ORDER VISCOSITY
548 * FOR ENTROPY VISCOSITY
554 * FOR FCT (flux and limited flux)
583 * NON-GHOSTED VECTORS
621 *
double solver_tolerance;
622 *
double entropy_normalization_factor;
630 * std::string ALGORITHM;
631 *
unsigned int TIME_INTEGRATION;
635 * std::map<types::global_dof_index, types::global_dof_index> map_from_Q1_to_Q2;
636 * std::map<types::global_dof_index, std::vector<types::global_dof_index> > sparsity_pattern;
640 * LevelSetSolver<dim>::LevelSetSolver (
const unsigned int degree_LS,
641 *
const unsigned int degree_U,
642 *
const double time_step,
645 *
const bool verbose,
646 * std::string ALGORITHM,
647 *
const unsigned int TIME_INTEGRATION,
651 * mpi_communicator (mpi_communicator),
652 * degree_LS(degree_LS),
655 * degree_U(degree_U),
658 * time_step(time_step),
662 * ALGORITHM(ALGORITHM),
663 * TIME_INTEGRATION(TIME_INTEGRATION),
666 * pcout <<
"********** LEVEL SET SETUP **********" << std::endl;
671 * LevelSetSolver<dim>::~LevelSetSolver ()
673 * dof_handler_LS.clear ();
674 * dof_handler_U.clear ();
693 * this->locally_relevant_solution_vx = locally_relevant_solution_vx;
694 * this->locally_relevant_solution_vy = locally_relevant_solution_vy;
697 * initialize old vectors with current solution,
this just happens the
first time
701 * locally_relevant_solution_vx_old = locally_relevant_solution_vx;
702 * locally_relevant_solution_vy_old = locally_relevant_solution_vy;
712 * this->locally_relevant_solution_vx = locally_relevant_solution_vx;
713 * this->locally_relevant_solution_vy = locally_relevant_solution_vy;
714 * this->locally_relevant_solution_vz = locally_relevant_solution_vz;
717 * initialize old vectors with current solution,
this just happens the
first time
721 * locally_relevant_solution_vx_old = locally_relevant_solution_vx;
722 * locally_relevant_solution_vy_old = locally_relevant_solution_vy;
723 * locally_relevant_solution_vz_old = locally_relevant_solution_vz;
734 *
void LevelSetSolver<dim>::set_boundary_conditions(std::vector<types::global_dof_index> &boundary_values_id_u,
735 * std::vector<double> boundary_values_u)
737 * this->boundary_values_id_u = boundary_values_id_u;
738 * this->boundary_values_u = boundary_values_u;
757 * save_old_vel_solution();
763 * this->locally_relevant_solution_vx=locally_relevant_solution_vx;
764 * this->locally_relevant_solution_vy=locally_relevant_solution_vy;
777 * save_old_vel_solution();
783 * this->locally_relevant_solution_vx=locally_relevant_solution_vx;
784 * this->locally_relevant_solution_vy=locally_relevant_solution_vy;
785 * this->locally_relevant_solution_vz=locally_relevant_solution_vz;
800 * -------------------------------------------------------------------------------
801 * ------------------------------ COMPUTE SOLUTIONS ------------------------------
802 * -------------------------------------------------------------------------------
806 *
void LevelSetSolver<dim>::nth_time_step()
808 * assemble_EntRes_Matrix();
814 *
if (TIME_INTEGRATION==FORWARD_EULER)
815 * compute_solution(unp1,un,ALGORITHM);
817 * compute_solution_SSP33(unp1,un,ALGORITHM);
820 * BOUNDARY CONDITIONS
823 * unp1.set(boundary_values_id_u,boundary_values_u);
827 * CHECK MAXIMUM PRINCIPLE
830 *
if (CHECK_MAX_PRINCIPLE)
832 * compute_bounds(un);
833 * check_max_principle(unp1);
837 * pcout <<
"*********************************************************************... "
838 * << unp1.min() <<
", " << unp1.max() << std::endl;
841 * save_old_solution();
846 * --------------------------------------------------------------------
847 * ------------------------------ SETUP ------------------------------
848 * --------------------------------------------------------------------
852 *
void LevelSetSolver<dim>::setup()
854 * solver_tolerance=1
E-6;
855 * degree_MAX =
std::max(degree_LS,degree_U);
859 * SETUP FOR DOF HANDLERS
864 * dof_handler_LS.distribute_dofs (fe_LS);
865 * locally_owned_dofs_LS = dof_handler_LS.locally_owned_dofs ();
872 * dof_handler_U.distribute_dofs (fe_U);
873 * locally_owned_dofs_U = dof_handler_U.locally_owned_dofs ();
882 * constraints.clear ();
883 * constraints.reinit (locally_relevant_dofs_LS);
885 * constraints.close ();
889 * NON-GHOSTED VECTORS
893 * MPP_uL_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
894 * NMPP_uH_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
895 * RHS.reinit(locally_owned_dofs_LS,mpi_communicator);
896 * uStage1_nonGhosted.reinit (locally_owned_dofs_LS,mpi_communicator);
897 * uStage2_nonGhosted.reinit (locally_owned_dofs_LS,mpi_communicator);
898 * unp1.reinit (locally_owned_dofs_LS,mpi_communicator);
899 * MPP_uH_solution.reinit (locally_owned_dofs_LS,mpi_communicator);
902 * vectors
for lumped mass
matrix
905 * ML_vector.reinit(locally_owned_dofs_LS,mpi_communicator);
906 * inverse_ML_vector.reinit(locally_owned_dofs_LS,mpi_communicator);
907 * ones_vector.reinit(locally_owned_dofs_LS,mpi_communicator);
911 * operators times solution
914 * K_times_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
915 * DL_times_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
916 * DH_times_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
922 * R_pos_vector_nonGhosted.reinit (locally_owned_dofs_LS,mpi_communicator);
923 * R_neg_vector_nonGhosted.reinit (locally_owned_dofs_LS,mpi_communicator);
924 * umin_vector.reinit (locally_owned_dofs_LS,mpi_communicator);
925 * umax_vector.reinit (locally_owned_dofs_LS,mpi_communicator);
929 * GHOSTED VECTORS (used within some assemble process)
933 * uStage1.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
934 * uStage2.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
935 * unm1.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
936 * un.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
937 * MPP_uL_solution_ghosted.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
938 * MPP_uLkp1_solution_ghosted.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
939 * NMPP_uH_solution_ghosted.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
942 * init vectors
for vx
945 * locally_relevant_solution_vx.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
946 * locally_relevant_solution_vx_old.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
949 * init vectors
for vy
952 * locally_relevant_solution_vy.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
953 * locally_relevant_solution_vy_old.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
956 * init vectors
for vz
959 * locally_relevant_solution_vz.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
960 * locally_relevant_solution_vz_old.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
966 * R_pos_vector.reinit(locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
967 * R_neg_vector.reinit(locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
979 * dof_handler_LS.locally_owned_dofs(),
981 * locally_relevant_dofs_LS);
982 * MC_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
983 * dof_handler_LS.locally_owned_dofs(),
984 * dsp, mpi_communicator);
985 * Cx_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
986 * dof_handler_LS.locally_owned_dofs(),
987 * dsp, mpi_communicator);
988 * CTx_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
989 * dof_handler_LS.locally_owned_dofs(),
990 * dsp, mpi_communicator);
991 * Cy_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
992 * dof_handler_LS.locally_owned_dofs(),
993 * dsp, mpi_communicator);
994 * CTy_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
995 * dof_handler_LS.locally_owned_dofs(),
996 * dsp, mpi_communicator);
999 * Cz_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1000 * dof_handler_LS.locally_owned_dofs(),
1001 * dsp, mpi_communicator);
1002 * CTz_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1003 * dof_handler_LS.locally_owned_dofs(),
1004 * dsp, mpi_communicator);
1006 * dLij_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1007 * dof_handler_LS.locally_owned_dofs(),
1008 * dsp, mpi_communicator);
1009 * EntRes_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1010 * dof_handler_LS.locally_owned_dofs(),
1011 * dsp, mpi_communicator);
1012 * SuppSize_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1013 * dof_handler_LS.locally_owned_dofs(),
1014 * dsp, mpi_communicator);
1015 * dCij_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1016 * dof_handler_LS.locally_owned_dofs(),
1017 * dsp, mpi_communicator);
1018 * A_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1019 * dof_handler_LS.locally_owned_dofs(),
1020 * dsp, mpi_communicator);
1021 * LxA_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1022 * dof_handler_LS.locally_owned_dofs(),
1023 * dsp, mpi_communicator);
1024 * Akp1_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1025 * dof_handler_LS.locally_owned_dofs(),
1026 * dsp, mpi_communicator);
1027 * LxAkp1_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1028 * dof_handler_LS.locally_owned_dofs(),
1029 * dsp, mpi_communicator);
1032 * COMPUTE MASS MATRICES (AND OTHERS) FOR FIRST TIME STEP
1038 * assemble_C_Matrix();
1041 * get mat
for DOFs between Q1 and Q2
1044 * get_map_from_Q1_to_Q2();
1045 * get_sparsity_pattern();
1050 * ----------------------------------------------------------------------------
1051 * ------------------------------ MASS MATRICES ------------------------------
1052 * ----------------------------------------------------------------------------
1056 *
void LevelSetSolver<dim>::assemble_ML()
1060 *
const QGauss<dim> quadrature_formula(degree_MAX+1);
1066 *
const unsigned int dofs_per_cell = fe_LS.dofs_per_cell;
1067 *
const unsigned int n_q_points = quadrature_formula.size();
1070 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1073 * cell_LS = dof_handler_LS.begin_active(),
1074 * endc_LS = dof_handler_LS.end();
1076 *
for (; cell_LS!=endc_LS; ++cell_LS)
1077 *
if (cell_LS->is_locally_owned())
1080 * fe_values_LS.reinit (cell_LS);
1081 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
1083 *
const double JxW = fe_values_LS.JxW(q_point);
1084 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1085 * cell_ML (i) += fe_values_LS.shape_value(i,q_point)*JxW;
1092 * cell_LS->get_dof_indices (local_dof_indices);
1093 * constraints.distribute_local_to_global (cell_ML,local_dof_indices,ML_vector);
1104 *
void LevelSetSolver<dim>::invert_ML()
1108 *
loop on locally owned i-DOFs (rows)
1112 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
1114 *
int gi = *idofs_iter;
1115 * inverse_ML_vector(gi) = 1./ML_vector(gi);
1121 *
void LevelSetSolver<dim>::assemble_MC()
1125 *
const QGauss<dim> quadrature_formula(degree_MAX+1);
1131 *
const unsigned int dofs_per_cell = fe_LS.dofs_per_cell;
1132 *
const unsigned int n_q_points = quadrature_formula.size();
1135 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1136 * std::vector<double> shape_values(dofs_per_cell);
1139 * cell_LS = dof_handler_LS.begin_active(),
1140 * endc_LS = dof_handler_LS.end();
1142 *
for (; cell_LS!=endc_LS; ++cell_LS)
1143 *
if (cell_LS->is_locally_owned())
1146 * fe_values_LS.reinit (cell_LS);
1147 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
1149 *
const double JxW = fe_values_LS.JxW(q_point);
1150 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1151 * shape_values[i] = fe_values_LS.shape_value(i,q_point);
1153 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1154 *
for (
unsigned int j=0; j<dofs_per_cell; ++j)
1155 * cell_MC(i,j) += shape_values[i]*shape_values[j]*JxW;
1162 * cell_LS->get_dof_indices (local_dof_indices);
1163 * constraints.distribute_local_to_global (cell_MC,local_dof_indices,MC_matrix);
1176 * ---------------------------------------------------------------------------------------
1177 * ------------------------------
LO METHOD (Dij Viscosity) ------------------------------
1178 * ---------------------------------------------------------------------------------------
1181 *
template <
int dim>
1182 *
void LevelSetSolver<dim>::assemble_C_Matrix ()
1191 *
const QGauss<dim> quadrature_formula(degree_MAX+1);
1197 *
const unsigned int dofs_per_cell_LS = fe_LS.dofs_per_cell;
1198 *
const unsigned int n_q_points = quadrature_formula.size();
1207 * std::vector<Tensor<1, dim> > shape_grads_LS(dofs_per_cell_LS);
1208 * std::vector<double> shape_values_LS(dofs_per_cell_LS);
1210 * std::vector<types::global_dof_index> local_dof_indices_LS (dofs_per_cell_LS);
1213 * cell_LS = dof_handler_LS.begin_active();
1214 * endc_LS = dof_handler_LS.end();
1216 *
for (; cell_LS!=endc_LS; ++cell_LS)
1217 *
if (cell_LS->is_locally_owned())
1229 * fe_values_LS.reinit (cell_LS);
1230 * cell_LS->get_dof_indices (local_dof_indices_LS);
1232 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
1234 *
const double JxW = fe_values_LS.JxW(q_point);
1235 *
for (
unsigned int i=0; i<dofs_per_cell_LS; ++i)
1237 * shape_values_LS[i] = fe_values_LS.shape_value(i,q_point);
1238 * shape_grads_LS [i] = fe_values_LS.shape_grad (i,q_point);
1241 *
for (
unsigned int i=0; i<dofs_per_cell_LS; ++i)
1242 *
for (
unsigned int j=0; j < dofs_per_cell_LS; ++j)
1244 * cell_Cij_x(i,j) += (shape_grads_LS[j][0])*shape_values_LS[i]*JxW;
1245 * cell_Cij_y(i,j) += (shape_grads_LS[j][1])*shape_values_LS[i]*JxW;
1246 * cell_Cji_x(i,j) += (shape_grads_LS[i][0])*shape_values_LS[j]*JxW;
1247 * cell_Cji_y(i,j) += (shape_grads_LS[i][1])*shape_values_LS[j]*JxW;
1250 * cell_Cij_z(i,j) += (shape_grads_LS[j][2])*shape_values_LS[i]*JxW;
1251 * cell_Cji_z(i,j) += (shape_grads_LS[i][2])*shape_values_LS[j]*JxW;
1260 * constraints.distribute_local_to_global(cell_Cij_x,local_dof_indices_LS,Cx_matrix);
1261 * constraints.distribute_local_to_global(cell_Cji_x,local_dof_indices_LS,CTx_matrix);
1262 * constraints.distribute_local_to_global(cell_Cij_y,local_dof_indices_LS,Cy_matrix);
1263 * constraints.distribute_local_to_global(cell_Cji_y,local_dof_indices_LS,CTy_matrix);
1266 * constraints.distribute_local_to_global(cell_Cij_z,local_dof_indices_LS,Cz_matrix);
1267 * constraints.distribute_local_to_global(cell_Cji_z,local_dof_indices_LS,CTz_matrix);
1289 * K_times_solution = 0;
1291 *
const QGauss<dim> quadrature_formula(degree_MAX+1);
1300 *
const unsigned int dofs_per_cell = fe_LS.dofs_per_cell;
1301 *
const unsigned int n_q_points = quadrature_formula.size();
1305 * std::vector<Tensor<1,dim> > un_grads (n_q_points);
1306 * std::vector<double> old_vx_values (n_q_points);
1307 * std::vector<double> old_vy_values (n_q_points);
1308 * std::vector<double> old_vz_values (n_q_points);
1310 * std::vector<double> shape_values(dofs_per_cell);
1311 * std::vector<Tensor<1,dim> > shape_grads(dofs_per_cell);
1315 * std::vector<types::global_dof_index> indices_LS (dofs_per_cell);
1323 * cell_LS = dof_handler_LS.begin_active(),
1324 * endc_LS = dof_handler_LS.end();
1326 * cell_U = dof_handler_U.begin_active();
1329 *
for (; cell_LS!=endc_LS; ++cell_U, ++cell_LS)
1330 *
if (cell_LS->is_locally_owned())
1332 * cell_K_times_solution=0;
1334 * fe_values_LS.reinit (cell_LS);
1335 * cell_LS->get_dof_indices (indices_LS);
1336 * fe_values_LS.get_function_gradients(solution,un_grads);
1338 * fe_values_U.reinit (cell_U);
1339 * fe_values_U.get_function_values(locally_relevant_solution_vx,old_vx_values);
1340 * fe_values_U.get_function_values(locally_relevant_solution_vy,old_vy_values);
1341 *
if (dim==3) fe_values_U.get_function_values(locally_relevant_solution_vz,old_vz_values);
1345 * compute cell_K_times_solution
1348 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
1350 * v[0] = old_vx_values[q_point];
1351 * v[1] = old_vy_values[q_point];
1352 *
if (dim==3) v[2] = old_vz_values[q_point];
1354 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1355 * cell_K_times_solution(i) += (v*un_grads[q_point])
1356 * *fe_values_LS.shape_value(i,q_point)*fe_values_LS.JxW(q_point);
1363 * constraints.distribute_local_to_global (cell_K_times_solution, indices_LS, K_times_solution);
1368 *
template <
int dim>
1369 *
void LevelSetSolver<dim>::assemble_K_DL_DH_times_vector
1374 * K_times_solution=0;
1377 * DL_times_solution=0;
1378 * DH_times_solution=0;
1382 * PetscInt ncolumns;
1383 *
const PetscInt *gj;
1384 *
const PetscScalar *Cxi, *Cyi, *Czi, *CTxi, *CTyi, *CTzi;
1385 *
const PetscScalar *EntResi, *SuppSizei, *MCi;
1392 *
loop on locally owned i-DOFs (rows)
1397 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
1399 * PetscInt gi = *idofs_iter;
1402 *
double ith_K_times_solution = 0;
1406 * read velocity of i-th DOF
1409 * vi[0] = locally_relevant_solution_vx(map_from_Q1_to_Q2[gi]);
1410 * vi[1] = locally_relevant_solution_vy(map_from_Q1_to_Q2[gi]);
1411 *
if (dim==3) vi[2] = locally_relevant_solution_vz(map_from_Q1_to_Q2[gi]);
1412 * solni = solution(gi);
1416 * get i-th row of
C matrices
1419 * MatGetRow(Cx_matrix,gi,&ncolumns,&gj,&Cxi);
1420 * MatGetRow(Cy_matrix,gi,&ncolumns,&gj,&Cyi);
1421 * MatGetRow(CTx_matrix,gi,&ncolumns,&gj,&CTxi);
1422 * MatGetRow(CTy_matrix,gi,&ncolumns,&gj,&CTyi);
1425 * MatGetRow(Cz_matrix,gi,&ncolumns,&gj,&Czi);
1426 * MatGetRow(CTz_matrix,gi,&ncolumns,&gj,&CTzi);
1428 * MatGetRow(EntRes_matrix,gi,&ncolumns,&gj,&EntResi);
1429 * MatGetRow(SuppSize_matrix,gi,&ncolumns,&gj,&SuppSizei);
1430 * MatGetRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
1434 * get vector values
for column indices
1437 *
const std::vector<types::global_dof_index> gj_indices (gj,gj+ncolumns);
1438 * std::vector<double> soln(ncolumns);
1439 * std::vector<double> vx(ncolumns);
1440 * std::vector<double> vy(ncolumns);
1441 * std::vector<double> vz(ncolumns);
1442 * get_vector_values(solution,gj_indices,soln);
1443 * get_vector_values(locally_relevant_solution_vx,gj_indices,map_from_Q1_to_Q2,vx);
1444 * get_vector_values(locally_relevant_solution_vy,gj_indices,map_from_Q1_to_Q2,vy);
1446 * get_vector_values(locally_relevant_solution_vz,gj_indices,map_from_Q1_to_Q2,vz);
1450 * Array
for i-th row of matrices
1453 * std::vector<double> dLi(ncolumns), dCi(ncolumns);
1454 *
double dLii = 0, dCii = 0;
1457 *
loop on sparsity pattern of i-th DOF
1460 *
for (
int j =0; j < ncolumns; ++j)
1477 * ith_K_times_solution += soln[j]*(vj*
C);
1484 * low order dissipative
matrix
1491 * high order dissipative
matrix (entropy viscosity)
1495 * cE*
std::abs(EntResi[j])/(entropy_normalization_factor*MCi[j]/SuppSizei[j]));
1498 * high order compression
matrix
1502 * dCi[j] = dEij*
std::max(1-Compij,0.0);
1508 * save
K times solution vector
1509 * K_times_solution(gi)=ith_K_times_solution;
1510 * save i-th row of matrices on global matrices
1513 * MatSetValuesRow(dLij_matrix,gi,&dLi[0]);
1514 * dLij_matrix.set(gi,gi,dLii);
1515 * MatSetValuesRow(dCij_matrix,gi,&dCi[0]);
1516 * dCij_matrix.set(gi,gi,dCii);
1520 * Restore matrices after reading rows
1523 * MatRestoreRow(Cx_matrix,gi,&ncolumns,&gj,&Cxi);
1524 * MatRestoreRow(Cy_matrix,gi,&ncolumns,&gj,&Cyi);
1525 * MatRestoreRow(CTx_matrix,gi,&ncolumns,&gj,&CTxi);
1526 * MatRestoreRow(CTy_matrix,gi,&ncolumns,&gj,&CTyi);
1529 * MatRestoreRow(Cz_matrix,gi,&ncolumns,&gj,&Czi);
1530 * MatRestoreRow(CTz_matrix,gi,&ncolumns,&gj,&CTzi);
1532 * MatRestoreRow(EntRes_matrix,gi,&ncolumns,&gj,&EntResi);
1533 * MatRestoreRow(SuppSize_matrix,gi,&ncolumns,&gj,&SuppSizei);
1534 * MatRestoreRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
1546 * get matrices times vector
1549 * dLij_matrix.vmult(DL_times_solution,solution);
1550 * dCij_matrix.vmult(DH_times_solution,solution);
1555 * --------------------------------------------------------------------------------------
1556 * ------------------------------ ENTROPY VISCOSITY ------------------------------
1557 * --------------------------------------------------------------------------------------
1560 *
template <
int dim>
1561 *
void LevelSetSolver<dim>::assemble_EntRes_Matrix ()
1564 * entropy_normalization_factor=0;
1565 * SuppSize_matrix=0;
1567 *
const QGauss<dim> quadrature_formula(degree_MAX+1);
1577 *
const unsigned int dofs_per_cell_LS = fe_LS.dofs_per_cell;
1578 *
const unsigned int n_q_points = quadrature_formula.size();
1580 * std::vector<double> uqn (n_q_points);
1581 * std::vector<double> uqnm1 (n_q_points);
1582 * std::vector<Tensor<1,dim> > guqn (n_q_points);
1583 * std::vector<Tensor<1,dim> > guqnm1 (n_q_points);
1585 * std::vector<double> vxqn (n_q_points);
1586 * std::vector<double> vyqn (n_q_points);
1587 * std::vector<double> vzqn (n_q_points);
1588 * std::vector<double> vxqnm1 (n_q_points);
1589 * std::vector<double> vyqnm1 (n_q_points);
1590 * std::vector<double> vzqnm1 (n_q_points);
1595 * std::vector<Tensor<1, dim> > shape_grads_LS(dofs_per_cell_LS);
1596 * std::vector<double> shape_values_LS(dofs_per_cell_LS);
1598 * std::vector<types::global_dof_index> local_dof_indices_LS (dofs_per_cell_LS);
1601 * cell_LS = dof_handler_LS.begin_active();
1602 * endc_LS = dof_handler_LS.end();
1604 * cell_U = dof_handler_U.begin_active();
1607 *
double max_entropy=-1E10, min_entropy=1E10;
1608 *
double cell_max_entropy, cell_min_entropy;
1609 *
double cell_entropy_mass, entropy_mass=0;
1610 *
double cell_volume_double,
volume=0;
1612 *
for (; cell_LS!=endc_LS; ++cell_LS, ++cell_U)
1613 *
if (cell_LS->is_locally_owned())
1615 * cell_entropy_mass = 0;
1616 * cell_volume_double = 0;
1617 * cell_max_entropy = -1E10;
1618 * cell_min_entropy = 1E10;
1624 * get solutions at quadrature points
1627 * fe_values_LS.reinit(cell_LS);
1628 * cell_LS->get_dof_indices (local_dof_indices_LS);
1629 * fe_values_LS.get_function_values(un,uqn);
1630 * fe_values_LS.get_function_values(unm1,uqnm1);
1631 * fe_values_LS.get_function_gradients(un,guqn);
1632 * fe_values_LS.get_function_gradients(unm1,guqnm1);
1634 * fe_values_U.reinit(cell_U);
1635 * fe_values_U.get_function_values(locally_relevant_solution_vx,vxqn);
1636 * fe_values_U.get_function_values(locally_relevant_solution_vy,vyqn);
1637 *
if (dim==3) fe_values_U.get_function_values(locally_relevant_solution_vz,vzqn);
1638 * fe_values_U.get_function_values(locally_relevant_solution_vx_old,vxqnm1);
1639 * fe_values_U.get_function_values(locally_relevant_solution_vy_old,vyqnm1);
1640 *
if (dim==3) fe_values_U.get_function_values(locally_relevant_solution_vz_old,vzqnm1);
1642 *
for (
unsigned int q=0; q<n_q_points; ++q)
1644 * Rk = 1./time_step*(ENTROPY(uqn[q])-ENTROPY(uqnm1[q]))
1645 * +(vxqn[q]*ENTROPY_GRAD(uqn[q],guqn[q][0])+vyqn[q]*ENTROPY_GRAD(uqn[q],guqn[q][1]))/2.
1646 * +(vxqnm1[q]*ENTROPY_GRAD(uqnm1[q],guqnm1[q][0])+vyqnm1[q]*ENTROPY_GRAD(uqnm1[q],guqnm1[q][1]))/2.;
1648 * Rk += 0.5*(vzqn[q]*ENTROPY_GRAD(uqn[q],guqn[q][2])+vzqnm1[q]*ENTROPY_GRAD(uqnm1[q],guqnm1[q][2]));
1650 *
const double JxW = fe_values_LS.JxW(q);
1651 *
for (
unsigned int i=0; i<dofs_per_cell_LS; ++i)
1653 * shape_values_LS[i] = fe_values_LS.shape_value(i,q);
1654 * shape_grads_LS [i] = fe_values_LS.shape_grad (i,q);
1657 *
for (
unsigned int i=0; i<dofs_per_cell_LS; ++i)
1658 *
for (
unsigned int j=0; j < dofs_per_cell_LS; ++j)
1660 * cell_EntRes (i,j) += Rk*shape_values_LS[i]*shape_values_LS[j]*JxW;
1661 * cell_volume (i,j) += JxW;
1663 * cell_entropy_mass += ENTROPY(uqn[q])*JxW;
1664 * cell_volume_double += JxW;
1666 * cell_min_entropy =
std::min(cell_min_entropy,ENTROPY(uqn[q]));
1667 * cell_max_entropy =
std::max(cell_max_entropy,ENTROPY(uqn[q]));
1669 * entropy_mass += cell_entropy_mass;
1670 *
volume += cell_volume_double;
1672 * min_entropy =
std::min(min_entropy,cell_min_entropy);
1673 * max_entropy =
std::max(max_entropy,cell_max_entropy);
1679 * constraints.distribute_local_to_global(cell_EntRes,local_dof_indices_LS,EntRes_matrix);
1680 * constraints.distribute_local_to_global(cell_volume,local_dof_indices_LS,SuppSize_matrix);
1686 * ENTROPY NORM FACTOR
1698 * ------------------------------------------------------------------------------------
1699 * ------------------------------ TO CHECK MAX PRINCIPLE ------------------------------
1700 * ------------------------------------------------------------------------------------
1710 *
loop on locally owned i-DOFs (rows)
1714 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
1716 *
int gi = *idofs_iter;
1720 * get solution at DOFs on the sparsity pattern of i-th DOF
1723 * std::vector<types::global_dof_index> gj_indices = sparsity_pattern[gi];
1724 * std::vector<double> soln(gj_indices.size());
1725 * get_vector_values(un_solution,gj_indices,soln);
1728 * compute bounds, ith row of flux
matrix, P vectors
1731 *
double mini=1E10, maxi=-1E10;
1732 *
for (
unsigned int j =0; j < gj_indices.size(); ++j)
1742 * umin_vector(gi) = mini;
1743 * umax_vector(gi) = maxi;
1754 * compute
min and
max vectors
1757 *
const unsigned int dofs_per_cell = fe_LS.dofs_per_cell;
1758 * std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1762 * cell_LS = dof_handler_LS.begin_active(),
1763 * endc_LS = dof_handler_LS.end();
1765 *
for (; cell_LS!=endc_LS; ++cell_LS)
1766 *
if (cell_LS->is_locally_owned() && !cell_LS->at_boundary())
1768 * cell_LS->get_dof_indices(local_dof_indices);
1769 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1770 *
if (locally_owned_dofs_LS.is_element(local_dof_indices[i]))
1772 *
double solni = unp1_solution(local_dof_indices[i]);
1773 *
if (solni - umin_vector(local_dof_indices[i]) < -tol || umax_vector(local_dof_indices[i]) - solni < -tol)
1775 * pcout <<
"MAX Principle violated" << std::endl;
1784 * -------------------------------------------------------------------------------
1785 * ------------------------------ COMPUTE SOLUTIONS ------------------------------
1786 * -------------------------------------------------------------------------------
1790 *
void LevelSetSolver<dim>::compute_MPP_uL_and_NMPP_uH
1797 * NON-GHOSTED VECTORS: MPP_uL_solution, NMPP_uH_solution
1798 * GHOSTED VECTORS: un_solution
1801 * MPP_uL_solution=un_solution;
1802 * NMPP_uH_solution=un_solution;
1808 * assemble_K_times_vector(un_solution);
1809 * assemble_K_DL_DH_times_vector(un_solution);
1813 * COMPUTE MPP u1 solution
1817 * MPP_uL_solution.scale(ML_vector);
1818 * MPP_uL_solution.add(-time_step,K_times_solution);
1819 * MPP_uL_solution.add(-time_step,DL_times_solution);
1820 * MPP_uL_solution.scale(inverse_ML_vector);
1824 * COMPUTE GALERKIN u2 solution
1828 * MC_matrix.vmult(RHS,un_solution);
1829 * RHS.add(-time_step,K_times_solution,-time_step,DH_times_solution);
1830 * solve(constraints,MC_matrix,MC_preconditioner,NMPP_uH_solution,RHS);
1833 *
template <
int dim>
1834 *
void LevelSetSolver<dim>::compute_MPP_uH
1840 * MPP_uH_solution=0;
1843 *
loop on locally owned i-DOFs (rows)
1848 * PetscInt ncolumns;
1849 *
const PetscInt *gj;
1850 *
const PetscScalar *MCi, *dLi, *dCi;
1851 *
double solni, mi, solLi, solHi;
1853 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
1855 *
int gi = *idofs_iter;
1858 * read vectors at i-th DOF
1861 * solni=solution(gi);
1862 * solHi=NMPP_uH_solution_ghosted(gi);
1863 * solLi=MPP_uL_solution_ghosted(gi);
1868 * get i-th row of matrices
1871 * MatGetRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
1872 * MatGetRow(dLij_matrix,gi,&ncolumns,&gj,&dLi);
1873 * MatGetRow(dCij_matrix,gi,&ncolumns,&gj,&dCi);
1877 * get vector values
for support of i-th DOF
1880 *
const std::vector<types::global_dof_index> gj_indices (gj,gj+ncolumns);
1881 * std::vector<double> soln(ncolumns);
1882 * std::vector<double> solH(ncolumns);
1883 * get_vector_values(solution,gj_indices,soln);
1884 * get_vector_values(NMPP_uH_solution_ghosted,gj_indices,solH);
1888 * Array
for i-th row of matrices
1891 * std::vector<double> Ai(ncolumns);
1894 * compute bounds, ith row of flux
matrix, P vectors
1897 *
double mini=1E10, maxi=-1E10;
1898 *
double Pposi=0 ,Pnegi=0;
1899 *
for (
int j =0; j < ncolumns; ++j)
1911 * i-th row of flux
matrix A
1914 * Ai[j] = (((gi==gj[j]) ? 1 : 0)*mi - MCi[j])*(solH[j]-soln[j] - (solHi-solni))
1915 * +time_step*(dLi[j]-dCi[j])*(soln[j]-solni);
1922 * Pposi += Ai[j]*((Ai[j] > 0) ? 1. : 0.);
1923 * Pnegi += Ai[j]*((Ai[j] < 0) ? 1. : 0.);
1927 * save i-th row of flux
matrix A
1930 * MatSetValuesRow(A_matrix,gi,&Ai[0]);
1937 *
double Qposi = mi*(maxi-solLi);
1938 *
double Qnegi = mi*(mini-solLi);
1945 * R_pos_vector_nonGhosted(gi) = ((Pposi==0) ? 1. :
std::
min(1.0,Qposi/Pposi));
1946 * R_neg_vector_nonGhosted(gi) = ((Pnegi==0) ? 1. :
std::
min(1.0,Qnegi/Pnegi));
1950 * Restore matrices after reading rows
1953 * MatRestoreRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
1954 * MatRestoreRow(dLij_matrix,gi,&ncolumns,&gj,&dLi);
1955 * MatRestoreRow(dCij_matrix,gi,&ncolumns,&gj,&dCi);
1972 * update ghost values
for R vectors
1975 * R_pos_vector = R_pos_vector_nonGhosted;
1976 * R_neg_vector = R_neg_vector_nonGhosted;
1980 * compute limiters. NOTE:
this is a different
loop due to need of i- and j-th entries of R vectors
1984 *
double Rposi, Rnegi;
1985 * idofs_iter=locally_owned_dofs_LS.begin();
1986 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
1988 *
int gi = *idofs_iter;
1989 * Rposi = R_pos_vector(gi);
1990 * Rnegi = R_neg_vector(gi);
1994 * get i-th row of A
matrix
1997 * MatGetRow(A_matrix,gi,&ncolumns,&gj,&Ai);
2001 * get vector values
for column indices
2004 *
const std::vector<types::global_dof_index> gj_indices (gj,gj+ncolumns);
2005 * std::vector<double> Rpos(ncolumns);
2006 * std::vector<double> Rneg(ncolumns);
2007 * get_vector_values(R_pos_vector,gj_indices,Rpos);
2008 * get_vector_values(R_neg_vector,gj_indices,Rneg);
2012 * Array
for i-th row of A_times_L
matrix
2015 * std::vector<double> LxAi(ncolumns);
2018 *
loop in sparsity pattern of i-th DOF
2021 *
for (
int j =0; j < ncolumns; ++j)
2022 * LxAi[j] = Ai[j] * ((Ai[j]>0) ?
std::min(Rposi,Rneg[j]) :
std::
min(Rnegi,Rpos[j]));
2026 * save i-th row of LxA
2029 * MatSetValuesRow(LxA_matrix,gi,&LxAi[0]);
2032 * restore A
matrix after reading it
2035 * MatRestoreRow(A_matrix,gi,&ncolumns,&gj,&Ai);
2038 * LxA_matrix.vmult(MPP_uH_solution,ones_vector);
2039 * MPP_uH_solution.scale(inverse_ML_vector);
2040 * MPP_uH_solution.add(1.0,MPP_uL_solution_ghosted);
2044 *
void LevelSetSolver<dim>::compute_MPP_uH_with_iterated_FCT
2050 * MPP_uH_solution=0;
2051 * compute_MPP_uH(MPP_uH_solution,MPP_uL_solution_ghosted,NMPP_uH_solution_ghosted,un_solution);
2055 * Akp1_matrix.copy_from(A_matrix);
2056 * LxAkp1_matrix.copy_from(LxA_matrix);
2060 *
loop in num of FCT iterations
2063 * PetscInt ncolumns;
2064 *
const PetscInt *gj;
2065 *
const PetscScalar *Akp1i;
2067 *
for (
int iter=0; iter<NUM_ITER; ++iter)
2069 * MPP_uLkp1_solution_ghosted = MPP_uH_solution;
2070 * Akp1_matrix.add(-1.0, LxAkp1_matrix);
2074 *
loop on locally owned i-DOFs (rows)
2078 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
2080 *
int gi = *idofs_iter;
2084 * read vectors at i-th DOF
2088 *
double solLi = MPP_uLkp1_solution_ghosted(gi);
2092 * get i-th row of matrices
2095 * MatGetRow(Akp1_matrix,gi,&ncolumns,&gj,&Akp1i);
2098 * get vector values
for support of i-th DOF
2101 *
const std::vector<types::global_dof_index> gj_indices (gj,gj+ncolumns);
2102 * std::vector<double> soln(ncolumns);
2103 * get_vector_values(un_solution,gj_indices,soln);
2107 * compute bounds, ith row of flux
matrix, P vectors
2110 *
double mini=1E10, maxi=-1E10;
2111 *
double Pposi=0 ,Pnegi=0;
2112 *
for (
int j =0; j < ncolumns; ++j)
2127 * Pposi += Akp1i[j]*((Akp1i[j] > 0) ? 1. : 0.);
2128 * Pnegi += Akp1i[j]*((Akp1i[j] < 0) ? 1. : 0.);
2135 *
double Qposi = mi*(maxi-solLi);
2136 *
double Qnegi = mi*(mini-solLi);
2143 * R_pos_vector_nonGhosted(gi) = ((Pposi==0) ? 1. :
std::
min(1.0,Qposi/Pposi));
2144 * R_neg_vector_nonGhosted(gi) = ((Pnegi==0) ? 1. :
std::
min(1.0,Qnegi/Pnegi));
2148 * Restore matrices after reading rows
2151 * MatRestoreRow(Akp1_matrix,gi,&ncolumns,&gj,&Akp1i);
2162 * update ghost values
for R vectors
2165 * R_pos_vector = R_pos_vector_nonGhosted;
2166 * R_neg_vector = R_neg_vector_nonGhosted;
2170 * compute limiters. NOTE:
this is a different
loop due to need of i- and j-th entries of R vectors
2173 *
double Rposi, Rnegi;
2174 * idofs_iter=locally_owned_dofs_LS.begin();
2175 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
2177 *
int gi = *idofs_iter;
2178 * Rposi = R_pos_vector(gi);
2179 * Rnegi = R_neg_vector(gi);
2183 * get i-th row of Akp1
matrix
2186 * MatGetRow(Akp1_matrix,gi,&ncolumns,&gj,&Akp1i);
2190 * get vector values
for column indices
2193 *
const std::vector<types::global_dof_index> gj_indices(gj,gj+ncolumns);
2194 * std::vector<double> Rpos(ncolumns);
2195 * std::vector<double> Rneg(ncolumns);
2196 * get_vector_values(R_pos_vector,gj_indices,Rpos);
2197 * get_vector_values(R_neg_vector,gj_indices,Rneg);
2201 * Array
for i-th row of LxAkp1
matrix
2204 * std::vector<double> LxAkp1i(ncolumns);
2205 *
for (
int j =0; j < ncolumns; ++j)
2206 * LxAkp1i[j] = Akp1i[j] * ((Akp1i[j]>0) ?
std::min(Rposi,Rneg[j]) :
std::
min(Rnegi,Rpos[j]));
2210 * save i-th row of LxA
2213 * MatSetValuesRow(LxAkp1_matrix,gi,&LxAkp1i[0]);
2216 * restore A
matrix after reading it
2219 * MatRestoreRow(Akp1_matrix,gi,&ncolumns,&gj,&Akp1i);
2222 * LxAkp1_matrix.vmult(MPP_uH_solution,ones_vector);
2223 * MPP_uH_solution.scale(inverse_ML_vector);
2224 * MPP_uH_solution.add(1.0,MPP_uLkp1_solution_ghosted);
2232 * std::string algorithm)
2237 * COMPUTE MPP LOW-ORDER SOLN and NMPP HIGH-ORDER SOLN
2240 * compute_MPP_uL_and_NMPP_uH(MPP_uL_solution,NMPP_uH_solution,un);
2242 *
if (algorithm.compare(
"MPP_u1")==0)
2243 * unp1=MPP_uL_solution;
2244 *
else if (algorithm.compare(
"NMPP_uH")==0)
2245 * unp1=NMPP_uH_solution;
2246 *
else if (algorithm.compare(
"MPP_uH")==0)
2248 * MPP_uL_solution_ghosted = MPP_uL_solution;
2249 * NMPP_uH_solution_ghosted=NMPP_uH_solution;
2250 * compute_MPP_uH_with_iterated_FCT(MPP_uH_solution,MPP_uL_solution_ghosted,NMPP_uH_solution_ghosted,un);
2251 * unp1=MPP_uH_solution;
2255 * pcout <<
"Error in algorithm" << std::endl;
2263 * std::string algorithm)
2267 * GHOSTED VECTORS: un
2268 * NON-GHOSTED VECTORS: unp1
2272 * uStage1=0., uStage2=0.;
2273 * uStage1_nonGhosted=0., uStage2_nonGhosted=0.;
2282 * compute_solution(uStage1_nonGhosted,un,algorithm);
2283 * uStage1=uStage1_nonGhosted;
2289 * u2=3/4*un+1/4*(u1-dt*RH*u1)
2292 * compute_solution(uStage2_nonGhosted,uStage1,algorithm);
2293 * uStage2_nonGhosted*=1./4;
2294 * uStage2_nonGhosted.add(3./4,un);
2295 * uStage2=uStage2_nonGhosted;
2301 * unp1=1/3*un+2/3*(u2-dt*RH*u2)
2304 * compute_solution(unp1,uStage2,algorithm);
2306 * unp1.add(1./3,un);
2311 * -----------------------------------------------------------------------
2312 * ------------------------------ UTILITIES ------------------------------
2313 * -----------------------------------------------------------------------
2317 *
void LevelSetSolver<dim>::get_sparsity_pattern()
2325 * PetscInt ncolumns;
2326 *
const PetscInt *gj;
2327 *
const PetscScalar *MCi;
2329 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
2331 * PetscInt gi = *idofs_iter;
2334 * get i-th row of mass
matrix (dummy, I just need the indices gj)
2337 * MatGetRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
2338 * sparsity_pattern[gi] = std::vector<types::global_dof_index>(gj,gj+ncolumns);
2339 * MatRestoreRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
2344 *
void LevelSetSolver<dim>::get_map_from_Q1_to_Q2()
2346 * map_from_Q1_to_Q2.clear();
2347 *
const unsigned int dofs_per_cell_LS = fe_LS.dofs_per_cell;
2348 * std::vector<types::global_dof_index> local_dof_indices_LS (dofs_per_cell_LS);
2349 *
const unsigned int dofs_per_cell_U = fe_U.dofs_per_cell;
2350 * std::vector<types::global_dof_index> local_dof_indices_U (dofs_per_cell_U);
2353 * cell_LS = dof_handler_LS.begin_active(),
2354 * endc_LS = dof_handler_LS.end();
2356 * cell_U = dof_handler_U.begin_active();
2358 *
for (; cell_LS!=endc_LS; ++cell_LS, ++cell_U)
2359 *
if (!cell_LS->is_artificial())
2361 * cell_LS->get_dof_indices(local_dof_indices_LS);
2362 * cell_U->get_dof_indices(local_dof_indices_U);
2363 *
for (
unsigned int i=0; i<dofs_per_cell_LS; ++i)
2364 * map_from_Q1_to_Q2[local_dof_indices_LS[i]] = local_dof_indices_U[i];
2368 *
template <
int dim>
2371 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
2377 * all vectors are NON-GHOSTED
2380 *
SolverControl solver_control (dof_handler_LS.n_dofs(), solver_tolerance);
2382 * constraints.distribute (completely_distributed_solution);
2383 * solver.solve (Matrix, completely_distributed_solution, rhs, *preconditioner);
2384 * constraints.distribute (completely_distributed_solution);
2385 *
if (verbose==
true) pcout <<
" Solved in " << solver_control.last_step() <<
" iterations." << std::endl;
2388 *
template <
int dim>
2389 *
void LevelSetSolver<dim>::save_old_solution()
2395 *
template <
int dim>
2396 *
void LevelSetSolver<dim>::save_old_vel_solution()
2398 * locally_relevant_solution_vx_old = locally_relevant_solution_vx;
2399 * locally_relevant_solution_vy_old = locally_relevant_solution_vy;
2401 * locally_relevant_solution_vz_old = locally_relevant_solution_vz;
2406 * -------------------------------------------------------------------------------
2407 * ------------------------------ MY PETSC WRAPPERS ------------------------------
2408 * -------------------------------------------------------------------------------
2413 *
const std::vector<types::global_dof_index> &indices,
2414 * std::vector<PetscScalar> &values)
2418 *
PETSc wrapper to get sets of values from a petsc vector.
2419 * we assume the vector is ghosted
2420 * We need to figure out which elements we
2421 * own locally. Then get a pointer to the
2422 * elements that are stored here (both the
2423 * ones we own as well as the ghost elements).
2424 * In
this array, the locally owned elements
2425 * come
first followed by the ghost elements whose
2426 * position we can get from an
index set
2432 *
IndexSet ghost_indices = locally_relevant_dofs_LS;
2436 * n_idx = indices.size();
2438 * VecGetOwnershipRange (vector, &begin, &end);
2440 * Vec solution_in_local_form =
nullptr;
2441 * VecGhostGetLocalForm(vector, &solution_in_local_form);
2443 * PetscScalar *soln;
2444 * VecGetArray(solution_in_local_form, &soln);
2446 *
for (i = 0; i < n_idx; ++i)
2448 *
int index = indices[i];
2449 *
if (index >= begin && index < end)
2453 *
const unsigned int ghostidx = ghost_indices.index_within_set(index);
2454 * values[i] = *(soln+ghostidx+
end-
begin);
2457 * VecRestoreArray(solution_in_local_form, &soln);
2458 * VecGhostRestoreLocalForm(vector, &solution_in_local_form);
2463 *
const std::vector<types::global_dof_index> &indices,
2464 * std::map<types::global_dof_index, types::global_dof_index> &map_from_Q1_to_Q2,
2465 * std::vector<PetscScalar> &values)
2469 * THIS IS MEANT TO BE USED WITH VELOCITY VECTORS
2470 *
PETSc wrapper to get sets of values from a petsc vector.
2471 * we assume the vector is ghosted
2472 * We need to figure out which elements we
2473 * own locally. Then get a pointer to the
2474 * elements that are stored here (both the
2475 * ones we own as well as the ghost elements).
2476 * In
this array, the locally owned elements
2477 * come
first followed by the ghost elements whose
2478 * position we can get from an
index set
2484 *
IndexSet ghost_indices = locally_relevant_dofs_U;
2488 * n_idx = indices.size();
2490 * VecGetOwnershipRange (vector, &begin, &end);
2492 * Vec solution_in_local_form =
nullptr;
2493 * VecGhostGetLocalForm(vector, &solution_in_local_form);
2495 * PetscScalar *soln;
2496 * VecGetArray(solution_in_local_form, &soln);
2498 *
for (i = 0; i < n_idx; ++i)
2500 *
int index = map_from_Q1_to_Q2[indices[i]];
2501 *
if (index >= begin && index < end)
2505 *
const unsigned int ghostidx = ghost_indices.index_within_set(index);
2506 * values[i] = *(soln+ghostidx+
end-
begin);
2509 * VecRestoreArray(solution_in_local_form, &soln);
2510 * VecGhostRestoreLocalForm(vector, &solution_in_local_form);
2516<a name=
"ann-MultiPhase.cc"></a>
2517<h1>Annotated version of MultiPhase.cc</h1>
2533 * #include <deal.II/base/quadrature_lib.h>
2534 * #include <deal.II/base/function.h>
2535 * #include <deal.II/lac/affine_constraints.h>
2536 * #include <deal.II/lac/vector.h>
2537 * #include <deal.II/lac/full_matrix.h>
2538 * #include <deal.II/lac/solver_cg.h>
2539 * #include <deal.II/lac/petsc_sparse_matrix.h>
2540 * #include <deal.II/lac/petsc_vector.h>
2541 * #include <deal.II/lac/petsc_solver.h>
2542 * #include <deal.II/lac/petsc_precondition.h>
2543 * #include <deal.II/grid/grid_generator.h>
2544 * #include <deal.II/grid/tria_accessor.h>
2545 * #include <deal.II/grid/tria_iterator.h>
2546 * #include <deal.II/dofs/dof_handler.h>
2547 * #include <deal.II/dofs/dof_accessor.h>
2548 * #include <deal.II/dofs/dof_tools.h>
2549 * #include <deal.II/fe/fe_values.h>
2550 * #include <deal.II/fe/fe_q.h>
2551 * #include <deal.II/numerics/vector_tools.h>
2552 * #include <deal.II/numerics/data_out.h>
2553 * #include <deal.II/numerics/error_estimator.h>
2554 * #include <deal.II/base/utilities.h>
2555 * #include <deal.II/base/conditional_ostream.h>
2556 * #include <deal.II/base/index_set.h>
2557 * #include <deal.II/lac/sparsity_tools.h>
2558 * #include <deal.II/distributed/tria.h>
2559 * #include <deal.II/distributed/grid_refinement.h>
2560 * #include <deal.II/base/convergence_table.h>
2561 * #include <deal.II/base/timer.h>
2562 * #include <deal.II/base/parameter_handler.h>
2563 * #include <fstream>
2564 * #include <iostream>
2565 * #include <deal.II/grid/grid_tools.h>
2566 * #include <deal.II/fe/mapping_q.h>
2568 *
using namespace dealii;
2573 * FOR TRANSPORT PROBLEM
2585 * #define FILLING_TANK 0
2586 * #define BREAKING_DAM 1
2587 * #define FALLING_DROP 2
2588 * #define SMALL_WAVE_PERTURBATION 3
2590 * #include
"NavierStokesSolver.cc"
2591 * #include
"LevelSetSolver.cc"
2592 * #include
"utilities.cc"
2601 *
template <
int dim>
2605 * MultiPhase (
const unsigned int degree_LS,
2606 *
const unsigned int degree_U);
2611 *
void set_boundary_inlet();
2612 *
void get_boundary_values_U();
2613 *
void get_boundary_values_phi(std::vector<types::global_dof_index> &boundary_values_id_phi,
2614 * std::vector<double> &boundary_values_phi);
2615 *
void output_results();
2616 *
void output_vectors();
2617 *
void output_rho();
2619 *
void initial_condition();
2620 *
void init_constraints();
2629 *
IndexSet locally_relevant_dofs_LS;
2636 *
IndexSet locally_relevant_dofs_U;
2641 *
IndexSet locally_relevant_dofs_P;
2663 * std::vector<types::global_dof_index> boundary_values_id_u;
2664 * std::vector<types::global_dof_index> boundary_values_id_v;
2665 * std::vector<types::global_dof_index> boundary_values_id_phi;
2666 * std::vector<double> boundary_values_u;
2667 * std::vector<double> boundary_values_v;
2668 * std::vector<double> boundary_values_phi;
2674 *
double final_time;
2675 *
unsigned int timestep_number;
2681 *
int sharpness_integer;
2683 *
unsigned int n_refinement;
2684 *
unsigned int output_number;
2685 *
double output_time;
2709 *
unsigned int TRANSPORT_TIME_INTEGRATION;
2710 * std::string ALGORITHM;
2711 *
unsigned int PROBLEM;
2714 *
template <
int dim>
2715 * MultiPhase<dim>::MultiPhase (
const unsigned int degree_LS,
2716 *
const unsigned int degree_U)
2718 * mpi_communicator (MPI_COMM_WORLD),
2723 * degree_LS(degree_LS),
2725 * fe_LS (degree_LS),
2726 * degree_U(degree_U),
2730 * fe_P (degree_U-1),
2734 *
template <
int dim>
2735 * MultiPhase<dim>::~MultiPhase ()
2737 * dof_handler_LS.
clear ();
2738 * dof_handler_U.clear ();
2739 * dof_handler_P.clear ();
2749 *
template <
int dim>
2750 *
void MultiPhase<dim>::setup()
2757 * dof_handler_LS.distribute_dofs (fe_LS);
2758 * locally_owned_dofs_LS = dof_handler_LS.locally_owned_dofs ();
2765 * dof_handler_U.distribute_dofs (fe_U);
2766 * locally_owned_dofs_U = dof_handler_U.locally_owned_dofs ();
2773 * dof_handler_P.distribute_dofs (fe_P);
2774 * locally_owned_dofs_P = dof_handler_P.locally_owned_dofs ();
2778 * init vectors
for phi
2781 * locally_relevant_solution_phi.reinit(locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
2782 * locally_relevant_solution_phi = 0;
2783 * completely_distributed_solution_phi.reinit (locally_owned_dofs_P,mpi_communicator);
2786 * init vectors
for u
2789 * locally_relevant_solution_u.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
2790 * locally_relevant_solution_u = 0;
2791 * completely_distributed_solution_u.reinit (locally_owned_dofs_U,mpi_communicator);
2794 * init vectors
for v
2797 * locally_relevant_solution_v.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
2798 * locally_relevant_solution_v = 0;
2799 * completely_distributed_solution_v.reinit (locally_owned_dofs_U,mpi_communicator);
2802 * init vectors
for p
2805 * locally_relevant_solution_p.reinit (locally_owned_dofs_P,locally_relevant_dofs_P,mpi_communicator);
2806 * locally_relevant_solution_p = 0;
2807 * completely_distributed_solution_p.reinit (locally_owned_dofs_P,mpi_communicator);
2813 * init_constraints();
2816 *
template <
int dim>
2817 *
void MultiPhase<dim>::initial_condition()
2822 * Initial conditions
2823 * init condition
for phi
2826 * completely_distributed_solution_phi = 0;
2828 * InitialPhi<dim>(PROBLEM, sharpness),
2829 * completely_distributed_solution_phi);
2830 * constraints.distribute (completely_distributed_solution_phi);
2831 * locally_relevant_solution_phi = completely_distributed_solution_phi;
2834 * init condition
for u=0
2837 * completely_distributed_solution_u = 0;
2840 * completely_distributed_solution_u);
2841 * constraints.distribute (completely_distributed_solution_u);
2842 * locally_relevant_solution_u = completely_distributed_solution_u;
2845 * init condition
for v
2848 * completely_distributed_solution_v = 0;
2851 * completely_distributed_solution_v);
2852 * constraints.distribute (completely_distributed_solution_v);
2853 * locally_relevant_solution_v = completely_distributed_solution_v;
2856 * init condition
for p
2859 * completely_distributed_solution_p = 0;
2862 * completely_distributed_solution_p);
2863 * constraints.distribute (completely_distributed_solution_p);
2864 * locally_relevant_solution_p = completely_distributed_solution_p;
2867 *
template <
int dim>
2868 *
void MultiPhase<dim>::init_constraints()
2870 * constraints.clear ();
2871 * constraints.reinit (locally_relevant_dofs_LS);
2873 * constraints.close ();
2876 *
template <
int dim>
2877 *
void MultiPhase<dim>::get_boundary_values_U()
2879 * std::map<types::global_dof_index, double> map_boundary_values_u;
2880 * std::map<types::global_dof_index, double> map_boundary_values_v;
2881 * std::map<types::global_dof_index, double> map_boundary_values_w;
2888 *
if (PROBLEM==BREAKING_DAM || PROBLEM==FALLING_DROP)
2919 *
else if (PROBLEM==SMALL_WAVE_PERTURBATION)
2923 * no slip in bottom and top and slip in left and right
2949 *
else if (PROBLEM==FILLING_TANK)
2953 * LEFT: entry in x, zero in y
2960 * RIGHT: no-slip condition
2974 * TOP: exit in y, zero in x
2982 * pcout <<
"Error in type of PROBLEM at Boundary Conditions" << std::endl;
2985 * boundary_values_id_u.resize(map_boundary_values_u.size());
2986 * boundary_values_id_v.resize(map_boundary_values_v.size());
2987 * boundary_values_u.resize(map_boundary_values_u.size());
2988 * boundary_values_v.resize(map_boundary_values_v.size());
2989 * std::map<types::global_dof_index,double>::const_iterator boundary_value_u =map_boundary_values_u.begin();
2990 * std::map<types::global_dof_index,double>::const_iterator boundary_value_v =map_boundary_values_v.begin();
2992 *
for (
int i=0; boundary_value_u !=map_boundary_values_u.end(); ++boundary_value_u, ++i)
2994 * boundary_values_id_u[i]=boundary_value_u->first;
2995 * boundary_values_u[i]=boundary_value_u->second;
2997 *
for (
int i=0; boundary_value_v !=map_boundary_values_v.end(); ++boundary_value_v, ++i)
2999 * boundary_values_id_v[i]=boundary_value_v->first;
3000 * boundary_values_v[i]=boundary_value_v->second;
3004 *
template <
int dim>
3005 *
void MultiPhase<dim>::set_boundary_inlet()
3007 *
const QGauss<dim-1> face_quadrature_formula(1);
3011 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
3012 * std::vector<double> u_value (n_face_q_points);
3013 * std::vector<double> v_value (n_face_q_points);
3016 * cell_U = dof_handler_U.begin_active(),
3017 * endc_U = dof_handler_U.end();
3020 *
for (; cell_U!=endc_U; ++cell_U)
3021 *
if (cell_U->is_locally_owned())
3022 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
3023 *
if (cell_U->face(face)->at_boundary())
3025 * fe_face_values.reinit(cell_U,face);
3026 * fe_face_values.get_function_values(locally_relevant_solution_u,u_value);
3027 * fe_face_values.get_function_values(locally_relevant_solution_v,v_value);
3030 *
if (fe_face_values.normal_vector(0)*u < -1e-14)
3031 * cell_U->face(face)->set_boundary_id(10);
3035 *
template <
int dim>
3036 *
void MultiPhase<dim>::get_boundary_values_phi(std::vector<types::global_dof_index> &boundary_values_id_phi,
3037 * std::vector<double> &boundary_values_phi)
3039 * std::map<types::global_dof_index, double> map_boundary_values_phi;
3042 * set_boundary_inlet();
3045 * boundary_values_id_phi.resize(map_boundary_values_phi.size());
3046 * boundary_values_phi.resize(map_boundary_values_phi.size());
3047 * std::map<types::global_dof_index,double>::const_iterator boundary_value_phi = map_boundary_values_phi.begin();
3048 *
for (
int i=0; boundary_value_phi !=map_boundary_values_phi.end(); ++boundary_value_phi, ++i)
3050 * boundary_values_id_phi[i]=boundary_value_phi->first;
3051 * boundary_values_phi[i]=boundary_value_phi->second;
3056 *
void MultiPhase<dim>::output_results()
3067 *
template <
int dim>
3068 *
void MultiPhase<dim>::output_vectors()
3072 * data_out.add_data_vector (locally_relevant_solution_phi,
"phi");
3073 * data_out.build_patches ();
3075 *
const std::string filename = (
"sol_vectors-" +
3080 * std::ofstream output ((filename +
".vtu").c_str());
3081 * data_out.write_vtu (output);
3085 * std::vector<std::string> filenames;
3086 *
for (
unsigned int i=0;
3089 * filenames.push_back (
"sol_vectors-" +
3095 * std::ofstream master_output ((filename +
".pvtu").c_str());
3096 * data_out.write_pvtu_record (master_output, filenames);
3100 *
template <
int dim>
3101 *
void MultiPhase<dim>::output_rho()
3103 * Postprocessor<dim> postprocessor(eps,rho_air,rho_fluid);
3106 * data_out.add_data_vector (locally_relevant_solution_phi, postprocessor);
3108 * data_out.build_patches ();
3110 *
const std::string filename = (
"sol_rho-" +
3115 * std::ofstream output ((filename +
".vtu").c_str());
3116 * data_out.write_vtu (output);
3120 * std::vector<std::string> filenames;
3121 *
for (
unsigned int i=0;
3124 * filenames.push_back (
"sol_rho-" +
3130 * std::ofstream master_output ((filename +
".pvtu").c_str());
3131 * data_out.write_pvtu_record (master_output, filenames);
3135 *
template <
int dim>
3136 *
void MultiPhase<dim>::run()
3141 * GENERAL PARAMETERS
3148 * get_output =
true;
3149 * output_number = 0;
3151 * output_time = 0.1;
3152 * final_time = 10.0;
3156 * PARAMETERS FOR THE NAVIER STOKES PROBLEM
3160 * rho_fluid = 1000.;
3164 * PROBLEM=BREAKING_DAM;
3167 * PROBLEM=FILLING_TANK;
3168 * PROBLEM=SMALL_WAVE_PERTURBATION;
3169 * PROBLEM=FALLING_DROP;
3175 * ForceTerms<dim> force_function(std::vector<double> {0.0,-1.0});
3179 * PARAMETERS FOR TRANSPORT PROBLEM
3185 * sharpness_integer=10;
3191 * TRANSPORT_TIME_INTEGRATION=SSP33;
3194 * ALGORITHM =
"MPP_u1";
3195 * ALGORITHM =
"NMPP_uH";
3198 * ALGORITHM =
"MPP_uH";
3202 * ADJUST PARAMETERS ACCORDING TO PROBLEM
3205 *
if (PROBLEM==FALLING_DROP)
3215 *
if (PROBLEM==FILLING_TANK)
3218 *
else if (PROBLEM==BREAKING_DAM || PROBLEM==SMALL_WAVE_PERTURBATION)
3220 * std::vector< unsigned int > repetitions;
3221 * repetitions.push_back(2);
3222 * repetitions.push_back(1);
3226 *
else if (PROBLEM==FALLING_DROP)
3228 * std::vector< unsigned int > repetitions;
3229 * repetitions.push_back(1);
3230 * repetitions.push_back(4);
3244 * PARAMETERS FOR TIME STEPPING
3248 * time_step = cfl*min_h/umax;
3250 * sharpness=sharpness_integer*min_h;
3254 * INITIAL CONDITIONS
3257 * initial_condition();
3262 * NAVIER STOKES SOLVER
3265 * NavierStokesSolver<dim> navier_stokes (degree_LS,degree_U,
3268 * rho_fluid,nu_fluid,
3274 * BOUNDARY CONDITIONS FOR NAVIER STOKES
3277 * get_boundary_values_U();
3278 * navier_stokes.set_boundary_conditions(boundary_values_id_u, boundary_values_id_v,
3279 * boundary_values_u, boundary_values_v);
3283 *
set INITIAL CONDITION within NAVIER STOKES
3286 * navier_stokes.initial_condition(locally_relevant_solution_phi,
3287 * locally_relevant_solution_u,
3288 * locally_relevant_solution_v,
3289 * locally_relevant_solution_p);
3295 * LevelSetSolver<dim> transport_solver (degree_LS,degree_U,
3299 * TRANSPORT_TIME_INTEGRATION,
3301 * mpi_communicator);
3304 * BOUNDARY CONDITIONS FOR PHI
3307 * get_boundary_values_phi(boundary_values_id_phi,boundary_values_phi);
3308 * transport_solver.set_boundary_conditions(boundary_values_id_phi,boundary_values_phi);
3312 *
set INITIAL CONDITION within TRANSPORT PROBLEM
3315 * transport_solver.initial_condition(locally_relevant_solution_phi,
3316 * locally_relevant_solution_u,
3317 * locally_relevant_solution_v);
3318 *
int dofs_U = 2*dof_handler_U.n_dofs();
3319 *
int dofs_P = 2*dof_handler_P.n_dofs();
3320 *
int dofs_LS = dof_handler_LS.n_dofs();
3321 *
int dofs_TOTAL = dofs_U+dofs_P+dofs_LS;
3325 * NO BOUNDARY CONDITIONS
for LEVEL SET
3328 * pcout <<
"Cfl: " << cfl <<
"; umax: " << umax <<
"; min h: " << min_h
3329 * <<
"; time step: " << time_step << std::endl;
3330 * pcout <<
" Number of active cells: "
3332 * <<
" Number of degrees of freedom: " << std::endl
3333 * <<
" U: " << dofs_U << std::endl
3334 * <<
" P: " << dofs_P << std::endl
3335 * <<
" LS: " << dofs_LS << std::endl
3336 * <<
" TOTAL: " << dofs_TOTAL
3344 *
for (timestep_number=1, time=time_step; time<=final_time;
3345 * time+=time_step,++timestep_number)
3347 * pcout <<
"Time step " << timestep_number
3348 * <<
" at t=" << time
3352 * GET NAVIER STOKES VELOCITY
3355 * navier_stokes.set_phi(locally_relevant_solution_phi);
3356 * navier_stokes.nth_time_step();
3357 * navier_stokes.get_velocity(locally_relevant_solution_u,locally_relevant_solution_v);
3358 * transport_solver.set_velocity(locally_relevant_solution_u,locally_relevant_solution_v);
3361 * GET LEVEL SET SOLUTION
3364 * transport_solver.nth_time_step();
3365 * transport_solver.get_unp1(locally_relevant_solution_phi);
3366 *
if (get_output && time-(output_number)*output_time>0)
3369 * navier_stokes.get_velocity(locally_relevant_solution_u, locally_relevant_solution_v);
3370 * transport_solver.get_unp1(locally_relevant_solution_phi);
3375 *
int main(
int argc,
char *argv[])
3379 *
using namespace dealii;
3381 * PetscInitialize(&argc, &argv,
nullptr,
nullptr);
3384 *
unsigned int degree_LS = 1;
3385 *
unsigned int degree_U = 2;
3386 * MultiPhase<2> multi_phase(degree_LS, degree_U);
3387 * multi_phase.run();
3392 *
catch (std::exception &exc)
3394 * std::cerr << std::endl << std::endl
3395 * <<
"----------------------------------------------------"
3397 * std::cerr <<
"Exception on processing: " << std::endl
3398 * << exc.what() << std::endl
3399 * <<
"Aborting!" << std::endl
3400 * <<
"----------------------------------------------------"
3406 * std::cerr << std::endl << std::endl
3407 * <<
"----------------------------------------------------"
3409 * std::cerr <<
"Unknown exception!" << std::endl
3410 * <<
"Aborting!" << std::endl
3411 * <<
"----------------------------------------------------"
3420<a name=
"ann-NavierStokesSolver.cc"></a>
3421<h1>Annotated version of NavierStokesSolver.cc</h1>
3437 * #include <deal.II/base/quadrature_lib.h>
3438 * #include <deal.II/base/function.h>
3439 * #include <deal.II/lac/affine_constraints.h>
3440 * #include <deal.II/lac/vector.h>
3441 * #include <deal.II/lac/full_matrix.h>
3442 * #include <deal.II/lac/solver_cg.h>
3443 * #include <deal.II/lac/petsc_sparse_matrix.h>
3444 * #include <deal.II/lac/petsc_vector.h>
3445 * #include <deal.II/lac/petsc_solver.h>
3446 * #include <deal.II/lac/petsc_precondition.h>
3447 * #include <deal.II/grid/grid_generator.h>
3448 * #include <deal.II/grid/tria_accessor.h>
3449 * #include <deal.II/grid/tria_iterator.h>
3450 * #include <deal.II/dofs/dof_handler.h>
3451 * #include <deal.II/dofs/dof_accessor.h>
3452 * #include <deal.II/dofs/dof_tools.h>
3453 * #include <deal.II/fe/fe_values.h>
3454 * #include <deal.II/fe/fe_q.h>
3455 * #include <deal.II/numerics/vector_tools.h>
3456 * #include <deal.II/numerics/data_out.h>
3457 * #include <deal.II/numerics/error_estimator.h>
3458 * #include <deal.II/base/utilities.h>
3459 * #include <deal.II/base/conditional_ostream.h>
3460 * #include <deal.II/base/index_set.h>
3461 * #include <deal.II/lac/sparsity_tools.h>
3462 * #include <deal.II/distributed/tria.h>
3463 * #include <deal.II/distributed/grid_refinement.h>
3464 * #include <deal.II/lac/petsc_vector.h>
3465 * #include <deal.II/base/convergence_table.h>
3466 * #include <deal.II/base/timer.h>
3467 * #include <deal.II/base/parameter_handler.h>
3468 * #include <deal.II/grid/grid_tools.h>
3469 * #include <deal.II/fe/mapping_q.h>
3471 * #include <fstream>
3472 * #include <iostream>
3475 *
using namespace dealii;
3477 * #define MAX_NUM_ITER_TO_RECOMPUTE_PRECONDITIONER 10
3487 *
class NavierStokesSolver
3492 * constructor
for using LEVEL SET
3495 * NavierStokesSolver(
const unsigned int degree_LS,
3496 *
const unsigned int degree_U,
3497 *
const double time_step,
3499 *
const double rho_air,
3500 *
const double nu_air,
3501 *
const double rho_fluid,
3502 *
const double nu_fluid,
3504 *
const bool verbose,
3509 * constructor
for NOT LEVEL SET
3512 * NavierStokesSolver(
const unsigned int degree_LS,
3513 *
const unsigned int degree_U,
3514 *
const double time_step,
3518 *
const bool verbose,
3527 *
void set_rho_and_nu_functions(
const Function<dim> &rho_function,
3545 * boundary conditions
3548 *
void set_boundary_conditions(std::vector<types::global_dof_index> boundary_values_id_u,
3549 * std::vector<types::global_dof_index> boundary_values_id_v, std::vector<double> boundary_values_u,
3550 * std::vector<double> boundary_values_v);
3551 *
void set_boundary_conditions(std::vector<types::global_dof_index> boundary_values_id_u,
3552 * std::vector<types::global_dof_index> boundary_values_id_v,
3553 * std::vector<types::global_dof_index> boundary_values_id_w, std::vector<double> boundary_values_u,
3554 * std::vector<double> boundary_values_v, std::vector<double> boundary_values_w);
3572 *
void nth_time_step();
3580 * ~NavierStokesSolver();
3585 * SETUP AND INITIAL CONDITION
3589 *
void setup_VECTORS();
3590 *
void init_constraints();
3596 *
void assemble_system_U();
3597 *
void assemble_system_dpsi_q();
3604 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
3608 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
3613 * GET DIFFERENT FIELDS
3616 *
void get_rho_and_nu(
double phi);
3617 *
void get_velocity();
3618 *
void get_pressure();
3624 *
void save_old_solution();
3633 *
IndexSet locally_relevant_dofs_LS;
3639 *
IndexSet locally_relevant_dofs_U;
3644 *
IndexSet locally_relevant_dofs_P;
3659 *
unsigned int LEVEL_SET;
3660 *
unsigned int RHO_TIMES_RHS;
3676 * std::vector<types::global_dof_index> boundary_values_id_u;
3677 * std::vector<types::global_dof_index> boundary_values_id_v;
3678 * std::vector<types::global_dof_index> boundary_values_id_w;
3679 * std::vector<double> boundary_values_u;
3680 * std::vector<double> boundary_values_v;
3681 * std::vector<double> boundary_values_w;
3686 *
bool rebuild_Matrix_U;
3687 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_Matrix_u;
3688 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_Matrix_v;
3689 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_Matrix_w;
3691 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_S;
3693 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_M;
3695 *
bool rebuild_Matrix_U_preconditioners;
3696 *
bool rebuild_S_M_preconditioners;
3724 * CONSTRUCTOR FOR LEVEL SET
3728 * NavierStokesSolver<dim>::NavierStokesSolver(
const unsigned int degree_LS,
3729 *
const unsigned int degree_U,
3730 *
const double time_step,
3732 *
const double rho_air,
3733 *
const double nu_air,
3734 *
const double rho_fluid,
3735 *
const double nu_fluid,
3737 *
const bool verbose,
3741 * mpi_communicator(mpi_communicator),
3743 * degree_LS(degree_LS),
3746 * degree_U(degree_U),
3751 * force_function(force_function),
3754 * This is dummy since rho and nu
functions won
't be used
3757 * rho_function(force_function),
3758 * nu_function(force_function),
3761 * rho_fluid(rho_fluid),
3762 * nu_fluid(nu_fluid),
3763 * time_step(time_step),
3768 * pcout(std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator)==0)),
3769 * rebuild_Matrix_U(true),
3770 * rebuild_S_M(true),
3771 * rebuild_Matrix_U_preconditioners(true),
3772 * rebuild_S_M_preconditioners(true)
3777 * CONSTRUCTOR NOT FOR LEVEL SET
3781 * NavierStokesSolver<dim>::NavierStokesSolver(const unsigned int degree_LS,
3782 * const unsigned int degree_U,
3783 * const double time_step,
3784 * Function<dim> &force_function,
3785 * Function<dim> &rho_function,
3786 * Function<dim> &nu_function,
3787 * const bool verbose,
3788 * parallel::distributed::Triangulation<dim> &triangulation,
3789 * MPI_Comm &mpi_communicator) :
3790 * mpi_communicator(mpi_communicator),
3791 * triangulation(triangulation),
3792 * degree_LS(degree_LS),
3793 * dof_handler_LS(triangulation),
3795 * degree_U(degree_U),
3796 * dof_handler_U(triangulation),
3798 * dof_handler_P(triangulation),
3800 * force_function(force_function),
3801 * rho_function(rho_function),
3802 * nu_function(nu_function),
3803 * time_step(time_step),
3807 * pcout(std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator)==0)),
3808 * rebuild_Matrix_U(true),
3809 * rebuild_S_M(true),
3810 * rebuild_Matrix_U_preconditioners(true),
3811 * rebuild_S_M_preconditioners(true)
3815 * NavierStokesSolver<dim>::~NavierStokesSolver()
3817 * dof_handler_LS.clear();
3818 * dof_handler_U.clear();
3819 * dof_handler_P.clear();
3824 * /////////////////////////////////////////////////////////
3825 * ////////////////// SETTERS AND GETTERS //////////////////
3826 * /////////////////////////////////////////////////////////
3830 * void NavierStokesSolver<dim>::set_rho_and_nu_functions(const Function<dim> &rho_function,
3831 * const Function<dim> &nu_function)
3833 * this->rho_function=rho_function;
3834 * this->nu_function=nu_function;
3838 * void NavierStokesSolver<dim>::initial_condition(PETScWrappers::MPI::Vector locally_relevant_solution_phi,
3839 * PETScWrappers::MPI::Vector locally_relevant_solution_u,
3840 * PETScWrappers::MPI::Vector locally_relevant_solution_v,
3841 * PETScWrappers::MPI::Vector locally_relevant_solution_p)
3843 * this->locally_relevant_solution_phi=locally_relevant_solution_phi;
3844 * this->locally_relevant_solution_u=locally_relevant_solution_u;
3845 * this->locally_relevant_solution_v=locally_relevant_solution_v;
3846 * this->locally_relevant_solution_p=locally_relevant_solution_p;
3849 * set old vectors to the initial condition (just for first time step)
3852 * save_old_solution();
3856 * void NavierStokesSolver<dim>::initial_condition(PETScWrappers::MPI::Vector locally_relevant_solution_phi,
3857 * PETScWrappers::MPI::Vector locally_relevant_solution_u,
3858 * PETScWrappers::MPI::Vector locally_relevant_solution_v,
3859 * PETScWrappers::MPI::Vector locally_relevant_solution_w,
3860 * PETScWrappers::MPI::Vector locally_relevant_solution_p)
3862 * this->locally_relevant_solution_phi=locally_relevant_solution_phi;
3863 * this->locally_relevant_solution_u=locally_relevant_solution_u;
3864 * this->locally_relevant_solution_v=locally_relevant_solution_v;
3865 * this->locally_relevant_solution_w=locally_relevant_solution_w;
3866 * this->locally_relevant_solution_p=locally_relevant_solution_p;
3869 * set old vectors to the initial condition (just for first time step)
3872 * save_old_solution();
3876 * void NavierStokesSolver<dim>::set_boundary_conditions(std::vector<types::global_dof_index> boundary_values_id_u,
3877 * std::vector<types::global_dof_index> boundary_values_id_v,
3878 * std::vector<double> boundary_values_u,
3879 * std::vector<double> boundary_values_v)
3881 * this->boundary_values_id_u=boundary_values_id_u;
3882 * this->boundary_values_id_v=boundary_values_id_v;
3883 * this->boundary_values_u=boundary_values_u;
3884 * this->boundary_values_v=boundary_values_v;
3888 * void NavierStokesSolver<dim>::set_boundary_conditions(std::vector<types::global_dof_index> boundary_values_id_u,
3889 * std::vector<types::global_dof_index> boundary_values_id_v,
3890 * std::vector<types::global_dof_index> boundary_values_id_w,
3891 * std::vector<double> boundary_values_u,
3892 * std::vector<double> boundary_values_v,
3893 * std::vector<double> boundary_values_w)
3895 * this->boundary_values_id_u=boundary_values_id_u;
3896 * this->boundary_values_id_v=boundary_values_id_v;
3897 * this->boundary_values_id_w=boundary_values_id_w;
3898 * this->boundary_values_u=boundary_values_u;
3899 * this->boundary_values_v=boundary_values_v;
3900 * this->boundary_values_w=boundary_values_w;
3904 * void NavierStokesSolver<dim>::set_velocity(PETScWrappers::MPI::Vector locally_relevant_solution_u,
3905 * PETScWrappers::MPI::Vector locally_relevant_solution_v)
3907 * this->locally_relevant_solution_u=locally_relevant_solution_u;
3908 * this->locally_relevant_solution_v=locally_relevant_solution_v;
3912 * void NavierStokesSolver<dim>::set_velocity(PETScWrappers::MPI::Vector locally_relevant_solution_u,
3913 * PETScWrappers::MPI::Vector locally_relevant_solution_v,
3914 * PETScWrappers::MPI::Vector locally_relevant_solution_w)
3916 * this->locally_relevant_solution_u=locally_relevant_solution_u;
3917 * this->locally_relevant_solution_v=locally_relevant_solution_v;
3918 * this->locally_relevant_solution_w=locally_relevant_solution_w;
3922 * void NavierStokesSolver<dim>::set_phi(PETScWrappers::MPI::Vector locally_relevant_solution_phi)
3924 * this->locally_relevant_solution_phi=locally_relevant_solution_phi;
3928 * void NavierStokesSolver<dim>::get_rho_and_nu(double phi)
3938 * else if (phi<-eps)
3942 * rho_value=rho_fluid*(1+H)/2.+rho_air*(1-H)/2.;
3943 * nu_value=nu_fluid*(1+H)/2.+nu_air*(1-H)/2.;
3946 * rho_value=rho_fluid*(1+phi)/2.+rho_air*(1-phi)/2.;
3947 * nu_value=nu_fluid*(1+phi)/2.+nu_air*(1-phi)/2.;
3953 * void NavierStokesSolver<dim>::get_pressure(PETScWrappers::MPI::Vector &locally_relevant_solution_p)
3955 * locally_relevant_solution_p=this->locally_relevant_solution_p;
3959 * void NavierStokesSolver<dim>::get_velocity(PETScWrappers::MPI::Vector &locally_relevant_solution_u,
3960 * PETScWrappers::MPI::Vector &locally_relevant_solution_v)
3962 * locally_relevant_solution_u=this->locally_relevant_solution_u;
3963 * locally_relevant_solution_v=this->locally_relevant_solution_v;
3967 * void NavierStokesSolver<dim>::get_velocity(PETScWrappers::MPI::Vector &locally_relevant_solution_u,
3968 * PETScWrappers::MPI::Vector &locally_relevant_solution_v,
3969 * PETScWrappers::MPI::Vector &locally_relevant_solution_w)
3971 * locally_relevant_solution_u=this->locally_relevant_solution_u;
3972 * locally_relevant_solution_v=this->locally_relevant_solution_v;
3973 * locally_relevant_solution_w=this->locally_relevant_solution_w;
3978 * ///////////////////////////////////////////////////
3979 * /////////// SETUP AND INITIAL CONDITION ///////////
3980 * ///////////////////////////////////////////////////
3984 * void NavierStokesSolver<dim>::setup()
3986 * pcout<<"***** SETUP IN NAVIER STOKES SOLVER *****"<<std::endl;
3988 * init_constraints();
3993 * void NavierStokesSolver<dim>::setup_DOF()
3996 * degree_MAX=std::max(degree_LS,degree_U);
4002 * dof_handler_LS.distribute_dofs(fe_LS);
4003 * locally_owned_dofs_LS = dof_handler_LS.locally_owned_dofs();
4004 * locally_relevant_dofs_LS = DoFTools::extract_locally_relevant_dofs(dof_handler_LS);
4010 * dof_handler_U.distribute_dofs(fe_U);
4011 * locally_owned_dofs_U = dof_handler_U.locally_owned_dofs();
4012 * locally_relevant_dofs_U = DoFTools::extract_locally_relevant_dofs(dof_handler_U);
4018 * dof_handler_P.distribute_dofs(fe_P);
4019 * locally_owned_dofs_P = dof_handler_P.locally_owned_dofs();
4020 * locally_relevant_dofs_P = DoFTools::extract_locally_relevant_dofs(dof_handler_P);
4024 * void NavierStokesSolver<dim>::setup_VECTORS()
4028 * init vectors for phi
4031 * locally_relevant_solution_phi.reinit(locally_owned_dofs_LS,locally_relevant_dofs_LS,
4032 * mpi_communicator);
4033 * locally_relevant_solution_phi=0;
4036 * init vectors for u
4039 * locally_relevant_solution_u.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4040 * mpi_communicator);
4041 * locally_relevant_solution_u=0;
4042 * completely_distributed_solution_u.reinit(locally_owned_dofs_U,mpi_communicator);
4043 * system_rhs_u.reinit(locally_owned_dofs_U,mpi_communicator);
4046 * init vectors for u_old
4049 * locally_relevant_solution_u_old.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4050 * mpi_communicator);
4051 * locally_relevant_solution_u_old=0;
4054 * init vectors for v
4057 * locally_relevant_solution_v.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4058 * mpi_communicator);
4059 * locally_relevant_solution_v=0;
4060 * completely_distributed_solution_v.reinit(locally_owned_dofs_U,mpi_communicator);
4061 * system_rhs_v.reinit(locally_owned_dofs_U,mpi_communicator);
4064 * init vectors for v_old
4067 * locally_relevant_solution_v_old.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4068 * mpi_communicator);
4069 * locally_relevant_solution_v_old=0;
4072 * init vectors for w
4075 * locally_relevant_solution_w.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4076 * mpi_communicator);
4077 * locally_relevant_solution_w=0;
4078 * completely_distributed_solution_w.reinit(locally_owned_dofs_U,mpi_communicator);
4079 * system_rhs_w.reinit(locally_owned_dofs_U,mpi_communicator);
4082 * init vectors for w_old
4085 * locally_relevant_solution_w_old.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4086 * mpi_communicator);
4087 * locally_relevant_solution_w_old=0;
4090 * init vectors for dpsi
4093 * locally_relevant_solution_psi.reinit(locally_owned_dofs_P,locally_relevant_dofs_P,
4094 * mpi_communicator);
4095 * locally_relevant_solution_psi=0;
4096 * system_rhs_psi.reinit(locally_owned_dofs_P,mpi_communicator);
4099 * init vectors for dpsi old
4102 * locally_relevant_solution_psi_old.reinit(locally_owned_dofs_P,locally_relevant_dofs_P,
4103 * mpi_communicator);
4104 * locally_relevant_solution_psi_old=0;
4107 * init vectors for q
4110 * completely_distributed_solution_q.reinit(locally_owned_dofs_P,mpi_communicator);
4111 * system_rhs_q.reinit(locally_owned_dofs_P,mpi_communicator);
4114 * init vectors for psi
4117 * completely_distributed_solution_psi.reinit(locally_owned_dofs_P,mpi_communicator);
4120 * init vectors for p
4123 * locally_relevant_solution_p.reinit(locally_owned_dofs_P,locally_relevant_dofs_P,
4124 * mpi_communicator);
4125 * locally_relevant_solution_p=0;
4126 * completely_distributed_solution_p.reinit(locally_owned_dofs_P,mpi_communicator);
4129 * ////////////////////////
4130 * Initialize constraints
4131 * ////////////////////////
4134 * init_constraints();
4137 * //////////////////
4139 * //////////////////
4140 * sparsity pattern for A
4143 * DynamicSparsityPattern dsp_Matrix(locally_relevant_dofs_U);
4144 * DoFTools::make_sparsity_pattern(dof_handler_U,dsp_Matrix,constraints,false);
4145 * SparsityTools::distribute_sparsity_pattern(dsp_Matrix,
4146 * dof_handler_U.locally_owned_dofs(),
4148 * locally_relevant_dofs_U);
4149 * system_Matrix_u.reinit(dof_handler_U.locally_owned_dofs(),
4150 * dof_handler_U.locally_owned_dofs(),
4152 * mpi_communicator);
4153 * system_Matrix_v.reinit(dof_handler_U.locally_owned_dofs(),
4154 * dof_handler_U.locally_owned_dofs(),
4156 * mpi_communicator);
4157 * system_Matrix_w.reinit(dof_handler_U.locally_owned_dofs(),
4158 * dof_handler_U.locally_owned_dofs(),
4160 * mpi_communicator);
4161 * rebuild_Matrix_U=true;
4164 * sparsity pattern for S
4167 * DynamicSparsityPattern dsp_S(locally_relevant_dofs_P);
4168 * DoFTools::make_sparsity_pattern(dof_handler_P,dsp_S,constraints_psi,false);
4169 * SparsityTools::distribute_sparsity_pattern(dsp_S,
4170 * dof_handler_P.locally_owned_dofs(),
4172 * locally_relevant_dofs_P);
4173 * system_S.reinit(dof_handler_P.locally_owned_dofs(),
4174 * dof_handler_P.locally_owned_dofs(),
4176 * mpi_communicator);
4179 * sparsity pattern for M
4182 * DynamicSparsityPattern dsp_M(locally_relevant_dofs_P);
4183 * DoFTools::make_sparsity_pattern(dof_handler_P,dsp_M,constraints_psi,false);
4184 * SparsityTools::distribute_sparsity_pattern(dsp_M,
4185 * dof_handler_P.locally_owned_dofs(),
4187 * locally_relevant_dofs_P);
4188 * system_M.reinit(dof_handler_P.locally_owned_dofs(),
4189 * dof_handler_P.locally_owned_dofs(),
4191 * mpi_communicator);
4196 * void NavierStokesSolver<dim>::init_constraints()
4203 * constraints.clear();
4204 * constraints.reinit(locally_relevant_dofs_U);
4205 * DoFTools::make_hanging_node_constraints(dof_handler_U,constraints);
4206 * constraints.close();
4209 * constraints for dpsi
4212 * constraints_psi.clear();
4213 * constraints_psi.reinit(locally_relevant_dofs_P);
4214 * DoFTools::make_hanging_node_constraints(dof_handler_P,constraints_psi);
4217 * if (constraints_psi.can_store_line(0))
4218 * constraints_psi.add_line(0); //constraint u0 = 0
4221 * constraints_psi.close();
4226 * ///////////////////////////////////////////////////
4227 * //////////////// ASSEMBLE SYSTEMS /////////////////
4228 * ///////////////////////////////////////////////////
4232 * void NavierStokesSolver<dim>::assemble_system_U()
4234 * if (rebuild_Matrix_U==true)
4236 * system_Matrix_u=0;
4237 * system_Matrix_v=0;
4238 * system_Matrix_w=0;
4244 * const QGauss<dim> quadrature_formula(degree_MAX+1);
4245 * FEValues<dim> fe_values_LS(fe_LS,quadrature_formula,
4246 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4247 * FEValues<dim> fe_values_U(fe_U,quadrature_formula,
4248 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4249 * FEValues<dim> fe_values_P(fe_P,quadrature_formula,
4250 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4252 * const unsigned int dofs_per_cell=fe_U.dofs_per_cell;
4253 * const unsigned int n_q_points=quadrature_formula.size();
4255 * FullMatrix<double> cell_A_u(dofs_per_cell,dofs_per_cell);
4256 * Vector<double> cell_rhs_u(dofs_per_cell);
4257 * Vector<double> cell_rhs_v(dofs_per_cell);
4258 * Vector<double> cell_rhs_w(dofs_per_cell);
4260 * std::vector<double> phiqnp1(n_q_points);
4262 * std::vector<double> uqn(n_q_points);
4263 * std::vector<double> uqnm1(n_q_points);
4264 * std::vector<double> vqn(n_q_points);
4265 * std::vector<double> vqnm1(n_q_points);
4266 * std::vector<double> wqn(n_q_points);
4267 * std::vector<double> wqnm1(n_q_points);
4271 * FOR Explicit nonlinearity
4272 * std::vector<Tensor<1, dim> > grad_un(n_q_points);
4273 * std::vector<Tensor<1, dim> > grad_vn(n_q_points);
4274 * std::vector<Tensor<1, dim> > grad_wn(n_q_points);
4275 * Tensor<1, dim> Un;
4281 * std::vector<Tensor<1, dim> > grad_pqn(n_q_points);
4282 * std::vector<Tensor<1, dim> > grad_psiqn(n_q_points);
4283 * std::vector<Tensor<1, dim> > grad_psiqnm1(n_q_points);
4285 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
4286 * std::vector<Tensor<1, dim> > shape_grad(dofs_per_cell);
4287 * std::vector<double> shape_value(dofs_per_cell);
4292 * double pressure_grad_u;
4293 * double pressure_grad_v;
4294 * double pressure_grad_w;
4300 * Vector<double> force_terms(dim);
4302 * typename DoFHandler<dim>::active_cell_iterator
4303 * cell_U=dof_handler_U.begin_active(), endc_U=dof_handler_U.end();
4304 * typename DoFHandler<dim>::active_cell_iterator cell_P=dof_handler_P.begin_active();
4305 * typename DoFHandler<dim>::active_cell_iterator cell_LS=dof_handler_LS.begin_active();
4307 * for (; cell_U!=endc_U; ++cell_U,++cell_P,++cell_LS)
4308 * if (cell_U->is_locally_owned())
4315 * fe_values_LS.reinit(cell_LS);
4316 * fe_values_U.reinit(cell_U);
4317 * fe_values_P.reinit(cell_P);
4321 * get function values for LS
4324 * fe_values_LS.get_function_values(locally_relevant_solution_phi,phiqnp1);
4327 * get function values for U
4330 * fe_values_U.get_function_values(locally_relevant_solution_u,uqn);
4331 * fe_values_U.get_function_values(locally_relevant_solution_u_old,uqnm1);
4332 * fe_values_U.get_function_values(locally_relevant_solution_v,vqn);
4333 * fe_values_U.get_function_values(locally_relevant_solution_v_old,vqnm1);
4336 * fe_values_U.get_function_values(locally_relevant_solution_w,wqn);
4337 * fe_values_U.get_function_values(locally_relevant_solution_w_old,wqnm1);
4341 * For explicit nonlinearity
4342 * get gradient values for U
4343 * fe_values_U.get_function_gradients(locally_relevant_solution_u,grad_un);
4344 * fe_values_U.get_function_gradients(locally_relevant_solution_v,grad_vn);
4346 * fe_values_U.get_function_gradients(locally_relevant_solution_w,grad_wn);
4350 * get values and gradients for p and dpsi
4353 * fe_values_P.get_function_gradients(locally_relevant_solution_p,grad_pqn);
4354 * fe_values_P.get_function_gradients(locally_relevant_solution_psi,grad_psiqn);
4355 * fe_values_P.get_function_gradients(locally_relevant_solution_psi_old,grad_psiqnm1);
4357 * for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
4359 * const double JxW=fe_values_U.JxW(q_point);
4360 * for (unsigned int i=0; i<dofs_per_cell; ++i)
4362 * shape_grad[i]=fe_values_U.shape_grad(i,q_point);
4363 * shape_value[i]=fe_values_U.shape_value(i,q_point);
4366 * pressure_grad_u=(grad_pqn[q_point][0]+4./3*grad_psiqn[q_point][0]-1./3*grad_psiqnm1[q_point][0]);
4367 * pressure_grad_v=(grad_pqn[q_point][1]+4./3*grad_psiqn[q_point][1]-1./3*grad_psiqnm1[q_point][1]);
4369 * pressure_grad_w=(grad_pqn[q_point][2]+4./3*grad_psiqn[q_point][2]-1./3*grad_psiqnm1[q_point][2]);
4371 * if (LEVEL_SET==1) // use level set to define rho and nu
4372 * get_rho_and_nu(phiqnp1[q_point]);
4373 * else // rho and nu are defined through functions
4375 * rho_value=rho_function.value(fe_values_U.quadrature_point(q_point));
4376 * nu_value=nu_function.value(fe_values_U.quadrature_point(q_point));
4381 * Non-linearity: for semi-implicit
4384 * u_star=2*uqn[q_point]-uqnm1[q_point];
4385 * v_star=2*vqn[q_point]-vqnm1[q_point];
4387 * w_star=2*wqn[q_point]-wqnm1[q_point];
4391 * for explicit nonlinearity
4392 * Un[0] = uqn[q_point];
4393 * Un[1] = vqn[q_point];
4395 * Un[2] = wqn[q_point];
4399 * double nonlinearity_u = Un*grad_un[q_point];
4400 * double nonlinearity_v = Un*grad_vn[q_point];
4401 * double nonlinearity_w = 0;
4403 * nonlinearity_w = Un*grad_wn[q_point];
4409 * rho_star=rho_value; // This is because we consider rho*u_t instead of (rho*u)_t
4417 * force_function.vector_value(fe_values_U.quadrature_point(q_point),force_terms);
4418 * force_u=force_terms[0];
4419 * force_v=force_terms[1];
4421 * force_w=force_terms[2];
4422 * if (RHO_TIMES_RHS==1)
4430 * for (unsigned int i=0; i<dofs_per_cell; ++i)
4432 * cell_rhs_u(i)+=((4./3*rho*uqn[q_point]-1./3*rho*uqnm1[q_point]
4433 * +2./3*time_step*(force_u-pressure_grad_u)
4436 * -2./3*time_step*rho*nonlinearity_u
4439 * )*shape_value[i])*JxW;
4440 * cell_rhs_v(i)+=((4./3*rho*vqn[q_point]-1./3*rho*vqnm1[q_point]
4441 * +2./3*time_step*(force_v-pressure_grad_v)
4444 * -2./3*time_step*rho*nonlinearity_v
4447 * )*shape_value[i])*JxW;
4449 * cell_rhs_w(i)+=((4./3*rho*wqn[q_point]-1./3*rho*wqnm1[q_point]
4450 * +2./3*time_step*(force_w-pressure_grad_w)
4453 * -2./3*time_step*rho*nonlinearity_w
4456 * )*shape_value[i])*JxW;
4457 * if (rebuild_Matrix_U==true)
4458 * for (unsigned int j=0; j<dofs_per_cell; ++j)
4461 * cell_A_u(i,j)+=(rho_star*shape_value[i]*shape_value[j]
4462 * +2./3*time_step*nu_value*(shape_grad[i]*shape_grad[j])
4463 * +2./3*time_step*rho*shape_value[i]
4464 * *(u_star*shape_grad[j][0]+v_star*shape_grad[j][1]) // semi-implicit NL
4467 * cell_A_u(i,j)+=(rho_star*shape_value[i]*shape_value[j]
4468 * +2./3*time_step*nu_value*(shape_grad[i]*shape_grad[j])
4469 * +2./3*time_step*rho*shape_value[i]
4470 * *(u_star*shape_grad[j][0]+v_star*shape_grad[j][1]+w_star*shape_grad[j][2]) // semi-implicit NL
4475 * cell_U->get_dof_indices(local_dof_indices);
4481 * if (rebuild_Matrix_U==true)
4482 * constraints.distribute_local_to_global(cell_A_u,local_dof_indices,system_Matrix_u);
4483 * constraints.distribute_local_to_global(cell_rhs_u,local_dof_indices,system_rhs_u);
4484 * constraints.distribute_local_to_global(cell_rhs_v,local_dof_indices,system_rhs_v);
4486 * constraints.distribute_local_to_global(cell_rhs_w,local_dof_indices,system_rhs_w);
4488 * system_rhs_u.compress(VectorOperation::add);
4489 * system_rhs_v.compress(VectorOperation::add);
4490 * if (dim==3) system_rhs_w.compress(VectorOperation::add);
4491 * if (rebuild_Matrix_U==true)
4493 * system_Matrix_u.compress(VectorOperation::add);
4494 * system_Matrix_v.copy_from(system_Matrix_u);
4496 * system_Matrix_w.copy_from(system_Matrix_u);
4500 * BOUNDARY CONDITIONS
4503 * system_rhs_u.set(boundary_values_id_u,boundary_values_u);
4504 * system_rhs_u.compress(VectorOperation::insert);
4505 * system_rhs_v.set(boundary_values_id_v,boundary_values_v);
4506 * system_rhs_v.compress(VectorOperation::insert);
4509 * system_rhs_w.set(boundary_values_id_w,boundary_values_w);
4510 * system_rhs_w.compress(VectorOperation::insert);
4512 * if (rebuild_Matrix_U)
4514 * system_Matrix_u.clear_rows(boundary_values_id_u,1);
4515 * system_Matrix_v.clear_rows(boundary_values_id_v,1);
4517 * system_Matrix_w.clear_rows(boundary_values_id_w,1);
4518 * if (rebuild_Matrix_U_preconditioners)
4525 * rebuild_Matrix_U_preconditioners=false;
4526 * preconditioner_Matrix_u.reset(new PETScWrappers::PreconditionBoomerAMG
4527 * (system_Matrix_u,PETScWrappers::PreconditionBoomerAMG::AdditionalData(false)));
4528 * preconditioner_Matrix_v.reset( new PETScWrappers::PreconditionBoomerAMG
4529 * (system_Matrix_v,PETScWrappers::PreconditionBoomerAMG::AdditionalData(false)));
4531 * preconditioner_Matrix_w.reset(new PETScWrappers::PreconditionBoomerAMG
4532 * (system_Matrix_w,PETScWrappers::PreconditionBoomerAMG::AdditionalData(false)));
4535 * rebuild_Matrix_U=true;
4539 * void NavierStokesSolver<dim>::assemble_system_dpsi_q()
4541 * if (rebuild_S_M==true)
4549 * const QGauss<dim> quadrature_formula(degree_MAX+1);
4551 * FEValues<dim> fe_values_U(fe_U,quadrature_formula,
4552 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4553 * FEValues<dim> fe_values_P(fe_P,quadrature_formula,
4554 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4555 * FEValues<dim> fe_values_LS(fe_LS,quadrature_formula,
4556 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4558 * const unsigned int dofs_per_cell=fe_P.dofs_per_cell;
4559 * const unsigned int n_q_points=quadrature_formula.size();
4561 * FullMatrix<double> cell_S(dofs_per_cell,dofs_per_cell);
4562 * FullMatrix<double> cell_M(dofs_per_cell,dofs_per_cell);
4563 * Vector<double> cell_rhs_psi(dofs_per_cell);
4564 * Vector<double> cell_rhs_q(dofs_per_cell);
4566 * std::vector<double> phiqnp1(n_q_points);
4567 * std::vector<Tensor<1, dim> > gunp1(n_q_points);
4568 * std::vector<Tensor<1, dim> > gvnp1(n_q_points);
4569 * std::vector<Tensor<1, dim> > gwnp1(n_q_points);
4571 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
4572 * std::vector<double> shape_value(dofs_per_cell);
4573 * std::vector<Tensor<1, dim> > shape_grad(dofs_per_cell);
4575 * typename DoFHandler<dim>::active_cell_iterator
4576 * cell_P=dof_handler_P.begin_active(), endc_P=dof_handler_P.end();
4577 * typename DoFHandler<dim>::active_cell_iterator cell_U=dof_handler_U.begin_active();
4578 * typename DoFHandler<dim>::active_cell_iterator cell_LS=dof_handler_LS.begin_active();
4580 * for (; cell_P!=endc_P; ++cell_P,++cell_U,++cell_LS)
4581 * if (cell_P->is_locally_owned())
4588 * fe_values_P.reinit(cell_P);
4589 * fe_values_U.reinit(cell_U);
4590 * fe_values_LS.reinit(cell_LS);
4594 * get function values for LS
4597 * fe_values_LS.get_function_values(locally_relevant_solution_phi,phiqnp1);
4601 * get function grads for u and v
4604 * fe_values_U.get_function_gradients(locally_relevant_solution_u,gunp1);
4605 * fe_values_U.get_function_gradients(locally_relevant_solution_v,gvnp1);
4607 * fe_values_U.get_function_gradients(locally_relevant_solution_w,gwnp1);
4609 * for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
4611 * const double JxW=fe_values_P.JxW(q_point);
4612 * double divU = gunp1[q_point][0]+gvnp1[q_point][1];
4613 * if (dim==3) divU += gwnp1[q_point][2];
4614 * for (unsigned int i=0; i<dofs_per_cell; ++i)
4616 * shape_value[i]=fe_values_P.shape_value(i,q_point);
4617 * shape_grad[i]=fe_values_P.shape_grad(i,q_point);
4619 * if (LEVEL_SET==1) // use level set to define rho and nu
4620 * get_rho_and_nu (phiqnp1[q_point]);
4621 * else // rho and nu are defined through functions
4622 * nu_value=nu_function.value(fe_values_U.quadrature_point(q_point));
4624 * for (unsigned int i=0; i<dofs_per_cell; ++i)
4626 * cell_rhs_psi(i)+=-3./2./time_step*rho_min*divU*shape_value[i]*JxW;
4627 * cell_rhs_q(i)-=nu_value*divU*shape_value[i]*JxW;
4628 * if (rebuild_S_M==true)
4630 * for (unsigned int j=0; j<dofs_per_cell; ++j)
4633 * cell_S(i,j)+=shape_grad[i]*shape_grad[j]*JxW+1E-10;
4634 * cell_M(i,j)+=shape_value[i]*shape_value[j]*JxW;
4638 * cell_S(i,j)+=shape_grad[i]*shape_grad[j]*JxW;
4639 * cell_M(i,j)+=shape_value[i]*shape_value[j]*JxW;
4644 * cell_P->get_dof_indices(local_dof_indices);
4650 * if (rebuild_S_M==true)
4652 * constraints_psi.distribute_local_to_global(cell_S,local_dof_indices,system_S);
4653 * constraints_psi.distribute_local_to_global(cell_M,local_dof_indices,system_M);
4655 * constraints_psi.distribute_local_to_global(cell_rhs_q,local_dof_indices,system_rhs_q);
4656 * constraints_psi.distribute_local_to_global(cell_rhs_psi,local_dof_indices,system_rhs_psi);
4658 * if (rebuild_S_M==true)
4660 * system_M.compress(VectorOperation::add);
4661 * system_S.compress(VectorOperation::add);
4662 * if (rebuild_S_M_preconditioners)
4664 * rebuild_S_M_preconditioners=false;
4665 * preconditioner_S.reset(new PETScWrappers::PreconditionBoomerAMG
4666 * (system_S,PETScWrappers::PreconditionBoomerAMG::AdditionalData(true)));
4667 * preconditioner_M.reset(new PETScWrappers::PreconditionBoomerAMG
4668 * (system_M,PETScWrappers::PreconditionBoomerAMG::AdditionalData(true)));
4671 * system_rhs_psi.compress(VectorOperation::add);
4672 * system_rhs_q.compress(VectorOperation::add);
4673 * rebuild_S_M=false;
4678 * ///////////////////////////////////////////////////
4679 * ///////////////////// SOLVERS /////////////////////
4680 * ///////////////////////////////////////////////////
4684 * void NavierStokesSolver<dim>::solve_U(const AffineConstraints<double> &constraints,
4685 * PETScWrappers::MPI::SparseMatrix &Matrix,
4686 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
4687 * PETScWrappers::MPI::Vector &completely_distributed_solution,
4688 * const PETScWrappers::MPI::Vector &rhs)
4690 * SolverControl solver_control(dof_handler_U.n_dofs(),1e-6);
4693 * PETScWrappers::SolverCG solver(solver_control, mpi_communicator);
4694 * PETScWrappers::SolverGMRES solver(solver_control, mpi_communicator);
4695 * PETScWrappers::SolverChebychev solver(solver_control, mpi_communicator);
4698 * PETScWrappers::SolverBicgstab solver(solver_control,mpi_communicator);
4699 * constraints.distribute(completely_distributed_solution);
4700 * solver.solve(Matrix,completely_distributed_solution,rhs,*preconditioner);
4701 * constraints.distribute(completely_distributed_solution);
4702 * if (solver_control.last_step() > MAX_NUM_ITER_TO_RECOMPUTE_PRECONDITIONER)
4703 * rebuild_Matrix_U_preconditioners=true;
4704 * if (verbose==true)
4705 * pcout<<" Solved U in "<<solver_control.last_step()<<" iterations."<<std::endl;
4709 * void NavierStokesSolver<dim>::solve_P(const AffineConstraints<double> &constraints,
4710 * PETScWrappers::MPI::SparseMatrix &Matrix,
4711 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
4712 * PETScWrappers::MPI::Vector &completely_distributed_solution,
4713 * const PETScWrappers::MPI::Vector &rhs)
4715 * SolverControl solver_control(dof_handler_P.n_dofs(),1e-6);
4716 * PETScWrappers::SolverCG solver(solver_control,mpi_communicator);
4719 * PETScWrappers::SolverGMRES solver(solver_control, mpi_communicator);
4722 * constraints.distribute(completely_distributed_solution);
4723 * solver.solve(Matrix,completely_distributed_solution,rhs,*preconditioner);
4724 * constraints.distribute(completely_distributed_solution);
4725 * if (solver_control.last_step() > MAX_NUM_ITER_TO_RECOMPUTE_PRECONDITIONER)
4726 * rebuild_S_M_preconditioners=true;
4727 * if (verbose==true)
4728 * pcout<<" Solved P in "<<solver_control.last_step()<<" iterations."<<std::endl;
4733 * ///////////////////////////////////////////////////
4734 * ////////////// get different fields ///////////////
4735 * ///////////////////////////////////////////////////
4739 * void NavierStokesSolver<dim>::get_velocity()
4741 * assemble_system_U();
4742 * save_old_solution();
4743 * solve_U(constraints,system_Matrix_u,preconditioner_Matrix_u,completely_distributed_solution_u,system_rhs_u);
4744 * locally_relevant_solution_u=completely_distributed_solution_u;
4745 * solve_U(constraints,system_Matrix_v,preconditioner_Matrix_v,completely_distributed_solution_v,system_rhs_v);
4746 * locally_relevant_solution_v=completely_distributed_solution_v;
4749 * solve_U(constraints,system_Matrix_w,preconditioner_Matrix_w,completely_distributed_solution_w,system_rhs_w);
4750 * locally_relevant_solution_w=completely_distributed_solution_w;
4755 * void NavierStokesSolver<dim>::get_pressure()
4762 * assemble_system_dpsi_q();
4763 * solve_P(constraints_psi,system_S,preconditioner_S,completely_distributed_solution_psi,system_rhs_psi);
4764 * locally_relevant_solution_psi=completely_distributed_solution_psi;
4770 * solve_P(constraints,system_M,preconditioner_M,completely_distributed_solution_q,system_rhs_q);
4773 * UPDATE THE PRESSURE
4776 * completely_distributed_solution_p.add(1,completely_distributed_solution_psi);
4777 * completely_distributed_solution_p.add(1,completely_distributed_solution_q);
4778 * locally_relevant_solution_p = completely_distributed_solution_p;
4783 * ///////////////////////////////////////////////////
4784 * ///////////////////// DO STEPS ////////////////////
4785 * ///////////////////////////////////////////////////
4789 * void NavierStokesSolver<dim>::nth_time_step()
4797 * ///////////////////////////////////////////////////
4798 * ////////////////////// OTHERS /////////////////////
4799 * ///////////////////////////////////////////////////
4803 * void NavierStokesSolver<dim>::save_old_solution()
4805 * locally_relevant_solution_u_old=locally_relevant_solution_u;
4806 * locally_relevant_solution_v_old=locally_relevant_solution_v;
4807 * locally_relevant_solution_w_old=locally_relevant_solution_w;
4808 * locally_relevant_solution_psi_old=locally_relevant_solution_psi;
4814<a name="ann-TestLevelSet.cc"></a>
4815<h1>Annotated version of TestLevelSet.cc</h1>
4821 * /* -----------------------------------------------------------------------------
4823 * * SPDX-License-Identifier: LGPL-2.1-or-later
4824 * * Copyright (C) 2016 Manuel Quezada de Luna
4826 * * This file is part of the deal.II code gallery.
4828 * * -----------------------------------------------------------------------------
4831 * #include <deal.II/base/quadrature_lib.h>
4832 * #include <deal.II/base/function.h>
4833 * #include <deal.II/lac/affine_constraints.h>
4834 * #include <deal.II/lac/vector.h>
4835 * #include <deal.II/lac/full_matrix.h>
4836 * #include <deal.II/lac/solver_cg.h>
4837 * #include <deal.II/lac/petsc_sparse_matrix.h>
4838 * #include <deal.II/lac/petsc_vector.h>
4839 * #include <deal.II/lac/petsc_solver.h>
4840 * #include <deal.II/lac/petsc_precondition.h>
4841 * #include <deal.II/grid/grid_generator.h>
4842 * #include <deal.II/grid/tria_accessor.h>
4843 * #include <deal.II/grid/tria_iterator.h>
4844 * #include <deal.II/dofs/dof_handler.h>
4845 * #include <deal.II/dofs/dof_accessor.h>
4846 * #include <deal.II/dofs/dof_tools.h>
4847 * #include <deal.II/fe/fe_values.h>
4848 * #include <deal.II/fe/fe_q.h>
4849 * #include <deal.II/numerics/vector_tools.h>
4850 * #include <deal.II/numerics/data_out.h>
4851 * #include <deal.II/numerics/error_estimator.h>
4852 * #include <deal.II/base/utilities.h>
4853 * #include <deal.II/base/conditional_ostream.h>
4854 * #include <deal.II/base/index_set.h>
4855 * #include <deal.II/lac/sparsity_tools.h>
4856 * #include <deal.II/distributed/tria.h>
4857 * #include <deal.II/distributed/grid_refinement.h>
4858 * #include <deal.II/lac/petsc_vector.h>
4859 * #include <deal.II/base/convergence_table.h>
4860 * #include <deal.II/base/timer.h>
4861 * #include <deal.II/base/parameter_handler.h>
4862 * #include <deal.II/grid/grid_tools.h>
4863 * #include <deal.II/fe/mapping_q.h>
4864 * #include <deal.II/fe/fe_system.h>
4866 * #include <fstream>
4867 * #include <iostream>
4870 * using namespace dealii;
4874 * ///////////////////////
4875 * FOR TRANSPORT PROBLEM
4876 * ///////////////////////
4880 * #define FORWARD_EULER 0
4887 * #define CIRCULAR_ROTATION 0
4888 * #define DIAGONAL_ADVECTION 1
4894 * #define VARIABLE_VELOCITY 0
4896 * #include "utilities_test_LS.cc"
4897 * #include "LevelSetSolver.cc"
4901 * ///////////////////////////////////////////////////
4902 * /////////////////// MAIN CLASS ////////////////////
4903 * ///////////////////////////////////////////////////
4906 * template <int dim>
4907 * class TestLevelSet
4910 * TestLevelSet (const unsigned int degree_LS,
4911 * const unsigned int degree_U);
4921 * void set_boundary_inlet();
4922 * void get_boundary_values_phi(std::vector<unsigned int> &boundary_values_id_phi,
4923 * std::vector<double> &boundary_values_phi);
4929 * void get_interpolated_velocity();
4932 * SETUP AND INIT CONDITIONS
4936 * void initial_condition();
4937 * void init_constraints();
4943 * void process_solution(parallel::distributed::Triangulation<dim> &triangulation,
4944 * DoFHandler<dim> &dof_handler_LS,
4945 * PETScWrappers::MPI::Vector &solution);
4946 * void output_results();
4947 * void output_solution();
4954 * PETScWrappers::MPI::Vector locally_relevant_solution_phi;
4955 * PETScWrappers::MPI::Vector locally_relevant_solution_u;
4956 * PETScWrappers::MPI::Vector locally_relevant_solution_v;
4957 * PETScWrappers::MPI::Vector locally_relevant_solution_w;
4958 * PETScWrappers::MPI::Vector completely_distributed_solution_phi;
4959 * PETScWrappers::MPI::Vector completely_distributed_solution_u;
4960 * PETScWrappers::MPI::Vector completely_distributed_solution_v;
4961 * PETScWrappers::MPI::Vector completely_distributed_solution_w;
4967 * std::vector<unsigned int> boundary_values_id_phi;
4968 * std::vector<double> boundary_values_phi;
4975 * MPI_Comm mpi_communicator;
4976 * parallel::distributed::Triangulation<dim> triangulation;
4980 * DoFHandler<dim> dof_handler_LS;
4982 * IndexSet locally_owned_dofs_LS;
4983 * IndexSet locally_relevant_dofs_LS;
4986 * DoFHandler<dim> dof_handler_U;
4988 * IndexSet locally_owned_dofs_U;
4989 * IndexSet locally_relevant_dofs_U;
4991 * DoFHandler<dim> dof_handler_U_disp_field;
4992 * FESystem<dim> fe_U_disp_field;
4993 * IndexSet locally_owned_dofs_U_disp_field;
4994 * IndexSet locally_relevant_dofs_U_disp_field;
4996 * AffineConstraints<double> constraints;
4997 * AffineConstraints<double> constraints_disp_field;
5001 * double final_time;
5002 * unsigned int timestep_number;
5007 * int sharpness_integer;
5009 * unsigned int n_refinement;
5010 * unsigned int output_number;
5011 * double output_time;
5015 * ConditionalOStream pcout;
5022 * double cK; //compression coeff
5023 * double cE; //entropy-visc coeff
5024 * unsigned int TRANSPORT_TIME_INTEGRATION;
5025 * std::string ALGORITHM;
5026 * unsigned int PROBLEM;
5030 * FOR RECONSTRUCTION OF MATERIAL FIELDS
5033 * double eps, rho_air, rho_fluid;
5040 * PETScWrappers::MPI::SparseMatrix matrix_MC, matrix_MC_tnm1;
5041 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_MC;
5045 * template <int dim>
5046 * TestLevelSet<dim>::TestLevelSet (const unsigned int degree_LS,
5047 * const unsigned int degree_U)
5049 * mpi_communicator (MPI_COMM_WORLD),
5050 * triangulation (mpi_communicator,
5051 * typename Triangulation<dim>::MeshSmoothing
5052 * (Triangulation<dim>::smoothing_on_refinement |
5053 * Triangulation<dim>::smoothing_on_coarsening)),
5054 * degree_LS(degree_LS),
5055 * dof_handler_LS (triangulation),
5056 * fe_LS (degree_LS),
5057 * degree_U(degree_U),
5058 * dof_handler_U (triangulation),
5060 * dof_handler_U_disp_field(triangulation),
5061 * fe_U_disp_field(FE_Q<dim>(degree_U),dim),
5062 * pcout (std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator)== 0))
5065 * template <int dim>
5066 * TestLevelSet<dim>::~TestLevelSet ()
5068 * dof_handler_U_disp_field.clear();
5069 * dof_handler_LS.clear ();
5070 * dof_handler_U.clear ();
5079 * template <int dim>
5080 * void TestLevelSet<dim>::get_interpolated_velocity()
5087 * completely_distributed_solution_u = 0;
5088 * VectorTools::interpolate(dof_handler_U,
5089 * ExactU<dim>(PROBLEM,time),
5090 * completely_distributed_solution_u);
5091 * constraints.distribute (completely_distributed_solution_u);
5092 * locally_relevant_solution_u = completely_distributed_solution_u;
5098 * completely_distributed_solution_v = 0;
5099 * VectorTools::interpolate(dof_handler_U,
5100 * ExactV<dim>(PROBLEM,time),
5101 * completely_distributed_solution_v);
5102 * constraints.distribute (completely_distributed_solution_v);
5103 * locally_relevant_solution_v = completely_distributed_solution_v;
5106 * completely_distributed_solution_w = 0;
5107 * VectorTools::interpolate(dof_handler_U,
5108 * ExactW<dim>(PROBLEM,time),
5109 * completely_distributed_solution_w);
5110 * constraints.distribute (completely_distributed_solution_w);
5111 * locally_relevant_solution_w = completely_distributed_solution_w;
5122 * template <int dim>
5123 * void TestLevelSet<dim>::set_boundary_inlet()
5125 * const QGauss<dim-1> face_quadrature_formula(1); // center of the face
5126 * FEFaceValues<dim> fe_face_values (fe_U,face_quadrature_formula,
5127 * update_values | update_quadrature_points |
5128 * update_normal_vectors);
5129 * const unsigned int n_face_q_points = face_quadrature_formula.size();
5130 * std::vector<double> u_value (n_face_q_points);
5131 * std::vector<double> v_value (n_face_q_points);
5132 * std::vector<double> w_value (n_face_q_points);
5134 * typename DoFHandler<dim>::active_cell_iterator
5135 * cell_U = dof_handler_U.begin_active(),
5136 * endc_U = dof_handler_U.end();
5139 * for (; cell_U!=endc_U; ++cell_U)
5140 * if (cell_U->is_locally_owned())
5141 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
5142 * if (cell_U->face(face)->at_boundary())
5144 * fe_face_values.reinit(cell_U,face);
5145 * fe_face_values.get_function_values(locally_relevant_solution_u,u_value);
5146 * fe_face_values.get_function_values(locally_relevant_solution_v,v_value);
5148 * fe_face_values.get_function_values(locally_relevant_solution_w,w_value);
5153 * if (fe_face_values.normal_vector(0)*u < -1e-14)
5154 * cell_U->face(face)->set_boundary_id(10);
5158 * template <int dim>
5159 * void TestLevelSet<dim>::get_boundary_values_phi(std::vector<unsigned int> &boundary_values_id_phi,
5160 * std::vector<double> &boundary_values_phi)
5162 * std::map<unsigned int, double> map_boundary_values_phi;
5163 * unsigned int boundary_id=0;
5165 * set_boundary_inlet();
5166 * boundary_id=10; // inlet
5167 * VectorTools::interpolate_boundary_values (dof_handler_LS,
5168 * boundary_id,BoundaryPhi<dim>(),
5169 * map_boundary_values_phi);
5171 * boundary_values_id_phi.resize(map_boundary_values_phi.size());
5172 * boundary_values_phi.resize(map_boundary_values_phi.size());
5173 * std::map<unsigned int,double>::const_iterator boundary_value_phi = map_boundary_values_phi.begin();
5174 * for (int i=0; boundary_value_phi !=map_boundary_values_phi.end(); ++boundary_value_phi, ++i)
5176 * boundary_values_id_phi[i]=boundary_value_phi->first;
5177 * boundary_values_phi[i]=boundary_value_phi->second;
5183 * ///////////////////////////////
5184 * SETUP AND INITIAL CONDITIONS
5185 * //////////////////////////////
5188 * template <int dim>
5189 * void TestLevelSet<dim>::setup()
5191 * degree = std::max(degree_LS,degree_U);
5197 * dof_handler_LS.distribute_dofs (fe_LS);
5198 * locally_owned_dofs_LS = dof_handler_LS.locally_owned_dofs ();
5199 * locally_relevant_dofs_LS = DoFTools::extract_locally_relevant_dofs (dof_handler_LS);
5205 * dof_handler_U.distribute_dofs (fe_U);
5206 * locally_owned_dofs_U = dof_handler_U.locally_owned_dofs ();
5207 * locally_relevant_dofs_U = DoFTools::extract_locally_relevant_dofs (dof_handler_U);
5210 * setup system U for disp field
5213 * dof_handler_U_disp_field.distribute_dofs (fe_U_disp_field);
5214 * locally_owned_dofs_U_disp_field = dof_handler_U_disp_field.locally_owned_dofs ();
5215 * locally_relevant_dofs_U_disp_field = DoFTools::extract_locally_relevant_dofs (dof_handler_U_disp_field);
5218 * init vectors for phi
5221 * locally_relevant_solution_phi.reinit(locally_owned_dofs_LS,
5222 * locally_relevant_dofs_LS,
5223 * mpi_communicator);
5224 * locally_relevant_solution_phi = 0;
5225 * completely_distributed_solution_phi.reinit(mpi_communicator,
5226 * dof_handler_LS.n_dofs(),
5227 * dof_handler_LS.n_locally_owned_dofs());
5230 * init vectors for u
5233 * locally_relevant_solution_u.reinit(locally_owned_dofs_U,
5234 * locally_relevant_dofs_U,
5235 * mpi_communicator);
5236 * locally_relevant_solution_u = 0;
5237 * completely_distributed_solution_u.reinit(mpi_communicator,
5238 * dof_handler_U.n_dofs(),
5239 * dof_handler_U.n_locally_owned_dofs());
5242 * init vectors for v
5245 * locally_relevant_solution_v.reinit(locally_owned_dofs_U,
5246 * locally_relevant_dofs_U,
5247 * mpi_communicator);
5248 * locally_relevant_solution_v = 0;
5249 * completely_distributed_solution_v.reinit(mpi_communicator,
5250 * dof_handler_U.n_dofs(),
5251 * dof_handler_U.n_locally_owned_dofs());
5254 * init vectors for w
5257 * locally_relevant_solution_w.reinit(locally_owned_dofs_U,
5258 * locally_relevant_dofs_U,
5259 * mpi_communicator);
5260 * locally_relevant_solution_w = 0;
5261 * completely_distributed_solution_w.reinit(mpi_communicator,
5262 * dof_handler_U.n_dofs(),
5263 * dof_handler_U.n_locally_owned_dofs());
5264 * init_constraints();
5270 * DynamicSparsityPattern dsp (locally_relevant_dofs_LS);
5271 * DoFTools::make_sparsity_pattern (dof_handler_LS,dsp,constraints,false);
5272 * SparsityTools::distribute_sparsity_pattern (dsp,
5273 * dof_handler_LS.n_locally_owned_dofs_per_processor(),
5275 * locally_relevant_dofs_LS);
5276 * matrix_MC.reinit (mpi_communicator,
5278 * dof_handler_LS.n_locally_owned_dofs_per_processor(),
5279 * dof_handler_LS.n_locally_owned_dofs_per_processor(),
5280 * Utilities::MPI::this_mpi_process(mpi_communicator));
5281 * matrix_MC_tnm1.reinit (mpi_communicator,
5283 * dof_handler_LS.n_locally_owned_dofs_per_processor(),
5284 * dof_handler_LS.n_locally_owned_dofs_per_processor(),
5285 * Utilities::MPI::this_mpi_process(mpi_communicator));
5288 * template <int dim>
5289 * void TestLevelSet<dim>::initial_condition()
5294 * Initial conditions
5295 * init condition for phi
5298 * completely_distributed_solution_phi = 0;
5299 * VectorTools::interpolate(dof_handler_LS,
5300 * InitialPhi<dim>(PROBLEM, sharpness),
5303 * Functions::ZeroFunction<dim>(),
5306 * completely_distributed_solution_phi);
5307 * constraints.distribute (completely_distributed_solution_phi);
5308 * locally_relevant_solution_phi = completely_distributed_solution_phi;
5311 * init condition for u=0
5314 * completely_distributed_solution_u = 0;
5315 * VectorTools::interpolate(dof_handler_U,
5316 * ExactU<dim>(PROBLEM,time),
5317 * completely_distributed_solution_u);
5318 * constraints.distribute (completely_distributed_solution_u);
5319 * locally_relevant_solution_u = completely_distributed_solution_u;
5322 * init condition for v
5325 * completely_distributed_solution_v = 0;
5326 * VectorTools::interpolate(dof_handler_U,
5327 * ExactV<dim>(PROBLEM,time),
5328 * completely_distributed_solution_v);
5329 * constraints.distribute (completely_distributed_solution_v);
5330 * locally_relevant_solution_v = completely_distributed_solution_v;
5333 * template <int dim>
5334 * void TestLevelSet<dim>::init_constraints()
5336 * constraints.clear ();
5337 * constraints.reinit (locally_relevant_dofs_LS);
5338 * DoFTools::make_hanging_node_constraints (dof_handler_LS, constraints);
5339 * constraints.close ();
5340 * constraints_disp_field.clear ();
5341 * constraints_disp_field.reinit (locally_relevant_dofs_LS);
5342 * DoFTools::make_hanging_node_constraints (dof_handler_LS, constraints_disp_field);
5343 * constraints_disp_field.close ();
5353 * template <int dim>
5354 * void TestLevelSet<dim>::process_solution(parallel::distributed::Triangulation<dim> &triangulation,
5355 * DoFHandler<dim> &dof_handler_LS,
5356 * PETScWrappers::MPI::Vector &solution)
5358 * Vector<double> difference_per_cell (triangulation.n_active_cells());
5364 * VectorTools::integrate_difference (dof_handler_LS,
5366 * InitialPhi<dim>(PROBLEM,sharpness),
5367 * difference_per_cell,
5368 * QGauss<dim>(degree_LS+3),
5369 * VectorTools::L1_norm);
5371 * double u_L1_error = difference_per_cell.l1_norm();
5372 * u_L1_error = std::sqrt(Utilities::MPI::sum(u_L1_error * u_L1_error, mpi_communicator));
5374 * VectorTools::integrate_difference (dof_handler_LS,
5376 * InitialPhi<dim>(PROBLEM,sharpness),
5377 * difference_per_cell,
5378 * QGauss<dim>(degree_LS+3),
5379 * VectorTools::L2_norm);
5380 * double u_L2_error = difference_per_cell.l2_norm();
5381 * u_L2_error = std::sqrt(Utilities::MPI::sum(u_L2_error * u_L2_error, mpi_communicator));
5383 * pcout << "L1 error: " << u_L1_error << std::endl;
5384 * pcout << "L2 error: " << u_L2_error << std::endl;
5388 * void TestLevelSet<dim>::output_results()
5390 * output_solution();
5394 * template <int dim>
5395 * void TestLevelSet<dim>::output_solution()
5397 * DataOut<dim> data_out;
5398 * data_out.attach_dof_handler(dof_handler_LS);
5399 * data_out.add_data_vector (locally_relevant_solution_phi, "phi");
5400 * data_out.build_patches();
5402 * const std::string filename = ("solution-" +
5403 * Utilities::int_to_string (output_number, 3) +
5405 * Utilities::int_to_string
5406 * (triangulation.locally_owned_subdomain(), 4));
5407 * std::ofstream output ((filename + ".vtu").c_str());
5408 * data_out.write_vtu (output);
5410 * if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
5412 * std::vector<std::string> filenames;
5413 * for (unsigned int i=0;
5414 * i<Utilities::MPI::n_mpi_processes(mpi_communicator);
5416 * filenames.push_back ("solution-" +
5417 * Utilities::int_to_string (output_number, 3) +
5419 * Utilities::int_to_string (i, 4) +
5422 * std::ofstream master_output ((filename + ".pvtu").c_str());
5423 * data_out.write_pvtu_record (master_output, filenames);
5427 * template <int dim>
5428 * void TestLevelSet<dim>::run()
5432 * ////////////////////
5433 * GENERAL PARAMETERS
5434 * ////////////////////
5439 * get_output = true;
5440 * output_number = 0;
5443 * output_time = 0.1;
5445 * PROBLEM=CIRCULAR_ROTATION;
5448 * PROBLEM=DIAGONAL_ADVECTION;
5452 * if (PROBLEM==CIRCULAR_ROTATION)
5453 * umax = std::sqrt(2)*numbers::PI;
5455 * umax = std::sqrt(2);
5459 * //////////////////////////////////
5460 * PARAMETERS FOR TRANSPORT PROBLEM
5461 * //////////////////////////////////
5464 * cK = 1.0; // compression constant
5465 * cE = 1.0; // entropy viscosity constant
5466 * sharpness_integer=1; //this will be multiplied by min_h
5469 * TRANSPORT_TIME_INTEGRATION=FORWARD_EULER;
5472 * TRANSPORT_TIME_INTEGRATION=SSP33;
5475 * ALGORITHM = "MPP_u1";
5478 * ALGORITHM = "NMPP_uH";
5481 * ALGORITHM = "MPP_uH";
5490 * if (PROBLEM==CIRCULAR_ROTATION || PROBLEM==DIAGONAL_ADVECTION)
5491 * GridGenerator::hyper_cube(triangulation);
5494 * GridGenerator::hyper_rectangle(triangulation, Point<dim>(0.0,0.0), Point<dim>(1.0,1.0), true);
5497 * triangulation.refine_global (n_refinement);
5510 * for Reconstruction of MATERIAL FIELDS
5513 * min_h = GridTools::minimal_cell_diameter(triangulation)/std::sqrt(dim)/degree;
5514 * eps=1*min_h; //For reconstruction of density in Navier Stokes
5515 * sharpness=sharpness_integer*min_h; //adjust value of sharpness (for init cond of phi)
5524 * time_step = cfl*min_h/umax;
5528 * //////////////////
5530 * //////////////////
5533 * LevelSetSolver<dim> level_set (degree_LS,degree_U,
5537 * TRANSPORT_TIME_INTEGRATION,
5539 * mpi_communicator);
5543 * ///////////////////
5545 * ///////////////////
5548 * initial_condition();
5551 * level_set.initial_condition(locally_relevant_solution_phi,
5552 * locally_relevant_solution_u,locally_relevant_solution_v);
5554 * level_set.initial_condition(locally_relevant_solution_phi,
5555 * locally_relevant_solution_u,locally_relevant_solution_v,locally_relevant_solution_w);
5559 * /////////////////////////////
5560 * BOUNDARY CONDITIONS FOR PHI
5561 * /////////////////////////////
5564 * get_boundary_values_phi(boundary_values_id_phi,boundary_values_phi);
5565 * level_set.set_boundary_conditions(boundary_values_id_phi,boundary_values_phi);
5569 * OUTPUT DATA REGARDING TIME STEPPING AND MESH
5572 * int dofs_LS = dof_handler_LS.n_dofs();
5573 * pcout << "Cfl: " << cfl << std::endl;
5574 * pcout << " Number of active cells: "
5575 * << triangulation.n_global_active_cells() << std::endl
5576 * << " Number of degrees of freedom: " << std::endl
5577 * << " LS: " << dofs_LS << std::endl;
5584 * timestep_number=0;
5586 * while (time<final_time)
5588 * timestep_number++;
5589 * if (time+time_step > final_time)
5591 * pcout << "FINAL TIME STEP... " << std::endl;
5592 * time_step = final_time-time;
5594 * pcout << "Time step " << timestep_number
5595 * << "\twith dt=" << time_step
5596 * << "\tat tn=" << time << std::endl;
5601 * GET VELOCITY // (NS or interpolate from a function) at current time tn
5605 * if (VARIABLE_VELOCITY)
5607 * get_interpolated_velocity();
5610 * SET VELOCITY TO LEVEL SET SOLVER
5613 * level_set.set_velocity(locally_relevant_solution_u,locally_relevant_solution_v);
5617 * ////////////////////////
5618 * GET LEVEL SET SOLUTION // (at tnp1)
5619 * ////////////////////////
5622 * level_set.nth_time_step();
5631 * time+=time_step; // time tnp1
5640 * if (get_output && time-(output_number)*output_time>=0)
5642 * level_set.get_unp1(locally_relevant_solution_phi);
5646 * pcout << "FINAL TIME T=" << time << std::endl;
5649 * int main(int argc, char *argv[])
5653 * using namespace dealii;
5654 * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
5655 * PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
5656 * deallog.depth_console (0);
5658 * unsigned int degree = 1;
5659 * TestLevelSet<2> multiphase(degree, degree);
5664 * catch (std::exception &exc)
5666 * std::cerr << std::endl << std::endl
5667 * << "----------------------------------------------------"
5669 * std::cerr << "Exception on processing: " << std::endl
5670 * << exc.what() << std::endl
5671 * << "Aborting!" << std::endl
5672 * << "----------------------------------------------------"
5678 * std::cerr << std::endl << std::endl
5679 * << "----------------------------------------------------"
5681 * std::cerr << "Unknown exception!" << std::endl
5682 * << "Aborting!" << std::endl
5683 * << "----------------------------------------------------"
5694<a name="ann-TestNavierStokes.cc"></a>
5695<h1>Annotated version of TestNavierStokes.cc</h1>
5701 * /* -----------------------------------------------------------------------------
5703 * * SPDX-License-Identifier: LGPL-2.1-or-later
5704 * * Copyright (C) 2016 Manuel Quezada de Luna
5706 * * This file is part of the deal.II code gallery.
5708 * * -----------------------------------------------------------------------------
5711 * #include <deal.II/base/quadrature_lib.h>
5712 * #include <deal.II/base/function.h>
5713 * #include <deal.II/lac/affine_constraints.h>
5714 * #include <deal.II/lac/vector.h>
5715 * #include <deal.II/lac/full_matrix.h>
5716 * #include <deal.II/lac/solver_cg.h>
5717 * #include <deal.II/lac/petsc_sparse_matrix.h>
5718 * #include <deal.II/lac/petsc_vector.h>
5719 * #include <deal.II/lac/petsc_solver.h>
5720 * #include <deal.II/lac/petsc_precondition.h>
5721 * #include <deal.II/grid/grid_generator.h>
5722 * #include <deal.II/grid/tria_accessor.h>
5723 * #include <deal.II/grid/tria_iterator.h>
5724 * #include <deal.II/dofs/dof_handler.h>
5725 * #include <deal.II/dofs/dof_accessor.h>
5726 * #include <deal.II/dofs/dof_tools.h>
5727 * #include <deal.II/fe/fe_values.h>
5728 * #include <deal.II/fe/fe_q.h>
5729 * #include <deal.II/numerics/vector_tools.h>
5730 * #include <deal.II/numerics/data_out.h>
5731 * #include <deal.II/numerics/error_estimator.h>
5732 * #include <deal.II/base/utilities.h>
5733 * #include <deal.II/base/conditional_ostream.h>
5734 * #include <deal.II/base/index_set.h>
5735 * #include <deal.II/lac/sparsity_tools.h>
5736 * #include <deal.II/distributed/tria.h>
5737 * #include <deal.II/distributed/grid_refinement.h>
5738 * #include <deal.II/lac/petsc_vector.h>
5739 * #include <deal.II/base/convergence_table.h>
5740 * #include <deal.II/base/timer.h>
5741 * #include <deal.II/base/parameter_handler.h>
5742 * #include <fstream>
5743 * #include <iostream>
5744 * #include <deal.II/grid/grid_tools.h>
5745 * #include <deal.II/fe/mapping_q.h>
5746 * #include <deal.II/base/function.h>
5748 * using namespace dealii;
5750 * #include "utilities_test_NS.cc"
5751 * #include "NavierStokesSolver.cc"
5755 * ///////////////////////////////////////////////////
5756 * /////////////////// MAIN CLASS ////////////////////
5757 * ///////////////////////////////////////////////////
5760 * template <int dim>
5761 * class TestNavierStokes
5764 * TestNavierStokes (const unsigned int degree_LS,
5765 * const unsigned int degree_U);
5766 * ~TestNavierStokes ();
5770 * void get_boundary_values_U(double t);
5771 * void fix_pressure();
5772 * void output_results();
5773 * void process_solution(const unsigned int cycle);
5775 * void initial_condition();
5776 * void init_constraints();
5778 * PETScWrappers::MPI::Vector locally_relevant_solution_rho;
5779 * PETScWrappers::MPI::Vector locally_relevant_solution_u;
5780 * PETScWrappers::MPI::Vector locally_relevant_solution_v;
5781 * PETScWrappers::MPI::Vector locally_relevant_solution_w;
5782 * PETScWrappers::MPI::Vector locally_relevant_solution_p;
5783 * PETScWrappers::MPI::Vector completely_distributed_solution_rho;
5784 * PETScWrappers::MPI::Vector completely_distributed_solution_u;
5785 * PETScWrappers::MPI::Vector completely_distributed_solution_v;
5786 * PETScWrappers::MPI::Vector completely_distributed_solution_w;
5787 * PETScWrappers::MPI::Vector completely_distributed_solution_p;
5789 * std::vector<unsigned int> boundary_values_id_u;
5790 * std::vector<unsigned int> boundary_values_id_v;
5791 * std::vector<unsigned int> boundary_values_id_w;
5792 * std::vector<double> boundary_values_u;
5793 * std::vector<double> boundary_values_v;
5794 * std::vector<double> boundary_values_w;
5801 * MPI_Comm mpi_communicator;
5802 * parallel::distributed::Triangulation<dim> triangulation;
5805 * DoFHandler<dim> dof_handler_LS;
5807 * IndexSet locally_owned_dofs_LS;
5808 * IndexSet locally_relevant_dofs_LS;
5811 * DoFHandler<dim> dof_handler_U;
5813 * IndexSet locally_owned_dofs_U;
5814 * IndexSet locally_relevant_dofs_U;
5816 * DoFHandler<dim> dof_handler_P;
5818 * IndexSet locally_owned_dofs_P;
5819 * IndexSet locally_relevant_dofs_P;
5821 * AffineConstraints<double> constraints;
5825 * TimerOutput timer;
5833 * double final_time;
5834 * unsigned int timestep_number;
5839 * unsigned int n_cycles;
5840 * unsigned int n_refinement;
5841 * unsigned int output_number;
5842 * double output_time;
5850 * ConditionalOStream pcout;
5851 * ConvergenceTable convergence_table;
5856 * template <int dim>
5857 * TestNavierStokes<dim>::TestNavierStokes (const unsigned int degree_LS,
5858 * const unsigned int degree_U)
5860 * mpi_communicator (MPI_COMM_WORLD),
5861 * triangulation (mpi_communicator,
5862 * typename Triangulation<dim>::MeshSmoothing
5863 * (Triangulation<dim>::smoothing_on_refinement |
5864 * Triangulation<dim>::smoothing_on_coarsening)),
5865 * degree_LS(degree_LS),
5866 * dof_handler_LS (triangulation),
5867 * fe_LS (degree_LS),
5868 * degree_U(degree_U),
5869 * dof_handler_U (triangulation),
5871 * dof_handler_P (triangulation),
5872 * fe_P (degree_U-1), //TODO: change this to be degree_Q-1
5875 * timer(std::cout, TimerOutput::summary, TimerOutput::wall_times),
5878 * pcout (std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator)== 0))
5881 * template <int dim>
5882 * TestNavierStokes<dim>::~TestNavierStokes ()
5884 * dof_handler_LS.clear ();
5885 * dof_handler_U.clear ();
5886 * dof_handler_P.clear ();
5891 * /////////////////////////////////////
5892 * /////////////// SETUP ///////////////
5893 * /////////////////////////////////////
5896 * template <int dim>
5897 * void TestNavierStokes<dim>::setup()
5904 * dof_handler_LS.distribute_dofs (fe_LS);
5905 * locally_owned_dofs_LS = dof_handler_LS.locally_owned_dofs ();
5906 * locally_relevant_dofs_LS = DoFTools::extract_locally_relevant_dofs (dof_handler_LS);
5912 * dof_handler_U.distribute_dofs (fe_U);
5913 * locally_owned_dofs_U = dof_handler_U.locally_owned_dofs ();
5914 * locally_relevant_dofs_U = DoFTools::extract_locally_relevant_dofs (dof_handler_U);
5920 * dof_handler_P.distribute_dofs (fe_P);
5921 * locally_owned_dofs_P = dof_handler_P.locally_owned_dofs ();
5922 * locally_relevant_dofs_P = DoFTools::extract_locally_relevant_dofs (dof_handler_P);
5923 * init_constraints();
5926 * init vectors for rho
5929 * locally_relevant_solution_rho.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
5930 * locally_relevant_solution_rho = 0;
5931 * completely_distributed_solution_rho.reinit(locally_owned_dofs_LS,mpi_communicator);
5934 * init vectors for u
5937 * locally_relevant_solution_u.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
5938 * locally_relevant_solution_u = 0;
5939 * completely_distributed_solution_u.reinit(locally_owned_dofs_U,mpi_communicator);
5942 * init vectors for v
5945 * locally_relevant_solution_v.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
5946 * locally_relevant_solution_v = 0;
5947 * completely_distributed_solution_v.reinit(locally_owned_dofs_U,mpi_communicator);
5950 * init vectors for w
5953 * locally_relevant_solution_w.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
5954 * locally_relevant_solution_w = 0;
5955 * completely_distributed_solution_w.reinit(locally_owned_dofs_U,mpi_communicator);
5958 * init vectors for p
5961 * locally_relevant_solution_p.reinit(locally_owned_dofs_P,locally_relevant_dofs_P,mpi_communicator);
5962 * locally_relevant_solution_p = 0;
5963 * completely_distributed_solution_p.reinit(locally_owned_dofs_P,mpi_communicator);
5966 * template <int dim>
5967 * void TestNavierStokes<dim>::initial_condition()
5972 * Initial conditions
5973 * init condition for rho
5976 * completely_distributed_solution_rho = 0;
5977 * VectorTools::interpolate(dof_handler_LS,
5978 * RhoFunction<dim>(0),
5979 * completely_distributed_solution_rho);
5980 * constraints.distribute (completely_distributed_solution_rho);
5981 * locally_relevant_solution_rho = completely_distributed_solution_rho;
5984 * init condition for u
5987 * completely_distributed_solution_u = 0;
5988 * VectorTools::interpolate(dof_handler_U,
5989 * ExactSolution_and_BC_U<dim>(0,0),
5990 * completely_distributed_solution_u);
5991 * constraints.distribute (completely_distributed_solution_u);
5992 * locally_relevant_solution_u = completely_distributed_solution_u;
5995 * init condition for v
5998 * completely_distributed_solution_v = 0;
5999 * VectorTools::interpolate(dof_handler_U,
6000 * ExactSolution_and_BC_U<dim>(0,1),
6001 * completely_distributed_solution_v);
6002 * constraints.distribute (completely_distributed_solution_v);
6003 * locally_relevant_solution_v = completely_distributed_solution_v;
6006 * init condition for w
6011 * completely_distributed_solution_w = 0;
6012 * VectorTools::interpolate(dof_handler_U,
6013 * ExactSolution_and_BC_U<dim>(0,2),
6014 * completely_distributed_solution_w);
6015 * constraints.distribute (completely_distributed_solution_w);
6016 * locally_relevant_solution_w = completely_distributed_solution_w;
6020 * init condition for p
6023 * completely_distributed_solution_p = 0;
6024 * VectorTools::interpolate(dof_handler_P,
6025 * ExactSolution_p<dim>(0),
6026 * completely_distributed_solution_p);
6027 * constraints.distribute (completely_distributed_solution_p);
6028 * locally_relevant_solution_p = completely_distributed_solution_p;
6031 * template <int dim>
6032 * void TestNavierStokes<dim>::init_constraints()
6034 * constraints.clear ();
6035 * constraints.reinit (locally_relevant_dofs_LS);
6036 * DoFTools::make_hanging_node_constraints (dof_handler_LS, constraints);
6037 * constraints.close ();
6041 * void TestNavierStokes<dim>::fix_pressure()
6045 * fix the constant in the pressure
6048 * completely_distributed_solution_p = locally_relevant_solution_p;
6049 * double mean_value = VectorTools::compute_mean_value(dof_handler_P,
6051 * locally_relevant_solution_p,
6054 * completely_distributed_solution_p.add(-mean_value+std::sin(1)*(std::cos(time)-cos(1+time)));
6056 * completely_distributed_solution_p.add(-mean_value+8*std::pow(std::sin(0.5),3)*std::sin(1.5+time));
6057 * locally_relevant_solution_p = completely_distributed_solution_p;
6060 * template <int dim>
6061 * void TestNavierStokes<dim>::output_results ()
6063 * DataOut<dim> data_out;
6064 * data_out.attach_dof_handler (dof_handler_U);
6065 * data_out.add_data_vector (locally_relevant_solution_u, "u");
6066 * data_out.add_data_vector (locally_relevant_solution_v, "v");
6067 * if (dim==3) data_out.add_data_vector (locally_relevant_solution_w, "w");
6069 * Vector<float> subdomain (triangulation.n_active_cells());
6070 * for (unsigned int i=0; i<subdomain.size(); ++i)
6071 * subdomain(i) = triangulation.locally_owned_subdomain();
6072 * data_out.add_data_vector (subdomain, "subdomain");
6074 * data_out.build_patches ();
6076 * const std::string filename = ("solution-" +
6077 * Utilities::int_to_string (output_number, 3) +
6079 * Utilities::int_to_string
6080 * (triangulation.locally_owned_subdomain(), 4));
6081 * std::ofstream output ((filename + ".vtu").c_str());
6082 * data_out.write_vtu (output);
6084 * if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
6086 * std::vector<std::string> filenames;
6087 * for (unsigned int i=0;
6088 * i<Utilities::MPI::n_mpi_processes(mpi_communicator);
6090 * filenames.push_back ("solution-" +
6091 * Utilities::int_to_string (output_number, 3) +
6093 * Utilities::int_to_string (i, 4) +
6096 * std::ofstream master_output ((filename + ".pvtu").c_str());
6097 * data_out.write_pvtu_record (master_output, filenames);
6102 * template <int dim>
6103 * void TestNavierStokes<dim>::process_solution(const unsigned int cycle)
6105 * Vector<double> difference_per_cell (triangulation.n_active_cells());
6111 * VectorTools::integrate_difference (dof_handler_U,
6112 * locally_relevant_solution_u,
6113 * ExactSolution_and_BC_U<dim>(time,0),
6114 * difference_per_cell,
6115 * QGauss<dim>(degree_U+1),
6116 * VectorTools::L2_norm);
6117 * double u_L2_error = difference_per_cell.l2_norm();
6119 * std::sqrt(Utilities::MPI::sum(u_L2_error * u_L2_error, mpi_communicator));
6120 * VectorTools::integrate_difference (dof_handler_U,
6121 * locally_relevant_solution_u,
6122 * ExactSolution_and_BC_U<dim>(time,0),
6123 * difference_per_cell,
6124 * QGauss<dim>(degree_U+1),
6125 * VectorTools::H1_norm);
6126 * double u_H1_error = difference_per_cell.l2_norm();
6128 * std::sqrt(Utilities::MPI::sum(u_H1_error * u_H1_error, mpi_communicator));
6134 * VectorTools::integrate_difference (dof_handler_U,
6135 * locally_relevant_solution_v,
6136 * ExactSolution_and_BC_U<dim>(time,1),
6137 * difference_per_cell,
6138 * QGauss<dim>(degree_U+1),
6139 * VectorTools::L2_norm);
6140 * double v_L2_error = difference_per_cell.l2_norm();
6142 * std::sqrt(Utilities::MPI::sum(v_L2_error * v_L2_error,
6143 * mpi_communicator));
6144 * VectorTools::integrate_difference (dof_handler_U,
6145 * locally_relevant_solution_v,
6146 * ExactSolution_and_BC_U<dim>(time,1),
6147 * difference_per_cell,
6148 * QGauss<dim>(degree_U+1),
6149 * VectorTools::H1_norm);
6150 * double v_H1_error = difference_per_cell.l2_norm();
6152 * std::sqrt(Utilities::MPI::sum(v_H1_error *
6153 * v_H1_error, mpi_communicator));
6159 * double w_L2_error = 0;
6160 * double w_H1_error = 0;
6163 * VectorTools::integrate_difference (dof_handler_U,
6164 * locally_relevant_solution_w,
6165 * ExactSolution_and_BC_U<dim>(time,2),
6166 * difference_per_cell,
6167 * QGauss<dim>(degree_U+1),
6168 * VectorTools::L2_norm);
6169 * w_L2_error = difference_per_cell.l2_norm();
6171 * std::sqrt(Utilities::MPI::sum(w_L2_error * w_L2_error,
6172 * mpi_communicator));
6173 * VectorTools::integrate_difference (dof_handler_U,
6174 * locally_relevant_solution_w,
6175 * ExactSolution_and_BC_U<dim>(time,2),
6176 * difference_per_cell,
6177 * QGauss<dim>(degree_U+1),
6178 * VectorTools::H1_norm);
6179 * w_H1_error = difference_per_cell.l2_norm();
6181 * std::sqrt(Utilities::MPI::sum(w_H1_error *
6182 * w_H1_error, mpi_communicator));
6189 * VectorTools::integrate_difference (dof_handler_P,
6190 * locally_relevant_solution_p,
6191 * ExactSolution_p<dim>(time),
6192 * difference_per_cell,
6193 * QGauss<dim>(degree_U+1),
6194 * VectorTools::L2_norm);
6195 * double p_L2_error = difference_per_cell.l2_norm();
6197 * std::sqrt(Utilities::MPI::sum(p_L2_error * p_L2_error,
6198 * mpi_communicator));
6199 * VectorTools::integrate_difference (dof_handler_P,
6200 * locally_relevant_solution_p,
6201 * ExactSolution_p<dim>(time),
6202 * difference_per_cell,
6203 * QGauss<dim>(degree_U+1),
6204 * VectorTools::H1_norm);
6205 * double p_H1_error = difference_per_cell.l2_norm();
6207 * std::sqrt(Utilities::MPI::sum(p_H1_error * p_H1_error,
6208 * mpi_communicator));
6210 * const unsigned int n_active_cells=triangulation.n_active_cells();
6211 * const unsigned int n_dofs_U=dof_handler_U.n_dofs();
6212 * const unsigned int n_dofs_P=dof_handler_P.n_dofs();
6214 * convergence_table.add_value("cycle", cycle);
6215 * convergence_table.add_value("cells", n_active_cells);
6216 * convergence_table.add_value("dofs_U", n_dofs_U);
6217 * convergence_table.add_value("dofs_P", n_dofs_P);
6218 * convergence_table.add_value("dt", time_step);
6219 * convergence_table.add_value("u L2", u_L2_error);
6220 * convergence_table.add_value("u H1", u_H1_error);
6221 * convergence_table.add_value("v L2", v_L2_error);
6222 * convergence_table.add_value("v H1", v_H1_error);
6225 * convergence_table.add_value("w L2", w_L2_error);
6226 * convergence_table.add_value("w H1", w_H1_error);
6228 * convergence_table.add_value("p L2", p_L2_error);
6229 * convergence_table.add_value("p H1", p_H1_error);
6232 * template <int dim>
6233 * void TestNavierStokes<dim>::get_boundary_values_U(double t)
6235 * std::map<unsigned int, double> map_boundary_values_u;
6236 * std::map<unsigned int, double> map_boundary_values_v;
6238 * VectorTools::interpolate_boundary_values (dof_handler_U,0,ExactSolution_and_BC_U<dim>(t,0),map_boundary_values_u);
6239 * VectorTools::interpolate_boundary_values (dof_handler_U,0,ExactSolution_and_BC_U<dim>(t,1),map_boundary_values_v);
6241 * boundary_values_id_u.resize(map_boundary_values_u.size());
6242 * boundary_values_id_v.resize(map_boundary_values_v.size());
6243 * boundary_values_u.resize(map_boundary_values_u.size());
6244 * boundary_values_v.resize(map_boundary_values_v.size());
6245 * std::map<unsigned int,double>::const_iterator boundary_value_u =map_boundary_values_u.begin();
6246 * std::map<unsigned int,double>::const_iterator boundary_value_v =map_boundary_values_v.begin();
6249 * std::map<unsigned int, double> map_boundary_values_w;
6250 * VectorTools::interpolate_boundary_values (dof_handler_U,0,ExactSolution_and_BC_U<dim>(t,2),map_boundary_values_w);
6251 * boundary_values_id_w.resize(map_boundary_values_w.size());
6252 * boundary_values_w.resize(map_boundary_values_w.size());
6253 * std::map<unsigned int,double>::const_iterator boundary_value_w =map_boundary_values_w.begin();
6254 * for (int i=0; boundary_value_w !=map_boundary_values_w.end(); ++boundary_value_w, ++i)
6256 * boundary_values_id_w[i]=boundary_value_w->first;
6257 * boundary_values_w[i]=boundary_value_w->second;
6260 * for (int i=0; boundary_value_u !=map_boundary_values_u.end(); ++boundary_value_u, ++i)
6262 * boundary_values_id_u[i]=boundary_value_u->first;
6263 * boundary_values_u[i]=boundary_value_u->second;
6265 * for (int i=0; boundary_value_v !=map_boundary_values_v.end(); ++boundary_value_v, ++i)
6267 * boundary_values_id_v[i]=boundary_value_v->first;
6268 * boundary_values_v[i]=boundary_value_v->second;
6272 * template <int dim>
6273 * void TestNavierStokes<dim>::run()
6275 * if (Utilities::MPI::this_mpi_process(mpi_communicator)== 0)
6277 * std::cout << "***** CONVERGENCE TEST FOR NS *****" << std::endl;
6278 * std::cout << "DEGREE LS: " << degree_LS << std::endl;
6279 * std::cout << "DEGREE U: " << degree_U << std::endl;
6283 * PARAMETERS FOR THE NAVIER STOKES PROBLEM
6290 * ForceTerms<dim> force_function;
6291 * RhoFunction<dim> rho_function;
6292 * NuFunction<dim> nu_function;
6296 * bool get_output = false;
6297 * bool get_error = true;
6300 * for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
6304 * GridGenerator::hyper_cube (triangulation);
6305 * triangulation.refine_global (n_refinement);
6307 * initial_condition();
6311 * triangulation.refine_global(1);
6313 * initial_condition();
6323 * NavierStokesSolver<dim> navier_stokes (degree_LS,
6331 * mpi_communicator);
6334 * set INITIAL CONDITION within TRANSPORT PROBLEM
6338 * navier_stokes.initial_condition(locally_relevant_solution_rho,
6339 * locally_relevant_solution_u,
6340 * locally_relevant_solution_v,
6341 * locally_relevant_solution_p);
6343 * navier_stokes.initial_condition(locally_relevant_solution_rho,
6344 * locally_relevant_solution_u,
6345 * locally_relevant_solution_v,
6346 * locally_relevant_solution_w,
6347 * locally_relevant_solution_p);
6349 * pcout << "Cycle " << cycle << ':
' << std::endl;
6350 * pcout << " Cycle " << cycle
6351 * << " Number of active cells: "
6352 * << triangulation.n_global_active_cells() << std::endl
6353 * << " Number of degrees of freedom (velocity): "
6354 * << dof_handler_U.n_dofs() << std::endl
6355 * << " min h=" << GridTools::minimal_cell_diameter(triangulation)/std::sqrt(2)/degree_U
6363 * timestep_number=0;
6365 * double time_step_backup=time_step;
6366 * while (time<final_time)
6368 * timestep_number++;
6376 * if (time+time_step > final_time-1E-10)
6378 * pcout << "FINAL TIME STEP..." << std::endl;
6379 * time_step_backup=time_step;
6380 * time_step=final_time-time;
6382 * pcout << "Time step " << timestep_number
6383 * << "\twith dt=" << time_step
6384 * << "\tat tn=" << time
6393 * force_function.set_time(time+time_step);
6396 * /////////////////////////////
6397 * DENSITY AND VISCOSITY FIELD
6398 * /////////////////////////////
6401 * rho_function.set_time(time+time_step);
6402 * nu_function.set_time(time+time_step);
6405 * /////////////////////
6406 * BOUNDARY CONDITIONS
6407 * /////////////////////
6410 * get_boundary_values_U(time+time_step);
6411 * if (dim==2) navier_stokes.set_boundary_conditions(boundary_values_id_u, boundary_values_id_v,
6412 * boundary_values_u, boundary_values_v);
6413 * else navier_stokes.set_boundary_conditions(boundary_values_id_u,
6414 * boundary_values_id_v,
6415 * boundary_values_id_w,
6416 * boundary_values_u, boundary_values_v, boundary_values_w);
6424 * navier_stokes.nth_time_step();
6426 * navier_stokes.get_velocity(locally_relevant_solution_u,locally_relevant_solution_v);
6428 * navier_stokes.get_velocity(locally_relevant_solution_u,
6429 * locally_relevant_solution_v,
6430 * locally_relevant_solution_w);
6431 * navier_stokes.get_pressure(locally_relevant_solution_p);
6458 * if (get_output && time-(output_number)*output_time>=1E-10)
6461 * pcout << "FINAL TIME: " << time << std::endl;
6462 * time_step=time_step_backup;
6464 * process_solution(cycle);
6468 * convergence_table.set_precision("u L2", 2);
6469 * convergence_table.set_precision("u H1", 2);
6470 * convergence_table.set_scientific("u L2",true);
6471 * convergence_table.set_scientific("u H1",true);
6473 * convergence_table.set_precision("v L2", 2);
6474 * convergence_table.set_precision("v H1", 2);
6475 * convergence_table.set_scientific("v L2",true);
6476 * convergence_table.set_scientific("v H1",true);
6480 * convergence_table.set_precision("w L2", 2);
6481 * convergence_table.set_precision("w H1", 2);
6482 * convergence_table.set_scientific("w L2",true);
6483 * convergence_table.set_scientific("w H1",true);
6486 * convergence_table.set_precision("p L2", 2);
6487 * convergence_table.set_precision("p H1", 2);
6488 * convergence_table.set_scientific("p L2",true);
6489 * convergence_table.set_scientific("p H1",true);
6491 * convergence_table.set_tex_format("cells","r");
6492 * convergence_table.set_tex_format("dofs_U","r");
6493 * convergence_table.set_tex_format("dofs_P","r");
6494 * convergence_table.set_tex_format("dt","r");
6496 * if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
6498 * std::cout << std::endl;
6499 * convergence_table.write_text(std::cout);
6505 * int main(int argc, char *argv[])
6509 * using namespace dealii;
6510 * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
6511 * PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
6512 * deallog.depth_console (0);
6515 * unsigned int degree_LS = 1;
6516 * unsigned int degree_U = 2;
6517 * TestNavierStokes<2> test_navier_stokes(degree_LS, degree_U);
6518 * test_navier_stokes.run();
6525 * catch (std::exception &exc)
6527 * std::cerr << std::endl << std::endl
6528 * << "----------------------------------------------------"
6530 * std::cerr << "Exception on processing: " << std::endl
6531 * << exc.what() << std::endl
6532 * << "Aborting!" << std::endl
6533 * << "----------------------------------------------------"
6540 * std::cerr << std::endl << std::endl
6541 * << "----------------------------------------------------"
6543 * std::cerr << "Unknown exception!" << std::endl
6544 * << "Aborting!" << std::endl
6545 * << "----------------------------------------------------"
6555<a name="ann-clean.sh"></a>
6556<h1>Annotated version of clean.sh</h1>
6558rm -rf CMakeFiles CMakeCache.txt Makefile cmake_install.cmake *~
6559rm -f MultiPhase TestLevelSet TestNavierStokes
6566<a name="ann-utilities.cc"></a>
6567<h1>Annotated version of utilities.cc</h1>
6573 * /* -----------------------------------------------------------------------------
6575 * * SPDX-License-Identifier: LGPL-2.1-or-later
6576 * * Copyright (C) 2016 Manuel Quezada de Luna
6578 * * This file is part of the deal.II code gallery.
6580 * * -----------------------------------------------------------------------------
6585 * /////////////////////////////////////////////////
6586 * ////////////////// INITIAL PHI //////////////////
6587 * /////////////////////////////////////////////////
6590 * template <int dim>
6591 * class InitialPhi : public Function <dim>
6594 * InitialPhi (unsigned int PROBLEM, double sharpness=0.005) : Function<dim>(),
6595 * sharpness(sharpness),
6596 * PROBLEM(PROBLEM) {}
6597 * virtual double value (const Point<dim> &p, const unsigned int component=0) const override;
6599 * unsigned int PROBLEM;
6601 * template <int dim>
6602 * double InitialPhi<dim>::value (const Point<dim> &p,
6603 * const unsigned int) const
6607 * double pi=numbers::PI;
6609 * if (PROBLEM==FILLING_TANK)
6610 * return 0.5*(-std::tanh((y-0.3)/sharpness)*std::tanh((y-0.35)/sharpness)+1)
6611 * *(-std::tanh((x-0.02)/sharpness)+1)-1;
6612 * else if (PROBLEM==BREAKING_DAM)
6613 * return 0.5*(-std::tanh((x-0.35)/sharpness)*std::tanh((x-0.65)/sharpness)+1)
6614 * *(1-std::tanh((y-0.35)/sharpness))-1;
6615 * else if (PROBLEM==FALLING_DROP)
6620 * double r = std::sqrt(std::pow(x-x0,2)+std::pow(y-y0,2));
6621 * return 1-(std::tanh((r-r0)/sharpness)+std::tanh((y-0.3)/sharpness));
6623 * else if (PROBLEM==SMALL_WAVE_PERTURBATION)
6625 * double wave = 0.1*std::sin(pi*x)+0.25;
6626 * return -std::tanh((y-wave)/sharpness);
6630 * std::cout << "Error in type of PROBLEM" << std::endl;
6637 * ///////////////////////////////////////////////////
6638 * ////////////////// FORCE TERMS ///// //////////////
6639 * ///////////////////////////////////////////////////
6642 * template <int dim>
6643 * class ForceTerms : public Functions::ConstantFunction <dim>
6646 * ForceTerms (const std::vector<double> values) : Functions::ConstantFunction<dim>(values) {}
6651 * /////////////////////////////////////////////////
6652 * ////////////////// BOUNDARY PHI /////////////////
6653 * /////////////////////////////////////////////////
6656 * template <int dim>
6657 * class BoundaryPhi : public Functions::ConstantFunction <dim>
6660 * BoundaryPhi (const double value, const unsigned int n_components=1) : Functions::ConstantFunction<dim>(value,n_components) {}
6665 * //////////////////////////////////////////////////////
6666 * ////////////////// BOUNDARY VELOCITY /////////////////
6667 * //////////////////////////////////////////////////////
6670 * template <int dim>
6671 * class BoundaryU : public Function <dim>
6674 * BoundaryU (unsigned int PROBLEM, double t=0) : Function<dim>(), PROBLEM(PROBLEM) {this->set_time(t);}
6675 * virtual double value (const Point<dim> &p, const unsigned int component=0) const override;
6678 * template <int dim>
6679 * double BoundaryU<dim>::value (const Point<dim> &p, const unsigned int) const
6683 * //////////////////
6685 * //////////////////
6686 * boundary for filling the tank (inlet)
6692 * if (PROBLEM==FILLING_TANK)
6694 * if (x==0 && y>=0.3 && y<=0.35)
6701 * std::cout << "Error in PROBLEM definition" << std::endl;
6706 * template <int dim>
6707 * class BoundaryV : public Function <dim>
6710 * BoundaryV (unsigned int PROBLEM, double t=0) : Function<dim>(), PROBLEM(PROBLEM) {this->set_time(t);}
6711 * virtual double value (const Point<dim> &p, const unsigned int component=0) const override;
6712 * unsigned int PROBLEM;
6714 * template <int dim>
6715 * double BoundaryV<dim>::value (const Point<dim> &p, const unsigned int) const
6719 * boundary for filling the tank (outlet)
6724 * double return_value = 0;
6726 * if (PROBLEM==FILLING_TANK)
6728 * if (y==0.4 && x>=0.3 && x<=0.35)
6729 * return_value = 0.25;
6731 * return return_value;
6736 * ///////////////////////////////////////////////////
6737 * ///////////////// POST-PROCESSING /////////////////
6738 * ///////////////////////////////////////////////////
6741 * template <int dim>
6742 * class Postprocessor : public DataPostprocessorScalar <dim>
6745 * Postprocessor(double eps, double rho_air, double rho_fluid)
6747 * DataPostprocessorScalar<dim>("Density",update_values)
6750 * this->rho_air=rho_air;
6751 * this->rho_fluid=rho_fluid;
6756 * evaluate_scalar_field (const DataPostprocessorInputs::Scalar<dim> &input_data,
6757 * std::vector<Vector<double> > &computed_quantities) const override;
6765 * template <int dim>
6767 * Postprocessor<dim>::
6768 * evaluate_scalar_field (const DataPostprocessorInputs::Scalar<dim> &input_data,
6769 * std::vector<Vector<double> > &computed_quantities) const
6771 * const unsigned int n_quadrature_points = input_data.solution_values.size();
6772 * for (unsigned int q=0; q<n_quadrature_points; ++q)
6776 * double phi_value=input_data.solution_values[q];
6777 * if (phi_value > eps)
6779 * else if (phi_value < -eps)
6783 * rho_value = rho_fluid*(1+H)/2. + rho_air*(1-H)/2.;
6784 * computed_quantities[q] = rho_value;
6791<a name="ann-utilities_test_LS.cc"></a>
6792<h1>Annotated version of utilities_test_LS.cc</h1>
6798 * /* -----------------------------------------------------------------------------
6800 * * SPDX-License-Identifier: LGPL-2.1-or-later
6801 * * Copyright (C) 2016 Manuel Quezada de Luna
6803 * * This file is part of the deal.II code gallery.
6805 * * -----------------------------------------------------------------------------
6810 * ///////////////////////////////////////////////////
6811 * ////////////////// INITIAL PHI //////////////////
6812 * ///////////////////////////////////////////////////
6815 * template <int dim>
6816 * class InitialPhi : public Function <dim>
6819 * InitialPhi (unsigned int PROBLEM, double sharpness=0.005) : Function<dim>(),
6820 * sharpness(sharpness),
6821 * PROBLEM(PROBLEM) {}
6822 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6824 * unsigned int PROBLEM;
6826 * template <int dim>
6827 * double InitialPhi<dim>::value (const Point<dim> &p,
6828 * const unsigned int) const
6832 * double return_value = -1.;
6834 * if (PROBLEM==CIRCULAR_ROTATION)
6839 * double r = std::sqrt(std::pow(x-x0,2)+std::pow(y-y0,2));
6840 * return_value = -std::tanh((r-r0)/sharpness);
6842 * else // (PROBLEM==DIAGONAL_ADVECTION)
6849 * r = std::sqrt(std::pow(x-x0,2)+std::pow(y-y0,2));
6854 * r = std::sqrt(std::pow(x-x0,2)+std::pow(y-y0,2)+std::pow(z-z0,2));
6856 * return_value = -std::tanh((r-r0)/sharpness);
6858 * return return_value;
6863 * /////////////////////////////////////////////////
6864 * ////////////////// BOUNDARY PHI /////////////////
6865 * /////////////////////////////////////////////////
6868 * template <int dim>
6869 * class BoundaryPhi : public Function <dim>
6872 * BoundaryPhi (double t=0)
6875 * {this->set_time(t);}
6876 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6879 * template <int dim>
6880 * double BoundaryPhi<dim>::value (const Point<dim> &, const unsigned int) const
6887 * ///////////////////////////////////////////////////
6888 * ////////////////// EXACT VELOCITY /////////////////
6889 * ///////////////////////////////////////////////////
6892 * template <int dim>
6893 * class ExactU : public Function <dim>
6896 * ExactU (unsigned int PROBLEM, double time=0) : Function<dim>(), PROBLEM(PROBLEM), time(time) {}
6897 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6898 * void set_time(double time) {this->time=time;};
6903 * template <int dim>
6904 * double ExactU<dim>::value (const Point<dim> &p, const unsigned int) const
6906 * if (PROBLEM==CIRCULAR_ROTATION)
6907 * return -2*numbers::PI*(p[1]-0.5);
6908 * else // (PROBLEM==DIAGONAL_ADVECTION)
6912 * template <int dim>
6913 * class ExactV : public Function <dim>
6916 * ExactV (unsigned int PROBLEM, double time=0) : Function<dim>(), PROBLEM(PROBLEM), time(time) {}
6917 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6918 * void set_time(double time) {this->time=time;};
6919 * unsigned int PROBLEM;
6923 * template <int dim>
6924 * double ExactV<dim>::value (const Point<dim> &p, const unsigned int) const
6926 * if (PROBLEM==CIRCULAR_ROTATION)
6927 * return 2*numbers::PI*(p[0]-0.5);
6928 * else // (PROBLEM==DIAGONAL_ADVECTION)
6932 * template <int dim>
6933 * class ExactW : public Function <dim>
6936 * ExactW (unsigned int PROBLEM, double time=0) : Function<dim>(), PROBLEM(PROBLEM), time(time) {}
6937 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6938 * void set_time(double time) {this->time=time;};
6939 * unsigned int PROBLEM;
6943 * template <int dim>
6944 * double ExactW<dim>::value (const Point<dim> &, const unsigned int) const
6948 * PROBLEM = 3D_DIAGONAL_ADVECTION
6957<a name="ann-utilities_test_NS.cc"></a>
6958<h1>Annotated version of utilities_test_NS.cc</h1>
6964 * /* -----------------------------------------------------------------------------
6966 * * SPDX-License-Identifier: LGPL-2.1-or-later
6967 * * Copyright (C) 2016 Manuel Quezada de Luna
6969 * * This file is part of the deal.II code gallery.
6971 * * -----------------------------------------------------------------------------
6976 * ///////////////////////////////////////////////////
6977 * ////////// EXACT SOLUTION RHO TO TEST NS //////////
6978 * ///////////////////////////////////////////////////
6981 * template <int dim>
6982 * class RhoFunction : public Function <dim>
6985 * RhoFunction (double t=0) : Function<dim>() {this->set_time(t);}
6986 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6988 * template <int dim>
6989 * double RhoFunction<dim>::value (const Point<dim> &p,
6990 * const unsigned int) const
6992 * double t = this->get_time();
6993 * double return_value = 0;
6995 * return_value = std::pow(std::sin(p[0]+p[1]+t),2)+1;
6997 * return_value = std::pow(std::sin(p[0]+p[1]+p[2]+t),2)+1;
6998 * return return_value;
7001 * template <int dim>
7002 * class NuFunction : public Function <dim>
7005 * NuFunction (double t=0) : Function<dim>() {this->set_time(t);}
7006 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
7008 * template <int dim>
7009 * double NuFunction<dim>::value (const Point<dim> &, const unsigned int) const
7016 * //////////////////////////////////////////////////////////////
7017 * ///////////////// EXACT SOLUTION U to TEST NS ////////////////
7018 * //////////////////////////////////////////////////////////////
7021 * template <int dim>
7022 * class ExactSolution_and_BC_U : public Function <dim>
7025 * ExactSolution_and_BC_U (double t=0, int field=0)
7030 * this->set_time(t);
7032 * virtual double value (const Point<dim> &p, const unsigned int component=1) const;
7033 * virtual Tensor<1,dim> gradient (const Point<dim> &p, const unsigned int component=1) const;
7034 * virtual void set_field(int field) {this->field=field;}
7036 * unsigned int type_simulation;
7038 * template <int dim>
7039 * double ExactSolution_and_BC_U<dim>::value (const Point<dim> &p,
7040 * const unsigned int) const
7042 * double t = this->get_time();
7043 * double return_value = 0;
7044 * double Pi = numbers::PI;
7051 * return_value = std::sin(x)*std::sin(y+t);
7053 * return_value = std::cos(x)*std::cos(y+t);
7058 * return_value = std::cos(t)*std::cos(Pi*y)*std::cos(Pi*z)*std::sin(Pi*x);
7059 * else if (field == 1)
7060 * return_value = std::cos(t)*std::cos(Pi*x)*std::cos(Pi*z)*std::sin(Pi*y);
7062 * return_value = -2*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::sin(Pi*z);
7064 * return return_value;
7066 * template <int dim>
7067 * Tensor<1,dim> ExactSolution_and_BC_U<dim>::gradient (const Point<dim> &p,
7068 * const unsigned int) const
7072 * THIS IS USED JUST FOR TESTING NS
7075 * Tensor<1,dim> return_value;
7076 * double t = this->get_time();
7077 * double Pi = numbers::PI;
7084 * return_value[0] = std::cos(x)*std::sin(y+t);
7085 * return_value[1] = std::sin(x)*std::cos(y+t);
7089 * return_value[0] = -std::sin(x)*std::cos(y+t);
7090 * return_value[1] = -std::cos(x)*std::sin(y+t);
7097 * return_value[0] = Pi*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::cos(Pi*z);
7098 * return_value[1] = -(Pi*std::cos(t)*std::cos(Pi*z)*std::sin(Pi*x)*std::sin(Pi*y));
7099 * return_value[2] = -(Pi*std::cos(t)*std::cos(Pi*y)*std::sin(Pi*x)*std::sin(Pi*z));
7101 * else if (field == 1)
7103 * return_value[0] = -(Pi*std::cos(t)*std::cos(Pi*z)*std::sin(Pi*x)*std::sin(Pi*y));
7104 * return_value[1] = Pi*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::cos(Pi*z);
7105 * return_value[2] = -(Pi*std::cos(t)*std::cos(Pi*x)*std::sin(Pi*y)*std::sin(Pi*z));
7109 * return_value[0] = 2*Pi*std::cos(t)*std::cos(Pi*y)*std::sin(Pi*x)*std::sin(Pi*z);
7110 * return_value[1] = 2*Pi*std::cos(t)*std::cos(Pi*x)*std::sin(Pi*y)*std::sin(Pi*z);
7111 * return_value[2] = -2*Pi*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::cos(Pi*z);
7114 * return return_value;
7119 * ///////////////////////////////////////////////////
7120 * ///////// EXACT SOLUTION FOR p TO TEST NS /////////
7121 * ///////////////////////////////////////////////////
7124 * template <int dim>
7125 * class ExactSolution_p : public Function <dim>
7128 * ExactSolution_p (double t=0) : Function<dim>() {this->set_time(t);}
7129 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
7130 * virtual Tensor<1,dim> gradient (const Point<dim> &p, const unsigned int component = 0) const;
7133 * template <int dim>
7134 * double ExactSolution_p<dim>::value (const Point<dim> &p, const unsigned int) const
7136 * double t = this->get_time();
7137 * double return_value = 0;
7139 * return_value = std::cos(p[0])*std::sin(p[1]+t);
7141 * return_value = std::sin(p[0]+p[1]+p[2]+t);
7142 * return return_value;
7145 * template <int dim>
7146 * Tensor<1,dim> ExactSolution_p<dim>::gradient (const Point<dim> &p, const unsigned int) const
7148 * Tensor<1,dim> return_value;
7149 * double t = this->get_time();
7152 * return_value[0] = -std::sin(p[0])*std::sin(p[1]+t);
7153 * return_value[1] = std::cos(p[0])*std::cos(p[1]+t);
7157 * return_value[0] = std::cos(t+p[0]+p[1]+p[2]);
7158 * return_value[1] = std::cos(t+p[0]+p[1]+p[2]);
7159 * return_value[2] = std::cos(t+p[0]+p[1]+p[2]);
7161 * return return_value;
7166 * //////////////////////////////////////////////////////////////
7167 * ////////////////// FORCE TERMS to TEST NS ////////////////////
7168 * //////////////////////////////////////////////////////////////
7171 * template <int dim>
7172 * class ForceTerms : public Function <dim>
7175 * ForceTerms (double t=0)
7179 * this->set_time(t);
7182 * virtual void vector_value (const Point<dim> &p, Vector<double> &values) const;
7186 * template <int dim>
7187 * void ForceTerms<dim>::vector_value (const Point<dim> &p, Vector<double> &values) const
7192 * double t = this->get_time();
7193 * double Pi = numbers::PI;
7202 * values[0] = std::cos(t+y)*std::sin(x)*(1+std::pow(std::sin(t+x+y),2)) // time derivative
7203 * +2*nu*std::sin(x)*std::sin(t+y) // viscosity
7204 * +std::cos(x)*std::sin(x)*(1+std::pow(std::sin(t+x+y),2)) // non-linearity
7205 * -std::sin(x)*std::sin(y+t); // pressure
7211 * values[1] = -(std::cos(x)*std::sin(t+y)*(1+std::pow(std::sin(t+x+y),2))) // time derivative
7212 * +2*nu*std::cos(x)*std::cos(t+y) // viscosity
7213 * -(std::sin(2*(t+y))*(1+std::pow(std::sin(t+x+y),2)))/2. // non-linearity
7214 * +std::cos(x)*std::cos(y+t); // pressure
7225 * -(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.
7226 * +3*std::pow(Pi,2)*std::cos(t)*std::cos(Pi*y)*std::cos(Pi*z)*std::sin(Pi*x) //viscosity
7227 * -(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
7228 * +std::cos(t+x+y+z); // pressure
7230 * -(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
7231 * +3*std::pow(Pi,2)*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*z)*std::sin(Pi*y) //viscosity
7232 * -(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
7233 * +std::cos(t+x+y+z); // pressure
7235 * 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
7236 * -6*std::pow(Pi,2)*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::sin(Pi*z) //viscosity
7237 * -(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
7238 * +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)
void refine_global(const unsigned int times=1)
virtual types::global_cell_index n_global_active_cells() const override
types::subdomain_id locally_owned_subdomain() const override
virtual void clear() override
__global__ void set(Number *val, const Number s, const size_type 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