179 * <a name=
"equation_data.h-Equationdata"></a>
193 *
static const unsigned int degree_p = 1;
208 *
Velocity(
const double initial_time = 0.0);
211 *
const unsigned int component = 0)
const override;
213 *
virtual void vector_value(
const Point<dim>& p,
225 *
if(component == 0) {
226 *
const double Um = 1.5;
227 *
const double H = 4.1;
229 *
return 4.0*
Um*p(1)*(
H - p(1))/(
H*
H);
240 *
for(
unsigned int i = 0; i < dim; ++i)
256 *
Pressure(
const double initial_time = 0.0);
259 *
const unsigned int component = 0)
const override;
272 *
return 22.0 - p(0);
279<a name=
"ann-navier_stokes_TRBDF2_DG.cc"></a>
365 *
#
include "runtime_parameters.h"
371 *
template<
int dim,
typename Number,
typename VectorizedArrayType>
376 *
const unsigned int&,
380 *
const unsigned int&,
384 *
const unsigned int&,
386 *
const unsigned int dof_no = 0) {
394 *
const unsigned int dummy = 0;
428 *
virtual std::vector<std::string> get_names()
const override;
430 *
virtual std::vector<DataComponentInterpretation::DataComponentInterpretation>
431 *
get_data_component_interpretation()
const override;
433 *
virtual UpdateFlags get_needed_update_flags()
const override;
447 *
const unsigned int n_quadrature_points =
inputs.solution_values.size();
464 *
for(
unsigned int q = 0;
q < n_quadrature_points; ++
q)
468 *
for(
unsigned int q = 0;
q < n_quadrature_points; ++
q) {
487 *
std::vector<std::string> names;
488 *
names.emplace_back(
"vorticity");
490 *
names.emplace_back(
"vorticity");
491 *
names.emplace_back(
"vorticity");
507 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
509 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
interpretation;
546 *
struct ScratchData {
551 *
ScratchData(
const ScratchData<dim>& scratch_data): fe_values(scratch_data.fe_values.get_fe(),
552 *
scratch_data.fe_values.get_quadrature(),
553 *
scratch_data.fe_values.get_update_flags()) {}
569 *
CopyData(
const CopyData &) =
default;
579 * <a name=
"navier_stokes_TRBDF2_DG.cc-codeNavierStokesProjectionOperatorNavierStokesProjectionOperatorcode"></a>
595 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
598 *
using Number =
typename Vec::value_type;
634 *
virtual void apply_add(Vec& dst,
const Vec& src)
const override;
639 *
const double a21 = 0.5;
640 *
const double a22 = 0.5;
656 *
const std::vector<Vec>& src,
657 *
const std::pair<unsigned int, unsigned int>&
cell_range)
const;
660 *
const std::vector<Vec>& src,
661 *
const std::pair<unsigned int, unsigned int>&
face_range)
const;
664 *
const std::vector<Vec>& src,
665 *
const std::pair<unsigned int, unsigned int>&
face_range)
const;
669 *
const std::vector<Vec>& src,
670 *
const std::pair<unsigned int, unsigned int>&
cell_range)
const;
673 *
const std::vector<Vec>& src,
674 *
const std::pair<unsigned int, unsigned int>&
face_range)
const;
677 *
const std::vector<Vec>& src,
678 *
const std::pair<unsigned int, unsigned int>&
face_range)
const;
683 *
const std::pair<unsigned int, unsigned int>&
cell_range)
const;
687 *
const std::pair<unsigned int, unsigned int>&
face_range)
const;
691 *
const std::pair<unsigned int, unsigned int>&
face_range)
const;
696 *
const std::pair<unsigned int, unsigned int>&
cell_range)
const;
700 *
const std::pair<unsigned int, unsigned int>&
face_range)
const;
704 *
const std::pair<unsigned int, unsigned int>&
face_range)
const;
709 *
const std::pair<unsigned int, unsigned int>&
cell_range)
const;
713 *
const std::pair<unsigned int, unsigned int>&
cell_range)
const;
717 *
const unsigned int& src,
718 *
const std::pair<unsigned int, unsigned int>&
cell_range)
const;
721 *
const unsigned int& src,
722 *
const std::pair<unsigned int, unsigned int>&
face_range)
const;
725 *
const unsigned int& src,
726 *
const std::pair<unsigned int, unsigned int>&
face_range)
const;
730 *
const unsigned int& src,
731 *
const std::pair<unsigned int, unsigned int>&
cell_range)
const;
734 *
const unsigned int& src,
735 *
const std::pair<unsigned int, unsigned int>&
face_range)
const;
738 *
const unsigned int& src,
739 *
const std::pair<unsigned int, unsigned int>&
face_range)
const;
752 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
767 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
784 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
800 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
819 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
837 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
841 *
u_extr.update_ghost_values();
855 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
859 *
const std::vector<Vec>& src,
860 *
const std::pair<unsigned int, unsigned int>&
cell_range)
const {
888 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
896 *
for(
unsigned int d = 0;
d < dim; ++
d)
899 *
phi.submit_value(1.0/(gamma*
dt)*
u_n,
q);
927 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
937 *
for(
unsigned int d = 0;
d < dim; ++
d)
958 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
962 *
const std::vector<Vec>& src,
963 *
const std::pair<unsigned int, unsigned int>&
face_range)
const {
991 *
phi_p.reinit(face);
992 *
phi_m.reinit(face);
995 *
for(
unsigned int q = 0;
q <
phi_p.n_q_points; ++
q) {
1037 *
phi_p.reinit(face);
1038 *
phi_m.reinit(face);
1041 *
for(
unsigned int q = 0;
q <
phi_p.n_q_points; ++
q) {
1072 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1076 *
const std::vector<Vec>& src,
1077 *
const std::pair<unsigned int, unsigned int>&
face_range)
const {
1101 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
1102 *
const auto&
n_plus =
phi.get_normal_vector(
q);
1109 *
if(boundary_id == 0) {
1114 *
for(
unsigned int d = 0;
d < dim; ++
d)
1116 *
for(
unsigned int d = 0;
d < dim; ++
d)
1157 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
1158 *
const auto&
n_plus =
phi.get_normal_vector(
q);
1167 *
if(boundary_id == 0) {
1170 *
for(
unsigned int d = 0;
d < dim; ++
d)
1172 *
for(
unsigned int d = 0;
d < dim; ++
d)
1199 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1202 *
for(
auto&
vec : src)
1203 *
vec.update_ghost_values();
1208 *
this, dst, src,
true,
1223 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1227 *
const std::vector<Vec>& src,
1228 *
const std::pair<unsigned int, unsigned int>&
cell_range)
const {
1248 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
1268 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1272 *
const std::vector<Vec>& src,
1273 *
const std::pair<unsigned int, unsigned int>&
face_range)
const {
1288 *
phi_p.reinit(face);
1289 *
phi_m.reinit(face);
1292 *
for(
unsigned int q = 0;
q <
phi_p.n_q_points; ++
q) {
1313 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1317 *
const std::vector<Vec>& src,
1318 *
const std::pair<unsigned int, unsigned int>&
face_range)
const {
1332 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
1333 *
const auto&
n_plus =
phi.get_normal_vector(
q);
1350 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1353 *
for(
auto&
vec : src)
1354 *
vec.update_ghost_values();
1359 *
this, dst, src,
true,
1374 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1379 *
const std::pair<unsigned int, unsigned int>&
cell_range)
const {
1394 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
1419 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
1442 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1447 *
const std::pair<unsigned int, unsigned int>&
face_range)
const {
1457 *
phi_p.reinit(face);
1459 *
phi_m.reinit(face);
1470 *
for(
unsigned int q = 0;
q <
phi_p.n_q_points; ++
q) {
1500 *
phi_p.reinit(face);
1502 *
phi_m.reinit(face);
1513 *
for(
unsigned int q = 0;
q <
phi_p.n_q_points; ++
q) {
1545 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1550 *
const std::pair<unsigned int, unsigned int>&
face_range)
const {
1568 *
if(boundary_id != 1) {
1572 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
1573 *
const auto&
n_plus =
phi.get_normal_vector(
q);
1587 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
1588 *
const auto&
n_plus =
phi.get_normal_vector(
q);
1598 *
for(
unsigned int d = 0;
d < dim; ++
d)
1631 *
if(boundary_id != 1) {
1635 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
1636 *
const auto&
n_plus =
phi.get_normal_vector(
q);
1638 *
const auto&
u =
phi.get_value(
q);
1650 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
1651 *
const auto&
n_plus =
phi.get_normal_vector(
q);
1653 *
const auto&
u =
phi.get_value(
q);
1661 *
for(
unsigned int d = 0;
d < dim; ++
d)
1664 *
u_m[1][v] = -
u_m[1][v];
1689 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1694 *
const std::pair<unsigned int, unsigned int>&
cell_range)
const {
1706 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
1707 *
phi.submit_gradient(
phi.get_gradient(
q),
q);
1724 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1729 *
const std::pair<unsigned int, unsigned int>&
face_range)
const {
1736 *
phi_p.reinit(face);
1738 *
phi_m.reinit(face);
1745 *
for(
unsigned int q = 0;
q <
phi_p.n_q_points; ++
q) {
1770 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1775 *
const std::pair<unsigned int, unsigned int>&
face_range)
const {
1781 *
if(boundary_id == 1) {
1787 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
1788 *
const auto&
n_plus =
phi.get_normal_vector(
q);
1791 *
const auto&
pres =
phi.get_value(
q);
1813 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1818 *
const std::pair<unsigned int, unsigned int>&
cell_range)
const {
1830 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q)
1846 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1850 *
this, dst, src,
true);
1862 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1867 *
const std::pair<unsigned int, unsigned int>&
cell_range)
const {
1876 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q)
1877 *
phi.submit_value(
phi.get_value(
q),
q);
1893 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1895 *
apply_add(Vec& dst,
const Vec& src)
const {
1900 *
this, dst, src,
false,
1908 *
this, dst, src,
false,
1914 *
this, dst, src,
false);
1932 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
1936 *
const unsigned int& ,
1937 *
const std::pair<unsigned int, unsigned int>&
cell_range)
const {
1946 *
for(
unsigned int d = 0;
d < dim; ++
d)
1956 *
for(
unsigned int i = 0; i <
phi.dofs_per_component; ++i) {
1957 *
for(
unsigned int j = 0;
j <
phi.dofs_per_component; ++
j)
1959 *
phi.submit_dof_value(tmp, i);
1963 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
1975 *
for(
unsigned int i = 0; i <
phi.dofs_per_component; ++i)
1976 *
phi.submit_dof_value(diagonal[i], i);
1977 *
phi.distribute_local_to_global(dst);
1986 *
for(
unsigned int d = 0;
d < dim; ++
d)
1996 *
for(
unsigned int i = 0; i <
phi.dofs_per_component; ++i) {
1997 *
for(
unsigned int j = 0;
j <
phi.dofs_per_component; ++
j)
1999 *
phi.submit_dof_value(tmp, i);
2003 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
2015 *
for(
unsigned int i = 0; i <
phi.dofs_per_component; ++i)
2016 *
phi.submit_dof_value(diagonal[i], i);
2017 *
phi.distribute_local_to_global(dst);
2031 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
2035 *
const unsigned int& ,
2036 *
const std::pair<unsigned int, unsigned int>&
face_range)
const {
2050 *
for(
unsigned int d = 0;
d < dim; ++
d)
2059 *
phi_p.reinit(face);
2060 *
phi_m.reinit(face);
2066 *
for(
unsigned int i = 0; i <
phi_p.dofs_per_component; ++i) {
2067 *
for(
unsigned int j = 0;
j <
phi_p.dofs_per_component; ++
j) {
2071 *
phi_p.submit_dof_value(tmp, i);
2073 *
phi_m.submit_dof_value(tmp, i);
2077 *
for(
unsigned int q = 0;
q <
phi_p.n_q_points; ++
q) {
2098 *
for(
unsigned int i = 0; i <
phi_p.dofs_per_component; ++i) {
2102 *
phi_p.distribute_local_to_global(dst);
2103 *
phi_m.distribute_local_to_global(dst);
2116 *
for(
unsigned int d = 0;
d < dim; ++
d)
2125 *
phi_p.reinit(face);
2126 *
phi_m.reinit(face);
2132 *
for(
unsigned int i = 0; i <
phi_p.dofs_per_component; ++i) {
2133 *
for(
unsigned int j = 0;
j <
phi_p.dofs_per_component; ++
j) {
2137 *
phi_p.submit_dof_value(tmp, i);
2139 *
phi_m.submit_dof_value(tmp, i);
2143 *
for(
unsigned int q = 0;
q <
phi_p.n_q_points; ++
q) {
2164 *
for(
unsigned int i = 0; i <
phi_p.dofs_per_component; ++i) {
2168 *
phi_p.distribute_local_to_global(dst);
2169 *
phi_m.distribute_local_to_global(dst);
2183 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
2187 *
const unsigned int& ,
2188 *
const std::pair<unsigned int, unsigned int>&
face_range)
const {
2195 *
for(
unsigned int d = 0;
d < dim; ++
d)
2207 *
if(boundary_id != 1) {
2211 *
for(
unsigned int i = 0; i <
phi.dofs_per_component; ++i) {
2212 *
for(
unsigned int j = 0;
j <
phi.dofs_per_component; ++
j)
2214 *
phi.submit_dof_value(tmp, i);
2218 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
2219 *
const auto&
n_plus =
phi.get_normal_vector(
q);
2232 *
for(
unsigned int i = 0; i <
phi.dofs_per_component; ++i)
2233 *
phi.submit_dof_value(diagonal[i], i);
2234 *
phi.distribute_local_to_global(dst);
2238 *
for(
unsigned int i = 0; i <
phi.dofs_per_component; ++i) {
2239 *
for(
unsigned int j = 0;
j <
phi.dofs_per_component; ++
j)
2241 *
phi.submit_dof_value(tmp, i);
2245 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
2246 *
const auto&
n_plus =
phi.get_normal_vector(
q);
2256 *
for(
unsigned int d = 0;
d < dim; ++
d)
2273 *
for(
unsigned int i = 0; i <
phi.dofs_per_component; ++i)
2274 *
phi.submit_dof_value(diagonal[i], i);
2275 *
phi.distribute_local_to_global(dst);
2285 *
for(
unsigned int d = 0;
d < dim; ++
d)
2297 *
if(boundary_id != 1) {
2301 *
for(
unsigned int i = 0; i <
phi.dofs_per_component; ++i) {
2302 *
for(
unsigned int j = 0;
j <
phi.dofs_per_component; ++
j)
2304 *
phi.submit_dof_value(tmp, i);
2308 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
2309 *
const auto&
n_plus =
phi.get_normal_vector(
q);
2311 *
const auto&
u =
phi.get_value(
q);
2322 *
for(
unsigned int i = 0; i <
phi.dofs_per_component; ++i)
2323 *
phi.submit_dof_value(diagonal[i], i);
2324 *
phi.distribute_local_to_global(dst);
2328 *
for(
unsigned int i = 0; i <
phi.dofs_per_component; ++i) {
2329 *
for(
unsigned int j = 0;
j <
phi.dofs_per_component; ++
j)
2331 *
phi.submit_dof_value(tmp, i);
2335 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
2336 *
const auto&
n_plus =
phi.get_normal_vector(
q);
2338 *
const auto&
u =
phi.get_value(
q);
2346 *
for(
unsigned int d = 0;
d < dim; ++
d)
2349 *
u_m[1][v] = -
u_m[1][v];
2363 *
for(
unsigned int i = 0; i <
phi.dofs_per_component; ++i)
2364 *
phi.submit_dof_value(diagonal[i], i);
2365 *
phi.distribute_local_to_global(dst);
2380 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
2384 *
const unsigned int& ,
2385 *
const std::pair<unsigned int, unsigned int>&
cell_range)
const {
2399 *
for(
unsigned int i = 0; i <
phi.dofs_per_component; ++i) {
2400 *
for(
unsigned int j = 0;
j <
phi.dofs_per_component; ++
j)
2408 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
2410 *
phi.submit_gradient(
phi.get_gradient(
q),
q);
2415 *
for(
unsigned int i = 0; i <
phi.dofs_per_component; ++i)
2416 *
phi.submit_dof_value(diagonal[i], i);
2418 *
phi.distribute_local_to_global(dst);
2431 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
2435 *
const unsigned int& ,
2436 *
const std::pair<unsigned int, unsigned int>&
face_range)
const {
2448 *
phi_p.reinit(face);
2449 *
phi_m.reinit(face);
2455 *
for(
unsigned int i = 0; i <
phi_p.dofs_per_component; ++i) {
2456 *
for(
unsigned int j = 0;
j <
phi_p.dofs_per_component; ++
j) {
2466 *
for(
unsigned int q = 0;
q <
phi_p.n_q_points; ++
q) {
2482 *
for(
unsigned int i = 0; i <
phi_p.dofs_per_component; ++i) {
2486 *
phi_p.distribute_local_to_global(dst);
2487 *
phi_m.distribute_local_to_global(dst);
2500 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
2504 *
const unsigned int& ,
2505 *
const std::pair<unsigned int, unsigned int>&
face_range)
const {
2513 *
if(boundary_id == 1) {
2518 *
for(
unsigned int i = 0; i <
phi.dofs_per_component; ++i) {
2519 *
for(
unsigned int j = 0;
j <
phi.dofs_per_component; ++
j)
2524 *
for(
unsigned int q = 0;
q <
phi.n_q_points; ++
q) {
2525 *
const auto&
n_plus =
phi.get_normal_vector(
q);
2528 *
const auto&
pres =
phi.get_value(
q);
2536 *
for(
unsigned int i = 0; i <
phi.dofs_per_component; ++i)
2537 *
phi.submit_dof_value(diagonal[i], i);
2538 *
phi.distribute_local_to_global(dst);
2555 *
template<
int dim,
int fe_degree_p,
int fe_degree_v,
int n_q_po
ints_1d_p,
int n_q_po
ints_1d_v,
typename Vec>
2564 * ::MatrixFreeTools::compute_diagonal<dim, Number, VectorizedArray<Number>>
2567 *
[&](
const auto&
data,
auto& dst,
const auto& src,
const auto&
cell_range) {
2570 *
[&](
const auto&
data,
auto& dst,
const auto& src,
const auto&
face_range) {
2579 * ::MatrixFreeTools::compute_diagonal<dim, Number, VectorizedArray<Number>>
2582 *
[&](
const auto&
data,
auto& dst,
const auto& src,
const auto&
cell_range) {
2585 *
[&](
const auto&
data,
auto& dst,
const auto& src,
const auto&
face_range) {
2594 *
for(
unsigned int i = 0; i <
inverse_diagonal.locally_owned_size(); ++i) {
2596 *
ExcMessage(
"No diagonal entry in a definite operator should be zero"));
2605 * <a name=
"navier_stokes_TRBDF2_DG.cc-ThecodeNavierStokesProjectioncodeclass"></a>
2627 *
const double gamma;
2671 *
<<
" The time step " <<
arg1 <<
" is out of range."
2673 *
<<
" The permitted range is (0," <<
arg2 <<
"]");
2679 *
void initialize();
2764 *
T(
data.final_time),
2797 *
<<
" WARNING: The chosen pair of finite element spaces is not stable."
2799 *
<<
" The obtained results will be nonsense" << std::endl;
2825 *
GridGenerator::plate_with_a_hole(
triangulation, 0.5, 1.0, 1.0, 1.1, 1.0, 19.0,
Point<2>(2.0, 2.0), 0, 1, 1.0, 2,
true);
2845 *
pcout <<
"Number of active cells: " <<
triangulation.n_global_active_cells() << std::endl;
2856 *
<<
"Re = " <<
Re << std::endl
2873 *
std::vector<const DoFHandler<dim>*> dof_handlers;
2883 *
std::vector<const AffineConstraints<double>*> constraints;
2943 *
const std::vector<unsigned int> tmp = {1};
2991 *
u_extr.equ(1.0 + gamma/(2.0*(1.0 - gamma)),
u_n);
3002 *
u_tmp.equ((1.0 - gamma)/gamma,
u_n);
3022 *
const std::vector<unsigned int> tmp = {0};
3074 *
const std::vector<unsigned int> tmp = {1};
3167 *
const std::vector<unsigned int> tmp = {0};
3172 *
else if(flag == 2)
3237 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
3239 *
u_n.update_ghost_values();
3241 *
pres_n.update_ghost_values();
3262 * <a name=
"navier_stokes_TRBDF2_DG.cc-codeNavierStokesProjectioncompute_lift_and_dragcode"></a>
3303 *
if(cell->is_locally_owned()) {
3305 *
if(cell->face(face)->at_boundary() && cell->face(face)->boundary_id() == 4) {
3312 *
for(
int q = 0;
q < n_q_points;
q++) {
3315 *
for(
unsigned int d = 0;
d < dim; ++
d) {
3317 *
for(
unsigned int k = 0;
k < dim; ++
k)
3347 * <a name=
"navier_stokes_TRBDF2_DG.cc-codeNavierStokesProjectionrefine_meshcode"></a>
3360 *
template <
int dim>
3376 *
const auto cell_worker = [&](
const Iterator& cell,
3378 *
CopyData& copy_data) {
3380 *
fe_values.reinit(cell);
3383 *
std::vector<std::vector<Tensor<1, dim>>>
gradients(fe_values.n_quadrature_points, std::vector<
Tensor<1, dim>>(dim));
3384 *
fe_values.get_function_gradients(
tmp_velocity, gradients);
3385 *
copy_data.cell_index = cell->active_cell_index();
3389 *
for(
unsigned k = 0;
k < fe_values.n_quadrature_points; ++
k) {
3398 *
const auto copier = [&](
const CopyData ©_data) {
3405 *
CopyData copy_data;
3417 *
for(
const auto& cell:
triangulation.active_cell_iterators()) {
3418 *
if(cell->refine_flag_set() &&
static_cast<unsigned int>(cell->level()) ==
max_loc_refinements)
3419 *
cell->clear_refine_flag();
3420 *
if(cell->coarsen_flag_set() &&
static_cast<unsigned int>(cell->level()) ==
min_loc_refinements)
3421 *
cell->clear_coarsen_flag();
3428 *
std::vector<const LinearAlgebra::distributed::Vector<double>*>
velocities;
3484 *
std::vector<const LinearAlgebra::distributed::Vector<double>*>
velocities;
3494 *
if(cell->is_locally_owned())
3495 *
cell->set_refine_flag();
3540 *
GridGenerator::plate_with_a_hole(
triangulation_tmp, 0.5, 1.0, 1.0, 1.1, 1.0, 19.0,
Point<2>(2.0, 2.0), 0, 1, 1.0, 2,
true);
3555 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
3558 *
u_n_tmp.update_ghost_values();
3567 *
const std::string output =
"./" +
saving_dir +
"/solution_max_res_end.vtu";
3575 * <a name=
"navier_stokes_TRBDF2_DG.cc-codeNavierStokesProjectionruncode"></a>
3599 *
double time =
t_0 +
dt;
3600 *
unsigned int n = 1;
3601 *
while(
std::abs(T - time) > 1e-10) {
3604 *
pcout <<
"Step = " << n <<
" Time = " << time << std::endl;
3612 *
verbose_cout <<
" Interpolating the velocity stage 1" << std::endl;
3615 *
verbose_cout <<
" Diffusion Step stage 1 " << std::endl;
3618 *
verbose_cout <<
" Projection Step stage 1" << std::endl;
3624 *
verbose_cout <<
" Updating the Velocity stage 1" << std::endl;
3638 *
verbose_cout <<
" Interpolating the velocity stage 2" << std::endl;
3641 *
verbose_cout <<
" Diffusion Step stage 2 " << std::endl;
3644 *
verbose_cout <<
" Projection Step stage 2" << std::endl;
3649 *
verbose_cout <<
" Updating the Velocity stage 2" << std::endl;
3656 *
pcout<<
"Maximal velocity = " <<
max_vel << std::endl;
3662 *
verbose_cout <<
"Plotting Solution final" << std::endl;
3679 *
verbose_cout <<
"Plotting Solution final" << std::endl;
3695 * <a name=
"navier_stokes_TRBDF2_DG.cc-Themainfunction"></a>
3713 *
data.read_data(
"parameter-file.prm");
3721 *
test.run(
data.verbose,
data.output_interval);
3724 *
std::cout <<
"----------------------------------------------------"
3726 *
<<
"Apparently everything went fine!" << std::endl
3727 *
<<
"Don't forget to brush your teeth :-)" << std::endl
3732 *
catch(std::exception &exc) {
3733 *
std::cerr << std::endl
3735 *
<<
"----------------------------------------------------"
3737 *
std::cerr <<
"Exception on processing: " << std::endl
3738 *
<< exc.what() << std::endl
3739 *
<<
"Aborting!" << std::endl
3740 *
<<
"----------------------------------------------------"
3745 *
std::cerr << std::endl
3747 *
<<
"----------------------------------------------------"
3749 *
std::cerr <<
"Unknown exception!" << std::endl
3750 *
<<
"Aborting!" << std::endl
3751 *
<<
"----------------------------------------------------"
3760<a name=
"ann-runtime_parameters.h"></a>
3790 * <a name=
"runtime_parameters.h-Runtimeparameters"></a>
3791 * <
h3>
Run time parameters</
h3>
3803 *
using namespace dealii;
3811 *
double initial_time;
3812 *
double final_time;
3823 *
unsigned int max_iterations;
3852 *
max_iterations(1000),
3857 *
prm.enter_subsection(
"Physical data");
3859 *
prm.declare_entry(
"initial_time",
3862 *
" The initial time of the simulation. ");
3863 *
prm.declare_entry(
"final_time",
3866 *
" The final time of the simulation. ");
3867 *
prm.declare_entry(
"Reynolds",
3870 *
" The Reynolds number. ");
3872 *
prm.leave_subsection();
3874 *
prm.enter_subsection(
"Time step data");
3876 *
prm.declare_entry(
"dt",
3879 *
" The time step size. ");
3881 *
prm.leave_subsection();
3883 *
prm.enter_subsection(
"Space discretization");
3885 *
prm.declare_entry(
"n_of_refines",
3888 *
" The number of cells we want on each direction of the mesh. ");
3889 *
prm.declare_entry(
"max_loc_refinements",
3892 *
" The number of maximum local refinements. ");
3893 *
prm.declare_entry(
"min_loc_refinements",
3896 *
" The number of minimum local refinements. ");
3898 *
prm.leave_subsection();
3900 *
prm.enter_subsection(
"Data solve");
3902 *
prm.declare_entry(
"max_iterations",
3905 *
" The maximal number of iterations linear solvers must make. ");
3906 *
prm.declare_entry(
"eps",
3909 *
" The stopping criterion. ");
3911 *
prm.leave_subsection();
3913 *
prm.declare_entry(
"refinement_iterations",
3916 *
" This number indicates how often we need to "
3917 *
"refine the mesh");
3919 *
prm.declare_entry(
"saving directory",
"SimTest");
3921 *
prm.declare_entry(
"verbose",
3924 *
" This indicates whether the output of the solution "
3925 *
"process should be verbose. ");
3927 *
prm.declare_entry(
"output_interval",
3930 *
" This indicates between how many time steps we print "
3931 *
"the solution. ");
3946 *
prm.parse_input(file);
3948 *
prm.enter_subsection(
"Physical data");
3950 *
initial_time = prm.get_double(
"initial_time");
3951 *
final_time = prm.get_double(
"final_time");
3952 *
Reynolds = prm.get_double(
"Reynolds");
3954 *
prm.leave_subsection();
3956 *
prm.enter_subsection(
"Time step data");
3958 *
dt = prm.get_double(
"dt");
3960 *
prm.leave_subsection();
3962 *
prm.enter_subsection(
"Space discretization");
3964 *
n_refines = prm.get_integer(
"n_of_refines");
3968 *
prm.leave_subsection();
3970 *
prm.enter_subsection(
"Data solve");
3972 *
max_iterations = prm.get_integer(
"max_iterations");
3973 *
eps = prm.get_double(
"eps");
3975 *
prm.leave_subsection();
3977 *
dir = prm.get(
"saving directory");
3981 *
verbose = prm.get_bool(
"verbose");
unsigned int depth_console(const unsigned int n)
Number linfty_norm(const Tensor< 2, dim, Number > &t)
static constexpr unsigned int dimension
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
static ::ExceptionBase & ExcFileNotOpen(std::string arg1)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define DeclException2(Exception2, type1, type2, outsequence)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void mesh_loop(const CellIteratorType &begin, const CellIteratorType &end, const CellWorkerFunctionType &cell_worker, const CopierType &copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const AssembleFlags flags=assemble_own_cells, const BoundaryWorkerFunctionType &boundary_worker=BoundaryWorkerFunctionType(), const FaceWorkerFunctionType &face_worker=FaceWorkerFunctionType(), const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
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())
@ 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.
std::vector< index_type > data
@ component_is_part_of_vector
void create_triangulation(Triangulation< dim, dim > &tria, const AdditionalData &additional_data=AdditionalData())
void plate_with_a_hole(Triangulation< dim > &tria, const double inner_radius=0.4, const double outer_radius=1., const double pad_bottom=2., const double pad_top=2., const double pad_left=1., const double pad_right=1., const Point< dim > ¢er=Point< dim >(), const types::manifold_id polar_manifold_id=0, const types::manifold_id tfi_manifold_id=1, const double L=1., const unsigned int n_slices=2, const bool colorize=false)
Rectangular plate with an (offset) cylindrical hole.
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
constexpr types::blas_int one
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
@ construct_multigrid_hierarchy
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
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)
long double gamma(const unsigned int n)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
constexpr unsigned int invalid_unsigned_int
void refine_and_coarsen_fixed_number(::Triangulation< dim, spacedim > &tria, const ::Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const types::global_cell_index max_n_cells=std::numeric_limits< types::global_cell_index >::max())
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
std_cxx20::type_identity< T > identity