1269 *
const unsigned int = 0)
const override
1274 *
virtual void vector_value(
const Point<dim> &p,
1277 *
for (
unsigned int c = 0; c < this->n_components; ++c)
1284 *
template <
int dim>
1297 *
ExcMessage(
"Invalid operation for a scalar function."));
1303 *
(dim == 2 ?
Point<dim>(.45, .1) :
Point<dim>(.45, .5, .1)),
1304 *
(dim == 2 ?
Point<dim>(.75, .1) :
Point<dim>(.75, .5, .1))};
1305 *
static const double source_radius = (dim == 2 ? 1. / 32 : 1. / 8);
1314 *
virtual void vector_value(
const Point<dim> &p,
1317 *
for (
unsigned int c = 0; c < this->n_components; ++c)
1328 * <a name=
"step_31-Linearsolversandpreconditioners"></a>
1340 *
that here we don't use the Schur complement to solve the Stokes
1341 * equations, though an approximate Schur complement (the mass matrix on the
1342 * pressure space) appears in the preconditioner.
1345 * namespace LinearSolvers
1350 * <a name="step_31-ThecodeInverseMatrixcodeclasstemplate"></a>
1351 * <h4>The <code>InverseMatrix</code> class template</h4>
1355 * This class is an interface to calculate the action of an "inverted"
1356 * matrix on a vector (using the <code>vmult</code> operation) in the same
1357 * way as the corresponding class in @ref step_22 "step-22": when the product of an
1358 * object of this class is requested, we solve a linear equation system
1359 * with that matrix using the CG method, accelerated by a preconditioner
1360 * of (templated) class <code>PreconditionerType</code>.
1364 * In a minor deviation from the implementation of the same class in
1365 * @ref step_22 "step-22", we make the <code>vmult</code> function take any
1366 * kind of vector type (it will yield compiler errors, however, if the
1367 * matrix does not allow a matrix-vector product with this kind of
1372 * Secondly, we catch any exceptions that the solver may have thrown. The
1373 * reason is as follows: When debugging a program like this one
1374 * occasionally makes a mistake of passing an indefinite or nonsymmetric
1375 * matrix or preconditioner to the current class. The solver will, in that
1376 * case, not converge and throw a run-time exception. If not caught here
1377 * it will propagate up the call stack and may end up in
1378 * <code>main()</code> where we output an error message that will say that
1379 * the CG solver failed. The question then becomes: Which CG solver? The
1380 * one that inverted the mass matrix? The one that inverted the top left
1381 * block with the Laplace operator? Or a CG solver in one of the several
1382 * other nested places where we use linear solvers in the current code? No
1383 * indication about this is present in a run-time exception because it
1408 *
const PreconditionerType &preconditioner);
1411 *
template <
typename VectorType>
1412 *
void vmult(VectorType &dst,
const VectorType &src)
const;
1416 *
const PreconditionerType &preconditioner;
1420 *
template <
class MatrixType,
class PreconditionerType>
1422 *
const MatrixType &m,
1423 *
const PreconditionerType &preconditioner)
1425 *
, preconditioner(preconditioner)
1430 *
template <
class MatrixType,
class PreconditionerType>
1431 *
template <
typename VectorType>
1434 *
const VectorType &src)
const
1436 *
SolverControl solver_control(src.size(), 1e-7 * src.l2_norm());
1443 *
cg.solve(*matrix, dst, src, preconditioner);
1445 *
catch (std::exception &e)
1454 * <a name=
"step_31-Schurcomplementpreconditioner"></a>
1467 *
Let's have a look at the ideal preconditioner matrix
1468 * @f$P=\left(\begin{array}{cc} A & 0 \\ B & -S \end{array}\right)@f$
1469 * described in the introduction. If we apply this matrix in the solution
1470 * of a linear system, convergence of an iterative GMRES solver will be
1471 * governed by the matrix @f{eqnarray*} P^{-1}\left(\begin{array}{cc} A &
1472 * B^T \\ B & 0 \end{array}\right) = \left(\begin{array}{cc} I & A^{-1}
1473 * B^T \\ 0 & I \end{array}\right), @f} which indeed is very simple. A
1474 * GMRES solver based on exact matrices would converge in one iteration,
1475 * since all eigenvalues are equal (any Krylov method takes at most as
1476 * many iterations as there are distinct eigenvalues). Such a
1477 * preconditioner for the blocked Stokes system has been proposed by
1478 * Silvester and Wathen ("Fast iterative solution of stabilised Stokes
1479 * systems part II. Using general block preconditioners", SIAM
1480 * J. Numer. Anal., 31 (1994), pp. 1352-1367).
1484 * Replacing @f$P@f$ by @f$\tilde{P}@f$ keeps that spirit alive: the product
1485 * @f$P^{-1} A@f$ will still be close to a matrix with eigenvalues 1 with a
1486 * distribution that does not depend on the problem size. This lets us
1487 * hope to be able to get a number of GMRES iterations that is
1488 * problem-size independent.
1492 * The deal.II users who have already gone through the @ref step_20 "step-20" and @ref step_22 "step-22"
1504 * (
using the AMG preconditioner).
1595 *
template <
class PreconditionerTypeA,
class PreconditionerTypeMp>
1602 *
stokes_matrix->block(1, 0).residual(tmp, dst.block(0), src.block(1));
1613 * <a name=
"step_31-ThecodeBoussinesqFlowProblemcodeclasstemplate"></a>
1645 *
template <
int dim>
1677 *
const double cell_diameter)
const;
1728 * <a name=
"step_31-BoussinesqFlowProblemclassimplementation"></a>
1734 * <a name=
"step_31-BoussinesqFlowProblemBoussinesqFlowProblem"></a>
1748 *
and preconditioning:
1751 *
template <
int dim>
1778 * <a name=
"step_31-BoussinesqFlowProblemget_maximal_velocity"></a>
1848 *
fe_values.
reinit(cell);
1852 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
1864 * <a name=
"step_31-BoussinesqFlowProblemget_extrapolated_temperature_range"></a>
1898 *
template <
int dim>
1899 *
std::pair<double, double>
1917 *
fe_values.
reinit(cell);
1923 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
1943 *
fe_values.
reinit(cell);
1947 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
1965 * <a name=
"step_31-BoussinesqFlowProblemcompute_viscosity"></a>
1997 *
template <
int dim>
2010 *
const double cell_diameter)
const
2012 *
constexpr double beta = 0.017 * dim;
2013 *
constexpr double alpha = 1.0;
2016 *
return 5
e-3 * cell_diameter;
2023 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
2028 *
const double dT_dt =
2038 *
const double residual =
2062 * <a name=
"step_31-BoussinesqFlowProblemsetup_dofs"></a>
2080 *
dependent on ILU's, whereas we use AMG here which is not sensitive to the
2081 * DoF numbering. The IC preconditioner for the inversion of the pressure
2082 * mass matrix would of course take advantage of a Cuthill-McKee like
2083 * renumbering, but its costs are low compared to the velocity portion, so
2084 * the additional work does not pay off.
2088 * We then proceed with the generation of the hanging node constraints that
2089 * arise from adaptive grid refinement for both DoFHandler objects. For the
2090 * velocity, we impose no-flux boundary conditions @f$\mathbf{u}\cdot
2091 * \mathbf{n}=0@f$ by adding constraints to the object that already stores the
2092 * hanging node constraints matrix. The second parameter in the function
2093 * describes the first of the velocity components in the total dof vector,
2094 * which is zero here. The variable <code>no_normal_flux_boundaries</code>
2095 * denotes the boundary indicators for which to set the no flux boundary
2096 * conditions; here, this is boundary indicator zero.
2100 * After having done so, we count the number of degrees of freedom in the
2104 * template <int dim>
2105 * void BoussinesqFlowProblem<dim>::setup_dofs()
2107 * std::vector<unsigned int> stokes_sub_blocks(dim + 1, 0);
2108 * stokes_sub_blocks[dim] = 1;
2111 * stokes_dof_handler.distribute_dofs(stokes_fe);
2112 * DoFRenumbering::component_wise(stokes_dof_handler, stokes_sub_blocks);
2114 * stokes_constraints.clear();
2115 * DoFTools::make_hanging_node_constraints(stokes_dof_handler,
2116 * stokes_constraints);
2117 * const std::set<types::boundary_id> no_normal_flux_boundaries = {0};
2118 * VectorTools::compute_no_normal_flux_constraints(stokes_dof_handler,
2120 * no_normal_flux_boundaries,
2121 * stokes_constraints);
2122 * stokes_constraints.close();
2125 * temperature_dof_handler.distribute_dofs(temperature_fe);
2127 * temperature_constraints.clear();
2128 * DoFTools::make_hanging_node_constraints(temperature_dof_handler,
2129 * temperature_constraints);
2130 * temperature_constraints.close();
2133 * const std::vector<types::global_dof_index> stokes_dofs_per_block =
2134 * DoFTools::count_dofs_per_fe_block(stokes_dof_handler, stokes_sub_blocks);
2136 * const types::global_dof_index n_u = stokes_dofs_per_block[0],
2137 * n_p = stokes_dofs_per_block[1],
2138 * n_T = temperature_dof_handler.n_dofs();
2140 * std::cout << "Number of active cells: " << triangulation.n_active_cells()
2141 * << " (on " << triangulation.n_levels() << " levels)" << std::endl
2142 * << "Number of degrees of freedom: " << n_u + n_p + n_T << " ("
2143 * << n_u << '+
' << n_p << '+
' << n_T << ')
' << std::endl
2148 * The next step is to create the sparsity pattern for the Stokes and
2149 * temperature system matrices as well as the preconditioner matrix from
2150 * which we build the Stokes preconditioner. As in @ref step_22 "step-22", we choose to
2151 * create the pattern by
2152 * using the blocked version of DynamicSparsityPattern.
2156 * So, we first release the memory stored in the matrices, then set up an
2157 * object of type BlockDynamicSparsityPattern consisting of
2158 * @f$2\times 2@f$ blocks (for the Stokes system matrix and preconditioner) or
2159 * DynamicSparsityPattern (for the temperature part). We then
2160 * fill these objects with the nonzero pattern, taking into account that
2161 * for the Stokes system matrix, there are no entries in the
2162 * pressure-pressure block (but all velocity vector components couple with
2163 * each other and with the pressure). Similarly, in the Stokes
2164 * preconditioner matrix, only the diagonal blocks are nonzero, since we
2165 * use the vector Laplacian as discussed in the introduction. This
2166 * operator only couples each vector component of the Laplacian with
2167 * itself, but not with the other vector components. (Application of the
2168 * constraints resulting from the no-flux boundary conditions will couple
2169 * vector components at the boundary again, however.)
2173 * When generating the sparsity pattern, we directly apply the constraints
2174 * from hanging nodes and no-flux boundary conditions. This approach was
2175 * already used in @ref step_27 "step-27", but is different from the one in early
2176 * tutorial programs where we first built the original sparsity pattern
2177 * and only then added the entries resulting from constraints. The reason
2178 * for doing so is that later during assembly we are going to distribute
2179 * the constraints immediately when transferring local to global
2180 * dofs. Consequently, there will be no data written at positions of
2181 * constrained degrees of freedom, so we can let the
2182 * DoFTools::make_sparsity_pattern function omit these entries by setting
2183 * the last Boolean flag to <code>false</code>. Once the sparsity pattern
2184 * is ready, we can use it to initialize the Trilinos matrices. Since the
2185 * Trilinos matrices store the sparsity pattern internally, there is no
2186 * need to keep the sparsity pattern around after the initialization of
2190 * stokes_partitioning.resize(2);
2191 * stokes_partitioning[0] = complete_index_set(n_u);
2192 * stokes_partitioning[1] = complete_index_set(n_p);
2194 * stokes_matrix.clear();
2196 * BlockDynamicSparsityPattern dsp(stokes_dofs_per_block,
2197 * stokes_dofs_per_block);
2199 * Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
2201 * for (unsigned int c = 0; c < dim + 1; ++c)
2202 * for (unsigned int d = 0; d < dim + 1; ++d)
2203 * if (!((c == dim) && (d == dim)))
2204 * coupling[c][d] = DoFTools::always;
2206 * coupling[c][d] = DoFTools::none;
2208 * DoFTools::make_sparsity_pattern(
2209 * stokes_dof_handler, coupling, dsp, stokes_constraints, false);
2211 * stokes_matrix.reinit(dsp);
2215 * Amg_preconditioner.reset();
2216 * Mp_preconditioner.reset();
2217 * stokes_preconditioner_matrix.clear();
2219 * BlockDynamicSparsityPattern dsp(stokes_dofs_per_block,
2220 * stokes_dofs_per_block);
2222 * Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
2223 * for (unsigned int c = 0; c < dim + 1; ++c)
2224 * for (unsigned int d = 0; d < dim + 1; ++d)
2226 * coupling[c][d] = DoFTools::always;
2228 * coupling[c][d] = DoFTools::none;
2230 * DoFTools::make_sparsity_pattern(
2231 * stokes_dof_handler, coupling, dsp, stokes_constraints, false);
2233 * stokes_preconditioner_matrix.reinit(dsp);
2238 * The creation of the temperature matrix (or, rather, matrices, since we
2239 * provide a temperature mass matrix and a temperature @ref GlossStiffnessMatrix "stiffness matrix",
2240 * that will be added together for time discretization) follows the
2241 * generation of the Stokes matrix – except that it is much easier
2242 * here since we do not need to take care of any blocks or coupling
2243 * between components. Note how we initialize the three temperature
2244 * matrices: We only use the sparsity pattern for reinitialization of the
2245 * first matrix, whereas we use the previously generated matrix for the
2246 * two remaining reinits. The reason for doing so is that reinitialization
2247 * from an already generated matrix allows Trilinos to reuse the sparsity
2248 * pattern instead of generating a new one for each copy. This saves both
2249 * some time and memory.
2253 * temperature_mass_matrix.clear();
2254 * temperature_stiffness_matrix.clear();
2255 * temperature_matrix.clear();
2257 * DynamicSparsityPattern dsp(n_T, n_T);
2258 * DoFTools::make_sparsity_pattern(temperature_dof_handler,
2260 * temperature_constraints,
2263 * temperature_matrix.reinit(dsp);
2264 * temperature_mass_matrix.reinit(temperature_matrix);
2265 * temperature_stiffness_matrix.reinit(temperature_matrix);
2270 * Lastly, we set the vectors for the Stokes solutions @f$\mathbf u^{n-1}@f$
2271 * and @f$\mathbf u^{n-2}@f$, as well as for the temperatures @f$T^{n}@f$,
2272 * @f$T^{n-1}@f$ and @f$T^{n-2}@f$ (required for time stepping) and all the system
2273 * right hand sides to their correct sizes and block structure:
2276 * IndexSet temperature_partitioning = complete_index_set(n_T);
2277 * stokes_solution.reinit(stokes_partitioning, MPI_COMM_WORLD);
2278 * old_stokes_solution.reinit(stokes_partitioning, MPI_COMM_WORLD);
2279 * stokes_rhs.reinit(stokes_partitioning, MPI_COMM_WORLD);
2281 * temperature_solution.reinit(temperature_partitioning, MPI_COMM_WORLD);
2282 * old_temperature_solution.reinit(temperature_partitioning, MPI_COMM_WORLD);
2283 * old_old_temperature_solution.reinit(temperature_partitioning,
2286 * temperature_rhs.reinit(temperature_partitioning, MPI_COMM_WORLD);
2294 * <a name="step_31-BoussinesqFlowProblemassemble_stokes_preconditioner"></a>
2295 * <h4>BoussinesqFlowProblem::assemble_stokes_preconditioner</h4>
2299 * This function assembles the matrix we use for preconditioning the Stokes
2300 * system. What we need are a vector Laplace matrix on the velocity
2301 * components and a mass matrix weighted by @f$\eta^{-1}@f$ on the pressure
2302 * component. We start by generating a quadrature object of appropriate
2303 * order, the FEValues object that can give values and gradients at the
2304 * quadrature points (together with quadrature weights). Next we create data
2305 * structures for the cell matrix and the relation between local and global
2306 * DoFs. The vectors <code>grad_phi_u</code> and <code>phi_p</code> are
2307 * going to hold the values of the basis functions in order to faster build
2308 * up the local matrices, as was already done in @ref step_22 "step-22". Before we start
2309 * the loop over all active cells, we have to specify which components are
2310 * pressure and which are velocity.
2313 * template <int dim>
2314 * void BoussinesqFlowProblem<dim>::assemble_stokes_preconditioner()
2316 * stokes_preconditioner_matrix = 0;
2318 * const QGauss<dim> quadrature_formula(stokes_degree + 2);
2319 * FEValues<dim> stokes_fe_values(stokes_fe,
2320 * quadrature_formula,
2321 * update_JxW_values | update_values |
2322 * update_gradients);
2324 * const unsigned int dofs_per_cell = stokes_fe.n_dofs_per_cell();
2325 * const unsigned int n_q_points = quadrature_formula.size();
2327 * FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
2328 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2330 * std::vector<Tensor<2, dim>> grad_phi_u(dofs_per_cell);
2331 * std::vector<double> phi_p(dofs_per_cell);
2333 * const FEValuesExtractors::Vector velocities(0);
2334 * const FEValuesExtractors::Scalar pressure(dim);
2336 * for (const auto &cell : stokes_dof_handler.active_cell_iterators())
2338 * stokes_fe_values.reinit(cell);
2343 * The creation of the local matrix is rather simple. There are only a
2344 * Laplace term (on the velocity) and a mass matrix weighted by
2345 * @f$\eta^{-1}@f$ to be generated, so the creation of the local matrix is
2346 * done in two lines. Once the local matrix is ready (loop over rows
2347 * and columns in the local matrix on each quadrature point), we get
2348 * the local DoF indices and write the local information into the
2349 * global matrix. We do this as in @ref step_27 "step-27", i.e., we directly apply the
2350 * constraints from hanging nodes locally. By doing so, we don't
have
2352 * matrix that will actually be set to zero again later when
2353 * eliminating constraints.
2356 * for (unsigned int q = 0; q < n_q_points; ++q)
2358 * for (unsigned int k = 0; k < dofs_per_cell; ++k)
2360 * grad_phi_u[k] = stokes_fe_values[velocities].gradient(k, q);
2361 * phi_p[k] = stokes_fe_values[pressure].value(k, q);
2364 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2365 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
2366 * local_matrix(i, j) +=
2367 * (EquationData::eta *
2368 * scalar_product(grad_phi_u[i], grad_phi_u[j]) +
2369 * (1. / EquationData::eta) * phi_p[i] * phi_p[j]) *
2370 * stokes_fe_values.JxW(q);
2373 * cell->get_dof_indices(local_dof_indices);
2374 * stokes_constraints.distribute_local_to_global(
2375 * local_matrix, local_dof_indices, stokes_preconditioner_matrix);
2384 * <a name="step_31-BoussinesqFlowProblembuild_stokes_preconditioner"></a>
2385 * <h4>BoussinesqFlowProblem::build_stokes_preconditioner</h4>
2389 * This function generates the inner preconditioners that are going to be
2390 * used for the Schur complement block preconditioner. Since the
2391 * preconditioners need only to be regenerated when the matrices change,
2392 * this function does not have to do anything in case the matrices have not
2393 * changed (i.e., the flag <code>rebuild_stokes_preconditioner</code> has
2394 * the value <code>false</code>). Otherwise its first task is to call
2395 * <code>assemble_stokes_preconditioner</code> to generate the
2396 * preconditioner matrices.
2400 * Next, we set up the preconditioner for the velocity-velocity matrix
2401 * @f$A@f$. As explained in the introduction, we are going to use an AMG
2402 * preconditioner based on a vector Laplace matrix @f$\hat{A}@f$ (which is
2403 * spectrally close to the Stokes matrix @f$A@f$). Usually, the
2404 * TrilinosWrappers::PreconditionAMG class can be seen as a good black-box
2405 * preconditioner which does not need any special knowledge. In this case,
2406 * however, we have to be careful: since we build an AMG for a vector
2407 * problem, we have to tell the preconditioner setup which dofs belong to
2408 * which vector component. We do this using the function
2409 * DoFTools::extract_constant_modes, a function that generates a set of
2410 * <code>dim</code> vectors, where each one has ones in the respective
2411 * component of the vector problem and zeros elsewhere. Hence, these are the
2412 * constant modes on each component, which explains the name of the
2416 * template <int dim>
2417 * void BoussinesqFlowProblem<dim>::build_stokes_preconditioner()
2419 * if (rebuild_stokes_preconditioner == false)
2422 * std::cout << " Rebuilding Stokes preconditioner..." << std::flush;
2424 * assemble_stokes_preconditioner();
2426 * Amg_preconditioner = std::make_shared<TrilinosWrappers::PreconditionAMG>();
2428 * const FEValuesExtractors::Vector velocity_components(0);
2429 * const std::vector<std::vector<bool>> constant_modes =
2430 * DoFTools::extract_constant_modes(
2431 * stokes_dof_handler, stokes_fe.component_mask(velocity_components));
2432 * TrilinosWrappers::PreconditionAMG::AdditionalData amg_data;
2433 * amg_data.constant_modes = constant_modes;
2437 * Next, we set some more options of the AMG preconditioner. In
2438 * particular, we need to tell the AMG setup that we use quadratic basis
2439 * functions for the velocity matrix (this implies more nonzero elements
2440 * in the matrix, so that a more robust algorithm needs to be chosen
2441 * internally). Moreover, we want to be able to control how the coarsening
2442 * structure is build up. The way the Trilinos smoothed aggregation AMG
2443 * does this is to look which matrix entries are of similar size as the
2444 * diagonal entry in order to algebraically build a coarse-grid
2445 * structure. By setting the parameter <code>aggregation_threshold</code>
2446 * to 0.02, we specify that all entries that are more than two percent of
2447 * size of some diagonal pivots in that row should form one coarse grid
2448 * point. This parameter is rather ad hoc, and some fine-tuning of it can
2449 * influence the performance of the preconditioner. As a rule of thumb,
2450 * larger values of <code>aggregation_threshold</code> will decrease the
2451 * number of iterations, but increase the costs per iteration. A look at
2452 * the Trilinos documentation will provide more information on these
2453 * parameters. With this data set, we then initialize the preconditioner
2454 * with the matrix we want it to apply to.
2458 * Finally, we also initialize the preconditioner for the inversion of the
2459 * pressure mass matrix. This matrix is symmetric and well-behaved, so we
2460 * can chose a simple preconditioner. We stick with an incomplete Cholesky
2461 * (IC) factorization preconditioner, which is designed for symmetric
2462 * matrices. We could have also chosen an SSOR preconditioner with
2463 * relaxation factor around 1.2, but IC is cheaper for our example. We
2464 * wrap the preconditioners into a <code>std::shared_ptr</code>
2465 * pointer, which makes it easier to recreate the preconditioner next time
2466 * around since we do not have to care about destroying the previously
2470 * amg_data.elliptic = true;
2471 * amg_data.higher_order_elements = true;
2472 * amg_data.smoother_sweeps = 2;
2473 * amg_data.aggregation_threshold = 0.02;
2474 * Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0, 0),
2477 * Mp_preconditioner = std::make_shared<TrilinosWrappers::PreconditionIC>();
2478 * Mp_preconditioner->initialize(stokes_preconditioner_matrix.block(1, 1));
2480 * std::cout << std::endl;
2482 * rebuild_stokes_preconditioner = false;
2490 * <a name="step_31-BoussinesqFlowProblemassemble_stokes_system"></a>
2491 * <h4>BoussinesqFlowProblem::assemble_stokes_system</h4>
2495 * The time lag scheme we use for advancing the coupled Stokes-temperature
2496 * system forces us to split up the assembly (and the solution of linear
2497 * systems) into two step. The first one is to create the Stokes system
2498 * matrix and right hand side, and the second is to create matrix and right
2499 * hand sides for the temperature dofs, which depends on the result of the
2500 * linear system for the velocity.
2504 * This function is called at the beginning of each time step. In the first
2505 * time step or if the mesh has changed, indicated by the
2506 * <code>rebuild_stokes_matrix</code>, we need to assemble the Stokes
2528 * use the temperature structures and set an update flag for the basis
2529 * function values which we need for evaluation of the temperature
2530 * solution. The only important part to remember here is that the same
2531 * quadrature formula is used for both FEValues objects to ensure that we
2532 * get matching information when we loop over the quadrature points of the
2537 * The declarations proceed with some shortcuts for array sizes, the
2538 * creation of the local matrix and right hand side as well as the vector
2539 * for the indices of the local dofs compared to the global system.
2542 * template <int dim>
2543 * void BoussinesqFlowProblem<dim>::assemble_stokes_system()
2545 * std::cout << " Assembling..." << std::flush;
2547 * if (rebuild_stokes_matrix == true)
2548 * stokes_matrix = 0;
2552 * const QGauss<dim> quadrature_formula(stokes_degree + 2);
2553 * FEValues<dim> stokes_fe_values(
2555 * quadrature_formula,
2556 * update_values | update_quadrature_points | update_JxW_values |
2557 * (rebuild_stokes_matrix == true ? update_gradients : UpdateFlags(0)));
2559 * FEValues<dim> temperature_fe_values(temperature_fe,
2560 * quadrature_formula,
2563 * const unsigned int dofs_per_cell = stokes_fe.n_dofs_per_cell();
2564 * const unsigned int n_q_points = quadrature_formula.size();
2566 * FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
2567 * Vector<double> local_rhs(dofs_per_cell);
2569 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2573 * Next we need a vector that will contain the values of the temperature
2574 * solution at the previous time level at the quadrature points to
2575 * assemble the source term in the right hand side of the momentum
2595 *
std::vector<Tensor<1, dim>>
phi_u(dofs_per_cell);
2596 *
std::vector<SymmetricTensor<2, dim>>
grads_phi_u(dofs_per_cell);
2597 *
std::vector<double>
div_phi_u(dofs_per_cell);
2598 *
std::vector<double>
phi_p(dofs_per_cell);
2610 * in sync. The first statements within the loop are again all very
2611 * familiar, doing the update of the finite element data as specified by
2612 * the update flags, zeroing out the local arrays and getting the values
2613 * of the old solution at the quadrature points. Then we are ready to loop
2614 * over the quadrature points on the cell.
2617 * auto cell = stokes_dof_handler.begin_active();
2618 * const auto endc = stokes_dof_handler.end();
2619 * auto temperature_cell = temperature_dof_handler.begin_active();
2621 * for (; cell != endc; ++cell, ++temperature_cell)
2623 * stokes_fe_values.reinit(cell);
2624 * temperature_fe_values.reinit(temperature_cell);
2629 * temperature_fe_values.get_function_values(old_temperature_solution,
2630 * old_temperature_values);
2632 * for (unsigned int q = 0; q < n_q_points; ++q)
2634 * const double old_temperature = old_temperature_values[q];
2638 * Next we extract the values and gradients of basis functions
2639 * relevant to the terms in the inner products. As shown in
2640 * @ref step_22 "step-22" this helps accelerate assembly.
2644 * Once this is done, we start the loop over the rows and columns
2645 * of the local matrix and feed the matrix with the relevant
2646 * products. The right hand side is filled with the forcing term
2647 * driven by temperature in direction of gravity (which is
2648 * vertical in our example). Note that the right hand side term
2649 * is always generated, whereas the matrix contributions are only
2650 * updated when it is requested by the
2651 * <code>rebuild_matrices</code> flag.
2654 * for (unsigned int k = 0; k < dofs_per_cell; ++k)
2656 * phi_u[k] = stokes_fe_values[velocities].value(k, q);
2657 * if (rebuild_stokes_matrix)
2660 * stokes_fe_values[velocities].symmetric_gradient(k, q);
2662 * stokes_fe_values[velocities].divergence(k, q);
2663 * phi_p[k] = stokes_fe_values[pressure].value(k, q);
2667 * if (rebuild_stokes_matrix)
2668 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2669 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
2670 * local_matrix(i, j) +=
2671 * (EquationData::eta * 2 * (grads_phi_u[i] * grads_phi_u[j]) -
2672 * div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) *
2673 * stokes_fe_values.JxW(q);
2675 * const Point<dim> gravity =
2676 * -((dim == 2) ? (Point<dim>(0, 1)) : (Point<dim>(0, 0, 1)));
2677 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
2678 * local_rhs(i) += (-EquationData::density * EquationData::beta *
2679 * gravity * phi_u[i] * old_temperature) *
2680 * stokes_fe_values.JxW(q);
2685 * The last step in the loop over all cells is to enter the local
2686 * contributions into the global matrix and vector structures to the
2687 * positions specified in <code>local_dof_indices</code>. Again, we
2688 * let the AffineConstraints class do the insertion of the cell
2689 * matrix elements to the global matrix, which already condenses the
2690 * hanging node constraints.
2693 * cell->get_dof_indices(local_dof_indices);
2695 * if (rebuild_stokes_matrix == true)
2696 * stokes_constraints.distribute_local_to_global(local_matrix,
2698 * local_dof_indices,
2702 * stokes_constraints.distribute_local_to_global(local_rhs,
2703 * local_dof_indices,
2707 * rebuild_stokes_matrix = false;
2709 * std::cout << std::endl;
2717 * <a name="step_31-BoussinesqFlowProblemassemble_temperature_matrix"></a>
2718 * <h4>BoussinesqFlowProblem::assemble_temperature_matrix</h4>
2722 * This function assembles the matrix in the temperature equation. The
2723 * temperature matrix consists of two parts, a mass matrix and the time step
2724 * size times a stiffness matrix given by a Laplace term times the amount of
2725 * diffusion. Since the matrix depends on the time step size (which varies
2726 * from one step to another), the temperature matrix needs to be updated
2727 * every time step. We could simply regenerate the matrices in every time
2728 * step, but this is not really efficient since mass and Laplace matrix do
2729 * only change when we change the mesh. Hence, we do this more efficiently
2730 * by generating two separate matrices in this function, one for the mass
2731 * matrix and one for the stiffness (diffusion) matrix. We will then sum up
2732 * the matrix plus the stiffness matrix times the time step size once we
2733 * know the actual time step.
2737 * So the details for this first step are very simple. In case we need to
2738 * rebuild the matrix (i.e., the mesh has changed), we zero the data
2739 * structures, get a quadrature formula and a FEValues object, and create
2740 * local matrices, local dof indices and evaluation structures for the basis
2744 * template <int dim>
2745 * void BoussinesqFlowProblem<dim>::assemble_temperature_matrix()
2747 * if (rebuild_temperature_matrices == false)
2750 * temperature_mass_matrix = 0;
2751 * temperature_stiffness_matrix = 0;
2753 * const QGauss<dim> quadrature_formula(temperature_degree + 2);
2754 * FEValues<dim> temperature_fe_values(temperature_fe,
2755 * quadrature_formula,
2756 * update_values | update_gradients |
2757 * update_JxW_values);
2759 * const unsigned int dofs_per_cell = temperature_fe.n_dofs_per_cell();
2760 * const unsigned int n_q_points = quadrature_formula.size();
2762 * FullMatrix<double> local_mass_matrix(dofs_per_cell, dofs_per_cell);
2763 * FullMatrix<double> local_stiffness_matrix(dofs_per_cell, dofs_per_cell);
2765 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2767 * std::vector<double> phi_T(dofs_per_cell);
2768 * std::vector<Tensor<1, dim>> grad_phi_T(dofs_per_cell);
2789 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
2791 *
for (
unsigned int k = 0;
k < dofs_per_cell; ++
k)
2797 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
2798 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
2808 *
cell->get_dof_indices(local_dof_indices);
2814 *
local_dof_indices,
2826 * <a name=
"step_31-BoussinesqFlowProblemassemble_temperature_system"></a>
2849 *
template <
int dim>
2881 *
const unsigned int dofs_per_cell =
temperature_fe.n_dofs_per_cell();
2886 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
2913 *
std::vector<double>
phi_T(dofs_per_cell);
2914 *
std::vector<Tensor<1, dim>>
grad_phi_T(dofs_per_cell);
2923 *
Now,
let's start the loop over all cells in the triangulation. Again,
2924 * we need two cell iterators that walk in parallel through the cells of
2925 * the two involved DoFHandler objects for the Stokes and temperature
2926 * part. Within the loop, we first set the local rhs to zero, and then get
2927 * the values and derivatives of the old solution functions at the
2928 * quadrature points, since they are going to be needed for the definition
2929 * of the stabilization parameters and as coefficients in the equation,
2930 * respectively. Note that since the temperature has its own DoFHandler
2931 * and FEValues object we get the entire solution at the quadrature point
2932 * (which is the scalar temperature field only anyway) whereas for the
2933 * Stokes part we restrict ourselves to extracting the velocity part (and
2934 * ignoring the pressure part) by using
2935 * <code>stokes_fe_values[velocities].get_function_values</code>.
2938 * auto cell = temperature_dof_handler.begin_active();
2939 * const auto endc = temperature_dof_handler.end();
2940 * auto stokes_cell = stokes_dof_handler.begin_active();
2942 * for (; cell != endc; ++cell, ++stokes_cell)
2946 * temperature_fe_values.reinit(cell);
2947 * stokes_fe_values.reinit(stokes_cell);
2949 * temperature_fe_values.get_function_values(old_temperature_solution,
2950 * old_temperature_values);
2951 * temperature_fe_values.get_function_values(old_old_temperature_solution,
2952 * old_old_temperature_values);
2954 * temperature_fe_values.get_function_gradients(old_temperature_solution,
2955 * old_temperature_grads);
2956 * temperature_fe_values.get_function_gradients(
2957 * old_old_temperature_solution, old_old_temperature_grads);
2959 * temperature_fe_values.get_function_laplacians(
2960 * old_temperature_solution, old_temperature_laplacians);
2961 * temperature_fe_values.get_function_laplacians(
2962 * old_old_temperature_solution, old_old_temperature_laplacians);
2964 * temperature_right_hand_side.value_list(
2965 * temperature_fe_values.get_quadrature_points(), gamma_values);
2967 * stokes_fe_values[velocities].get_function_values(stokes_solution,
2968 * old_velocity_values);
2969 * stokes_fe_values[velocities].get_function_values(
2970 * old_stokes_solution, old_old_velocity_values);
2974 * Next, we calculate the artificial viscosity for stabilization
2975 * according to the discussion in the introduction using the dedicated
2976 * function. With that at hand, we can get into the loop over
2977 * quadrature points and local rhs vector components. The terms here
2978 * are quite lengthy, but their definition follows the time-discrete
2979 * system developed in the introduction of this program. The BDF-2
2980 * scheme needs one more term from the old time step (and involves
2981 * more complicated factors) than the backward Euler scheme that is
2982 * used for the first time step. When all this is done, we distribute
2983 * the local vector into the global one (including hanging node
2988 * compute_viscosity(old_temperature_values,
2989 * old_old_temperature_values,
2990 * old_temperature_grads,
2991 * old_old_temperature_grads,
2992 * old_temperature_laplacians,
2993 * old_old_temperature_laplacians,
2994 * old_velocity_values,
2995 * old_old_velocity_values,
2998 * global_T_range.second - global_T_range.first,
2999 * cell->diameter());
3001 * for (unsigned int q = 0; q < n_q_points; ++q)
3003 * for (unsigned int k = 0; k < dofs_per_cell; ++k)
3005 * grad_phi_T[k] = temperature_fe_values.shape_grad(k, q);
3006 * phi_T[k] = temperature_fe_values.shape_value(k, q);
3009 * const double T_term_for_rhs =
3010 * (use_bdf2_scheme ?
3011 * (old_temperature_values[q] * (1 + time_step / old_time_step) -
3012 * old_old_temperature_values[q] * (time_step * time_step) /
3013 * (old_time_step * (time_step + old_time_step))) :
3014 * old_temperature_values[q]);
3016 * const Tensor<1, dim> ext_grad_T =
3017 * (use_bdf2_scheme ?
3018 * (old_temperature_grads[q] * (1 + time_step / old_time_step) -
3019 * old_old_temperature_grads[q] * time_step / old_time_step) :
3020 * old_temperature_grads[q]);
3022 * const Tensor<1, dim> extrapolated_u =
3023 * (use_bdf2_scheme ?
3024 * (old_velocity_values[q] * (1 + time_step / old_time_step) -
3025 * old_old_velocity_values[q] * time_step / old_time_step) :
3026 * old_velocity_values[q]);
3028 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
3030 * (T_term_for_rhs * phi_T[i] -
3031 * time_step * extrapolated_u * ext_grad_T * phi_T[i] -
3032 * time_step * nu * ext_grad_T * grad_phi_T[i] +
3033 * time_step * gamma_values[q] * phi_T[i]) *
3034 * temperature_fe_values.JxW(q);
3037 * cell->get_dof_indices(local_dof_indices);
3038 * temperature_constraints.distribute_local_to_global(local_rhs,
3039 * local_dof_indices,
3049 * <a name="step_31-BoussinesqFlowProblemsolve"></a>
3050 * <h4>BoussinesqFlowProblem::solve</h4>
3054 * This function solves the linear systems of equations. Following the
3055 * introduction, we start with the Stokes system, where we need to generate
3056 * our block Schur preconditioner. Since all the relevant actions are
3057 * implemented in the class <code>BlockSchurPreconditioner</code>, all we
3058 * have to do is to initialize the class appropriately. What we need to pass
3059 * down is an <code>InverseMatrix</code> object for the pressure mass
3060 * matrix, which we set up using the respective class together with the IC
3061 * preconditioner we already generated, and the AMG preconditioner for the
3062 * velocity-velocity matrix. Note that both <code>Mp_preconditioner</code>
3063 * and <code>Amg_preconditioner</code> are only pointers, so we use
3064 * <code>*</code> to pass down the actual preconditioner objects.
3068 * Once the preconditioner is ready, we create a GMRES solver for the block
3069 * system. Since we are working with Trilinos data structures, we have to
3070 * set the respective template argument in the solver. GMRES needs to
3071 * internally store temporary vectors for each iteration (see the discussion
3072 * in the results section of @ref step_22 "step-22") – the more vectors it can use,
3073 * the better it will generally perform. To keep memory demands in check, we
3074 * set the number of vectors to 100. This means that up to 100 solver
3075 * iterations, every temporary vector can be stored. If the solver needs to
3076 * iterate more often to get the specified tolerance, it will work on a
3077 * reduced set of vectors by restarting at every 100 iterations.
3081 * With this all set up, we solve the system and distribute the constraints
3082 * in the Stokes system, i.e., hanging nodes and no-flux boundary condition,
3083 * in order to have the appropriate solution values even at constrained
3084 * dofs. Finally, we write the number of iterations to the screen.
3087 * template <int dim>
3088 * void BoussinesqFlowProblem<dim>::solve()
3090 * std::cout << " Solving..." << std::endl;
3093 * const LinearSolvers::InverseMatrix<TrilinosWrappers::SparseMatrix,
3094 * TrilinosWrappers::PreconditionIC>
3095 * mp_inverse(stokes_preconditioner_matrix.block(1, 1),
3096 * *Mp_preconditioner);
3098 * const LinearSolvers::BlockSchurPreconditioner<
3099 * TrilinosWrappers::PreconditionAMG,
3100 * TrilinosWrappers::PreconditionIC>
3101 * preconditioner(stokes_matrix, mp_inverse, *Amg_preconditioner);
3103 * SolverControl solver_control(stokes_matrix.m(),
3104 * 1e-6 * stokes_rhs.l2_norm());
3106 * SolverGMRES<TrilinosWrappers::MPI::BlockVector> gmres(
3108 * SolverGMRES<TrilinosWrappers::MPI::BlockVector>::AdditionalData(100));
3110 * for (unsigned int i = 0; i < stokes_solution.size(); ++i)
3111 * if (stokes_constraints.is_constrained(i))
3112 * stokes_solution(i) = 0;
3114 * gmres.solve(stokes_matrix, stokes_solution, stokes_rhs, preconditioner);
3116 * stokes_constraints.distribute(stokes_solution);
3118 * std::cout << " " << solver_control.last_step()
3119 * << " GMRES iterations for Stokes subsystem." << std::endl;
3124 * Once we know the Stokes solution, we can determine the new time step
3125 * from the maximal velocity. We have to do this to satisfy the CFL
3126 * condition since convection terms are treated explicitly in the
3127 * temperature equation, as discussed in the introduction. The exact form
3128 * of the formula used here for the time step is discussed in the results
3129 * section of this program.
3133 * There is a snatch here. The formula contains a division by the maximum
3134 * value of the velocity. However, at the start of the computation, we
3135 * have a constant temperature field (we start with a constant
3136 * temperature, and it will be nonconstant only after the first time step
3137 * during which the source acts). Constant temperature means that no
3138 * buoyancy acts, and so the velocity is zero. Dividing by it will not
3139 * likely lead to anything good.
3143 * To avoid the resulting infinite time step, we ask whether the maximal
3144 * velocity is very small (in particular smaller than the values we
3145 * encounter during any of the following time steps) and if so rather than
3146 * dividing by zero we just divide by a small value, resulting in a large
3147 * but finite time step.
3150 * old_time_step = time_step;
3151 * const double maximal_velocity = get_maximal_velocity();
3153 * if (maximal_velocity >= 0.01)
3154 * time_step = 1. / (1.7 * dim * std::sqrt(1. * dim)) / temperature_degree *
3155 * GridTools::minimal_cell_diameter(triangulation) /
3158 * time_step = 1. / (1.7 * dim * std::sqrt(1. * dim)) / temperature_degree *
3159 * GridTools::minimal_cell_diameter(triangulation) / .01;
3162 * << "Time step: " << time_step << std::endl;
3164 * temperature_solution = old_temperature_solution;
3168 * Next we set up the temperature system and the right hand side using the
3169 * function <code>assemble_temperature_system()</code>. Knowing the
3170 * matrix and right hand side of the temperature equation, we set up a
3171 * preconditioner and a solver. The temperature matrix is a mass matrix
3172 * (with eigenvalues around one) plus a Laplace matrix (with eigenvalues
3173 * between zero and @f$ch^{-2}@f$) times a small number proportional to the
3174 * time step @f$k_n@f$. Hence, the resulting symmetric and positive definite
3175 * matrix has eigenvalues in the range @f$[1,1+k_nh^{-2}]@f$ (up to
3176 * constants). This matrix is only moderately ill conditioned even for
3177 * small mesh sizes and we get a reasonably good preconditioner by simple
3178 * means, for example with an incomplete Cholesky decomposition
3179 * preconditioner (IC) as we also use for preconditioning the pressure
3180 * mass matrix solver. As a solver, we choose the conjugate gradient
3181 * method CG. As before, we tell the solver to use Trilinos vectors via
3182 * the template argument <code>TrilinosWrappers::MPI::Vector</code>.
3183 * Finally, we solve, distribute the hanging node constraints and write out
3184 * the number of iterations.
3187 * assemble_temperature_system(maximal_velocity);
3189 * SolverControl solver_control(temperature_matrix.m(),
3190 * 1e-8 * temperature_rhs.l2_norm());
3191 * SolverCG<TrilinosWrappers::MPI::Vector> cg(solver_control);
3193 * TrilinosWrappers::PreconditionIC preconditioner;
3194 * preconditioner.initialize(temperature_matrix);
3196 * cg.solve(temperature_matrix,
3197 * temperature_solution,
3201 * temperature_constraints.distribute(temperature_solution);
3203 * std::cout << " " << solver_control.last_step()
3204 * << " CG iterations for temperature." << std::endl;
3208 * At the end of this function, we step through the vector and read out
3209 * the maximum and minimum temperature value, which we also want to
3210 * output. This will come in handy when determining the correct constant
3211 * in the choice of time step as discuss in the results section of this
3215 * double min_temperature = temperature_solution(0),
3216 * max_temperature = temperature_solution(0);
3217 * for (unsigned int i = 0; i < temperature_solution.size(); ++i)
3220 * std::min<double>(min_temperature, temperature_solution(i));
3222 * std::max<double>(max_temperature, temperature_solution(i));
3225 * std::cout << " Temperature range: " << min_temperature << ' '
3226 * << max_temperature << std::endl;
3235 * <a name="step_31-BoussinesqFlowProblemoutput_results"></a>
3236 * <h4>BoussinesqFlowProblem::output_results</h4>
3240 * This function writes the solution to a VTK output file for visualization,
3241 * which is done every tenth time step. This is usually quite a simple task,
3242 * since the deal.II library provides functions that do almost all the job
3243 * for us. There is one new function compared to previous examples: We want
3244 * to visualize both the Stokes solution and the temperature as one data
3245 * set, but we have done all the calculations based on two different
3246 * DoFHandler objects. Luckily, the DataOut class is prepared to deal with
3247 * it. All we have to do is to not attach one single DoFHandler at the
3248 * beginning and then use that for all added vector, but specify the
3249 * DoFHandler to each vector separately. The rest is done as in @ref step_22 "step-22". We
3250 * create solution names (that are going to appear in the visualization
3251 * program for the individual components). The first <code>dim</code>
3252 * components are the vector velocity, and then we have pressure for the
3253 * Stokes part, whereas temperature is scalar. This information is read out
3254 * using the DataComponentInterpretation helper class. Next, we actually
3255 * attach the data vectors with their DoFHandler objects, build patches
3256 * according to the degree of freedom, which are (sub-) elements that
3257 * describe the data for visualization programs. Finally, we open a file
3258 * (that includes the time step number) and write the vtk data into it.
3261 * template <int dim>
3262 * void BoussinesqFlowProblem<dim>::output_results() const
3264 * if (timestep_number % 10 != 0)
3267 * std::vector<std::string> stokes_names(dim, "velocity");
3268 * stokes_names.emplace_back("p");
3269 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
3270 * stokes_component_interpretation(
3271 * dim + 1, DataComponentInterpretation::component_is_scalar);
3272 * for (unsigned int i = 0; i < dim; ++i)
3273 * stokes_component_interpretation[i] =
3274 * DataComponentInterpretation::component_is_part_of_vector;
3276 * DataOut<dim> data_out;
3277 * data_out.add_data_vector(stokes_dof_handler,
3280 * stokes_component_interpretation);
3281 * data_out.add_data_vector(temperature_dof_handler,
3282 * temperature_solution,
3284 * data_out.build_patches(std::min(stokes_degree, temperature_degree));
3286 * std::ofstream output("solution-" +
3287 * Utilities::int_to_string(timestep_number, 4) + ".vtk");
3288 * data_out.write_vtk(output);
3296 * <a name="step_31-BoussinesqFlowProblemrefine_mesh"></a>
3297 * <h4>BoussinesqFlowProblem::refine_mesh</h4>
3301 * This function takes care of the adaptive mesh refinement. The three tasks
3302 * this function performs is to first find out which cells to
3303 * refine/coarsen, then to actually do the refinement and eventually
3304 * transfer the solution vectors between the two different grids. The first
3305 * task is simply achieved by using the well-established Kelly error
3306 * estimator on the temperature (it is the temperature we're mainly
3358 *
cell->clear_refine_flag();
3383 *
const std::vector<TrilinosWrappers::MPI::Vector>
x_temperature = {
3455 * <a name=
"step_31-BoussinesqFlowProblemrun"></a>
3456 * <
h4>BoussinesqFlowProblem::run</
h4>
3477 * initialize time step number
and time step
and start
the time
loop.
3483 * const unsigned int initial_refinement = (dim == 2 ? 4 : 2);
3484 * const unsigned int n_pre_refinement_steps = (dim == 2 ? 4 : 3);
3487 * GridGenerator::hyper_cube(triangulation);
3488 * global_Omega_diameter = GridTools::diameter(triangulation);
3490 * triangulation.refine_global(initial_refinement);
3494 * unsigned int pre_refinement_step = 0;
3496 * start_time_iteration:
3498 * VectorTools::project(temperature_dof_handler,
3499 * temperature_constraints,
3500 * QGauss<dim>(temperature_degree + 2),
3501 * EquationData::TemperatureInitialValues<dim>(),
3502 * old_temperature_solution);
3504 * timestep_number = 0;
3505 * time_step = old_time_step = 0;
3511 * std::cout <<
"Timestep " << timestep_number <<
": t=" << time
3516 * The first steps in the time loop are all obvious – we
3517 * assemble the Stokes system, the preconditioner, the temperature
3518 * matrix (matrices and preconditioner do actually only change in case
3519 * we
've remeshed before), and then do the solve. Before going on with
3520 * the next time step, we have to check whether we should first finish
3521 * the pre-refinement steps or if we should remesh (every fifth time
3522 * step), refining up to a level that is consistent with initial
3523 * refinement and pre-refinement steps. Last in the loop is to advance
3524 * the solutions, i.e., to copy the solutions to the next "older" time
3528 * assemble_stokes_system();
3529 * build_stokes_preconditioner();
3530 * assemble_temperature_matrix();
3536 * std::cout << std::endl;
3538 * if ((timestep_number == 0) &&
3539 * (pre_refinement_step < n_pre_refinement_steps))
3541 * refine_mesh(initial_refinement + n_pre_refinement_steps);
3542 * ++pre_refinement_step;
3543 * goto start_time_iteration;
3545 * else if ((timestep_number > 0) && (timestep_number % 5 == 0))
3546 * refine_mesh(initial_refinement + n_pre_refinement_steps);
3548 * time += time_step;
3549 * ++timestep_number;
3551 * old_stokes_solution = stokes_solution;
3552 * old_old_temperature_solution = old_temperature_solution;
3553 * old_temperature_solution = temperature_solution;
3557 * Do all the above until we arrive at time 100.
3560 * while (time <= 100);
3562 * } // namespace Step31
3569 * <a name="step_31-Thecodemaincodefunction"></a>
3570 * <h3>The <code>main</code> function</h3>
3574 * The main function looks almost the same as in all other programs.
3578 * There is one difference we have to be careful about. This program uses
3579 * Trilinos and, typically, Trilinos is configured so that it can run in
3580 * %parallel using MPI. This doesn't mean that it <i>has</i> to run in
3581 * %parallel, and in fact this program (unlike @ref step_32
"step-32") makes no attempt at
3582 * all to do anything in %parallel using MPI. Nevertheless, Trilinos wants the
3583 * MPI system to be initialized. We do that be creating an object of type
3584 * Utilities::MPI::MPI_InitFinalize that initializes MPI (if available) using
3585 * the arguments given to main() (i.e., <code>argc</code> and
3586 * <code>argv</code>) and de-initializes it again when the object goes out of
3590 * int main(int argc, char *argv[])
3594 * using namespace dealii;
3595 * using namespace Step31;
3597 * Utilities::MPI::MPI_InitFinalize mpi_initialization(
3598 * argc, argv, numbers::invalid_unsigned_int);
3602 * This program can only be run in serial. Otherwise, throw an exception.
3605 * AssertThrow(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) == 1,
3607 *
"This program can only be run in serial, use ./step-31"));
3609 * BoussinesqFlowProblem<2> flow_problem;
3610 * flow_problem.run();
3612 * catch (std::exception &exc)
3614 * std::cerr << std::endl
3616 * <<
"----------------------------------------------------"
3618 * std::cerr <<
"Exception on processing: " << std::endl
3619 * << exc.what() << std::endl
3620 * <<
"Aborting!" << std::endl
3621 * <<
"----------------------------------------------------"
3628 * std::cerr << std::endl
3630 * <<
"----------------------------------------------------"
3632 * std::cerr <<
"Unknown exception!" << std::endl
3633 * <<
"Aborting!" << std::endl
3634 * <<
"----------------------------------------------------"
3642<a name=
"step_31-Results"></a><h1>Results</h1>
3645<a name=
"step_31-Resultsin2d"></a><h3> Results in 2d </h3>
3648When you run the program in 2d, the output will look something like
3652Number of active cells: 256 (on 5 levels)
3653Number of degrees of freedom: 3556 (2178+289+1089)
3657 Rebuilding Stokes preconditioner...
3659 0 GMRES iterations for Stokes subsystem.
3661 9 CG iterations for temperature.
3662 Temperature range: -0.16687 1.30011
3664Number of active cells: 280 (on 6 levels)
3665Number of degrees of freedom: 4062 (2490+327+1245)
3669 Rebuilding Stokes preconditioner...
3671 0 GMRES iterations for Stokes subsystem.
3673 9 CG iterations for temperature.
3674 Temperature range: -0.0982971 0.598503
3676Number of active cells: 520 (on 7 levels)
3677Number of degrees of freedom: 7432 (4562+589+2281)
3681 Rebuilding Stokes preconditioner...
3683 0 GMRES iterations for Stokes subsystem.
3685 9 CG iterations for temperature.
3686 Temperature range: -0.0551098 0.294493
3688Number of active cells: 1072 (on 8 levels)
3689Number of degrees of freedom: 15294 (9398+1197+4699)
3693 Rebuilding Stokes preconditioner...
3695 0 GMRES iterations for Stokes subsystem.
3697 9 CG iterations for temperature.
3698 Temperature range: -0.0273524 0.156861
3700Number of active cells: 2116 (on 9 levels)
3701Number of degrees of freedom: 30114 (18518+2337+9259)
3705 Rebuilding Stokes preconditioner...
3707 0 GMRES iterations for Stokes subsystem.
3708 Time step: 0.0574449
3709 9 CG iterations for temperature.
3710 Temperature range: -0.014993 0.0738328
3712Timestep 1: t=0.0574449
3715 56 GMRES iterations for Stokes subsystem.
3716 Time step: 0.0574449
3717 9 CG iterations for temperature.
3718 Temperature range: -0.0273934 0.14488
3724In the beginning we refine the mesh several times adaptively and
3725always return to time step zero to restart on the newly refined
3726mesh. Only then do we start the actual time iteration.
3728The program runs for a while. The temperature field for time steps 0,
3729500, 1000, 1500, 2000, 3000, 4000, and 5000 looks like this (note that
3730the color scale used for the temperature is not always the same):
3732<table align=
"center" class=
"doxtable">
3735 <img src=
"https://www.dealii.org/images/steps/developer/step-31.2d.solution.00.png" alt=
"">
3738 <img src=
"https://www.dealii.org/images/steps/developer/step-31.2d.solution.01.png" alt=
"">
3741 <img src=
"https://www.dealii.org/images/steps/developer/step-31.2d.solution.02.png" alt=
"">
3744 <img src=
"https://www.dealii.org/images/steps/developer/step-31.2d.solution.03.png" alt=
"">
3749 <img src=
"https://www.dealii.org/images/steps/developer/step-31.2d.solution.04.png" alt=
"">
3752 <img src=
"https://www.dealii.org/images/steps/developer/step-31.2d.solution.05.png" alt=
"">
3755 <img src=
"https://www.dealii.org/images/steps/developer/step-31.2d.solution.06.png" alt=
"">
3758 <img src=
"https://www.dealii.org/images/steps/developer/step-31.2d.solution.07.png" alt=
"">
3763The visualizations shown here were generated using a version of the example
3764which did not enforce the constraints after transferring the mesh.
3766As can be seen, we have three heat sources that heat fluid and
3767therefore produce a buoyancy effect that lets hots pockets of fluid
3768rise up and swirl around. By a chimney effect, the three streams are
3769pressed together by fluid that comes from the outside and wants to
3770join the updraft party. Note that because the fluid is initially at
3771rest, those parts of the fluid that were initially over the sources
3772receive a longer heating time than that fluid that is later dragged
3773over the source by the fully developed flow field. It is therefore
3774hotter, a fact that can be seen in the red tips of the three
3775plumes. Note also the relatively fine features of the flow field, a
3776result of the sophisticated transport stabilization of the temperature
3777equation we have chosen.
3779In addition to the pictures above, the following ones show the
3780adaptive mesh and the flow field at the same time steps:
3782<table align=
"center" class=
"doxtable">
3785 <img src=
"https://www.dealii.org/images/steps/developer/step-31.2d.grid.00.png" alt=
"">
3788 <img src=
"https://www.dealii.org/images/steps/developer/step-31.2d.grid.01.png" alt=
"">
3791 <img src=
"https://www.dealii.org/images/steps/developer/step-31.2d.grid.02.png" alt=
"">
3794 <img src=
"https://www.dealii.org/images/steps/developer/step-31.2d.grid.03.png" alt=
"">
3799 <img src=
"https://www.dealii.org/images/steps/developer/step-31.2d.grid.04.png" alt=
"">
3802 <img src=
"https://www.dealii.org/images/steps/developer/step-31.2d.grid.05.png" alt=
"">
3805 <img src=
"https://www.dealii.org/images/steps/developer/step-31.2d.grid.06.png" alt=
"">
3808 <img src=
"https://www.dealii.org/images/steps/developer/step-31.2d.grid.07.png" alt=
"">
3814<a name=
"step_31-Resultsin3d"></a><h3> Results in 3d </h3>
3817The same thing can of course be done in 3d by changing the template
3818parameter to the BoussinesqFlowProblem object in <code>main()</code>
3819from 2 to 3, so that the output now looks like follows:
3823Number of active cells: 64 (on 3 levels)
3824Number of degrees of freedom: 3041 (2187+125+729)
3828 Rebuilding Stokes preconditioner...
3830 0 GMRES iterations for Stokes subsystem.
3832 9 CG iterations for temperature.
3833 Temperature range: -0.675683 4.94725
3835Number of active cells: 288 (on 4 levels)
3836Number of degrees of freedom: 12379 (8943+455+2981)
3840 Rebuilding Stokes preconditioner...
3842 0 GMRES iterations for Stokes subsystem.
3844 9 CG iterations for temperature.
3845 Temperature range: -0.527701 2.25764
3847Number of active cells: 1296 (on 5 levels)
3848Number of degrees of freedom: 51497 (37305+1757+12435)
3852 Rebuilding Stokes preconditioner...
3854 0 GMRES iterations for Stokes subsystem.
3856 10 CG iterations for temperature.
3857 Temperature range: -0.496942 0.847395
3859Number of active cells: 5048 (on 6 levels)
3860Number of degrees of freedom: 192425 (139569+6333+46523)
3864 Rebuilding Stokes preconditioner...
3866 0 GMRES iterations for Stokes subsystem.
3868 10 CG iterations for temperature.
3869 Temperature range: -0.267683 0.497739
3871Timestep 1: t=0.306373
3874 27 GMRES iterations for Stokes subsystem.
3876 10 CG iterations for temperature.
3877 Temperature range: -0.461787 0.958679
3883Visualizing the temperature isocontours at time steps 0,
388450, 100, 150, 200, 300, 400, 500, 600, 700, and 800 yields the
3887<table align=
"center" class=
"doxtable">
3890 <img src=
"https://www.dealii.org/images/steps/developer/step-31.3d.solution.00.png" alt=
"">
3893 <img src=
"https://www.dealii.org/images/steps/developer/step-31.3d.solution.01.png" alt=
"">
3896 <img src=
"https://www.dealii.org/images/steps/developer/step-31.3d.solution.02.png" alt=
"">
3899 <img src=
"https://www.dealii.org/images/steps/developer/step-31.3d.solution.03.png" alt=
"">
3904 <img src=
"https://www.dealii.org/images/steps/developer/step-31.3d.solution.04.png" alt=
"">
3907 <img src=
"https://www.dealii.org/images/steps/developer/step-31.3d.solution.05.png" alt=
"">
3910 <img src=
"https://www.dealii.org/images/steps/developer/step-31.3d.solution.06.png" alt=
"">
3913 <img src=
"https://www.dealii.org/images/steps/developer/step-31.3d.solution.07.png" alt=
"">
3918 <img src=
"https://www.dealii.org/images/steps/developer/step-31.3d.solution.08.png" alt=
"">
3921 <img src=
"https://www.dealii.org/images/steps/developer/step-31.3d.solution.09.png" alt=
"">
3924 <img src=
"https://www.dealii.org/images/steps/developer/step-31.3d.solution.10.png" alt=
"">
3931That the first picture looks like three hedgehogs stems from the fact that our
3932scheme essentially projects the source times the first time step size onto the
3933mesh to obtain the temperature field in the first time step. Since the source
3934function is discontinuous, we need to expect over- and undershoots from this
3935project. This is in fact what happens (it
's easier to check this in 2d) and
3936leads to the crumpled appearance of the isosurfaces. The visualizations shown
3937here were generated using a version of the example which did not enforce the
3938constraints after transferring the mesh.
3942<a name="step_31-Numericalexperimentstodetermineoptimalparameters"></a><h3> Numerical experiments to determine optimal parameters </h3>
3945The program as is has three parameters that we don't have much of a
3946theoretical handle on how to choose in an optimal way. These are:
3948 <li>The time step must satisfy a CFL condition
3949 @f$k\le \min_K \frac{c_kh_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$. Here, @f$c_k@f$ is
3950 dimensionless, but what is the right value?
3951 <li>In the computation of the artificial viscosity,
3956 \|\mathbf{u}\|_{L^\infty(K)}
3960 \frac{\|R_\alpha(T)\|_{L^\infty(K)}}{c(\mathbf{u},T)}
3963 with @f$c(\mathbf{u},T) =
3964 c_R\ \|\mathbf{u}\|_{L^\infty(\Omega)} \ \mathrm{var}(T)
3965 \ |\mathrm{diam}(\Omega)|^{\alpha-2}@f$.
3966 Here, the choice of the dimensionless %numbers @f$\beta,c_R@f$ is of
3969In all of these cases, we will have to expect that the correct choice of each
3970value depends on that of the others, and most likely also on the space
3971dimension and polynomial degree of the finite element used for the
3972temperature. Below we
'll discuss a few numerical experiments to choose
3973constants @f$c_k@f$ and @f$\beta@f$.
3975Below, we will not discuss the choice of @f$c_R@f$. In the program, we set
3976it to @f$c_R=2^{\frac{4-2\alpha}{d}}@f$. The reason for this value is a
3977bit complicated and has more to do with the history of the program
3978than reasoning: while the correct formula for the global scaling
3979parameter @f$c(\mathbf{u},T)@f$ is shown above, the program (including the
3980version shipped with deal.II 6.2) initially had a bug in that we
3983 \|\mathbf{u}\|_{L^\infty(\Omega)} \ \mathrm{var}(T)
3984 \ \frac{1}{|\mathrm{diam}(\Omega)|^{\alpha-2}}@f$ instead, where
3985we had set the scaling parameter to one. Since we only computed on the
3986unit square/cube where @f$\mathrm{diam}(\Omega)=2^{1/d}@f$, this was
3987entirely equivalent to using the correct formula with
3988@f$c_R=\left(2^{1/d}\right)^{4-2\alpha}=2^{\frac{4-2\alpha}{d}}@f$. Since
3989this value for @f$c_R@f$ appears to work just fine for the current
3990program, we corrected the formula in the program and set @f$c_R@f$ to a
3991value that reproduces exactly the results we had before. We will,
3992however, revisit this issue again in @ref step_32 "step-32".
3994Now, however, back to the discussion of what values of @f$c_k@f$ and
3995@f$\beta@f$ to choose:
3998<a name="step_31-Choosingicsubksubiandbeta"></a><h4> Choosing <i>c<sub>k</sub></i> and beta </h4>
4001These two constants are definitely linked in some way. The reason is easy to
4002see: In the case of a pure advection problem,
4003@f$\frac{\partial T}{\partial t} + \mathbf{u}\cdot\nabla T = \gamma@f$, any
4004explicit scheme has to satisfy a CFL condition of the form
4005@f$k\le \min_K \frac{c_k^a h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$. On the other hand,
4006for a pure diffusion problem,
4007@f$\frac{\partial T}{\partial t} + \nu \Delta T = \gamma@f$,
4008explicit schemes need to satisfy a condition
4009@f$k\le \min_K \frac{c_k^d h_K^2}{\nu}@f$. So given the form of @f$\nu@f$ above, an
4010advection diffusion problem like the one we have to solve here will result in
4011a condition of the form
4013k\le \min_K \min \left\{
4014 \frac{c_k^a h_K}{\|\mathbf{u}\|_{L^\infty(K)}},
4015 \frac{c_k^d h_K^2}{\beta \|\mathbf{u}\|_{L^\infty(K)} h_K}\right\}
4017 \min_K \left( \min \left\{
4019 \frac{c_k^d}{\beta}\right\}
4020 \frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}} \right)
4022It follows that we have to face the fact that we might want to choose @f$\beta@f$
4023larger to improve the stability of the numerical scheme (by increasing the
4024amount of artificial diffusion), but we have to pay a price in the form of
4025smaller, and consequently more time steps. In practice, one would therefore
4026like to choose @f$\beta@f$ as small as possible to keep the transport problem
4027sufficiently stabilized while at the same time trying to choose the time step
4028as large as possible to reduce the overall amount of work.
4030The find the right balance, the only way is to do a few computational
4031experiments. Here's what we did: We modified the program slightly to allow
4032less mesh refinement (so we don
't always have to wait that long) and to choose
4037 \|\mathbf{u}\|_{L^\infty(K)} h_K
4038@f$ to eliminate the effect of the constant @f$c_R@f$ (we know that
4039solutions are stable by using this version of @f$\nu(T)@f$ as an artificial
4040viscosity, but that we can improve things -- i.e. make the solution
4041sharper -- by using the more complicated formula for this artificial
4042viscosity). We then run the program
4043for different values @f$c_k,\beta@f$ and observe maximal and minimal temperatures
4044in the domain. What we expect to see is this: If we choose the time step too
4045big (i.e. choose a @f$c_k@f$ bigger than theoretically allowed) then we will get
4046exponential growth of the temperature. If we choose @f$\beta@f$ too small, then
4047the transport stabilization becomes insufficient and the solution will show
4048significant oscillations but not exponential growth.
4051<a name="step_31-ResultsforQsub1subelements"></a><h5>Results for Q<sub>1</sub> elements</h5>
4054Here is what we get for
4055@f$\beta=0.01, \beta=0.1@f$, and @f$\beta=0.5@f$, different choices of @f$c_k@f$, and
4056bilinear elements (<code>temperature_degree=1</code>) in 2d:
4058<table align="center" class="doxtable">
4061 <img src="https://www.dealii.org/images/steps/developer/step-31.timestep.q1.beta=0.01.png" alt="">
4064 <img src="https://www.dealii.org/images/steps/developer/step-31.timestep.q1.beta=0.03.png" alt="">
4069 <img src="https://www.dealii.org/images/steps/developer/step-31.timestep.q1.beta=0.1.png" alt="">
4072 <img src="https://www.dealii.org/images/steps/developer/step-31.timestep.q1.beta=0.5.png" alt="">
4077The way to interpret these graphs goes like this: for @f$\beta=0.01@f$ and
4078@f$c_k=\frac 12,\frac 14@f$, we see exponential growth or at least large
4079variations, but if we choose
4080@f$k=\frac 18\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$
4081or smaller, then the scheme is
4082stable though a bit wobbly. For more artificial diffusion, we can choose
4083@f$k=\frac 14\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$
4084or smaller for @f$\beta=0.03@f$,
4085@f$k=\frac 13\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$
4086or smaller for @f$\beta=0.1@f$, and again need
4087@f$k=\frac 1{15}\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$
4088for @f$\beta=0.5@f$ (this time because much diffusion requires a small time
4091So how to choose? If we were simply interested in a large time step, then we
4092would go with @f$\beta=0.1@f$ and
4093@f$k=\frac 13\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$.
4094On the other hand, we're also interested in accuracy and here it may be of
4095interest to actually investigate what these curves show. To this end note that
4096we start with a zero temperature and that our sources are positive — so
4097we would intuitively expect that the temperature can never drop below
4098zero. But it does, a consequence of Gibb
's phenomenon when using continuous
4099elements to approximate a discontinuous solution. We can therefore see that
4100choosing @f$\beta@f$ too small is bad: too little artificial diffusion leads to
4101over- and undershoots that aren't diffused away. On the other hand, for large
4102@f$\beta@f$, the minimum temperature drops below zero at the beginning but then
4103quickly diffuses back to zero.
4105On the other hand, let
's also look at the maximum temperature. Watching the
4106movie of the solution, we see that initially the fluid is at rest. The source
4107keeps heating the same volume of fluid whose temperature increases linearly at
4108the beginning until its buoyancy is able to move it upwards. The hottest part
4109of the fluid is therefore transported away from the solution and fluid taking
4110its place is heated for only a short time before being moved out of the source
4111region, therefore remaining cooler than the initial bubble. If @f$\kappa=0@f$
4112(in the program it is nonzero but very small) then the hottest part of the
4113fluid should be advected along with the flow with its temperature
4114constant. That's what we can see in the graphs with the smallest @f$\beta@f$: Once
4115the maximum temperature is reached, it hardly changes any more. On the other
4116hand, the larger the artificial diffusion, the more the hot spot is
4117diffused. Note that for this criterion, the time step size does not play a
4120So to sum up, likely the best choice would appear to be @f$\beta=0.03@f$
4121and @f$k=\frac 14\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$. The curve is
4122a bit wobbly, but overall pictures looks pretty reasonable with the
4123exception of some over and undershoots close to the start time due to
4127<a name="step_31-ResultsforQsub2subelements"></a><h5>Results for Q<sub>2</sub> elements</h5>
4130One can repeat the same sequence of experiments for higher order
4131elements as well. Here are the graphs for bi-quadratic shape functions
4132(<code>temperature_degree=2</code>) for the temperature, while we
4133retain the @f$Q_2/Q_1@f$ stable Taylor-Hood element for the Stokes system:
4135<table align="center" class="doxtable">
4138 <img src="https://www.dealii.org/images/steps/developer/step-31.timestep.q2.beta=0.01.png" alt="">
4141 <img src="https://www.dealii.org/images/steps/developer/step-31.timestep.q2.beta=0.03.png" alt="">
4146 <img src="https://www.dealii.org/images/steps/developer/step-31.timestep.q2.beta=0.1.png" alt="">
4151Again, small values of @f$\beta@f$ lead to less diffusion but we have to
4152choose the time step very small to keep things under control. Too
4153large values of @f$\beta@f$ make for more diffusion, but again require
4154small time steps. The best value would appear to be @f$\beta=0.03@f$, as
4155for the @f$Q_1@f$ element, and then we have to choose
4156@f$k=\frac 18\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$ — exactly
4157half the size for the @f$Q_1@f$ element, a fact that may not be surprising
4158if we state the CFL condition as the requirement that the time step be
4159small enough so that the distance transport advects in each time step
4160is no longer than one <i>grid point</i> away (which for @f$Q_1@f$ elements
4161is @f$h_K@f$, but for @f$Q_2@f$ elements is @f$h_K/2@f$). It turns out that @f$\beta@f$
4162needs to be slightly larger for obtaining stable results also late in
4163the simulation at times larger than 60, so we actually choose it as
4164@f$\beta = 0.034@f$ in the code.
4167<a name="step_31-Resultsfor3d"></a><h5>Results for 3d</h5>
4170One can repeat these experiments in 3d and find the optimal time step
4171for each value of @f$\beta@f$ and find the best value of @f$\beta@f$. What one
4172finds is that for the same @f$\beta@f$ already used in 2d, the time steps
4173needs to be a bit smaller, by around a factor of 1.2 or so. This is
4174easily explained: the time step restriction is
4175@f$k=\min_K \frac{ch_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$ where @f$h_K@f$ is
4176the <i>diameter</i> of the cell. However, what is really needed is the
4177distance between mesh points, which is @f$\frac{h_K}{\sqrt{d}}@f$. So a
4178more appropriate form would be
4179@f$k=\min_K \frac{ch_K}{\|\mathbf{u}\|_{L^\infty(K)}\sqrt{d}}@f$.
4181The second find is that one needs to choose @f$\beta@f$ slightly bigger
4182(about @f$\beta=0.05@f$ or so). This then again reduces the time step we
4188<a name="step_31-Conclusions"></a><h5>Conclusions</h5>
4191Concluding, from the simple computations above, @f$\beta=0.034@f$ appears to be a
4192good choice for the stabilization parameter in 2d, and @f$\beta=0.05@f$ in 3d. In
4193a dimension independent way, we can model this as @f$\beta=0.017d@f$. If one does
4194longer computations (several thousand time steps) on finer meshes, one
4195realizes that the time step size is not quite small enough and that for
4196stability one will have to reduce the above values a bit more (by about a
4197factor of @f$\frac 78@f$).
4199As a consequence, a formula that reconciles 2d, 3d, and variable polynomial
4200degree and takes all factors in account reads as follows:
4203 \frac 1{2 \cdot 1.7} \frac 1{\sqrt{d}}
4206 \frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}
4208 \frac 1{1.7 d\sqrt{d}}
4210 \frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}.
4212In the first form (in the center of the equation), @f$\frac
42131{2 \cdot 1.7}@f$ is a universal constant, @f$\frac 1{\sqrt{d}}@f$
4214is the factor that accounts for the difference between cell diameter
4215and grid point separation,
4216@f$\frac 2d@f$ accounts for the increase in @f$\beta@f$ with space dimension,
4217@f$\frac 1{q_T}@f$ accounts for the distance between grid points for
4218higher order elements, and @f$\frac{h_K}{\|\mathbf{u}\|_{L^\infty(K)}}@f$
4219for the local speed of transport relative to the cell size. This is
4220the formula that we use in the program.
4222As for the question of whether to use @f$Q_1@f$ or @f$Q_2@f$ elements for the
4223temperature, the following considerations may be useful: First,
4224solving the temperature equation is hardly a factor in the overall
4225scheme since almost the entire compute time goes into solving the
4226Stokes system in each time step. Higher order elements for the
4227temperature equation are therefore not a significant drawback. On the
4228other hand, if one compares the size of the over- and undershoots the
4229solution produces due to the discontinuous source description, one
4230notices that for the choice of @f$\beta@f$ and @f$k@f$ as above, the @f$Q_1@f$
4231solution dips down to around @f$-0.47@f$, whereas the @f$Q_2@f$ solution only
4232goes to @f$-0.13@f$ (remember that the exact solution should never become
4233negative at all. This means that the @f$Q_2@f$ solution is significantly
4234more accurate; the program therefore uses these higher order elements,
4235despite the penalty we pay in terms of smaller time steps.
4238<a name="step_31-Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
4241There are various ways to extend the current program. Of particular interest
4242is, of course, to make it faster and/or increase the resolution of the
4243program, in particular in 3d. This is the topic of the @ref step_32 "step-32"
4244tutorial program which will implement strategies to solve this problem in
4245%parallel on a cluster. It is also the basis of the much larger open
4246source code ASPECT (see https://aspect.geodynamics.org/ ) that can solve realistic
4247problems and that constitutes the further development of @ref step_32 "step-32".
4249Another direction would be to make the fluid flow more realistic. The program
4250was initially written to simulate various cases simulating the convection of
4251material in the earth's mantle, i.e. the zone between the outer earth core and
4252the solid earth crust: there, material is heated from below and cooled from
4253above, leading to thermal convection. The physics of this fluid are much more
4254complicated than shown in this program, however: The viscosity of mantle
4255material is strongly dependent on the temperature, i.e. @f$\eta=\eta(T)@f$, with
4256the dependency frequently modeled as a viscosity that is reduced exponentially
4257with rising temperature. Secondly, much of the dynamics of the mantle is
4258determined by chemical reactions, primarily phase changes of the various
4259crystals that make up the mantle; the buoyancy term on the right hand side of
4260the Stokes equations then depends not only on the temperature, but also on the
4261chemical composition at a given location which is advected by the flow field
4262but also changes as a function of pressure and temperature. We will
4263investigate some of these effects in later tutorial programs as well.
4266<a name=
"step_31-PlainProg"></a>
4267<h1> The plain program</h1>
4268@include
"step-31.cc"
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
const std::vector< Point< dim > > & get_unit_support_points() const
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
numbers::NumberTraits< Number >::real_type norm() const
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 & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
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_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
IndexSet complete_index_set(const IndexSet::size_type N)
std::vector< index_type > data
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
Expression sign(const Expression &x)
void downstream(DoFHandler< dim, spacedim > &dof_handler, const Tensor< 1, spacedim > &direction, const bool dof_wise_renumbering=false)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void refine_and_coarsen_fixed_fraction(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double top_fraction, const double bottom_fraction, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max(), const VectorTools::NormType norm_type=VectorTools::L1_norm)
@ matrix
Contents is actually a matrix.
constexpr types::blas_int zero
constexpr types::blas_int one
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
VectorType::value_type * end(VectorType &V)
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)
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation