246 *
#include <deal.II/base/quadrature_lib.h>
247 *
#include <deal.II/base/function.h>
248 *
#include <deal.II/lac/affine_constraints.h>
249 *
#include <deal.II/lac/vector.h>
250 *
#include <deal.II/lac/full_matrix.h>
251 *
#include <deal.II/lac/solver_cg.h>
252 *
#include <deal.II/lac/petsc_sparse_matrix.h>
253 *
#include <deal.II/lac/petsc_sparse_matrix.h>
254 *
#include <deal.II/lac/petsc_vector.h>
255 *
#include <deal.II/lac/petsc_solver.h>
256 *
#include <deal.II/lac/petsc_precondition.h>
257 *
#include <deal.II/grid/grid_generator.h>
258 *
#include <deal.II/grid/tria_accessor.h>
259 *
#include <deal.II/grid/tria_iterator.h>
260 *
#include <deal.II/dofs/dof_handler.h>
261 *
#include <deal.II/dofs/dof_accessor.h>
262 *
#include <deal.II/dofs/dof_tools.h>
263 *
#include <deal.II/fe/fe_values.h>
264 *
#include <deal.II/fe/fe_q.h>
265 *
#include <deal.II/numerics/vector_tools.h>
266 *
#include <deal.II/numerics/data_out.h>
267 *
#include <deal.II/numerics/error_estimator.h>
268 *
#include <deal.II/base/utilities.h>
269 *
#include <deal.II/base/conditional_ostream.h>
270 *
#include <deal.II/base/index_set.h>
271 *
#include <deal.II/lac/sparsity_tools.h>
272 *
#include <deal.II/distributed/tria.h>
273 *
#include <deal.II/distributed/grid_refinement.h>
274 *
#include <deal.II/lac/vector.h>
275 *
#include <deal.II/base/convergence_table.h>
276 *
#include <deal.II/base/timer.h>
277 *
#include <deal.II/base/parameter_handler.h>
278 *
#include <deal.II/grid/grid_tools.h>
279 *
#include <deal.II/fe/mapping_q.h>
284 *
#include <iostream>
295 *
#define CHECK_MAX_PRINCIPLE 0
299 * LOG FOR LEVEL SET FROM -1 to 1
303 *
#define ENTROPY_GRAD(phi,phix) 2*phi*phix*((1-phi*phi>=0) ? -1 : 1)/(
std::
abs(1-phi*phi)+1
E-14)
310 * This is a solver
for the transpor solver.
311 * We assume the velocity is divergence free
312 * and solve the equation in conservation form.
314 * ---------- NOTATION ----------
316 * We use notation popular in the literature of conservation laws.
317 * For this reason the solution is denoted as u, unm1, unp1, etc.
318 * and the velocity is treated as vx, vy and vz.
322 *
class LevelSetSolver
342 * BOUNDARY CONDITIONS
346 *
void set_boundary_conditions(std::vector<types::global_dof_index> &boundary_values_id_u,
347 *
std::vector<double> boundary_values_u);
375 *
void nth_time_step();
385 *
LevelSetSolver (
const unsigned int degree_LS,
386 *
const unsigned int degree_U,
387 *
const double time_step,
390 *
const bool verbose,
391 *
std::string ALGORITHM,
392 *
const unsigned int TIME_INTEGRATION,
401 * ASSEMBLE MASS (and other) MATRICES
405 *
void assemble_ML();
407 *
void assemble_MC();
411 * LOW ORDER METHOD (DiJ Viscosity)
415 *
void assemble_C_Matrix();
425 *
void assemble_EntRes_Matrix();
429 * FOR MAXIMUM PRINCIPLE
455 *
std::string algorithm);
458 *
std::string algorithm);
466 *
void get_sparsity_pattern();
467 *
void get_map_from_Q1_to_Q2();
470 *
std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
473 *
void save_old_solution();
474 *
void save_old_vel_solution();
483 *
const std::vector<types::global_dof_index> &indices,
484 *
std::vector<PetscScalar> &values);
486 *
const std::vector<types::global_dof_index> &indices,
487 *
std::map<types::global_dof_index, types::global_dof_index> &map_from_Q1_to_Q2,
488 *
std::vector<PetscScalar> &values);
494 * FINITE ELEMENT SPACE
502 *
IndexSet locally_relevant_dofs_LS;
512 * OPERATORS times SOLUTION VECTOR
525 *
std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> MC_preconditioner;
532 *
std::vector<types::global_dof_index> boundary_values_id_u;
533 *
std::vector<double> boundary_values_u;
540 * FOR FIRST ORDER VISCOSITY
547 * FOR ENTROPY VISCOSITY
553 * FOR FCT (flux and limited flux)
582 * NON-GHOSTED VECTORS
620 *
double solver_tolerance;
621 *
double entropy_normalization_factor;
629 *
std::string ALGORITHM;
630 *
unsigned int TIME_INTEGRATION;
634 *
std::map<types::global_dof_index, types::global_dof_index> map_from_Q1_to_Q2;
635 *
std::map<types::global_dof_index, std::vector<types::global_dof_index> > sparsity_pattern;
639 *
LevelSetSolver<dim>::LevelSetSolver (
const unsigned int degree_LS,
640 *
const unsigned int degree_U,
641 *
const double time_step,
644 *
const bool verbose,
645 *
std::string ALGORITHM,
646 *
const unsigned int TIME_INTEGRATION,
650 *
mpi_communicator (mpi_communicator),
651 *
degree_LS(degree_LS),
652 *
dof_handler_LS (triangulation),
654 *
degree_U(degree_U),
655 *
dof_handler_U (triangulation),
657 *
time_step(time_step),
661 *
ALGORITHM(ALGORITHM),
662 *
TIME_INTEGRATION(TIME_INTEGRATION),
665 *
pcout <<
"********** LEVEL SET SETUP **********" << std::endl;
670 *
LevelSetSolver<dim>::~LevelSetSolver ()
672 *
dof_handler_LS.clear ();
673 *
dof_handler_U.clear ();
692 *
this->locally_relevant_solution_vx = locally_relevant_solution_vx;
693 *
this->locally_relevant_solution_vy = locally_relevant_solution_vy;
696 * initialize old vectors with current solution,
this just happens the
first time
700 *
locally_relevant_solution_vx_old = locally_relevant_solution_vx;
701 *
locally_relevant_solution_vy_old = locally_relevant_solution_vy;
711 *
this->locally_relevant_solution_vx = locally_relevant_solution_vx;
712 *
this->locally_relevant_solution_vy = locally_relevant_solution_vy;
713 *
this->locally_relevant_solution_vz = locally_relevant_solution_vz;
716 * initialize old vectors with current solution,
this just happens the
first time
720 *
locally_relevant_solution_vx_old = locally_relevant_solution_vx;
721 *
locally_relevant_solution_vy_old = locally_relevant_solution_vy;
722 *
locally_relevant_solution_vz_old = locally_relevant_solution_vz;
733 *
void LevelSetSolver<dim>::set_boundary_conditions(std::vector<types::global_dof_index> &boundary_values_id_u,
734 *
std::vector<double> boundary_values_u)
736 *
this->boundary_values_id_u = boundary_values_id_u;
737 *
this->boundary_values_u = boundary_values_u;
756 *
save_old_vel_solution();
762 *
this->locally_relevant_solution_vx=locally_relevant_solution_vx;
763 *
this->locally_relevant_solution_vy=locally_relevant_solution_vy;
776 *
save_old_vel_solution();
782 *
this->locally_relevant_solution_vx=locally_relevant_solution_vx;
783 *
this->locally_relevant_solution_vy=locally_relevant_solution_vy;
784 *
this->locally_relevant_solution_vz=locally_relevant_solution_vz;
799 * -------------------------------------------------------------------------------
800 * ------------------------------ COMPUTE SOLUTIONS ------------------------------
801 * -------------------------------------------------------------------------------
805 *
void LevelSetSolver<dim>::nth_time_step()
807 *
assemble_EntRes_Matrix();
813 *
if (TIME_INTEGRATION==FORWARD_EULER)
814 *
compute_solution(unp1,un,ALGORITHM);
816 *
compute_solution_SSP33(unp1,un,ALGORITHM);
819 * BOUNDARY CONDITIONS
822 *
unp1.set(boundary_values_id_u,boundary_values_u);
826 * CHECK MAXIMUM PRINCIPLE
829 *
if (CHECK_MAX_PRINCIPLE)
831 *
compute_bounds(un);
832 *
check_max_principle(unp1);
836 * pcout <<
"*********************************************************************... "
837 * << unp1.min() <<
", " << unp1.max() << std::endl;
840 *
save_old_solution();
845 * --------------------------------------------------------------------
846 * ------------------------------ SETUP ------------------------------
847 * --------------------------------------------------------------------
851 *
void LevelSetSolver<dim>::setup()
853 *
solver_tolerance=1
E-6;
854 *
degree_MAX =
std::max(degree_LS,degree_U);
858 * SETUP FOR DOF HANDLERS
863 *
dof_handler_LS.distribute_dofs (fe_LS);
864 *
locally_owned_dofs_LS = dof_handler_LS.locally_owned_dofs ();
871 *
dof_handler_U.distribute_dofs (fe_U);
872 *
locally_owned_dofs_U = dof_handler_U.locally_owned_dofs ();
881 *
constraints.clear ();
882 *
constraints.reinit (locally_owned_dofs_LS, locally_relevant_dofs_LS);
884 *
constraints.close ();
888 * NON-GHOSTED VECTORS
892 *
MPP_uL_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
893 *
NMPP_uH_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
894 *
RHS.reinit(locally_owned_dofs_LS,mpi_communicator);
895 *
uStage1_nonGhosted.reinit (locally_owned_dofs_LS,mpi_communicator);
896 *
uStage2_nonGhosted.reinit (locally_owned_dofs_LS,mpi_communicator);
897 *
unp1.reinit (locally_owned_dofs_LS,mpi_communicator);
898 *
MPP_uH_solution.reinit (locally_owned_dofs_LS,mpi_communicator);
901 * vectors
for lumped mass
matrix
904 *
ML_vector.reinit(locally_owned_dofs_LS,mpi_communicator);
905 *
inverse_ML_vector.reinit(locally_owned_dofs_LS,mpi_communicator);
906 *
ones_vector.reinit(locally_owned_dofs_LS,mpi_communicator);
910 * operators times solution
913 *
K_times_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
914 *
DL_times_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
915 *
DH_times_solution.reinit(locally_owned_dofs_LS,mpi_communicator);
921 *
R_pos_vector_nonGhosted.reinit (locally_owned_dofs_LS,mpi_communicator);
922 *
R_neg_vector_nonGhosted.reinit (locally_owned_dofs_LS,mpi_communicator);
923 *
umin_vector.reinit (locally_owned_dofs_LS,mpi_communicator);
924 *
umax_vector.reinit (locally_owned_dofs_LS,mpi_communicator);
928 * GHOSTED VECTORS (used within some assemble process)
932 *
uStage1.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
933 *
uStage2.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
934 *
unm1.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
935 *
un.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
936 *
MPP_uL_solution_ghosted.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
937 *
MPP_uLkp1_solution_ghosted.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
938 *
NMPP_uH_solution_ghosted.reinit (locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
941 * init vectors
for vx
944 *
locally_relevant_solution_vx.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
945 *
locally_relevant_solution_vx_old.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
948 * init vectors
for vy
951 *
locally_relevant_solution_vy.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
952 *
locally_relevant_solution_vy_old.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
955 * init vectors
for vz
958 *
locally_relevant_solution_vz.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
959 *
locally_relevant_solution_vz_old.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
965 *
R_pos_vector.reinit(locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
966 *
R_neg_vector.reinit(locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
978 *
dof_handler_LS.locally_owned_dofs(),
980 *
locally_relevant_dofs_LS);
981 *
MC_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
982 *
dof_handler_LS.locally_owned_dofs(),
983 *
dsp, mpi_communicator);
984 *
Cx_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
985 *
dof_handler_LS.locally_owned_dofs(),
986 *
dsp, mpi_communicator);
987 *
CTx_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
988 *
dof_handler_LS.locally_owned_dofs(),
989 *
dsp, mpi_communicator);
990 *
Cy_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
991 *
dof_handler_LS.locally_owned_dofs(),
992 *
dsp, mpi_communicator);
993 *
CTy_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
994 *
dof_handler_LS.locally_owned_dofs(),
995 *
dsp, mpi_communicator);
998 *
Cz_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
999 *
dof_handler_LS.locally_owned_dofs(),
1000 *
dsp, mpi_communicator);
1001 *
CTz_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1002 *
dof_handler_LS.locally_owned_dofs(),
1003 *
dsp, mpi_communicator);
1005 *
dLij_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1006 *
dof_handler_LS.locally_owned_dofs(),
1007 *
dsp, mpi_communicator);
1008 *
EntRes_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1009 *
dof_handler_LS.locally_owned_dofs(),
1010 *
dsp, mpi_communicator);
1011 *
SuppSize_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1012 *
dof_handler_LS.locally_owned_dofs(),
1013 *
dsp, mpi_communicator);
1014 *
dCij_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1015 *
dof_handler_LS.locally_owned_dofs(),
1016 *
dsp, mpi_communicator);
1017 *
A_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1018 *
dof_handler_LS.locally_owned_dofs(),
1019 *
dsp, mpi_communicator);
1020 *
LxA_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1021 *
dof_handler_LS.locally_owned_dofs(),
1022 *
dsp, mpi_communicator);
1023 *
Akp1_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1024 *
dof_handler_LS.locally_owned_dofs(),
1025 *
dsp, mpi_communicator);
1026 *
LxAkp1_matrix.reinit (dof_handler_LS.locally_owned_dofs(),
1027 *
dof_handler_LS.locally_owned_dofs(),
1028 *
dsp, mpi_communicator);
1031 * COMPUTE MASS MATRICES (AND OTHERS) FOR FIRST TIME STEP
1037 *
assemble_C_Matrix();
1040 * get mat
for DOFs between Q1 and Q2
1043 *
get_map_from_Q1_to_Q2();
1044 *
get_sparsity_pattern();
1049 * ----------------------------------------------------------------------------
1050 * ------------------------------ MASS MATRICES ------------------------------
1051 * ----------------------------------------------------------------------------
1055 *
void LevelSetSolver<dim>::assemble_ML()
1059 *
const QGauss<dim> quadrature_formula(degree_MAX+1);
1065 *
const unsigned int dofs_per_cell = fe_LS.dofs_per_cell;
1066 *
const unsigned int n_q_points = quadrature_formula.size();
1069 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1072 *
cell_LS = dof_handler_LS.begin_active(),
1073 *
endc_LS = dof_handler_LS.end();
1075 *
for (; cell_LS!=endc_LS; ++cell_LS)
1076 *
if (cell_LS->is_locally_owned())
1079 *
fe_values_LS.reinit (cell_LS);
1080 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
1082 *
const double JxW = fe_values_LS.JxW(q_point);
1083 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1084 *
cell_ML (i) += fe_values_LS.shape_value(i,q_point)*JxW;
1091 *
cell_LS->get_dof_indices (local_dof_indices);
1092 *
constraints.distribute_local_to_global (cell_ML,local_dof_indices,ML_vector);
1103 *
void LevelSetSolver<dim>::invert_ML()
1107 *
loop on locally owned i-DOFs (rows)
1111 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
1113 *
int gi = *idofs_iter;
1114 *
inverse_ML_vector(gi) = 1./ML_vector(gi);
1120 *
void LevelSetSolver<dim>::assemble_MC()
1124 *
const QGauss<dim> quadrature_formula(degree_MAX+1);
1130 *
const unsigned int dofs_per_cell = fe_LS.dofs_per_cell;
1131 *
const unsigned int n_q_points = quadrature_formula.size();
1134 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1135 *
std::vector<double> shape_values(dofs_per_cell);
1138 *
cell_LS = dof_handler_LS.begin_active(),
1139 *
endc_LS = dof_handler_LS.end();
1141 *
for (; cell_LS!=endc_LS; ++cell_LS)
1142 *
if (cell_LS->is_locally_owned())
1145 *
fe_values_LS.reinit (cell_LS);
1146 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
1148 *
const double JxW = fe_values_LS.JxW(q_point);
1149 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1150 *
shape_values[i] = fe_values_LS.shape_value(i,q_point);
1152 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1153 *
for (
unsigned int j=0; j<dofs_per_cell; ++j)
1154 *
cell_MC(i,j) += shape_values[i]*shape_values[j]*JxW;
1161 *
cell_LS->get_dof_indices (local_dof_indices);
1162 *
constraints.distribute_local_to_global (cell_MC,local_dof_indices,MC_matrix);
1175 * ---------------------------------------------------------------------------------------
1176 * ------------------------------
LO METHOD (Dij Viscosity) ------------------------------
1177 * ---------------------------------------------------------------------------------------
1180 *
template <
int dim>
1181 *
void LevelSetSolver<dim>::assemble_C_Matrix ()
1190 *
const QGauss<dim> quadrature_formula(degree_MAX+1);
1196 *
const unsigned int dofs_per_cell_LS = fe_LS.dofs_per_cell;
1197 *
const unsigned int n_q_points = quadrature_formula.size();
1206 *
std::vector<Tensor<1, dim> > shape_grads_LS(dofs_per_cell_LS);
1207 *
std::vector<double> shape_values_LS(dofs_per_cell_LS);
1209 *
std::vector<types::global_dof_index> local_dof_indices_LS (dofs_per_cell_LS);
1212 *
cell_LS = dof_handler_LS.begin_active();
1213 *
endc_LS = dof_handler_LS.end();
1215 *
for (; cell_LS!=endc_LS; ++cell_LS)
1216 *
if (cell_LS->is_locally_owned())
1228 *
fe_values_LS.reinit (cell_LS);
1229 *
cell_LS->get_dof_indices (local_dof_indices_LS);
1231 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
1233 *
const double JxW = fe_values_LS.JxW(q_point);
1234 *
for (
unsigned int i=0; i<dofs_per_cell_LS; ++i)
1236 *
shape_values_LS[i] = fe_values_LS.shape_value(i,q_point);
1237 *
shape_grads_LS [i] = fe_values_LS.shape_grad (i,q_point);
1240 *
for (
unsigned int i=0; i<dofs_per_cell_LS; ++i)
1241 *
for (
unsigned int j=0; j < dofs_per_cell_LS; ++j)
1243 *
cell_Cij_x(i,j) += (shape_grads_LS[j][0])*shape_values_LS[i]*JxW;
1244 *
cell_Cij_y(i,j) += (shape_grads_LS[j][1])*shape_values_LS[i]*JxW;
1245 *
cell_Cji_x(i,j) += (shape_grads_LS[i][0])*shape_values_LS[j]*JxW;
1246 *
cell_Cji_y(i,j) += (shape_grads_LS[i][1])*shape_values_LS[j]*JxW;
1249 *
cell_Cij_z(i,j) += (shape_grads_LS[j][2])*shape_values_LS[i]*JxW;
1250 *
cell_Cji_z(i,j) += (shape_grads_LS[i][2])*shape_values_LS[j]*JxW;
1259 *
constraints.distribute_local_to_global(cell_Cij_x,local_dof_indices_LS,Cx_matrix);
1260 *
constraints.distribute_local_to_global(cell_Cji_x,local_dof_indices_LS,CTx_matrix);
1261 *
constraints.distribute_local_to_global(cell_Cij_y,local_dof_indices_LS,Cy_matrix);
1262 *
constraints.distribute_local_to_global(cell_Cji_y,local_dof_indices_LS,CTy_matrix);
1265 *
constraints.distribute_local_to_global(cell_Cij_z,local_dof_indices_LS,Cz_matrix);
1266 *
constraints.distribute_local_to_global(cell_Cji_z,local_dof_indices_LS,CTz_matrix);
1288 *
K_times_solution = 0;
1290 *
const QGauss<dim> quadrature_formula(degree_MAX+1);
1299 *
const unsigned int dofs_per_cell = fe_LS.dofs_per_cell;
1300 *
const unsigned int n_q_points = quadrature_formula.size();
1304 *
std::vector<Tensor<1,dim> > un_grads (n_q_points);
1305 *
std::vector<double> old_vx_values (n_q_points);
1306 *
std::vector<double> old_vy_values (n_q_points);
1307 *
std::vector<double> old_vz_values (n_q_points);
1309 *
std::vector<double> shape_values(dofs_per_cell);
1310 *
std::vector<Tensor<1,dim> > shape_grads(dofs_per_cell);
1314 *
std::vector<types::global_dof_index> indices_LS (dofs_per_cell);
1322 *
cell_LS = dof_handler_LS.begin_active(),
1323 *
endc_LS = dof_handler_LS.end();
1325 *
cell_U = dof_handler_U.begin_active();
1328 *
for (; cell_LS!=endc_LS; ++cell_U, ++cell_LS)
1329 *
if (cell_LS->is_locally_owned())
1331 *
cell_K_times_solution=0;
1333 *
fe_values_LS.reinit (cell_LS);
1334 *
cell_LS->get_dof_indices (indices_LS);
1335 *
fe_values_LS.get_function_gradients(solution,un_grads);
1337 *
fe_values_U.reinit (cell_U);
1338 *
fe_values_U.get_function_values(locally_relevant_solution_vx,old_vx_values);
1339 *
fe_values_U.get_function_values(locally_relevant_solution_vy,old_vy_values);
1340 *
if (dim==3) fe_values_U.get_function_values(locally_relevant_solution_vz,old_vz_values);
1344 * compute cell_K_times_solution
1347 *
for (
unsigned int q_point=0; q_point<n_q_points; ++q_point)
1349 *
v[0] = old_vx_values[q_point];
1350 *
v[1] = old_vy_values[q_point];
1351 *
if (dim==3) v[2] = old_vz_values[q_point];
1353 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1354 *
cell_K_times_solution(i) += (v*un_grads[q_point])
1355 *
*fe_values_LS.shape_value(i,q_point)*fe_values_LS.JxW(q_point);
1362 *
constraints.distribute_local_to_global (cell_K_times_solution, indices_LS, K_times_solution);
1367 *
template <
int dim>
1368 *
void LevelSetSolver<dim>::assemble_K_DL_DH_times_vector
1373 * K_times_solution=0;
1376 *
DL_times_solution=0;
1377 *
DH_times_solution=0;
1381 *
PetscInt ncolumns;
1382 *
const PetscInt *gj;
1383 *
const PetscScalar *Cxi, *Cyi, *Czi, *CTxi, *CTyi, *CTzi;
1384 *
const PetscScalar *EntResi, *SuppSizei, *MCi;
1391 *
loop on locally owned i-DOFs (rows)
1396 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
1398 *
PetscInt gi = *idofs_iter;
1401 *
double ith_K_times_solution = 0;
1405 * read velocity of i-th DOF
1408 *
vi[0] = locally_relevant_solution_vx(map_from_Q1_to_Q2[gi]);
1409 *
vi[1] = locally_relevant_solution_vy(map_from_Q1_to_Q2[gi]);
1410 *
if (dim==3) vi[2] = locally_relevant_solution_vz(map_from_Q1_to_Q2[gi]);
1411 *
solni = solution(gi);
1415 * get i-th row of
C matrices
1418 *
MatGetRow(Cx_matrix,gi,&ncolumns,&gj,&Cxi);
1419 *
MatGetRow(Cy_matrix,gi,&ncolumns,&gj,&Cyi);
1420 *
MatGetRow(CTx_matrix,gi,&ncolumns,&gj,&CTxi);
1421 *
MatGetRow(CTy_matrix,gi,&ncolumns,&gj,&CTyi);
1424 *
MatGetRow(Cz_matrix,gi,&ncolumns,&gj,&Czi);
1425 *
MatGetRow(CTz_matrix,gi,&ncolumns,&gj,&CTzi);
1427 *
MatGetRow(EntRes_matrix,gi,&ncolumns,&gj,&EntResi);
1428 *
MatGetRow(SuppSize_matrix,gi,&ncolumns,&gj,&SuppSizei);
1429 *
MatGetRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
1433 * get vector
values for column indices
1436 *
const std::vector<types::global_dof_index> gj_indices (gj,gj+ncolumns);
1437 *
std::vector<double> soln(ncolumns);
1438 *
std::vector<double> vx(ncolumns);
1439 *
std::vector<double> vy(ncolumns);
1440 *
std::vector<double> vz(ncolumns);
1441 *
get_vector_values(solution,gj_indices,soln);
1442 *
get_vector_values(locally_relevant_solution_vx,gj_indices,map_from_Q1_to_Q2,vx);
1443 *
get_vector_values(locally_relevant_solution_vy,gj_indices,map_from_Q1_to_Q2,vy);
1445 *
get_vector_values(locally_relevant_solution_vz,gj_indices,map_from_Q1_to_Q2,vz);
1449 * Array
for i-th row of matrices
1452 *
std::vector<double> dLi(ncolumns), dCi(ncolumns);
1453 *
double dLii = 0, dCii = 0;
1456 *
loop on sparsity pattern of i-th DOF
1459 *
for (
int j =0; j < ncolumns; ++j)
1476 * ith_K_times_solution += soln[j]*(vj*
C);
1483 * low order dissipative
matrix
1490 * high order dissipative
matrix (entropy viscosity)
1494 *
cE*
std::abs(EntResi[j])/(entropy_normalization_factor*MCi[j]/SuppSizei[j]));
1497 * high order compression
matrix
1501 *
dCi[j] = dEij*
std::max(1-Compij,0.0);
1507 *
save K times solution vector
1508 * K_times_solution(gi)=ith_K_times_solution;
1509 *
save i-th row of matrices on global matrices
1512 *
MatSetValuesRow(dLij_matrix,gi,&dLi[0]);
1513 *
dLij_matrix.set(gi,gi,dLii);
1514 *
MatSetValuesRow(dCij_matrix,gi,&dCi[0]);
1515 *
dCij_matrix.set(gi,gi,dCii);
1519 * Restore matrices after reading rows
1522 *
MatRestoreRow(Cx_matrix,gi,&ncolumns,&gj,&Cxi);
1523 *
MatRestoreRow(Cy_matrix,gi,&ncolumns,&gj,&Cyi);
1524 *
MatRestoreRow(CTx_matrix,gi,&ncolumns,&gj,&CTxi);
1525 *
MatRestoreRow(CTy_matrix,gi,&ncolumns,&gj,&CTyi);
1528 *
MatRestoreRow(Cz_matrix,gi,&ncolumns,&gj,&Czi);
1529 *
MatRestoreRow(CTz_matrix,gi,&ncolumns,&gj,&CTzi);
1531 *
MatRestoreRow(EntRes_matrix,gi,&ncolumns,&gj,&EntResi);
1532 *
MatRestoreRow(SuppSize_matrix,gi,&ncolumns,&gj,&SuppSizei);
1533 *
MatRestoreRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
1545 * get matrices times vector
1548 *
dLij_matrix.vmult(DL_times_solution,solution);
1549 *
dCij_matrix.vmult(DH_times_solution,solution);
1554 * --------------------------------------------------------------------------------------
1555 * ------------------------------ ENTROPY VISCOSITY ------------------------------
1556 * --------------------------------------------------------------------------------------
1559 *
template <
int dim>
1560 *
void LevelSetSolver<dim>::assemble_EntRes_Matrix ()
1563 *
entropy_normalization_factor=0;
1564 *
SuppSize_matrix=0;
1566 *
const QGauss<dim> quadrature_formula(degree_MAX+1);
1576 *
const unsigned int dofs_per_cell_LS = fe_LS.dofs_per_cell;
1577 *
const unsigned int n_q_points = quadrature_formula.size();
1579 *
std::vector<double> uqn (n_q_points);
1580 *
std::vector<double> uqnm1 (n_q_points);
1581 *
std::vector<Tensor<1,dim> > guqn (n_q_points);
1582 *
std::vector<Tensor<1,dim> > guqnm1 (n_q_points);
1584 *
std::vector<double> vxqn (n_q_points);
1585 *
std::vector<double> vyqn (n_q_points);
1586 *
std::vector<double> vzqn (n_q_points);
1587 *
std::vector<double> vxqnm1 (n_q_points);
1588 *
std::vector<double> vyqnm1 (n_q_points);
1589 *
std::vector<double> vzqnm1 (n_q_points);
1594 *
std::vector<Tensor<1, dim> > shape_grads_LS(dofs_per_cell_LS);
1595 *
std::vector<double> shape_values_LS(dofs_per_cell_LS);
1597 *
std::vector<types::global_dof_index> local_dof_indices_LS (dofs_per_cell_LS);
1600 *
cell_LS = dof_handler_LS.begin_active();
1601 *
endc_LS = dof_handler_LS.end();
1603 *
cell_U = dof_handler_U.begin_active();
1606 *
double max_entropy=-1E10, min_entropy=1E10;
1607 *
double cell_max_entropy, cell_min_entropy;
1608 *
double cell_entropy_mass, entropy_mass=0;
1609 *
double cell_volume_double,
volume=0;
1611 *
for (; cell_LS!=endc_LS; ++cell_LS, ++cell_U)
1612 *
if (cell_LS->is_locally_owned())
1614 *
cell_entropy_mass = 0;
1615 *
cell_volume_double = 0;
1616 *
cell_max_entropy = -1E10;
1617 *
cell_min_entropy = 1E10;
1623 * get solutions at quadrature points
1626 *
fe_values_LS.reinit(cell_LS);
1627 *
cell_LS->get_dof_indices (local_dof_indices_LS);
1628 *
fe_values_LS.get_function_values(un,uqn);
1629 *
fe_values_LS.get_function_values(unm1,uqnm1);
1630 *
fe_values_LS.get_function_gradients(un,guqn);
1631 *
fe_values_LS.get_function_gradients(unm1,guqnm1);
1633 *
fe_values_U.reinit(cell_U);
1634 *
fe_values_U.get_function_values(locally_relevant_solution_vx,vxqn);
1635 *
fe_values_U.get_function_values(locally_relevant_solution_vy,vyqn);
1636 *
if (dim==3) fe_values_U.get_function_values(locally_relevant_solution_vz,vzqn);
1637 *
fe_values_U.get_function_values(locally_relevant_solution_vx_old,vxqnm1);
1638 *
fe_values_U.get_function_values(locally_relevant_solution_vy_old,vyqnm1);
1639 *
if (dim==3) fe_values_U.get_function_values(locally_relevant_solution_vz_old,vzqnm1);
1641 *
for (
unsigned int q=0; q<n_q_points; ++q)
1643 *
Rk = 1./time_step*(ENTROPY(uqn[q])-ENTROPY(uqnm1[q]))
1644 *
+(vxqn[q]*ENTROPY_GRAD(uqn[q],guqn[q][0])+vyqn[q]*ENTROPY_GRAD(uqn[q],guqn[q][1]))/2.
1645 *
+(vxqnm1[q]*ENTROPY_GRAD(uqnm1[q],guqnm1[q][0])+vyqnm1[q]*ENTROPY_GRAD(uqnm1[q],guqnm1[q][1]))/2.;
1647 *
Rk += 0.5*(vzqn[q]*ENTROPY_GRAD(uqn[q],guqn[q][2])+vzqnm1[q]*ENTROPY_GRAD(uqnm1[q],guqnm1[q][2]));
1649 *
const double JxW = fe_values_LS.JxW(q);
1650 *
for (
unsigned int i=0; i<dofs_per_cell_LS; ++i)
1652 *
shape_values_LS[i] = fe_values_LS.shape_value(i,q);
1653 *
shape_grads_LS [i] = fe_values_LS.shape_grad (i,q);
1656 *
for (
unsigned int i=0; i<dofs_per_cell_LS; ++i)
1657 *
for (
unsigned int j=0; j < dofs_per_cell_LS; ++j)
1659 *
cell_EntRes (i,j) += Rk*shape_values_LS[i]*shape_values_LS[j]*JxW;
1660 *
cell_volume (i,j) += JxW;
1662 *
cell_entropy_mass += ENTROPY(uqn[q])*JxW;
1663 *
cell_volume_double += JxW;
1665 *
cell_min_entropy =
std::min(cell_min_entropy,ENTROPY(uqn[q]));
1666 *
cell_max_entropy =
std::max(cell_max_entropy,ENTROPY(uqn[q]));
1668 *
entropy_mass += cell_entropy_mass;
1669 *
volume += cell_volume_double;
1671 *
min_entropy =
std::min(min_entropy,cell_min_entropy);
1672 *
max_entropy =
std::max(max_entropy,cell_max_entropy);
1678 *
constraints.distribute_local_to_global(cell_EntRes,local_dof_indices_LS,EntRes_matrix);
1679 *
constraints.distribute_local_to_global(cell_volume,local_dof_indices_LS,SuppSize_matrix);
1685 * ENTROPY NORM FACTOR
1697 * ------------------------------------------------------------------------------------
1698 * ------------------------------ TO CHECK MAX PRINCIPLE ------------------------------
1699 * ------------------------------------------------------------------------------------
1709 *
loop on locally owned i-DOFs (rows)
1713 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
1715 *
int gi = *idofs_iter;
1719 * get solution at DOFs on the sparsity pattern of i-th DOF
1722 *
std::vector<types::global_dof_index> gj_indices = sparsity_pattern[gi];
1723 *
std::vector<double> soln(gj_indices.size());
1724 *
get_vector_values(un_solution,gj_indices,soln);
1727 * compute bounds, ith row of flux
matrix, P vectors
1730 *
double mini=1E10, maxi=-1E10;
1731 *
for (
unsigned int j =0; j < gj_indices.size(); ++j)
1741 *
umin_vector(gi) = mini;
1742 *
umax_vector(gi) = maxi;
1753 * compute
min and
max vectors
1756 *
const unsigned int dofs_per_cell = fe_LS.dofs_per_cell;
1757 *
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1761 *
cell_LS = dof_handler_LS.begin_active(),
1762 *
endc_LS = dof_handler_LS.end();
1764 *
for (; cell_LS!=endc_LS; ++cell_LS)
1765 *
if (cell_LS->is_locally_owned() && !cell_LS->at_boundary())
1767 *
cell_LS->get_dof_indices(local_dof_indices);
1768 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
1769 *
if (locally_owned_dofs_LS.is_element(local_dof_indices[i]))
1771 *
double solni = unp1_solution(local_dof_indices[i]);
1772 *
if (solni - umin_vector(local_dof_indices[i]) < -tol || umax_vector(local_dof_indices[i]) - solni < -tol)
1774 *
pcout <<
"MAX Principle violated" << std::endl;
1783 * -------------------------------------------------------------------------------
1784 * ------------------------------ COMPUTE SOLUTIONS ------------------------------
1785 * -------------------------------------------------------------------------------
1789 *
void LevelSetSolver<dim>::compute_MPP_uL_and_NMPP_uH
1796 * NON-GHOSTED VECTORS: MPP_uL_solution, NMPP_uH_solution
1797 * GHOSTED VECTORS: un_solution
1800 *
MPP_uL_solution=un_solution;
1801 *
NMPP_uH_solution=un_solution;
1807 *
assemble_K_times_vector(un_solution);
1808 *
assemble_K_DL_DH_times_vector(un_solution);
1812 * COMPUTE MPP u1 solution
1816 *
MPP_uL_solution.scale(ML_vector);
1817 *
MPP_uL_solution.add(-time_step,K_times_solution);
1818 *
MPP_uL_solution.add(-time_step,DL_times_solution);
1819 *
MPP_uL_solution.scale(inverse_ML_vector);
1823 * COMPUTE GALERKIN u2 solution
1827 *
MC_matrix.vmult(RHS,un_solution);
1828 *
RHS.add(-time_step,K_times_solution,-time_step,DH_times_solution);
1829 *
solve(constraints,MC_matrix,MC_preconditioner,NMPP_uH_solution,RHS);
1832 *
template <
int dim>
1833 *
void LevelSetSolver<dim>::compute_MPP_uH
1839 *
MPP_uH_solution=0;
1842 *
loop on locally owned i-DOFs (rows)
1847 *
PetscInt ncolumns;
1848 *
const PetscInt *gj;
1849 *
const PetscScalar *MCi, *dLi, *dCi;
1850 *
double solni, mi, solLi, solHi;
1852 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
1854 *
int gi = *idofs_iter;
1857 * read vectors at i-th DOF
1860 *
solni=solution(gi);
1861 *
solHi=NMPP_uH_solution_ghosted(gi);
1862 *
solLi=MPP_uL_solution_ghosted(gi);
1867 * get i-th row of matrices
1870 *
MatGetRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
1871 *
MatGetRow(dLij_matrix,gi,&ncolumns,&gj,&dLi);
1872 *
MatGetRow(dCij_matrix,gi,&ncolumns,&gj,&dCi);
1876 * get vector
values for support of i-th DOF
1879 *
const std::vector<types::global_dof_index> gj_indices (gj,gj+ncolumns);
1880 *
std::vector<double> soln(ncolumns);
1881 *
std::vector<double> solH(ncolumns);
1882 *
get_vector_values(solution,gj_indices,soln);
1883 *
get_vector_values(NMPP_uH_solution_ghosted,gj_indices,solH);
1887 * Array
for i-th row of matrices
1890 *
std::vector<double> Ai(ncolumns);
1893 * compute bounds, ith row of flux
matrix, P vectors
1896 *
double mini=1E10, maxi=-1E10;
1897 *
double Pposi=0 ,Pnegi=0;
1898 *
for (
int j =0; j < ncolumns; ++j)
1913 *
Ai[j] = (((gi==gj[j]) ? 1 : 0)*mi - MCi[j])*(solH[j]-soln[j] - (solHi-solni))
1914 *
+time_step*(dLi[j]-dCi[j])*(soln[j]-solni);
1921 *
Pposi += Ai[j]*((Ai[j] > 0) ? 1. : 0.);
1922 *
Pnegi += Ai[j]*((Ai[j] < 0) ? 1. : 0.);
1929 *
MatSetValuesRow(A_matrix,gi,&Ai[0]);
1936 *
double Qposi = mi*(maxi-solLi);
1937 *
double Qnegi = mi*(mini-solLi);
1944 *
R_pos_vector_nonGhosted(gi) = ((Pposi==0) ? 1. :
std::
min(1.0,Qposi/Pposi));
1945 *
R_neg_vector_nonGhosted(gi) = ((Pnegi==0) ? 1. :
std::
min(1.0,Qnegi/Pnegi));
1949 * Restore matrices after reading rows
1952 *
MatRestoreRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
1953 *
MatRestoreRow(dLij_matrix,gi,&ncolumns,&gj,&dLi);
1954 *
MatRestoreRow(dCij_matrix,gi,&ncolumns,&gj,&dCi);
1971 * update ghost
values for R vectors
1974 *
R_pos_vector = R_pos_vector_nonGhosted;
1975 *
R_neg_vector = R_neg_vector_nonGhosted;
1979 * compute limiters. NOTE:
this is a different
loop due to need of i- and j-th entries of R vectors
1983 *
double Rposi, Rnegi;
1984 *
idofs_iter=locally_owned_dofs_LS.begin();
1985 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
1987 *
int gi = *idofs_iter;
1988 *
Rposi = R_pos_vector(gi);
1989 *
Rnegi = R_neg_vector(gi);
1996 *
MatGetRow(A_matrix,gi,&ncolumns,&gj,&Ai);
2000 * get vector
values for column indices
2003 *
const std::vector<types::global_dof_index> gj_indices (gj,gj+ncolumns);
2004 *
std::vector<double> Rpos(ncolumns);
2005 *
std::vector<double> Rneg(ncolumns);
2006 *
get_vector_values(R_pos_vector,gj_indices,Rpos);
2007 *
get_vector_values(R_neg_vector,gj_indices,Rneg);
2011 * Array
for i-th row of A_times_L
matrix
2014 *
std::vector<double> LxAi(ncolumns);
2017 *
loop in sparsity pattern of i-th DOF
2020 *
for (
int j =0; j < ncolumns; ++j)
2021 *
LxAi[j] = Ai[j] * ((Ai[j]>0) ?
std::min(Rposi,Rneg[j]) :
std::
min(Rnegi,Rpos[j]));
2025 *
save i-th row of LxA
2028 *
MatSetValuesRow(LxA_matrix,gi,&LxAi[0]);
2031 * restore
A matrix after reading it
2034 *
MatRestoreRow(A_matrix,gi,&ncolumns,&gj,&Ai);
2037 *
LxA_matrix.vmult(MPP_uH_solution,ones_vector);
2038 *
MPP_uH_solution.scale(inverse_ML_vector);
2039 *
MPP_uH_solution.add(1.0,MPP_uL_solution_ghosted);
2043 *
void LevelSetSolver<dim>::compute_MPP_uH_with_iterated_FCT
2049 *
MPP_uH_solution=0;
2050 *
compute_MPP_uH(MPP_uH_solution,MPP_uL_solution_ghosted,NMPP_uH_solution_ghosted,un_solution);
2054 *
Akp1_matrix.copy_from(A_matrix);
2055 *
LxAkp1_matrix.copy_from(LxA_matrix);
2059 *
loop in num of FCT iterations
2062 *
PetscInt ncolumns;
2063 *
const PetscInt *gj;
2064 *
const PetscScalar *Akp1i;
2066 *
for (
int iter=0; iter<NUM_ITER; ++iter)
2068 *
MPP_uLkp1_solution_ghosted = MPP_uH_solution;
2069 *
Akp1_matrix.add(-1.0, LxAkp1_matrix);
2073 *
loop on locally owned i-DOFs (rows)
2077 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
2079 *
int gi = *idofs_iter;
2083 * read vectors at i-th DOF
2087 *
double solLi = MPP_uLkp1_solution_ghosted(gi);
2091 * get i-th row of matrices
2094 *
MatGetRow(Akp1_matrix,gi,&ncolumns,&gj,&Akp1i);
2097 * get vector
values for support of i-th DOF
2100 *
const std::vector<types::global_dof_index> gj_indices (gj,gj+ncolumns);
2101 *
std::vector<double> soln(ncolumns);
2102 *
get_vector_values(un_solution,gj_indices,soln);
2106 * compute bounds, ith row of flux
matrix, P vectors
2109 *
double mini=1E10, maxi=-1E10;
2110 *
double Pposi=0 ,Pnegi=0;
2111 *
for (
int j =0; j < ncolumns; ++j)
2126 *
Pposi += Akp1i[j]*((Akp1i[j] > 0) ? 1. : 0.);
2127 *
Pnegi += Akp1i[j]*((Akp1i[j] < 0) ? 1. : 0.);
2134 *
double Qposi = mi*(maxi-solLi);
2135 *
double Qnegi = mi*(mini-solLi);
2142 *
R_pos_vector_nonGhosted(gi) = ((Pposi==0) ? 1. :
std::
min(1.0,Qposi/Pposi));
2143 *
R_neg_vector_nonGhosted(gi) = ((Pnegi==0) ? 1. :
std::
min(1.0,Qnegi/Pnegi));
2147 * Restore matrices after reading rows
2150 *
MatRestoreRow(Akp1_matrix,gi,&ncolumns,&gj,&Akp1i);
2161 * update ghost
values for R vectors
2164 *
R_pos_vector = R_pos_vector_nonGhosted;
2165 *
R_neg_vector = R_neg_vector_nonGhosted;
2169 * compute limiters. NOTE:
this is a different
loop due to need of i- and j-th entries of R vectors
2172 *
double Rposi, Rnegi;
2173 *
idofs_iter=locally_owned_dofs_LS.begin();
2174 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
2176 *
int gi = *idofs_iter;
2177 *
Rposi = R_pos_vector(gi);
2178 *
Rnegi = R_neg_vector(gi);
2182 * get i-th row of Akp1
matrix
2185 *
MatGetRow(Akp1_matrix,gi,&ncolumns,&gj,&Akp1i);
2189 * get vector
values for column indices
2192 *
const std::vector<types::global_dof_index> gj_indices(gj,gj+ncolumns);
2193 *
std::vector<double> Rpos(ncolumns);
2194 *
std::vector<double> Rneg(ncolumns);
2195 *
get_vector_values(R_pos_vector,gj_indices,Rpos);
2196 *
get_vector_values(R_neg_vector,gj_indices,Rneg);
2200 * Array
for i-th row of LxAkp1
matrix
2203 *
std::vector<double> LxAkp1i(ncolumns);
2204 *
for (
int j =0; j < ncolumns; ++j)
2205 *
LxAkp1i[j] = Akp1i[j] * ((Akp1i[j]>0) ?
std::min(Rposi,Rneg[j]) :
std::
min(Rnegi,Rpos[j]));
2209 *
save i-th row of LxA
2212 *
MatSetValuesRow(LxAkp1_matrix,gi,&LxAkp1i[0]);
2215 * restore
A matrix after reading it
2218 *
MatRestoreRow(Akp1_matrix,gi,&ncolumns,&gj,&Akp1i);
2221 *
LxAkp1_matrix.vmult(MPP_uH_solution,ones_vector);
2222 *
MPP_uH_solution.scale(inverse_ML_vector);
2223 *
MPP_uH_solution.add(1.0,MPP_uLkp1_solution_ghosted);
2231 *
std::string algorithm)
2236 * COMPUTE MPP LOW-ORDER SOLN and NMPP HIGH-ORDER SOLN
2239 *
compute_MPP_uL_and_NMPP_uH(MPP_uL_solution,NMPP_uH_solution,un);
2241 *
if (algorithm.compare(
"MPP_u1")==0)
2242 *
unp1=MPP_uL_solution;
2243 *
else if (algorithm.compare(
"NMPP_uH")==0)
2244 *
unp1=NMPP_uH_solution;
2245 *
else if (algorithm.compare(
"MPP_uH")==0)
2247 *
MPP_uL_solution_ghosted = MPP_uL_solution;
2248 *
NMPP_uH_solution_ghosted=NMPP_uH_solution;
2249 *
compute_MPP_uH_with_iterated_FCT(MPP_uH_solution,MPP_uL_solution_ghosted,NMPP_uH_solution_ghosted,un);
2250 *
unp1=MPP_uH_solution;
2254 *
pcout <<
"Error in algorithm" << std::endl;
2262 *
std::string algorithm)
2266 * GHOSTED VECTORS: un
2267 * NON-GHOSTED VECTORS: unp1
2271 *
uStage1=0., uStage2=0.;
2272 *
uStage1_nonGhosted=0., uStage2_nonGhosted=0.;
2281 *
compute_solution(uStage1_nonGhosted,un,algorithm);
2282 *
uStage1=uStage1_nonGhosted;
2288 * u2=3/4*un+1/4*(u1-dt*RH*u1)
2291 *
compute_solution(uStage2_nonGhosted,uStage1,algorithm);
2292 *
uStage2_nonGhosted*=1./4;
2293 *
uStage2_nonGhosted.add(3./4,un);
2294 *
uStage2=uStage2_nonGhosted;
2300 * unp1=1/3*un+2/3*(u2-dt*RH*u2)
2303 *
compute_solution(unp1,uStage2,algorithm);
2305 *
unp1.add(1./3,un);
2310 * -----------------------------------------------------------------------
2311 * ------------------------------ UTILITIES ------------------------------
2312 * -----------------------------------------------------------------------
2316 *
void LevelSetSolver<dim>::get_sparsity_pattern()
2324 *
PetscInt ncolumns;
2325 *
const PetscInt *gj;
2326 *
const PetscScalar *MCi;
2328 *
for (; idofs_iter!=locally_owned_dofs_LS.end(); ++idofs_iter)
2330 *
PetscInt gi = *idofs_iter;
2333 * get i-th row of mass
matrix (dummy, I just need the indices gj)
2336 *
MatGetRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
2337 *
sparsity_pattern[gi] = std::vector<types::global_dof_index>(gj,gj+ncolumns);
2338 *
MatRestoreRow(MC_matrix,gi,&ncolumns,&gj,&MCi);
2343 *
void LevelSetSolver<dim>::get_map_from_Q1_to_Q2()
2345 *
map_from_Q1_to_Q2.clear();
2346 *
const unsigned int dofs_per_cell_LS = fe_LS.dofs_per_cell;
2347 *
std::vector<types::global_dof_index> local_dof_indices_LS (dofs_per_cell_LS);
2348 *
const unsigned int dofs_per_cell_U = fe_U.dofs_per_cell;
2349 *
std::vector<types::global_dof_index> local_dof_indices_U (dofs_per_cell_U);
2352 *
cell_LS = dof_handler_LS.begin_active(),
2353 *
endc_LS = dof_handler_LS.end();
2355 *
cell_U = dof_handler_U.begin_active();
2357 *
for (; cell_LS!=endc_LS; ++cell_LS, ++cell_U)
2358 *
if (!cell_LS->is_artificial())
2360 *
cell_LS->get_dof_indices(local_dof_indices_LS);
2361 *
cell_U->get_dof_indices(local_dof_indices_U);
2362 *
for (
unsigned int i=0; i<dofs_per_cell_LS; ++i)
2363 *
map_from_Q1_to_Q2[local_dof_indices_LS[i]] = local_dof_indices_U[i];
2367 *
template <
int dim>
2370 *
std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
2376 * all vectors are NON-GHOSTED
2379 *
SolverControl solver_control (dof_handler_LS.n_dofs(), solver_tolerance);
2381 *
constraints.distribute (completely_distributed_solution);
2382 *
solver.solve (Matrix, completely_distributed_solution, rhs, *preconditioner);
2383 *
constraints.distribute (completely_distributed_solution);
2384 *
if (verbose==
true) pcout <<
" Solved in " << solver_control.last_step() <<
" iterations." << std::endl;
2387 *
template <
int dim>
2388 *
void LevelSetSolver<dim>::save_old_solution()
2394 *
template <
int dim>
2395 *
void LevelSetSolver<dim>::save_old_vel_solution()
2397 *
locally_relevant_solution_vx_old = locally_relevant_solution_vx;
2398 *
locally_relevant_solution_vy_old = locally_relevant_solution_vy;
2400 *
locally_relevant_solution_vz_old = locally_relevant_solution_vz;
2405 * -------------------------------------------------------------------------------
2406 * ------------------------------ MY PETSC WRAPPERS ------------------------------
2407 * -------------------------------------------------------------------------------
2412 *
const std::vector<types::global_dof_index> &indices,
2413 *
std::vector<PetscScalar> &values)
2417 *
PETSc wrapper to get sets of
values from a petsc vector.
2418 * we assume the vector is ghosted
2419 * We need to figure out which elements we
2420 * own locally. Then get a
pointer to the
2421 * elements that are stored here (both the
2422 * ones we own as well as the ghost elements).
2423 * In
this array, the locally owned elements
2424 * come
first followed by the ghost elements whose
2425 * position we can get from an
index set
2431 *
IndexSet ghost_indices = locally_relevant_dofs_LS;
2435 *
n_idx = indices.size();
2437 *
VecGetOwnershipRange (vector, &
begin, &
end);
2439 *
Vec solution_in_local_form =
nullptr;
2440 *
VecGhostGetLocalForm(vector, &solution_in_local_form);
2442 *
PetscScalar *soln;
2443 *
VecGetArray(solution_in_local_form, &soln);
2445 *
for (i = 0; i < n_idx; ++i)
2447 *
int index = indices[i];
2448 *
if (index >=
begin && index <
end)
2452 *
const unsigned int ghostidx = ghost_indices.index_within_set(index);
2456 *
VecRestoreArray(solution_in_local_form, &soln);
2457 *
VecGhostRestoreLocalForm(vector, &solution_in_local_form);
2462 *
const std::vector<types::global_dof_index> &indices,
2463 *
std::map<types::global_dof_index, types::global_dof_index> &map_from_Q1_to_Q2,
2464 *
std::vector<PetscScalar> &values)
2468 * THIS IS MEANT TO BE USED WITH VELOCITY VECTORS
2469 *
PETSc wrapper to get sets of
values from a petsc vector.
2470 * we assume the vector is ghosted
2471 * We need to figure out which elements we
2472 * own locally. Then get a
pointer to the
2473 * elements that are stored here (both the
2474 * ones we own as well as the ghost elements).
2475 * In
this array, the locally owned elements
2476 * come
first followed by the ghost elements whose
2477 * position we can get from an
index set
2483 *
IndexSet ghost_indices = locally_relevant_dofs_U;
2487 *
n_idx = indices.size();
2489 *
VecGetOwnershipRange (vector, &
begin, &
end);
2491 *
Vec solution_in_local_form =
nullptr;
2492 *
VecGhostGetLocalForm(vector, &solution_in_local_form);
2494 *
PetscScalar *soln;
2495 *
VecGetArray(solution_in_local_form, &soln);
2497 *
for (i = 0; i < n_idx; ++i)
2499 *
int index = map_from_Q1_to_Q2[indices[i]];
2500 *
if (index >=
begin && index <
end)
2504 *
const unsigned int ghostidx = ghost_indices.index_within_set(index);
2508 *
VecRestoreArray(solution_in_local_form, &soln);
2509 *
VecGhostRestoreLocalForm(vector, &solution_in_local_form);
2515<a name=
"ann-MultiPhase.cc"></a>
2516<h1>Annotated version of MultiPhase.cc</h1>
2532 *
#include <deal.II/base/quadrature_lib.h>
2533 *
#include <deal.II/base/function.h>
2534 *
#include <deal.II/lac/affine_constraints.h>
2535 *
#include <deal.II/lac/vector.h>
2536 *
#include <deal.II/lac/full_matrix.h>
2537 *
#include <deal.II/lac/solver_cg.h>
2538 *
#include <deal.II/lac/petsc_sparse_matrix.h>
2539 *
#include <deal.II/lac/petsc_vector.h>
2540 *
#include <deal.II/lac/petsc_solver.h>
2541 *
#include <deal.II/lac/petsc_precondition.h>
2542 *
#include <deal.II/grid/grid_generator.h>
2543 *
#include <deal.II/grid/tria_accessor.h>
2544 *
#include <deal.II/grid/tria_iterator.h>
2545 *
#include <deal.II/dofs/dof_handler.h>
2546 *
#include <deal.II/dofs/dof_accessor.h>
2547 *
#include <deal.II/dofs/dof_tools.h>
2548 *
#include <deal.II/fe/fe_values.h>
2549 *
#include <deal.II/fe/fe_q.h>
2550 *
#include <deal.II/numerics/vector_tools.h>
2551 *
#include <deal.II/numerics/data_out.h>
2552 *
#include <deal.II/numerics/error_estimator.h>
2553 *
#include <deal.II/base/utilities.h>
2554 *
#include <deal.II/base/conditional_ostream.h>
2555 *
#include <deal.II/base/index_set.h>
2556 *
#include <deal.II/lac/sparsity_tools.h>
2557 *
#include <deal.II/distributed/tria.h>
2558 *
#include <deal.II/distributed/grid_refinement.h>
2559 *
#include <deal.II/base/convergence_table.h>
2560 *
#include <deal.II/base/timer.h>
2561 *
#include <deal.II/base/parameter_handler.h>
2562 *
#include <fstream>
2563 *
#include <iostream>
2564 *
#include <deal.II/grid/grid_tools.h>
2565 *
#include <deal.II/fe/mapping_q.h>
2567 *
using namespace dealii;
2572 * FOR TRANSPORT PROBLEM
2584 *
#define FILLING_TANK 0
2585 *
#define BREAKING_DAM 1
2586 *
#define FALLING_DROP 2
2587 *
#define SMALL_WAVE_PERTURBATION 3
2589 *
#include
"NavierStokesSolver.cc"
2590 *
#include
"LevelSetSolver.cc"
2591 *
#include
"utilities.cc"
2600 *
template <
int dim>
2604 *
MultiPhase (
const unsigned int degree_LS,
2605 *
const unsigned int degree_U);
2610 *
void set_boundary_inlet();
2611 *
void get_boundary_values_U();
2612 *
void get_boundary_values_phi(std::vector<types::global_dof_index> &boundary_values_id_phi,
2613 *
std::vector<double> &boundary_values_phi);
2614 *
void output_results();
2615 *
void output_vectors();
2616 *
void output_rho();
2618 *
void initial_condition();
2619 *
void init_constraints();
2628 *
IndexSet locally_relevant_dofs_LS;
2635 *
IndexSet locally_relevant_dofs_U;
2640 *
IndexSet locally_relevant_dofs_P;
2662 *
std::vector<types::global_dof_index> boundary_values_id_u;
2663 *
std::vector<types::global_dof_index> boundary_values_id_v;
2664 *
std::vector<types::global_dof_index> boundary_values_id_phi;
2665 *
std::vector<double> boundary_values_u;
2666 *
std::vector<double> boundary_values_v;
2667 *
std::vector<double> boundary_values_phi;
2673 *
double final_time;
2674 *
unsigned int timestep_number;
2680 *
int sharpness_integer;
2682 *
unsigned int n_refinement;
2683 *
unsigned int output_number;
2684 *
double output_time;
2708 *
unsigned int TRANSPORT_TIME_INTEGRATION;
2709 *
std::string ALGORITHM;
2710 *
unsigned int PROBLEM;
2713 *
template <
int dim>
2714 *
MultiPhase<dim>::MultiPhase (
const unsigned int degree_LS,
2715 *
const unsigned int degree_U)
2717 *
mpi_communicator (MPI_COMM_WORLD),
2718 *
triangulation (mpi_communicator,
2722 *
degree_LS(degree_LS),
2723 *
dof_handler_LS (triangulation),
2724 *
fe_LS (degree_LS),
2725 *
degree_U(degree_U),
2726 *
dof_handler_U (triangulation),
2728 *
dof_handler_P (triangulation),
2729 *
fe_P (degree_U-1),
2733 *
template <
int dim>
2734 *
MultiPhase<dim>::~MultiPhase ()
2736 *
dof_handler_LS.clear ();
2737 *
dof_handler_U.clear ();
2738 *
dof_handler_P.clear ();
2748 *
template <
int dim>
2749 *
void MultiPhase<dim>::setup()
2756 *
dof_handler_LS.distribute_dofs (fe_LS);
2757 *
locally_owned_dofs_LS = dof_handler_LS.locally_owned_dofs ();
2764 *
dof_handler_U.distribute_dofs (fe_U);
2765 *
locally_owned_dofs_U = dof_handler_U.locally_owned_dofs ();
2772 *
dof_handler_P.distribute_dofs (fe_P);
2773 *
locally_owned_dofs_P = dof_handler_P.locally_owned_dofs ();
2777 * init vectors
for phi
2780 *
locally_relevant_solution_phi.reinit(locally_owned_dofs_LS,locally_relevant_dofs_LS,mpi_communicator);
2781 *
locally_relevant_solution_phi = 0;
2782 *
completely_distributed_solution_phi.reinit (locally_owned_dofs_P,mpi_communicator);
2785 * init vectors
for u
2788 *
locally_relevant_solution_u.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
2789 *
locally_relevant_solution_u = 0;
2790 *
completely_distributed_solution_u.reinit (locally_owned_dofs_U,mpi_communicator);
2793 * init vectors
for v
2796 *
locally_relevant_solution_v.reinit (locally_owned_dofs_U,locally_relevant_dofs_U,mpi_communicator);
2797 *
locally_relevant_solution_v = 0;
2798 *
completely_distributed_solution_v.reinit (locally_owned_dofs_U,mpi_communicator);
2801 * init vectors
for p
2804 *
locally_relevant_solution_p.reinit (locally_owned_dofs_P,locally_relevant_dofs_P,mpi_communicator);
2805 *
locally_relevant_solution_p = 0;
2806 *
completely_distributed_solution_p.reinit (locally_owned_dofs_P,mpi_communicator);
2812 *
init_constraints();
2815 *
template <
int dim>
2816 *
void MultiPhase<dim>::initial_condition()
2821 * Initial conditions
2822 * init condition
for phi
2825 *
completely_distributed_solution_phi = 0;
2827 *
InitialPhi<dim>(PROBLEM, sharpness),
2828 *
completely_distributed_solution_phi);
2829 *
constraints.distribute (completely_distributed_solution_phi);
2830 *
locally_relevant_solution_phi = completely_distributed_solution_phi;
2833 * init condition
for u=0
2836 *
completely_distributed_solution_u = 0;
2839 *
completely_distributed_solution_u);
2840 *
constraints.distribute (completely_distributed_solution_u);
2841 *
locally_relevant_solution_u = completely_distributed_solution_u;
2844 * init condition
for v
2847 *
completely_distributed_solution_v = 0;
2850 *
completely_distributed_solution_v);
2851 *
constraints.distribute (completely_distributed_solution_v);
2852 *
locally_relevant_solution_v = completely_distributed_solution_v;
2855 * init condition
for p
2858 *
completely_distributed_solution_p = 0;
2861 *
completely_distributed_solution_p);
2862 *
constraints.distribute (completely_distributed_solution_p);
2863 *
locally_relevant_solution_p = completely_distributed_solution_p;
2866 *
template <
int dim>
2867 *
void MultiPhase<dim>::init_constraints()
2869 *
constraints.clear ();
2870 *
constraints.reinit (locally_owned_dofs_LS, locally_relevant_dofs_LS);
2872 *
constraints.close ();
2875 *
template <
int dim>
2876 *
void MultiPhase<dim>::get_boundary_values_U()
2878 *
std::map<types::global_dof_index, double> map_boundary_values_u;
2879 *
std::map<types::global_dof_index, double> map_boundary_values_v;
2880 *
std::map<types::global_dof_index, double> map_boundary_values_w;
2887 *
if (PROBLEM==BREAKING_DAM || PROBLEM==FALLING_DROP)
2918 *
else if (PROBLEM==SMALL_WAVE_PERTURBATION)
2922 * no slip in bottom and top and slip in left and right
2948 *
else if (PROBLEM==FILLING_TANK)
2952 * LEFT: entry in x,
zero in y
2959 * RIGHT: no-slip condition
2973 * TOP: exit in y,
zero in x
2981 *
pcout <<
"Error in type of PROBLEM at Boundary Conditions" << std::endl;
2984 *
boundary_values_id_u.resize(map_boundary_values_u.size());
2985 *
boundary_values_id_v.resize(map_boundary_values_v.size());
2986 *
boundary_values_u.resize(map_boundary_values_u.size());
2987 *
boundary_values_v.resize(map_boundary_values_v.size());
2988 *
std::map<types::global_dof_index,double>::const_iterator boundary_value_u =map_boundary_values_u.begin();
2989 *
std::map<types::global_dof_index,double>::const_iterator boundary_value_v =map_boundary_values_v.begin();
2991 *
for (
int i=0; boundary_value_u !=map_boundary_values_u.end(); ++boundary_value_u, ++i)
2993 *
boundary_values_id_u[i]=boundary_value_u->first;
2994 *
boundary_values_u[i]=boundary_value_u->second;
2996 *
for (
int i=0; boundary_value_v !=map_boundary_values_v.end(); ++boundary_value_v, ++i)
2998 *
boundary_values_id_v[i]=boundary_value_v->first;
2999 *
boundary_values_v[i]=boundary_value_v->second;
3003 *
template <
int dim>
3004 *
void MultiPhase<dim>::set_boundary_inlet()
3006 *
const QGauss<dim-1> face_quadrature_formula(1);
3010 *
const unsigned int n_face_q_points = face_quadrature_formula.size();
3011 *
std::vector<double> u_value (n_face_q_points);
3012 *
std::vector<double> v_value (n_face_q_points);
3015 *
cell_U = dof_handler_U.begin_active(),
3016 *
endc_U = dof_handler_U.end();
3019 *
for (; cell_U!=endc_U; ++cell_U)
3020 *
if (cell_U->is_locally_owned())
3021 *
for (
unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
3022 *
if (cell_U->face(face)->at_boundary())
3024 *
fe_face_values.reinit(cell_U,face);
3025 *
fe_face_values.get_function_values(locally_relevant_solution_u,u_value);
3026 *
fe_face_values.get_function_values(locally_relevant_solution_v,v_value);
3029 *
if (fe_face_values.normal_vector(0)*u < -1e-14)
3030 *
cell_U->face(face)->set_boundary_id(10);
3034 *
template <
int dim>
3035 *
void MultiPhase<dim>::get_boundary_values_phi(std::vector<types::global_dof_index> &boundary_values_id_phi,
3036 *
std::vector<double> &boundary_values_phi)
3038 *
std::map<types::global_dof_index, double> map_boundary_values_phi;
3041 *
set_boundary_inlet();
3044 *
boundary_values_id_phi.resize(map_boundary_values_phi.size());
3045 *
boundary_values_phi.resize(map_boundary_values_phi.size());
3046 *
std::map<types::global_dof_index,double>::const_iterator boundary_value_phi = map_boundary_values_phi.begin();
3047 *
for (
int i=0; boundary_value_phi !=map_boundary_values_phi.end(); ++boundary_value_phi, ++i)
3049 *
boundary_values_id_phi[i]=boundary_value_phi->first;
3050 *
boundary_values_phi[i]=boundary_value_phi->second;
3055 *
void MultiPhase<dim>::output_results()
3066 *
template <
int dim>
3067 *
void MultiPhase<dim>::output_vectors()
3071 *
data_out.add_data_vector (locally_relevant_solution_phi,
"phi");
3072 *
data_out.build_patches ();
3074 *
const std::string filename = (
"sol_vectors-" +
3078 *
(triangulation.locally_owned_subdomain(), 4));
3079 *
std::ofstream output ((filename +
".vtu").c_str());
3080 *
data_out.write_vtu (output);
3084 *
std::vector<std::string> filenames;
3085 *
for (
unsigned int i=0;
3088 *
filenames.push_back (
"sol_vectors-" +
3094 *
std::ofstream master_output ((filename +
".pvtu").c_str());
3095 *
data_out.write_pvtu_record (master_output, filenames);
3099 *
template <
int dim>
3100 *
void MultiPhase<dim>::output_rho()
3102 *
Postprocessor<dim> postprocessor(eps,rho_air,rho_fluid);
3105 *
data_out.add_data_vector (locally_relevant_solution_phi, postprocessor);
3107 *
data_out.build_patches ();
3109 *
const std::string filename = (
"sol_rho-" +
3113 *
(triangulation.locally_owned_subdomain(), 4));
3114 *
std::ofstream output ((filename +
".vtu").c_str());
3115 *
data_out.write_vtu (output);
3119 *
std::vector<std::string> filenames;
3120 *
for (
unsigned int i=0;
3123 *
filenames.push_back (
"sol_rho-" +
3129 *
std::ofstream master_output ((filename +
".pvtu").c_str());
3130 *
data_out.write_pvtu_record (master_output, filenames);
3134 *
template <
int dim>
3135 *
void MultiPhase<dim>::run()
3140 * GENERAL PARAMETERS
3147 *
get_output =
true;
3148 *
output_number = 0;
3150 *
output_time = 0.1;
3151 *
final_time = 10.0;
3155 * PARAMETERS FOR THE NAVIER STOKES PROBLEM
3159 *
rho_fluid = 1000.;
3163 *
PROBLEM=BREAKING_DAM;
3166 * PROBLEM=FILLING_TANK;
3167 * PROBLEM=SMALL_WAVE_PERTURBATION;
3168 * PROBLEM=FALLING_DROP;
3174 *
ForceTerms<dim> force_function(std::vector<double> {0.0,-1.0});
3178 * PARAMETERS FOR TRANSPORT PROBLEM
3184 *
sharpness_integer=10;
3190 *
TRANSPORT_TIME_INTEGRATION=SSP33;
3193 * ALGORITHM =
"MPP_u1";
3194 * ALGORITHM =
"NMPP_uH";
3197 *
ALGORITHM =
"MPP_uH";
3201 * ADJUST PARAMETERS ACCORDING TO PROBLEM
3204 *
if (PROBLEM==FALLING_DROP)
3214 *
if (PROBLEM==FILLING_TANK)
3217 *
else if (PROBLEM==BREAKING_DAM || PROBLEM==SMALL_WAVE_PERTURBATION)
3219 *
std::vector< unsigned int > repetitions;
3220 *
repetitions.push_back(2);
3221 *
repetitions.push_back(1);
3225 *
else if (PROBLEM==FALLING_DROP)
3227 *
std::vector< unsigned int > repetitions;
3228 *
repetitions.push_back(1);
3229 *
repetitions.push_back(4);
3233 *
triangulation.refine_global (n_refinement);
3243 * PARAMETERS FOR TIME STEPPING
3247 *
time_step = cfl*min_h/umax;
3249 *
sharpness=sharpness_integer*min_h;
3253 * INITIAL CONDITIONS
3256 *
initial_condition();
3261 * NAVIER STOKES SOLVER
3264 *
NavierStokesSolver<dim> navier_stokes (degree_LS,degree_U,
3267 *
rho_fluid,nu_fluid,
3270 *
triangulation,mpi_communicator);
3273 * BOUNDARY CONDITIONS FOR NAVIER STOKES
3276 *
get_boundary_values_U();
3277 *
navier_stokes.set_boundary_conditions(boundary_values_id_u, boundary_values_id_v,
3278 *
boundary_values_u, boundary_values_v);
3282 *
set INITIAL CONDITION within NAVIER STOKES
3285 *
navier_stokes.initial_condition(locally_relevant_solution_phi,
3286 *
locally_relevant_solution_u,
3287 *
locally_relevant_solution_v,
3288 *
locally_relevant_solution_p);
3294 *
LevelSetSolver<dim> transport_solver (degree_LS,degree_U,
3298 *
TRANSPORT_TIME_INTEGRATION,
3300 *
mpi_communicator);
3303 * BOUNDARY CONDITIONS FOR PHI
3306 *
get_boundary_values_phi(boundary_values_id_phi,boundary_values_phi);
3307 *
transport_solver.set_boundary_conditions(boundary_values_id_phi,boundary_values_phi);
3311 *
set INITIAL CONDITION within TRANSPORT PROBLEM
3314 *
transport_solver.initial_condition(locally_relevant_solution_phi,
3315 *
locally_relevant_solution_u,
3316 *
locally_relevant_solution_v);
3317 *
int dofs_U = 2*dof_handler_U.n_dofs();
3318 *
int dofs_P = 2*dof_handler_P.n_dofs();
3319 *
int dofs_LS = dof_handler_LS.n_dofs();
3320 *
int dofs_TOTAL = dofs_U+dofs_P+dofs_LS;
3324 * NO BOUNDARY CONDITIONS
for LEVEL SET
3327 *
pcout <<
"Cfl: " << cfl <<
"; umax: " << umax <<
"; min h: " << min_h
3328 *
<<
"; time step: " << time_step << std::endl;
3329 *
pcout <<
" Number of active cells: "
3330 *
<< triangulation.n_global_active_cells() << std::endl
3331 *
<<
" Number of degrees of freedom: " << std::endl
3332 *
<<
" U: " << dofs_U << std::endl
3333 *
<<
" P: " << dofs_P << std::endl
3334 *
<<
" LS: " << dofs_LS << std::endl
3335 *
<<
" TOTAL: " << dofs_TOTAL
3343 *
for (timestep_number=1, time=time_step; time<=final_time;
3344 *
time+=time_step,++timestep_number)
3346 *
pcout <<
"Time step " << timestep_number
3347 *
<<
" at t=" << time
3351 * GET NAVIER STOKES VELOCITY
3354 *
navier_stokes.set_phi(locally_relevant_solution_phi);
3355 *
navier_stokes.nth_time_step();
3356 *
navier_stokes.get_velocity(locally_relevant_solution_u,locally_relevant_solution_v);
3357 *
transport_solver.set_velocity(locally_relevant_solution_u,locally_relevant_solution_v);
3360 * GET LEVEL SET SOLUTION
3363 *
transport_solver.nth_time_step();
3364 *
transport_solver.get_unp1(locally_relevant_solution_phi);
3365 *
if (get_output && time-(output_number)*output_time>0)
3368 *
navier_stokes.get_velocity(locally_relevant_solution_u, locally_relevant_solution_v);
3369 *
transport_solver.get_unp1(locally_relevant_solution_phi);
3374 *
int main(
int argc,
char *argv[])
3378 *
using namespace dealii;
3382 *
unsigned int degree_LS = 1;
3383 *
unsigned int degree_U = 2;
3384 *
MultiPhase<2> multi_phase(degree_LS, degree_U);
3385 *
multi_phase.run();
3389 *
catch (std::exception &exc)
3391 *
std::cerr << std::endl << std::endl
3392 *
<<
"----------------------------------------------------"
3394 *
std::cerr <<
"Exception on processing: " << std::endl
3395 *
<< exc.what() << std::endl
3396 *
<<
"Aborting!" << std::endl
3397 *
<<
"----------------------------------------------------"
3403 *
std::cerr << std::endl << std::endl
3404 *
<<
"----------------------------------------------------"
3406 *
std::cerr <<
"Unknown exception!" << std::endl
3407 *
<<
"Aborting!" << std::endl
3408 *
<<
"----------------------------------------------------"
3417<a name=
"ann-NavierStokesSolver.cc"></a>
3418<h1>Annotated version of NavierStokesSolver.cc</h1>
3434 *
#include <deal.II/base/quadrature_lib.h>
3435 *
#include <deal.II/base/function.h>
3436 *
#include <deal.II/lac/affine_constraints.h>
3437 *
#include <deal.II/lac/vector.h>
3438 *
#include <deal.II/lac/full_matrix.h>
3439 *
#include <deal.II/lac/solver_cg.h>
3440 *
#include <deal.II/lac/petsc_sparse_matrix.h>
3441 *
#include <deal.II/lac/petsc_vector.h>
3442 *
#include <deal.II/lac/petsc_solver.h>
3443 *
#include <deal.II/lac/petsc_precondition.h>
3444 *
#include <deal.II/grid/grid_generator.h>
3445 *
#include <deal.II/grid/tria_accessor.h>
3446 *
#include <deal.II/grid/tria_iterator.h>
3447 *
#include <deal.II/dofs/dof_handler.h>
3448 *
#include <deal.II/dofs/dof_accessor.h>
3449 *
#include <deal.II/dofs/dof_tools.h>
3450 *
#include <deal.II/fe/fe_values.h>
3451 *
#include <deal.II/fe/fe_q.h>
3452 *
#include <deal.II/numerics/vector_tools.h>
3453 *
#include <deal.II/numerics/data_out.h>
3454 *
#include <deal.II/numerics/error_estimator.h>
3455 *
#include <deal.II/base/utilities.h>
3456 *
#include <deal.II/base/conditional_ostream.h>
3457 *
#include <deal.II/base/index_set.h>
3458 *
#include <deal.II/lac/sparsity_tools.h>
3459 *
#include <deal.II/distributed/tria.h>
3460 *
#include <deal.II/distributed/grid_refinement.h>
3461 *
#include <deal.II/lac/petsc_vector.h>
3462 *
#include <deal.II/base/convergence_table.h>
3463 *
#include <deal.II/base/timer.h>
3464 *
#include <deal.II/base/parameter_handler.h>
3465 *
#include <deal.II/grid/grid_tools.h>
3466 *
#include <deal.II/fe/mapping_q.h>
3468 *
#include <fstream>
3469 *
#include <iostream>
3472 *
using namespace dealii;
3474 *
#define MAX_NUM_ITER_TO_RECOMPUTE_PRECONDITIONER 10
3484 *
class NavierStokesSolver
3489 * constructor
for using LEVEL SET
3492 *
NavierStokesSolver(
const unsigned int degree_LS,
3493 *
const unsigned int degree_U,
3494 *
const double time_step,
3496 *
const double rho_air,
3497 *
const double nu_air,
3498 *
const double rho_fluid,
3499 *
const double nu_fluid,
3501 *
const bool verbose,
3506 * constructor
for NOT LEVEL SET
3509 *
NavierStokesSolver(
const unsigned int degree_LS,
3510 *
const unsigned int degree_U,
3511 *
const double time_step,
3515 *
const bool verbose,
3524 *
void set_rho_and_nu_functions(
const Function<dim> &rho_function,
3542 * boundary conditions
3545 *
void set_boundary_conditions(std::vector<types::global_dof_index> boundary_values_id_u,
3546 *
std::vector<types::global_dof_index> boundary_values_id_v, std::vector<double> boundary_values_u,
3547 *
std::vector<double> boundary_values_v);
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,
3550 *
std::vector<types::global_dof_index> boundary_values_id_w, std::vector<double> boundary_values_u,
3551 *
std::vector<double> boundary_values_v, std::vector<double> boundary_values_w);
3569 *
void nth_time_step();
3577 *
~NavierStokesSolver();
3582 * SETUP AND INITIAL CONDITION
3586 *
void setup_VECTORS();
3587 *
void init_constraints();
3593 *
void assemble_system_U();
3594 *
void assemble_system_dpsi_q();
3601 *
std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
3605 *
std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
3610 * GET DIFFERENT FIELDS
3613 *
void get_rho_and_nu(
double phi);
3615 *
void get_pressure();
3621 *
void save_old_solution();
3630 *
IndexSet locally_relevant_dofs_LS;
3636 *
IndexSet locally_relevant_dofs_U;
3641 *
IndexSet locally_relevant_dofs_P;
3656 *
unsigned int LEVEL_SET;
3657 *
unsigned int RHO_TIMES_RHS;
3673 *
std::vector<types::global_dof_index> boundary_values_id_u;
3674 *
std::vector<types::global_dof_index> boundary_values_id_v;
3675 *
std::vector<types::global_dof_index> boundary_values_id_w;
3676 *
std::vector<double> boundary_values_u;
3677 *
std::vector<double> boundary_values_v;
3678 *
std::vector<double> boundary_values_w;
3683 *
bool rebuild_Matrix_U;
3684 *
std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_Matrix_u;
3685 *
std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_Matrix_v;
3686 *
std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_Matrix_w;
3688 *
std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_S;
3690 *
std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_M;
3692 *
bool rebuild_Matrix_U_preconditioners;
3693 *
bool rebuild_S_M_preconditioners;
3721 * CONSTRUCTOR FOR LEVEL SET
3725 *
NavierStokesSolver<dim>::NavierStokesSolver(
const unsigned int degree_LS,
3726 *
const unsigned int degree_U,
3727 *
const double time_step,
3729 *
const double rho_air,
3730 *
const double nu_air,
3731 *
const double rho_fluid,
3732 *
const double nu_fluid,
3734 *
const bool verbose,
3738 *
mpi_communicator(mpi_communicator),
3739 *
triangulation(triangulation),
3740 *
degree_LS(degree_LS),
3741 *
dof_handler_LS(triangulation),
3743 *
degree_U(degree_U),
3744 *
dof_handler_U(triangulation),
3746 *
dof_handler_P(triangulation),
3748 *
force_function(force_function),
3751 * This is dummy since rho and nu
functions won
't be used
3754 * rho_function(force_function),
3755 * nu_function(force_function),
3758 * rho_fluid(rho_fluid),
3759 * nu_fluid(nu_fluid),
3760 * time_step(time_step),
3765 * pcout(std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator)==0)),
3766 * rebuild_Matrix_U(true),
3767 * rebuild_S_M(true),
3768 * rebuild_Matrix_U_preconditioners(true),
3769 * rebuild_S_M_preconditioners(true)
3774 * CONSTRUCTOR NOT FOR LEVEL SET
3778 * NavierStokesSolver<dim>::NavierStokesSolver(const unsigned int degree_LS,
3779 * const unsigned int degree_U,
3780 * const double time_step,
3781 * Function<dim> &force_function,
3782 * Function<dim> &rho_function,
3783 * Function<dim> &nu_function,
3784 * const bool verbose,
3785 * parallel::distributed::Triangulation<dim> &triangulation,
3786 * MPI_Comm &mpi_communicator) :
3787 * mpi_communicator(mpi_communicator),
3788 * triangulation(triangulation),
3789 * degree_LS(degree_LS),
3790 * dof_handler_LS(triangulation),
3792 * degree_U(degree_U),
3793 * dof_handler_U(triangulation),
3795 * dof_handler_P(triangulation),
3797 * force_function(force_function),
3798 * rho_function(rho_function),
3799 * nu_function(nu_function),
3800 * time_step(time_step),
3804 * pcout(std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator)==0)),
3805 * rebuild_Matrix_U(true),
3806 * rebuild_S_M(true),
3807 * rebuild_Matrix_U_preconditioners(true),
3808 * rebuild_S_M_preconditioners(true)
3812 * NavierStokesSolver<dim>::~NavierStokesSolver()
3814 * dof_handler_LS.clear();
3815 * dof_handler_U.clear();
3816 * dof_handler_P.clear();
3821 * /////////////////////////////////////////////////////////
3822 * ////////////////// SETTERS AND GETTERS //////////////////
3823 * /////////////////////////////////////////////////////////
3827 * void NavierStokesSolver<dim>::set_rho_and_nu_functions(const Function<dim> &rho_function,
3828 * const Function<dim> &nu_function)
3830 * this->rho_function=rho_function;
3831 * this->nu_function=nu_function;
3835 * void NavierStokesSolver<dim>::initial_condition(PETScWrappers::MPI::Vector locally_relevant_solution_phi,
3836 * PETScWrappers::MPI::Vector locally_relevant_solution_u,
3837 * PETScWrappers::MPI::Vector locally_relevant_solution_v,
3838 * PETScWrappers::MPI::Vector locally_relevant_solution_p)
3840 * this->locally_relevant_solution_phi=locally_relevant_solution_phi;
3841 * this->locally_relevant_solution_u=locally_relevant_solution_u;
3842 * this->locally_relevant_solution_v=locally_relevant_solution_v;
3843 * this->locally_relevant_solution_p=locally_relevant_solution_p;
3846 * set old vectors to the initial condition (just for first time step)
3849 * save_old_solution();
3853 * void NavierStokesSolver<dim>::initial_condition(PETScWrappers::MPI::Vector locally_relevant_solution_phi,
3854 * PETScWrappers::MPI::Vector locally_relevant_solution_u,
3855 * PETScWrappers::MPI::Vector locally_relevant_solution_v,
3856 * PETScWrappers::MPI::Vector locally_relevant_solution_w,
3857 * PETScWrappers::MPI::Vector locally_relevant_solution_p)
3859 * this->locally_relevant_solution_phi=locally_relevant_solution_phi;
3860 * this->locally_relevant_solution_u=locally_relevant_solution_u;
3861 * this->locally_relevant_solution_v=locally_relevant_solution_v;
3862 * this->locally_relevant_solution_w=locally_relevant_solution_w;
3863 * this->locally_relevant_solution_p=locally_relevant_solution_p;
3866 * set old vectors to the initial condition (just for first time step)
3869 * save_old_solution();
3873 * void NavierStokesSolver<dim>::set_boundary_conditions(std::vector<types::global_dof_index> boundary_values_id_u,
3874 * std::vector<types::global_dof_index> boundary_values_id_v,
3875 * std::vector<double> boundary_values_u,
3876 * std::vector<double> boundary_values_v)
3878 * this->boundary_values_id_u=boundary_values_id_u;
3879 * this->boundary_values_id_v=boundary_values_id_v;
3880 * this->boundary_values_u=boundary_values_u;
3881 * this->boundary_values_v=boundary_values_v;
3885 * void NavierStokesSolver<dim>::set_boundary_conditions(std::vector<types::global_dof_index> boundary_values_id_u,
3886 * std::vector<types::global_dof_index> boundary_values_id_v,
3887 * std::vector<types::global_dof_index> boundary_values_id_w,
3888 * std::vector<double> boundary_values_u,
3889 * std::vector<double> boundary_values_v,
3890 * std::vector<double> boundary_values_w)
3892 * this->boundary_values_id_u=boundary_values_id_u;
3893 * this->boundary_values_id_v=boundary_values_id_v;
3894 * this->boundary_values_id_w=boundary_values_id_w;
3895 * this->boundary_values_u=boundary_values_u;
3896 * this->boundary_values_v=boundary_values_v;
3897 * this->boundary_values_w=boundary_values_w;
3901 * void NavierStokesSolver<dim>::set_velocity(PETScWrappers::MPI::Vector locally_relevant_solution_u,
3902 * PETScWrappers::MPI::Vector locally_relevant_solution_v)
3904 * this->locally_relevant_solution_u=locally_relevant_solution_u;
3905 * this->locally_relevant_solution_v=locally_relevant_solution_v;
3909 * void NavierStokesSolver<dim>::set_velocity(PETScWrappers::MPI::Vector locally_relevant_solution_u,
3910 * PETScWrappers::MPI::Vector locally_relevant_solution_v,
3911 * PETScWrappers::MPI::Vector locally_relevant_solution_w)
3913 * this->locally_relevant_solution_u=locally_relevant_solution_u;
3914 * this->locally_relevant_solution_v=locally_relevant_solution_v;
3915 * this->locally_relevant_solution_w=locally_relevant_solution_w;
3919 * void NavierStokesSolver<dim>::set_phi(PETScWrappers::MPI::Vector locally_relevant_solution_phi)
3921 * this->locally_relevant_solution_phi=locally_relevant_solution_phi;
3925 * void NavierStokesSolver<dim>::get_rho_and_nu(double phi)
3935 * else if (phi<-eps)
3939 * rho_value=rho_fluid*(1+H)/2.+rho_air*(1-H)/2.;
3940 * nu_value=nu_fluid*(1+H)/2.+nu_air*(1-H)/2.;
3943 * rho_value=rho_fluid*(1+phi)/2.+rho_air*(1-phi)/2.;
3944 * nu_value=nu_fluid*(1+phi)/2.+nu_air*(1-phi)/2.;
3950 * void NavierStokesSolver<dim>::get_pressure(PETScWrappers::MPI::Vector &locally_relevant_solution_p)
3952 * locally_relevant_solution_p=this->locally_relevant_solution_p;
3956 * void NavierStokesSolver<dim>::get_velocity(PETScWrappers::MPI::Vector &locally_relevant_solution_u,
3957 * PETScWrappers::MPI::Vector &locally_relevant_solution_v)
3959 * locally_relevant_solution_u=this->locally_relevant_solution_u;
3960 * locally_relevant_solution_v=this->locally_relevant_solution_v;
3964 * void NavierStokesSolver<dim>::get_velocity(PETScWrappers::MPI::Vector &locally_relevant_solution_u,
3965 * PETScWrappers::MPI::Vector &locally_relevant_solution_v,
3966 * PETScWrappers::MPI::Vector &locally_relevant_solution_w)
3968 * locally_relevant_solution_u=this->locally_relevant_solution_u;
3969 * locally_relevant_solution_v=this->locally_relevant_solution_v;
3970 * locally_relevant_solution_w=this->locally_relevant_solution_w;
3975 * ///////////////////////////////////////////////////
3976 * /////////// SETUP AND INITIAL CONDITION ///////////
3977 * ///////////////////////////////////////////////////
3981 * void NavierStokesSolver<dim>::setup()
3983 * pcout<<"***** SETUP IN NAVIER STOKES SOLVER *****"<<std::endl;
3985 * init_constraints();
3990 * void NavierStokesSolver<dim>::setup_DOF()
3993 * degree_MAX=std::max(degree_LS,degree_U);
3999 * dof_handler_LS.distribute_dofs(fe_LS);
4000 * locally_owned_dofs_LS = dof_handler_LS.locally_owned_dofs();
4001 * locally_relevant_dofs_LS = DoFTools::extract_locally_relevant_dofs(dof_handler_LS);
4007 * dof_handler_U.distribute_dofs(fe_U);
4008 * locally_owned_dofs_U = dof_handler_U.locally_owned_dofs();
4009 * locally_relevant_dofs_U = DoFTools::extract_locally_relevant_dofs(dof_handler_U);
4015 * dof_handler_P.distribute_dofs(fe_P);
4016 * locally_owned_dofs_P = dof_handler_P.locally_owned_dofs();
4017 * locally_relevant_dofs_P = DoFTools::extract_locally_relevant_dofs(dof_handler_P);
4021 * void NavierStokesSolver<dim>::setup_VECTORS()
4025 * init vectors for phi
4028 * locally_relevant_solution_phi.reinit(locally_owned_dofs_LS,locally_relevant_dofs_LS,
4029 * mpi_communicator);
4030 * locally_relevant_solution_phi=0;
4033 * init vectors for u
4036 * locally_relevant_solution_u.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4037 * mpi_communicator);
4038 * locally_relevant_solution_u=0;
4039 * completely_distributed_solution_u.reinit(locally_owned_dofs_U,mpi_communicator);
4040 * system_rhs_u.reinit(locally_owned_dofs_U,mpi_communicator);
4043 * init vectors for u_old
4046 * locally_relevant_solution_u_old.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4047 * mpi_communicator);
4048 * locally_relevant_solution_u_old=0;
4051 * init vectors for v
4054 * locally_relevant_solution_v.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4055 * mpi_communicator);
4056 * locally_relevant_solution_v=0;
4057 * completely_distributed_solution_v.reinit(locally_owned_dofs_U,mpi_communicator);
4058 * system_rhs_v.reinit(locally_owned_dofs_U,mpi_communicator);
4061 * init vectors for v_old
4064 * locally_relevant_solution_v_old.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4065 * mpi_communicator);
4066 * locally_relevant_solution_v_old=0;
4069 * init vectors for w
4072 * locally_relevant_solution_w.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4073 * mpi_communicator);
4074 * locally_relevant_solution_w=0;
4075 * completely_distributed_solution_w.reinit(locally_owned_dofs_U,mpi_communicator);
4076 * system_rhs_w.reinit(locally_owned_dofs_U,mpi_communicator);
4079 * init vectors for w_old
4082 * locally_relevant_solution_w_old.reinit(locally_owned_dofs_U,locally_relevant_dofs_U,
4083 * mpi_communicator);
4084 * locally_relevant_solution_w_old=0;
4087 * init vectors for dpsi
4090 * locally_relevant_solution_psi.reinit(locally_owned_dofs_P,locally_relevant_dofs_P,
4091 * mpi_communicator);
4092 * locally_relevant_solution_psi=0;
4093 * system_rhs_psi.reinit(locally_owned_dofs_P,mpi_communicator);
4096 * init vectors for dpsi old
4099 * locally_relevant_solution_psi_old.reinit(locally_owned_dofs_P,locally_relevant_dofs_P,
4100 * mpi_communicator);
4101 * locally_relevant_solution_psi_old=0;
4104 * init vectors for q
4107 * completely_distributed_solution_q.reinit(locally_owned_dofs_P,mpi_communicator);
4108 * system_rhs_q.reinit(locally_owned_dofs_P,mpi_communicator);
4111 * init vectors for psi
4114 * completely_distributed_solution_psi.reinit(locally_owned_dofs_P,mpi_communicator);
4117 * init vectors for p
4120 * locally_relevant_solution_p.reinit(locally_owned_dofs_P,locally_relevant_dofs_P,
4121 * mpi_communicator);
4122 * locally_relevant_solution_p=0;
4123 * completely_distributed_solution_p.reinit(locally_owned_dofs_P,mpi_communicator);
4126 * ////////////////////////
4127 * Initialize constraints
4128 * ////////////////////////
4131 * init_constraints();
4134 * //////////////////
4136 * //////////////////
4137 * sparsity pattern for A
4140 * DynamicSparsityPattern dsp_Matrix(locally_relevant_dofs_U);
4141 * DoFTools::make_sparsity_pattern(dof_handler_U,dsp_Matrix,constraints,false);
4142 * SparsityTools::distribute_sparsity_pattern(dsp_Matrix,
4143 * dof_handler_U.locally_owned_dofs(),
4145 * locally_relevant_dofs_U);
4146 * system_Matrix_u.reinit(dof_handler_U.locally_owned_dofs(),
4147 * dof_handler_U.locally_owned_dofs(),
4149 * mpi_communicator);
4150 * system_Matrix_v.reinit(dof_handler_U.locally_owned_dofs(),
4151 * dof_handler_U.locally_owned_dofs(),
4153 * mpi_communicator);
4154 * system_Matrix_w.reinit(dof_handler_U.locally_owned_dofs(),
4155 * dof_handler_U.locally_owned_dofs(),
4157 * mpi_communicator);
4158 * rebuild_Matrix_U=true;
4161 * sparsity pattern for S
4164 * DynamicSparsityPattern dsp_S(locally_relevant_dofs_P);
4165 * DoFTools::make_sparsity_pattern(dof_handler_P,dsp_S,constraints_psi,false);
4166 * SparsityTools::distribute_sparsity_pattern(dsp_S,
4167 * dof_handler_P.locally_owned_dofs(),
4169 * locally_relevant_dofs_P);
4170 * system_S.reinit(dof_handler_P.locally_owned_dofs(),
4171 * dof_handler_P.locally_owned_dofs(),
4173 * mpi_communicator);
4176 * sparsity pattern for M
4179 * DynamicSparsityPattern dsp_M(locally_relevant_dofs_P);
4180 * DoFTools::make_sparsity_pattern(dof_handler_P,dsp_M,constraints_psi,false);
4181 * SparsityTools::distribute_sparsity_pattern(dsp_M,
4182 * dof_handler_P.locally_owned_dofs(),
4184 * locally_relevant_dofs_P);
4185 * system_M.reinit(dof_handler_P.locally_owned_dofs(),
4186 * dof_handler_P.locally_owned_dofs(),
4188 * mpi_communicator);
4193 * void NavierStokesSolver<dim>::init_constraints()
4200 * constraints.clear();
4201 * constraints.reinit(locally_owned_dofs_U, locally_relevant_dofs_U);
4202 * DoFTools::make_hanging_node_constraints(dof_handler_U,constraints);
4203 * constraints.close();
4206 * constraints for dpsi
4209 * constraints_psi.clear();
4210 * constraints_psi.reinit(locally_owned_dofs_P, locally_relevant_dofs_P);
4211 * DoFTools::make_hanging_node_constraints(dof_handler_P,constraints_psi);
4214 * if (constraints_psi.can_store_line(0))
4215 * constraints_psi.add_line(0); //constraint u0 = 0
4218 * constraints_psi.close();
4223 * ///////////////////////////////////////////////////
4224 * //////////////// ASSEMBLE SYSTEMS /////////////////
4225 * ///////////////////////////////////////////////////
4229 * void NavierStokesSolver<dim>::assemble_system_U()
4231 * if (rebuild_Matrix_U==true)
4233 * system_Matrix_u=0;
4234 * system_Matrix_v=0;
4235 * system_Matrix_w=0;
4241 * const QGauss<dim> quadrature_formula(degree_MAX+1);
4242 * FEValues<dim> fe_values_LS(fe_LS,quadrature_formula,
4243 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4244 * FEValues<dim> fe_values_U(fe_U,quadrature_formula,
4245 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4246 * FEValues<dim> fe_values_P(fe_P,quadrature_formula,
4247 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4249 * const unsigned int dofs_per_cell=fe_U.dofs_per_cell;
4250 * const unsigned int n_q_points=quadrature_formula.size();
4252 * FullMatrix<double> cell_A_u(dofs_per_cell,dofs_per_cell);
4253 * Vector<double> cell_rhs_u(dofs_per_cell);
4254 * Vector<double> cell_rhs_v(dofs_per_cell);
4255 * Vector<double> cell_rhs_w(dofs_per_cell);
4257 * std::vector<double> phiqnp1(n_q_points);
4259 * std::vector<double> uqn(n_q_points);
4260 * std::vector<double> uqnm1(n_q_points);
4261 * std::vector<double> vqn(n_q_points);
4262 * std::vector<double> vqnm1(n_q_points);
4263 * std::vector<double> wqn(n_q_points);
4264 * std::vector<double> wqnm1(n_q_points);
4268 * FOR Explicit nonlinearity
4269 * std::vector<Tensor<1, dim> > grad_un(n_q_points);
4270 * std::vector<Tensor<1, dim> > grad_vn(n_q_points);
4271 * std::vector<Tensor<1, dim> > grad_wn(n_q_points);
4272 * Tensor<1, dim> Un;
4278 * std::vector<Tensor<1, dim> > grad_pqn(n_q_points);
4279 * std::vector<Tensor<1, dim> > grad_psiqn(n_q_points);
4280 * std::vector<Tensor<1, dim> > grad_psiqnm1(n_q_points);
4282 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
4283 * std::vector<Tensor<1, dim> > shape_grad(dofs_per_cell);
4284 * std::vector<double> shape_value(dofs_per_cell);
4289 * double pressure_grad_u;
4290 * double pressure_grad_v;
4291 * double pressure_grad_w;
4297 * Vector<double> force_terms(dim);
4299 * typename DoFHandler<dim>::active_cell_iterator
4300 * cell_U=dof_handler_U.begin_active(), endc_U=dof_handler_U.end();
4301 * typename DoFHandler<dim>::active_cell_iterator cell_P=dof_handler_P.begin_active();
4302 * typename DoFHandler<dim>::active_cell_iterator cell_LS=dof_handler_LS.begin_active();
4304 * for (; cell_U!=endc_U; ++cell_U,++cell_P,++cell_LS)
4305 * if (cell_U->is_locally_owned())
4312 * fe_values_LS.reinit(cell_LS);
4313 * fe_values_U.reinit(cell_U);
4314 * fe_values_P.reinit(cell_P);
4318 * get function values for LS
4321 * fe_values_LS.get_function_values(locally_relevant_solution_phi,phiqnp1);
4324 * get function values for U
4327 * fe_values_U.get_function_values(locally_relevant_solution_u,uqn);
4328 * fe_values_U.get_function_values(locally_relevant_solution_u_old,uqnm1);
4329 * fe_values_U.get_function_values(locally_relevant_solution_v,vqn);
4330 * fe_values_U.get_function_values(locally_relevant_solution_v_old,vqnm1);
4333 * fe_values_U.get_function_values(locally_relevant_solution_w,wqn);
4334 * fe_values_U.get_function_values(locally_relevant_solution_w_old,wqnm1);
4338 * For explicit nonlinearity
4339 * get gradient values for U
4340 * fe_values_U.get_function_gradients(locally_relevant_solution_u,grad_un);
4341 * fe_values_U.get_function_gradients(locally_relevant_solution_v,grad_vn);
4343 * fe_values_U.get_function_gradients(locally_relevant_solution_w,grad_wn);
4347 * get values and gradients for p and dpsi
4350 * fe_values_P.get_function_gradients(locally_relevant_solution_p,grad_pqn);
4351 * fe_values_P.get_function_gradients(locally_relevant_solution_psi,grad_psiqn);
4352 * fe_values_P.get_function_gradients(locally_relevant_solution_psi_old,grad_psiqnm1);
4354 * for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
4356 * const double JxW=fe_values_U.JxW(q_point);
4357 * for (unsigned int i=0; i<dofs_per_cell; ++i)
4359 * shape_grad[i]=fe_values_U.shape_grad(i,q_point);
4360 * shape_value[i]=fe_values_U.shape_value(i,q_point);
4363 * pressure_grad_u=(grad_pqn[q_point][0]+4./3*grad_psiqn[q_point][0]-1./3*grad_psiqnm1[q_point][0]);
4364 * pressure_grad_v=(grad_pqn[q_point][1]+4./3*grad_psiqn[q_point][1]-1./3*grad_psiqnm1[q_point][1]);
4366 * pressure_grad_w=(grad_pqn[q_point][2]+4./3*grad_psiqn[q_point][2]-1./3*grad_psiqnm1[q_point][2]);
4368 * if (LEVEL_SET==1) // use level set to define rho and nu
4369 * get_rho_and_nu(phiqnp1[q_point]);
4370 * else // rho and nu are defined through functions
4372 * rho_value=rho_function.value(fe_values_U.quadrature_point(q_point));
4373 * nu_value=nu_function.value(fe_values_U.quadrature_point(q_point));
4378 * Non-linearity: for semi-implicit
4381 * u_star=2*uqn[q_point]-uqnm1[q_point];
4382 * v_star=2*vqn[q_point]-vqnm1[q_point];
4384 * w_star=2*wqn[q_point]-wqnm1[q_point];
4388 * for explicit nonlinearity
4389 * Un[0] = uqn[q_point];
4390 * Un[1] = vqn[q_point];
4392 * Un[2] = wqn[q_point];
4396 * double nonlinearity_u = Un*grad_un[q_point];
4397 * double nonlinearity_v = Un*grad_vn[q_point];
4398 * double nonlinearity_w = 0;
4400 * nonlinearity_w = Un*grad_wn[q_point];
4406 * rho_star=rho_value; // This is because we consider rho*u_t instead of (rho*u)_t
4414 * force_function.vector_value(fe_values_U.quadrature_point(q_point),force_terms);
4415 * force_u=force_terms[0];
4416 * force_v=force_terms[1];
4418 * force_w=force_terms[2];
4419 * if (RHO_TIMES_RHS==1)
4427 * for (unsigned int i=0; i<dofs_per_cell; ++i)
4429 * cell_rhs_u(i)+=((4./3*rho*uqn[q_point]-1./3*rho*uqnm1[q_point]
4430 * +2./3*time_step*(force_u-pressure_grad_u)
4433 * -2./3*time_step*rho*nonlinearity_u
4436 * )*shape_value[i])*JxW;
4437 * cell_rhs_v(i)+=((4./3*rho*vqn[q_point]-1./3*rho*vqnm1[q_point]
4438 * +2./3*time_step*(force_v-pressure_grad_v)
4441 * -2./3*time_step*rho*nonlinearity_v
4444 * )*shape_value[i])*JxW;
4446 * cell_rhs_w(i)+=((4./3*rho*wqn[q_point]-1./3*rho*wqnm1[q_point]
4447 * +2./3*time_step*(force_w-pressure_grad_w)
4450 * -2./3*time_step*rho*nonlinearity_w
4453 * )*shape_value[i])*JxW;
4454 * if (rebuild_Matrix_U==true)
4455 * for (unsigned int j=0; j<dofs_per_cell; ++j)
4458 * cell_A_u(i,j)+=(rho_star*shape_value[i]*shape_value[j]
4459 * +2./3*time_step*nu_value*(shape_grad[i]*shape_grad[j])
4460 * +2./3*time_step*rho*shape_value[i]
4461 * *(u_star*shape_grad[j][0]+v_star*shape_grad[j][1]) // semi-implicit NL
4464 * cell_A_u(i,j)+=(rho_star*shape_value[i]*shape_value[j]
4465 * +2./3*time_step*nu_value*(shape_grad[i]*shape_grad[j])
4466 * +2./3*time_step*rho*shape_value[i]
4467 * *(u_star*shape_grad[j][0]+v_star*shape_grad[j][1]+w_star*shape_grad[j][2]) // semi-implicit NL
4472 * cell_U->get_dof_indices(local_dof_indices);
4478 * if (rebuild_Matrix_U==true)
4479 * constraints.distribute_local_to_global(cell_A_u,local_dof_indices,system_Matrix_u);
4480 * constraints.distribute_local_to_global(cell_rhs_u,local_dof_indices,system_rhs_u);
4481 * constraints.distribute_local_to_global(cell_rhs_v,local_dof_indices,system_rhs_v);
4483 * constraints.distribute_local_to_global(cell_rhs_w,local_dof_indices,system_rhs_w);
4485 * system_rhs_u.compress(VectorOperation::add);
4486 * system_rhs_v.compress(VectorOperation::add);
4487 * if (dim==3) system_rhs_w.compress(VectorOperation::add);
4488 * if (rebuild_Matrix_U==true)
4490 * system_Matrix_u.compress(VectorOperation::add);
4491 * system_Matrix_v.copy_from(system_Matrix_u);
4493 * system_Matrix_w.copy_from(system_Matrix_u);
4497 * BOUNDARY CONDITIONS
4500 * system_rhs_u.set(boundary_values_id_u,boundary_values_u);
4501 * system_rhs_u.compress(VectorOperation::insert);
4502 * system_rhs_v.set(boundary_values_id_v,boundary_values_v);
4503 * system_rhs_v.compress(VectorOperation::insert);
4506 * system_rhs_w.set(boundary_values_id_w,boundary_values_w);
4507 * system_rhs_w.compress(VectorOperation::insert);
4509 * if (rebuild_Matrix_U)
4511 * system_Matrix_u.clear_rows(boundary_values_id_u,1);
4512 * system_Matrix_v.clear_rows(boundary_values_id_v,1);
4514 * system_Matrix_w.clear_rows(boundary_values_id_w,1);
4515 * if (rebuild_Matrix_U_preconditioners)
4522 * rebuild_Matrix_U_preconditioners=false;
4523 * preconditioner_Matrix_u.reset(new PETScWrappers::PreconditionBoomerAMG
4524 * (system_Matrix_u,PETScWrappers::PreconditionBoomerAMG::AdditionalData(false)));
4525 * preconditioner_Matrix_v.reset( new PETScWrappers::PreconditionBoomerAMG
4526 * (system_Matrix_v,PETScWrappers::PreconditionBoomerAMG::AdditionalData(false)));
4528 * preconditioner_Matrix_w.reset(new PETScWrappers::PreconditionBoomerAMG
4529 * (system_Matrix_w,PETScWrappers::PreconditionBoomerAMG::AdditionalData(false)));
4532 * rebuild_Matrix_U=true;
4533 * rebuild_Matrix_U_preconditioners=true;
4537 * void NavierStokesSolver<dim>::assemble_system_dpsi_q()
4539 * if (rebuild_S_M==true)
4547 * const QGauss<dim> quadrature_formula(degree_MAX+1);
4549 * FEValues<dim> fe_values_U(fe_U,quadrature_formula,
4550 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4551 * FEValues<dim> fe_values_P(fe_P,quadrature_formula,
4552 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4553 * FEValues<dim> fe_values_LS(fe_LS,quadrature_formula,
4554 * update_values|update_gradients|update_quadrature_points|update_JxW_values);
4556 * const unsigned int dofs_per_cell=fe_P.dofs_per_cell;
4557 * const unsigned int n_q_points=quadrature_formula.size();
4559 * FullMatrix<double> cell_S(dofs_per_cell,dofs_per_cell);
4560 * FullMatrix<double> cell_M(dofs_per_cell,dofs_per_cell);
4561 * Vector<double> cell_rhs_psi(dofs_per_cell);
4562 * Vector<double> cell_rhs_q(dofs_per_cell);
4564 * std::vector<double> phiqnp1(n_q_points);
4565 * std::vector<Tensor<1, dim> > gunp1(n_q_points);
4566 * std::vector<Tensor<1, dim> > gvnp1(n_q_points);
4567 * std::vector<Tensor<1, dim> > gwnp1(n_q_points);
4569 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
4570 * std::vector<double> shape_value(dofs_per_cell);
4571 * std::vector<Tensor<1, dim> > shape_grad(dofs_per_cell);
4573 * typename DoFHandler<dim>::active_cell_iterator
4574 * cell_P=dof_handler_P.begin_active(), endc_P=dof_handler_P.end();
4575 * typename DoFHandler<dim>::active_cell_iterator cell_U=dof_handler_U.begin_active();
4576 * typename DoFHandler<dim>::active_cell_iterator cell_LS=dof_handler_LS.begin_active();
4578 * for (; cell_P!=endc_P; ++cell_P,++cell_U,++cell_LS)
4579 * if (cell_P->is_locally_owned())
4586 * fe_values_P.reinit(cell_P);
4587 * fe_values_U.reinit(cell_U);
4588 * fe_values_LS.reinit(cell_LS);
4592 * get function values for LS
4595 * fe_values_LS.get_function_values(locally_relevant_solution_phi,phiqnp1);
4599 * get function grads for u and v
4602 * fe_values_U.get_function_gradients(locally_relevant_solution_u,gunp1);
4603 * fe_values_U.get_function_gradients(locally_relevant_solution_v,gvnp1);
4605 * fe_values_U.get_function_gradients(locally_relevant_solution_w,gwnp1);
4607 * for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
4609 * const double JxW=fe_values_P.JxW(q_point);
4610 * double divU = gunp1[q_point][0]+gvnp1[q_point][1];
4611 * if (dim==3) divU += gwnp1[q_point][2];
4612 * for (unsigned int i=0; i<dofs_per_cell; ++i)
4614 * shape_value[i]=fe_values_P.shape_value(i,q_point);
4615 * shape_grad[i]=fe_values_P.shape_grad(i,q_point);
4617 * if (LEVEL_SET==1) // use level set to define rho and nu
4618 * get_rho_and_nu (phiqnp1[q_point]);
4619 * else // rho and nu are defined through functions
4620 * nu_value=nu_function.value(fe_values_U.quadrature_point(q_point));
4622 * for (unsigned int i=0; i<dofs_per_cell; ++i)
4624 * cell_rhs_psi(i)+=-3./2./time_step*rho_min*divU*shape_value[i]*JxW;
4625 * cell_rhs_q(i)-=nu_value*divU*shape_value[i]*JxW;
4626 * if (rebuild_S_M==true)
4628 * for (unsigned int j=0; j<dofs_per_cell; ++j)
4631 * cell_S(i,j)+=shape_grad[i]*shape_grad[j]*JxW+1E-10;
4632 * cell_M(i,j)+=shape_value[i]*shape_value[j]*JxW;
4636 * cell_S(i,j)+=shape_grad[i]*shape_grad[j]*JxW;
4637 * cell_M(i,j)+=shape_value[i]*shape_value[j]*JxW;
4642 * cell_P->get_dof_indices(local_dof_indices);
4648 * if (rebuild_S_M==true)
4650 * constraints_psi.distribute_local_to_global(cell_S,local_dof_indices,system_S);
4651 * constraints_psi.distribute_local_to_global(cell_M,local_dof_indices,system_M);
4653 * constraints_psi.distribute_local_to_global(cell_rhs_q,local_dof_indices,system_rhs_q);
4654 * constraints_psi.distribute_local_to_global(cell_rhs_psi,local_dof_indices,system_rhs_psi);
4656 * if (rebuild_S_M==true)
4658 * system_M.compress(VectorOperation::add);
4659 * system_S.compress(VectorOperation::add);
4660 * if (rebuild_S_M_preconditioners)
4662 * rebuild_S_M_preconditioners=false;
4663 * preconditioner_S.reset(new PETScWrappers::PreconditionBoomerAMG
4664 * (system_S,PETScWrappers::PreconditionBoomerAMG::AdditionalData(true)));
4665 * preconditioner_M.reset(new PETScWrappers::PreconditionBoomerAMG
4666 * (system_M,PETScWrappers::PreconditionBoomerAMG::AdditionalData(true)));
4669 * system_rhs_psi.compress(VectorOperation::add);
4670 * system_rhs_q.compress(VectorOperation::add);
4671 * rebuild_S_M=false;
4676 * ///////////////////////////////////////////////////
4677 * ///////////////////// SOLVERS /////////////////////
4678 * ///////////////////////////////////////////////////
4682 * void NavierStokesSolver<dim>::solve_U(const AffineConstraints<double> &constraints,
4683 * PETScWrappers::MPI::SparseMatrix &Matrix,
4684 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
4685 * PETScWrappers::MPI::Vector &completely_distributed_solution,
4686 * const PETScWrappers::MPI::Vector &rhs)
4688 * SolverControl solver_control(dof_handler_U.n_dofs(),1e-6);
4691 * PETScWrappers::SolverCG solver(solver_control);
4692 * PETScWrappers::SolverGMRES solver(solver_control);
4693 * PETScWrappers::SolverChebychev solver(solver_control);
4696 * PETScWrappers::SolverBicgstab solver(solver_control);
4697 * constraints.distribute(completely_distributed_solution);
4698 * solver.solve(Matrix,completely_distributed_solution,rhs,*preconditioner);
4699 * constraints.distribute(completely_distributed_solution);
4700 * if (solver_control.last_step() > MAX_NUM_ITER_TO_RECOMPUTE_PRECONDITIONER)
4701 * rebuild_Matrix_U_preconditioners=true;
4702 * if (verbose==true)
4703 * pcout<<" Solved U in "<<solver_control.last_step()<<" iterations."<<std::endl;
4707 * void NavierStokesSolver<dim>::solve_P(const AffineConstraints<double> &constraints,
4708 * PETScWrappers::MPI::SparseMatrix &Matrix,
4709 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner,
4710 * PETScWrappers::MPI::Vector &completely_distributed_solution,
4711 * const PETScWrappers::MPI::Vector &rhs)
4713 * SolverControl solver_control(dof_handler_P.n_dofs(),1e-6);
4714 * PETScWrappers::SolverCG solver(solver_control);
4717 * PETScWrappers::SolverGMRES solver(solver_control);
4720 * constraints.distribute(completely_distributed_solution);
4721 * solver.solve(Matrix,completely_distributed_solution,rhs,*preconditioner);
4722 * constraints.distribute(completely_distributed_solution);
4723 * if (solver_control.last_step() > MAX_NUM_ITER_TO_RECOMPUTE_PRECONDITIONER)
4724 * rebuild_S_M_preconditioners=true;
4725 * if (verbose==true)
4726 * pcout<<" Solved P in "<<solver_control.last_step()<<" iterations."<<std::endl;
4731 * ///////////////////////////////////////////////////
4732 * ////////////// get different fields ///////////////
4733 * ///////////////////////////////////////////////////
4737 * void NavierStokesSolver<dim>::get_velocity()
4739 * assemble_system_U();
4740 * save_old_solution();
4741 * solve_U(constraints,system_Matrix_u,preconditioner_Matrix_u,completely_distributed_solution_u,system_rhs_u);
4742 * locally_relevant_solution_u=completely_distributed_solution_u;
4743 * solve_U(constraints,system_Matrix_v,preconditioner_Matrix_v,completely_distributed_solution_v,system_rhs_v);
4744 * locally_relevant_solution_v=completely_distributed_solution_v;
4747 * solve_U(constraints,system_Matrix_w,preconditioner_Matrix_w,completely_distributed_solution_w,system_rhs_w);
4748 * locally_relevant_solution_w=completely_distributed_solution_w;
4753 * void NavierStokesSolver<dim>::get_pressure()
4760 * assemble_system_dpsi_q();
4761 * solve_P(constraints_psi,system_S,preconditioner_S,completely_distributed_solution_psi,system_rhs_psi);
4762 * locally_relevant_solution_psi=completely_distributed_solution_psi;
4768 * solve_P(constraints,system_M,preconditioner_M,completely_distributed_solution_q,system_rhs_q);
4771 * UPDATE THE PRESSURE
4774 * completely_distributed_solution_p.add(1,completely_distributed_solution_psi);
4775 * completely_distributed_solution_p.add(1,completely_distributed_solution_q);
4776 * locally_relevant_solution_p = completely_distributed_solution_p;
4781 * ///////////////////////////////////////////////////
4782 * ///////////////////// DO STEPS ////////////////////
4783 * ///////////////////////////////////////////////////
4787 * void NavierStokesSolver<dim>::nth_time_step()
4795 * ///////////////////////////////////////////////////
4796 * ////////////////////// OTHERS /////////////////////
4797 * ///////////////////////////////////////////////////
4801 * void NavierStokesSolver<dim>::save_old_solution()
4803 * locally_relevant_solution_u_old=locally_relevant_solution_u;
4804 * locally_relevant_solution_v_old=locally_relevant_solution_v;
4805 * locally_relevant_solution_w_old=locally_relevant_solution_w;
4806 * locally_relevant_solution_psi_old=locally_relevant_solution_psi;
4812<a name="ann-TestLevelSet.cc"></a>
4813<h1>Annotated version of TestLevelSet.cc</h1>
4819 * /* -----------------------------------------------------------------------------
4821 * * SPDX-License-Identifier: LGPL-2.1-or-later
4822 * * Copyright (C) 2016 Manuel Quezada de Luna
4824 * * This file is part of the deal.II code gallery.
4826 * * -----------------------------------------------------------------------------
4829 * #include <deal.II/base/quadrature_lib.h>
4830 * #include <deal.II/base/function.h>
4831 * #include <deal.II/lac/affine_constraints.h>
4832 * #include <deal.II/lac/vector.h>
4833 * #include <deal.II/lac/full_matrix.h>
4834 * #include <deal.II/lac/solver_cg.h>
4835 * #include <deal.II/lac/petsc_sparse_matrix.h>
4836 * #include <deal.II/lac/petsc_vector.h>
4837 * #include <deal.II/lac/petsc_solver.h>
4838 * #include <deal.II/lac/petsc_precondition.h>
4839 * #include <deal.II/grid/grid_generator.h>
4840 * #include <deal.II/grid/tria_accessor.h>
4841 * #include <deal.II/grid/tria_iterator.h>
4842 * #include <deal.II/dofs/dof_handler.h>
4843 * #include <deal.II/dofs/dof_accessor.h>
4844 * #include <deal.II/dofs/dof_tools.h>
4845 * #include <deal.II/fe/fe_values.h>
4846 * #include <deal.II/fe/fe_q.h>
4847 * #include <deal.II/numerics/vector_tools.h>
4848 * #include <deal.II/numerics/data_out.h>
4849 * #include <deal.II/numerics/error_estimator.h>
4850 * #include <deal.II/base/utilities.h>
4851 * #include <deal.II/base/conditional_ostream.h>
4852 * #include <deal.II/base/index_set.h>
4853 * #include <deal.II/lac/sparsity_tools.h>
4854 * #include <deal.II/distributed/tria.h>
4855 * #include <deal.II/distributed/grid_refinement.h>
4856 * #include <deal.II/lac/petsc_vector.h>
4857 * #include <deal.II/base/convergence_table.h>
4858 * #include <deal.II/base/timer.h>
4859 * #include <deal.II/base/parameter_handler.h>
4860 * #include <deal.II/grid/grid_tools.h>
4861 * #include <deal.II/fe/mapping_q.h>
4862 * #include <deal.II/fe/fe_system.h>
4864 * #include <fstream>
4865 * #include <iostream>
4868 * using namespace dealii;
4872 * ///////////////////////
4873 * FOR TRANSPORT PROBLEM
4874 * ///////////////////////
4878 * #define FORWARD_EULER 0
4885 * #define CIRCULAR_ROTATION 0
4886 * #define DIAGONAL_ADVECTION 1
4892 * #define VARIABLE_VELOCITY 0
4894 * #include "utilities_test_LS.cc"
4895 * #include "LevelSetSolver.cc"
4899 * ///////////////////////////////////////////////////
4900 * /////////////////// MAIN CLASS ////////////////////
4901 * ///////////////////////////////////////////////////
4904 * template <int dim>
4905 * class TestLevelSet
4908 * TestLevelSet (const unsigned int degree_LS,
4909 * const unsigned int degree_U);
4919 * void set_boundary_inlet();
4920 * void get_boundary_values_phi(std::vector<unsigned int> &boundary_values_id_phi,
4921 * std::vector<double> &boundary_values_phi);
4927 * void get_interpolated_velocity();
4930 * SETUP AND INIT CONDITIONS
4934 * void initial_condition();
4935 * void init_constraints();
4941 * void process_solution(parallel::distributed::Triangulation<dim> &triangulation,
4942 * DoFHandler<dim> &dof_handler_LS,
4943 * PETScWrappers::MPI::Vector &solution);
4944 * void output_results();
4945 * void output_solution();
4952 * PETScWrappers::MPI::Vector locally_relevant_solution_phi;
4953 * PETScWrappers::MPI::Vector locally_relevant_solution_u;
4954 * PETScWrappers::MPI::Vector locally_relevant_solution_v;
4955 * PETScWrappers::MPI::Vector locally_relevant_solution_w;
4956 * PETScWrappers::MPI::Vector completely_distributed_solution_phi;
4957 * PETScWrappers::MPI::Vector completely_distributed_solution_u;
4958 * PETScWrappers::MPI::Vector completely_distributed_solution_v;
4959 * PETScWrappers::MPI::Vector completely_distributed_solution_w;
4965 * std::vector<unsigned int> boundary_values_id_phi;
4966 * std::vector<double> boundary_values_phi;
4973 * MPI_Comm mpi_communicator;
4974 * parallel::distributed::Triangulation<dim> triangulation;
4978 * DoFHandler<dim> dof_handler_LS;
4980 * IndexSet locally_owned_dofs_LS;
4981 * IndexSet locally_relevant_dofs_LS;
4984 * DoFHandler<dim> dof_handler_U;
4986 * IndexSet locally_owned_dofs_U;
4987 * IndexSet locally_relevant_dofs_U;
4989 * DoFHandler<dim> dof_handler_U_disp_field;
4990 * FESystem<dim> fe_U_disp_field;
4991 * IndexSet locally_owned_dofs_U_disp_field;
4992 * IndexSet locally_relevant_dofs_U_disp_field;
4994 * AffineConstraints<double> constraints;
4995 * AffineConstraints<double> constraints_disp_field;
4999 * double final_time;
5000 * unsigned int timestep_number;
5005 * int sharpness_integer;
5007 * unsigned int n_refinement;
5008 * unsigned int output_number;
5009 * double output_time;
5013 * ConditionalOStream pcout;
5020 * double cK; //compression coeff
5021 * double cE; //entropy-visc coeff
5022 * unsigned int TRANSPORT_TIME_INTEGRATION;
5023 * std::string ALGORITHM;
5024 * unsigned int PROBLEM;
5028 * FOR RECONSTRUCTION OF MATERIAL FIELDS
5031 * double eps, rho_air, rho_fluid;
5038 * PETScWrappers::MPI::SparseMatrix matrix_MC, matrix_MC_tnm1;
5039 * std::shared_ptr<PETScWrappers::PreconditionBoomerAMG> preconditioner_MC;
5043 * template <int dim>
5044 * TestLevelSet<dim>::TestLevelSet (const unsigned int degree_LS,
5045 * const unsigned int degree_U)
5047 * mpi_communicator (MPI_COMM_WORLD),
5048 * triangulation (mpi_communicator,
5049 * typename Triangulation<dim>::MeshSmoothing
5050 * (Triangulation<dim>::smoothing_on_refinement |
5051 * Triangulation<dim>::smoothing_on_coarsening)),
5052 * degree_LS(degree_LS),
5053 * dof_handler_LS (triangulation),
5054 * fe_LS (degree_LS),
5055 * degree_U(degree_U),
5056 * dof_handler_U (triangulation),
5058 * dof_handler_U_disp_field(triangulation),
5059 * fe_U_disp_field(FE_Q<dim>(degree_U),dim),
5060 * pcout (std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator)== 0))
5063 * template <int dim>
5064 * TestLevelSet<dim>::~TestLevelSet ()
5066 * dof_handler_U_disp_field.clear();
5067 * dof_handler_LS.clear ();
5068 * dof_handler_U.clear ();
5077 * template <int dim>
5078 * void TestLevelSet<dim>::get_interpolated_velocity()
5085 * completely_distributed_solution_u = 0;
5086 * VectorTools::interpolate(dof_handler_U,
5087 * ExactU<dim>(PROBLEM,time),
5088 * completely_distributed_solution_u);
5089 * constraints.distribute (completely_distributed_solution_u);
5090 * locally_relevant_solution_u = completely_distributed_solution_u;
5096 * completely_distributed_solution_v = 0;
5097 * VectorTools::interpolate(dof_handler_U,
5098 * ExactV<dim>(PROBLEM,time),
5099 * completely_distributed_solution_v);
5100 * constraints.distribute (completely_distributed_solution_v);
5101 * locally_relevant_solution_v = completely_distributed_solution_v;
5104 * completely_distributed_solution_w = 0;
5105 * VectorTools::interpolate(dof_handler_U,
5106 * ExactW<dim>(PROBLEM,time),
5107 * completely_distributed_solution_w);
5108 * constraints.distribute (completely_distributed_solution_w);
5109 * locally_relevant_solution_w = completely_distributed_solution_w;
5120 * template <int dim>
5121 * void TestLevelSet<dim>::set_boundary_inlet()
5123 * const QGauss<dim-1> face_quadrature_formula(1); // center of the face
5124 * FEFaceValues<dim> fe_face_values (fe_U,face_quadrature_formula,
5125 * update_values | update_quadrature_points |
5126 * update_normal_vectors);
5127 * const unsigned int n_face_q_points = face_quadrature_formula.size();
5128 * std::vector<double> u_value (n_face_q_points);
5129 * std::vector<double> v_value (n_face_q_points);
5130 * std::vector<double> w_value (n_face_q_points);
5132 * typename DoFHandler<dim>::active_cell_iterator
5133 * cell_U = dof_handler_U.begin_active(),
5134 * endc_U = dof_handler_U.end();
5137 * for (; cell_U!=endc_U; ++cell_U)
5138 * if (cell_U->is_locally_owned())
5139 * for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
5140 * if (cell_U->face(face)->at_boundary())
5142 * fe_face_values.reinit(cell_U,face);
5143 * fe_face_values.get_function_values(locally_relevant_solution_u,u_value);
5144 * fe_face_values.get_function_values(locally_relevant_solution_v,v_value);
5146 * fe_face_values.get_function_values(locally_relevant_solution_w,w_value);
5151 * if (fe_face_values.normal_vector(0)*u < -1e-14)
5152 * cell_U->face(face)->set_boundary_id(10);
5156 * template <int dim>
5157 * void TestLevelSet<dim>::get_boundary_values_phi(std::vector<unsigned int> &boundary_values_id_phi,
5158 * std::vector<double> &boundary_values_phi)
5160 * std::map<unsigned int, double> map_boundary_values_phi;
5161 * unsigned int boundary_id=0;
5163 * set_boundary_inlet();
5164 * boundary_id=10; // inlet
5165 * VectorTools::interpolate_boundary_values (dof_handler_LS,
5166 * boundary_id,BoundaryPhi<dim>(),
5167 * map_boundary_values_phi);
5169 * boundary_values_id_phi.resize(map_boundary_values_phi.size());
5170 * boundary_values_phi.resize(map_boundary_values_phi.size());
5171 * std::map<unsigned int,double>::const_iterator boundary_value_phi = map_boundary_values_phi.begin();
5172 * for (int i=0; boundary_value_phi !=map_boundary_values_phi.end(); ++boundary_value_phi, ++i)
5174 * boundary_values_id_phi[i]=boundary_value_phi->first;
5175 * boundary_values_phi[i]=boundary_value_phi->second;
5181 * ///////////////////////////////
5182 * SETUP AND INITIAL CONDITIONS
5183 * //////////////////////////////
5186 * template <int dim>
5187 * void TestLevelSet<dim>::setup()
5189 * degree = std::max(degree_LS,degree_U);
5195 * dof_handler_LS.distribute_dofs (fe_LS);
5196 * locally_owned_dofs_LS = dof_handler_LS.locally_owned_dofs ();
5197 * locally_relevant_dofs_LS = DoFTools::extract_locally_relevant_dofs (dof_handler_LS);
5203 * dof_handler_U.distribute_dofs (fe_U);
5204 * locally_owned_dofs_U = dof_handler_U.locally_owned_dofs ();
5205 * locally_relevant_dofs_U = DoFTools::extract_locally_relevant_dofs (dof_handler_U);
5208 * setup system U for disp field
5211 * dof_handler_U_disp_field.distribute_dofs (fe_U_disp_field);
5212 * locally_owned_dofs_U_disp_field = dof_handler_U_disp_field.locally_owned_dofs ();
5213 * locally_relevant_dofs_U_disp_field = DoFTools::extract_locally_relevant_dofs (dof_handler_U_disp_field);
5216 * init vectors for phi
5219 * locally_relevant_solution_phi.reinit(locally_owned_dofs_LS,
5220 * locally_relevant_dofs_LS,
5221 * mpi_communicator);
5222 * locally_relevant_solution_phi = 0;
5223 * completely_distributed_solution_phi.reinit(mpi_communicator,
5224 * dof_handler_LS.n_dofs(),
5225 * dof_handler_LS.n_locally_owned_dofs());
5228 * init vectors for u
5231 * locally_relevant_solution_u.reinit(locally_owned_dofs_U,
5232 * locally_relevant_dofs_U,
5233 * mpi_communicator);
5234 * locally_relevant_solution_u = 0;
5235 * completely_distributed_solution_u.reinit(mpi_communicator,
5236 * dof_handler_U.n_dofs(),
5237 * dof_handler_U.n_locally_owned_dofs());
5240 * init vectors for v
5243 * locally_relevant_solution_v.reinit(locally_owned_dofs_U,
5244 * locally_relevant_dofs_U,
5245 * mpi_communicator);
5246 * locally_relevant_solution_v = 0;
5247 * completely_distributed_solution_v.reinit(mpi_communicator,
5248 * dof_handler_U.n_dofs(),
5249 * dof_handler_U.n_locally_owned_dofs());
5252 * init vectors for w
5255 * locally_relevant_solution_w.reinit(locally_owned_dofs_U,
5256 * locally_relevant_dofs_U,
5257 * mpi_communicator);
5258 * locally_relevant_solution_w = 0;
5259 * completely_distributed_solution_w.reinit(mpi_communicator,
5260 * dof_handler_U.n_dofs(),
5261 * dof_handler_U.n_locally_owned_dofs());
5262 * init_constraints();
5268 * DynamicSparsityPattern dsp (locally_relevant_dofs_LS);
5269 * DoFTools::make_sparsity_pattern (dof_handler_LS,dsp,constraints,false);
5271 * const std::vector<types::global_dof_index> n_locally_owned_dofs_per_processor =
5272 * Utilities::MPI::all_gather(mpi_communicator, dof_handler_LS.n_locally_owned_dofs());
5274 * SparsityTools::distribute_sparsity_pattern (dsp,
5275 * n_locally_owned_dofs_per_processor,
5277 * locally_relevant_dofs_LS);
5278 * matrix_MC.reinit (mpi_communicator,
5280 * n_locally_owned_dofs_per_processor,
5281 * n_locally_owned_dofs_per_processor,
5282 * Utilities::MPI::this_mpi_process(mpi_communicator));
5283 * matrix_MC_tnm1.reinit (mpi_communicator,
5285 * n_locally_owned_dofs_per_processor,
5286 * n_locally_owned_dofs_per_processor,
5287 * Utilities::MPI::this_mpi_process(mpi_communicator));
5290 * template <int dim>
5291 * void TestLevelSet<dim>::initial_condition()
5296 * Initial conditions
5297 * init condition for phi
5300 * completely_distributed_solution_phi = 0;
5301 * VectorTools::interpolate(dof_handler_LS,
5302 * InitialPhi<dim>(PROBLEM, sharpness),
5305 * Functions::ZeroFunction<dim>(),
5308 * completely_distributed_solution_phi);
5309 * constraints.distribute (completely_distributed_solution_phi);
5310 * locally_relevant_solution_phi = completely_distributed_solution_phi;
5313 * init condition for u=0
5316 * completely_distributed_solution_u = 0;
5317 * VectorTools::interpolate(dof_handler_U,
5318 * ExactU<dim>(PROBLEM,time),
5319 * completely_distributed_solution_u);
5320 * constraints.distribute (completely_distributed_solution_u);
5321 * locally_relevant_solution_u = completely_distributed_solution_u;
5324 * init condition for v
5327 * completely_distributed_solution_v = 0;
5328 * VectorTools::interpolate(dof_handler_U,
5329 * ExactV<dim>(PROBLEM,time),
5330 * completely_distributed_solution_v);
5331 * constraints.distribute (completely_distributed_solution_v);
5332 * locally_relevant_solution_v = completely_distributed_solution_v;
5335 * template <int dim>
5336 * void TestLevelSet<dim>::init_constraints()
5338 * constraints.clear ();
5339 * constraints.reinit (locally_owned_dofs_LS, locally_relevant_dofs_LS);
5340 * DoFTools::make_hanging_node_constraints (dof_handler_LS, constraints);
5341 * constraints.close ();
5342 * constraints_disp_field.clear ();
5343 * constraints_disp_field.reinit (locally_owned_dofs_LS, locally_relevant_dofs_LS);
5344 * DoFTools::make_hanging_node_constraints (dof_handler_LS, constraints_disp_field);
5345 * constraints_disp_field.close ();
5355 * template <int dim>
5356 * void TestLevelSet<dim>::process_solution(parallel::distributed::Triangulation<dim> &triangulation,
5357 * DoFHandler<dim> &dof_handler_LS,
5358 * PETScWrappers::MPI::Vector &solution)
5360 * Vector<double> difference_per_cell (triangulation.n_active_cells());
5366 * VectorTools::integrate_difference (dof_handler_LS,
5368 * InitialPhi<dim>(PROBLEM,sharpness),
5369 * difference_per_cell,
5370 * QGauss<dim>(degree_LS+3),
5371 * VectorTools::L1_norm);
5373 * double u_L1_error = difference_per_cell.l1_norm();
5374 * u_L1_error = std::sqrt(Utilities::MPI::sum(u_L1_error * u_L1_error, mpi_communicator));
5376 * VectorTools::integrate_difference (dof_handler_LS,
5378 * InitialPhi<dim>(PROBLEM,sharpness),
5379 * difference_per_cell,
5380 * QGauss<dim>(degree_LS+3),
5381 * VectorTools::L2_norm);
5382 * double u_L2_error = difference_per_cell.l2_norm();
5383 * u_L2_error = std::sqrt(Utilities::MPI::sum(u_L2_error * u_L2_error, mpi_communicator));
5385 * pcout << "L1 error: " << u_L1_error << std::endl;
5386 * pcout << "L2 error: " << u_L2_error << std::endl;
5390 * void TestLevelSet<dim>::output_results()
5392 * output_solution();
5396 * template <int dim>
5397 * void TestLevelSet<dim>::output_solution()
5399 * DataOut<dim> data_out;
5400 * data_out.attach_dof_handler(dof_handler_LS);
5401 * data_out.add_data_vector (locally_relevant_solution_phi, "phi");
5402 * data_out.build_patches();
5404 * const std::string filename = ("solution-" +
5405 * Utilities::int_to_string (output_number, 3) +
5407 * Utilities::int_to_string
5408 * (triangulation.locally_owned_subdomain(), 4));
5409 * std::ofstream output ((filename + ".vtu").c_str());
5410 * data_out.write_vtu (output);
5412 * if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
5414 * std::vector<std::string> filenames;
5415 * for (unsigned int i=0;
5416 * i<Utilities::MPI::n_mpi_processes(mpi_communicator);
5418 * filenames.push_back ("solution-" +
5419 * Utilities::int_to_string (output_number, 3) +
5421 * Utilities::int_to_string (i, 4) +
5424 * std::ofstream master_output ((filename + ".pvtu").c_str());
5425 * data_out.write_pvtu_record (master_output, filenames);
5429 * template <int dim>
5430 * void TestLevelSet<dim>::run()
5434 * ////////////////////
5435 * GENERAL PARAMETERS
5436 * ////////////////////
5441 * get_output = true;
5442 * output_number = 0;
5445 * output_time = 0.1;
5447 * PROBLEM=CIRCULAR_ROTATION;
5450 * PROBLEM=DIAGONAL_ADVECTION;
5454 * if (PROBLEM==CIRCULAR_ROTATION)
5455 * umax = std::sqrt(2)*numbers::PI;
5457 * umax = std::sqrt(2);
5461 * //////////////////////////////////
5462 * PARAMETERS FOR TRANSPORT PROBLEM
5463 * //////////////////////////////////
5466 * cK = 1.0; // compression constant
5467 * cE = 1.0; // entropy viscosity constant
5468 * sharpness_integer=1; //this will be multiplied by min_h
5471 * TRANSPORT_TIME_INTEGRATION=FORWARD_EULER;
5474 * TRANSPORT_TIME_INTEGRATION=SSP33;
5477 * ALGORITHM = "MPP_u1";
5480 * ALGORITHM = "NMPP_uH";
5483 * ALGORITHM = "MPP_uH";
5492 * if (PROBLEM==CIRCULAR_ROTATION || PROBLEM==DIAGONAL_ADVECTION)
5493 * GridGenerator::hyper_cube(triangulation);
5496 * GridGenerator::hyper_rectangle(triangulation, Point<dim>(0.0,0.0), Point<dim>(1.0,1.0), true);
5499 * triangulation.refine_global (n_refinement);
5512 * for Reconstruction of MATERIAL FIELDS
5515 * min_h = GridTools::minimal_cell_diameter(triangulation)/std::sqrt(dim)/degree;
5516 * eps=1*min_h; //For reconstruction of density in Navier Stokes
5517 * sharpness=sharpness_integer*min_h; //adjust value of sharpness (for init cond of phi)
5526 * time_step = cfl*min_h/umax;
5530 * //////////////////
5532 * //////////////////
5535 * LevelSetSolver<dim> level_set (degree_LS,degree_U,
5539 * TRANSPORT_TIME_INTEGRATION,
5541 * mpi_communicator);
5545 * ///////////////////
5547 * ///////////////////
5550 * initial_condition();
5553 * level_set.initial_condition(locally_relevant_solution_phi,
5554 * locally_relevant_solution_u,locally_relevant_solution_v);
5556 * level_set.initial_condition(locally_relevant_solution_phi,
5557 * locally_relevant_solution_u,locally_relevant_solution_v,locally_relevant_solution_w);
5561 * /////////////////////////////
5562 * BOUNDARY CONDITIONS FOR PHI
5563 * /////////////////////////////
5566 * get_boundary_values_phi(boundary_values_id_phi,boundary_values_phi);
5567 * level_set.set_boundary_conditions(boundary_values_id_phi,boundary_values_phi);
5571 * OUTPUT DATA REGARDING TIME STEPPING AND MESH
5574 * int dofs_LS = dof_handler_LS.n_dofs();
5575 * pcout << "Cfl: " << cfl << std::endl;
5576 * pcout << " Number of active cells: "
5577 * << triangulation.n_global_active_cells() << std::endl
5578 * << " Number of degrees of freedom: " << std::endl
5579 * << " LS: " << dofs_LS << std::endl;
5586 * timestep_number=0;
5588 * while (time<final_time)
5590 * timestep_number++;
5591 * if (time+time_step > final_time)
5593 * pcout << "FINAL TIME STEP... " << std::endl;
5594 * time_step = final_time-time;
5596 * pcout << "Time step " << timestep_number
5597 * << "\twith dt=" << time_step
5598 * << "\tat tn=" << time << std::endl;
5603 * GET VELOCITY // (NS or interpolate from a function) at current time tn
5607 * if (VARIABLE_VELOCITY)
5609 * get_interpolated_velocity();
5612 * SET VELOCITY TO LEVEL SET SOLVER
5615 * level_set.set_velocity(locally_relevant_solution_u,locally_relevant_solution_v);
5619 * ////////////////////////
5620 * GET LEVEL SET SOLUTION // (at tnp1)
5621 * ////////////////////////
5624 * level_set.nth_time_step();
5633 * time+=time_step; // time tnp1
5642 * if (get_output && time-(output_number)*output_time>=0)
5644 * level_set.get_unp1(locally_relevant_solution_phi);
5648 * pcout << "FINAL TIME T=" << time << std::endl;
5651 * int main(int argc, char *argv[])
5655 * using namespace dealii;
5656 * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
5657 * deallog.depth_console (0);
5659 * unsigned int degree = 1;
5660 * 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_owned_dofs_LS, 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 * deallog.depth_console (0);
6513 * unsigned int degree_LS = 1;
6514 * unsigned int degree_U = 2;
6515 * TestNavierStokes<2> test_navier_stokes(degree_LS, degree_U);
6516 * test_navier_stokes.run();
6520 * catch (std::exception &exc)
6522 * std::cerr << std::endl << std::endl
6523 * << "----------------------------------------------------"
6525 * std::cerr << "Exception on processing: " << std::endl
6526 * << exc.what() << std::endl
6527 * << "Aborting!" << std::endl
6528 * << "----------------------------------------------------"
6535 * std::cerr << std::endl << std::endl
6536 * << "----------------------------------------------------"
6538 * std::cerr << "Unknown exception!" << std::endl
6539 * << "Aborting!" << std::endl
6540 * << "----------------------------------------------------"
6550<a name="ann-clean.sh"></a>
6551<h1>Annotated version of clean.sh</h1>
6553rm -rf CMakeFiles CMakeCache.txt Makefile cmake_install.cmake *~
6554rm -f MultiPhase TestLevelSet TestNavierStokes
6561<a name="ann-utilities.cc"></a>
6562<h1>Annotated version of utilities.cc</h1>
6568 * /* -----------------------------------------------------------------------------
6570 * * SPDX-License-Identifier: LGPL-2.1-or-later
6571 * * Copyright (C) 2016 Manuel Quezada de Luna
6573 * * This file is part of the deal.II code gallery.
6575 * * -----------------------------------------------------------------------------
6580 * /////////////////////////////////////////////////
6581 * ////////////////// INITIAL PHI //////////////////
6582 * /////////////////////////////////////////////////
6585 * template <int dim>
6586 * class InitialPhi : public Function <dim>
6589 * InitialPhi (unsigned int PROBLEM, double sharpness=0.005) : Function<dim>(),
6590 * sharpness(sharpness),
6591 * PROBLEM(PROBLEM) {}
6592 * virtual double value (const Point<dim> &p, const unsigned int component=0) const override;
6594 * unsigned int PROBLEM;
6596 * template <int dim>
6597 * double InitialPhi<dim>::value (const Point<dim> &p,
6598 * const unsigned int) const
6602 * double pi=numbers::PI;
6604 * if (PROBLEM==FILLING_TANK)
6605 * return 0.5*(-std::tanh((y-0.3)/sharpness)*std::tanh((y-0.35)/sharpness)+1)
6606 * *(-std::tanh((x-0.02)/sharpness)+1)-1;
6607 * else if (PROBLEM==BREAKING_DAM)
6608 * return 0.5*(-std::tanh((x-0.35)/sharpness)*std::tanh((x-0.65)/sharpness)+1)
6609 * *(1-std::tanh((y-0.35)/sharpness))-1;
6610 * else if (PROBLEM==FALLING_DROP)
6615 * double r = std::sqrt(std::pow(x-x0,2)+std::pow(y-y0,2));
6616 * return 1-(std::tanh((r-r0)/sharpness)+std::tanh((y-0.3)/sharpness));
6618 * else if (PROBLEM==SMALL_WAVE_PERTURBATION)
6620 * double wave = 0.1*std::sin(pi*x)+0.25;
6621 * return -std::tanh((y-wave)/sharpness);
6625 * std::cout << "Error in type of PROBLEM" << std::endl;
6632 * ///////////////////////////////////////////////////
6633 * ////////////////// FORCE TERMS ///// //////////////
6634 * ///////////////////////////////////////////////////
6637 * template <int dim>
6638 * class ForceTerms : public Functions::ConstantFunction <dim>
6641 * ForceTerms (const std::vector<double> values) : Functions::ConstantFunction<dim>(values) {}
6646 * /////////////////////////////////////////////////
6647 * ////////////////// BOUNDARY PHI /////////////////
6648 * /////////////////////////////////////////////////
6651 * template <int dim>
6652 * class BoundaryPhi : public Functions::ConstantFunction <dim>
6655 * BoundaryPhi (const double value, const unsigned int n_components=1) : Functions::ConstantFunction<dim>(value,n_components) {}
6660 * //////////////////////////////////////////////////////
6661 * ////////////////// BOUNDARY VELOCITY /////////////////
6662 * //////////////////////////////////////////////////////
6665 * template <int dim>
6666 * class BoundaryU : public Function <dim>
6669 * BoundaryU (unsigned int PROBLEM, double t=0) : Function<dim>(), PROBLEM(PROBLEM) {this->set_time(t);}
6670 * virtual double value (const Point<dim> &p, const unsigned int component=0) const override;
6673 * template <int dim>
6674 * double BoundaryU<dim>::value (const Point<dim> &p, const unsigned int) const
6678 * //////////////////
6680 * //////////////////
6681 * boundary for filling the tank (inlet)
6687 * if (PROBLEM==FILLING_TANK)
6689 * if (x==0 && y>=0.3 && y<=0.35)
6696 * std::cout << "Error in PROBLEM definition" << std::endl;
6701 * template <int dim>
6702 * class BoundaryV : public Function <dim>
6705 * BoundaryV (unsigned int PROBLEM, double t=0) : Function<dim>(), PROBLEM(PROBLEM) {this->set_time(t);}
6706 * virtual double value (const Point<dim> &p, const unsigned int component=0) const override;
6707 * unsigned int PROBLEM;
6709 * template <int dim>
6710 * double BoundaryV<dim>::value (const Point<dim> &p, const unsigned int) const
6714 * boundary for filling the tank (outlet)
6719 * double return_value = 0;
6721 * if (PROBLEM==FILLING_TANK)
6723 * if (y==0.4 && x>=0.3 && x<=0.35)
6724 * return_value = 0.25;
6726 * return return_value;
6731 * ///////////////////////////////////////////////////
6732 * ///////////////// POST-PROCESSING /////////////////
6733 * ///////////////////////////////////////////////////
6736 * template <int dim>
6737 * class Postprocessor : public DataPostprocessorScalar <dim>
6740 * Postprocessor(double eps, double rho_air, double rho_fluid)
6742 * DataPostprocessorScalar<dim>("Density",update_values)
6745 * this->rho_air=rho_air;
6746 * this->rho_fluid=rho_fluid;
6751 * evaluate_scalar_field (const DataPostprocessorInputs::Scalar<dim> &input_data,
6752 * std::vector<Vector<double> > &computed_quantities) const override;
6760 * template <int dim>
6762 * Postprocessor<dim>::
6763 * evaluate_scalar_field (const DataPostprocessorInputs::Scalar<dim> &input_data,
6764 * std::vector<Vector<double> > &computed_quantities) const
6766 * const unsigned int n_quadrature_points = input_data.solution_values.size();
6767 * for (unsigned int q=0; q<n_quadrature_points; ++q)
6771 * double phi_value=input_data.solution_values[q];
6772 * if (phi_value > eps)
6774 * else if (phi_value < -eps)
6778 * rho_value = rho_fluid*(1+H)/2. + rho_air*(1-H)/2.;
6779 * computed_quantities[q] = rho_value;
6786<a name="ann-utilities_test_LS.cc"></a>
6787<h1>Annotated version of utilities_test_LS.cc</h1>
6793 * /* -----------------------------------------------------------------------------
6795 * * SPDX-License-Identifier: LGPL-2.1-or-later
6796 * * Copyright (C) 2016 Manuel Quezada de Luna
6798 * * This file is part of the deal.II code gallery.
6800 * * -----------------------------------------------------------------------------
6805 * ///////////////////////////////////////////////////
6806 * ////////////////// INITIAL PHI //////////////////
6807 * ///////////////////////////////////////////////////
6810 * template <int dim>
6811 * class InitialPhi : public Function <dim>
6814 * InitialPhi (unsigned int PROBLEM, double sharpness=0.005) : Function<dim>(),
6815 * sharpness(sharpness),
6816 * PROBLEM(PROBLEM) {}
6817 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6819 * unsigned int PROBLEM;
6821 * template <int dim>
6822 * double InitialPhi<dim>::value (const Point<dim> &p,
6823 * const unsigned int) const
6827 * double return_value = -1.;
6829 * if (PROBLEM==CIRCULAR_ROTATION)
6834 * double r = std::sqrt(std::pow(x-x0,2)+std::pow(y-y0,2));
6835 * return_value = -std::tanh((r-r0)/sharpness);
6837 * else // (PROBLEM==DIAGONAL_ADVECTION)
6844 * r = std::sqrt(std::pow(x-x0,2)+std::pow(y-y0,2));
6849 * r = std::sqrt(std::pow(x-x0,2)+std::pow(y-y0,2)+std::pow(z-z0,2));
6851 * return_value = -std::tanh((r-r0)/sharpness);
6853 * return return_value;
6858 * /////////////////////////////////////////////////
6859 * ////////////////// BOUNDARY PHI /////////////////
6860 * /////////////////////////////////////////////////
6863 * template <int dim>
6864 * class BoundaryPhi : public Function <dim>
6867 * BoundaryPhi (double t=0)
6870 * {this->set_time(t);}
6871 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6874 * template <int dim>
6875 * double BoundaryPhi<dim>::value (const Point<dim> &, const unsigned int) const
6882 * ///////////////////////////////////////////////////
6883 * ////////////////// EXACT VELOCITY /////////////////
6884 * ///////////////////////////////////////////////////
6887 * template <int dim>
6888 * class ExactU : public Function <dim>
6891 * ExactU (unsigned int PROBLEM, double time=0) : Function<dim>(), PROBLEM(PROBLEM), time(time) {}
6892 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6893 * void set_time(double time) {this->time=time;};
6898 * template <int dim>
6899 * double ExactU<dim>::value (const Point<dim> &p, const unsigned int) const
6901 * if (PROBLEM==CIRCULAR_ROTATION)
6902 * return -2*numbers::PI*(p[1]-0.5);
6903 * else // (PROBLEM==DIAGONAL_ADVECTION)
6907 * template <int dim>
6908 * class ExactV : public Function <dim>
6911 * ExactV (unsigned int PROBLEM, double time=0) : Function<dim>(), PROBLEM(PROBLEM), time(time) {}
6912 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6913 * void set_time(double time) {this->time=time;};
6914 * unsigned int PROBLEM;
6918 * template <int dim>
6919 * double ExactV<dim>::value (const Point<dim> &p, const unsigned int) const
6921 * if (PROBLEM==CIRCULAR_ROTATION)
6922 * return 2*numbers::PI*(p[0]-0.5);
6923 * else // (PROBLEM==DIAGONAL_ADVECTION)
6927 * template <int dim>
6928 * class ExactW : public Function <dim>
6931 * ExactW (unsigned int PROBLEM, double time=0) : Function<dim>(), PROBLEM(PROBLEM), time(time) {}
6932 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6933 * void set_time(double time) {this->time=time;};
6934 * unsigned int PROBLEM;
6938 * template <int dim>
6939 * double ExactW<dim>::value (const Point<dim> &, const unsigned int) const
6943 * PROBLEM = 3D_DIAGONAL_ADVECTION
6952<a name="ann-utilities_test_NS.cc"></a>
6953<h1>Annotated version of utilities_test_NS.cc</h1>
6959 * /* -----------------------------------------------------------------------------
6961 * * SPDX-License-Identifier: LGPL-2.1-or-later
6962 * * Copyright (C) 2016 Manuel Quezada de Luna
6964 * * This file is part of the deal.II code gallery.
6966 * * -----------------------------------------------------------------------------
6971 * ///////////////////////////////////////////////////
6972 * ////////// EXACT SOLUTION RHO TO TEST NS //////////
6973 * ///////////////////////////////////////////////////
6976 * template <int dim>
6977 * class RhoFunction : public Function <dim>
6980 * RhoFunction (double t=0) : Function<dim>() {this->set_time(t);}
6981 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
6983 * template <int dim>
6984 * double RhoFunction<dim>::value (const Point<dim> &p,
6985 * const unsigned int) const
6987 * double t = this->get_time();
6988 * double return_value = 0;
6990 * return_value = std::pow(std::sin(p[0]+p[1]+t),2)+1;
6992 * return_value = std::pow(std::sin(p[0]+p[1]+p[2]+t),2)+1;
6993 * return return_value;
6996 * template <int dim>
6997 * class NuFunction : public Function <dim>
7000 * NuFunction (double t=0) : Function<dim>() {this->set_time(t);}
7001 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
7003 * template <int dim>
7004 * double NuFunction<dim>::value (const Point<dim> &, const unsigned int) const
7011 * //////////////////////////////////////////////////////////////
7012 * ///////////////// EXACT SOLUTION U to TEST NS ////////////////
7013 * //////////////////////////////////////////////////////////////
7016 * template <int dim>
7017 * class ExactSolution_and_BC_U : public Function <dim>
7020 * ExactSolution_and_BC_U (double t=0, int field=0)
7025 * this->set_time(t);
7027 * virtual double value (const Point<dim> &p, const unsigned int component=1) const;
7028 * virtual Tensor<1,dim> gradient (const Point<dim> &p, const unsigned int component=1) const;
7029 * virtual void set_field(int field) {this->field=field;}
7031 * unsigned int type_simulation;
7033 * template <int dim>
7034 * double ExactSolution_and_BC_U<dim>::value (const Point<dim> &p,
7035 * const unsigned int) const
7037 * double t = this->get_time();
7038 * double return_value = 0;
7039 * double Pi = numbers::PI;
7046 * return_value = std::sin(x)*std::sin(y+t);
7048 * return_value = std::cos(x)*std::cos(y+t);
7053 * return_value = std::cos(t)*std::cos(Pi*y)*std::cos(Pi*z)*std::sin(Pi*x);
7054 * else if (field == 1)
7055 * return_value = std::cos(t)*std::cos(Pi*x)*std::cos(Pi*z)*std::sin(Pi*y);
7057 * return_value = -2*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::sin(Pi*z);
7059 * return return_value;
7061 * template <int dim>
7062 * Tensor<1,dim> ExactSolution_and_BC_U<dim>::gradient (const Point<dim> &p,
7063 * const unsigned int) const
7067 * THIS IS USED JUST FOR TESTING NS
7070 * Tensor<1,dim> return_value;
7071 * double t = this->get_time();
7072 * double Pi = numbers::PI;
7079 * return_value[0] = std::cos(x)*std::sin(y+t);
7080 * return_value[1] = std::sin(x)*std::cos(y+t);
7084 * return_value[0] = -std::sin(x)*std::cos(y+t);
7085 * return_value[1] = -std::cos(x)*std::sin(y+t);
7092 * return_value[0] = Pi*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::cos(Pi*z);
7093 * return_value[1] = -(Pi*std::cos(t)*std::cos(Pi*z)*std::sin(Pi*x)*std::sin(Pi*y));
7094 * return_value[2] = -(Pi*std::cos(t)*std::cos(Pi*y)*std::sin(Pi*x)*std::sin(Pi*z));
7096 * else if (field == 1)
7098 * return_value[0] = -(Pi*std::cos(t)*std::cos(Pi*z)*std::sin(Pi*x)*std::sin(Pi*y));
7099 * return_value[1] = Pi*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::cos(Pi*z);
7100 * return_value[2] = -(Pi*std::cos(t)*std::cos(Pi*x)*std::sin(Pi*y)*std::sin(Pi*z));
7104 * return_value[0] = 2*Pi*std::cos(t)*std::cos(Pi*y)*std::sin(Pi*x)*std::sin(Pi*z);
7105 * return_value[1] = 2*Pi*std::cos(t)*std::cos(Pi*x)*std::sin(Pi*y)*std::sin(Pi*z);
7106 * return_value[2] = -2*Pi*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::cos(Pi*z);
7109 * return return_value;
7114 * ///////////////////////////////////////////////////
7115 * ///////// EXACT SOLUTION FOR p TO TEST NS /////////
7116 * ///////////////////////////////////////////////////
7119 * template <int dim>
7120 * class ExactSolution_p : public Function <dim>
7123 * ExactSolution_p (double t=0) : Function<dim>() {this->set_time(t);}
7124 * virtual double value (const Point<dim> &p, const unsigned int component=0) const;
7125 * virtual Tensor<1,dim> gradient (const Point<dim> &p, const unsigned int component = 0) const;
7128 * template <int dim>
7129 * double ExactSolution_p<dim>::value (const Point<dim> &p, const unsigned int) const
7131 * double t = this->get_time();
7132 * double return_value = 0;
7134 * return_value = std::cos(p[0])*std::sin(p[1]+t);
7136 * return_value = std::sin(p[0]+p[1]+p[2]+t);
7137 * return return_value;
7140 * template <int dim>
7141 * Tensor<1,dim> ExactSolution_p<dim>::gradient (const Point<dim> &p, const unsigned int) const
7143 * Tensor<1,dim> return_value;
7144 * double t = this->get_time();
7147 * return_value[0] = -std::sin(p[0])*std::sin(p[1]+t);
7148 * return_value[1] = std::cos(p[0])*std::cos(p[1]+t);
7152 * return_value[0] = std::cos(t+p[0]+p[1]+p[2]);
7153 * return_value[1] = std::cos(t+p[0]+p[1]+p[2]);
7154 * return_value[2] = std::cos(t+p[0]+p[1]+p[2]);
7156 * return return_value;
7161 * //////////////////////////////////////////////////////////////
7162 * ////////////////// FORCE TERMS to TEST NS ////////////////////
7163 * //////////////////////////////////////////////////////////////
7166 * template <int dim>
7167 * class ForceTerms : public Function <dim>
7170 * ForceTerms (double t=0)
7174 * this->set_time(t);
7177 * virtual void vector_value (const Point<dim> &p, Vector<double> &values) const;
7181 * template <int dim>
7182 * void ForceTerms<dim>::vector_value (const Point<dim> &p, Vector<double> &values) const
7187 * double t = this->get_time();
7188 * double Pi = numbers::PI;
7197 * values[0] = std::cos(t+y)*std::sin(x)*(1+std::pow(std::sin(t+x+y),2)) // time derivative
7198 * +2*nu*std::sin(x)*std::sin(t+y) // viscosity
7199 * +std::cos(x)*std::sin(x)*(1+std::pow(std::sin(t+x+y),2)) // non-linearity
7200 * -std::sin(x)*std::sin(y+t); // pressure
7206 * values[1] = -(std::cos(x)*std::sin(t+y)*(1+std::pow(std::sin(t+x+y),2))) // time derivative
7207 * +2*nu*std::cos(x)*std::cos(t+y) // viscosity
7208 * -(std::sin(2*(t+y))*(1+std::pow(std::sin(t+x+y),2)))/2. // non-linearity
7209 * +std::cos(x)*std::cos(y+t); // pressure
7220 * -(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.
7221 * +3*std::pow(Pi,2)*std::cos(t)*std::cos(Pi*y)*std::cos(Pi*z)*std::sin(Pi*x) //viscosity
7222 * -(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
7223 * +std::cos(t+x+y+z); // pressure
7225 * -(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
7226 * +3*std::pow(Pi,2)*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*z)*std::sin(Pi*y) //viscosity
7227 * -(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
7228 * +std::cos(t+x+y+z); // pressure
7230 * 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
7231 * -6*std::pow(Pi,2)*std::cos(t)*std::cos(Pi*x)*std::cos(Pi*y)*std::sin(Pi*z) //viscosity
7232 * -(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
7233 * +std::cos(t+x+y+z); // pressure
* * for(const auto &cell :triangulation.active_cell_iterators())
* * int main(int argc, char **argv)
* x_component_mask set(0, true)
* * * struct InterferenceTaperTransform *
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
void subtract_set(const IndexSet &other)
unsigned int depth_console(const unsigned int n)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
@ matrix
Contents is actually a matrix.
constexpr types::blas_int zero
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
* * * * ::Tensor< 1, dim, Number > TimeRateRequest< ValueType, dim, Number > get_velocity() const
* * * * void TimeRateRequest< ValueType, dim, Number > set_velocity(const ::Tensor< 1, dim, Number > &velocity)
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 save(Archive &ar, const ::std_cxx26::inplace_vector< T, N > &vec, const unsigned int)
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 > &)
1.9.8