Skip to content

Commit

Permalink
Merge pull request #1976 from gassmoeller/improve_compressible_precon…
Browse files Browse the repository at this point in the history 
…ditioner

Improve compressible preconditioner
  • Loading branch information
@gassmoeller
gassmoeller committed Oct 30, 2017
2 parents 65ea661 + bb87b60 commit d06f1f1
Show file tree
Hide file tree
Showing 108 changed files with 10,727 additions and 10,556 deletions.
5 changes: 3 additions & 2 deletions doc/modules/changes/20171027_bangerth
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
Fixed: A forgotten factor of 2 in assembling the preconditioner for
the Stokes equation led to a suboptimal performance of the solver. The
the Stokes equation, and a forgotten term for compressible models,
led to a suboptimal performance of the solver. The
answers were always correct, but it took more iterations than
necessary. This is now fixed, and results in significantly fewer
GMRES iterations when solving the Stokes equations.
<br>
(Menno Fraters, Wolfgang Bangerth, 2017/10/27)
(Menno Fraters, Wolfgang Bangerth, Rene Gassmoeller, Timo Heister, 2017/10/27)
10 changes: 9 additions & 1 deletion include/aspect/assembly.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,7 @@ namespace aspect
std::vector<types::global_dof_index> local_dof_indices;
std::vector<unsigned int> dof_component_indices;
std::vector<SymmetricTensor<2,dim> > grads_phi_u;
std::vector<double> div_phi_u;
std::vector<double> phi_p;
std::vector<double> phi_p_c;
std::vector<Tensor<1,dim> > grad_phi_p;
Expand Down Expand Up @@ -103,7 +104,6 @@ namespace aspect

std::vector<Tensor<1,dim> > phi_u;
std::vector<SymmetricTensor<2,dim> > grads_phi_u;
std::vector<double> div_phi_u;
std::vector<Tensor<1,dim> > velocity_values;
std::vector<double> velocity_divergence;

Expand Down Expand Up @@ -736,6 +736,14 @@ namespace aspect
internal::Assembly::Scratch::StokesPreconditioner<dim> &scratch,
internal::Assembly::CopyData::StokesPreconditioner<dim> &data) const;

/**
* This function assembles the terms of the Stokes preconditioner matrix for the current cell.
*/
void
compressible_preconditioner (const double pressure_scaling,
internal::Assembly::Scratch::StokesPreconditioner<dim> &scratch,
internal::Assembly::CopyData::StokesPreconditioner<dim> &data) const;

/**
* This function assembles the terms for the matrix and right-hand-side of the incompressible
* Stokes equation for the current cell.
Expand Down
64 changes: 64 additions & 0 deletions source/simulator/assemblers/stokes.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,6 +122,70 @@ namespace aspect



template <int dim>
void
StokesAssembler<dim>::
compressible_preconditioner (const double /*pressure_scaling*/,
internal::Assembly::Scratch::StokesPreconditioner<dim> &scratch,
internal::Assembly::CopyData::StokesPreconditioner<dim> &data) const
{
const Introspection<dim> &introspection = this->introspection();
const FiniteElement<dim> &fe = this->get_fe();
const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points;

// First loop over all dofs and find those that are in the Stokes system
// save the component (pressure and dim velocities) each belongs to.
for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/)
{
if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
{
scratch.dof_component_indices[i_stokes] = fe.system_to_component_index(i).first;
++i_stokes;
}
++i;
}

// Loop over all quadrature points and assemble their contributions to
// the preconditioner matrix
for (unsigned int q = 0; q < n_q_points; ++q)
{
for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/)
{
if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
{
scratch.grads_phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].symmetric_gradient(i,q);
scratch.div_phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].divergence (i, q);

++i_stokes;
}
++i;
}

const double eta_two_thirds = scratch.material_model_outputs.viscosities[q] * 2.0 / 3.0;

const SymmetricTensor<4, dim> &stress_strain_director = scratch
.material_model_outputs.stress_strain_directors[q];
const bool use_tensor = (stress_strain_director
!= dealii::identity_tensor<dim>());

const double JxW = scratch.finite_element_values.JxW(q);

for (unsigned int i = 0; i < stokes_dofs_per_cell; ++i)
for (unsigned int j = 0; j < stokes_dofs_per_cell; ++j)
if (scratch.dof_component_indices[i] ==
scratch.dof_component_indices[j])
data.local_matrix(i, j) += (- (use_tensor ?
eta_two_thirds * (scratch.div_phi_u[i] * trace(stress_strain_director * scratch.grads_phi_u[j]))
:
eta_two_thirds * (scratch.div_phi_u[i] * scratch.div_phi_u[j])
))
* JxW;
}
}



template <int dim>
void
StokesAssembler<dim>::
Expand Down
20 changes: 14 additions & 6 deletions source/simulator/assembly.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,7 @@ namespace aspect
local_dof_indices (finite_element.dofs_per_cell),
dof_component_indices(stokes_dofs_per_cell),
grads_phi_u (stokes_dofs_per_cell, numbers::signaling_nan<SymmetricTensor<2,dim> >()),
div_phi_u (stokes_dofs_per_cell, numbers::signaling_nan<double>()),
phi_p (stokes_dofs_per_cell, numbers::signaling_nan<double>()),
phi_p_c (add_compaction_pressure ? stokes_dofs_per_cell : 0, numbers::signaling_nan<double>()),
grad_phi_p (add_compaction_pressure ? stokes_dofs_per_cell : 0, numbers::signaling_nan<Tensor<1,dim> >()),
Expand All @@ -84,6 +85,7 @@ namespace aspect
local_dof_indices (scratch.local_dof_indices),
dof_component_indices( scratch.dof_component_indices),
grads_phi_u (scratch.grads_phi_u),
div_phi_u (scratch.div_phi_u),
phi_p (scratch.phi_p),
phi_p_c (scratch.phi_p_c),
grad_phi_p(scratch.grad_phi_p),
Expand Down Expand Up @@ -127,7 +129,6 @@ namespace aspect

phi_u (stokes_dofs_per_cell, numbers::signaling_nan<Tensor<1,dim> >()),
grads_phi_u (stokes_dofs_per_cell, numbers::signaling_nan<SymmetricTensor<2,dim> >()),
div_phi_u (stokes_dofs_per_cell, numbers::signaling_nan<double>()),
velocity_values (quadrature.size(), numbers::signaling_nan<Tensor<1,dim> >()),
velocity_divergence (quadrature.size(), numbers::signaling_nan<double>()),
face_material_model_inputs(face_quadrature.size(), n_compositional_fields),
Expand All @@ -151,7 +152,6 @@ namespace aspect

phi_u (scratch.phi_u),
grads_phi_u (scratch.grads_phi_u),
div_phi_u (scratch.div_phi_u),
velocity_values (scratch.velocity_values),
velocity_divergence (scratch.velocity_divergence),
face_material_model_inputs(scratch.face_material_model_inputs),
Expand Down Expand Up @@ -667,10 +667,18 @@ namespace aspect
std_cxx11::_2,
std_cxx11::_3));
else
assemblers->local_assemble_stokes_preconditioner
.connect (std_cxx11::bind(&aspect::Assemblers::StokesAssembler<dim>::preconditioner,
std_cxx11::cref (*stokes_assembler),
std_cxx11::_1, std_cxx11::_2, std_cxx11::_3));
{
assemblers->local_assemble_stokes_preconditioner
.connect (std_cxx11::bind(&aspect::Assemblers::StokesAssembler<dim>::preconditioner,
std_cxx11::cref (*stokes_assembler),
std_cxx11::_1, std_cxx11::_2, std_cxx11::_3));

if (material_model->is_compressible())
assemblers->local_assemble_stokes_preconditioner
.connect (std_cxx11::bind(&aspect::Assemblers::StokesAssembler<dim>::compressible_preconditioner,
std_cxx11::cref (*stokes_assembler),
std_cxx11::_1, std_cxx11::_2, std_cxx11::_3));
}


if (parameters.include_melt_transport)
Expand Down
11 changes: 9 additions & 2 deletions source/simulator/melt.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,6 +152,7 @@ namespace aspect
if (is_velocity_or_pressures(introspection,p_c_component_index,p_f_component_index,component_index_i))
{
scratch.grads_phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].symmetric_gradient(i,q);
scratch.div_phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].divergence (i, q);
scratch.phi_p[i_stokes] = scratch.finite_element_values[ex_p_f].value (i, q);
scratch.phi_p_c[i_stokes] = scratch.finite_element_values[ex_p_c].value (i, q);
scratch.grad_phi_p[i_stokes] = scratch.finite_element_values[ex_p_f].gradient (i, q);
Expand All @@ -162,6 +163,7 @@ namespace aspect

const double eta = scratch.material_model_outputs.viscosities[q];
const double one_over_eta = 1. / eta;
const double eta_two_thirds = scratch.material_model_outputs.viscosities[q] * 2.0 / 3.0;

/*
- R = 1/eta M_p + K_D L_p for p
Expand All @@ -187,9 +189,14 @@ namespace aspect
for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
if (scratch.dof_component_indices[i] == scratch.dof_component_indices[j])
data.local_matrix(i,j) += ((use_tensor ?
eta * (scratch.grads_phi_u[i] * stress_strain_director * scratch.grads_phi_u[j])
2.0 * eta * (scratch.grads_phi_u[i] * stress_strain_director * scratch.grads_phi_u[j])
:
eta * (scratch.grads_phi_u[i] * scratch.grads_phi_u[j]))
2.0 * eta * (scratch.grads_phi_u[i] * scratch.grads_phi_u[j]))
-
(use_tensor ?
eta_two_thirds * (scratch.div_phi_u[i] * trace(stress_strain_director * scratch.grads_phi_u[j]))
:
eta_two_thirds * (scratch.div_phi_u[i] * scratch.div_phi_u[j]))
+
(one_over_eta *
pressure_scaling *
Expand Down
22 changes: 11 additions & 11 deletions tests/adiabatic_heating_with_melt/screen-output
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,20 +7,20 @@ Number of degrees of freedom: 9,270 (2,898+810+2,898+405+1,449+405+405)
Solving porosity system ... 0 iterations.
Skipping peridotite composition solve because RHS is zero.
Rebuilding Stokes preconditioner...
Solving Stokes system... 57+0 iterations.
Solving Stokes system... 54+0 iterations.
Solving for u_f in 1 iterations.
Relative nonlinear residuals: 1.80908e-16, 9.58741e-17, 0, 1
Relative nonlinear residuals: 1.85418e-16, 9.58741e-17, 0, 1
Total relative residual after nonlinear iteration 1: 1


Solving temperature system... 0 iterations.
Solving porosity system ... 0 iterations.
Skipping peridotite composition solve because RHS is zero.
Rebuilding Stokes preconditioner...
Solving Stokes system... 4+0 iterations.
Solving Stokes system... 5+0 iterations.
Solving for u_f in 1 iterations.
Relative nonlinear residuals: 1.80908e-16, 9.58741e-17, 0, 6.74528e-13
Total relative residual after nonlinear iteration 2: 6.74528e-13
Relative nonlinear residuals: 1.85418e-16, 9.58741e-17, 0, 9.29021e-13
Total relative residual after nonlinear iteration 2: 9.29021e-13



Expand All @@ -33,9 +33,9 @@ Number of degrees of freedom: 9,270 (2,898+810+2,898+405+1,449+405+405)
Solving porosity system ... 0 iterations.
Skipping peridotite composition solve because RHS is zero.
Rebuilding Stokes preconditioner...
Solving Stokes system... 18+0 iterations.
Solving Stokes system... 17+0 iterations.
Solving for u_f in 1 iterations.
Relative nonlinear residuals: 2.0391e-06, 1.55118e-16, 0, 4.87025e-09
Relative nonlinear residuals: 2.0391e-06, 1.53408e-16, 0, 4.87025e-09
Total relative residual after nonlinear iteration 1: 2.0391e-06


Expand All @@ -49,9 +49,9 @@ Number of degrees of freedom: 9,270 (2,898+810+2,898+405+1,449+405+405)
Solving porosity system ... 0 iterations.
Skipping peridotite composition solve because RHS is zero.
Rebuilding Stokes preconditioner...
Solving Stokes system... 21+0 iterations.
Solving Stokes system... 20+0 iterations.
Solving for u_f in 1 iterations.
Relative nonlinear residuals: 2.21934e-06, 2.29702e-16, 0, 4.20136e-09
Relative nonlinear residuals: 2.21934e-06, 1.21242e-16, 0, 4.20135e-09
Total relative residual after nonlinear iteration 1: 2.21934e-06


Expand All @@ -65,9 +65,9 @@ Number of degrees of freedom: 9,270 (2,898+810+2,898+405+1,449+405+405)
Solving porosity system ... 0 iterations.
Skipping peridotite composition solve because RHS is zero.
Rebuilding Stokes preconditioner...
Solving Stokes system... 17+0 iterations.
Solving Stokes system... 16+0 iterations.
Solving for u_f in 1 iterations.
Relative nonlinear residuals: 4.7975e-08, 1.75107e-16, 0, 8.5155e-10
Relative nonlinear residuals: 4.7975e-08, 1.56341e-16, 0, 8.51551e-10
Total relative residual after nonlinear iteration 1: 4.7975e-08


Expand Down
8 changes: 4 additions & 4 deletions tests/adiabatic_heating_with_melt/statistics
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@
# 20: Total adiabatic heating rate (W)
# 21: Average adiabatic heating of melt rate (W/kg)
# 22: Total adiabatic heating of melt rate (W)
0 0.000000000000e+00 0.000000000000e+00 320 3303 1449 810 2 0 0 0 61 64 64 1.72315000e+03 1.75516666e+03 1.78757747e+03 8.48573063e-03 -3.19121211e-11 -1.77378538e+02 3.19121211e-11 1.77378538e+02
1 1.147124999998e+13 1.147124999998e+13 320 3303 1449 810 1 8 0 0 18 19 19 1.72319196e+03 1.75516670e+03 1.78757747e+03 8.47462138e-03 -3.19121219e-11 -1.77378542e+02 -4.62946794e-26 -2.57321740e-13
2 2.294249984881e+13 1.147124984883e+13 320 3303 1449 810 1 13 0 0 21 22 22 1.72316615e+03 1.75516665e+03 1.78757747e+03 8.48144718e-03 -3.19121210e-11 -1.77378537e+02 -1.78871215e-26 -9.94227693e-14
3 2.500000000000e+13 2.057500151191e+12 320 3303 1449 810 1 6 0 0 17 18 18 1.72315936e+03 1.75516664e+03 1.78757747e+03 8.48324430e-03 -3.19121208e-11 -1.77378536e+02 -8.33979653e-26 -4.63554556e-13
0 0.000000000000e+00 0.000000000000e+00 320 3303 1449 810 2 0 0 0 59 62 62 1.72315000e+03 1.75516666e+03 1.78757747e+03 8.48573063e-03 -3.19121211e-11 -1.77378538e+02 3.19121211e-11 1.77378538e+02
1 1.147124999994e+13 1.147124999994e+13 320 3303 1449 810 1 8 0 0 17 18 18 1.72319196e+03 1.75516670e+03 1.78757747e+03 8.47462138e-03 -3.19121219e-11 -1.77378542e+02 -1.02517773e-26 -5.69829019e-14
2 2.294249984878e+13 1.147124984883e+13 320 3303 1449 810 1 13 0 0 20 21 21 1.72316615e+03 1.75516665e+03 1.78757747e+03 8.48144718e-03 -3.19121210e-11 -1.77378537e+02 2.35199516e-27 1.30731975e-14
3 2.500000000000e+13 2.057500151224e+12 320 3303 1449 810 1 6 0 0 16 17 17 1.72315936e+03 1.75516664e+03 1.78757747e+03 8.48324430e-03 -3.19121208e-11 -1.77378536e+02 4.13359881e-26 2.29759630e-13

0 comments on commit d06f1f1

Please sign in to comment.