Fix two places where we write into vectors with ghosts.

source/postprocess/heat_flux_map.cc

@@ -94,8 +94,10 @@ namespace aspect
         LinearAlgebra::BlockVector rhs_vector(simulator_access.introspection().index_sets.system_partitioning,
                                               simulator_access.get_mpi_communicator());
 
+        // Set both of these vectors to zero in a way that respects that
+        // heat_flux_vector has ghost entries:
         distributed_heat_flux_vector = 0.;
-        heat_flux_vector = 0.;
+        heat_flux_vector = distributed_heat_flux_vector;
 
         typename MaterialModel::Interface<dim>::MaterialModelInputs
         in(fe_volume_values.n_quadrature_points, simulator_access.n_compositional_fields());

source/simulator/helper_functions.cc

@@ -2599,10 +2599,17 @@ namespace aspect
   Simulator<dim>::compute_initial_newton_residual()
   {
     // Store the values of current_linearization_point to be able to restore it later.
-    LinearAlgebra::BlockVector temp_linearization_point = current_linearization_point;
-
-    // Set the velocity initial guess to zero.
-    current_linearization_point.block(introspection.block_indices.velocities) = 0;
+    // We will want to set velocity component of the current_linearization_points
+    // to a zero vector until the end of the function,
+    // which is most efficiently done by creating the
+    // temp_linearization_point as a zero vector, swapping contents, and at the end
+    // of the function, swapping things back:
+    LinearAlgebra::Vector temp_velocity_linearization_point
+    (introspection.index_sets.system_partitioning[introspection.block_indices.velocities],
+     introspection.index_sets.system_relevant_partitioning[introspection.block_indices.velocities],
+     mpi_communicator);
+    temp_velocity_linearization_point
+    .swap(current_linearization_point.block(introspection.block_indices.velocities));
 
     // Rebuild the whole system to compute the rhs.
     assemble_newton_stokes_system = true;
@@ -2623,7 +2630,9 @@ namespace aspect
     const double initial_newton_residual_p = system_rhs.block(introspection.block_indices.pressure).l2_norm();
     const double initial_newton_residual = std::sqrt(initial_newton_residual_vel * initial_newton_residual_vel + initial_newton_residual_p * initial_newton_residual_p);
 
-    current_linearization_point = temp_linearization_point;
+    // Swap old content back in:
+    current_linearization_point.block(introspection.block_indices.velocities)
+    .swap(temp_velocity_linearization_point);
 
     pcout << "   Initial Newton Stokes residual = " << initial_newton_residual << ", v = " << initial_newton_residual_vel << ", p = " << initial_newton_residual_p << std::endl << std::endl;
     return initial_newton_residual;

source/simulator/solver.cc

@@ -422,7 +422,11 @@ namespace aspect
     if (system_rhs.block(block_idx).l2_norm() <= std::numeric_limits<double>::min())
       {
         pcout << "   Skipping " + field_name + " solve because RHS is zero." << std::endl;
-        solution.block(block_idx) = 0;
+
+        const LinearAlgebra::Vector zero_vector_no_ghosts (
+          introspection.index_sets.system_partitioning[block_idx],
+          mpi_communicator);
+        solution.block(block_idx) = zero_vector_no_ghosts;
 
         // signal successful solver and signal residual of zero
         solver_control.check(0, 0.0);