diff --git a/source/mesh_deformation/fastscape.cc b/source/mesh_deformation/fastscape.cc index 29cf61dd8a9..60ef3daa6ce 100644 --- a/source/mesh_deformation/fastscape.cc +++ b/source/mesh_deformation/fastscape.cc @@ -499,7 +499,7 @@ namespace aspect &sand_transport_coefficient, &silt_transport_coefficient); - // generate a combined array for kf kd both onshore and offshore + // generate a combined array for kf and kd both onshore and offshore std::vector combined_kd(fastscape_array_size); std::vector combined_kf(fastscape_array_size); for (unsigned int i = 0; i < fastscape_array_size; ++i) @@ -515,7 +515,7 @@ namespace aspect { combined_kf[i] = sediment_river_incision_rate; combined_kd[i] = sediment_transport_coefficient; - }; + } } // select additional output for Fastscape vtu @@ -558,11 +558,13 @@ namespace aspect mesh_velocity_z[i] = (elevation[i] - elevation_old[i])/aspect_timestep_in_years; } else + // For ranks other than the root: { - // add a declare for kd - std::vector bedrock_transport_coefficient_array(fastscape_array_size); for (unsigned int i=0; iget_mpi_communicator()); + MPI_Ssend(&local_aspect_values[i][0], local_aspect_values[i].size(), MPI_DOUBLE, + /* destination is root= */ 0, + /* tag= */ 42, + this->get_mpi_communicator()); // Check whether the FastScape mesh was filled with data. const bool fastscape_mesh_filled = Utilities::MPI::broadcast (this->get_mpi_communicator(), true, 0); @@ -570,12 +572,13 @@ namespace aspect throw aspect::QuietException(); // This is called solely so we can set the timer and will return immediately. + std::vector dummy(fastscape_array_size); execute_fastscape(mesh_velocity_z, mesh_velocity_z, mesh_velocity_z, mesh_velocity_z, mesh_velocity_z, - bedrock_transport_coefficient_array, + dummy, aspect_timestep_in_years, fastscape_steps_per_aspect_step); } @@ -585,14 +588,15 @@ namespace aspect // that will be interpolated back to ASPECT. We need this table on all // processes, and can achieve this goal by first filling it on the root process, // and then replicating it on all processes (if possible using shared memory). - TableIndices size_idx; - for (unsigned int d=0; d velocity_table; if (Utilities::MPI::this_mpi_process(this->get_mpi_communicator()) == 0) - velocity_table = fill_data_table(mesh_velocity_z, size_idx, fastscape_nx, fastscape_ny); + { + TableIndices size_idx; + for (unsigned int d=0; dget_mpi_communicator(), /*root_process=*/0); @@ -608,9 +612,9 @@ namespace aspect // Then create a function that can be interpolated from the data table. // Use move semantics to ensure that we keep using the replicated // table: - Functions::InterpolatedUniformGridData velocities (std::move(interpolation_extent), - std::move(table_intervals), - std::move(velocity_table)); + const Functions::InterpolatedUniformGridData velocities (std::move(interpolation_extent), + std::move(table_intervals), + std::move(velocity_table)); VectorFunctionFromScalarFunctionObject vector_function_object( [&](const Point &p) -> double @@ -620,10 +624,11 @@ namespace aspect dim-1, dim); - VectorTools::interpolate_boundary_values (mesh_deformation_dof_handler, - *boundary_ids.begin(), - vector_function_object, - mesh_velocity_constraints); + for (const types::boundary_id boundary : boundary_ids) + VectorTools::interpolate_boundary_values (mesh_deformation_dof_handler, + boundary, + vector_function_object, + mesh_velocity_constraints); this->get_computing_timer().leave_subsection("FastScape plugin"); }