diff --git a/include/aspect/simulator.h b/include/aspect/simulator.h
index 4d5c99386d4..c28cbf9d5b7 100644
--- a/include/aspect/simulator.h
+++ b/include/aspect/simulator.h
@@ -1300,14 +1300,6 @@ namespace aspect
* "physical" pressure so that all following postprocessing
* steps can use the latter.
*
- * In the case of the surface average, whether a face is part of
- * the surface is determined by asking whether its depth of its
- * midpoint (as determined by the geometry model) is less than
- * 1/3*1/sqrt(dim-1)*diameter of the face. For reasonably curved
- * boundaries, this rules out side faces that are perpendicular
- * to the surface boundary but includes those faces that are
- * along the boundary even if the real boundary is curved.
- *
* Whether the pressure should be normalized based on the
* surface or volume average is decided by a parameter in the
* input file.
@@ -1352,12 +1344,7 @@ namespace aspect
* come out of GMRES, namely the one on which we later called
* normalize_pressure().
*
- * This function modifies @p vector in-place. In some cases, we need
- * locally_relevant values of the pressure. To avoid creating a new vector
- * and transferring data, this function uses a second vector with relevant
- * dofs (@p relevant_vector) for accessing these pressure values. Both
- * @p vector and @p relevant_vector are expected to already contain
- * the correct pressure values.
+ * This function modifies @p vector in-place.
*
* @note The adjustment made in this function is done using the
* negative of the @p pressure_adjustment function argument that
@@ -1371,8 +1358,7 @@ namespace aspect
* source/simulator/helper_functions.cc.
*/
void denormalize_pressure(const double pressure_adjustment,
- LinearAlgebra::BlockVector &vector,
- const LinearAlgebra::BlockVector &relevant_vector) const;
+ LinearAlgebra::BlockVector &vector) const;
/**
* Apply the bound preserving limiter to the discontinuous Galerkin solutions:
@@ -1780,7 +1766,7 @@ namespace aspect
* Computes the initial Newton residual.
*/
double
- compute_initial_newton_residual (const LinearAlgebra::BlockVector &linearized_stokes_initial_guess);
+ compute_initial_newton_residual ();
/**
* This function computes the Eisenstat Walker linear tolerance used for the Newton iterations
diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc
index 9cebb8efe48..0671bf7ad9b 100644
--- a/source/simulator/helper_functions.cc
+++ b/source/simulator/helper_functions.cc
@@ -1085,8 +1085,7 @@ namespace aspect
void
Simulator::
denormalize_pressure (const double pressure_adjustment,
- LinearAlgebra::BlockVector &vector,
- const LinearAlgebra::BlockVector &relevant_vector) const
+ LinearAlgebra::BlockVector &vector) const
{
if (parameters.pressure_normalization == "no")
return;
@@ -1106,6 +1105,17 @@ namespace aspect
finite_element.base_element(introspection.variable("fluid pressure").base_index).dofs_per_cell
: finite_element.base_element(introspection.base_elements.pressure).dofs_per_cell);
+ // We may touch the same DoF multiple times, so we need to copy the
+ // vector before modifying it to have access to the original value.
+ LinearAlgebra::BlockVector vector_backup;
+ vector_backup.reinit(vector, /* omit_zeroing_entries = */ true);
+ const unsigned int pressure_block_index =
+ parameters.include_melt_transport ?
+ introspection.variable("fluid pressure").block_index
+ :
+ introspection.block_indices.pressure;
+ vector_backup.block(pressure_block_index) = vector.block(pressure_block_index);
+
std::vector local_dof_indices (finite_element.dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
@@ -1117,11 +1127,15 @@ namespace aspect
= finite_element.component_to_system_index(pressure_component,
/*dof index within component=*/ j);
- // then adjust its value. Note that because we end up touching
+ // Then adjust its value. vector could be a vector with ghost elements
+ // or a fully distributed vector. In the latter case only access dofs that are
+ // locally owned. Note that because we end up touching
// entries more than once, we are not simply incrementing
- // distributed_vector but copy from the unchanged vector.
- vector(local_dof_indices[local_dof_index])
- = relevant_vector(local_dof_indices[local_dof_index]) - pressure_adjustment;
+ // vector but copy from the vector_backup copy.
+ if (vector.has_ghost_elements() ||
+ dof_handler.locally_owned_dofs().is_element(local_dof_indices[local_dof_index]))
+ vector(local_dof_indices[local_dof_index])
+ = vector_backup(local_dof_indices[local_dof_index]) - pressure_adjustment;
}
}
vector.compress(VectorOperation::insert);
@@ -1144,7 +1158,20 @@ namespace aspect
ExcInternalError());
Assert(!parameters.include_melt_transport, ExcNotImplemented());
const unsigned int pressure_component = introspection.component_indices.pressure;
+
+ // We may touch the same DoF multiple times, so we need to copy the
+ // vector before modifying it to have access to the original value.
+ LinearAlgebra::BlockVector vector_backup;
+ vector_backup.reinit(vector, /* omit_zeroing_entries = */ true);
+ const unsigned int pressure_block_index =
+ parameters.include_melt_transport ?
+ introspection.variable("fluid pressure").block_index
+ :
+ introspection.block_indices.pressure;
+ vector_backup.block(pressure_block_index) = vector.block(pressure_block_index);
+
std::vector local_dof_indices (finite_element.dofs_per_cell);
+
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
@@ -1154,12 +1181,11 @@ namespace aspect
= finite_element.component_to_system_index (pressure_component, 0);
// make sure that this DoF is really owned by the current processor
- // and that it is in fact a pressure dof
Assert (dof_handler.locally_owned_dofs().is_element(local_dof_indices[first_pressure_dof]),
ExcInternalError());
// then adjust its value
- vector (local_dof_indices[first_pressure_dof]) = relevant_vector(local_dof_indices[first_pressure_dof])
+ vector (local_dof_indices[first_pressure_dof]) = vector_backup(local_dof_indices[first_pressure_dof])
- pressure_adjustment;
}
@@ -1341,15 +1367,7 @@ namespace aspect
else
linearized_stokes_variables.block (pressure_block_index) = current_linearization_point.block (pressure_block_index);
- // TODO: we don't have .stokes_relevant_partitioning so I am creating a much
- // bigger vector here, oh well.
- LinearAlgebra::BlockVector ghosted (introspection.index_sets.system_partitioning,
- introspection.index_sets.system_relevant_partitioning,
- mpi_communicator);
- // TODO for Timo: can we create the ghost vector inside of denormalize_pressure
- // (only in cases where we need it)
- ghosted.block(pressure_block_index) = linearized_stokes_variables.block(pressure_block_index);
- denormalize_pressure (this->last_pressure_normalization_adjustment, linearized_stokes_variables, ghosted);
+ denormalize_pressure (this->last_pressure_normalization_adjustment, linearized_stokes_variables);
current_constraints.set_zero (linearized_stokes_variables);
linearized_stokes_variables.block (pressure_block_index) /= pressure_scaling;
@@ -2586,23 +2604,22 @@ namespace aspect
template
double
- Simulator::compute_initial_newton_residual(const LinearAlgebra::BlockVector &linearized_stokes_initial_guess)
+ Simulator::compute_initial_newton_residual()
{
- // Store the values of the current_linearization_point and linearized_stokes_initial_guess so we can reset them again.
+ // Store the values of current_linearization_point to be able to restore it later.
LinearAlgebra::BlockVector temp_linearization_point = current_linearization_point;
- LinearAlgebra::BlockVector temp_linearized_stokes_initial_guess = linearized_stokes_initial_guess;
- // Set the velocity initial guess to zero, but we use the initial guess for the pressure.
+ // Set the velocity initial guess to zero.
current_linearization_point.block(introspection.block_indices.velocities) = 0;
- temp_linearized_stokes_initial_guess.block (introspection.block_indices.velocities) = 0;
- denormalize_pressure (last_pressure_normalization_adjustment,
- temp_linearized_stokes_initial_guess,
- current_linearization_point);
-
- // rebuild the whole system to compute the rhs.
+ // Rebuild the whole system to compute the rhs.
assemble_newton_stokes_system = true;
rebuild_stokes_preconditioner = false;
+
+ // Technically we only need the rhs, but we have asserts in place that check if
+ // the system is assembled correctly when boundary conditions are prescribed, so we assemble the whole system.
+ // TODO: This is a waste of time in the first nonlinear iteration. Check if we can modify the asserts in the
+ // assemble_stokes_system() function to only assemble the RHS.
rebuild_stokes_matrix = boundary_velocity_manager.get_active_boundary_velocity_conditions().size()!=0;
assemble_newton_stokes_matrix = boundary_velocity_manager.get_active_boundary_velocity_conditions().size()!=0;
@@ -2740,8 +2757,7 @@ namespace aspect
template struct Simulator::AdvectionField; \
template double Simulator::normalize_pressure(LinearAlgebra::BlockVector &vector) const; \
template void Simulator::denormalize_pressure(const double pressure_adjustment, \
- LinearAlgebra::BlockVector &vector, \
- const LinearAlgebra::BlockVector &relevant_vector) const; \
+ LinearAlgebra::BlockVector &vector) const; \
template double Simulator::compute_pressure_scaling_factor () const; \
template double Simulator::get_maximal_velocity (const LinearAlgebra::BlockVector &solution) const; \
template std::pair Simulator::get_extrapolated_advection_field_range (const AdvectionField &advection_field) const; \
@@ -2766,7 +2782,7 @@ namespace aspect
template void Simulator::replace_outflow_boundary_ids(const unsigned int boundary_id_offset); \
template void Simulator::restore_outflow_boundary_ids(const unsigned int boundary_id_offset); \
template void Simulator::check_consistency_of_boundary_conditions() const; \
- template double Simulator::compute_initial_newton_residual(const LinearAlgebra::BlockVector &linearized_stokes_initial_guess); \
+ template double Simulator::compute_initial_newton_residual(); \
template double Simulator::compute_Eisenstat_Walker_linear_tolerance(const bool EisenstatWalkerChoiceOne, \
const double maximum_linear_stokes_solver_tolerance, \
const double linear_stokes_solver_tolerance, \
diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc
index 219a5c262ff..f6a58fe8fe0 100644
--- a/source/simulator/solver.cc
+++ b/source/simulator/solver.cc
@@ -855,8 +855,7 @@ namespace aspect
// initial residual we could skip all of this stuff.
distributed_stokes_solution.block(velocity_and_pressure_block) = solution.block(velocity_and_pressure_block);
denormalize_pressure (this->last_pressure_normalization_adjustment,
- distributed_stokes_solution,
- solution);
+ distributed_stokes_solution);
current_stokes_constraints.set_zero (distributed_stokes_solution);
// Undo the pressure scaling:
@@ -973,8 +972,7 @@ namespace aspect
linearized_stokes_initial_guess.block (pressure_block_index) = current_linearization_point.block (pressure_block_index);
denormalize_pressure (this->last_pressure_normalization_adjustment,
- linearized_stokes_initial_guess,
- current_linearization_point);
+ linearized_stokes_initial_guess);
}
else
{
diff --git a/source/simulator/solver_schemes.cc b/source/simulator/solver_schemes.cc
index be2268b99e0..d46f7a64752 100644
--- a/source/simulator/solver_schemes.cc
+++ b/source/simulator/solver_schemes.cc
@@ -437,19 +437,13 @@ namespace aspect
: introspection.block_indices.pressure;
Assert(introspection.block_indices.velocities == 0, ExcNotImplemented());
Assert(pressure_block_index == 1, ExcNotImplemented());
+ (void) pressure_block_index;
Assert(!parameters.include_melt_transport
|| introspection.variable("compaction pressure").block_index == 1, ExcNotImplemented());
- // create a completely distributed vector that will be used for
- // the scaled and denormalized solution and later used as a
- // starting guess for the linear solver
- LinearAlgebra::BlockVector linearized_stokes_initial_guess(introspection.index_sets.stokes_partitioning, mpi_communicator);
- linearized_stokes_initial_guess.block(introspection.block_indices.velocities) = current_linearization_point.block(introspection.block_indices.velocities);
- linearized_stokes_initial_guess.block(pressure_block_index) = current_linearization_point.block(pressure_block_index);
-
if (nonlinear_iteration == 0)
{
- dcr.initial_residual = compute_initial_newton_residual(linearized_stokes_initial_guess);
+ dcr.initial_residual = compute_initial_newton_residual();
dcr.switch_initial_residual = dcr.initial_residual;
dcr.residual_old = dcr.initial_residual;
dcr.residual = dcr.initial_residual;
@@ -461,12 +455,6 @@ namespace aspect
{
assemble_newton_stokes_system = assemble_newton_stokes_matrix = false;
}
- else
- {
- denormalize_pressure (last_pressure_normalization_adjustment,
- linearized_stokes_initial_guess,
- current_linearization_point);
- }
if (nonlinear_iteration <= 1)
{
diff --git a/source/simulator/stokes_matrix_free.cc b/source/simulator/stokes_matrix_free.cc
index 7a893f2555e..6ade5af78e7 100644
--- a/source/simulator/stokes_matrix_free.cc
+++ b/source/simulator/stokes_matrix_free.cc
@@ -1919,8 +1919,7 @@ namespace aspect
linearized_stokes_initial_guess.block (block_p) = sim.current_linearization_point.block (block_p);
sim.denormalize_pressure (sim.last_pressure_normalization_adjustment,
- linearized_stokes_initial_guess,
- sim.current_linearization_point);
+ linearized_stokes_initial_guess);
}
else
{