problem. In the following table, the residuals at each Newton's iteration on
every mesh is shown. The data in the table shows that Newton's iteration converges quadratically.
-<table align="center" border="1">
- <tr>
- <th> </th>
+<table align="center" class="doxtable">
+<tr>
+ <th>Re=400</th>
<th colspan="2">Mesh0</th>
<th colspan="2">Mesh1</th>
<th colspan="2">Mesh2</th>
<th colspan="2">Mesh3</th>
<th colspan="2">Mesh4</th>
- </tr>
- <tr>
+</tr>
+<tr>
<th>Newton iter </th>
<th>Residual </th>
<th>FGMRES </th>
<th>FGMRES </th>
<th>Residual </th>
<th>FGMRES </th>
- </tr>
- <tr>
- <td>1 </td>
- <td>7.40396e-3</td>
- <td>3 </td>
- <td>1.05562e-3 </td>
- <td>3 </td>
- <td>4.94796e-4 </td>
- <td>3 </td>
- <td>2.5624e-4 </td>
- <td>2 </td>
- <td>1.26733e-4 </td>
- <td>2 </td>
- </tr>
- <tr>
- <td>2 </td>
- <td>3.86766e-3 </td>
- <td>4 </td>
- <td>1.3549e-5 </td>
- <td>3 </td>
- <td>1.41981e-6 </td>
- <td>3 </td>
- <td>1.29108e-6 </td>
- <td>4 </td>
- <td>6.14794e-7 </td>
- <td>4 </td>
- </tr>
- <tr>
- <td>3 </td>
- <td>1.60421e-3</td>
- <td>4 </td>
- <td>1.24836e-9 </td>
- <td>3 </td>
- <td>9.11557e-11 </td>
- <td>3 </td>
- <td>3.35933e-11 </td>
- <td>3 </td>
- <td>5.86734e-11 </td>
- <td>2 </td>
- </tr>
- <tr>
- <td>4 </td>
- <td>9.26748e-4 </td>
- <td>4 </td>
- <td>2.75537e-14 </td>
- <td>4 </td>
- <td>1.39986e-14 </td>
- <td>5 </td>
- <td>2.18864e-14 </td>
- <td>5 </td>
- <td>3.38787e-14 </td>
- <td>5 </td>
- </tr>
- <tr>
- <td>5 </td>
- <td>1.34601e-5</td>
- <td>4 </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- </tr>
- <tr>
- <td>6 </td>
- <td>2.5235e-8 </td>
- <td>5 </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- </tr>
- <tr>
- <td>7 </td>
- <td>1.38899e-12 </td>
- <td>4 </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- </tr>
- <tr>
- <td>8 </td>
- <td>4.68224e-15 </td>
- <td>4 </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- </tr>
+</tr>
+<tr>
+ <td>1</td>
+ <td>3.7112e-03</td>
+ <td>5</td>
+ <td>6.4189e-03</td>
+ <td>3</td>
+ <td>2.4338e-03</td>
+ <td>3</td>
+ <td>1.0570e-03</td>
+ <td>3</td>
+ <td>4.9499e-04</td>
+ <td>3</td>
+</tr>
+<tr>
+ <td>2</td>
+ <td>7.0849e-04</td>
+ <td>5.0000e+00</td>
+ <td>9.9458e-04</td>
+ <td>5</td>
+ <td>1.1409e-04</td>
+ <td>6</td>
+ <td>1.3544e-05</td>
+ <td>6</td>
+ <td>1.4171e-06</td>
+ <td>6</td>
+</tr>
+<tr>
+ <td>3</td>
+ <td>1.9980e-05</td>
+ <td>5.0000e+00</td>
+ <td>4.5007e-05</td>
+ <td>5</td>
+ <td>2.9020e-08</td>
+ <td>5</td>
+ <td>4.4021e-10</td>
+ <td>6</td>
+ <td>6.3435e-11</td>
+ <td>6</td>
+</tr>
+<tr>
+ <td>4</td>
+ <td>2.3165e-09</td>
+ <td>6.0000e+00</td>
+ <td>1.6891e-07</td>
+ <td>5</td>
+ <td>1.2338e-14</td>
+ <td>7</td>
+ <td>1.8506e-14</td>
+ <td>8</td>
+ <td>8.8563e-15</td>
+ <td>8</td>
+</tr>
+<tr>
+ <td>5</td>
+ <td>1.2585e-13</td>
+ <td>7.0000e+00</td>
+ <td>1.4520e-11</td>
+ <td>6</td>
+ <td>1.9044e-13</td>
+ <td>8</td>
+ <td></td>
+ <td></td>
+ <td></td>
+ <td></td>
+</tr>
+<tr>
+ <td>6</td>
+ <td></td>
+ <td></td>
+ <td>1.3998e-15</td>
+ <td>8</td>
+ <td></td>
+ <td></td>
+ <td></td>
+ <td></td>
+ <td></td>
+ <td></td>
+</tr>
</table>
+
+
+
+
+
The following figures show the sequence of generated grids. For the case
of Re=400, the initial guess is obtained by solving Stokes on an $8 \times 8$
mesh, and the mesh is refined adaptively. Between meshes, the solution from
We also show the residual from each step of Newton's iteration on every
mesh. The quadratic convergence is clearly visible in the table.
-<table align="center" border="1">
+<table align="center" class="doxtable">
<tr>
- <th> </th>
+ <th>Re=7500</th>
<th colspan="2">Mesh0</th>
<th colspan="2">Mesh1</th>
<th colspan="2">Mesh2</th>
<th>Residual </th>
<th>FGMRES </th>
</tr>
- <tr>
- <td>1 </td>
- <td>1.89223e-6 </td>
- <td>6 </td>
- <td>4.2506e-3 </td>
- <td>3 </td>
- <td>1.42993e-3 </td>
- <td>3 </td>
- <td>4.87932e-4 </td>
- <td>2 </td>
- <td>1.89981e-04 </td>
- <td>2 </td>
- </tr>
- <tr>
- <td>2 </td>
- <td>3.16439e-9</td>
- <td>8 </td>
- <td>1.3732e-3 </td>
- <td>7 </td>
- <td>4.15062e-4 </td>
- <td>7 </td>
- <td>9.11191e-5 </td>
- <td>8 </td>
- <td>1.35553e-5</td>
- <td>8 </td>
- </tr>
- <tr>
- <td>3 </td>
- <td>1.7628e-14</td>
- <td>9 </td>
- <td>2.19455e-4 </td>
- <td>6 </td>
- <td>1.78805e-5 </td>
- <td>6 </td>
- <td>5.26782e-7 </td>
- <td>7 </td>
- <td>9.37391e-9 </td>
- <td>7 </td>
- </tr>
- <tr>
- <td>4 </td>
- <td> </td>
- <td> </td>
- <td>8.82693e-6 </td>
- <td>6 </td>
- <td>6.82096e-9 </td>
- <td>7 </td>
- <td>2.27696e-11 </td>
- <td>8 </td>
- <td>1.25899e-13</td>
- <td>9 </td>
- </tr>
- <tr>
- <td>5 </td>
- <td> </td>
- <td> </td>
- <td>1.29739e-7</td>
- <td>7 </td>
- <td>1.25167e-13 </td>
- <td>9 </td>
- <td>1.76128e-14 </td>
- <td>10 </td>
- <td> </td>
- <td> </td>
- </tr>
- <tr>
- <td>6 </td>
- <td> </td>
- <td> </td>
- <td>4.43518e-11</td>
- <td>7 </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- </tr>
- <tr>
- <td>7 </td>
- <td> </td>
- <td> </td>
- <td>6.42323e-15 </td>
- <td>9 </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- <td> </td>
- </tr>
+<tr>
+ <td>1</td>
+ <td>1.8922e-06</td>
+ <td>6</td>
+ <td>4.2506e-03</td>
+ <td>3</td>
+ <td>1.4299e-03</td>
+ <td>3</td>
+ <td>4.8793e-04</td>
+ <td>2</td>
+ <td>1.8998e-04</td>
+ <td>2</td>
+</tr>
+<tr>
+ <td>2</td>
+ <td>3.1644e-09</td>
+ <td>8</td>
+ <td>1.3732e-03</td>
+ <td>7</td>
+ <td>4.1506e-04</td>
+ <td>7</td>
+ <td>9.1119e-05</td>
+ <td>8</td>
+ <td>1.3555e-05</td>
+ <td>8</td>
+</tr>
+<tr>
+ <td>3</td>
+ <td>1.7611e-14</td>
+ <td>9</td>
+ <td>2.1946e-04</td>
+ <td>6</td>
+ <td>1.7881e-05</td>
+ <td>6</td>
+ <td>5.2678e-07</td>
+ <td>7</td>
+ <td>9.3739e-09</td>
+ <td>7</td>
+</tr>
+<tr>
+ <td>4</td>
+ <td></td>
+ <td></td>
+ <td>8.8269e-06</td>
+ <td>6</td>
+ <td>6.8210e-09</td>
+ <td>7</td>
+ <td>2.2770e-11</td>
+ <td>8</td>
+ <td>1.2588e-13</td>
+ <td>9</td>
+</tr>
+<tr>
+ <td>5</td>
+ <td></td>
+ <td></td>
+ <td>1.2974e-07</td>
+ <td>7</td>
+ <td>1.2515e-13</td>
+ <td>9</td>
+ <td>1.7801e-14</td>
+ <td>1</td>
+ <td></td>
+ <td></td>
+</tr>
+<tr>
+ <td>6</td>
+ <td></td>
+ <td></td>
+ <td>4.4352e-11</td>
+ <td>7</td>
+ <td></td>
+ <td></td>
+ <td></td>
+ <td></td>
+ <td></td>
+ <td></td>
+</tr>
+<tr>
+ <td>7</td>
+ <td></td>
+ <td></td>
+ <td>6.2863e-15</td>
+ <td>9</td>
+ <td></td>
+ <td></td>
+ <td></td>
+ <td></td>
+ <td></td>
+ <td></td>
+</tr>
</table>
+
+
+
+
+
+
The sequence of generated grids looks like this:
<TABLE ALIGN="center">
<tr>
private:
void setup_dofs();
void initialize_system();
- void assemble_system(const bool initial_step,
- const bool assemble_matrix,
- const bool assemble_rhs);
- void assemble_matrix(const bool initial_step);
+ void assemble(const bool initial_step,
+ const bool assemble_matrix);
+ void assemble_system(const bool initial_step);
void assemble_rhs(const bool initial_step);
void solve(const bool initial_step);
void refine_mesh();
system_rhs.reinit (dofs_per_block);
}
- // @sect4{StationaryNavierStokes::assemble_system}
+ // @sect4{StationaryNavierStokes::assemble}
// This function builds the system matrix and right hand side that we
- // actually work on. "initial_step" is given for applying different
- // constraints (nonzero for the initial step and zero for the others). The
- // other two flags are to determine whether to assemble the system matrix
- // or the right hand side vector, respectively.
+ // currently work on. The @p initial_step argument is used to determine
+ // which set of constraints we apply (nonzero for the initial step and zero
+ // for the others). The @p assemble_matrix flag determines whether to
+ // assemble the whole system or only the right hand side vector,
+ // respectively.
template <int dim>
- void StationaryNavierStokes<dim>::assemble_system(const bool initial_step,
- const bool assemble_matrix,
- const bool assemble_rhs)
+ void StationaryNavierStokes<dim>::assemble(const bool initial_step,
+ const bool assemble_matrix)
{
if (assemble_matrix)
system_matrix = 0;
- if (assemble_rhs)
- system_rhs = 0;
+ system_rhs = 0;
QGauss<dim> quadrature_formula(degree+2);
fe_values[pressure].get_function_values(evaluation_point,
present_pressure_values);
- // The assembly is similar to step-22. An additional term with gamma as a coefficient
- // is the Augmented Lagrangian (AL), which is assembled via grad-div stabilization.
- // As we discussed in the introduction, the bottom right block of the system matrix should be
- // zero. Since the pressure mass matrix is used while creating the preconditioner,
- // we assemble it here and then move it into a separate SparseMatrix at the end (same as in step-22).
-
+ // The assembly is similar to step-22. An additional term with gamma
+ // as a coefficient is the Augmented Lagrangian (AL), which is
+ // assembled via grad-div stabilization. As we discussed in the
+ // introduction, the bottom right block of the system matrix should be
+ // zero. Since the pressure mass matrix is used while creating the
+ // preconditioner, we assemble it here and then move it into a
+ // separate SparseMatrix at the end (same as in step-22).
for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int k=0; k<dofs_per_cell; ++k)
}
}
- if (assemble_rhs)
- {
- double present_velocity_divergence = trace(present_velocity_gradients[q]);
- local_rhs(i) += ( - viscosity*scalar_product(present_velocity_gradients[q],grad_phi_u[i])
- - present_velocity_gradients[q]*present_velocity_values[q]*phi_u[i]
- + present_pressure_values[q]*div_phi_u[i]
- + present_velocity_divergence*phi_p[i]
- - gamma*present_velocity_divergence*div_phi_u[i])
- * fe_values.JxW(q);
- }
+ double present_velocity_divergence = trace(present_velocity_gradients[q]);
+ local_rhs(i) += ( - viscosity*scalar_product(present_velocity_gradients[q],grad_phi_u[i])
+ - present_velocity_gradients[q]*present_velocity_values[q]*phi_u[i]
+ + present_pressure_values[q]*div_phi_u[i]
+ + present_velocity_divergence*phi_p[i]
+ - gamma*present_velocity_divergence*div_phi_u[i])
+ * fe_values.JxW(q);
}
}
cell->get_dof_indices (local_dof_indices);
const ConstraintMatrix &constraints_used = initial_step ? nonzero_constraints : zero_constraints;
- // Finally we move pressure mass matrix into a separate matrix:
if (assemble_matrix)
{
constraints_used.distribute_local_to_global(local_matrix,
+ local_rhs,
local_dof_indices,
- system_matrix);
+ system_matrix,
+ system_rhs);
}
-
- if (assemble_rhs)
+ else
{
constraints_used.distribute_local_to_global(local_rhs,
local_dof_indices,
if (assemble_matrix)
{
+ // Finally we move pressure mass matrix into a separate matrix:
+
pressure_mass_matrix.reinit(sparsity_pattern.block(1,1));
pressure_mass_matrix.copy_from(system_matrix.block(1,1));
}
template <int dim>
- void StationaryNavierStokes<dim>::assemble_matrix(const bool initial_step)
+ void StationaryNavierStokes<dim>::assemble_system(const bool initial_step)
{
- assemble_system(initial_step, true, false);
+ assemble(initial_step, true);
}
template <int dim>
void StationaryNavierStokes<dim>::assemble_rhs(const bool initial_step)
{
- assemble_system(initial_step, false, true);
+ assemble(initial_step, false);
}
+
// @sect4{StationaryNavierStokes::solve}
// In this function, we use FGMRES together with the block preconditioner,
// which is defined at the beginning of the program, to solve the linear
setup_dofs();
initialize_system();
evaluation_point = present_solution;
- assemble_matrix(first_step);
- assemble_rhs(first_step);
+ assemble_system(first_step);
solve(first_step);
present_solution = newton_update;
nonzero_constraints.distribute(present_solution);
else
{
evaluation_point = present_solution;
- assemble_matrix(first_step);
- if (outer_iteration == 0)
- assemble_rhs(first_step);
+ assemble_system(first_step);
solve(first_step);
// To make sure our solution is getting close to the exact solution, we
projects / dealii.git / commitdiff